Stephen R Novascone

Benjamin W Spencer

Jason D Hales

Richard L Williamson

Evaluation of coupling approaches for thermomechanical simulations

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The processes involved in capturing energy from nuclear reactions and converting that to usable form can involve extreme thermal environments. To characterize the thermal and mechanical response of nuclear power plant components subjected to those conditions, one must consider the physics driving both the thermal and mechanical response, as well as the interactions between the two.

Numerical methods for the implicit solution of the partial differential equations that describe physical phenomena typically lead to the solution of a system of discretized equations. If multiple coupled physics are included in a model, the set of equations to be solved includes degrees of freedom from all of these physics. The strategies used to solve coupled sets of physics equations can be generally categorized as loose coupling and tight coupling. In loose coupling, the individual physics in a coupled problem are solved individually, keeping the solutions for the other physics fixed. After a solution is obtained for an individual physics, it is transferred to other physics that depend on it, and solutions are obtained for those physics.

These fixed-point iterations are repeated until convergence is obtained. If there is not a strong two-way feedback between the physics involved, convergence can be obtained quickly with a minimal number of loose-coupling iterations. An advantage of this approach is that it allows for independent codes to be coupled with relatively minor modifications to those codes, and they can each use their own solution strategies that are tailored for their solution domain. The disadvantage of loose coupling is that if there is strong two-way feedback between the physics, that approach can have an unacceptably slow convergence rate and may encounter convergence difficulty.

In tight coupling solution methods, a single system of equations is assembled and solved for the full set of coupled physics. The nonlinear iterations operate on the full system of equations simultaneously, taking into account the interactions between the equations for the coupled physics in each iteration. In cases where there is strong coupling between the physics, this approach can have faster convergence rates than loose coupling. The primary disadvantage of this approach is that it necessitates tighter coordination between the codes to solve the individual physics.

In thermomechanical problems, in the absence of evolving contact between components, the coupling between the heat conduction and solid mechanics equations is often primarily one-way. The temperatures obtained from the heat conduction equations cause thermal strains, which result in displacement of the mechanical model. These displacements typically have a negligible effect on the thermal model, and such problems can be readily solved using loose coupling strategies, or even by transferring data from a thermal code to a solid mechanics code and completely neglecting the effect of the mechanical solution on the thermal solution. The simulation of the response of a reactor pressure vessel to pressurized thermal shock conditions is a good example of such a problem in nuclear power generation. During an accident, the vessel could be subjected to rapidly decreasing temperature and pressure, potentially followed by a rapid repressurization. High tensile stresses can occur on the interior of the vessel due to the combined effects of thermal gradients and internal pressure. In this problem, temperature changes lead to thermally-induced strains, but the displacement caused by those thermal strains has a negligible effect on the thermal response.

Introducing evolving mechanical and thermal contact to thermomechanical problems transforms thermomechanical problems from being essentially one-way coupled problems to strongly twoway coupled problems. This is because the heat conductance across gaps between adjacent bodies is highly dependent on the distance between those bodies, which is a function of the mechanical deformation. A good example of this type of problem is the simulation of the performance of a light water reactor (LWR) fuel rod. Heat generated by fission is transferred through the fuel pellet, across the gap between the fuel and cladding, and through the cladding to the coolant. The conductance across the gap is strongly dependent on the composition of the gas in that gap and the size of that gap, which is driven by the mechanical response of the fuel system. The fission gas released from the fuel has a strong influence on the composition, and thus, the conductivity of the gas in the gap. The mechanical effects driving the evolving gap size include thermal expansion, swelling, densification and relocation of the fuel, and cladding creep. Because the gap conductance has a strong effect on the thermal response of the fuel system, the ability to efficiently and robustly solve the strongly coupled thermal and mechanical equations in the presence of evolving contact conditions is critical for a successful fuel performance modeling code.

This paper describes the solution environment used to enable tightly and loosely coupled simulations of thermomechanical problems, provides a review of the equations governing thermal and mechanical response, and demonstrates the performance of loose and tight coupling strategies on simple thermomechanical problems with varying degrees of feedback between the two systems. Following these simple demonstrations, the performance of these solution strategies is demonstrated on real-world nuclear engineering problems, first on a simulation of reactor pressure vessel response during pressurized thermal shock conditions and then on a fuel performance simulation. This work extends similar studies presented by the authors in Novascone et al. (2013) and Novascone et al. (2013).

Multiphysics solution environment

The work performed in this paper was done using codes built on the open source Multiphysics Object-Oriented Simulation Environment (MOOSE) [START_REF] Gaston | MOOSE: a parallel computational framework for coupled systems of nonlinear equations[END_REF], developed at Idaho National Laboratory (INL). MOOSE is a parallel computing environment for solving general systems of coupled partial differential equations based on the finite element method. MOOSE provides the framework for rapid development of physics simulation codes as well as access to solvers appropriate for nonlinear multiphysics problems.

To solve a nonlinear system of equations, it is common to begin with a residual statement [START_REF] Hales | Solving nonlinear solid mechanics problems with the Jacobian-free Newton Krylov method[END_REF]

r(x) = 0 (1)
where r is the residual with x as the unknown solution. The Jacobian is written as

J(x) = ∂r(x) ∂x . (2)
Newton's method is then

Compute J(x k ), r(x k ) (3) Solve J(x k )s = -r(x k ) for s (4)
x k+1 = x k + s (5)
which is continued until the update is sufficiently small or some other criterion is met. The Jacobian-free Newton Krylov (JFNK) method evaluates the action of the Jacobian through a finite difference approximation,

J(x k )v ≈ r(x k + v) -r(x k ) . (6)
This is an attractive form since neither the full Jacobian nor its element-by-element contributions are required. Despite not requiring the analytic Jacobian, the effect of the full Jacobian is seen from the first iteration of the iterative solver, unlike modified Newton or quasi-Newton algorithms. Thus, with only GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving linear systems[END_REF]) (which does not require J but Jv) and a function that computes the residual, JFNK finds solutions to nonlinear coupled equations with the convergence rate of a traditional Newton algorithm.

JFNK, like Newton's method, is a general technique for solving nonlinear equations. As such, it provides flexibility in selecting which phenomena are active in a given simulation as well as flexibility in choosing which model to activate for a given phenomenon. This generality is due to the fact that the approach is based on the evaluation of the residual, which is done according to whatever models and options are active for a given analysis.

Efficient solves using iterative methods require good preconditioners. The purpose of preconditioning is to decrease the condition number of the system being solved. In JFNK, it is common to use right preconditioning,

J(x k )M -1 (Ms) = -f (x k ) ( 7
)
where M is the preconditioner or preconditioning process. In this form, the solution approach involves two steps. First, solve J(x k )M -1 w = -f (x k ) for w. Then, compute s = M -1 w. Note that if M -1 = J -1 the iterative solve will converge in one iteration. However, computing J -1 is equivalent to solving the original system and so is not advantageous. It is necessary, therefore, to choose a preconditioner that reflects the character of J( x k ) in order to accelerate the iterative method but also one that is inexpensive to compute and apply. Simple approximations to sub-blocks of the Jacobian along the diagonal may suffice.

