AIAA-2003-1119 Implicit Approaches for Moving Boundaries in a 3-D Cartesian Method Scott M. Murman* Michael J. Aftosmis t ELORET NASA Ames Research Center Moffett Field, CA 94035 MS T27B Moffett Field, CA 94035 smurman_nas.nasa.gov Marsha J. Berger* Courant Institute 251 Mercer St. New York, NY 10012 1 Introduction This work considers numerical simulation of three-dimensional flows with time- evolving boundaries. Such problems pose a variety of challenges for numerical schemes, and have received a substantial amount of attention in the recent literature. Since such simulations are unsteady, time-accurate solution of the governing equations is required. In special cases, the body motion can be treated by a uniform rigid mo- tion of the computational domain. For the more general situation of relative-body motion, however, this simplification is unavailable and tile simulations require a mech- anism for ensuring that the mesh evolves with the moving boundaries. This involves a "remeshing" of the computational domain (either locali:,ed or global) at each physical timestep, and places a premium on both the speed and robustness of the remeshing algorithms. This work presents a method which includes unsteady flow simulation, rigid domain motion, and relative body motion using a time-evolving Cartesian grid system in three dimensions. *Member AIAA tSenior Member AIAA Copyright @2003 by the American Institute of Aeronautics and Astronautics, Inc. No copyright is asserted in the United States under Title 17, U. S. Code. The U. S. Government has a royalty-free license to exercise all rights under the copyright claimed herein for Governmental purposes. All other rights are reserved by the copyright owner. While 3-D moving-boundary simulations have been performed on body-fitted structured and unstructured grid systems for some timer1 4], the literature on Carte- sian methods has been largely restricted to two dimensions and/or simplified configu- rations. Early approaches to the 3-D Cartesian moving-body problem included both the Volume of Fluid (VOF)[5] and level-set approaehe,'_[6]. NIore recently, there has been interest in the immersed-boundary[7-9] and cut-cell Cartesian approaches[10- 13]. Non-body-fitted, Cartesian approaches like these are particularly interesting since they can be made both extremely fast and robust[14]. Moreover, they are comparatively insensitive to the complexity of the input geometry since the surface description is decoupled from the volmne mesh, and can therefore handle complex geometries with relative ease. This work builds upon the inviscid, Cartesian cut-cell solver in Ref. [151. In this method, the cells which are cut by the boundary geometry can be arbitrarily small, making explicit update schemes overly restrictive for time-dependent problems. Ap- proaches to overcome this restriction usually either extend the difference stencil of the spatial terms[16], or use a cell-merging approach [10], so that cut-cells carl be advanced with the explicit timestep of a full urtcut cell. Both approaches have been successful for two-dimensional simulations with movinf; boundaries[10, 13, 161, how- ever, the coupling of cell size and allowable timestep with explicit methods implies that tile boundary motion will be restricted by the size of the finest boundary inter- secting Cartesian cell during a single timestep. In a Cartesian moving-boundary scheme, the most delicate operation is the re- intersection of the body with Cartesian grid at each time step. Not only is this the most computationally expensive part of Cartesian me'_h generation, but it requires special procedures to ensure that the floating-point intersection calculations are al- ways robust[14]. In a two dimensional mesh with O(N 2) ceils, the boundary geometry only intersects O(N) cells. In three dimensions, however, the boundary is a surface instead of a line and the number of intersection calculat!ons is squared with respect to the 2-D case. Moreover, these intersection calculation,'_ are higher-dimensional, and therefore each is more expensive and more difficult to compute robustly. Arguments for both efficiency and robustness in three dimensions weigh heavily in favor of moving the body as rarely as possible to minimize the re-cutting of Cartesian cells. These ar- guments are amplified in parallel-computing environments, where mesh modifications often imply rebalancing the load distributed to the var!ous processors by exchanging ceils across subdomain boundaries. Following this line of reasoning, the present work adopts a fully-implicit temporal operator to decouple the timestep from the local mest, scale. It leverages the same multigrid smoother used by the steady-state solver[15] by embedding it in a dual- time framework. The implicit approach means that finer mesh simulations do not automatically require more remeshings than coarse simulations. Simulations proceed at a timestep dictated by the appropriate physics rather than stability constraints. As noted earlier, special subclasses of arbitrary body motion can be treated by rigid motion of tile entire domain. The current work implements a rigid domain motionwithintheArbitraryLagrangian-Euleria(nALEJformulationofHirt etal.[17]. In this Lagrangian framework, a single transformation matrix can be applied to all the faces in a Cartesian mesh, and since these faces all have the same orientations, the transformation carl be inexpensively precomputed, stored and applied. Relative boundary motion is superimposed on top of the rigid domain motion using an adaptive re-meshing strategy which does not rely on explicit cell merging. Since the geometry moves relative to the volume cells, this aspect of the implemen- tation is Eulerian. A space-time analysis is used to ensure conservation, and to develop a hierarchy of approximations to the moving boundary. The basic analysis is presented in detail in I-D, with the complexity which arises when extending to mul- tiple dimensions also discussed. The final section presents initial numerical results in one, two, and three dimensions using the implicit, relative-motion sclmme, including comparison with analytic solutions and experimental d_ta. 2 Dual-time formulation In order to leverage the infrastructure of the steady-state flow solver outlined in Ref. [15], a dual-time formulation (cf. Refs. [i8, i9]) was developed for the time- dependent scheme, dQ d--T+ R*(Q) = o (1) c_Q R*(Q)=_+R(Q) where 7-is referred to here as "pseudo-time", and is the iterative parameter, and t is the physical time. Q is the vector of conserved variables, and R (Q) is an appropriate 1 numerical quadrature of the flux divergence, V fs f' ndS. As _ _ 0 the time- dependent formulation is recovered. The multi-grid smoother described in [15] is used to converge the inner pseudo-time integration. An explicit, multi-stage, pseudo- time-integration scheme is utilized to converge the "irner loop" in Eqn. 1. This is similar to the scheme outlined by Jameson[20], however, the semi-implicit approach of Melson et a1.[21] is used here for the physical time-derivative term. Various time-dependent schemes can be constructed for Eqn. 1 by appropriately discretizing the time derivative. As noted in the Intrcduction, it's desirable to uti- lize an unconditionally-stable, implicit scheme to allow a large timestep to be chosen based upon physical considerations rather than a potentially smaller stability-limited timestep. Beam and Warming[22] outline a family Df consistent time-dependent schemes that utilize three time levels, Q_-I Q,_, Qn+l These are given by fii*(Q) -- (1 + ()Q,_+I _ (1 + 2_)Q,_ + _Qn-1 2st (2) + OR (q,,+l) + (1 - 0 + ¢)R (qn) _ OR (q,_-l) I 0 I _ I _ I Method " I Order I 1 0 I 0 Backward Euler 1 1/2 0 0 Trapezoidal 2 1 t 1/2 0 2nd-order backward 2 3/4 0 -1/4 Adam's type 2 1/3 -1/2 -1/3 Lee's type 2 1/2 -1/2 -1/2 Two-step trapezoidal 2 5/9 -1/6 -2/9 A-contractive 2 Table 1: Partial listofA-stable, two-and three-level methods (Eqn. 2) fromBeam and Warming[22]. Table 1 contains a partial list of the A-stable, three time level methods that can be formulated using Eqn. 2. 3 ALE formulation In many applications with moving geometry, the motion of the geometry can be decomposed into a "uniform" rigid-bodjy ....,,,,°_,,.,.,.,._..; ,,1o,,;+r_elative motion confined to a subset of the domain. Examples include rotating airframes with dithering canards[23], rotoreraft, or stage separation from space vehicles. It's desirable to simulate such a motion using a rigid-body motion of the entire comput._Ltional domain, and treat the relative motion within the moving domain separately. This again limits the amount of computational work that is required to process the moving geometry. An ALE formulation (eft Hirtet al. [17]) was utilized in order to account for the rigid-body motion of the computational domain. This is accomplished by modifying the flux through a boundary to account for the motic)n of the boundary. For the inviscid flux vector used here, this becomes (3) f . n = pu,_u + pn pzt,_ } pu,_e + pu •n where is the velocity relative to the moving boundary, and un is the velocity of the moving domain. Hence the convective part of the flux is modified to account for the motion of the boundary, compared to the treatment for a fixed domain. A modified form of van Leer's flux-vector splitting (FVS) [24J is used with the ALE formulation. This modification generalizes the scheme by using the Mach number relative to the moving boundary when determining the characteristic speeds of the system. If the relative Mach number is not used the scheme will not be Galilean invariant; simulations with a static domain will not provide the same numerical results as those which use a moving domain to compute the same physical problem. The subsonicflux forvanLeer'sFVScanbewrittenasacombinationofconvectiveand acousticterms(forastaticdomain)as {0}v t*-•n = pcu 4-M,_ (2:FM,,) Pnn pcH o} :F (2:FM,,)2 0 (4) "7-7 " t pc where M,,c = u. n and M_ = _¼ (1 + M,,) 2. The numerical flux across a boundary which separates two domains (left and right) then becomes f(QL, QR) = f+(Qc) + With an ALE formulation, the Mach number relative to the moving boundary is M,_c = (u - u_) •n. If this relative Ivlach number is simply substituted into van Leer's FVS when computing on a moving domain, the energy equation will not be consistent, in the sense of f(Q, Q) = f(Q). This is due to combining the total energy and pressure work terms into the total enthalpy H in the numerical energy flux in Eqn. 4. Physically, the pressure work due to the domain motion (pua) is identically zero, and hence this term does not appear in the energy flux in Eqn. 3. A Inodified form of van Leer's FVS was developed for the ALE scheme. The mass and momentmn flux are unmodified from the original form, only the energy flux is changed. Examining the energy" flux in Eqn. 4, it's seen that the last term disappears when QL = QR. Neglecting this last term in the energy flux, the total enthalpy is split into the total energy and pressure work, and these terms are treated separately in the same manner as the momentum equations. The modified scheme then becomes {0} f±. n = M,_ pcu :LM_ (2 :F M. ) pn (5) pce pu. n where M,_, is now based upon the relative Mach number. The convective terms and acoustic terms are treated consistently in all 3equations of the system. This numerical flux is similar to the WPS scheme developed by Agarwal and Halt[25] The modified scheme for the energy flux maintains t:_esmoothness property of van Leer's original FVS. It is continuously-differentiable through M,, = El, as can be seen in Fig. 1, where the normalized energy flux (o@) is plotted against the Mach number. Both FVS schemes smoothly approach the asymptotic limit of pure upwinding. The geometry of the domain is expressed in the moving coordinate system, and must be transformed to the inertial system where the equations of motion are spec- ified. With a Cartesian scheme, the faces of each cmnputational cell have normals pointing in one of the Cartesian directions. These normsls all transform to the inertial system similarly, and are simply pre-computed and stcred prior to each timestep. [- v nLee, " /"J - -#-'2 [-- " van Leer F- _" -- l-- O_ _-- ' Mod. van Leer F _,_ = ! ..... 0 , 1 L I _ , L -0.5 0 0.5 Mach Number Figure 1: Normalized energy flux for tile FVS schemes, Eqns. 4 and 5. The ALE formulation was utilized to sinmlate a transonic NACA 0012 pitch- ing airfoil (cf. AGARD Repor_ 702 [26]). The 2nd-order backward time-integration scheme was used to compute the unsteady motion star_;ing from a converged steady- state solution. The sinmlation uses 100 timesteps per complete cycle of the airfoil, and an isotropic mesh with a finest resolution of 0.001 chord. As the airfoil oscillates, the shock transitions between the upper and lower surfaces of the airfoil, and the fluid state is path-dependent. This hysteresis is evident in the normal force history plotted in Fig. 2. After an initial transient of approximately 1/2 cycle, the solution is periodic. At a given angle of attack, the normal force is multi-valued depending upon whether the airfoil is pitching up or down. The computed variation of normal force is in good agreement with the experimental data. 4 Relative motion As an introduction to computing moving boundaries with a non-body-fitted Carte- sian scheme, Fig..9? contrasts three methods of moving geometry during a timestep: overset meshes, deforming meshes, and the current Cartesian approach. Both the overset and deforming mesh approaches use body-fitted volume meshes. Overset ap- proaches (cf. [1, 27, 28]) are relatively easy to implement and efficient, as the volume mesh processing is done a priori; only (lower-dimensional) boundary interpolations are updated at each step. Overset methods cannot maintain conservation across composite grid boundaries without complex re-meshing of the overlap region. It also can be difficult, and labor-intensive to maintain a compatible spatial resolution across composite grids, and to process the boundaries when components are close to- 0.5 ' I ' I T I ' 1 l O.25 0 D • Exp. (AGARD 702) m Comp. -Rigid-Domain Motion _0.51 , -I2 , -1I , 0I , 1I , 2l , Angle of Attack (deg) Figure 2: Variation of normal force coeflicient with angle of attack for oscillating NACA 0012. (3/I_ = 0.755, a(t) = 0.016 + 2.51sin (27r62.5t)). The simulation uses 100 timesteps per complete cycle of the airfoil, and an isotropic mesh with a finest resolution of 0.001 chord. Experimental data from [26]. gether. With deforming meshes (of. [3.29, 30]), the volume mesh deforms in response to the surface motion, and after large deformations e. volume re-meshing and (of- ten non-conservative) interpolation to the new mesh is. performed. Deforming-mesh approaches can be conservative over a timestep, making them attractive for small de- formations, however for large timesteps (which lead to large boundary movements), the quality of the cell stencil can severely degrade due to the deformations. In the current Cartesian approach, the moving boundary "sweeps" through a fixed Eulerian mesh over a timestep. The cells within the swept region of the domain change volume and shape over the timestep, and cells can appear or disappear (or both) as well. The space-time geometry associated with these swept cells is complicated, however it is possible to formulate conservative schemes as the boundary' is impermeable, hence the flux through the moving boundary is always known to some order of accuracy. Formulating a non-body-fitted Cartesian scheme which maintains conservation for large timesteps (large motions) of a moving boundary, while still retaining efficiency and robustness is the challenge for the current (and fllture) research. 4.1 Governing Equations The motion of a solid body through an inviscid fluid discretized by a fixed Carte- sian mesh is governed by the same ALE set of conservation equations (a Lagrangian body moves through an Eulerian mesh) as the rigid-domain motion of the previous ! section, S = OV dt "(0 (t) } (G) f. n = p_zTu, + pn pu_ p_ne + pu •n = (u- w).n Here w is the velocity of the moving boundary with respect to the Eulerian frame, and is used to differentiate from the rigid domain velocity tin. For the current discussion the rigid-body motion of the domain will be ignored, anJ only the relative motion will be formulated. No changes to the current scheme are required when the rigid-domain motion is superposed. For a cell in a Cartesian mesh swept by a moving boundary, Eqn. 6 can be simpli- fied as the boundaries of the cell are fixed, hence w •n = 0, except for the motion of the solid surface through the volume, for which u. n = w. n holds (i.e. the convective portion of the flux is zero through an impermeable surface). Eqn. 6 is preferred over a simplified form however, as it emphasizes the deforming-cell nature of the problem. Applying Liebniz' rule to the left side of Eqn. 6 givss d QdV = --_ + Qw ndS (7) dt (t) (t) (t) and the right-hand-side flux through the boundary sur[ace can be decomposed as -_s f'ndS=-_s [Q(u-w)+P].ndS (8) (t] (t) (0) where P = p6 is the acoustic portion of the flux. Balancing terms gives pu /v(t) _t dV + _(t) [Qu + P] .n:tS = O (9) Equation 9 can be recast into a space-time divergence form by applying Gauss's thereom to the spatial volume integral. Defining the 4-D space-time normal as _ = {[, n}, where [ is the normal in the time direction, n is the unit spatial normal, and Q = {Q, Qu + P}, the space-time conservation equation is Q. : o 0o) where f_ is the boundary of the space-time volume. Expanding Q for clarity gives = pu t, puu + p6 n (11) pe peu + pu ......... 0 ..7 -I u< ,oso t t ,@!Qu + P Qu+P ,,Pl I-+n J"_X Figure 3: One-dimensional space-time cell. The impermeable moving boundary (shown in red) crosses the fixed left face of the cell at time tc. The impermeable boundary has a normal with elements inboth space and time. Note that the velocity of the moving solid surface no longer explicitly appears in Eqn. 10 as the motion of the boundaries is accounted for in the direction and "area" of the space-time boundary. Fig. 3 shows a 1-D space-time cell volume for a Cartesian cell as a moving-boundary crosses the left cell face. The impermeable portion of tile space-time cell boundary (red boundary in Fig. 3 with slope l/w) has a normal with elements in both space and time, and is analogous to the role of the boundary velocity in Eqn. 6. In essence, the mixed Lagrangian/Eulerian Jbrmulation of Eqn. 6 has been converted to an Eulerian formulation in space-time. Equations 6 and l0 explicitly satisfy the so-called Geometric Conservation Law (GCL) (cf. Thomas and Lombard[31]), and no supplementary information is required. The GCL is usually presented as dS dV t; w ndS (12) dt (t) (t) which states that the change in cell volume is equivalent to the area swept by the moving boundary. The supplementary information the.t the cell must close, f ndS = 0, is also required. The space-time analog to the GCL is _dfi 0 (13) which states that the space-time cell must close. Whether numerical schemes are built for Eqn. 6 or Eqn. 10 is largely a matter of convenience, as both are mathematically equivalent. In the current work, Eqn. 6 is discretized directly, since it maintains a similar formulation as the governing equations for both static and rigid-domain motions. This leverages the existing flow solver infrastructure for multigrid, dual-time, etc., for solving the relative motion equations. 9 (a) Behind (Emerging) 1 t= n+l _ . i I Jo 1 t=/'/ 1 J r Ahead (dissappearing) (b) k i AJ' t = n+l -" im -" tc '. _l t'_ i J k Figure 4: One-dimensional motions of an impermeable boundary for motion less than the uncut hexahedron spacing (CFL_ < 1). and spatial directions. The temporal (convective) contribution is identically zero, as u. n = w. n for an impermeable boundary. This leaves only the pressure contribution to the moving-boundary flux, which in the current 1-D example becomes tn + l _in_r l (16) L /?d_ = pfwdntJt where Pw represents an approximation for the wall pressure to some order of accuracy. In other words, pressure acts only in the spatial directions, and hence the space-time area over which the surface pressure acts is the projection along the spatial axes of the moving wall area. In the vector notation of the previous section = (p. n)w (lr) An example in multiple dimensions is shown in Fig. 5 where a boundary moves through a group of hexahedra in 2-D. Examining the space-time cell k, the pres- sure acts on the projection of the moving front along the -z and +y axes. The space-time areas over which the pressure acts are thei_ the triangular projections in Fig. 5b. 11