Background Thermoelasticity Theory
When an elastic rod is stretched reversibly and adiabatically its temperature drops. The drop in temperature compensates for the increase in entropy caused by the stress in the rod (since the process is reversible the entropy remains constant). The equations of thermoelasticity incorporate this phenomenon, together with the irreversible processes that occur in a vibrating rod (which result in net entropy increases). When a structure vibrates in a mode with both local compression and expansion there are always some losses associated with the irreversible heat conduction between the expanding areas that cool and the contracting areas that heat. These losses result in thermoelastic damping, which can be important, particularly in small vibrating structures, such as resonant MEMS devices.
References 1 to 5 in the section References for the Thermoelasticity Interface provide useful background information.
The equations of thermoelasticity are derived from the first law of thermodynamics, which can be stated as follows:
(5-7)
where dU is the change in internal energy, dQ' is the heat flow into the system (the prime indicates an inexact differential in this case) and dW' is the work done on the system. For a small part of a solid (sufficiently small that the stresses and strains are uniform), with an initial reference density, ρ0, the first law can be rewritten in the following form (assuming that the differential changes occur between equilibrium states):
(5-8)
where Ta is the absolute temperature, s is the entropy per unit mass, σ is the elastic part of the second Piola-Kirchhoff stress (in general a rank 2 tensor), ε is the material strain (also a tensor). In general the second Piola-Kirchhoff stress tensor, p, must be split into elastic (σ) and inelastic (τ) parts such that:
The elastic part of the stress tensor, σ, does work σ:dε during a change in the strain. The inelastic part of the stress tensor, τ, generates heat at a rate τ:(dε/dτ) when the strain is changing and is identified with internal or material damping. These internal damping mechanisms are associated with microscopic phenomena such as dislocation movement.
From Equation 5-8 it is possible to make the following identifications for Ta and σ:
Next the entropy balance equation must be derived. Since thermoelasticity involves irreversible processes, the assumption of equilibrium required to derive Equation 5-8 is no longer valid. Instead an assumption of “local” equilibrium is made. It is assumed that although the system is not in equilibrium, there exists within small elements a state of local equilibrium, for which the local entropy per unit mass, s, is the same function of the internal energy, strain, and particle number as it was in equilibrium. This assumption is commonly employed in the modeling of transport phenomena and is justified only by the validity of conclusions derived from it and by results obtained from specific microscopic models, for near-equilibrium situations. For a small volume element in the material frame Equation 5-8 can then be written as:
The rate of change of entropy can then be written as:
(5-9)
From the first law (Equation 5-7) the rate of change of internal energy is given by:
where w is the work done per unit volume and q is the heat accumulated per unit volume. The heat accumulated can be written as the sum of the heat sources and the divergence in the material frame heat flux:
where Q represents the heat source per unit volume and τ is the inelastic part of the stress tensor. The rate of doing work (per unit reference volume) by a linear elastic material is given by the elastic part of the second Piola–Kirchhoff stress contracted with the rate of material strain. Per unit volume the following equation is obtained:
so Equation 5-9 reduces to:
(5-10)
The definition of the material thermal conductivity gives:
where κ is the thermal conductivity, defined in the material frame.
Therefore the equation is:
(5-11)
It is now necessary to derive an expression for the rate of change of entropy with respect to time. In order to do this an assumption of local equilibrium is used once again. Using Equation 5-8 the equation is written:
which defines a new thermodynamic potential, the Gibbs free energy per unit mass, given by:
Changes in the Gibbs free energy per unit mass take the form:
which leads to the relations:
By differentiating each of the above equations a second time, it is possible to derive the following Maxwell relation:
(5-12)
It is now possible to derive an expression for the entropy of the solid. Assuming that the elastic stress is an invertible function of the strain, we can write s=s(σ,Ta). Thus:
Using the Maxwell relation in Equation 5-12:
Therefore:
By definition the heat capacity of the solid at constant stress is given by:
Thus:
(5-13)
Substituting Equation 5-13 into Equation 5-11 gives the following equation for thermoelasticity:
(5-14)
An additional heat source term is present in Equation 5-14, compared to the standard heat transfer equations in solids. This term couples the structural problem with the heat transfer problem. In turn the heat transfer equation couples back into the structural problem through the constitutive relationship.
In the particular case of a linear elastic material (in the absence of damping) the stress and strain are related by Duhamel–Hooke’s law:
(5-15)
where C is the elasticity tensor, σi is the initial stress, εi is the initial strain, and Tref is the reference temperature at which the strain and stresses take the initial values.
This equation couples the heat transfer equation to the structural problem. Given a temperature independent thermal expansivity, and no material damping, Equation 5-14 takes the form:
This is the usual form of the equation for linear thermoelasticity. In the COMSOL Multiphysics Material Library, a slightly unusual definition of the thermal expansivity is used, so that α is not temperature independent. COMSOL Multiphysics therefore implements a version of Equation 5-14, which automatically includes these effects. In brief, an expression for ε is obtained by solving Equation 5-15 and the expression is differentiated with respect to Ta while keeping σ constant.
Linearized Equations and Frequency Domain Formulation
In the time domain it is convenient to linearize Equation 5-14, and this is necessary to solve problems in the frequency domain.
The independent variables in the system are Ta and u (where u is the displacement vector). Since the displacements drive the thermal equations it is natural to write the displacement vector as consisting of small deviations from its equilibrium value such that:
where a is a constant that is subsequently set to 1 and that is used to keep track of the order of the term in the expansion.
The corresponding temperatures, stresses, and strains in the device are given by:
where a is a constant that is subsequently set to 1. a can be used to track the order of the deviation from T0.
Substituting these relations into Equation 5-14 and keeping terms up to the second order gives:
Note that the heat source term, Q, is assumed to be time independent (Q = Q0). The final term in the above equation results from consideration of the Taylor expansion of the strain at constant stress, which to first order can be represented by a Taylor series of the form:
Consequently:
Note that for the particular case of a linear elastic material:
where the material damping has been introduced into the equation system by means of an anisotropic loss factor, η.
Grouping together terms of the same order in a gives the following equations:
The linear equation can be written in the form:
(5-16)
Where a=1 and the subscript has been dropped from the dependent variable T (previously T1) and the stress deviation σ (previously σ1). Equation 5-16 is solved by the Thermoelasticity interface. Writing the equation in the frequency domain gives:
This is the frequency domain form of the equations solved by COMSOL Multiphysics.
Considering the quadratic terms in Equation 5-12, there are four second-order heat source terms in the equation system:
In the frequency domain, the first three terms involve the products of small harmonic deviations with other similar terms (that is, terms that vary in time in a similar manner to cos2ω t. The trigonometric identity cos2ω t=(1cos 2ω t)/2 makes it explicit that these terms have a time-independent heating effect in addition to the heating that occurs at double the frequency of the driving signal). These second-order terms have a small effect on the time-varying solution and lead to a second harmonic term in the frequency domain. However, if the time average is nonzero, this can, over many oscillatory periods, result in a change in the background temperature. In COMSOL Multiphysics the (nonzero) time-averaged heat sources are available as variables, so that the self-heating of the structure can be computed.
First note that the time average of the term QTED3 must be zero, since although the second-order stress term σ2 has a nonzero time average, its time derivative cannot. The other terms in the expression for QTED3 are all time independent. Each of the remaining terms can have a nonzero time average as these involve the multiplication of two separate first-order terms (potentially with phase differences between them). The following result for the multiplication of two general oscillatory complex numbers, a1ei(ωt+ϕ1) and, a2ei(ωt+ϕ2) are useful:
Taking the time average of the above expression gives:
The time average heat sources above can be written in the following form in the frequency domain:
The derivation of the last term <QTED2> strongly relies on the prerequisite that the temperature dependence of the thermal expansion coefficient alpha is thermodynamically consistent. Since this prerequisite is usually not satisfied due to measurement or formulation errors, we have decided to remove the variable from the physics interface.
COMSOL Multiphysics also defines a total heat source term (solid.Qh) of the form:
The realdot operator is used to evaluate these quantities in a manner that produces the correct terms in the Jacobian for sensitivity analyses. To compute the slow time dependence of the background temperature, the steady-state equation can be modified to allow for the time dependence of the background temperature and solved together with the frequency domain equations. In such a frequency-transient formulation, the following equations are solved:
Operators, Functions, and Constants and Studies and Solvers in the COMSOL Multiphysics Reference Manual