The software described here solves the system of equations using preconditioned JFNK using PETSc [START_REF] Balay | PETSc Users Manual[END_REF], which permits the use of a variety of preconditioning methods. In this environment, it is common to use algebraic multigrid preconditioner provided by the hypre [START_REF] Falgout | HYPRE: a library of high performance preconditioners[END_REF] project for large problems running in parallel. Computing factors with LU is common for single-processor cases or with the SuperLU [START_REF] Li | An overview of SuperLU: algorithms, implementation, and user interface[END_REF] distributed direct solver for smaller parallel problems (2-8 processors). The performance and scalability of the preconditioner varies depending on the problem being solved. For the solution of typical thermo-mechanical problems in the MOOSE environment, algebraic multigrid preconditioners often deliver the best performance and scalability. For small to moderately sized problems involving mechanical contact, direct solvers are often employed in this environment because of improved robustness in the presence of contact constraints. A review of the JFNK approach, and its application to solid mechanics, can be found in [START_REF] Hales | Solving nonlinear solid mechanics problems with the Jacobian-free Newton Krylov method[END_REF].

The equations to compute updates for multiphysics systems may be described as

J aa J ab J ba J bb s a s b = - r a r b (8)
where subscripts a and b represent different phenomena. In the tight coupling approach, JFNK solves for s a and s b simultaneously and includes the coupling terms in the matrix.

For loose coupling, s a and s b are computed separately. First s a is computed, with the values of s a held constant during the subsequent calculation of s b . Then, s a is again computed while holding s b constant. This pattern is repeated until the change in solution is acceptably small. This is referred to as fixed-point or Picard iteration. The FALCON fuel performance code uses a fixed number of nonlinear iterations to solve the thermal equations, a fixed number of iterations to solve the mechanical equations, and a fixed number of loose coupling iterations [START_REF] Lyon | Fuel Performance Analysis Capability in FALCON[END_REF].

In cases where the coupling terms J ab and J ba are significant, a tight coupling approach is advantageous. Given a sufficiently good initial guess, the quadratic convergence rate of JFNK will provide a distinct advantage over the linear convergence rate of the Picard approach [START_REF] Knoll | An efficient nonlinear solution method for non-equilibrium radiation diffusion[END_REF][START_REF] Lemieux | Implementation of the Jacobian-free Newton-Krylov method for solving the first-order ice sheet momentum balance[END_REF]. However, for a general nonlinear problem, a poor initial guess may lead to nonconvergence for either approach [START_REF] Mehl | Use of Picard and Newton iteration for solving nonlinear ground water flow equations[END_REF]. [START_REF] Danowski | A monolithic computational approach to thermo-structure interaction[END_REF] present results of loose and tight thermomechanical coupling applied to analysis of rocket nozzles as well as a discussion of possible convergence difficulty that may arise due to ranges of material properties, solution algorithm choices, and mesh density. For a review of acceleration techniques for the loosely coupled scheme, see [START_REF] Erbts | Acceleration of partitioned coupling schemes for problems of thermoelasticity[END_REF].

In the MOOSE environment, multiphysics models are typically solved by applying JFNK to solve Eq. ( 8) as a single tightly coupled set of equations. MOOSE has recently been extended to permit a physics simulation to run other nested physics simulations and pass a variety of types of data to and from those nested simulations. Multiple levels of recursive nesting can be used to permit simulation of phenomena occurring on widely disparate length scales, as demonstrated in [START_REF] Gaston | Physics-based multiscale coupling for full core nuclear reactor simulation[END_REF]. The mechanisms for running and transferring data to other simulations can also be used to partition a multiphysics simulation into separate models for the individual physics, transfer data between those models, and perform fixed-point iterations until all of the individual models are converged. The individual physics models are solved using preconditioned JFNK. This technique has been used in the work presented here to explore the relative performance of tight and loose coupling on thermomechanics problems.

To perform an objective comparison between loose and tight coupling strategies, it is important to use consistent convergence criteria for both approaches. In the work presented here, magnitudes of vector (l 2 ) norms of the residual vectors of the individual solution variables are used as convergence criteria. The model is considered to be converged when the norm of the residual is below a target quantity for each of the primary solution variables. The temperature and the x, y, and z components of the displacement vector are all treated as separate variables for this purpose. This ensures that the partitioned systems in the loose coupling strategy use the same convergence limits that are used for the tightly coupled solution. The fixed-point iterations in the loosely coupled solutions are considered to be converged when the initial residual for an individual physics solution, after applying the solution from the other physics, meets the same convergence criterion used for the nonlinear solutions of the individual physics models.

Governing equations for coupled thermal/mechanical simulations

The emphasis of this work is on methods for solving the equations of energy and momentum conservation. A brief overview of the governing equations for these physics is given here. The energy balance is given in terms of the heat conduction equation:

C p ∂T ∂t + ∇ • q -Q = 0, (9)
where T, and C p are the temperature, density and specific heat, respectively, and Q is the volumetric heating rate. The heat flux is given as

q = -k∇T, (10)
where k denotes the thermal conductivity of the material. Momentum conservation is prescribed assuming static equilibrium at each time increment using Cauchy's equation,

∇ • + f = 0 (11)
where is the Cauchy stress tensor and f is the body force per unit mass (e.g. gravity). The displacement vector u, which is the primary solution variable, is connected to the stress field via the strain, through kinematic and constitutive relations.

