Analysis of Nonlinear Thermoelastic Dissipation in Euler-Bernoulli Beam Resonators

The linear theory of thermoelastic damping (TED) has been extensively developed over the past eight decades, but relatively little is known about the different types of nonlinearities that are associated with this fundamental mechanism of material damping. Here, we initiate the study of a dissipative nonlinearity (also called thermomechanical nonlinearity) whose origins reside at the heart of the thermomechanical coupling that gives rise to TED. The finite difference method is used to solve the nonlinear governing equation and estimate nonlinear TED in Euler-Bernoulli beams. The maximum difference between the nonlinear and linear estimates ranges from 0.06% for quartz and 0.3% for silicon to 7% for aluminum and 28% for zinc.


Introduction
Thermoelastic damping (TED) refers to energy dissipation due to irreversible heat conduction across thermoelastic temperature gradients in vibrating structures [1,2]. The foundations of the linear theory of TED were established in the 1930s, and developed extensively over the past thirty years (see, for instance [3][4][5] for reviews). The linear analysis is now widely used to gain insight into a fundamental mechanism of material damping, calibrate measurements of internal friction in thin films, and guide the design of miniaturized resonators used in microelectromechanical systems (MEMS).
In stark contrast, nonlinear thermoelastic dissipation has received relatively little attention. The literature is sparse and focuses mainly on the following three types of nonlinearities: (i) geometric nonlinearities due to large deformations, (ii) material nonlinearity due to the temperature dependence of mechanical properties, and (iii) transduction nonlinearities in oscillators caused by electrostatic actuation [6][7][8][9][10][11][12][13][14][15]. In this paper, we address a fourth type of nonlinearity whose origins reside at the heart of the thermoelastic coupling that gives rise to dissipation. For this reason, it may be termed a dissipative nonlinearity (or thermoelastic nonlinearity).
To understand the source of the dissipative nonlinearity, let us consider the time-harmonic oscillations of a thermoelastic structure that is initially at equilibrium at temperature T 0 . Due to thermoelastic coupling, the bending vibrations create a time-dependent temperature field T within the structure. In turn, the temperature gradients lead to irreversible heat conduction, entropy generation, and energy dissipation [1].
For diffusive heat conduction according to Fourier's law, the temperature is governed by a nonlinear partial differential equation (PDE) given by [16] C @y @t where θ = T−T 0 is the excess temperature, C is the specific heat per unit volume, t is time, k is the thermal conductivity, E is Young's modulus, α is the coefficient of thermal expansion, υ is Poisson's ratio, and ε is the total strain (that is, the sum of elastic strain and thermal strain). The summation convention is implied for the index n; hence, ε nn denotes the trace of the strain tensor.
The second term on the right-hand side of Eq (1) gives rise to the dissipative nonlinearity. Thus far, this effect has been largely ignored in the literature, and the governing equation is usually linearized by replacing T with T 0 . The error caused by linearization is expected to be small, but this important assumption has yet to be rigorously examined and quantified. In the remaining sections of the paper, we present a detailed numerical analysis of nonlinear TED and quantify the error incurred by ignoring the dissipative nonlinearity.

