- Email: [email protected]
Mathematical Formulation for Preconditioned Compressible Flow Solver
Chapter 2
Mathematical Formulation for Preconditioned Compressible Flow Solver 2.1 BASIC GOVERNING EQUATIONS FOR FLUID FLOWS The 3D unsteady NavierStokes equations may be written in a number of forms. One common form is as follows: Continuity: @ρ @ðρuÞ @ðρvÞ @ðρwÞ 1 1 1 50 @t @x @y @z
ð2:1Þ
u-Momentum: @ðρuÞ @ðρu2 1 pÞ @ðρuvÞ @ðρuwÞ @τ xx @τ xy @τ xz 1 1 1 5 1 1 @t @x @y @z @x @y @z
ð2:2Þ
v-Momentum: @ðρvÞ @ðρuvÞ @ðρv2 1 pÞ @ðρvwÞ @τ yx @τ yy @τ yz 1 1 1 5 1 1 @t @x @y @z @x @y @z
ð2:3Þ
w-Momentum: @ðρwÞ @ðρuwÞ @ðρvwÞ @ðρw2 1 pÞ @τ zx @τ zy @τ zz 1 1 1 5 1 1 @t @x @y @z @x @y @z
ð2:4Þ
Energy equation: @ðρet Þ @ðρet uÞ @ðρet vÞ ðρet wÞ 1 1 1 @t @x @y @z 5
@ðκð@T=@xÞÞ @ðκð@T=@yÞÞ @ðκð@T=@zÞÞ @ðupÞ @ðvpÞ @ðwpÞ 1 1 2 2 2 @x @y @z @x @y @z
1
@ðuτ xx 1 vτ xy 1 wτ xz Þ @ðuτ yx 1 vτ yy 1 wτ yz Þ @ðuτ zx 1 vτ zy 1 wτ zz Þ 1 1 @x @y @z ð2:5Þ
Computational Fluid-Structure Interaction. DOI: https://doi.org/10.1016/B978-0-12-814770-2.00002-7 © 2019 Elsevier Inc. All rights reserved.
11
12
Computational Fluid-Structure Interaction
where the volumetric heating, such as absorption or emission of radiation and body source terms, such as gravity terms are not included. ρ is density; u, v, and w are the Cartesian components of velocity along x-, y-, and z-axes; p is pressure; and T is temperature. They can be related by the state equation of perfect gas p 5 ρRT, where R is the gas constant, which is 287 m2 =ðs2 KÞ for air at standard conditions. Furthermore, et is the total energy defined as et 5 e 1 ð1=2Þðu2 1 v2 1 w2 Þ, where e is internal energy per unit mass of the fluid and e 5 Cv T, where Cv is the specific heat at constant volume and Cv 5 R=ðγ 2 1Þ. γ is the ratio of specific heats and γ 5 1:4. If we assume that the fluid considered is a calorically perfect gas, then the NavierStokes equations are closed by the state equation for perfect gas: p 5 ðγ 2 1Þρe T5
ðγ 2 1Þe R
H 5e1
p 1 1 ðu2 1 v2 1 w2 Þ ρ 2
Finally, the viscous stresses are related to relations: @u 2 @u @v @w 1 1 τ xx 5 2μ 2 μ ; @x 3 @x @y @z @v 2 @u @v @w 1 1 τ yy 5 2μ 2 μ ; @y 3 @x @y @z @w 2 @u @v @w 2 μ 1 1 τ zz 5 2μ ; @z 3 @x @y @z
the velocity field by the Stokes @v @u 1 τ xy 5 τ yx 5 μ @x @y @u @w 1 τ xz 5 τ zx 5 μ @z @x @w @v 1 τ yz 5 τ zy 5 μ @y @z
The dynamic viscosity μ and conductivity κ are properties of the fluid and are functions of temperature. These two quantities are related by the Prandtl number: Pr 5
μCp κ
For air, the Prandtl number is around 0.72 at room temperature. Equations of fluid motion may be nondimensionalized to achieve certain objectives. First, it would provide conditions upon which dynamic and energetic similarity may be obtained for geometrically similar situations. Second, the solution of such equations would usually provide values within limits between zero and one, which can enhance accuracy for numerical computation.
Mathematical Formulation Chapter | 2
13
The nondimensional variables used in this study are defined as follows: x y z u v w ðx ; y ; z Þ 5 ; ; ; ; ðu ; v ; w Þ 5 L L L UN UN UN t 5 ρ 5
t p 2 pin et H ; p 5 ; et 5 ; H 5 2 2 L=U N ρN ðU N Þ ðU N Þ ðU N Þ2
ρ T μ ρN U N L γμ Cp r T ; γ5 ; T5 ; μ 5 ; Re 5 ; q 52 N 2 ρN μN μN Pr Cv ðU N Þ =Cv
where pin is the inlet or reference pressure; L is the characteristic length of the computed model; U N is the inflow velocity; and μ is the dynamic viscosity of the fluid. The variables with a superscript asterisk sign ( ) are nondimensional parameters and the asterisk sign will be dropped in subsequent equations for sake of convenience. It should be noted that we subtract a constant value (the reference pressure) from the pressure term to control the round-off errors for low-speed flows, which is found to be critical in controlling computational errors in the momentum equations for low-speed compressible flows. The use of gauge pressure is a common treatment for incompressible solvers because only pressure gradients are needed for all calculations. For compressible solvers, the absolute value of pressure must be used when dealing with the energy equation and state equation of gases. But when we use nondimensional absolute pressure at low Mach numbers, it becomes extremely large although the pressure gradients in momentum equations are small. The use of gauge pressure can avoid performing the addition and subtraction operations between two extraordinarily large values of nondimensional absolute pressures. In our experience, the result obtained with gauge pressure is more accurate than without using it. And the relative error can be up to 4%5% of the result. It should be noted that the original absolute pressure (designated by p0 hereafter) should be used instead of gauge pressure whenever the state equation of gas is involved. The nondimensional equations of fluid motion may be expressed in vector form as @W c @Ei @Fi @Gi @Ev @Fv @Gv 1 1 1 5 1 1 @t @x @y @z @x @y @z 2 3 ρ 6 ρu 7 6 7 7 Wc 5 6 6 ρv 7 4 ρw 5 ρet
ð2:6Þ
ð2:7Þ
14
Computational Fluid-Structure Interaction
2
3 ρu 6 ρu2 1 p 7 6 7 7 Ei 5 6 6 ρuv 7 4 ρuw 5 ðρet 1 pÞu 3 2 0 7 6 τ xx 7 1 6 7 6 τ xy Ev 5 7 6 ReN 4 τ 5 xz uτ xx 1 vτ xy 1 wτ xz 2 qx 2 3 ρv 6 ρvu 7 6 2 7 7 Fi 5 6 ρv 1 p 6 7 4 ρvw 5 ðρet 1 pÞv 3 2 0 7 6 τ yx 7 1 6 7 6 τ yy Fv 5 7 6 ReN 4 τ 5 yz uτ yx 1 vτ yy 1 wτ yz 2 qy 3 2 ρw 7 6 ρwu 7 6 7 6 Gi 5 6 ρwv 7 4 ρw2 1 p 5 ðρet 1 pÞw 2 3 0 6 τ zx 7 7 1 6 6 7 τ zy Gv 5 7 ReN 6 4 τ zz 5 uτ zx 1 vτ zy 1 wτ zz 2 qz
ð2:8Þ
ð2:9Þ
ð2:10Þ
ð2:11Þ
ð2:12Þ
ð2:13Þ
2.2 LOW-SPEED PRECONDITIONING FORMULATION One difficulty with compressible NavierStokes solvers is their slow convergence rates and even unstable solutions for low-Mach-number flows. This difficulty can be traced to a disparity between the acoustic and convective speeds [110], and can be addressed by a preconditioning algorithm. Previous work in this area has been reported by Venkateswaran and
Mathematical Formulation Chapter | 2
15
Merkle [47], Turkel [8], Van Leer et al. [9], and Weiss and Smith [10]. The applications of the preconditioning methods have been found in the computation of steady flows without considering arbitrarily moving objects in the flow field. The preconditioned NavierStokes equations for 3D compressible unsteady flows can be given in vector form, explicitly expressing the conservation laws of mass, momentum, and energy. We also introduce, in the equations, pseudo time terms to provide pseudo time marching for their numerical solutions: Γ1 where
@W p @W c 1 1 rUFi 5 rUFv @τ @t 3 ρ 6 ρu 7 6 7 7 Wc 5 6 6 ρv 7 4 ρw 5 ρet 2 3 p 6u 7 6 7 7 Wp 5 6 6v 7 4w 5 T0 3 2 ρU 6 ρuU 1 pi 7 7 6 7 Fi 5 6 6 ρvU 1 pj 7 4 ρwU 1 pk 5 ρHU 2 3 0 6 τx 7 7 1 6 6 τy 7 Fv 5 6 7 ReN 4 τ 5 z τ UU 2 q
ð2:14Þ
2
ð2:15Þ
ð2:16Þ
ð2:17Þ
ð2:18Þ
τ is the pseudo time and Γ1 is the preconditioning matrix in the pseudo time terms for low-Mach-number flows, which is defined in the appendix. W c and W p are the vectors of conservative and primitive dependent variables, respectively; Fi and Fv are the inviscid convective flux and viscous flux vectors. Furthermore, we have the following formulas: U 5 ui 1 vj 1 wk
16
Computational Fluid-Structure Interaction
τ 5 τ xi 1 τ yj 1 τ zk τ i 5 τ ix i 1 τ iy j 1 τ iz k @ui @uj 2 τ ij 5 μ 1 2 δij rUU 3 @xj @xi q 5 qx i 1 qy j 1 qz k T 0 5 p0 =ρ 5 c2 =γ i, j, and k are the three unit vectors in three Cartesian directions; τ ix , τ iy , and τ iz are the viscous stresses.
2.3 DISCRETIZATION METHODS The 3D equations (2.14) is transformed into an integral form and discretized on an unstructured grid. A cell-vertex finite-volume scheme is adopted here. For every vertex, a control volume is constructed using the median duals of the tetrahedral cells by connecting all the neighboring centroids of edges, surfaces, and tetrahedra. As shown in Fig. 2.1 is a part of the control volume associated with vertex P within tetrahedron PABC. Spatial discretization is performed by using the integral form of the conservation equations over the control volume surrounding node P: ððð ððð ððð ððð @Q0 1 @W c dV 1 dV 1 rUFi dV 2 rUFv dV 5 0 ð2:19Þ cv @τ cv @t cv cv
FIG. 2.1 Construction of control volume within a tetrahedron for a node P.
Mathematical Formulation Chapter | 2
17
Noted that a new variable Q0 1 has arisen as @Q0 1 =ð@τÞ 5 Γ1 @W p =ð@τÞ, and the Jacobian Γ1 5 @Q0 1 =ð@W p Þ. So that we have: @Q0 1 @Q0 1 @W p @W p 5 5 Γ1 @τ @W p @τ @τ
2.3.1 Edge-Based Method and Cell-Based Method The convective term in Eq. (2.19) is transformed into a summation: ððð ð nbseg X rUFi dV 5 Fi UndS 5 ½ðFi Þij UnΔSk cv
Scv
ð2:20Þ
k51
where nbseg is the number of the edges associated with node P; ðFi Þkij is the inviscid flux through the part of control volume surface associated with edge k; and n is the unit normal vector of the control volume surface. Finally, ΔSk is a part of the control volume surface associated with edge k (as shown in Fig. 2.1, if k is edge PC, then 1O2c is ΔSk ). Therefore all the fluxes are calculated for the edges and then collected at the two ends of each edge for updating of flow variables in time marching, which is the so-called edgebased data structure for fast inviscid flux computation and efficient data storage and retrieval. On the other hand, the viscous term is calculated using a cell-based method for the same reason of high computational efficiency: ððð ð ncell X ½Fv UnΔSc i rUFv dV 5 Fv UndS 5 ð2:21Þ cv
Scv
i51
where ncell is the number of tetrahedral elements associated with node P and ΔSci is the part of control volume surface in cell i. By using the following relation: ð dS 5 0: Scv
The total vector surface of the control volume in a cell i becomes nΔSci 5 ð1=3ÞðnΔSpi Þ. Thus the calculation of viscous terms can be simplified as ncell X i51
½Fv UnΔSc i 5
ncell 1X ½Fv UnΔSp i 3 i51
ð2:22Þ
where nΔSpi is the surface vector of the face opposite node P of the tetrahedron under consideration, where the ðFv Þi is calculated at the center of the tetrahedron with a node P, and can be obtained by using the Green’s theorem based on the variables at the four vertices of the tetrahedron. Similar to the
18
Computational Fluid-Structure Interaction
Galerkin type of formulation, the gradient of a flow variable φ at the center of a tetrahedron is evaluated as follows: P4 P 1 4i51 φi Si i51 φi 9Si 52 ð2:23Þ gradφc 5 2 3 27V V where φi is the flow variable at a vertex i of the tetrahedron; Si is the surface area that is opposite to node i; and V is the volume of the tetrahedron. Gradients at the vertices are obtained by a volume averaging of the gradients at the centers of cells associated with the vertex under consideration.
2.3.2 Roe’s TVD Scheme for Inviscid Flux in Edge-Based Method The Riemann problem is an initial problem with piecewise constant initial data. Its exact solution represents the real physical characteristics of a flow with several families of waves and their propagation. When the wave is a shock, the equations for this initial problem cannot be explicitly solved and an iterative solution has to be used. In order to reduce the amount of calculation, some approximate treatments have been employed in numerical solutions. In this work, a high-order Roe’s TVD scheme for compressible flow for arbitrary unstructured 3D grids has been adopted. In Roe’s approach [11], averaged values of density, velocity, and enthalpy at the interface between the two initial states at two neighboring points are calculated by using a special averaging procedure, called Roe’s averaging: pffiffiffiffiffipffiffiffiffiffi ρ 5 ρL ρR pffiffiffiffiffi pffiffiffiffiffi ðu; v; wÞL ρL 1 ðu; v; wÞR ρR pffiffiffiffiffi pffiffiffiffiffi ðu; v; wÞ 5 ð2:24Þ ρL 1 ρR pffiffiffiffiffi pffiffiffiffiffi H L ρL 1 H R ρR p p ffiffiffiffiffi ffiffiffiffiffi H5 ð2:25Þ ρL 1 ρR where L and R represent the two neighboring points on an edge, and H denotes total enthalpy which has the following relation: H5e1
p0 1 1 u2 1 v 2 1 w 2 2 ρ
The speed of sound can be derived as: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 2 c5 H 2 ðu 1 v 1 w Þ ðγ 2 1Þ 2
ð2:26Þ
ð2:27Þ
Mathematical Formulation Chapter | 2
19
The inviscid flux ðFi Þij is evaluated based on Roe’s approximate Riemann solver (or the flux difference splitting scheme): 1
ðFi Þij 5 FLi 1 FRi 2 A ðW R 2 W L Þ 2 1 L 5 Fi 1 FRi 2 jΔFj ð2:28Þ 2 5 FLi 1 ΔF2 5 FRi 2 ΔF1 where
FLi 5 Fi W Lc FRi 5 Fi W Rc
and W Lc and W Rc are the left and right conserved state vectors on the two neighboring points of an edge. The flux difference is: ð2:29Þ jΔFj 5 ΔF1 2 ΔF2 8 3 39 2 2 1 0 > > > > > > > 7 7> 6 6 > > u ΔU Δu2n > > 7 7 6 6 x n = < 7 7 6 6 Δp 6 6 7 7 6 6 1ρ ΔF 5λ1 Δρ2 2 6 v Δv2n ΔU y n 7 7> 6 > ðcÞ 6 > 7 7> 6 > > > 5 5> 4w 4 Δw2nz ΔUn > > > > ; : 2 ðU n Þ =2 uΔu1vΔv1wΔw2U n ΔUn 3 2 1 7 6 6 u1nx c 7 7 6 Δp1ρcΔU n 6 v1ny c 7 1λ46 7 6 2 2ðcÞ 7 6 4 w1nz c 5 H 1U n c 3 2 1 7 6 6 u6nx c 7 7 6 Δp2ρcΔUn 6 1λ56 v6ny c 7 2 7 6 2ðcÞ 7 6 4 w6nz c 5 H 6U n c ð2:30Þ where ðU n Þ2 5 u2 1 v2 1 w2 ΔUn 5 Δunx 1 Δvny 1 Δwnz U n 5 unx 1 vny 1 wnz
20
Computational Fluid-Structure Interaction
Furthermore, the eigenvalues of A are λ1 5 U n , λ4;5 5 U n 6 c. Because of the introduction of preconditioning matrix Γ1 , the inviscid fluxes, ðFi Þijk , through the face k is now reformulated as 1 1 @Fi 0 k ðFi Þij ððFi Þi 1ðFi Þj Þk 2 @Q0 ðδQ1 Þk 2 2 1 k ! ð2:31Þ 1 1 @Fi @W p @Q01 5 ððFi Þi 1ðFi Þj Þk 2 @W p @Q0 @W p ððW p Þj 2ðW p Þi Þk 2 2 1 k
k
Note that we have retained the variable, Q0 1 , in computing this flux. Defining the Jacobian in the normal direction as @Fi ðHp Þk 5 @W p k And using the previously defined Jacobian Γ1 5 @Q0 1 =@W p , then the above expression becomes 1 1 ðFi Þkij ððFi Þi 1ðFi Þj Þk 2 Hp Γ21 1 k Γ1k ððW p Þj 2ðW p Þi Þk 2 2
Drop the subscript k in the flux vector and the Jacobian with the assumption that the fluxes and Jacobians all correspond to conditions in the normal direction on the given control volume surface. And after some simple algebraic derivations, we have 1 1 ðFi Þij ððFi Þi 1 ðFi Þj Þ 2 Γ1 Γ21 1 Hp ððW p Þj 2 ðW p Þi Þ 2 2
ð2:32Þ
2.3.3 Upwind-Biased Interpolation Combined with the third-order Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) interpolation, it can produce accurate and stable solution on unstructured grids. The left and right state vectors W L and W R at a control volume surface are evaluated using a nominally third-order upwind-biased interpolation scheme: 1
1 ð1 2 κÞΔ2 i 1 ð1 1 κÞΔi 4 i 1h 2 W R 5 W j 2 ð1 2 κÞΔ1 j 1 ð1 1 κÞΔj 4 WL 5 Wi 1
where 2 Δ1 i 5 Δj 5 W j 2 W i 1 Δ2 i 5 W i 2 W i21 5 2ijUrW i 2 ðW j 2 W i Þ 5 2ijUrW i 2 Δi 2 Δ1 j 5 W j11 2 W j 5 2ijUrW j 2 ðW j 2 W i Þ 5 2ijUrW j 2 Δj
ð2:33aÞ ð2:33bÞ
Mathematical Formulation Chapter | 2
21
Therefore substituting the above equations into Eqs. (2.33a) and (2.33b), the final equations based on upwind-biased interpolation scheme is shown as follows: 1
ð1 2 κÞijUrW i 1 κΔ1 i 2 i 1h W R 5 W j 2 ð1 2 κÞijUrW j 1 κΔ2 j 2 WL 5 Wi 1
ð2:34aÞ ð2:34bÞ
where κ is set to 1/3, which corresponds to a nominally third-order accuracy. ij is the vector representing the edge, which points from node P to its neighboring node under consideration. The gradients of W at i and j are calculated by volume averaging the gradients of the cells that surround i and j. Finally, for a given node P, the spatially discretized Eq. (2.19) form a system of coupled ordinary differential equations. This can be reformulated as: @Q0 1 @W c ΔVcv 1 ΔVcv @τ @t 8 9 ncell <1 nbseg = X X 1 52 ½ððFi Þi1ðFi Þj Þ2Γ1 jΓ21 ½Fv UnΔSP i 1 Hp jððW p Þj2ðW p Þi Þk UðnΔSÞk 2 :2 k51 ; 3 i51 52R ð2:35Þ
where R represents the residual, which includes the convective and diffusive fluxes and ΔVcv is the control volume of node P. the overbar in Eq. (2.35) denotes the cell-averaging value. An implicit scheme is adopted for Eq. (2.35) and the time-dependent term is discretized using a second-orderaccurate backward differencing scheme: n11 n11 n n21 n21 @Q0 1 1:5ΔVcv W c 2 2:0ΔVcv W nc 1 0:5ΔVcv Wc n11 ~ n11 5R ΔVcv 5 2 R 2 @τ Δt ð2:36Þ where the superscript (n 1 1) denotes the physical time level (n 1 1)Δt and ~ n11 Þ is the new modiall the variables are evaluated at this time level, RðW p fied residual which contains both the time derivative and flux vectors. The derivative with respect to a pseudo time τ is discretized as n11 ΔVcv Γ1
W n11;m11 2 W n11;m n11;m11 p p 5 R~ dτ
ð2:37Þ
whose solution is sought by marching to a pseudo steady state in τ, where m and (m 1 1) denote the initial and final pseudo time levels. Once the artificial steady state is reached, the derivative of Wp with respect to τ becomes zero, ~ n11 5 0. Hence, the original unsteady and the solution satisfies R
22
Computational Fluid-Structure Interaction
NavierStokes equations are fully recovered. Therefore instead of solving the equations in each time-step in the physical time domain (t), the problem is transformed into a sequence of steady-state computations in the artificial time domain ðτÞ. Eq. (2.37) can be integrated in pseudo time by an explicit five-stage RungeKutta scheme. However the pseudo time-step size may be severely restricted if the physical time-step size is very small. Hence, a fully implicit dual-time-stepping is adopted here.
2.3.4 Matrix-Free, Implicit, Dual-Time-Stepping Algorithm A Taylor series expansion is performed for the residual in Eq. (2.37) with respect to the pseudo time for node i, m11 ~m 1 R~ i 5 R i
nbseg X @R ~i ~i @R ΔðW p Þi 1 ΔðW p Þj @ðW p Þi @ðW p Þj j51
ð2:38Þ
where nbseg is the number of edges connected to i, which is also equal to the number of neighboring points connected to point i through the edges. An approximate flux function is introduced here to simplify the implicit timestepping calculation. The total flux (including both convective and viscous fluxes) across a control volume surface associated with a certain edge ij can be approximated as Fij
1
ðFi Þi Un 1 ðFi Þj Un 2 jλij jððW p Þj 2 ðW p Þi Þ 2
where λij is the spectral radius based on the preconditioned system which is associated with edge ij, and is given in the appendix. Then in all of the R(W) terms (Eq. (2.35)), we have " #
nbseg nbseg X 1 @ðFi Þij X 1 @ðFi Þij @Ri @Ri 5 1 jλij j and 5 2 jλij j : 2 @ðW p Þi 2 @ðW p Þj @ðW p Þi @ðW p Þj j51 i51 And for the physical time-dependent terms, we have ðW c Þm11 5 ðW c Þm i i 1
@W c ΔðW p Þi @W p
After combining all the residual terms at every node in the flow field into a vector and dropping the third term of the right-hand side of Eq. (2.38) as an approximation, we have n11 ~ n11;m11 5 R ~ n11;m 2 @Ri ΔðW p Þi 2 1:5 ΔVcv @W c ΔðW p Þi R i i @ðW p Þi Δt @W p n11;m 5 R~ i 2
nbseg X j51
ðHp;j ÞΔðW p Þi 2 1:5
n11 ΔVcv MΔðW p Þi Δt
Mathematical Formulation Chapter | 2
23
where Hp;j 5 ð1=2Þ½ð@ðFi Þij =@ðW p Þi Þ 1 jλij j, M 5 @W c =@W p . And the wholefield equivalent of Eq. (2.36) can then be rewritten as 0 1 nbseg X W n11;m11 2 W n11;m 1:5Δτ Δτ p 21 n11 @I 1 A p Γ21 ΔVcv ðΓ H Þ p;j 1 M1 1 n11 Δt ΔVcv Δτ j51 0 1 n11;m n n21 n11 n n21 1:5W ΔV 2 2:0W ΔV 1 0:5W ΔV cv cv cv c c c @2 2 Rn11;m A 5 Γ21 1 Δt ð2:39Þ that is, A~
W n11;m11 2 W n11;m p p n11 ~ n11;m ΔVcv 5 Γ21 1 R Δτ
thus W n11;m11 2 W n11;m p p ~~ n11;m n11 ΔVcv 5R Δτ where ~~ R
n11;m
X 21 ~ n11;m and A~ 5 I 1 1:5Δτ Γ121 M 1 Δτ R 5 A~ Γ21 Γ21 1 1 Hp;j n11 Δt ΔVcv j51 nbseg
Therefore n21 n11 n11;m 2 2:0ΔSncv W n 1 0:5ΔSn21 cv W ~Rn11;m 5 2 Rn11;m 2 K 1:5ΔScv W Δt Further approximation can be introduced in order to achieve matrix-free computation. If we employ point implicit treatment to the preceding equations, then only the diagonal terms in A~ are used in the pseudo timestepping. As a result, the equation for every node can now be written as W n11;m11 2 W n11;m p p ~~ n11;m n11 ΔVcv 5R Δτ where ~~ n11;m 5 A~21 Γ 21 R ~ n11;m R 1 ii i i
ð2:40Þ
24
Computational Fluid-Structure Interaction
and
2
21 A~ii
nbseg 1:5Δτ 21 Δτ X 21 Γ1 M1 5 diag4 I1 Γ1 Hp;j n11 Δt ΔVcv j51
!21 3 5:
Pseudo time-stepping is then performed on Eq. (2.40).
2.3.5 Time Integration In this work, a RungeKutta time-integration algorithm is used between each physical time-step to iterate the numerical solution in an artificial time τ until convergence is reached. For a five-stage RungeKutta scheme, the stage coefficients are α1 5
1 1 3 1 ; α2 5 ; α3 5 ; α4 5 ; α5 5 1 4 6 8 2
and we have m W ð0Þ c;p 5 W c;p
Δτ ~~ ð0Þ R W c;p ΔVcv Δτ ~~ ð1Þ ð0Þ R W c;p W ð2Þ c;p 5 W c;p 2 α2 C ΔVcv Δτ ~~ ð2Þ ð0Þ R W c;p W ð3Þ c;p 5 W c;p 2 α3 C ΔVcv Δτ ~~ ð3Þ ð0Þ R Wc;p W ð4Þ c;p 5 W c;p 2 α4 C ΔVcv Δτ ~~ ð4Þ ð0Þ R Wc;p 5 W 2 α C W ð5Þ 5 c;p c;p ΔVcv
ð0Þ W ð1Þ c;p 5 W c;p 2 α1 C
ð2:41Þ
W ðm11Þ 5 W ð5Þ c;p c;p
2.4 CONVERGENCE ACCELERATION TECHNIQUES To speed up the convergence rate, apart from the matrix-free, implicit, dualtime-stepping algorithm and local time-stepping scheme introduced in the previous sections, an implicit residual smoothing scheme developed for unstructured grids is employed. The smoothing equation for a vertex k can be expressed as follows: Rk 5 Rk 1 ε r2 Rk
ð2:42Þ
Mathematical Formulation Chapter | 2
25
where R is the original residual; R is smoothed residual; and ε is the smoothing coefficient, which can be defined as ( " # ) 1 CFL 2 21 ; 0 ε 5 max 4 CFL where CFL is the maximum CFL number of the basic scheme. The solution to the above equations can be obtained on an unstructured grid by using the Jacobi iterative method as follows:
ðmÞ
Rk 5
Rð0Þ k 1ε
numnodðkÞ P i51
ðm21;mÞ
Ri
1 1 ε UnumnodðkÞ
ð2:43Þ
where numnodðkÞ is the number of neighboring nodes of vertex k. Furthermore, the multigrid method can also be employed for convergence acceleration, which will be introduced and discussed in details in Chapter 5, The Multigrid Method.
2.5 BOUNDARY CONDITIONS Physically correct numerical boundary conditions are used to achieve timeaccurate solutions for both steady and unsteady flows. This kind of boundary conditions should be compatible with modern low-dissipation and lowdispersion numerical algorithms, which have very low-dispersion errors and require precise boundary conditions to avoid numerical instabilities and to control spurious nonphysical wave reflections at the computational boundaries. A high-order and high-resolution scheme as one of those numerical algorithms gives results sensitive to the boundary values in the computation, thus the quality of solutions even on interior nodes are dependent on the accuracy and correctness of boundary conditions [12]. Thompson [13,14] decomposed hyperbolic equations into wave modes of definite velocity and then specified characteristic boundary conditions for incoming waves. The starting point of his analysis was nonlinear Euler equations. The idea of his approach was that 1D characteristic analysis could be performed by consideration of the transverse terms as a constant source term. The amplitudes of outward propagating waves are defined entirely from the variables inside computational domain, while those of inward propagating waves are specified as the characteristic boundary conditions. In our studies, we choose to let the state variables of the far-field boundary be floating rather than fixed; in other word, we solve the characteristic equations on the boundary nodes. This method can make the boundary treatment nonreflective to the maximum extent and our solver more stable [15].
26
Computational Fluid-Structure Interaction
2.5.1 Far-Field Conditions In a cell-centered scheme, the Roe TVD upwind scheme described above can be extended to the boundaries conveniently. Since the flux-formula automatically selects the proper information from either side of the boundary face. One may add in an extra cell across the boundary in which the full state is prescribed at any time, even for subsonic inflow and outflow; while in our cellvertex scheme, the far-field boundary conditions for the flow equations are based on the flux vector splitting technique for the inviscid fluxes [16]. For a boundary node k with its control volume boundary face ΔSk , the inviscid flux across ΔSk is evaluated by using the Stegger and Warming splitting [16]. At the point upstream of the inflow boundary (denoted by “in”), ρin ; uin ; vin ; win and free stream Mach number MN are prescribed. Then, the inflow pressure, temperature, and total energy are determined by: 1 γðMN Þ2
ð2:44Þ
1 γðγ 2 1ÞðMN Þ2
ð2:45Þ
pin 5 Tin 5 ðρet Þin 5
1 2 pin ρin uin 1 v2in 1 w2in 1 2 γ21
ð2:46Þ
At the point downstream of the outflow boundary, the only one known state is the pressure. All of the other states must be extrapolated from interior flow field using characteristic variables [9,15]. For an arbitrary boundary normal direction n, we can write the variation of characteristic variables as 0
@ρ 2
1 @p c2
1
C 1 B C B @W 1λn C B sU@U C B B @W 2 C B C λn C B tU@U C B 3 C 5 @W λn 5 B @W 1 C B λn C B C B @p nU@U 1 @ @W 4 A B C ρc λn C B @W 5λn C B 1 A @ 2 nU@U 1 @p ρc 0
ð2:47Þ
The first characteristic variation is proportional to entropy variations, the second and third are related to variations in shear and tangent velocity, and the last two represent acoustic disturbances. Thus from @W 1λn 1=γ pout ρout 5 ρb ð2:48Þ pb
Mathematical Formulation Chapter | 2
27
where subscript “b” denotes the variable from the boundary node and “out” denotes flow state at the point downstream of the outflow boundary. From @W 4λn , @W 5λn , we have ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi γpout =ρout b 2 2 2 2 2 2 Un;out 5 un;out 1 vn;out 1 wn;out 5 un;b 1 vn;b 1 wn;b 1 γ21 From @W 3λn U t;out 5 U t;b 5 U b 2 Un;b n Therefore we have Uout 5 U t;out 1 Un;out n 5 U b 1 ðUn;out 2 Un;b Þn
ð2:49Þ
uout 5 ub 1 ðUn;out 2 Un;b Þnx
ð2:50Þ
vout 5 vb 1 ðUn;out 2 Un;b Þny
ð2:51Þ
wout 5 wb 1 ðUn;out 2 Un;b Þnz
ð2:52Þ
pout ρ jU out j2 1 out γ21 2
ð2:53Þ
ðρet Þout 5
As an example, the far-field boundary conditions treatment of a strong converged nozzle is shown in Fig. 2.2. The boundary fluxes are calculated using Stegger and Warming’s flux vector splitting formula [16,17]: Fi UnSk 5 ðFi 1 1 Fi 2 ÞUΔSk
Fi 6
3 2 ρ α 7 6 2γ 7 6 7 6 ρ 7 6 ðαu 1 βcnx Þ 7 6 2γ 7 6 7 6 ρ 7 6 ðαv 1 βcny Þ 7 6 2γ 56 7 7 6 ρ 7 6 ðαw 1 βcn Þ 7 6 z 7 6 2γ 6 0 17 7 6 2 2 2 7 6 ρ 4 @α u 1 v 1 w 1 2βcU 1 ϖ c2 A 5 n 2γ γ21 2
where Un 5 UUn λ1 5 UUn
ð2:54Þ
28
Computational Fluid-Structure Interaction
FIGURE 2.2 The inflow and outflow boundary conditions for a strong converged nozzle.
λ4 5 UUn 1 c λ5 5 UUn 2 c α 5 2ðγ 2 1Þλ16 1 λ46 1 λ56 β 5 λ46 2 λ56 ϖ 5 λ46 1 λ56 λ1 5 maxfλ; 0:0g λ2 5 minfλ; 0:0g It should be noted that at the inflow, the boundary nodal states are used for calculating Fi 1 , the preimposed inflow states are used for Fi 2 . Similarly, at the outflow, the boundary nodal states are used while calculating Fi 1 , the preimposed and extrapolated outflow states are used for Fi 2 .
2.5.2 Solid Wall/Slip Conditions At the solid wall, slip (for inviscid flows) or no-slip (for viscous flows) and no-injection boundary conditions are imposed, that is, the zero normal fluxes of mass, momentum, and energy are imposed. In addition, the solid surface is assumed to be adiabatic, that is: rUT 5 0
2.5.3 Far-Field Conditions for Preconditioned Systems The far-field calculations are based on characteristic variables (Reimann invariants). Thus at inflow, the incoming variables corresponding to positive eigenvalues are specified while the outgoing variables corresponding to negative eigenvalues are extrapolated. Once we change the time-dependent equations, we also change the characteristics of the system (though not the signs
Mathematical Formulation Chapter | 2
29
of the eigenvalues). Hence, it is also necessary to modify the boundary conditions for the preconditioned system. Here, the flux at a boundary is defined as Fi UnSk 5 ðFi 1 1 Fi 2 ÞUΔSk where Fi 6 has been redefined as Fi 6 5 XHp ;R Λ 6 X Hp ;L W p , where Λ 6 5 λi 6 jλi j=ð2Þ, and λi represents the eigenvalues of Hp (see appendix for details of λi ). X Hp ;R and XHp ;L are the right and left eigenvectors of Hp (see appendix). We first calculate 0 1 ð1 2 γÞnx p=ðγρÞ 1 nz v 1 nx T 0 2 ny w B ð1 2 γÞny p=ðγρÞ 2 nz u 1 nx w 1 ny T 0 C B C 0 B C z 5 XHp ;L W p 5 B ð1 2 γÞnz p=ðγρÞ 1 ny u 2 nx v 1 nz T C @ p 1 nx ρðλ1 2 βUÞu 1 ny ρðλ1 2 βUÞv 1 nz ρðλ1 2 βUÞw A p 1 nx ρðλ2 2 βUÞu 1 ny ρðλ2 2 βUÞv 1 nz ρðλ2 2 βUÞw ð2:55Þ And then,
Fi 6 5 X Hp ;R Λ 6 XHp ;L W p 5 X Hp ;R Λ 6 z 1 0 ððβU 2 λ2 Þλ1 6 z4 1 λ1 2 βUÞλ2 6 z5 =S C B ny λ0 6 z3 2 nz λ0 6 z2 1 nx λ1 6 z4 2 nx λ2 6 z5 =ðρSÞ C B B nz λ0 6 z1 2 nx λ0 6 z3 1 ny λ1 6 z4 2 ny 6 λ2 6 z5 =ðρSÞ C C 5B C B nx λ0 6 z2 2 ny λ0 6 z1 1 nz λ1 6 z4 2 nz λ2 6 z5 =ðρSÞ C B A @ nx λ0 6 z1 1 ny λ0 6 z2 1 nz λ0 6 z3 6 6 1 ðγ 2 1Þ ðβU 2 λ2 Þλ1 z4 1 ðλ1 2 βUÞλ2 z5 =ðγρSÞ
where
ð2:56Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S 5 U 2 ðβ21Þ2 1 4βc2 :
REFERENCES [1] N.A. Pierce and M.B. Giles, Preconditioning Compressible Flow Calculations on Stretched Meshes, AIAA Paper 96-0889, 1996. [2] P. Moinier, M.B. Giles, Stability analysis of preconditioned approximations of the Euler equations on unstructured meshes, J. Comput. Phys. 178 (2002) 498519. [3] P. Moinier and M.B. Giles. Preconditioned Euler and Navier-Stokes Calculations on Unstructured Grids, 6th ICFD Conference on Numerical Methods for Fluid Dynamics, Oxford, UK, 1998. [4] S. Venkateswaran, D. Li and C.L. Merkle, Influence of Stagnation Regions on Preconditioned Solutions at Low Speeds, AIAA-2003-0435, 41st Aerospace Sciences Meeting and Exhibit, Reno, NV, January 69, 2003. [5] S. Venkateswaran, and C.L. Merkle, Dual-Time stepping and Preconditioning for Unsteady Computations, AIAA Paper 95-0078, 33rd Aerospace Sciences Meeting and Exhibit, January, 912, 1995.
30
Computational Fluid-Structure Interaction
[6] S. Venkateswaran, P.E.O. Buelow and C.L. Merkle, Development of Linearized Preconditioning Methods for Enhancing Robustness and Efficiency of Euler and NavierStokes Computations, AIAA Paper 97-2030, 1997. [7] S. Venkateswaran, C.L. Merkle, Analysis of Time-Derivative Preconditioning for the Navier-Stokes Equations, 6th International Symposium on Computational Fluid Dynamics, pp. 13231328, 1995. [8] E. Turkel, Preconditioned methods for solving the incompressible and low speed compressible equations, J. Comput. Phys. 72 (1987). [9] B. Van Leer, W.T. Lee, and P.L. Roe, Characteristic Time-Stepping or Local Preconditioning of the Euler Equations, AIAA Paper 91-1552, 1991. [10] J.M. Weiss, W.A. Smith, Preconditioning applied to variable and constant density flows, AIAA J. 33 (11) (1995) 2050. [11] P.L. Roe, Approximate Riemann solvers, parameter vectors and difference schemes, J. Comput. Phys. 43 (1981) 357372. [12] Jae Wook Kim and Duck Joo Lee, Generalized Formulation and Application of Characteristic Boundary Conditions, AIAA Paper, 98-222. [13] J.W. Kim, D.J. Lee, Implementation of boundary conditions for optimized high-order compact schemes, J. Comput. Acous. 5 (2) (1997) 177191. [14] K.W. Thompson, Time dependent boundary conditions for hyperbolic systems, J. Comput. Phys. 68 (1987) 124. [15] K.W. Thompson, Time dependent boundary conditions for hyperbolic systems II, J. Comput. Phys. 89 (1990) 439461. [16] T.J. Poinsot, S.K. Lele, Boundary conditions for direct simulations of compressible viscous flows, J. Comput. Phys. 101 (1992) 104129. [17] J.L. Steger, R.F. Warming, Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods, J. Comput. Phys. 40 (1981) 263293.
FURTHER READING Joseph L. Steger and R. F. Warming, Flux Vector Splitting of the Inviscid Gasdynamic Equations with Application to Finite DifferenceMethods, NASA contractor report, HC A04/ MFA01, N79-28950.
APPENDIX PRECONDITIONING MATRIX DEFINITION 0
1=β 0 γT 0 B u=β 0 γT 0 B 0 0 Γ1 5 B B v=β 0 γT 0 @ w=β γT H=β 0 γT 0 2 1
0 ρ 0 0 ρu
0 0 ρ 0 ρv
0 0 0 ρ ρw
1 2 ρ=T 0 C 2 ρu=T 0 C 0 C 2 ρv=T C 0 A 2 ρw=T 0 ρ½γ=ðγ 2 1Þ 2 H=T
where: γ 5 Cp =Cv c2 5 γp=ρ T 0 5 p=ρ 5 c2 =γ H 5 ðρet 1 pÞ=ρ β 0 5 β=½1 1 ðγ 2 1Þβ
Mathematical Formulation Chapter | 2
31
Noted that Γ1 is a rank one modification of M (which is defined in the main body) allows one to easily compute the matrix products Γ21 1 M and M 21 Γ1 (which will be needed) via the ShermanMorrisonWoodbury formula. The general preconditioned Jacobian matrix is: 21 @Fc H~ p 5 Γ21 1 H p 5 Γ1 @W p 2 3 βU ρβγT 0 nx ρβγT 0 ny ρβγT 0 nz 0 6 nx =ρ U 0 0 07 6 7 6 ny =ρ 0 U 0 07 56 7 4 nz =ρ 0 0 U 05 Uðγ 2 1Þðβ 2 1Þ=ðργÞ ðγ 2 1ÞβT 0 nx ðγ 2 1ÞβT 0 ny ðγ 2 1ÞβT 0 nz U
where U 5 nx u 1 ny v 1 nz w, with nx ; ny ; nz being the x, y, and z components of unit vector n. The eigenvalues of H~ p are: ðβ 1 1ÞU 6 S ~ λðH p Þ 5 λ0 5 U; λ0 5 U; λ0 5 U; λ1;2 5 2 where: S5
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U 2 ðβ21Þ2 1 4βc2
The right eigenvectors of H~ p are: 3 2 ðλ1 2βUÞ=S 0 0 0 ðβU 2λ2 Þ=S 7 6 0 2n n nx =ðρSÞ 2nx =ðρSÞ 7 6 z y 7 6 7 X H~ p ;R 5 6 n 0 2n n =ðρSÞ 2n =ðρSÞ z x y y 7 6 7 6 0 nz =ðρSÞ 2nz =ðρSÞ 5 4 2ny nx nx ny nz ðγ 2 1ÞðβU 2 λ2 Þ=ðγρSÞ ðγ 2 1Þðλ1 2βUÞ=ðγρSÞ And the corresponding left eigenvectors matrix is given by: X H~ p ;L 3 0 nz 2ny nx 2ðγ 2 1Þnx =ðγρÞ 6 2ðγ 2 1Þn =ðγρÞ 2nz 0 nx ny 7 7 6 y 7 6 7 2ðγ 2 1Þn =ðγρÞ n 2n 0 n 56 z y x z 7 6 7 6 1 ρðλ1 2 βUÞnx ρðλ1 2 βUÞny ρðλ1 2 βUÞnz 0 5 4 1 ρðλ2 2 βUÞnx ρðλ2 2 βUÞny ρðλ2 2 βUÞnz 0 2
Mathematical Formulation for Preconditioned Compressible Flow Solver 2.1 BASIC GOVERNING EQUATIONS FOR FLUID FLOWS The 3D unsteady NavierStokes equations may be written in a number of forms. One common form is as follows: Continuity: @ρ @ðρuÞ @ðρvÞ @ðρwÞ 1 1 1 50 @t @x @y @z
ð2:1Þ
u-Momentum: @ðρuÞ @ðρu2 1 pÞ @ðρuvÞ @ðρuwÞ @τ xx @τ xy @τ xz 1 1 1 5 1 1 @t @x @y @z @x @y @z
ð2:2Þ
v-Momentum: @ðρvÞ @ðρuvÞ @ðρv2 1 pÞ @ðρvwÞ @τ yx @τ yy @τ yz 1 1 1 5 1 1 @t @x @y @z @x @y @z
ð2:3Þ
w-Momentum: @ðρwÞ @ðρuwÞ @ðρvwÞ @ðρw2 1 pÞ @τ zx @τ zy @τ zz 1 1 1 5 1 1 @t @x @y @z @x @y @z
ð2:4Þ
Energy equation: @ðρet Þ @ðρet uÞ @ðρet vÞ ðρet wÞ 1 1 1 @t @x @y @z 5
@ðκð@T=@xÞÞ @ðκð@T=@yÞÞ @ðκð@T=@zÞÞ @ðupÞ @ðvpÞ @ðwpÞ 1 1 2 2 2 @x @y @z @x @y @z
1
@ðuτ xx 1 vτ xy 1 wτ xz Þ @ðuτ yx 1 vτ yy 1 wτ yz Þ @ðuτ zx 1 vτ zy 1 wτ zz Þ 1 1 @x @y @z ð2:5Þ
Computational Fluid-Structure Interaction. DOI: https://doi.org/10.1016/B978-0-12-814770-2.00002-7 © 2019 Elsevier Inc. All rights reserved.
11
12
Computational Fluid-Structure Interaction
where the volumetric heating, such as absorption or emission of radiation and body source terms, such as gravity terms are not included. ρ is density; u, v, and w are the Cartesian components of velocity along x-, y-, and z-axes; p is pressure; and T is temperature. They can be related by the state equation of perfect gas p 5 ρRT, where R is the gas constant, which is 287 m2 =ðs2 KÞ for air at standard conditions. Furthermore, et is the total energy defined as et 5 e 1 ð1=2Þðu2 1 v2 1 w2 Þ, where e is internal energy per unit mass of the fluid and e 5 Cv T, where Cv is the specific heat at constant volume and Cv 5 R=ðγ 2 1Þ. γ is the ratio of specific heats and γ 5 1:4. If we assume that the fluid considered is a calorically perfect gas, then the NavierStokes equations are closed by the state equation for perfect gas: p 5 ðγ 2 1Þρe T5
ðγ 2 1Þe R
H 5e1
p 1 1 ðu2 1 v2 1 w2 Þ ρ 2
Finally, the viscous stresses are related to relations: @u 2 @u @v @w 1 1 τ xx 5 2μ 2 μ ; @x 3 @x @y @z @v 2 @u @v @w 1 1 τ yy 5 2μ 2 μ ; @y 3 @x @y @z @w 2 @u @v @w 2 μ 1 1 τ zz 5 2μ ; @z 3 @x @y @z
the velocity field by the Stokes @v @u 1 τ xy 5 τ yx 5 μ @x @y @u @w 1 τ xz 5 τ zx 5 μ @z @x @w @v 1 τ yz 5 τ zy 5 μ @y @z
The dynamic viscosity μ and conductivity κ are properties of the fluid and are functions of temperature. These two quantities are related by the Prandtl number: Pr 5
μCp κ
For air, the Prandtl number is around 0.72 at room temperature. Equations of fluid motion may be nondimensionalized to achieve certain objectives. First, it would provide conditions upon which dynamic and energetic similarity may be obtained for geometrically similar situations. Second, the solution of such equations would usually provide values within limits between zero and one, which can enhance accuracy for numerical computation.
Mathematical Formulation Chapter | 2
13
The nondimensional variables used in this study are defined as follows: x y z u v w ðx ; y ; z Þ 5 ; ; ; ; ðu ; v ; w Þ 5 L L L UN UN UN t 5 ρ 5
t p 2 pin et H ; p 5 ; et 5 ; H 5 2 2 L=U N ρN ðU N Þ ðU N Þ ðU N Þ2
ρ T μ ρN U N L γμ Cp r T ; γ5 ; T5 ; μ 5 ; Re 5 ; q 52 N 2 ρN μN μN Pr Cv ðU N Þ =Cv
where pin is the inlet or reference pressure; L is the characteristic length of the computed model; U N is the inflow velocity; and μ is the dynamic viscosity of the fluid. The variables with a superscript asterisk sign ( ) are nondimensional parameters and the asterisk sign will be dropped in subsequent equations for sake of convenience. It should be noted that we subtract a constant value (the reference pressure) from the pressure term to control the round-off errors for low-speed flows, which is found to be critical in controlling computational errors in the momentum equations for low-speed compressible flows. The use of gauge pressure is a common treatment for incompressible solvers because only pressure gradients are needed for all calculations. For compressible solvers, the absolute value of pressure must be used when dealing with the energy equation and state equation of gases. But when we use nondimensional absolute pressure at low Mach numbers, it becomes extremely large although the pressure gradients in momentum equations are small. The use of gauge pressure can avoid performing the addition and subtraction operations between two extraordinarily large values of nondimensional absolute pressures. In our experience, the result obtained with gauge pressure is more accurate than without using it. And the relative error can be up to 4%5% of the result. It should be noted that the original absolute pressure (designated by p0 hereafter) should be used instead of gauge pressure whenever the state equation of gas is involved. The nondimensional equations of fluid motion may be expressed in vector form as @W c @Ei @Fi @Gi @Ev @Fv @Gv 1 1 1 5 1 1 @t @x @y @z @x @y @z 2 3 ρ 6 ρu 7 6 7 7 Wc 5 6 6 ρv 7 4 ρw 5 ρet
ð2:6Þ
ð2:7Þ
14
Computational Fluid-Structure Interaction
2
3 ρu 6 ρu2 1 p 7 6 7 7 Ei 5 6 6 ρuv 7 4 ρuw 5 ðρet 1 pÞu 3 2 0 7 6 τ xx 7 1 6 7 6 τ xy Ev 5 7 6 ReN 4 τ 5 xz uτ xx 1 vτ xy 1 wτ xz 2 qx 2 3 ρv 6 ρvu 7 6 2 7 7 Fi 5 6 ρv 1 p 6 7 4 ρvw 5 ðρet 1 pÞv 3 2 0 7 6 τ yx 7 1 6 7 6 τ yy Fv 5 7 6 ReN 4 τ 5 yz uτ yx 1 vτ yy 1 wτ yz 2 qy 3 2 ρw 7 6 ρwu 7 6 7 6 Gi 5 6 ρwv 7 4 ρw2 1 p 5 ðρet 1 pÞw 2 3 0 6 τ zx 7 7 1 6 6 7 τ zy Gv 5 7 ReN 6 4 τ zz 5 uτ zx 1 vτ zy 1 wτ zz 2 qz
ð2:8Þ
ð2:9Þ
ð2:10Þ
ð2:11Þ
ð2:12Þ
ð2:13Þ
2.2 LOW-SPEED PRECONDITIONING FORMULATION One difficulty with compressible NavierStokes solvers is their slow convergence rates and even unstable solutions for low-Mach-number flows. This difficulty can be traced to a disparity between the acoustic and convective speeds [110], and can be addressed by a preconditioning algorithm. Previous work in this area has been reported by Venkateswaran and
Mathematical Formulation Chapter | 2
15
Merkle [47], Turkel [8], Van Leer et al. [9], and Weiss and Smith [10]. The applications of the preconditioning methods have been found in the computation of steady flows without considering arbitrarily moving objects in the flow field. The preconditioned NavierStokes equations for 3D compressible unsteady flows can be given in vector form, explicitly expressing the conservation laws of mass, momentum, and energy. We also introduce, in the equations, pseudo time terms to provide pseudo time marching for their numerical solutions: Γ1 where
@W p @W c 1 1 rUFi 5 rUFv @τ @t 3 ρ 6 ρu 7 6 7 7 Wc 5 6 6 ρv 7 4 ρw 5 ρet 2 3 p 6u 7 6 7 7 Wp 5 6 6v 7 4w 5 T0 3 2 ρU 6 ρuU 1 pi 7 7 6 7 Fi 5 6 6 ρvU 1 pj 7 4 ρwU 1 pk 5 ρHU 2 3 0 6 τx 7 7 1 6 6 τy 7 Fv 5 6 7 ReN 4 τ 5 z τ UU 2 q
ð2:14Þ
2
ð2:15Þ
ð2:16Þ
ð2:17Þ
ð2:18Þ
τ is the pseudo time and Γ1 is the preconditioning matrix in the pseudo time terms for low-Mach-number flows, which is defined in the appendix. W c and W p are the vectors of conservative and primitive dependent variables, respectively; Fi and Fv are the inviscid convective flux and viscous flux vectors. Furthermore, we have the following formulas: U 5 ui 1 vj 1 wk
16
Computational Fluid-Structure Interaction
τ 5 τ xi 1 τ yj 1 τ zk τ i 5 τ ix i 1 τ iy j 1 τ iz k @ui @uj 2 τ ij 5 μ 1 2 δij rUU 3 @xj @xi q 5 qx i 1 qy j 1 qz k T 0 5 p0 =ρ 5 c2 =γ i, j, and k are the three unit vectors in three Cartesian directions; τ ix , τ iy , and τ iz are the viscous stresses.
2.3 DISCRETIZATION METHODS The 3D equations (2.14) is transformed into an integral form and discretized on an unstructured grid. A cell-vertex finite-volume scheme is adopted here. For every vertex, a control volume is constructed using the median duals of the tetrahedral cells by connecting all the neighboring centroids of edges, surfaces, and tetrahedra. As shown in Fig. 2.1 is a part of the control volume associated with vertex P within tetrahedron PABC. Spatial discretization is performed by using the integral form of the conservation equations over the control volume surrounding node P: ððð ððð ððð ððð @Q0 1 @W c dV 1 dV 1 rUFi dV 2 rUFv dV 5 0 ð2:19Þ cv @τ cv @t cv cv
FIG. 2.1 Construction of control volume within a tetrahedron for a node P.
Mathematical Formulation Chapter | 2
17
Noted that a new variable Q0 1 has arisen as @Q0 1 =ð@τÞ 5 Γ1 @W p =ð@τÞ, and the Jacobian Γ1 5 @Q0 1 =ð@W p Þ. So that we have: @Q0 1 @Q0 1 @W p @W p 5 5 Γ1 @τ @W p @τ @τ
2.3.1 Edge-Based Method and Cell-Based Method The convective term in Eq. (2.19) is transformed into a summation: ððð ð nbseg X rUFi dV 5 Fi UndS 5 ½ðFi Þij UnΔSk cv
Scv
ð2:20Þ
k51
where nbseg is the number of the edges associated with node P; ðFi Þkij is the inviscid flux through the part of control volume surface associated with edge k; and n is the unit normal vector of the control volume surface. Finally, ΔSk is a part of the control volume surface associated with edge k (as shown in Fig. 2.1, if k is edge PC, then 1O2c is ΔSk ). Therefore all the fluxes are calculated for the edges and then collected at the two ends of each edge for updating of flow variables in time marching, which is the so-called edgebased data structure for fast inviscid flux computation and efficient data storage and retrieval. On the other hand, the viscous term is calculated using a cell-based method for the same reason of high computational efficiency: ððð ð ncell X ½Fv UnΔSc i rUFv dV 5 Fv UndS 5 ð2:21Þ cv
Scv
i51
where ncell is the number of tetrahedral elements associated with node P and ΔSci is the part of control volume surface in cell i. By using the following relation: ð dS 5 0: Scv
The total vector surface of the control volume in a cell i becomes nΔSci 5 ð1=3ÞðnΔSpi Þ. Thus the calculation of viscous terms can be simplified as ncell X i51
½Fv UnΔSc i 5
ncell 1X ½Fv UnΔSp i 3 i51
ð2:22Þ
where nΔSpi is the surface vector of the face opposite node P of the tetrahedron under consideration, where the ðFv Þi is calculated at the center of the tetrahedron with a node P, and can be obtained by using the Green’s theorem based on the variables at the four vertices of the tetrahedron. Similar to the
18
Computational Fluid-Structure Interaction
Galerkin type of formulation, the gradient of a flow variable φ at the center of a tetrahedron is evaluated as follows: P4 P 1 4i51 φi Si i51 φi 9Si 52 ð2:23Þ gradφc 5 2 3 27V V where φi is the flow variable at a vertex i of the tetrahedron; Si is the surface area that is opposite to node i; and V is the volume of the tetrahedron. Gradients at the vertices are obtained by a volume averaging of the gradients at the centers of cells associated with the vertex under consideration.
2.3.2 Roe’s TVD Scheme for Inviscid Flux in Edge-Based Method The Riemann problem is an initial problem with piecewise constant initial data. Its exact solution represents the real physical characteristics of a flow with several families of waves and their propagation. When the wave is a shock, the equations for this initial problem cannot be explicitly solved and an iterative solution has to be used. In order to reduce the amount of calculation, some approximate treatments have been employed in numerical solutions. In this work, a high-order Roe’s TVD scheme for compressible flow for arbitrary unstructured 3D grids has been adopted. In Roe’s approach [11], averaged values of density, velocity, and enthalpy at the interface between the two initial states at two neighboring points are calculated by using a special averaging procedure, called Roe’s averaging: pffiffiffiffiffipffiffiffiffiffi ρ 5 ρL ρR pffiffiffiffiffi pffiffiffiffiffi ðu; v; wÞL ρL 1 ðu; v; wÞR ρR pffiffiffiffiffi pffiffiffiffiffi ðu; v; wÞ 5 ð2:24Þ ρL 1 ρR pffiffiffiffiffi pffiffiffiffiffi H L ρL 1 H R ρR p p ffiffiffiffiffi ffiffiffiffiffi H5 ð2:25Þ ρL 1 ρR where L and R represent the two neighboring points on an edge, and H denotes total enthalpy which has the following relation: H5e1
p0 1 1 u2 1 v 2 1 w 2 2 ρ
The speed of sound can be derived as: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 2 c5 H 2 ðu 1 v 1 w Þ ðγ 2 1Þ 2
ð2:26Þ
ð2:27Þ
Mathematical Formulation Chapter | 2
19
The inviscid flux ðFi Þij is evaluated based on Roe’s approximate Riemann solver (or the flux difference splitting scheme): 1
ðFi Þij 5 FLi 1 FRi 2 A ðW R 2 W L Þ 2 1 L 5 Fi 1 FRi 2 jΔFj ð2:28Þ 2 5 FLi 1 ΔF2 5 FRi 2 ΔF1 where
FLi 5 Fi W Lc FRi 5 Fi W Rc
and W Lc and W Rc are the left and right conserved state vectors on the two neighboring points of an edge. The flux difference is: ð2:29Þ jΔFj 5 ΔF1 2 ΔF2 8 3 39 2 2 1 0 > > > > > > > 7 7> 6 6 > > u ΔU Δu2n > > 7 7 6 6 x n = < 7 7 6 6 Δp 6 6 7 7 6 6 1ρ ΔF 5λ1 Δρ2 2 6 v Δv2n ΔU y n 7 7> 6 > ðcÞ 6 > 7 7> 6 > > > 5 5> 4w 4 Δw2nz ΔUn > > > > ; : 2 ðU n Þ =2 uΔu1vΔv1wΔw2U n ΔUn 3 2 1 7 6 6 u1nx c 7 7 6 Δp1ρcΔU n 6 v1ny c 7 1λ46 7 6 2 2ðcÞ 7 6 4 w1nz c 5 H 1U n c 3 2 1 7 6 6 u6nx c 7 7 6 Δp2ρcΔUn 6 1λ56 v6ny c 7 2 7 6 2ðcÞ 7 6 4 w6nz c 5 H 6U n c ð2:30Þ where ðU n Þ2 5 u2 1 v2 1 w2 ΔUn 5 Δunx 1 Δvny 1 Δwnz U n 5 unx 1 vny 1 wnz
20
Computational Fluid-Structure Interaction
Furthermore, the eigenvalues of A are λ1 5 U n , λ4;5 5 U n 6 c. Because of the introduction of preconditioning matrix Γ1 , the inviscid fluxes, ðFi Þijk , through the face k is now reformulated as 1 1 @Fi 0 k ðFi Þij ððFi Þi 1ðFi Þj Þk 2 @Q0 ðδQ1 Þk 2 2 1 k ! ð2:31Þ 1 1 @Fi @W p @Q01 5 ððFi Þi 1ðFi Þj Þk 2 @W p @Q0 @W p ððW p Þj 2ðW p Þi Þk 2 2 1 k
k
Note that we have retained the variable, Q0 1 , in computing this flux. Defining the Jacobian in the normal direction as @Fi ðHp Þk 5 @W p k And using the previously defined Jacobian Γ1 5 @Q0 1 =@W p , then the above expression becomes 1 1 ðFi Þkij ððFi Þi 1ðFi Þj Þk 2 Hp Γ21 1 k Γ1k ððW p Þj 2ðW p Þi Þk 2 2
Drop the subscript k in the flux vector and the Jacobian with the assumption that the fluxes and Jacobians all correspond to conditions in the normal direction on the given control volume surface. And after some simple algebraic derivations, we have 1 1 ðFi Þij ððFi Þi 1 ðFi Þj Þ 2 Γ1 Γ21 1 Hp ððW p Þj 2 ðW p Þi Þ 2 2
ð2:32Þ
2.3.3 Upwind-Biased Interpolation Combined with the third-order Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) interpolation, it can produce accurate and stable solution on unstructured grids. The left and right state vectors W L and W R at a control volume surface are evaluated using a nominally third-order upwind-biased interpolation scheme: 1
1 ð1 2 κÞΔ2 i 1 ð1 1 κÞΔi 4 i 1h 2 W R 5 W j 2 ð1 2 κÞΔ1 j 1 ð1 1 κÞΔj 4 WL 5 Wi 1
where 2 Δ1 i 5 Δj 5 W j 2 W i 1 Δ2 i 5 W i 2 W i21 5 2ijUrW i 2 ðW j 2 W i Þ 5 2ijUrW i 2 Δi 2 Δ1 j 5 W j11 2 W j 5 2ijUrW j 2 ðW j 2 W i Þ 5 2ijUrW j 2 Δj
ð2:33aÞ ð2:33bÞ
Mathematical Formulation Chapter | 2
21
Therefore substituting the above equations into Eqs. (2.33a) and (2.33b), the final equations based on upwind-biased interpolation scheme is shown as follows: 1
ð1 2 κÞijUrW i 1 κΔ1 i 2 i 1h W R 5 W j 2 ð1 2 κÞijUrW j 1 κΔ2 j 2 WL 5 Wi 1
ð2:34aÞ ð2:34bÞ
where κ is set to 1/3, which corresponds to a nominally third-order accuracy. ij is the vector representing the edge, which points from node P to its neighboring node under consideration. The gradients of W at i and j are calculated by volume averaging the gradients of the cells that surround i and j. Finally, for a given node P, the spatially discretized Eq. (2.19) form a system of coupled ordinary differential equations. This can be reformulated as: @Q0 1 @W c ΔVcv 1 ΔVcv @τ @t 8 9 ncell <1 nbseg = X X 1 52 ½ððFi Þi1ðFi Þj Þ2Γ1 jΓ21 ½Fv UnΔSP i 1 Hp jððW p Þj2ðW p Þi Þk UðnΔSÞk 2 :2 k51 ; 3 i51 52R ð2:35Þ
where R represents the residual, which includes the convective and diffusive fluxes and ΔVcv is the control volume of node P. the overbar in Eq. (2.35) denotes the cell-averaging value. An implicit scheme is adopted for Eq. (2.35) and the time-dependent term is discretized using a second-orderaccurate backward differencing scheme: n11 n11 n n21 n21 @Q0 1 1:5ΔVcv W c 2 2:0ΔVcv W nc 1 0:5ΔVcv Wc n11 ~ n11 5R ΔVcv 5 2 R 2 @τ Δt ð2:36Þ where the superscript (n 1 1) denotes the physical time level (n 1 1)Δt and ~ n11 Þ is the new modiall the variables are evaluated at this time level, RðW p fied residual which contains both the time derivative and flux vectors. The derivative with respect to a pseudo time τ is discretized as n11 ΔVcv Γ1
W n11;m11 2 W n11;m n11;m11 p p 5 R~ dτ
ð2:37Þ
whose solution is sought by marching to a pseudo steady state in τ, where m and (m 1 1) denote the initial and final pseudo time levels. Once the artificial steady state is reached, the derivative of Wp with respect to τ becomes zero, ~ n11 5 0. Hence, the original unsteady and the solution satisfies R
22
Computational Fluid-Structure Interaction
NavierStokes equations are fully recovered. Therefore instead of solving the equations in each time-step in the physical time domain (t), the problem is transformed into a sequence of steady-state computations in the artificial time domain ðτÞ. Eq. (2.37) can be integrated in pseudo time by an explicit five-stage RungeKutta scheme. However the pseudo time-step size may be severely restricted if the physical time-step size is very small. Hence, a fully implicit dual-time-stepping is adopted here.
2.3.4 Matrix-Free, Implicit, Dual-Time-Stepping Algorithm A Taylor series expansion is performed for the residual in Eq. (2.37) with respect to the pseudo time for node i, m11 ~m 1 R~ i 5 R i
nbseg X @R ~i ~i @R ΔðW p Þi 1 ΔðW p Þj @ðW p Þi @ðW p Þj j51
ð2:38Þ
where nbseg is the number of edges connected to i, which is also equal to the number of neighboring points connected to point i through the edges. An approximate flux function is introduced here to simplify the implicit timestepping calculation. The total flux (including both convective and viscous fluxes) across a control volume surface associated with a certain edge ij can be approximated as Fij
1
ðFi Þi Un 1 ðFi Þj Un 2 jλij jððW p Þj 2 ðW p Þi Þ 2
where λij is the spectral radius based on the preconditioned system which is associated with edge ij, and is given in the appendix. Then in all of the R(W) terms (Eq. (2.35)), we have " #
nbseg nbseg X 1 @ðFi Þij X 1 @ðFi Þij @Ri @Ri 5 1 jλij j and 5 2 jλij j : 2 @ðW p Þi 2 @ðW p Þj @ðW p Þi @ðW p Þj j51 i51 And for the physical time-dependent terms, we have ðW c Þm11 5 ðW c Þm i i 1
@W c ΔðW p Þi @W p
After combining all the residual terms at every node in the flow field into a vector and dropping the third term of the right-hand side of Eq. (2.38) as an approximation, we have n11 ~ n11;m11 5 R ~ n11;m 2 @Ri ΔðW p Þi 2 1:5 ΔVcv @W c ΔðW p Þi R i i @ðW p Þi Δt @W p n11;m 5 R~ i 2
nbseg X j51
ðHp;j ÞΔðW p Þi 2 1:5
n11 ΔVcv MΔðW p Þi Δt
Mathematical Formulation Chapter | 2
23
where Hp;j 5 ð1=2Þ½ð@ðFi Þij =@ðW p Þi Þ 1 jλij j, M 5 @W c =@W p . And the wholefield equivalent of Eq. (2.36) can then be rewritten as 0 1 nbseg X W n11;m11 2 W n11;m 1:5Δτ Δτ p 21 n11 @I 1 A p Γ21 ΔVcv ðΓ H Þ p;j 1 M1 1 n11 Δt ΔVcv Δτ j51 0 1 n11;m n n21 n11 n n21 1:5W ΔV 2 2:0W ΔV 1 0:5W ΔV cv cv cv c c c @2 2 Rn11;m A 5 Γ21 1 Δt ð2:39Þ that is, A~
W n11;m11 2 W n11;m p p n11 ~ n11;m ΔVcv 5 Γ21 1 R Δτ
thus W n11;m11 2 W n11;m p p ~~ n11;m n11 ΔVcv 5R Δτ where ~~ R
n11;m
X 21 ~ n11;m and A~ 5 I 1 1:5Δτ Γ121 M 1 Δτ R 5 A~ Γ21 Γ21 1 1 Hp;j n11 Δt ΔVcv j51 nbseg
Therefore n21 n11 n11;m 2 2:0ΔSncv W n 1 0:5ΔSn21 cv W ~Rn11;m 5 2 Rn11;m 2 K 1:5ΔScv W Δt Further approximation can be introduced in order to achieve matrix-free computation. If we employ point implicit treatment to the preceding equations, then only the diagonal terms in A~ are used in the pseudo timestepping. As a result, the equation for every node can now be written as W n11;m11 2 W n11;m p p ~~ n11;m n11 ΔVcv 5R Δτ where ~~ n11;m 5 A~21 Γ 21 R ~ n11;m R 1 ii i i
ð2:40Þ
24
Computational Fluid-Structure Interaction
and
2
21 A~ii
nbseg 1:5Δτ 21 Δτ X 21 Γ1 M1 5 diag4 I1 Γ1 Hp;j n11 Δt ΔVcv j51
!21 3 5:
Pseudo time-stepping is then performed on Eq. (2.40).
2.3.5 Time Integration In this work, a RungeKutta time-integration algorithm is used between each physical time-step to iterate the numerical solution in an artificial time τ until convergence is reached. For a five-stage RungeKutta scheme, the stage coefficients are α1 5
1 1 3 1 ; α2 5 ; α3 5 ; α4 5 ; α5 5 1 4 6 8 2
and we have m W ð0Þ c;p 5 W c;p
Δτ ~~ ð0Þ R W c;p ΔVcv Δτ ~~ ð1Þ ð0Þ R W c;p W ð2Þ c;p 5 W c;p 2 α2 C ΔVcv Δτ ~~ ð2Þ ð0Þ R W c;p W ð3Þ c;p 5 W c;p 2 α3 C ΔVcv Δτ ~~ ð3Þ ð0Þ R Wc;p W ð4Þ c;p 5 W c;p 2 α4 C ΔVcv Δτ ~~ ð4Þ ð0Þ R Wc;p 5 W 2 α C W ð5Þ 5 c;p c;p ΔVcv
ð0Þ W ð1Þ c;p 5 W c;p 2 α1 C
ð2:41Þ
W ðm11Þ 5 W ð5Þ c;p c;p
2.4 CONVERGENCE ACCELERATION TECHNIQUES To speed up the convergence rate, apart from the matrix-free, implicit, dualtime-stepping algorithm and local time-stepping scheme introduced in the previous sections, an implicit residual smoothing scheme developed for unstructured grids is employed. The smoothing equation for a vertex k can be expressed as follows: Rk 5 Rk 1 ε r2 Rk
ð2:42Þ
Mathematical Formulation Chapter | 2
25
where R is the original residual; R is smoothed residual; and ε is the smoothing coefficient, which can be defined as ( " # ) 1 CFL 2 21 ; 0 ε 5 max 4 CFL where CFL is the maximum CFL number of the basic scheme. The solution to the above equations can be obtained on an unstructured grid by using the Jacobi iterative method as follows:
ðmÞ
Rk 5
Rð0Þ k 1ε
numnodðkÞ P i51
ðm21;mÞ
Ri
1 1 ε UnumnodðkÞ
ð2:43Þ
where numnodðkÞ is the number of neighboring nodes of vertex k. Furthermore, the multigrid method can also be employed for convergence acceleration, which will be introduced and discussed in details in Chapter 5, The Multigrid Method.
2.5 BOUNDARY CONDITIONS Physically correct numerical boundary conditions are used to achieve timeaccurate solutions for both steady and unsteady flows. This kind of boundary conditions should be compatible with modern low-dissipation and lowdispersion numerical algorithms, which have very low-dispersion errors and require precise boundary conditions to avoid numerical instabilities and to control spurious nonphysical wave reflections at the computational boundaries. A high-order and high-resolution scheme as one of those numerical algorithms gives results sensitive to the boundary values in the computation, thus the quality of solutions even on interior nodes are dependent on the accuracy and correctness of boundary conditions [12]. Thompson [13,14] decomposed hyperbolic equations into wave modes of definite velocity and then specified characteristic boundary conditions for incoming waves. The starting point of his analysis was nonlinear Euler equations. The idea of his approach was that 1D characteristic analysis could be performed by consideration of the transverse terms as a constant source term. The amplitudes of outward propagating waves are defined entirely from the variables inside computational domain, while those of inward propagating waves are specified as the characteristic boundary conditions. In our studies, we choose to let the state variables of the far-field boundary be floating rather than fixed; in other word, we solve the characteristic equations on the boundary nodes. This method can make the boundary treatment nonreflective to the maximum extent and our solver more stable [15].
26
Computational Fluid-Structure Interaction
2.5.1 Far-Field Conditions In a cell-centered scheme, the Roe TVD upwind scheme described above can be extended to the boundaries conveniently. Since the flux-formula automatically selects the proper information from either side of the boundary face. One may add in an extra cell across the boundary in which the full state is prescribed at any time, even for subsonic inflow and outflow; while in our cellvertex scheme, the far-field boundary conditions for the flow equations are based on the flux vector splitting technique for the inviscid fluxes [16]. For a boundary node k with its control volume boundary face ΔSk , the inviscid flux across ΔSk is evaluated by using the Stegger and Warming splitting [16]. At the point upstream of the inflow boundary (denoted by “in”), ρin ; uin ; vin ; win and free stream Mach number MN are prescribed. Then, the inflow pressure, temperature, and total energy are determined by: 1 γðMN Þ2
ð2:44Þ
1 γðγ 2 1ÞðMN Þ2
ð2:45Þ
pin 5 Tin 5 ðρet Þin 5
1 2 pin ρin uin 1 v2in 1 w2in 1 2 γ21
ð2:46Þ
At the point downstream of the outflow boundary, the only one known state is the pressure. All of the other states must be extrapolated from interior flow field using characteristic variables [9,15]. For an arbitrary boundary normal direction n, we can write the variation of characteristic variables as 0
@ρ 2
1 @p c2
1
C 1 B C B @W 1λn C B sU@U C B B @W 2 C B C λn C B tU@U C B 3 C 5 @W λn 5 B @W 1 C B λn C B C B @p nU@U 1 @ @W 4 A B C ρc λn C B @W 5λn C B 1 A @ 2 nU@U 1 @p ρc 0
ð2:47Þ
The first characteristic variation is proportional to entropy variations, the second and third are related to variations in shear and tangent velocity, and the last two represent acoustic disturbances. Thus from @W 1λn 1=γ pout ρout 5 ρb ð2:48Þ pb
Mathematical Formulation Chapter | 2
27
where subscript “b” denotes the variable from the boundary node and “out” denotes flow state at the point downstream of the outflow boundary. From @W 4λn , @W 5λn , we have ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi γpout =ρout b 2 2 2 2 2 2 Un;out 5 un;out 1 vn;out 1 wn;out 5 un;b 1 vn;b 1 wn;b 1 γ21 From @W 3λn U t;out 5 U t;b 5 U b 2 Un;b n Therefore we have Uout 5 U t;out 1 Un;out n 5 U b 1 ðUn;out 2 Un;b Þn
ð2:49Þ
uout 5 ub 1 ðUn;out 2 Un;b Þnx
ð2:50Þ
vout 5 vb 1 ðUn;out 2 Un;b Þny
ð2:51Þ
wout 5 wb 1 ðUn;out 2 Un;b Þnz
ð2:52Þ
pout ρ jU out j2 1 out γ21 2
ð2:53Þ
ðρet Þout 5
As an example, the far-field boundary conditions treatment of a strong converged nozzle is shown in Fig. 2.2. The boundary fluxes are calculated using Stegger and Warming’s flux vector splitting formula [16,17]: Fi UnSk 5 ðFi 1 1 Fi 2 ÞUΔSk
Fi 6
3 2 ρ α 7 6 2γ 7 6 7 6 ρ 7 6 ðαu 1 βcnx Þ 7 6 2γ 7 6 7 6 ρ 7 6 ðαv 1 βcny Þ 7 6 2γ 56 7 7 6 ρ 7 6 ðαw 1 βcn Þ 7 6 z 7 6 2γ 6 0 17 7 6 2 2 2 7 6 ρ 4 @α u 1 v 1 w 1 2βcU 1 ϖ c2 A 5 n 2γ γ21 2
where Un 5 UUn λ1 5 UUn
ð2:54Þ
28
Computational Fluid-Structure Interaction
FIGURE 2.2 The inflow and outflow boundary conditions for a strong converged nozzle.
λ4 5 UUn 1 c λ5 5 UUn 2 c α 5 2ðγ 2 1Þλ16 1 λ46 1 λ56 β 5 λ46 2 λ56 ϖ 5 λ46 1 λ56 λ1 5 maxfλ; 0:0g λ2 5 minfλ; 0:0g It should be noted that at the inflow, the boundary nodal states are used for calculating Fi 1 , the preimposed inflow states are used for Fi 2 . Similarly, at the outflow, the boundary nodal states are used while calculating Fi 1 , the preimposed and extrapolated outflow states are used for Fi 2 .
2.5.2 Solid Wall/Slip Conditions At the solid wall, slip (for inviscid flows) or no-slip (for viscous flows) and no-injection boundary conditions are imposed, that is, the zero normal fluxes of mass, momentum, and energy are imposed. In addition, the solid surface is assumed to be adiabatic, that is: rUT 5 0
2.5.3 Far-Field Conditions for Preconditioned Systems The far-field calculations are based on characteristic variables (Reimann invariants). Thus at inflow, the incoming variables corresponding to positive eigenvalues are specified while the outgoing variables corresponding to negative eigenvalues are extrapolated. Once we change the time-dependent equations, we also change the characteristics of the system (though not the signs
Mathematical Formulation Chapter | 2
29
of the eigenvalues). Hence, it is also necessary to modify the boundary conditions for the preconditioned system. Here, the flux at a boundary is defined as Fi UnSk 5 ðFi 1 1 Fi 2 ÞUΔSk where Fi 6 has been redefined as Fi 6 5 XHp ;R Λ 6 X Hp ;L W p , where Λ 6 5 λi 6 jλi j=ð2Þ, and λi represents the eigenvalues of Hp (see appendix for details of λi ). X Hp ;R and XHp ;L are the right and left eigenvectors of Hp (see appendix). We first calculate 0 1 ð1 2 γÞnx p=ðγρÞ 1 nz v 1 nx T 0 2 ny w B ð1 2 γÞny p=ðγρÞ 2 nz u 1 nx w 1 ny T 0 C B C 0 B C z 5 XHp ;L W p 5 B ð1 2 γÞnz p=ðγρÞ 1 ny u 2 nx v 1 nz T C @ p 1 nx ρðλ1 2 βUÞu 1 ny ρðλ1 2 βUÞv 1 nz ρðλ1 2 βUÞw A p 1 nx ρðλ2 2 βUÞu 1 ny ρðλ2 2 βUÞv 1 nz ρðλ2 2 βUÞw ð2:55Þ And then,
Fi 6 5 X Hp ;R Λ 6 XHp ;L W p 5 X Hp ;R Λ 6 z 1 0 ððβU 2 λ2 Þλ1 6 z4 1 λ1 2 βUÞλ2 6 z5 =S C B ny λ0 6 z3 2 nz λ0 6 z2 1 nx λ1 6 z4 2 nx λ2 6 z5 =ðρSÞ C B B nz λ0 6 z1 2 nx λ0 6 z3 1 ny λ1 6 z4 2 ny 6 λ2 6 z5 =ðρSÞ C C 5B C B nx λ0 6 z2 2 ny λ0 6 z1 1 nz λ1 6 z4 2 nz λ2 6 z5 =ðρSÞ C B A @ nx λ0 6 z1 1 ny λ0 6 z2 1 nz λ0 6 z3 6 6 1 ðγ 2 1Þ ðβU 2 λ2 Þλ1 z4 1 ðλ1 2 βUÞλ2 z5 =ðγρSÞ
where
ð2:56Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S 5 U 2 ðβ21Þ2 1 4βc2 :
REFERENCES [1] N.A. Pierce and M.B. Giles, Preconditioning Compressible Flow Calculations on Stretched Meshes, AIAA Paper 96-0889, 1996. [2] P. Moinier, M.B. Giles, Stability analysis of preconditioned approximations of the Euler equations on unstructured meshes, J. Comput. Phys. 178 (2002) 498519. [3] P. Moinier and M.B. Giles. Preconditioned Euler and Navier-Stokes Calculations on Unstructured Grids, 6th ICFD Conference on Numerical Methods for Fluid Dynamics, Oxford, UK, 1998. [4] S. Venkateswaran, D. Li and C.L. Merkle, Influence of Stagnation Regions on Preconditioned Solutions at Low Speeds, AIAA-2003-0435, 41st Aerospace Sciences Meeting and Exhibit, Reno, NV, January 69, 2003. [5] S. Venkateswaran, and C.L. Merkle, Dual-Time stepping and Preconditioning for Unsteady Computations, AIAA Paper 95-0078, 33rd Aerospace Sciences Meeting and Exhibit, January, 912, 1995.
30
Computational Fluid-Structure Interaction
[6] S. Venkateswaran, P.E.O. Buelow and C.L. Merkle, Development of Linearized Preconditioning Methods for Enhancing Robustness and Efficiency of Euler and NavierStokes Computations, AIAA Paper 97-2030, 1997. [7] S. Venkateswaran, C.L. Merkle, Analysis of Time-Derivative Preconditioning for the Navier-Stokes Equations, 6th International Symposium on Computational Fluid Dynamics, pp. 13231328, 1995. [8] E. Turkel, Preconditioned methods for solving the incompressible and low speed compressible equations, J. Comput. Phys. 72 (1987). [9] B. Van Leer, W.T. Lee, and P.L. Roe, Characteristic Time-Stepping or Local Preconditioning of the Euler Equations, AIAA Paper 91-1552, 1991. [10] J.M. Weiss, W.A. Smith, Preconditioning applied to variable and constant density flows, AIAA J. 33 (11) (1995) 2050. [11] P.L. Roe, Approximate Riemann solvers, parameter vectors and difference schemes, J. Comput. Phys. 43 (1981) 357372. [12] Jae Wook Kim and Duck Joo Lee, Generalized Formulation and Application of Characteristic Boundary Conditions, AIAA Paper, 98-222. [13] J.W. Kim, D.J. Lee, Implementation of boundary conditions for optimized high-order compact schemes, J. Comput. Acous. 5 (2) (1997) 177191. [14] K.W. Thompson, Time dependent boundary conditions for hyperbolic systems, J. Comput. Phys. 68 (1987) 124. [15] K.W. Thompson, Time dependent boundary conditions for hyperbolic systems II, J. Comput. Phys. 89 (1990) 439461. [16] T.J. Poinsot, S.K. Lele, Boundary conditions for direct simulations of compressible viscous flows, J. Comput. Phys. 101 (1992) 104129. [17] J.L. Steger, R.F. Warming, Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods, J. Comput. Phys. 40 (1981) 263293.
FURTHER READING Joseph L. Steger and R. F. Warming, Flux Vector Splitting of the Inviscid Gasdynamic Equations with Application to Finite DifferenceMethods, NASA contractor report, HC A04/ MFA01, N79-28950.
APPENDIX PRECONDITIONING MATRIX DEFINITION 0
1=β 0 γT 0 B u=β 0 γT 0 B 0 0 Γ1 5 B B v=β 0 γT 0 @ w=β γT H=β 0 γT 0 2 1
0 ρ 0 0 ρu
0 0 ρ 0 ρv
0 0 0 ρ ρw
1 2 ρ=T 0 C 2 ρu=T 0 C 0 C 2 ρv=T C 0 A 2 ρw=T 0 ρ½γ=ðγ 2 1Þ 2 H=T
where: γ 5 Cp =Cv c2 5 γp=ρ T 0 5 p=ρ 5 c2 =γ H 5 ðρet 1 pÞ=ρ β 0 5 β=½1 1 ðγ 2 1Þβ
Mathematical Formulation Chapter | 2
31
Noted that Γ1 is a rank one modification of M (which is defined in the main body) allows one to easily compute the matrix products Γ21 1 M and M 21 Γ1 (which will be needed) via the ShermanMorrisonWoodbury formula. The general preconditioned Jacobian matrix is: 21 @Fc H~ p 5 Γ21 1 H p 5 Γ1 @W p 2 3 βU ρβγT 0 nx ρβγT 0 ny ρβγT 0 nz 0 6 nx =ρ U 0 0 07 6 7 6 ny =ρ 0 U 0 07 56 7 4 nz =ρ 0 0 U 05 Uðγ 2 1Þðβ 2 1Þ=ðργÞ ðγ 2 1ÞβT 0 nx ðγ 2 1ÞβT 0 ny ðγ 2 1ÞβT 0 nz U
where U 5 nx u 1 ny v 1 nz w, with nx ; ny ; nz being the x, y, and z components of unit vector n. The eigenvalues of H~ p are: ðβ 1 1ÞU 6 S ~ λðH p Þ 5 λ0 5 U; λ0 5 U; λ0 5 U; λ1;2 5 2 where: S5
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U 2 ðβ21Þ2 1 4βc2
The right eigenvectors of H~ p are: 3 2 ðλ1 2βUÞ=S 0 0 0 ðβU 2λ2 Þ=S 7 6 0 2n n nx =ðρSÞ 2nx =ðρSÞ 7 6 z y 7 6 7 X H~ p ;R 5 6 n 0 2n n =ðρSÞ 2n =ðρSÞ z x y y 7 6 7 6 0 nz =ðρSÞ 2nz =ðρSÞ 5 4 2ny nx nx ny nz ðγ 2 1ÞðβU 2 λ2 Þ=ðγρSÞ ðγ 2 1Þðλ1 2βUÞ=ðγρSÞ And the corresponding left eigenvectors matrix is given by: X H~ p ;L 3 0 nz 2ny nx 2ðγ 2 1Þnx =ðγρÞ 6 2ðγ 2 1Þn =ðγρÞ 2nz 0 nx ny 7 7 6 y 7 6 7 2ðγ 2 1Þn =ðγρÞ n 2n 0 n 56 z y x z 7 6 7 6 1 ρðλ1 2 βUÞnx ρðλ1 2 βUÞny ρðλ1 2 βUÞnz 0 5 4 1 ρðλ2 2 βUÞnx ρðλ2 2 βUÞny ρðλ2 2 βUÞnz 0 2