On interfaces between solid bodies, MOOSE enforces thermal and mechanical contact using a node on face, master-slave algorithm. Mechanical contact is enforced between interacting faces using a kinematic enforcement algorithm, which strictly enforces a non-penetration condition between slave nodes and master faces when the contact force between them is positive, and which releases the node otherwise. Alternatively, a penalty enforcement algorithm can optionally be used for mechanical contact. Details on enforcement of mechanical contact in the finite element method may be found in [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. For thermal contact, a geometric search is performed to find the face on the master surface onto which each node on the slave surface is projected. Heat transfer across the gap is based on the gap conductance, which can be calculated using a variety of models.

In the simplest model, the gap conductance, h gap , is inversely proportional to the gap distance, g:

h gap = k g g , (12)
where k g is the thermal conductivity of the material in the gap. The gap distance is taken as the greater of the actual gap and a userspecified minimum gap value, which prevents a singularity as the gap closes.

A series of thermomechanics problems of increasing complexity are provided in the following sections to evaluate the performance of tight and loose coupling strategies for solving thermomechanics problems.

Single-block thermomechanics example

Problem description

The first example studied here is a simple rectangular block considered under two different heating conditions. The domain is 5 units in length, has a 1 unit square cross section, and is meshed with fully integrated 8-noded hexahedral elements. The domain is assigned an elastic modulus of 1 × 10 6 , Poisson's ratio of 0.3, coefficient of thermal expansion of 1 × 10 -5 , and thermal conductivity of 1. The model is constrained such that only axial displacements are permitted. Note that the geometry and material properties are contrived and meant only to serve as a demonstration problem.

Two loading scenarios are considered. In the first, the block temperature is raised uniformly from 100 to 1100 over a single step.

In the second, one end of the block is held at a temperature of 100 while the opposite end is raised to 1100, again over a single step. Thus the second case provides a significant thermal gradient not present in the first case. Both scenarios were run using both loose and tight coupling, using a direct solver for application of the preconditioner and a single processor.

Results

The final temperature and axial displacements for the two loading cases are shown in Fig. 1. In each case, temperature and displacement results from the two coupling approaches were effectively identical. Note that in the first case, a uniform temperature field (Fig. 1(a)) results in a linearly varying displacement field along the axis (Fig. 1(b)) while in the second case, a linearly varying temperature field (Fig. 1(c)) results in a nonlinear displacement field (Fig. 1(d)).

Significant numerical solution differences were observed between the two cases, both in terms of iteration counts and run times. With tight coupling, the solution behavior was essentially independent of the load case, requiring four nonlinear iterations and having a run time of 1.26 s in both cases. With loose coupling, the uniform temperature case required two thermal and three mechanical nonlinear iterations, converged with a single fixedpoint iteration, and had a run time of 0.8 s. For the thermal gradient case, four fixed-point iterations were required, with a total of seven thermal and seven mechanical nonlinear iterations and a run time of 2.15 s.

Feedback between the thermal and mechanical solutions is clearly more significant in the second case, and the tightly coupled solution was significantly more efficient in that case. This increased feedback is due to the nonlinear displacement field induced by the change in the thermal gradient over a solution step. In the case with constant temperature, the displacement field does not affect the thermal solution. In the case with the temperature gradient, the change in the thermal gradient results in nonuniform incremental thermal strains across the beam. There are larger relative displacement increments between neighboring nodes in regions with larger incremental temperature changes. As a result, the nodes are redistributed, resulting in an incremental displacement field that varies nonlinearly as a function of the axial position. When those nonlinear incremental displacements are transferred from the mechanical solution, they render the prior thermal solution invalid, requiring multiple fixed-point iterations for convergence in a loose coupling solution strategy. This effect is obviously more pronounced with larger incremental thermal strains.

Thermomechanical contact example

Problem description

The second example studied here is a thermomechanics problem with thermal and mechanical contact between two domains. This problem was developed to compare the performance of the coupling approaches when applied to a difficult coupled thermal and mechanical contact problem with an evolving gap between two domains. Each domain is the same rectangular block used in the problem described in the prior section. One block is situated at a distance of 0.25 above and offset length-wise by 1.25 from the other block. The material properties are as described in the prior section, except that the lower half of the top block is assigned a coefficient of thermal expansion of 50 × 10 -6 . In other words, the lower half of the top block should have 5 times the thermal strain when compared to all other points in the domain given an equal change in temperature. These problems were run using one processor. Again, the unitless geometry and material properties are for demonstration purposes and are not representative of a real problem.

The upper block is cantilevered at the far end while the bottom surface of the lower block is fixed from translation in any direction. The simplified gap conductance model was used with a thermal conductivity of 4 × 10 -2 , and a minimum gap size of 5 × 10 -4 .

Each block has an initial temperature of 100. The temperature of the top surface of the upper block is held constant at 100 for the entire simulation. The temperature of the bottom surface of the lower block is increased from 100 to 1000 during the first 8 quasistatic load steps of this simulation, held at 1000 through load step 16 then decreased to 100 from load step 16 to load step 24, which is the final load step. All other surfaces are thermally insulated. The right side of the top block is initially held fixed, and is then monotonically displaced downward toward the bottom block (thereby closing the gap completely), from load step 8 to load step 16, at which point the displacement on that boundary is -0.26.

As the temperature of the lower block increases, heat is transferred to the upper block via the gap, which causes the upper block to thermally expand. The higher coefficient of thermal expansion in the lower half of the upper block causes the block to bend away from the lower block. After the lower block reaches a temperature of 1000, the upper block is lowered, so that the gap decreases, causing increased heat transfer from the lower block to the upper block, which results in more bending. After the upper block makes contact with the lower block, the temperature of the lower block is decreased, which reverses the bending of the upper beam, and thus decreases the gap size. As for the single-block problem, these cases were all run on a single processor using a direct solver for preconditioning.

Results

Temperature contour plots and the displaced mesh are shown in Fig. 2 at load steps 0, 8, 16, and 24.

The results for the primary solution variables (displacements and temperatures) for the tightly and loosely coupled simulations were identical within the convergence limits.

Fig. 3 shows a comparison of iteration counts for each step for the tight and loose coupling approaches for this model, as well as these same iteration counts plotted cumulatively. It shows the nonlinear iterations taken in the tight coupling approach, the number of fixed-point iterations in the loose coupling approach, and the total number of nonlinear iterations taken by the individual physics solutions in the loose coupling approach, accumulated for all fixed-point iterations in each step.

These plots of iteration counts show that the tightly coupled solution took a roughly five nonlinear iterations for each step, except at the steps when mechanical contact was first being established. The loose coupling case took about five fixed-point iterations per step early in the analysis, when the gap was large, but this number increased as the beams came into contact. The loose coupling approach overall took roughly twice as long as the tight coupling approach. This increased computational cost is largely due to the increased cost when the beams were in contact, as evident by the increased slope of the cumulative run time plot for the loosely coupled approach. The slope of the cumulative run time plot is more constant for the tight coupling approach.

Gap material conductivity sensitivity study

This problem clearly demonstrates that the behavior of a closing gap has a significant effect on solution performance. To better understand the effect of parameters controlling the behavior of the gap, two additional variants of this model were run with different values for the conductivity of the gap material, k g . Conductivities of 2 × 10 -2 and 1 × 10 -2 were used, which were 2× and 4× lower the conductivity of 4 × 10 -2 used in the baseline model. Comparisons of the iteration counts and run times (wall clock time) for these cases are shown in Fig. 4. These results demonstrate a moderate increase in iterations required for both solution approaches as the thermal conductivity of the gap is decreased. This behavior is expected because gap conductance is more sensitive to changes in gap size as the conductivity of the gap material is decreased. This effect is particularly pronounced in load step 16, when the two blocks come into contact.

Reactor pressure vessel simulation example

To ensure safe operation of LWR nuclear power plants, it is critical that the reactor pressure vessel (RPV) containing the reactor core maintain integrity during a wide variety of operational and postulated accident conditions. The regions of the RPV near the core are subject to significant neutron flux, which over time leads to embrittlement of the vessel. During an accident, the RPV could be subjected to a transient loading event known as a pressurized thermal shock (PTS), in which the interior of the vessel is rapidly cooled, followed by a rapid increase in pressure. The rapid cooling causes large temperature gradients through the wall of the vessel, and the accompanying thermal strains can lead to high tensile stresses on the interior of the vessel. The pressure increase can lead to a further elevation of the tensile stresses. In the unlikely event that a critical flaw is present in a vessel subjected to those conditions, the combination of tensile stresses, decreased temperature, and embrittled material could lead to unstable crack propagation through the wall of the vessel [START_REF] Odette | Embrittlement of nuclear reactor pressure vessels[END_REF]. A simulation of a RPV subjected to PTS conditions is used here as an example of an engineering thermomechanics application.

Grizzly component aging code

An analysis code called Grizzly is being developed at INL to address a variety of concerns related to aging of nuclear power plant systems, structures, and components. Grizzly is intended to model both aging processes in these components, as well as their response, in an aged condition, to the full set of loading conditions that they must safely withstand. The Grizzly application is built on the MOOSE framework, which provides a parallel, multiphysics equation solving environment, as well as standard physics models, including those for heat conduction and solid mechanics. Grizzly then builds on those base capabilities by adding applicationspecific models in a modular fashion.

The analysis of PTS in an embrittled RPV is an initial application targeted by Grizzly. The intent is to provide a capability similar to that provided by the FAVOR code [START_REF] Williams | Fracture Analysis of Vessels -Oak Ridge, FAVOR, v12.1[END_REF][START_REF] Dickson | Fracture Analysis of Vessels -Oak Ridge, FAVOR[END_REF], but in the framework of a code that can run 1D, 2D, or 3D models to enable the exploration of effects that can only be observed in 3D. To model the RPV fracture problem, Grizzly currently provides a general 2D axisymmetric or 3D thermomechanics capability to model the global response of the vessel to PTS conditions. It also provides a capability to transfer results to a 3D submodel of the local region in the vicinity of a postulated flaw, and provides fracture domain integral capabilities to evaluate mode-I stress intensity factors along the front of a meshed crack. The embrittling effects of neutron irradiation and thermal aging are accounted for by using the physically-based correlations in the EONY model [START_REF] Eason | A physically-based correlation of irradiation-induced transition temperature shifts for RPV steels[END_REF] to calculate the transition temperature shift.

RPV Simulation model

An analysis of a 3D model of a typical RPV for a pressurized water reactor (PWR) subjected to PTS loading conditions is used here for a comparative study of the performance of tight and loose coupling strategies for this type of problem. This model, shown in Fig. 6, uses two symmetry planes, with appropriate mechanical boundary conditions, to represent one fourth of the vessel. This model is representative of a typical PWR RPV, with an inner radius of the base metal in the beltline region of 2.034 m, a wall thickness of 0.199 m, and an overall height of 12.88 m. The inner surface of the vessel is lined with a 3.75 mm thick stainless steel cladding. The finite element model uses 8-node hexahedral elements and has a total of 659,133 elements and 734,771 nodes. In the beltline region of the vessel, there are 10 elements through the thickness: one to represent the liner and 9 for the base metal.

The base metal and stainless steel liner are represented using linearly elastic material models, with temperature-dependent coefficients of thermal expansion, elastic moduli, Poisson's ratios, and thermal conductivities. The coefficient of thermal expansion of the liner is roughly 20% higher than that of the base metal, which can cause significant thermally-induced stresses in the vicinity of the liner. For this demonstration, the PTS event is simulated by applying the spatially uniform coolant temperature and pressure histories shown in Fig. 5 to the entire inner surface of the vessel. The model starts out with a uniform temperature of 559.37 K throughout the vessel, to simulate normal operating conditions. A convective flux boundary condition is used with a time-varying convection coefficient to apply the thermal boundary conditions. These uniform boundary conditions result in a response that is essentially axisymmetric in the beltline region of the vessel, but this model can be used with nonuniform boundary conditions to simulate the effects of nonuniformities in both the boundary conditions and the vessel.

This model was run with uniform 100 s time steps through the loading history to obtain the coupled thermal and mechanical response to this event. The model was run using both a tightly coupled and a loosely coupled solution strategy. In both cases, it was run on 160 cores, spread across 10 nodes of a distributed memory cluster. For this model, a block-diagonal preconditioning matrix was used in conjunction with an algebraic multigrid method to apply the preconditioner. As for the other problems studied here, good convergence of the linear iterations was obtained, so the nonlinear iterations of the system were effectively the same as those that would be taken with a Newton method on either the full or partitioned system of equations.

RPV simulation results

Fig. 6 shows representative contour plots of the temperature and von Mises stress solutions obtained at the time of the rapid repressurization of the vessel. Fig. 7 shows time histories of the hoop stress at three locations at the beltline: an element on the outer diameter of the base metal, an element on the inner diameter of the base metal, and an element in the liner. The results for the tight and loose coupling cases overlay each other.

The thermomechanical response of an RPV under PTS conditions is a good example of a problem with relatively minimal two-way feedback between the thermal and mechanical response. Early on in the event, when the coolant temperature is rapidly decreasing, there are significant gradients in the incremental thermal response, which cause the incremental mechanical solution to affect the thermal solution, but as the temperature equilibrates in the vessel, this effect is diminished. for each step for the tight and loose coupling approaches for this model. The same metrics shown previously for the beam problem in Fig. 3 are used here.

From these iteration counts, it is evident that this model requires a very small number of fixed-point iterations when solved loosely coupled. Because of this, the loose coupling approach is very competitive, as can be seen in the comparative plot of the cumulative wall clock run times for the two approaches in Fig. 8. Early in the analysis, when there is more two-way feedback between the physics, the two approaches require almost identical run time, but later on, when only one fixed-point iteration per step is required, the loose coupling approach requires less time to solve per step, and in the end, the loose coupling approach requires less overall time to run this model.

LWR fuel performance simulation example

Thermomechanical simulations are an important tool used to understand the performance of nuclear fuel in the reactor. Nuclear fuel operates in an environment that induces complex coupling between physics, occurring over distances ranging from inter-atomic spacing to meters and time scales ranging from microseconds to years. Many important aspects of fuel response are inherently multidimensional. An important defining aspect of LWR fuel performance simulations is that thermomechanical contact between fuel and cladding plays a significant role in their response. For that reason, a fuel performance simulation is used in this study to compare the performance of solution strategies on an engineering application dominated by contact.

BISON fuel performance code

BISON is a fuel performance simulation code based on the MOOSE framework [START_REF] Williamson | Multidimensional multiphysics simulation of nuclear fuel behavior[END_REF]. It builds on the thermomechanics simulation capabilities provided by MOOSE and adds models specific to the materials and regimes experienced in nuclear fuel systems. Because it is based on a finite element-based framework, BISON can be used for modeling a wide variety of fuel types. More detailed information on the models in BISON and its application to a variety of fuel types can be found in [START_REF] Williamson | Multidimensional multiphysics simulation of nuclear fuel behavior[END_REF], [START_REF] Hales | Multidimensional multiphysics simulation of TRISO particle fuel[END_REF][START_REF] Hales | Verification of the BISON fuel performance code[END_REF]. BISON is used here to compare solution strategies for fuel performance simulations.

In BISON fuel simulations, a gap conductance model more complete than the one described in Eq. ( 12) is typically used:

h gap = h g + h s + h r (13)
where the gap conductance is calculated as a summation of components due to gas conductance, h g , solid-to-solid contact, h s , and radiative heat transfer, h r . The details of how these individual

LWR fuel performance simulation model

Tightly and loosely coupled strategies are evaluated using a 2D axisymmetric simulation of a short section of an LWR fuel pin. The details of this model are documented in [START_REF] Williamson | Multidimensional multiphysics simulation of nuclear fuel behavior[END_REF]. The assumed geometry, shown in Fig. 9, includes ten individual UO 2 pellets, Zr-4 cladding, an initial 80 m pellet-clad gap, and an open region to simulate the upper plenum. The clad outer surface temperature is held fixed at 605 K. The rod power is assumed to rise linearly over 3 h and then held constant for 8 × 10 7 s (about 30 months). An axially varying power profile is applied over the length of the fuel column. Although this variation is unrealistic for a short rodlet, the resulting axial variation in fuel temperature ensures that once pellet-clad contact begins, some portion of the rodlet is coming into contact over the remaining irradiation period. The finite element mesh for a single fuel pellet and neighboring clad is also shown in Fig. 9. This mesh has 2448 linear quadrilateral elements and 2898 nodes.

The thermal conductivity of the fuel is simulated using the temperature and burnup dependent Fink-Lucuta approach [START_REF] Fink | Thermophysical properties of uranium dioxide[END_REF][START_REF] Lucuta | A pragmatic approach to modelling thermal conductivity of irradiated UO2 fuel: review and recommendations[END_REF]. The Young's modulus, Poisson's ratio, and thermal expansion coefficient for the fuel were assumed constant at 219 GPa, 0.345, and 10 -6 (1/K) respectively, as given in [START_REF] Olander | Fundamental Aspects of Nuclear Reactor Fuel Elements[END_REF]. The thermal conductivity, Young's modulus, Poisson's ratio, and thermal expansion coefficient of Zr-4 were assumed constant at 16 W/m K, 75 GPa, 0.3, and 5 × 10 -6 (1/K) respectively, using typical values from MATPRO [START_REF] Allison | SCDAP/RELAP5/MOD3[END_REF]. Frictionless mechanical contact is enforced between the fuel and clad. A sophisticated model is used to calculate the fission gas released into the plenum [START_REF] Pastore | Physics-based modelling of fission gas swelling and release in UO2 applied to integral fuel rod analysis[END_REF], and the conductivity of the gas in the gap is a function of the composition of the gas. The rod is initially filled with 100% He, but as fission gas is released, the gas gap conductivity is reduced. The fuel performance simulations were run with 8 cores using a distributed direct solver to apply the preconditioning.

LWR Fuel performance simulation results

The fuel performance simulation results obtained from the tightly and loosely coupled solution strategies are effectively identical. Fig. 10 shows representative temperature contour plots at the final time step. The axial variation in the temperature is reflective of the nonuniform axial profile applied to this model. At this point in time, sufficient fuel swelling and cladding creep have occurred to completely close the gap between the fuel and cladding.

Fig. 11 shows a comparison of iteration counts for each step for the tight and loose coupling approaches. Early in the analysis, the loose coupling approach required a small number of fixed-point iterations, but as the gap between the fuel and cladding closed, the number of fixed-point iterations required per step increases considerably. The tightly coupled solution, on the other hand, required a fairly constant number of nonlinear iterations per step. Later in the analysis, after the gap completely closed, the number of fixedpoint iterations decreased for the loose coupling case. The overall run time for the loose coupling case was significantly higher than for the tight coupling case.

Fill gas composition sensitivity study

The presence of fission gas in the plenum has a significant effect on the behavior of LWR fuel rods because of the corresponding decrease in the gas conductivity. To investigate the influence this has on the efficiency of the solution strategy, the composition of the initial fill gas in the LWR fuel rod model was modified to include a mixture of He and Xe. Two additional variants of the model were run: one where the gas was composed of 25% Xe, and one where it was composed of 50% Xe. The remainder of the gas in both cases is pure He. These results are compared with the baseline case, which had 100% He.

As would be expected, increasing the amount of Xe in the fill gas results in increased temperatures across the fuel. Fig. 12 shows time histories of the fuel temperature at the fuel centerline and fuel surface for these three cases. These plots show only the results obtained from the tightly coupled solution strategy because they are exactly overlaid by the results from the loose coupling strategy.

Fig. 13 shows a comparison of the performance of the loose and tight coupling solution strategies for the fuel model with varying fill gas. The tightly coupled solution required roughly the same number of nonlinear iterations per step regardless of the fill gas mix, but as the Xe fraction in the gas increased, the solution with loose coupling required significantly more fixed-point iterations, which resulted in a significantly increased run time. It is interesting to note that this increase in solution effort occurred over the period when the gap was closing, as evidenced by the increase in the slopes of the cumulative iteration count and solution time plots. Once the gap fully closed, the slopes for the various models once again assumed relatively similar values.

Summary of performance of loose and tight coupling

A summary of run times for all of the problems presented in this paper is given in Table 1. The table also gives the ratio of loosely coupled run time to tightly coupled run time for each problem. A ratio of less than one (loosely coupled required a shorter wall clock time than tightly coupled) is seen for two problems: the one block problem with constant temperature of Section 4 and the reactor pressure vessel problem from Section 6. With the exception of the case with a large thermal gradient, neither of those problems exhibited a strong two-way feedback, with the result that the loosely coupled approach was more efficient. For the two block problem of Section 5 and the LWR fuel rod problem of Section 7, the ratio of loosely coupled wall clock time to tightly coupled time increases with decreasing gap conductivity, particularly for the LWR fuel rod problem. With a lower gap conductivity, the thermal solution is more sensitive to changing gap distance and therefore has stronger two-way feedback. The tightly coupled approach is expected to be increasingly more efficient in this scenario, and the data in Table 1 is in agreement with this anticipated trend.

Conclusion

To investigate the effectiveness of loose and tight coupling strategies for solving thermomechanics problems, a series of simulations has been conducted on a set of models representative of behavioral regimes of engineering interest. A flexible computational framework that permits running simulations in both a tightly and loosely coupled fashion has been used to run these models both ways. Performance metrics on these problems have been gathered to permit a quantitative comparison of tight and loose coupling strategies. The following observations can been made from this study:

• A loose coupling strategy that uses a series of fixed-point iterations to repeatedly solve the individual physics models and transfer results between them generates results that are effectively identical to those obtained from a tightly coupled solution of the full coupled system of equations. This has been demonstrated on a variety of thermomechanical problems. For an objective comparison of the two approaches, care must be taken to use consistent convergence criteria. • For problems that involve essentially one-way feedback between the thermal and mechanical response, a small number of fixedpoint iterations are required for convergence in the loose coupling approach, and loose coupling gives better overall performance than tight coupling. In the models of this class studied here, overall run times for tight coupling were about 30% to 60% longer than those for loose coupling. This stands to reason, as it is more efficient to solve two smaller systems of equations than one large system of equations. • As the amount of two-way feedback between the thermal and mechanical response increases, loose coupling requires an increasing number of fixed-point iterations, while the number of nonlinear iterations for a tightly coupled solution of the full system remains relatively constant. In the models studied here with strong two-way feedback, the run times for loose coupling ranged from about 70% to 275% longer than those for tight coupling. • Two phenomena observed here that can lead to strong twoway feedback between the thermal and mechanical solutions are incremental thermal gradients and thermal contact. In both cases, a change in the configuration due to the mechanical response leads to a change in the thermal solution. This effect is typically much more pronounced for thermal contact than for incremental thermal gradients. • For the nuclear engineering applications studied here, the RPV model is effectively solved using either approach, while fuel performance simulations appear to be much better solved using a tightly coupled approach.

The models considered in this study are a mix of relatively small problems run on a single processor and larger problems run using many processors. In addition, multiple methods were used for preconditioning, with an algebraic multigrid preconditioner employed for the RPV problem and a direct solver for the remaining problems. The preconditioning strategy and parallelism have a strong effect on the run time due to the overhead costs and scalability of the preconditioner being employed. However, within the JFNK solver, as long as sufficiently tight convergence is obtained in the linear iterations, the preconditioning and parallelism do not affect the convergence of the Newton iterations or the Picard iterations (for the loosely coupled cases).

For the problems considered in this study, the relative solution times between loosely and tightly coupled solution strategies for a given problem type appear to be independent of preconditioning and parallelism. Because the tightly coupled approach requires preconditioning on a larger equation system, preconditioning that system generally requires more computational effort than the combined effort required to precondition the smaller single-physics equation systems used in the loosely coupled approach. If the preconditioning method being used scales well for problem size being considered, the relative solution costs of loose and tight coupling should be fairly independent of the preconditioner and number of processors. A study to quantify the effects of preconditioning and parallelism on the relative performance of loose and tight coupling strategies would be a logical follow-on to this work. Scalable parallel performance is becoming increasingly important for practical engineering applications as these problems become increasingly large (3D) and complex (multiphysics).

It is important to note that there are other factors besides solution efficiency that could motivate the decision to use tight or loose coupling. For instance, loose coupling permits more flexibility in model definition, because dissimilar meshes can be used for the individual physics and preconditioning could be tailored for individual physics, which could lead to improved convergence robustness. Conversely, model setup is much more convenient for tightly coupled solutions in the multiphysics solution framework used in this study. This is especially true if a larger set of coupled physics equations is to be solved.

Fig. 1 .

 1 Fig. 1. Temperature and displacement results for the two thermal load cases considered. (a) and (b) are for a uniform temperature loading and (c) and (d) are for the case having an axial temperature gradient.

Fig. 2 .

 2 Fig. 2. Thermomechanics simulation with gap at several steps.

Fig. 3 .

 3 Fig. 3. Comparison of solver performance for loose and tight coupling approaches for beam model. (top) Iterations per step. (middle) Cumulative iterations up to current step. (bottom) Cumulative wall clock run time.

Fig. 4 .

 4 Fig. 4. Comparison of solver performance for loose and tight coupling approaches for beam model gap material conductivity sensitivity study. (top) Iterations per step, showing fixed-point iterations for loose coupling and nonlinear iterations for tight coupling. (middle) Cumulative iterations up to current step. (bottom) Cumulative wall clock run time.

Fig. 5 .

 5 Fig. 5. Time histories of coolant pressure and temperature uniformly applied to interior surface of RPV to simulate PTS event.

Fig. 6 .

 6 Fig. 6. Solution results for RPV model showing (a) temperature contours and (b) von Mises stress contours during pressure spike at time 4800 s. Results for loosely and tightly coupled cases are the same.

FigFig. 7 .Fig. 8 .

 78 Fig. 7.Hoop stress time histories for tight and loose coupling runs at three locations at the beltline: an element on the outer diameter of the base metal, an element on the inner diameter of the base metal, and an element in the liner.

Fig. 9 .

 9 Fig. 9. Geometry, materials and mesh used to simulate an axisymmetric discretepellet fuel rodlet.

Fig. 10 .

 10 Fig. 10. Contour plot of temperature at the end of simulation. (a) Tightly coupled and (b) loosely coupled. Coordinates are scaled by 0.5 in the axial direction.

Fig. 11 .

 11 Fig. 11. Comparison of solver performance for loose and tight coupling approaches for fuel model. (top) Iterations per step. (middle) Cumulative iterations up to current step. (bottom) Cumulative wall clock run time.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Temperature histories with varying percentages of Xe in the initial fill gas.

Table 1

 1 Wall clock time comparison.

		Loose (s)	Tight (s)	Ratio
	One block constant temp	0.800	1.26	0.635
	One block temp gradient	2.15	1.26	1.71
	Two block kg=4 × 10 -2	1627	780	2.086
	Two block kg=2 × 10 -2	1555	726	2.14
	Two block kg=1 × 10 -2	1452	647	2.24
	Reactor pressure vessel	3.31 × 10 4	4.23 × 10 4	0.783
	LWR Fuel He gap	4.22 × 10 3	2.50 × 10 3	1.69
	LWR Fuel 25%Xe gap	5.70 × 10 3	2.56 × 10 3	2.23
	LWR Fuel 50%Xe gap	1.23 × 10 4	3.28 × 10 3	3.75

Acknowledgements

The authors thank Derek Gaston, David Andrs, and Cody Permann from INL for support developing coupling methods in MOOSE, and Marie Backman from the University of Tennessee and William Hoffman from the University of Idaho for support developing the RPV model. The submitted manuscript has been authored by a contractor of the U.S. Government under Contract DE-AC07-05ID14517. Accordingly, the U.S. Government retains a non-exclusive, royalty free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes.


Amy Duwel

John Gorman

Marcie Weinstein

Jeff Borenstein

Paul Ward

Experimental study of thermoelastic damping in MEMS gyros

Keywords: Quality factor, Damping, Q, Thermoelastic

We present new experimental data illustrating the importance of thermoelastic damping (TED) in MEMS resonant sensors. MEMS gyroscopes have been used to demonstrate that both the choice of materials and variations in device design can lead to signi®cant differences in the measured quality (Q) factors of the device. These differences in the Q-factor can be explained by including the contribution of thermoelastic damping, which varies strongly between the different silicon etch-stop compositions used in this study. Known damping mechanisms, such as ¯uid damping, anchor damping, and electronics damping are minimized and held ®xed in this experiment so that materials effects can be isolated.

Introduction

Thermoelastic damping (TED) has been identi®ed as an important loss mechanism in MEMS resonators [1±3]. The impact of device geometry on the level of damping has been considered both experimentally and theoretically [3±5]. Although most work till date has focused on analytically tractable beam systems, the need for high Q resonators extends to a wide variety of applications, including accelerometers [START_REF] Gibbons | The design, fabrication, and testing of a micromechanical silicon oscillating accelerometer[END_REF], Coriolis rate sensors [START_REF] Kourepenis | Performance of MEMS inertial sensors[END_REF], chemical sensors [START_REF] Weinberg | Modeling flexural plate wave devices[END_REF], and RF ®lters [START_REF] Wang | VHF free±free beam high-Q micromechanical resonators[END_REF]. In most applications, maintaining a high quality factor results in reduced readout errors, lower power requirements, improved stability, and increased sensitivity.

The resonator devices used in this study were MEMS gyros fabricated at Draper Laboratory. The excitation and readout mechanism is capacitive, and the measurements are made in a vacuum better than 1 mTorr. We ®nd that the simple calculation by Zener [1,2] for thermoelastic damping in ¯exural mode beam resonators work well to describe the effects of beam width and material properties on the Q-factor of the resonant-mode gyros. In this paper, the effects of alloying on thermal conductivity and the resulting effect on Q-factor are explored. We compare the measurements on devices made from silicon±germanium alloys with Ge concentrations of up to 30%, and with boron doping of up to 2 Â 10 20 cm À3 . Using thermal conductivity data from Dismukes et al. [START_REF] Dismukes | [END_REF] on Boron doped SiGe materials, we show that the differences in thermal conductivity can account for our observed variations in Q factors.

In Section 2, we introduce the MEMS resonator being studied and the experimental setup. We also discuss the impact of our measurement setup on the measured Q of the resonator. In Section 3, we review the thermoelastic damping and discuss some of the important material parameters that must be known to quantify this effect. Section 4 compares a theory for total Q-factor to measured data.

Experimental setup

The MEMS gyro is a tuning fork resonator which senses the angular rate. A photograph of the gyro (top view) is shown in Fig. 1. The masses are driven in a tuning fork resonance mode in the plane of the wafer. In response to a ẑ-directed angular rate, the two proof masses move out of the wafer plane (AEŷ) in opposite directions. This motion is sensed capacitively and the amplitude of this motion is the desired signal proportional to angular rate. The resonating proof mass always serves as one plate of the sense capacitor, while the other plate is ®xed to the substrate. The gyro has two possible readoutsÐone ®xed capacitor plate is beneath the proof mass to sense y-motion, and the other senses the tuning fork motion in the x-direction. In this experiment, we focus on measurements of the tuning fork resonance in the x-direction for comparing the Q values of devices made from various materials.

In the design of resonant sensors, the readout circuit must balance the need for low insertion loss with high Q-factor. Any output signal current, for example, results from a coupling to the resonating sensor and can lead to dissipation of the mechanical resonance [START_REF] Varghese | Resistive Damping of pulsesensed capacitive position sensors[END_REF]. In the presence of a voltage bias, the current is de®ned by I _ CV _ VC, and for the case of a dc bias is simply proportional to the proof mass motion through _ C _ u dC=du, where u is the amplitude of motion in the x-direction. This current is ampli®ed and converted to a voltage, and the transfer function of the output to the input voltage is used to measure Q. Fig. 2 shows the circuit schematic used to drive the sensor and measure the Q-value. For Q measurements, the input drive voltage is the sum:

V drive V dc V in sin ot;
(1) so that the force on the sensor is:

F drive V 2 dc V dc V in sin ot V 2 in sin 2 ot (2)
The transfer function is taken by comparing V out at o to V in . We use the half power bandwidth relative to the resonant frequency to measure the Q. This works well as long as the drive coupling to other modes is minimized and the spectrum looks like that of an isolated resonator. When using this method, it is also important to drive the resonator at low amplitudes, so that nonlinear effects are negligible. Bandwidth and ringdown testing1 typically agreed to better than 10% even without taking the different bias conditions and corresponding electrostatic spring effects into account. The resonance peaks were very symmetric, and the standard deviation of the repeated bandwidth measurements is 1%.

All measurements were conducted at the wafer level under a vacuum of approximately 1 mTorr. Fig. 3 shows measurements of Q-factor as the pressure in the wafer probe station is reduced. In the range of 1 mTorr, the contribution of ¯uidic damping is reduced so that measurements of Q are found to be independent of small changes in vacuum level. We expect any residual ¯uidic damping to be constant for all measured devices due to similar device geometry.

The measurement itself can in¯uence the quality factor of the gyro. In order to produce a signal out, power has to be coupled from the resonator, and this power can be dissipated. We include in the Appendix A, a simple model for the gyro dynamics and the role of the applied voltages. We use this model to calculate the electronics contribution to the Q-factor and include it in the total calculated Q. Since all of our gyro designs use almost identical input and output coupling capacitors, the contribution of the electronics damping does not vary signi®cantly from device to device.

Thermoelastic damping

Thermoelastic damping arises from the coupling of the stress±strain state equation to heat ¯ow in the material. The familiar effects of the coef®cient of thermal expansion (a), where stresses (s) and strains (e) can be thermally induced, can be derived from the dependence of the Hemholtz free energy on both strain and temperature. From this starting point, however, one also ®nds that the temperature gradient in a material can be driven by changes in the strain. The coupled (linearized) state equations for an isotropic medium are [START_REF] Lifshitz | Theory of Elasticity[END_REF][START_REF] Nowacki | Thermoelasticity[END_REF]:

e s E aT (3)
dT dt k C V r 2 T À EaT 0 C V de dt ( 4
)
where E is the Young's modulus (Poisson's effect is neglected here), T 0 the nominal equilibrium temperature, and T the temperature offset from equilibrium. The coef®cient of thermal expansion is denoted by a, the thermal conductivity by k, and the heat capacity at constant volume by C V .