Methods
Consider an isotropic, homogeneous, Euler-Bernoulli beam of length L, thickness h, and width b. A Cartesian coordinate system (x,y,z) is attached to the structure so that the beam occupies the domain defined by 0 x L, 0 y b, and 0 z h. The axial stress σ xx = σ 0 exp(iωt) is the only nonzero component of the stress tensor; here, ω is the oscillation frequency (in units of radians per second) and i ¼ . The oscillatory curvature of the beam under pure bending is κ xx = z 0 exp(iωt), where z 0 is the amplitude of oscillation. The normal strains are given by [16] and z = h/2 is the position of the neutral axis. Using the expressions for the stresses and strains, and for one-dimensional heat conduction across the thickness of the beam, Eq (1) can be expressed as where C = Eα 2 T 0 /C is the dimensionless Zener modulus. The beam is initially at equilibrium at temperature T 0 and the boundaries along the z-direction are adiabatic; hence, the initial and boundary conditions are given by We employ the method of finite difference to solve the nonlinear partial differential equation, Eq (3), subject to the initial and boundary conditions expressed in Eq (4). The first step is to cast the PDE in dimensionless form by defining the following variables Here, " y is the normalized temperature, " t is the normalized time, " z is the normalized coordinate, and O is the normalized frequency. Hence, Eq (3) can be expressed in dimensionless form as Eq (6) was solved using the implicit Crank-Nicholson finite-difference scheme with discretization grids that are uniform in space (" z) and time ( " t) [17]. All the derivatives are replaced by the corresponding finite difference approximation to get where N ¼ 0:5D " t=D" z 2 and M ¼ i OD " tð" z À 0:5ÞexpðiO " tÞ. The superscripts n + 1 and n denote the temperature at the current and previous iteration, respectively, and the subscripts refer to the nodes of the finite-difference mesh. This equation was expressed in matrix form and solved iteratively using the tri-diagonal matrix algorithm to obtain the temperature field. Finally, by definition, the nonlinear thermoelastic damping is given by [3,18] Here, V is the volume of the beam, H denotes the integral over one cycle of oscillation, ε th = αθ is the thermal strain, and the peak elastic energy in the beam is W ¼ Ez 2 0 h 3 b=24. Taken together, the aforementioned equations constitute a framework for computing nonlinear thermoelastic dissipation in Euler-Bernoulli beam resonators. The framework was implemented using MATLAB to analyze nonlinear TED and the results are presented in the next section.

Results
Nonlinear TED was analyzed in Euler-Bernoulli beams made of five different materials: quartz, silicon, diamond-like carbon (DLC), aluminum, and zinc. The materials were chosen to span the range of thermomechanical properties pertinent to TED. Table 1 lists the effective isotropic material properties at room temperature.
In all cases, the nonlinear estimates were compared with the prediction of the linear theory of TED in Euler-Bernoulli beams given by [20] The main results are presented in Fig 1 on a  , the nonlinear analysis predicts a higher value for TED. The maximum difference between the nonlinear and linear estimates is 0.06% for quartz, 0.3% for both silicon and diamond-like carbon, 7% for aluminum, and 28% for zinc. Taken together, the results indicate that the error caused by ignoring the dissipative nonlinearity scales with the Zener modulus of the material.

Discussion
Thermoelastic damping arises as a consequence of the nonlinear coupling between thermal and elastic fields in vibrating solids. In this paper, we presented the first detailed numerical analysis of the effects of the dissipative nonlinearity on thermoelastic damping in Euler-Bernoulli beams. All other types of nonlinearities (geometric, material, and transduction) were eliminated from the analysis, and the effects of the dissipative nonlinearity were quantified by comparison with a linear model for thermoelastic damping. Results were presented for five representative materials that span the full range of the Zener modulus encountered in metals, alloys, glasses, and ceramics. The dissipative nonlinearity does not introduce any new spectral features, and the dissipation spectra display a single peak. The error caused by ignoring the dissipative nonlinearity is negligible at below the peak frequency for all materials. Thus, the linear approximation can be used with confidence for design and analysis in this regime. However, beyond the peak frequency, the error is proportional to the Zener modulus of the material. The maximum error can be as high as 7% for aluminum and 28% for zinc.
To our knowledge, the present study is the first to quantify the effects of the dissipative nonlinearity on thermoelastic damping. The analysis focused on isotropic and homogeneous Euler-Bernoulli beams, and ignored the effects of material and geometric nonlinearities. Our results suggest several fruitful and intriguing topics for future investigations including the study of nonlinear thermoelastic dissipation in other structures (plates, hollow tubes, rings, shells, and layered composites) and anisotropic materials, and exploring the interactions between material, geometric, and dissipative nonlinearities. A useful starting point for such studies is to replace Eq (1) with a fully nonlinear equation for thermoelastic heat conduction in anisotropic materials [21].