An effective transfer function for stress versus strain can be found by applying Fourier transform techniques to (3) and [START_REF] Lifshitz | [END_REF]. By assuming sinusoidal steady-state solutions with periodic or discrete spatial dependence [et; x D ẽs; k], temperature can be eliminated in (3). This yields a complex, frequency dependent Young's modulus that includes both mechanical and thermal material properties [START_REF] Lakes | Viscoelastic Solids[END_REF]. This effective modulus exhibits a zero followed closely by a pole in the frequency domain, which results in a resonance (called a Debye peak) in the imaginary component. With a complex modulus, the strain in the material can be out of phase with the applied sinusoidal stress, and the integral of work done (energy lost per cycle) will be a function of the phase shift.

The ratio of the imaginary component to the real component of the modulus de®nes a frequency dependant quality factor for the material. This approach reproduces with the quality factor derived by Zener [1,2] for a ¯exural mode beam resonator:

Q À1 TED Ea 2 T 0 C V ot 1 ot 2 (5) with t C V b 2 kp 2 (6)
where o is the mechanical resonance frequency of the beam and b represents the width of the resonator in the plane of the resonance motion. In deriving expression ( 5) from ( 3) and (4) using Fourier transforms, a spatial dependence, such as k p=b (the ®rst eigen value of the uncoupled thermal equation of a ¯exural mode beam) must be applied, which is analogous to the type of approximation used by Zener [1,2]. While this approach does not offer a rigorous solution to the problem (which can be found in [START_REF] Lifshitz | [END_REF]), it does connect thermoelastic effects with standard inelastic damping theory and offers intuition for the loss mechanism. Different material relaxation mechanisms, such as point defect motion, dislocation movement, and thermoelastic damping will each have an associated relaxation time constant, as in Eq. ( 6). The total amount of damping will depend on the proximity of the device operating frequency to the Debye peak frequencies.

In the resonators tested, several different materials were used. We especially focused on Boron-doped SiGe epitaxial materials since the presence of Ge reduces strain in the epi layer and the resulting devices can be extremely ¯at and well-machined. In addition, the use of Ge as an etch stop eliminates the need for EDP in processing [START_REF] Borenstein | A new ultra-hard etch-stop layer for high precision micromachining[END_REF], which is a huge advantage. However, considering the issue of thermoelastic damping, the use of an alloy can signi®cantly reduce the thermal conductivity due to phonon scattering. Based on the data of Dismukes [START_REF] Dismukes | [END_REF], we summarize in Table 1 the The notation ``Int'' refers to a linear interpolation between the silicon and germanium values, based on the alloy composition.

material parameters used for the calculation of Q TED . We also measured devices made from Boron diffused silicon (noted as B-diff in the plots), where the Boron concentration was approximately 10 20 . Although the Debye resonance frequency (1/t in Table 1) is order of magnitude above our operating frequency of about 12 kHz, the importance of the thermal conductivity is apparent when we compare the theory to our measurements.

Data versus theory

The materials tested, as summarized in Table 1, are made from: (1) silicon boron epi, (2) SiGeB epi, with 2 at.% germanium, (3) SiGeB epi, with 23 at.% germanium, (4) SiGeB epi, with 30 at.% germanium, and (5) boron diffused silicon. We measured 9±10 devices from each wafer, and took the average Q-value. The standard deviation of the Q measurements across a wafer was typically about 12%. We compared the measured (average) Q values for each material to a calculation based on the following formula:

1 Q total 1 Q electronics 1 Q TED 1 Q other ( 7
)
While Q TED and Q electronics can be calculated for each material, the Q other represents as yet unmodeled effects such as anchor damping, any residual gas damping, and other material losses. We reason that Q other should be the same for all devices and materials if the calculated Q TED accounts for all differences between the materials. Thus, we ®nd a single value of Q other that minimizes the error between our set of measurements and the theory expressed by Eq. [START_REF] Kourepenis | Performance of MEMS inertial sensors[END_REF]. The extent to which the data and theory do not match can be helpful in identifying additional material damping mechanisms. For the particular device design used in this study, Q other 250,000 worked well for all the materials measured.

The electronics damping, calculated in the Appendix A, was in the range of Q electronics 3:5 Â 10 11 but varies slightly between materials since the undercut (which we measure to be typically less than 0.5 mm) can affect the readout capacitor size. The thermoelastic Q-value for each material is listed in Table 1. The plot in Fig. 4 compares the data to the theory.

In addition to the data presented in Fig. 4, we also compared the measured Q values of various device designs. In several iterations of the tuning fork gyro, the beam widths were varied to optimize the modes of the devices. Based only on the beam width and resonance frequency data, we have been able to use thermoelastic damping to explain why some designs exhibited much higher Q values than others. The match of thermoelastic damping theory to our data over a range of device designs indicates that the anchor damping does not vary much between these iterations. However, the relatively low best ®t values of Q other 250,000 suggest that anchor damping may be an important limiting factor in all designs.

Conclusions

We have measured the quality factors of resonators made from several different semiconductor materials. Our data can be explained well by thermoelastic damping and electronics damping calculations. Changes in Q TED can affect our ®nal measured Q values by factors of two and three. Because the MEMS devices studied operate at frequencies well below the thermoelastic damping resonance, the value of Q TED can be optimized by moving the thermoelastic Debye peak to higher frequencies. In the gyros studied, this corresponds to choosing materials with higher thermal conductivity and designing gyros with narrower ¯exures. In contrast, for MEMS devices, such as RF resonators and ®lters, the optimization of Q TED may require materials and designs that lower the thermoelastic Debye resonance, to move it further from the frequency of operation.

The simple ¯exural beam models provided by Zener [1,2], and Lifshitz and Roukes [START_REF] Lifshitz | [END_REF] were surprisingly successful in estimating thermoelastic Q values and in identifying the material and design parameters that optimize the device Q. More precise optimization of complex geometries can be achieved using numerical approaches to solving the coupled thermoelastic equations. Future work in modeling additional loss mechanisms, such as anchor damping will further improve MEMS resonator designs.

Fig. 1 .

 1 Fig. 1. Photograph of a Draper tuning fork gyro.

Fig. 2 .

 2 Fig. 2. Circuit schematic of capacitive readout for MEMS resonator. The mode of interest results in a movement of the gyro mass in the x±z plane. The gyro resonance is excited on one side of the mass by an electrostatic force and the change in the readout capacitance is sensed. Comb finger capacitors are used so that the readout is linear in x.

Fig. 3 .

 3 Fig. 3. Pressure versus measured Q-value of a MEMS gyro. The measurements are taken to verify that the influence of gas damping is minimized at 1 mTorr.

Fig. 4 .

 4 Fig. 4. Data versus calculated Q values. In the theory, Q other 250,000 is used as a fitting parameter.

Table 1

 1 Material and device parameters used in thermoelastic damping calculations

	Quantity	Symbol	Units	SiB epi	SiGeB 2%	SiGeB 23% SiGeB 30% B-diff	Si value	Ge value
	Thermal coefficient	a	1/8C	2.57E À 06 2.65E À 06 3.33E À 06 3.55E À 06 2.59E À 06 2.59E À 06 5.80E À 06
	of expansion
	Reference			Si value	Int	Int	Int	[20]	[17,18]	[18]
	Modulus	E	N/m 2	1.69E 11 1.68E 11 1.54E 11 1.49E 11 169E 11 1.69E 11 1.03E 11
	Reference			Si value	Int	Int	Int	Si value	[18]	[18]
	Density	r	kg/m 3	2.33E 03 2.39E 03 3.02E 03 3.24E 03 2.33E 03 2.33E 03 5.35E 03
	Reference			Si value	Int	Int	Int	Si value	[17]	[18]
	Specific heat	C sp C V /r J/(K kg)	7.00E 02 7.05E 02 6.23E 02 5.95E 02 7.00E 02 7.13E 02 3.20E 02
	Reference			Si value	Int	Int	Int	Si value	[19]
	Thermal conductivity	K	J/(K s m) 80	10.67	5.71	5.52	80.00	1.56E 02 58.61
	Reference			[16,21]	[16,21]	[16,21]	[16,21]	[16,21]	[18,19]	[18]
	Tire constant	t	s	9.00E À 08 7.62E À 07 1.40E À 06 1.54E À 06 7.44E À 08
	Beam width	b	m	6.60E À 06 6.90E À 06 6.48E À 06 6.00E À 06 6.00E À 06
	Resonant frequency	f o/2p	Hz	12,851	13,327	10,357	10,785	12,361
	Thermoelastic Q-value	Q TED		6.70E 05 7.48E 04 4.07E 04 3.30E 04 8.29E 05

Note that both bandwidth and ringdown tests for measuring Q will be affected by the presence of nearby modes. According to theory (Nayfeh), ringdown tests should not be affected by amplitude nonlinearities. For small drive amplitudes, the effects of amplitude nonlinearities on the resonance curve measurements can also be minimized.

Acknowledgements

The authors wish to thank Professor Gene Fitzgerald and Mathew Curie for SiGe ®lm growth, and thanks to Draper laboratory technicians for fabrication assistance. Thanks to Robert White for his initial measurements of Q factors in gyros and to Bill Sawyer for assistance with the fabrication. We especially thank Mark Mescher for valuable conversations, as well as Marc Weinberg, Chris Dube  and Richard Elliott for valuable advice and support. Finally, we thank Neil Barbour, Jim Connelly, Tony Kourepenis, and Jim Sitomer for their support of this work.

Appendix A. Calculation of electronics damping

The tuning fork resonance is excited by applying an ac voltage to the outer capacitor (left capacitor in Fig. 2). The electrostatic force on the gyro is proportional to the square of the drive voltage. At the output, the voltage bias needed to generate a readout signal can also place a force on the sensor. To model the system dynamics, we balance these electrostatic forces on the sensor with the inertia and resorting force on the masses. The force balance becomes simplest when the motion of the resonator is already decomposed into linear eigen modes. For the tuning fork mode of interest, the masses move in the wafer plane, in the x-direction of Fig. 2. The modal amplitude is represented by u(t), and the force balance becomes

The electrostatic forces in [START_REF] Gibbons | The design, fabrication, and testing of a micromechanical silicon oscillating accelerometer[END_REF] are proportional to the change in capacitance with sensor mode amplitude. For these gyros, the drive and sense capacitors on the motor axis are comb designs, and to ®rst order the term dC/du is a constant that depends on geometry but is independent of u. In Eq. ( 6), m is the modal mass and b is the modal damping from sources other than the electronics.

If the ampli®er gain is in®nite, the node voltage V n is equal to V B . The gain is assumed to be purely imaginary at the operating frequency. For ®nite gain, we write

I s is the sensor current, R in the input resistance of the op-amp, and A the ampli®er gain. These equations combine to give

For large ampli®er gain, R in is small and a Taylor expansion can be applied to [START_REF] Dismukes | [END_REF]. In the limit of large ampli®er gain, Eq. ( 6) becomes

Thus, the applied bias results in a term proportional to velocity. By taking the ratio of stored energy to the energy dissipated per cycle, the contribution of the electronics to the Q of the system is

Note that designing for large readout gains (high V B and dC/du) can degrade the resonator Q-value. However, the readout ampli®er gain can compensate for this effect.