- Email: [email protected]
Discrete conservation and coupling strategies in nonlinear aeroelasticity
Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102 www.elsevier.com/locate/cma
Discrete conservation and coupling strategies in nonlinear aeroelasticity Ralf Massjung
1
Institut fu¨r Geometrie und Praktische Mathematik, RWTH Aachen, Templergraben 55, 52056 Aachen, Germany Received 30 December 2004; accepted 30 January 2006
Abstract As a model problem of aeroelasticity the panel flutter problem is considered to study the discretization of nonlinear aeroelastic problems and the behaviour of the corresponding solution strategies. A fluid–structure discretization that mimics the energy budget of the continuous problem is derived. The set of equations occurring in each time step are implicit in time over fluid and structure. To solve these equations approximately, we employ predictor strategies and fixed-point iterations. We compare the method with conventional fluid–structure methods, addressing applications from transonic aeroelasticity, concerning in particular the numerical prediction of bifurcations. 2006 Elsevier B.V. All rights reserved. Keywords: Computational aeroelasticity; Fluid structure discretization; Coupling algorithms
1. Introduction The numerical simulation of large scale fluid–structure interaction problems is typically achieved by merging a given fluid code and a given structural code: due to some coupling strategy, a master program advances the discrete fluid–structure system in time by calling the time stepping modules of fluid and structure and modules that accomplish the data transfer at the interface. The discretization techniques chosen in fluid and structure, the data transfer formulas and the coupling strategy, define the underlying discrete fluid–structure system and its properties. In practice, after having chosen established fluid and structure codes, data transfer formulas and coupling strategy are often made up from adhoc considerations that hopefully result in a good method, having low memory requirements and being accurate, fast and manageable for the program1 Supported by the Deutsche Forschungsgemeinschaft (DFG) within the SFB 401 ‘‘Stro¨mungsbeeinflussung und Stro¨mungs-Struktur-Wechselwirkung an Tragflu¨geln’’ at the RWTH Aachen. E-mail address: [email protected]
0045-7825/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.01.011
mer. Here we have set out from the other perspective, namely to built a method that results in a discrete fluid– structure system with some determined property, trying to keep an eye on the practical issue, which mainly is for now, the feasibility with established fluid and structure codes. As a nonlinear model problem, we have chosen the panel flutter problem, modelled by the 2D Euler equations of gas dynamics and a structure that is a strip of a von-Karman plate. For the model problem we have designed a fluid– structure discretization that mimics the energy balance of the continuous aeroelastic system on the discrete level. The construction makes use of standard finite volume flow methods and finite element structure methods and allows non-matching grids at the fluid–structure interface. Apart from the construction of the scheme, we discuss its behaviour in numerical experiments in comparison with other coupling strategies. We document how the quality of the discretization can influence the numerical predictions, in particular when Dt varies. We also consider the constructed scheme as a reference scheme to check, to which extend cheaper schemes can reproduce its results, and can thus be believed to possess the same properties in practice.
92
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
In the literature, various approaches can be found: a comparison of popular partitioned procedures, also referred to as staggered procedures, is given in [8]. These are the loose coupling [11], the coupling of Lesoinne and Farhat [14], subcycling [17] and predictor strategies [17]. Amongst these, only the scheme of Lesoinne and Farhat has a design criterion, namely the continuity in both displacement and velocity at the fluid–structure interface without violating the discrete geometric conservation law (DGCL). The DGCL is an important compatibility condition to be met for flow computations on dynamic grids and is in fact related to the numerical stability of such computations [7]. Its satisfaction is thus another important factor in fluid– structure calculations. Other design criteria are given in [4,9]. Both references address the handling of non-matching grids at the fluid–structure interface. In [4] a conservative load transfer from fluid to structure is discussed, but still the proposed scheme is a loose coupling scheme: the fluid pressure present at the beginning of the time step is transferred conservatively from the fluid grid to the structure grid and used as a load to advance the structure in time, irrespectively of the time stepping schemes employed in fluid and structure. This does not necessarily result in a set of discrete fluid–structure equations, which possess a conservation property that holds in the discrete time evolution. Nevertheless, clear improvements over load transfers based on interpolation are observed in [4], when the load contains a discontinuity. In [9] a conservative energy transfer is proposed, but also here, a thorough treatment of the time discretization is circumvented. An idea similar to ours can be found in [16]. In [16], the exact energy conservation on the discrete level is perturbed by an OðDtÞ term, which is due to fluid discretization issues. The methodology is applied to a fluid structure problem with incompressible flow and in a pure finite element context for the space discretizations. Predictor schemes try to remove the time lags present in the loose coupling, by employing time extrapolations of the structural deflection [17] or the fluid pressure [12]. They intend to achieve high time orders and a correct energy exchange at the fluid–structure interface. This is also the case for the scheme of Bendiksen [2,5], which is an explicit five-stage Runge–Kutta time integrator, that exchanges data between fluid and structure at each Runge–Kutta stage. Other schemes, often referred to as strong coupling schemes, try to prevent losses in conservation and accuracy, by requiring a convergence criterion to be met in each time step by some iteration process on the discrete fluid–structure equations. Such schemes have been proposed in [1,15] without attributing conservation properties to the discrete equations in case the iteration is converged. In [1] an airfoil section with plunge and pitch degree-of-freedom in transonic flow is considered. The iteration method is Jameson’s dual time-stepping, where a data-transfer between fluid and structure takes place at each pseudo time step. In [15] a low Mach number flow interacting with an elastically mounted cylinder, where vortex shedding plays an important role, is simulated and considerable improve-
ments over the loose coupling are obtained for large time step sizes. The iteration is a straight forward fixed-point iteration, alternating several times between fluid and structural solver within one time step. In Section 2, we discuss the continuous model problem. In Section 3, the discretization of fluid and structure will be given and the construction of the aeroelastic scheme that satisfies our design criteria will be derived. The solution techniques for the corresponding equations will be discussed in Section 4. Finally, in Section 5, the proposed scheme will be compared in aeroelastic cases with other schemes from the literature. 2. The continuous problem Our aeroelastic model problem is the panel flutter problem [6]. The panel flutter problem consists of a plate (panel) over which a compressible fluid (air) flows. We assume the panel as infinitely long in the spanwise direction (z-coordinate), so that a 2D flow passes over a strip of the panel, compare Fig. 1. The flow is modelled by the 2D Euler equations of gas dynamics and the structure by a strip of a von-Karman plate. The panel is simply supported at its ends. It is placed on the x-axis between solid walls that stretch along the x-axis, as drawn in Fig. 1. This aeroelastic system exhibits self-excited oscillations (limit cycle oscillations) caused by the aerodynamic coupling of structural modes. This phenomena is known as flutter. Other observed long time behaviours are steady states with either deflected panel (divergence) or an undeflected panel (damped or stable behaviour). Furthermore, the aeroelastic system exhibits a transonic dip. This means that the dynamic pressure (kinetic energy of the inflow) needed to obtain unstable behaviour (flutter or divergence) drops considerably to a minimum at an inflow Mach number near unity. Another important feature of panel flutter is that shock movements play a significant role in the destabilizing mechanisms when the flow is transonic. The fluid equations are 0 1 0 Z Z Z B pn C d B 1 C T _ n ds ¼ U dx dy þ U ðv xÞ B Cds @ A dt XðtÞ pn oXðtÞ oXðtÞ 2 pvT n ð1Þ for arbitrarily moving domains X(t). The conservative state p T Þ is made up of the denvector U = (q, qu, qv, 12 qjvj2 þ c1
Fig. 1. Geometry of the 2D panel flutter problem.
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
sity q, the velocity vector v = (u, v)T and the pressure p. On a point of the boundary oX(t), n = (n1, n2)T is the outward unit normal of X(t) and x_ is the velocity of a boundary point x(t). For the formulation of the panel equation we introduce the space V :¼ H 2 ð0; lÞ \ H 10 ð0; lÞ. The variational formulation is to find the panel deflection w(t, x) satisfying Dðwxx ; uxx Þ þ N ðwÞ ðwx ; ux Þ þ mð€ w; uÞ ¼ ðp2 p1 ; uÞ
8u 2 V ;
ð2Þ
where the nonlinear term models a restoring force due to the panel stretching with Z Eh l 2 w dx. N ðwÞ ¼ 2l 0 x The constants in the panel equation are the stiffness D = E h3/12(1 m2), the mass per unit area m = qsh, the panel thickness h and length l, the density qs, Young’s modulus E and the Poisson ratio m. Now the aeroelastic problem can be formulated. For convenience, all interface and boundary conditions are given in the strong sense. • The fluid domain at time t is defined as the space above the deflecting panel and solid walls and is denoted by XF(t), i.e. for each fluid interface point x(t) associated with a panel point n 2 [0, l] T
xð0Þ ¼ ðn; 0Þ ;
T
xðtÞ ¼ ðn; wðt; nÞÞ .
ð3Þ
• For all moving domains X(t) XF(t) (1) holds. • (2) holds with the r.h.s. determined from the fluid pressure on top of the panel and p1 on the bottom p1 ðt; nÞ pðt; n; wðt; nÞÞ;
n 2 ½0; l;
p2 p1 .
ð4Þ
• At infinity we have q = q1, u = u1, v = 0, p = p1. • On the remaining boundary of XF(t) the fluid velocity in normal direction n equals the velocity of the boundary in direction n, in particular on the fluid–structure interface, _ 2. vT n ¼ x_ T n ¼ wn
ð5Þ
2.1. Energy equation of the aeroelastic problem The energy equations for fluid and structure are given by d EF þ dt Z ¼
Z
1 2 p _ T n ds qv þ ðv xÞ c1 oXF ðtÞ 2 pvT n ds;
oXF ðtÞ
d ES ¼ dt
Z
0
l
ðp2 p1 Þw_ dx;
93
where we have denoted the energy in fluid and structure by Z 1 2 p qv þ dx; EF ðtÞ 2 c 1 XF ðtÞ Z 2 Z l m 2 D 2 Eh 1 l 2 ES ðtÞ w dx . w_ þ wxx dx þ 2 2l 2 0 x 0 2 The energy equation of the fluid is the fourth component of (1) and the energy equation of the structure is derived by plugging u ¼ w_ into (2). Summation of the two energy equations yields the energy equation of the aeroelastic system. Here we restrict the fluid solution to a finite domain XF(t) with non-moving inflow and outflow boundary. Due to the boundary and coupling conditions (3)–(5), we end up with d ðEF þ ES Þ ¼ dt
Z l 1 2 pc _ qv þ p1 wdx. vT nds þ c1 0 oXin=outflow 2
Z
ð6Þ 3. Discretization 3.1. Flow discretization The flow is discretized in space by the finite volume MUSCL scheme on polygonal grids, applying a leastsquares linear reconstruction and the Venkatakrishnan limiter on the primitive fluid state vector. As a numerical flux function F we use the van-Leer flux vector splitting for moving edges. Let us write down the discrete flow equations for an inner cell, when employing the mid-point rule in time X nþ1=2 jXinþ1 jU inþ1 jXni jU ni þ Dt jeij j F
eij oXi nþ1=2 nþ1=2 nþ1=2 nþ1=2 U ij ; U ji ; x_ ij ; nij
¼ 0.
Here, the fluxes are evaluated at the half time grid position xn+1/2 = (xn+1 + xn)/2, with the corresponding normal nþ1=2 nþ1=2 nþ1=2 and length jeij j of the edge eij , which connects nij the cell Xi to its neighbour cell Xj. The edge’s velocity is given by 1 xlnþ1 xnl xnþ1 xnr þ r x_ nþ1=2 ð7Þ 2 Dt Dt and corresponds to the constant velocity by which the edge’s mid-point can be advanced from its position at time level n to its position at time level n + 1. Here, xl and xr are the two end vertices of an edge. The half time geometry xn+1/2 is used to apply reconstruction and limiting on the averaged fluid state Un+1/2 (Un+1 + Un)/2 and as a result nþ1=2 one obtains the state U ij inside the ith cell at the midpoint of the edge eij. With this choice the scheme satisfies its discrete geometric conservation law (DGCL) [7]. At inflow and outflow boundaries non-reflecting conditions according to [20] are implemented. At impermeable boundaries the flux
94
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
0
0
1
B nij;1 C B C jeij j F ðU ij ; ; x_ ij ; nij Þ ¼ jeij jpij B C @ nij;2 A x_ Tij nij
ð8Þ
is prescribed, which results from the slip-condition vT n ¼ x_ T n. Another time stepping scheme, used in the numerical experiments later, is the important second order backward differencing scheme. For this scheme the DGCL is satisfied with the choice of edge normals and velocities as given in [13]. 3.2. Structural discretization Conformal elements in one dimension for fourth order operators are constructed from the Hermite shape functions on the unit interval, n 2 [0, 1], s0;l ðnÞ ¼ ð1 þ 2nÞð1 nÞ2 ; s0;r ðnÞ ¼ ð3 2nÞn2 ; s1;l ðnÞ ¼ nð1 nÞ2 ; s1;r ðnÞ ¼ ð1 nÞn2 . We have n 1 inner nodes, so that 2(n 1) + 2 basis functions /l(x) are constructed from these shape functions when the panel is simply supported. We plug X wðt; xÞ ¼ al ðtÞ/l ðxÞ l
into (2) and test with the basis functions. We introduce bl a_ l to obtain a first order system of ODEs X Dal ð/l ; /k Þ2 þ al N ðwx Þð/l ; /k Þ1 þ mb_ l ð/l ; /k Þ0 l
¼ ðp2 p1 ; /k Þ0 a_ k ¼ bk
8k;
8k.
Here we have set ( )n+1/2 (( )n+1 + ( )n)/2. In particular, the nonlinearity N is evaluated with a trapezoidal rule. The evaluation of the load, ðp1 p2 Þjtnþ1=2 , is not made explicit now, but will be specified P below for the aeroelastic application. DefiningPwn ðxÞ l anl /l ðxÞ and introducing the velocity w_ n ðxÞ l bnl /l ðxÞ, the system reads 9 N nþ1=2 ¼ ðN ðwnþ1 Þ þ N ðwn ÞÞ=2; > > > > = DtDðwnþ1=2 ; /k Þ2 þ DtN nþ1=2 ðwnþ1=2 ; /k Þ1 ð9Þ nþ1 n þ mðw_ w_ ; /k Þ0 ¼ Dtððp2 p1 Þjtnþ1=2 ; /k Þ0 ; > > > > ; nþ1=2 anþ1 ank ¼ Dtbk : k 3.3. Mimicing the energy budget of the continuous problem on the discrete level By definition of the finite volume method as a conservative scheme, the energy conservation properties are carried over from the continuous to the discrete level in the fluid. In the structure, the choice of the time integrator decides if the energy conservation is reflected correctly in the discrete scheme. Here, after the spacial FEM discretization, the time integration scheme (9) is a scheme with the desired property. Furthermore, the use of non-matching interface representations has to be taken care of, so that the energy exchange at the interface is conservative. Here, due to the use of fluid cells with polygonal boundary, a piecewise linear representation of the interface is given for the fluid. Whereas the panel deflection is a piecewise polynomial of third order due to the finite element space employed. And in general, the fluid vertices and the structural nodes do not coincide on the interface, compare Fig. 2. These points will be addressed in the following construction of the aeroelastic scheme. This scheme is constructed such that the energy conservation properties of the continuous problem, compare (6), are mimiced on the discrete level. Using again the mid-point rule in time, the energy balance for a cell Xi can be written as
Here, we have used abbreviations for the bilinear forms ð/; wÞ0 ð/; wÞ; ð/; wÞ1 ð/x ; wx Þ; ð/; wÞ2 ð/xx ; wxx Þ. The set of 4(n 1) + 4 ODEs may now be integrated in time with any kind of Runge–Kutta scheme to approximately solve for discrete values of the coefficients al(t), bl(t). The following second order implicit rule will be used below in the construction of an aeroelastic scheme with correct energy transfer: X nþ1=2 nþ1=2 nþ1=2 Dal ð/l ; /k Þ2 þ al N ð/l ; /k Þ1 l
blnþ1 bnl ð/l ; /k Þ0 ¼ ððp2 p1 Þjtnþ1=2 ; /k Þ0 ; Dt aknþ1 ank nþ1=2 ¼ bk . Dt þm
Fig. 2. Representation of the interface in fluid and structure, the Nj are panel nodes.
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
X
nþ1 EF;i EnF;i ¼ Dt
nþ1=2
jeij
j
eij oXi nþ1=2
F 4 ðU ij
nþ1=2
; U ji
nþ1=2
; x_ ij
nþ1=2
; nij
Þ.
Here F4 is the fourth component of the numerical flux F and all its arguments are taken at the edge center. The energy contained in the ith cell at time level n is EnF;i ¼ n n 2 n jXni j ðq i jvi j =2 þ p i =ðc 1ÞÞ. The discrete fluid energy is P n n EF i EF;i . We can divide the energy production into nþ1=2 the contribution oEFS coming from the interface and nþ1=2
the contribution oEF Thus we have nþ1=2
EFnþ1 EnF ¼ oEFS
from the remaining boundary. nþ1=2
þ oEF
.
The energy flux at impermeable non-moving boundary parts vanishes due to (8). Thus the only contribution nþ1=2 remaining in oEF comes from inflow and outflow boundaries, so that we can write nþ1=2
EFnþ1 EnF ¼ oEFS
nþ1=2
þ oEin=outflow .
ð10Þ
The energy flux at an interface edge, due to (8), is given by nþ1=2
F 4 ðU ij
nþ1=2
; ; x_ ij
nþ1=2
; nij
nþ1=2
Þ pij
nþ1=2 T nþ1=2 Þ nij
ðx_ ij
Dt ðp2 p1 Þjtnþ1=2 ; /k 0 ¼ Dt p1 pðn; wnþ1=2 ðnÞÞ; /k ðnÞ 0 X nþ1=2
Dt Dnij ðp1 pij Þ /k ðnl;ij Þ þ /k ðnr;ij Þ =2. eij CFS
ð12Þ Note that most summands do not contribute anything since in most cases nl,ij, nr,ij 62 supp(/k). For the evaluation of the bilinear products on the l.h.s. of (9) we may use a different quadrature rule, let us assume that this rule is exact. Plugging the quadrature formula (12) into (9), multiplying nþ1=2 the resulting equation with bk ¼ ðaknþ1 ank Þ=Dt and summing over all k, we obtain DtDðwnþ1=2 ; w_ nþ1=2 Þ2 þ DtN nþ1=2 ðwnþ1=2 ; w_ nþ1=2 Þ1 þ mðw_ nþ1 w_ n ; w_ nþ1=2 Þ0
X Dnij ðp1 pijnþ1=2 Þ wnþ1 ðnl;ij Þ wn ðnl;ij Þ ¼ Dt 2 Dt eij CFS nþ1 n w ðnr;ij Þ w ðnr;ij Þ þ . Dt From the third equation of (9) we have wnþ1 wn ¼ Dtw_ nþ1=2 .
and summing over all edges on the interface, we obtain X nþ1=2 nþ1=2 nþ1=2 nþ1=2 nþ1=2 jeij j pij ðx_ ij ÞT nij . oEFS Dt eij CFS
Let us look at a single edge e xl xr CFS with end vertices xl, xr, omitting the ij-indices for the moment, compare Fig. 2. Then x_ nþ1=2 is given by (7). On the interface, we lock the fluid vertices on material points of the panel, demanding that T
xnl ¼ ðnl ; wnl Þ where wnl wn ðnl Þ for some fixed nl 2 ½0; l and correspondingly for the right point xnr and some fixed nr 2 [0, l]. This is the discretization of (3). Thus we get 1 wlnþ1 wnl wrnþ1 wnr nþ1=2 nþ1=2 T nþ1=2 þ Þ n ¼ ðx_ n2 2 Dt Dt on the interface. Looking again at specific edges, adding nþ1=2 nþ1=2 and the indices ij, using jeij j ¼ ðnr;ij nl;ij Þ=n2;ij denoting Dnij nr,ij nl,ij, we obtain X nþ1=2 ¼ Dt Dnij pijnþ1=2 oEFS eij CFS
95
1 wnþ1 ðnl;ij Þ wn ðnl;ij Þ wnþ1 ðnr;ij Þ wn ðnr;ij Þ þ . 2 Dt Dt ð11Þ
We turn to the discrete panel equations (9). We have to specify the load term due to the coupling condition (4). We choose the following quadrature rule to evaluate the load term in (9):
Differentiating s times with respect to x, multiplying by nþ1=2 ak /k , respectively its derivatives, summing over all k and integrating, gives 1 nþ1 ðw wn ; wnþ1 þ wn Þs ¼ Dtðw_ nþ1=2 ; wnþ1=2 Þs 2 for s ¼ 0; 1; 2 or 1 fðwnþ1 ; wnþ1 Þs ðwn ; wn Þs g ¼ Dtðw_ nþ1=2 ; wnþ1=2 Þs 2 for s ¼ 0; 1; 2. fðwnþ1 ; wnþ1 Þ1 þ ðwn ; wn Þ1 g=2 and Note that N nþ1=2 ¼ Eh 2l thus DtN nþ1=2 ðwnþ1=2 ; w_ nþ1=2 Þ1 ¼
Eh ðwnþ1 ; wnþ1 Þ1 þ ðwn ; wn Þ1 ðwnþ1 ; wnþ1 Þ1 ðwn ; wn Þ1 2l 2 2 2
2
Eh ðwnþ1 ; wnþ1 Þ1 ðwn ; wn Þ1 . 2l 4 With the discrete panel energy 2 m n n Eh ðwn ; wn Þ1 D n ES ðw_ ; w_ Þ0 þ þ ðwn ; wn Þ2 2 2l 2 2 ¼
the energy change on the discrete level is X Dnij ðp1 pijnþ1=2 Þ 2 eij CFS
nþ1 w ðnl;ij Þ wn ðnl;ij Þ þ wnþ1 ðnr;ij Þ wn ðnr;ij Þ .
Enþ1 EnS ¼ S
ð13Þ
96
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
Now, from (10), (11), and (13) we see that in the energy change of the discrete fluid–structure system, the contributions from the interface cancel, and EFnþ1
þ
Enþ1 S
EnF
EnS
X Dnij p1
nþ1=2 wnþ1 ðnl;ij Þ wn ðnl;ij Þ ¼ oEin=outflow þ 2 eij CFS nþ1 þ w ðnr;ij Þ wn ðnr;ij Þ .
ð14Þ
results. The following theorem holds: Theorem 1. The numerical scheme P nþ1=2 jXnþ1 jU inþ1 jXni jU ni þ Dt jeij j i
9 > > > > > > > > =
eij oXi
nþ1=2 nþ1=2 nþ1=2 nþ1=2 F ðU ij ; U ji ; x_ ij ; nij Þ
¼ 0;
with the grid velocity (7), the boundary treatment (8), and > > > T T > > > xnl;ij ¼ nl;ij ; wn ðnl;ij Þ ; xnr;ij ¼ nr;ij ; wn ðnr;ij Þ > > ; for all flow edges eij on the panel, 9 > DtDðw ; /k Þ2 þ DtN ðw ; /k Þ1 > > > P nþ1=2 > nþ1 n þ mðw_ w_ ; /k Þ0 ¼ Dt Dnij ðp1 pij Þ > = nþ1=2
nþ1=2
eij CFS
/k ðnl;ij Þ þ /k ðnr;ij Þ =2; nþ1=2
aknþ1 ank ¼ Dtbk
;
ð15Þ
nþ1=2
> > > > > > ; ð16Þ
approximating the aeroelastic initial-boundary value problem of Section 2, carries over the energy balance from the continuous level (6) to the discrete level (14). The discretization satisfies the DGCL in the fluid. Remark 2. The discretization given in the theorem employs, even on the interface, the mid-point rule for the time integration of the fluid equations and a mixed trapezoidal-midpoint rule for the time integration of the structural equations; the nonlinearity N is evaluated with a trapezoidal rule, whereas for the fluid load the mid-point rule is used. We may expect that this scheme is a second order time approximation of our aeroelastic problem. This possible increase in order may be another beneficial point of the scheme apart from its conservation property. It may be the dominating reason for the improvement over conventional fluid–structure methods shown in the numerical experiments below.
4. Algorithms An algorithm that solves the discrete flow equations (15) with a prescribed boundary movement, may serve as a fluid time integrator. Similarly, a structural time integrator is given if the Eq. (16) can be solved when a load is prescribed. In both cases we have employed a matrix-free Newton-GMRes method to solve the discrete equations. With the resulting fluid time integrator and structural time integrator, the aeroelastic problem can be tackled with all partitioned schemes. Amongst the partitioned schemes we have implemented, are the loose coupling, the coupling of Lesoinne and Farhat (denoted from now as the LFcoupling), a predictor scheme and a fixed-point iteration that solves the coupled set (15) and (16). The latter method applies the fluid integrator and the structural integrator several times within one time step. They are successively applied in the sequence given by the loose coupling with the purpose of diminishing the residual of the coupled system (15) and (16). In the sequel we denote this coupling strategy based on the fixed-point iteration as the FPI-coupling. The loose coupling and the LF-coupling are schematically depicted in Fig. 3. In the figure the time axis goes from left to right and the first three time levels are shown. The encircled numbers give the order in which structural time integration (bottom line), fluid grid movement (vertical lines) and fluid time integration (top line) are performed. And we have indicated from which time level the fluid pressure is taken as the load to perform a structural time integration. Furthermore, we can see how the grid position is transferred at the interface from structure (w) to fluid (x). Using again the mid-point rule as the time integrator, we have indicated from which time level the fluid flux is taken. In the loose coupling and the LF-coupling the quadrature rule defining the load transfer from fluid to structure is taken to be exact. This means that the product of the piecewise linear fluid pressure and a structural test function is integrated by placing adequately Gauss quadrature points. We make a remark at this point: Remark 3. We have seen in numerical experiments that when altering the FPI-coupling by replacing the load quadrature in (16) by exact quadrature, the altered FPIcoupling brings results identical to the ones obtained with the original FPI-coupling. It seems that the spacial resolution, necessary to resolve the flow details sufficiently
Fig. 3. Loose coupling (left) and LF-coupling (right). LF-coupling: in step 2, for time level 1, the fluid grid position at the interface is determined by x1 ¼ x0 þ Dtð0; w_ 1=2 ÞT and accordingly, in step 5, for time level 2, by x2 ¼ x1 þ Dtð0; w_ 3=2 ÞT .
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
(e.g. the shock movements, compare Section 5), is too fine to let any differences occur. We stress here that thus the differences we will observe in Section 5 stem from the distinct time discretization used in (15) and (16). Our predictor scheme uses the idea of [12], where the fluid pressure on the interface is predicted in time by a linear extrapolation. This means that our predictor scheme aims to approximate the solution of (15) and (16), i.e. it is precisely defined, when we replace pn+1/2 by 1.5pn 0.5pn1 in (16). In the numerical experiments we denote the scheme as PRED. Another predictor scheme, which is not related to (15) and (16) and is very close to the one proposed in [12], can be obtained by using the second order BDF time integrator in the fluid in conjunction with (16) for the structure, but again replacing pn+1/2 by 1.5pn 0.5pn1. This is a predictor scheme that is not related to a scheme with a design criterion such as PRED is. In the numerical experiments we denote the scheme as PREDBDF. The matrix-free Newton-GMRes solves a nonlinear system G(x) = 0 by applying a Newton iteration and employing GMRes [19] as a linear solver. Since the GMRes iteration only needs to know the evaluation of matrix vec0 0 tor products G (x*)x, with the Jacobian G , the determina0 tion and storage of G may be avoided by approximating G0 ðx Þx
Gðx þ xÞ Gðx Þ .
Here should be chosen big enough to avoid cancellation and small enough to keep the approximation accurate. In all matrix-free Newton-GMRes iterations we have used ¼
eps1=2 jjxjj2
where eps is the machine-accuracy,
which was also recommended in [18]. To accelerate convergence in the linear solvers, a block symmetric Gauss–Seidel preconditioner based on a first order spacial discretization is used in the fluid, whereas a block symmetric Jacobi preconditioner is used in the structure.
97
in the initial guess, then a stagnation of the fluid residual in the following iteration steps can be accepted. 4.2. Grid movement In general applications, the grid movement algorithm, which assigns the vertex positions according to the given boundary movement, may require the solution of a system of equations as-well. In our case, the computational domain has a simple geometry, also when it deforms in the aeroelastic simulation. A simple algebraic rule to stretch the grid suffices and is implemented here. 5. Numerical experiments We mainly examine the behaviour of different coupling schemes with respect to Dt. A sufficiently accurate spacial resolution of 80 fluid cells on the panel and 26 panel nodes was chosen for all computations presented, compare [3,5, 10,12]. In total the fluid grid consists of 28,000 quadrilateral cells, reaching the in- and outflow boundaries at 20 panel lengths away from the panel. Away from the panel, the cells are stretched towards the in- and outflow boundaries. Throughout, we use aluminium as panel material, thus fixing the panel parameters qs = 2700 kg/m3, E = 7.1 · 1010 Pa and m = 0.34. In the dimensionless form all coefficients of the panel equation are determined, when the ratio h/l is specified, which will be done in each of the following cases. The fluid parameters may be defined in two ways: either we assign the inflow Mach number M1, the fluid structure mass ratio l q1l/(qsh) and the dimensionless dynamic inflow pressure k q1 u21 l3 =D or we assign M1 and all remaining flow conditions as those found in the atmosphere at a certain altitude. The panel is taken as undeflected at the initial time, but with a velocity profile, which would cause the panel to oscillate in its first fundamental mode with an amplitude of 0.2 h, if the nonlinearity and the load were not present. Only in Section 5.5 a different initial velocity distribution is used. 5.1. Transonic and supersonic panel flutter
4.1. Convergence criterion for the solution of (15) and (16) As an initial guess in the fixed-point iteration to solve (15) and (16), the result of a loose coupling step can be used. The loose coupling first solves the structural equations with the wrong fluid load and then solves the fluid equations due to the obtained structural deflections. Thus the structural equations of the loose coupling and of (16) do not coincide, whereas the fluid equations of the loose coupling and of (15) do coincide. This means that, when plugging the initial guess into the system (15) and (16), the error is only visible in the residual of the structural equations. Thus, the convergence criterion is to reduce the structural residual below a given threshold. If the residual threshold in the fluid equation is already small enough
Let us first envision some qualitative properties of the flutter solutions. We look at a transonic flutter case with M1 = 0.95, h/l = 0.002, l = 0.1, k = 2392 and a supersonic flutter case defined through h/l = 0.004 and an inflow of M1 = 1.2 at 20,000 feet altitude. The latter case and cases similar to the first one were also considered in [3,5]. The periodic panel deflections resulting in these cases are visualized in Fig. 4 by the deflection at 25%, at 50% and at 75% chord of the panel. Higher harmonics in the panel deflection are more pronounced in the transonic case in contrast to the supersonic case. In Fig. 5, the development of the Mach number distribution on the panel over a flutter cycle is shown. We can observe that the Mach number on the panel has a moving shock in the limit cycle in the transonic case,
98
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102 25 % 50 % 75 %
M = 0.95
10 8
6
25 % 50 % 75 %
M = 1.2
4
6 4
2 w/h
w/h
2 0
0
-2 -2
-4 -6
-4
-8 -10
0
50
100 t / t*
-6 260
270
280
290
300
t / t*
Fig. 4. The panel deflections at various points on the panel.
Fig. 5. Mach numbers along the panel during one period of the limit cycle. M1 = 0.95 on the left, M1 = 1.2 on the right.
whereas in the supersonic case the fluid quantities are smooth along the panel. In the supersonic case shocks periodically build up and disappear at the trailing and the leading edge. In the transonic case we can relate the shock speed to the CFL number in the fluid, and for the given grid, at a CFL number of roughly 100, the shock moves about one cell during a time step, furthermore at this CFL number, we have roughly 60 time steps per flutter period. 5.2. Flutter at M1 = 0.95 Except for k, which we take now as k = 3000, the transonic flutter case considered here is defined with the same data as the one in the previous section. Here we compare the behaviour of the solutions obtained with the LF-coupling and the FPI-coupling, as Dt ! 0. Increasing the fluid CFL number from 10, 30, 50 to 70, we can observe in Fig. 6 how time accuracy is lost in the LF-coupling. We even obtain a ‘‘numerical bifurcation’’ from flutter to divergence with the LF-coupling, when changing from CFL = 50 to CFL = 70. At the latter CFL number the run with the LF-coupling fails at t/t* 50. In contrast, the results obtained with the FPI-coupling remain very accurate even at a CFL number of 70. In comparison with Fig. 7 we see that the predictor method also works quite accurate, but
less accurate then the FPI-coupling. The loose coupling fails to detect a flutter solution at CFL = 50 already. 5.3. Flutter at M1 = 1.2 In this case, for the given fluid grid, a fluid CFL number of 20 corresponds to 56 time steps per flutter period. In Figs. 8 and 9, we see that the loose coupling and the LFcoupling have the tendency to build up the deflection amplitude much earlier than the predictor or the FPIcoupling. A similar tendency can be observed in the loose coupling, the LF-coupling and the predictor coupling when we increase Dt. 5.4. A bifurcation at M1 = 0.95 Another numerical experiment in the transonic regime addresses the determination of a bifurcation point with different coupling schemes and time step sizes. Increasing the dynamic pressure at fixed Mach number a bifurcation from divergence to flutter occurs, where strong shocks move along the panel in the flutter solution. Again an aluminium panel with h/l = 0.002 is considered and we fix M1 = 0.95 and l = 0.1. The non-dimensional dynamic pressure is varied in the range 2000 6 k 6 4000. The bifurcation points
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102 CFL = 10 CFL = 30 CFL = 50 CFL = 70
8
8
6
6
4
4
2
2
0
0
-2
-2
-4
-4
-6
-6
-8
-8
-10
0
50
100
150
CFL = 10 CFL = 30 CFL = 50 CFL = 70
10
w/h
w/h
10
99
-10
0
50
t / t*
100
150
t / t*
Fig. 6. Panel mid-point deflection for LF-coupling and FPI-coupling at various Dt; M1 = 0.95.
CFL = 10 CFL = 30 CFL = 50 CFL = 70
8
8
6
6
4
4
2
2
0
0
-2
-2
-4
-4
-6
-6
-8
-8
-10
-10
0
50
100
150
CFL = 10 CFL = 30 CFL = 50 CFL = 70
10
w/h
w/h
10
0
50
t / t*
100
150
t / t*
Fig. 7. Panel mid-point deflection for loose coupling and predictor coupling at various Dt; M1 = 0.95.
CFL = 10 CFL = 20 CFL = 40
4
4
2
2
0
0
-2
-2
-4
-4
-6 50
100
150
200
t / t*
CFL = 10 CFL = 20 CFL = 40
6
w/h
w/h
6
-6 50
100
150
200
t / t*
Fig. 8. Panel mid-point deflection for loose coupling and predictor coupling at various Dt; M1 = 1.2.
obtained are shown in Fig. 10 and reveal a strong variation depending on the coupling scheme, when using moderate CFL numbers. As we know from Section 5.1, a CFL num-
ber of 50 corresponds to 120 time steps per flutter period for flutter solutions in the considered range of dynamic pressures. Clearly, the loose coupling and the LF-coupling
100
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102 CFL = 10 CFL = 20 CFL = 40
4
4
2
2
0
0
-2
-2
-4
-4
-6 50
100
150
200
CFL = 10 CFL = 20 CFL = 40
6
w/h
w/h
6
-6 50
100
t / t*
150
200
t / t*
Fig. 9. Panel mid-point deflection for LF-coupling and FPI-coupling at various Dt; M1 = 1.2.
4000 LOOSEBDF LOOSE LF PREDBDF PRED FPI
3800 3600 3400
λ
3200 3000 2800 2600 2400 2200 2000 0
10
20
30
40
50
detected around k = 2500. We obtain the diagram of Fig. 11. In contrast to Fig. 10 all schemes now underestimate the bifurcation point more or less at higher CFL numbers. At CFL = 10 there is no entry for the FPI scheme. In this singular case we observe that the fixedpoint iteration does not converge in some of the first 50 time steps. The result is, that the FPI produces only flutter solutions in the k-range around the bifurcation point, although it performs very well at higher CFL numbers. A theoretical argument for such failure in convergence at small Dt is given in [16] in a similar context. As a cure, a relaxation strategy is proposed in [16], which we should also test in the future.
CFL
overestimate the dynamic pressure at the bifurcation point immensely at CFL = 50. The relative error in the dynamic pressure at the bifurcation point is 10% for the loose coupling and 9% for the LF-coupling, both at CFL = 10, whereas the relative errors are 3% for the FPI-coupling and 15% for the predictor scheme, both at CFL = 50. Further, we see that at CFL = 50 the schemes LOOSEBDF and PREDBDF are less accurate than their counterparts LOOSE and PRED. The latter ones can be interpreted as approximations to the conservative method (15) and (16) by employing an estimate for the pressure pn+1/2, thereas the former ones deviate even further from 15 and 16 by using additionally BDF as the time integrator in the fluid instead of the mid-point rule. 5.5. Change of the panel’s initial velocity distribution We now change the initial amplitude of the panel velocity to 0.00005 m/s and keep the other parameters of the above bifurcation problem fixed. In [10] the bifurcation is
5.6. Convergence of FPI We have picked the case with M1 = 0.95 from Section 5.2 at CFL = 50 to exhibit the convergence behaviour of the fixed-point iteration for the solution of (15) and (16). In Fig. 12, the behaviour of the structural residual over a
2600
LOOSE LF PRED FPI
2500
2400
λ
Fig. 10. Bifurcation at M1 = 0.95. For each coupling scheme we find divergence below the corresponding line and flutter above.
2300
2200 0
10
20
CFL
30
40
50
Fig. 11. Bifurcation at M1 = 0.95 with initial panel velocity 0.00005 m/s. For each coupling scheme we find divergence below the corresponding line and flutter above.
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
the fixed-point iteration. From the fluid residual we demand that it is diminished due to the assigned threshold given in the fluid solver. We do not need a decrease of the fluid residual after each fixed-point iteration step. Here, a sufficient threshold in the fluid solver is to reduce the residual below 108. We observe in Fig. 13 that the FPI always keeps the fluid residual below 108.
Log (structural residual)
1
0
-1
-2
6. Conclusions
-3
-4 0
50
100
150 200 Time step
250
300
Fig. 12. Behaviour of structural residual.
range of 300 time steps can be seen. This covers the dimensionless time interval from 0 to 84 in the right plot of Fig. 6, i.e. during these time steps, the solution has already gone through one flutter cycle. After each fluid–structure loop the residual is plotted in the diagram. At each time step we have at least two entries, the one from the initial guess and the one from the following iteration loop. The structural residual of the initial guess, which stems from a loose coupling step, lies around 100.2. This can be red off from the accumulation of dots producing a wavy shape on top of the diagram, whereas after a single FPI iteration the residual is around 101.7, as red off from the wavy shape on the bottom. Throughout, the FPI reduces with one iteration step the structural residual by more than one order of magnitude, which was the convergence criterion for all computations shown for the FPI-coupling. Thus, together with the determination of the initial guess, a numerical simulation with the FPI-coupling needs roughly twice the computational work as the numerical simulation with the loose coupling at the same CFL number. Analogous plots for the fluid residuals can be found in Fig. 13. As explained in Section 4.1, the structural residual is the crucial one, when we judge about the convergence of
The panel flutter problem was considered as a model problem to investigate upon numerical methods for nonlinear aeroelastic problems. We have proposed a discretization of the aeroelastic initial-boundary value problem that is designed such that the energy conservation property of the continuous model is carried over to the discrete model. It employs trapezoidal and mid-point rule as time integrators and deals with non-matching grids at the interface. The DGCL is satisfied in the fluid. The resulting scheme is implicit over the fluid and structure unknowns. We have found that a fixed point iteration (FPI-coupling) solves the system sufficiently, resulting in a fluid–structure algorithm with roughly twice the computational cost of a loose coupling scheme. A singular case, where we encounter convergence problems at small Dt may be dealt with by introducing a relaxation strategy in the future, an idea that is theoretically supported by [16] in a similar context. A predictor scheme derived from the proposed discretization, the scheme proposed in [14] and the loose coupling all have roughly the same computational cost and are compared to the FPI-coupling in transonic aeroelastic cases. We find that the FPI-coupling is still very accurate at moderate to large Dt, i.e. the range of 120 time steps per flutter cycle. We can observe a loss of accuracy at these time step sizes for the predictor scheme, and we observe a severe loss of accuracy at even smaller Dt for the loose coupling and the scheme of [14]. These observations were made by numerical experiments when predicting a bifurcation in the transonic regime, where shock movements play a significant role. References
-6
Log (fluid residual)
101
-8
-10
0
50
100
150 200 Time step
250
Fig. 13. Behaviour of fluid residual.
300
[1] J.J. Alonso, A. Jameson, Fully-implicit time-marching aeroelastic solutions, AIAA-paper 94-0056, 1994. [2] O. Bendiksen, A new approach to computational aeroelasticity, AIAA-paper 91-0939-CP, 1991. [3] O. Bendiksen, G. Davis, Nonlinear travelling wave flutter of panels in transonic flow, AIAA-paper 95-1496, 1995. [4] J.R. Cebral, R. Lo¨hner, Conservative load projection and tracking for fluid–structure problems, AIAA J. 35 (4) (1997) 678–692. [5] G. Davis, O. Bendiksen, Transonic panel flutter, AIAA-paper 931476, 1993. [6] E.H. Dowell, Aeroelasticity of Plates and Shells, Noordhoff International Publishing, 1975. [7] C. Farhat, P. Geuzaine, C. Grandmont, The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids, J. Comput. Phys. 174 (2001) 669–694.
102
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
[8] C. Farhat, M. Lesoinne, On the accuracy, stability, and performance of the three-dimensional nonlinear transient aeroelastic problems by partitioned procedures, AIAA-paper 96-1388-CP, 1996. [9] C. Farhat, M. Lesoinne, P. LeTallec, Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: momentum and energy conservation, optimal discretization and application to aeroelasticity, Comput. Methods Appl. Mech. Engrg. 157 (1998) 95–114. [10] R.E. Gordnier, M.R. Visbal, Development of a three-dimensional viscous aeroelastic solver for nonlinear panel flutter, AIAA-paper 2000-2337, 2000. [11] G. Guruswamy, Coupled finite-difference/finite-element approach for wing-body aeroelasticity, AIAA-paper 92-4680-CP, 1992. [12] J. Hurka, J. Ballmann, Elastic panels in transonic flow, AIAA-paper 2001-2722, 2001. [13] B. Koobus, C. Farhat, Second-order time-accurate and geometrically conservative implicit schemes for flow computations on unstructured dynamic meshes, Comput. Methods Appl. Mech. Engrg. 170 (1999) 103–129.
[14] M. Lesoinne, C. Farhat, Geometric conservation laws for flow problems with moving boundaries, and their impact on aeroelastic computations, Comput. Methods Appl. Mech. Engrg. 134 (1996) 71–80. [15] S. Morton, R. Melville, M. Visbal, Accuracy and coupling issues of aeroelastic Navier–Stokes solutions on deforming meshes, J. Aircraft 35 (5) (1998) 798–805. [16] P. LeTallec, J. Mouro, Fluid structure interaction with large structural displacements, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3039–3067. [17] S. Piperno, Explicit/implicit fluid/structure staggered procedures with a structural predictor and fluid subcycling for 2D inviscid aeroelastic simulations, Int. J. Numer. Methods Fluids 25 (1997) 1207–1226. [18] N. Qin, D. Ludlow, S. Shaw, A matrix-free preconditioned Newton/ GMRES method for unsteady Navier–Stokes solutions, Int. J. Numer. Methods Fluids 33 (2000) 223–248. [19] Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Publishing, New York, 1996. [20] K. Thompson, Time dependent boundary conditions for hyperbolic systems, J. Comput. Phys. 68 (1987) 1–24.
Discrete conservation and coupling strategies in nonlinear aeroelasticity Ralf Massjung
1
Institut fu¨r Geometrie und Praktische Mathematik, RWTH Aachen, Templergraben 55, 52056 Aachen, Germany Received 30 December 2004; accepted 30 January 2006
Abstract As a model problem of aeroelasticity the panel flutter problem is considered to study the discretization of nonlinear aeroelastic problems and the behaviour of the corresponding solution strategies. A fluid–structure discretization that mimics the energy budget of the continuous problem is derived. The set of equations occurring in each time step are implicit in time over fluid and structure. To solve these equations approximately, we employ predictor strategies and fixed-point iterations. We compare the method with conventional fluid–structure methods, addressing applications from transonic aeroelasticity, concerning in particular the numerical prediction of bifurcations. 2006 Elsevier B.V. All rights reserved. Keywords: Computational aeroelasticity; Fluid structure discretization; Coupling algorithms
1. Introduction The numerical simulation of large scale fluid–structure interaction problems is typically achieved by merging a given fluid code and a given structural code: due to some coupling strategy, a master program advances the discrete fluid–structure system in time by calling the time stepping modules of fluid and structure and modules that accomplish the data transfer at the interface. The discretization techniques chosen in fluid and structure, the data transfer formulas and the coupling strategy, define the underlying discrete fluid–structure system and its properties. In practice, after having chosen established fluid and structure codes, data transfer formulas and coupling strategy are often made up from adhoc considerations that hopefully result in a good method, having low memory requirements and being accurate, fast and manageable for the program1 Supported by the Deutsche Forschungsgemeinschaft (DFG) within the SFB 401 ‘‘Stro¨mungsbeeinflussung und Stro¨mungs-Struktur-Wechselwirkung an Tragflu¨geln’’ at the RWTH Aachen. E-mail address: [email protected]
0045-7825/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.01.011
mer. Here we have set out from the other perspective, namely to built a method that results in a discrete fluid– structure system with some determined property, trying to keep an eye on the practical issue, which mainly is for now, the feasibility with established fluid and structure codes. As a nonlinear model problem, we have chosen the panel flutter problem, modelled by the 2D Euler equations of gas dynamics and a structure that is a strip of a von-Karman plate. For the model problem we have designed a fluid– structure discretization that mimics the energy balance of the continuous aeroelastic system on the discrete level. The construction makes use of standard finite volume flow methods and finite element structure methods and allows non-matching grids at the fluid–structure interface. Apart from the construction of the scheme, we discuss its behaviour in numerical experiments in comparison with other coupling strategies. We document how the quality of the discretization can influence the numerical predictions, in particular when Dt varies. We also consider the constructed scheme as a reference scheme to check, to which extend cheaper schemes can reproduce its results, and can thus be believed to possess the same properties in practice.
92
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
In the literature, various approaches can be found: a comparison of popular partitioned procedures, also referred to as staggered procedures, is given in [8]. These are the loose coupling [11], the coupling of Lesoinne and Farhat [14], subcycling [17] and predictor strategies [17]. Amongst these, only the scheme of Lesoinne and Farhat has a design criterion, namely the continuity in both displacement and velocity at the fluid–structure interface without violating the discrete geometric conservation law (DGCL). The DGCL is an important compatibility condition to be met for flow computations on dynamic grids and is in fact related to the numerical stability of such computations [7]. Its satisfaction is thus another important factor in fluid– structure calculations. Other design criteria are given in [4,9]. Both references address the handling of non-matching grids at the fluid–structure interface. In [4] a conservative load transfer from fluid to structure is discussed, but still the proposed scheme is a loose coupling scheme: the fluid pressure present at the beginning of the time step is transferred conservatively from the fluid grid to the structure grid and used as a load to advance the structure in time, irrespectively of the time stepping schemes employed in fluid and structure. This does not necessarily result in a set of discrete fluid–structure equations, which possess a conservation property that holds in the discrete time evolution. Nevertheless, clear improvements over load transfers based on interpolation are observed in [4], when the load contains a discontinuity. In [9] a conservative energy transfer is proposed, but also here, a thorough treatment of the time discretization is circumvented. An idea similar to ours can be found in [16]. In [16], the exact energy conservation on the discrete level is perturbed by an OðDtÞ term, which is due to fluid discretization issues. The methodology is applied to a fluid structure problem with incompressible flow and in a pure finite element context for the space discretizations. Predictor schemes try to remove the time lags present in the loose coupling, by employing time extrapolations of the structural deflection [17] or the fluid pressure [12]. They intend to achieve high time orders and a correct energy exchange at the fluid–structure interface. This is also the case for the scheme of Bendiksen [2,5], which is an explicit five-stage Runge–Kutta time integrator, that exchanges data between fluid and structure at each Runge–Kutta stage. Other schemes, often referred to as strong coupling schemes, try to prevent losses in conservation and accuracy, by requiring a convergence criterion to be met in each time step by some iteration process on the discrete fluid–structure equations. Such schemes have been proposed in [1,15] without attributing conservation properties to the discrete equations in case the iteration is converged. In [1] an airfoil section with plunge and pitch degree-of-freedom in transonic flow is considered. The iteration method is Jameson’s dual time-stepping, where a data-transfer between fluid and structure takes place at each pseudo time step. In [15] a low Mach number flow interacting with an elastically mounted cylinder, where vortex shedding plays an important role, is simulated and considerable improve-
ments over the loose coupling are obtained for large time step sizes. The iteration is a straight forward fixed-point iteration, alternating several times between fluid and structural solver within one time step. In Section 2, we discuss the continuous model problem. In Section 3, the discretization of fluid and structure will be given and the construction of the aeroelastic scheme that satisfies our design criteria will be derived. The solution techniques for the corresponding equations will be discussed in Section 4. Finally, in Section 5, the proposed scheme will be compared in aeroelastic cases with other schemes from the literature. 2. The continuous problem Our aeroelastic model problem is the panel flutter problem [6]. The panel flutter problem consists of a plate (panel) over which a compressible fluid (air) flows. We assume the panel as infinitely long in the spanwise direction (z-coordinate), so that a 2D flow passes over a strip of the panel, compare Fig. 1. The flow is modelled by the 2D Euler equations of gas dynamics and the structure by a strip of a von-Karman plate. The panel is simply supported at its ends. It is placed on the x-axis between solid walls that stretch along the x-axis, as drawn in Fig. 1. This aeroelastic system exhibits self-excited oscillations (limit cycle oscillations) caused by the aerodynamic coupling of structural modes. This phenomena is known as flutter. Other observed long time behaviours are steady states with either deflected panel (divergence) or an undeflected panel (damped or stable behaviour). Furthermore, the aeroelastic system exhibits a transonic dip. This means that the dynamic pressure (kinetic energy of the inflow) needed to obtain unstable behaviour (flutter or divergence) drops considerably to a minimum at an inflow Mach number near unity. Another important feature of panel flutter is that shock movements play a significant role in the destabilizing mechanisms when the flow is transonic. The fluid equations are 0 1 0 Z Z Z B pn C d B 1 C T _ n ds ¼ U dx dy þ U ðv xÞ B Cds @ A dt XðtÞ pn oXðtÞ oXðtÞ 2 pvT n ð1Þ for arbitrarily moving domains X(t). The conservative state p T Þ is made up of the denvector U = (q, qu, qv, 12 qjvj2 þ c1
Fig. 1. Geometry of the 2D panel flutter problem.
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
sity q, the velocity vector v = (u, v)T and the pressure p. On a point of the boundary oX(t), n = (n1, n2)T is the outward unit normal of X(t) and x_ is the velocity of a boundary point x(t). For the formulation of the panel equation we introduce the space V :¼ H 2 ð0; lÞ \ H 10 ð0; lÞ. The variational formulation is to find the panel deflection w(t, x) satisfying Dðwxx ; uxx Þ þ N ðwÞ ðwx ; ux Þ þ mð€ w; uÞ ¼ ðp2 p1 ; uÞ
8u 2 V ;
ð2Þ
where the nonlinear term models a restoring force due to the panel stretching with Z Eh l 2 w dx. N ðwÞ ¼ 2l 0 x The constants in the panel equation are the stiffness D = E h3/12(1 m2), the mass per unit area m = qsh, the panel thickness h and length l, the density qs, Young’s modulus E and the Poisson ratio m. Now the aeroelastic problem can be formulated. For convenience, all interface and boundary conditions are given in the strong sense. • The fluid domain at time t is defined as the space above the deflecting panel and solid walls and is denoted by XF(t), i.e. for each fluid interface point x(t) associated with a panel point n 2 [0, l] T
xð0Þ ¼ ðn; 0Þ ;
T
xðtÞ ¼ ðn; wðt; nÞÞ .
ð3Þ
• For all moving domains X(t) XF(t) (1) holds. • (2) holds with the r.h.s. determined from the fluid pressure on top of the panel and p1 on the bottom p1 ðt; nÞ pðt; n; wðt; nÞÞ;
n 2 ½0; l;
p2 p1 .
ð4Þ
• At infinity we have q = q1, u = u1, v = 0, p = p1. • On the remaining boundary of XF(t) the fluid velocity in normal direction n equals the velocity of the boundary in direction n, in particular on the fluid–structure interface, _ 2. vT n ¼ x_ T n ¼ wn
ð5Þ
2.1. Energy equation of the aeroelastic problem The energy equations for fluid and structure are given by d EF þ dt Z ¼
Z
1 2 p _ T n ds qv þ ðv xÞ c1 oXF ðtÞ 2 pvT n ds;
oXF ðtÞ
d ES ¼ dt
Z
0
l
ðp2 p1 Þw_ dx;
93
where we have denoted the energy in fluid and structure by Z 1 2 p qv þ dx; EF ðtÞ 2 c 1 XF ðtÞ Z 2 Z l m 2 D 2 Eh 1 l 2 ES ðtÞ w dx . w_ þ wxx dx þ 2 2l 2 0 x 0 2 The energy equation of the fluid is the fourth component of (1) and the energy equation of the structure is derived by plugging u ¼ w_ into (2). Summation of the two energy equations yields the energy equation of the aeroelastic system. Here we restrict the fluid solution to a finite domain XF(t) with non-moving inflow and outflow boundary. Due to the boundary and coupling conditions (3)–(5), we end up with d ðEF þ ES Þ ¼ dt
Z l 1 2 pc _ qv þ p1 wdx. vT nds þ c1 0 oXin=outflow 2
Z
ð6Þ 3. Discretization 3.1. Flow discretization The flow is discretized in space by the finite volume MUSCL scheme on polygonal grids, applying a leastsquares linear reconstruction and the Venkatakrishnan limiter on the primitive fluid state vector. As a numerical flux function F we use the van-Leer flux vector splitting for moving edges. Let us write down the discrete flow equations for an inner cell, when employing the mid-point rule in time X nþ1=2 jXinþ1 jU inþ1 jXni jU ni þ Dt jeij j F
eij oXi nþ1=2 nþ1=2 nþ1=2 nþ1=2 U ij ; U ji ; x_ ij ; nij
¼ 0.
Here, the fluxes are evaluated at the half time grid position xn+1/2 = (xn+1 + xn)/2, with the corresponding normal nþ1=2 nþ1=2 nþ1=2 and length jeij j of the edge eij , which connects nij the cell Xi to its neighbour cell Xj. The edge’s velocity is given by 1 xlnþ1 xnl xnþ1 xnr þ r x_ nþ1=2 ð7Þ 2 Dt Dt and corresponds to the constant velocity by which the edge’s mid-point can be advanced from its position at time level n to its position at time level n + 1. Here, xl and xr are the two end vertices of an edge. The half time geometry xn+1/2 is used to apply reconstruction and limiting on the averaged fluid state Un+1/2 (Un+1 + Un)/2 and as a result nþ1=2 one obtains the state U ij inside the ith cell at the midpoint of the edge eij. With this choice the scheme satisfies its discrete geometric conservation law (DGCL) [7]. At inflow and outflow boundaries non-reflecting conditions according to [20] are implemented. At impermeable boundaries the flux
94
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
0
0
1
B nij;1 C B C jeij j F ðU ij ; ; x_ ij ; nij Þ ¼ jeij jpij B C @ nij;2 A x_ Tij nij
ð8Þ
is prescribed, which results from the slip-condition vT n ¼ x_ T n. Another time stepping scheme, used in the numerical experiments later, is the important second order backward differencing scheme. For this scheme the DGCL is satisfied with the choice of edge normals and velocities as given in [13]. 3.2. Structural discretization Conformal elements in one dimension for fourth order operators are constructed from the Hermite shape functions on the unit interval, n 2 [0, 1], s0;l ðnÞ ¼ ð1 þ 2nÞð1 nÞ2 ; s0;r ðnÞ ¼ ð3 2nÞn2 ; s1;l ðnÞ ¼ nð1 nÞ2 ; s1;r ðnÞ ¼ ð1 nÞn2 . We have n 1 inner nodes, so that 2(n 1) + 2 basis functions /l(x) are constructed from these shape functions when the panel is simply supported. We plug X wðt; xÞ ¼ al ðtÞ/l ðxÞ l
into (2) and test with the basis functions. We introduce bl a_ l to obtain a first order system of ODEs X Dal ð/l ; /k Þ2 þ al N ðwx Þð/l ; /k Þ1 þ mb_ l ð/l ; /k Þ0 l
¼ ðp2 p1 ; /k Þ0 a_ k ¼ bk
8k;
8k.
Here we have set ( )n+1/2 (( )n+1 + ( )n)/2. In particular, the nonlinearity N is evaluated with a trapezoidal rule. The evaluation of the load, ðp1 p2 Þjtnþ1=2 , is not made explicit now, but will be specified P below for the aeroelastic application. DefiningPwn ðxÞ l anl /l ðxÞ and introducing the velocity w_ n ðxÞ l bnl /l ðxÞ, the system reads 9 N nþ1=2 ¼ ðN ðwnþ1 Þ þ N ðwn ÞÞ=2; > > > > = DtDðwnþ1=2 ; /k Þ2 þ DtN nþ1=2 ðwnþ1=2 ; /k Þ1 ð9Þ nþ1 n þ mðw_ w_ ; /k Þ0 ¼ Dtððp2 p1 Þjtnþ1=2 ; /k Þ0 ; > > > > ; nþ1=2 anþ1 ank ¼ Dtbk : k 3.3. Mimicing the energy budget of the continuous problem on the discrete level By definition of the finite volume method as a conservative scheme, the energy conservation properties are carried over from the continuous to the discrete level in the fluid. In the structure, the choice of the time integrator decides if the energy conservation is reflected correctly in the discrete scheme. Here, after the spacial FEM discretization, the time integration scheme (9) is a scheme with the desired property. Furthermore, the use of non-matching interface representations has to be taken care of, so that the energy exchange at the interface is conservative. Here, due to the use of fluid cells with polygonal boundary, a piecewise linear representation of the interface is given for the fluid. Whereas the panel deflection is a piecewise polynomial of third order due to the finite element space employed. And in general, the fluid vertices and the structural nodes do not coincide on the interface, compare Fig. 2. These points will be addressed in the following construction of the aeroelastic scheme. This scheme is constructed such that the energy conservation properties of the continuous problem, compare (6), are mimiced on the discrete level. Using again the mid-point rule in time, the energy balance for a cell Xi can be written as
Here, we have used abbreviations for the bilinear forms ð/; wÞ0 ð/; wÞ; ð/; wÞ1 ð/x ; wx Þ; ð/; wÞ2 ð/xx ; wxx Þ. The set of 4(n 1) + 4 ODEs may now be integrated in time with any kind of Runge–Kutta scheme to approximately solve for discrete values of the coefficients al(t), bl(t). The following second order implicit rule will be used below in the construction of an aeroelastic scheme with correct energy transfer: X nþ1=2 nþ1=2 nþ1=2 Dal ð/l ; /k Þ2 þ al N ð/l ; /k Þ1 l
blnþ1 bnl ð/l ; /k Þ0 ¼ ððp2 p1 Þjtnþ1=2 ; /k Þ0 ; Dt aknþ1 ank nþ1=2 ¼ bk . Dt þm
Fig. 2. Representation of the interface in fluid and structure, the Nj are panel nodes.
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
X
nþ1 EF;i EnF;i ¼ Dt
nþ1=2
jeij
j
eij oXi nþ1=2
F 4 ðU ij
nþ1=2
; U ji
nþ1=2
; x_ ij
nþ1=2
; nij
Þ.
Here F4 is the fourth component of the numerical flux F and all its arguments are taken at the edge center. The energy contained in the ith cell at time level n is EnF;i ¼ n n 2 n jXni j ðq i jvi j =2 þ p i =ðc 1ÞÞ. The discrete fluid energy is P n n EF i EF;i . We can divide the energy production into nþ1=2 the contribution oEFS coming from the interface and nþ1=2
the contribution oEF Thus we have nþ1=2
EFnþ1 EnF ¼ oEFS
from the remaining boundary. nþ1=2
þ oEF
.
The energy flux at impermeable non-moving boundary parts vanishes due to (8). Thus the only contribution nþ1=2 remaining in oEF comes from inflow and outflow boundaries, so that we can write nþ1=2
EFnþ1 EnF ¼ oEFS
nþ1=2
þ oEin=outflow .
ð10Þ
The energy flux at an interface edge, due to (8), is given by nþ1=2
F 4 ðU ij
nþ1=2
; ; x_ ij
nþ1=2
; nij
nþ1=2
Þ pij
nþ1=2 T nþ1=2 Þ nij
ðx_ ij
Dt ðp2 p1 Þjtnþ1=2 ; /k 0 ¼ Dt p1 pðn; wnþ1=2 ðnÞÞ; /k ðnÞ 0 X nþ1=2
Dt Dnij ðp1 pij Þ /k ðnl;ij Þ þ /k ðnr;ij Þ =2. eij CFS
ð12Þ Note that most summands do not contribute anything since in most cases nl,ij, nr,ij 62 supp(/k). For the evaluation of the bilinear products on the l.h.s. of (9) we may use a different quadrature rule, let us assume that this rule is exact. Plugging the quadrature formula (12) into (9), multiplying nþ1=2 the resulting equation with bk ¼ ðaknþ1 ank Þ=Dt and summing over all k, we obtain DtDðwnþ1=2 ; w_ nþ1=2 Þ2 þ DtN nþ1=2 ðwnþ1=2 ; w_ nþ1=2 Þ1 þ mðw_ nþ1 w_ n ; w_ nþ1=2 Þ0
X Dnij ðp1 pijnþ1=2 Þ wnþ1 ðnl;ij Þ wn ðnl;ij Þ ¼ Dt 2 Dt eij CFS nþ1 n w ðnr;ij Þ w ðnr;ij Þ þ . Dt From the third equation of (9) we have wnþ1 wn ¼ Dtw_ nþ1=2 .
and summing over all edges on the interface, we obtain X nþ1=2 nþ1=2 nþ1=2 nþ1=2 nþ1=2 jeij j pij ðx_ ij ÞT nij . oEFS Dt eij CFS
Let us look at a single edge e xl xr CFS with end vertices xl, xr, omitting the ij-indices for the moment, compare Fig. 2. Then x_ nþ1=2 is given by (7). On the interface, we lock the fluid vertices on material points of the panel, demanding that T
xnl ¼ ðnl ; wnl Þ where wnl wn ðnl Þ for some fixed nl 2 ½0; l and correspondingly for the right point xnr and some fixed nr 2 [0, l]. This is the discretization of (3). Thus we get 1 wlnþ1 wnl wrnþ1 wnr nþ1=2 nþ1=2 T nþ1=2 þ Þ n ¼ ðx_ n2 2 Dt Dt on the interface. Looking again at specific edges, adding nþ1=2 nþ1=2 and the indices ij, using jeij j ¼ ðnr;ij nl;ij Þ=n2;ij denoting Dnij nr,ij nl,ij, we obtain X nþ1=2 ¼ Dt Dnij pijnþ1=2 oEFS eij CFS
95
1 wnþ1 ðnl;ij Þ wn ðnl;ij Þ wnþ1 ðnr;ij Þ wn ðnr;ij Þ þ . 2 Dt Dt ð11Þ
We turn to the discrete panel equations (9). We have to specify the load term due to the coupling condition (4). We choose the following quadrature rule to evaluate the load term in (9):
Differentiating s times with respect to x, multiplying by nþ1=2 ak /k , respectively its derivatives, summing over all k and integrating, gives 1 nþ1 ðw wn ; wnþ1 þ wn Þs ¼ Dtðw_ nþ1=2 ; wnþ1=2 Þs 2 for s ¼ 0; 1; 2 or 1 fðwnþ1 ; wnþ1 Þs ðwn ; wn Þs g ¼ Dtðw_ nþ1=2 ; wnþ1=2 Þs 2 for s ¼ 0; 1; 2. fðwnþ1 ; wnþ1 Þ1 þ ðwn ; wn Þ1 g=2 and Note that N nþ1=2 ¼ Eh 2l thus DtN nþ1=2 ðwnþ1=2 ; w_ nþ1=2 Þ1 ¼
Eh ðwnþ1 ; wnþ1 Þ1 þ ðwn ; wn Þ1 ðwnþ1 ; wnþ1 Þ1 ðwn ; wn Þ1 2l 2 2 2
2
Eh ðwnþ1 ; wnþ1 Þ1 ðwn ; wn Þ1 . 2l 4 With the discrete panel energy 2 m n n Eh ðwn ; wn Þ1 D n ES ðw_ ; w_ Þ0 þ þ ðwn ; wn Þ2 2 2l 2 2 ¼
the energy change on the discrete level is X Dnij ðp1 pijnþ1=2 Þ 2 eij CFS
nþ1 w ðnl;ij Þ wn ðnl;ij Þ þ wnþ1 ðnr;ij Þ wn ðnr;ij Þ .
Enþ1 EnS ¼ S
ð13Þ
96
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
Now, from (10), (11), and (13) we see that in the energy change of the discrete fluid–structure system, the contributions from the interface cancel, and EFnþ1
þ
Enþ1 S
EnF
EnS
X Dnij p1
nþ1=2 wnþ1 ðnl;ij Þ wn ðnl;ij Þ ¼ oEin=outflow þ 2 eij CFS nþ1 þ w ðnr;ij Þ wn ðnr;ij Þ .
ð14Þ
results. The following theorem holds: Theorem 1. The numerical scheme P nþ1=2 jXnþ1 jU inþ1 jXni jU ni þ Dt jeij j i
9 > > > > > > > > =
eij oXi
nþ1=2 nþ1=2 nþ1=2 nþ1=2 F ðU ij ; U ji ; x_ ij ; nij Þ
¼ 0;
with the grid velocity (7), the boundary treatment (8), and > > > T T > > > xnl;ij ¼ nl;ij ; wn ðnl;ij Þ ; xnr;ij ¼ nr;ij ; wn ðnr;ij Þ > > ; for all flow edges eij on the panel, 9 > DtDðw ; /k Þ2 þ DtN ðw ; /k Þ1 > > > P nþ1=2 > nþ1 n þ mðw_ w_ ; /k Þ0 ¼ Dt Dnij ðp1 pij Þ > = nþ1=2
nþ1=2
eij CFS
/k ðnl;ij Þ þ /k ðnr;ij Þ =2; nþ1=2
aknþ1 ank ¼ Dtbk
;
ð15Þ
nþ1=2
> > > > > > ; ð16Þ
approximating the aeroelastic initial-boundary value problem of Section 2, carries over the energy balance from the continuous level (6) to the discrete level (14). The discretization satisfies the DGCL in the fluid. Remark 2. The discretization given in the theorem employs, even on the interface, the mid-point rule for the time integration of the fluid equations and a mixed trapezoidal-midpoint rule for the time integration of the structural equations; the nonlinearity N is evaluated with a trapezoidal rule, whereas for the fluid load the mid-point rule is used. We may expect that this scheme is a second order time approximation of our aeroelastic problem. This possible increase in order may be another beneficial point of the scheme apart from its conservation property. It may be the dominating reason for the improvement over conventional fluid–structure methods shown in the numerical experiments below.
4. Algorithms An algorithm that solves the discrete flow equations (15) with a prescribed boundary movement, may serve as a fluid time integrator. Similarly, a structural time integrator is given if the Eq. (16) can be solved when a load is prescribed. In both cases we have employed a matrix-free Newton-GMRes method to solve the discrete equations. With the resulting fluid time integrator and structural time integrator, the aeroelastic problem can be tackled with all partitioned schemes. Amongst the partitioned schemes we have implemented, are the loose coupling, the coupling of Lesoinne and Farhat (denoted from now as the LFcoupling), a predictor scheme and a fixed-point iteration that solves the coupled set (15) and (16). The latter method applies the fluid integrator and the structural integrator several times within one time step. They are successively applied in the sequence given by the loose coupling with the purpose of diminishing the residual of the coupled system (15) and (16). In the sequel we denote this coupling strategy based on the fixed-point iteration as the FPI-coupling. The loose coupling and the LF-coupling are schematically depicted in Fig. 3. In the figure the time axis goes from left to right and the first three time levels are shown. The encircled numbers give the order in which structural time integration (bottom line), fluid grid movement (vertical lines) and fluid time integration (top line) are performed. And we have indicated from which time level the fluid pressure is taken as the load to perform a structural time integration. Furthermore, we can see how the grid position is transferred at the interface from structure (w) to fluid (x). Using again the mid-point rule as the time integrator, we have indicated from which time level the fluid flux is taken. In the loose coupling and the LF-coupling the quadrature rule defining the load transfer from fluid to structure is taken to be exact. This means that the product of the piecewise linear fluid pressure and a structural test function is integrated by placing adequately Gauss quadrature points. We make a remark at this point: Remark 3. We have seen in numerical experiments that when altering the FPI-coupling by replacing the load quadrature in (16) by exact quadrature, the altered FPIcoupling brings results identical to the ones obtained with the original FPI-coupling. It seems that the spacial resolution, necessary to resolve the flow details sufficiently
Fig. 3. Loose coupling (left) and LF-coupling (right). LF-coupling: in step 2, for time level 1, the fluid grid position at the interface is determined by x1 ¼ x0 þ Dtð0; w_ 1=2 ÞT and accordingly, in step 5, for time level 2, by x2 ¼ x1 þ Dtð0; w_ 3=2 ÞT .
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
(e.g. the shock movements, compare Section 5), is too fine to let any differences occur. We stress here that thus the differences we will observe in Section 5 stem from the distinct time discretization used in (15) and (16). Our predictor scheme uses the idea of [12], where the fluid pressure on the interface is predicted in time by a linear extrapolation. This means that our predictor scheme aims to approximate the solution of (15) and (16), i.e. it is precisely defined, when we replace pn+1/2 by 1.5pn 0.5pn1 in (16). In the numerical experiments we denote the scheme as PRED. Another predictor scheme, which is not related to (15) and (16) and is very close to the one proposed in [12], can be obtained by using the second order BDF time integrator in the fluid in conjunction with (16) for the structure, but again replacing pn+1/2 by 1.5pn 0.5pn1. This is a predictor scheme that is not related to a scheme with a design criterion such as PRED is. In the numerical experiments we denote the scheme as PREDBDF. The matrix-free Newton-GMRes solves a nonlinear system G(x) = 0 by applying a Newton iteration and employing GMRes [19] as a linear solver. Since the GMRes iteration only needs to know the evaluation of matrix vec0 0 tor products G (x*)x, with the Jacobian G , the determina0 tion and storage of G may be avoided by approximating G0 ðx Þx
Gðx þ xÞ Gðx Þ .
Here should be chosen big enough to avoid cancellation and small enough to keep the approximation accurate. In all matrix-free Newton-GMRes iterations we have used ¼
eps1=2 jjxjj2
where eps is the machine-accuracy,
which was also recommended in [18]. To accelerate convergence in the linear solvers, a block symmetric Gauss–Seidel preconditioner based on a first order spacial discretization is used in the fluid, whereas a block symmetric Jacobi preconditioner is used in the structure.
97
in the initial guess, then a stagnation of the fluid residual in the following iteration steps can be accepted. 4.2. Grid movement In general applications, the grid movement algorithm, which assigns the vertex positions according to the given boundary movement, may require the solution of a system of equations as-well. In our case, the computational domain has a simple geometry, also when it deforms in the aeroelastic simulation. A simple algebraic rule to stretch the grid suffices and is implemented here. 5. Numerical experiments We mainly examine the behaviour of different coupling schemes with respect to Dt. A sufficiently accurate spacial resolution of 80 fluid cells on the panel and 26 panel nodes was chosen for all computations presented, compare [3,5, 10,12]. In total the fluid grid consists of 28,000 quadrilateral cells, reaching the in- and outflow boundaries at 20 panel lengths away from the panel. Away from the panel, the cells are stretched towards the in- and outflow boundaries. Throughout, we use aluminium as panel material, thus fixing the panel parameters qs = 2700 kg/m3, E = 7.1 · 1010 Pa and m = 0.34. In the dimensionless form all coefficients of the panel equation are determined, when the ratio h/l is specified, which will be done in each of the following cases. The fluid parameters may be defined in two ways: either we assign the inflow Mach number M1, the fluid structure mass ratio l q1l/(qsh) and the dimensionless dynamic inflow pressure k q1 u21 l3 =D or we assign M1 and all remaining flow conditions as those found in the atmosphere at a certain altitude. The panel is taken as undeflected at the initial time, but with a velocity profile, which would cause the panel to oscillate in its first fundamental mode with an amplitude of 0.2 h, if the nonlinearity and the load were not present. Only in Section 5.5 a different initial velocity distribution is used. 5.1. Transonic and supersonic panel flutter
4.1. Convergence criterion for the solution of (15) and (16) As an initial guess in the fixed-point iteration to solve (15) and (16), the result of a loose coupling step can be used. The loose coupling first solves the structural equations with the wrong fluid load and then solves the fluid equations due to the obtained structural deflections. Thus the structural equations of the loose coupling and of (16) do not coincide, whereas the fluid equations of the loose coupling and of (15) do coincide. This means that, when plugging the initial guess into the system (15) and (16), the error is only visible in the residual of the structural equations. Thus, the convergence criterion is to reduce the structural residual below a given threshold. If the residual threshold in the fluid equation is already small enough
Let us first envision some qualitative properties of the flutter solutions. We look at a transonic flutter case with M1 = 0.95, h/l = 0.002, l = 0.1, k = 2392 and a supersonic flutter case defined through h/l = 0.004 and an inflow of M1 = 1.2 at 20,000 feet altitude. The latter case and cases similar to the first one were also considered in [3,5]. The periodic panel deflections resulting in these cases are visualized in Fig. 4 by the deflection at 25%, at 50% and at 75% chord of the panel. Higher harmonics in the panel deflection are more pronounced in the transonic case in contrast to the supersonic case. In Fig. 5, the development of the Mach number distribution on the panel over a flutter cycle is shown. We can observe that the Mach number on the panel has a moving shock in the limit cycle in the transonic case,
98
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102 25 % 50 % 75 %
M = 0.95
10 8
6
25 % 50 % 75 %
M = 1.2
4
6 4
2 w/h
w/h
2 0
0
-2 -2
-4 -6
-4
-8 -10
0
50
100 t / t*
-6 260
270
280
290
300
t / t*
Fig. 4. The panel deflections at various points on the panel.
Fig. 5. Mach numbers along the panel during one period of the limit cycle. M1 = 0.95 on the left, M1 = 1.2 on the right.
whereas in the supersonic case the fluid quantities are smooth along the panel. In the supersonic case shocks periodically build up and disappear at the trailing and the leading edge. In the transonic case we can relate the shock speed to the CFL number in the fluid, and for the given grid, at a CFL number of roughly 100, the shock moves about one cell during a time step, furthermore at this CFL number, we have roughly 60 time steps per flutter period. 5.2. Flutter at M1 = 0.95 Except for k, which we take now as k = 3000, the transonic flutter case considered here is defined with the same data as the one in the previous section. Here we compare the behaviour of the solutions obtained with the LF-coupling and the FPI-coupling, as Dt ! 0. Increasing the fluid CFL number from 10, 30, 50 to 70, we can observe in Fig. 6 how time accuracy is lost in the LF-coupling. We even obtain a ‘‘numerical bifurcation’’ from flutter to divergence with the LF-coupling, when changing from CFL = 50 to CFL = 70. At the latter CFL number the run with the LF-coupling fails at t/t* 50. In contrast, the results obtained with the FPI-coupling remain very accurate even at a CFL number of 70. In comparison with Fig. 7 we see that the predictor method also works quite accurate, but
less accurate then the FPI-coupling. The loose coupling fails to detect a flutter solution at CFL = 50 already. 5.3. Flutter at M1 = 1.2 In this case, for the given fluid grid, a fluid CFL number of 20 corresponds to 56 time steps per flutter period. In Figs. 8 and 9, we see that the loose coupling and the LFcoupling have the tendency to build up the deflection amplitude much earlier than the predictor or the FPIcoupling. A similar tendency can be observed in the loose coupling, the LF-coupling and the predictor coupling when we increase Dt. 5.4. A bifurcation at M1 = 0.95 Another numerical experiment in the transonic regime addresses the determination of a bifurcation point with different coupling schemes and time step sizes. Increasing the dynamic pressure at fixed Mach number a bifurcation from divergence to flutter occurs, where strong shocks move along the panel in the flutter solution. Again an aluminium panel with h/l = 0.002 is considered and we fix M1 = 0.95 and l = 0.1. The non-dimensional dynamic pressure is varied in the range 2000 6 k 6 4000. The bifurcation points
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102 CFL = 10 CFL = 30 CFL = 50 CFL = 70
8
8
6
6
4
4
2
2
0
0
-2
-2
-4
-4
-6
-6
-8
-8
-10
0
50
100
150
CFL = 10 CFL = 30 CFL = 50 CFL = 70
10
w/h
w/h
10
99
-10
0
50
t / t*
100
150
t / t*
Fig. 6. Panel mid-point deflection for LF-coupling and FPI-coupling at various Dt; M1 = 0.95.
CFL = 10 CFL = 30 CFL = 50 CFL = 70
8
8
6
6
4
4
2
2
0
0
-2
-2
-4
-4
-6
-6
-8
-8
-10
-10
0
50
100
150
CFL = 10 CFL = 30 CFL = 50 CFL = 70
10
w/h
w/h
10
0
50
t / t*
100
150
t / t*
Fig. 7. Panel mid-point deflection for loose coupling and predictor coupling at various Dt; M1 = 0.95.
CFL = 10 CFL = 20 CFL = 40
4
4
2
2
0
0
-2
-2
-4
-4
-6 50
100
150
200
t / t*
CFL = 10 CFL = 20 CFL = 40
6
w/h
w/h
6
-6 50
100
150
200
t / t*
Fig. 8. Panel mid-point deflection for loose coupling and predictor coupling at various Dt; M1 = 1.2.
obtained are shown in Fig. 10 and reveal a strong variation depending on the coupling scheme, when using moderate CFL numbers. As we know from Section 5.1, a CFL num-
ber of 50 corresponds to 120 time steps per flutter period for flutter solutions in the considered range of dynamic pressures. Clearly, the loose coupling and the LF-coupling
100
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102 CFL = 10 CFL = 20 CFL = 40
4
4
2
2
0
0
-2
-2
-4
-4
-6 50
100
150
200
CFL = 10 CFL = 20 CFL = 40
6
w/h
w/h
6
-6 50
100
t / t*
150
200
t / t*
Fig. 9. Panel mid-point deflection for LF-coupling and FPI-coupling at various Dt; M1 = 1.2.
4000 LOOSEBDF LOOSE LF PREDBDF PRED FPI
3800 3600 3400
λ
3200 3000 2800 2600 2400 2200 2000 0
10
20
30
40
50
detected around k = 2500. We obtain the diagram of Fig. 11. In contrast to Fig. 10 all schemes now underestimate the bifurcation point more or less at higher CFL numbers. At CFL = 10 there is no entry for the FPI scheme. In this singular case we observe that the fixedpoint iteration does not converge in some of the first 50 time steps. The result is, that the FPI produces only flutter solutions in the k-range around the bifurcation point, although it performs very well at higher CFL numbers. A theoretical argument for such failure in convergence at small Dt is given in [16] in a similar context. As a cure, a relaxation strategy is proposed in [16], which we should also test in the future.
CFL
overestimate the dynamic pressure at the bifurcation point immensely at CFL = 50. The relative error in the dynamic pressure at the bifurcation point is 10% for the loose coupling and 9% for the LF-coupling, both at CFL = 10, whereas the relative errors are 3% for the FPI-coupling and 15% for the predictor scheme, both at CFL = 50. Further, we see that at CFL = 50 the schemes LOOSEBDF and PREDBDF are less accurate than their counterparts LOOSE and PRED. The latter ones can be interpreted as approximations to the conservative method (15) and (16) by employing an estimate for the pressure pn+1/2, thereas the former ones deviate even further from 15 and 16 by using additionally BDF as the time integrator in the fluid instead of the mid-point rule. 5.5. Change of the panel’s initial velocity distribution We now change the initial amplitude of the panel velocity to 0.00005 m/s and keep the other parameters of the above bifurcation problem fixed. In [10] the bifurcation is
5.6. Convergence of FPI We have picked the case with M1 = 0.95 from Section 5.2 at CFL = 50 to exhibit the convergence behaviour of the fixed-point iteration for the solution of (15) and (16). In Fig. 12, the behaviour of the structural residual over a
2600
LOOSE LF PRED FPI
2500
2400
λ
Fig. 10. Bifurcation at M1 = 0.95. For each coupling scheme we find divergence below the corresponding line and flutter above.
2300
2200 0
10
20
CFL
30
40
50
Fig. 11. Bifurcation at M1 = 0.95 with initial panel velocity 0.00005 m/s. For each coupling scheme we find divergence below the corresponding line and flutter above.
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
the fixed-point iteration. From the fluid residual we demand that it is diminished due to the assigned threshold given in the fluid solver. We do not need a decrease of the fluid residual after each fixed-point iteration step. Here, a sufficient threshold in the fluid solver is to reduce the residual below 108. We observe in Fig. 13 that the FPI always keeps the fluid residual below 108.
Log (structural residual)
1
0
-1
-2
6. Conclusions
-3
-4 0
50
100
150 200 Time step
250
300
Fig. 12. Behaviour of structural residual.
range of 300 time steps can be seen. This covers the dimensionless time interval from 0 to 84 in the right plot of Fig. 6, i.e. during these time steps, the solution has already gone through one flutter cycle. After each fluid–structure loop the residual is plotted in the diagram. At each time step we have at least two entries, the one from the initial guess and the one from the following iteration loop. The structural residual of the initial guess, which stems from a loose coupling step, lies around 100.2. This can be red off from the accumulation of dots producing a wavy shape on top of the diagram, whereas after a single FPI iteration the residual is around 101.7, as red off from the wavy shape on the bottom. Throughout, the FPI reduces with one iteration step the structural residual by more than one order of magnitude, which was the convergence criterion for all computations shown for the FPI-coupling. Thus, together with the determination of the initial guess, a numerical simulation with the FPI-coupling needs roughly twice the computational work as the numerical simulation with the loose coupling at the same CFL number. Analogous plots for the fluid residuals can be found in Fig. 13. As explained in Section 4.1, the structural residual is the crucial one, when we judge about the convergence of
The panel flutter problem was considered as a model problem to investigate upon numerical methods for nonlinear aeroelastic problems. We have proposed a discretization of the aeroelastic initial-boundary value problem that is designed such that the energy conservation property of the continuous model is carried over to the discrete model. It employs trapezoidal and mid-point rule as time integrators and deals with non-matching grids at the interface. The DGCL is satisfied in the fluid. The resulting scheme is implicit over the fluid and structure unknowns. We have found that a fixed point iteration (FPI-coupling) solves the system sufficiently, resulting in a fluid–structure algorithm with roughly twice the computational cost of a loose coupling scheme. A singular case, where we encounter convergence problems at small Dt may be dealt with by introducing a relaxation strategy in the future, an idea that is theoretically supported by [16] in a similar context. A predictor scheme derived from the proposed discretization, the scheme proposed in [14] and the loose coupling all have roughly the same computational cost and are compared to the FPI-coupling in transonic aeroelastic cases. We find that the FPI-coupling is still very accurate at moderate to large Dt, i.e. the range of 120 time steps per flutter cycle. We can observe a loss of accuracy at these time step sizes for the predictor scheme, and we observe a severe loss of accuracy at even smaller Dt for the loose coupling and the scheme of [14]. These observations were made by numerical experiments when predicting a bifurcation in the transonic regime, where shock movements play a significant role. References
-6
Log (fluid residual)
101
-8
-10
0
50
100
150 200 Time step
250
Fig. 13. Behaviour of fluid residual.
300
[1] J.J. Alonso, A. Jameson, Fully-implicit time-marching aeroelastic solutions, AIAA-paper 94-0056, 1994. [2] O. Bendiksen, A new approach to computational aeroelasticity, AIAA-paper 91-0939-CP, 1991. [3] O. Bendiksen, G. Davis, Nonlinear travelling wave flutter of panels in transonic flow, AIAA-paper 95-1496, 1995. [4] J.R. Cebral, R. Lo¨hner, Conservative load projection and tracking for fluid–structure problems, AIAA J. 35 (4) (1997) 678–692. [5] G. Davis, O. Bendiksen, Transonic panel flutter, AIAA-paper 931476, 1993. [6] E.H. Dowell, Aeroelasticity of Plates and Shells, Noordhoff International Publishing, 1975. [7] C. Farhat, P. Geuzaine, C. Grandmont, The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids, J. Comput. Phys. 174 (2001) 669–694.
102
R. Massjung / Comput. Methods Appl. Mech. Engrg. 196 (2006) 91–102
[8] C. Farhat, M. Lesoinne, On the accuracy, stability, and performance of the three-dimensional nonlinear transient aeroelastic problems by partitioned procedures, AIAA-paper 96-1388-CP, 1996. [9] C. Farhat, M. Lesoinne, P. LeTallec, Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: momentum and energy conservation, optimal discretization and application to aeroelasticity, Comput. Methods Appl. Mech. Engrg. 157 (1998) 95–114. [10] R.E. Gordnier, M.R. Visbal, Development of a three-dimensional viscous aeroelastic solver for nonlinear panel flutter, AIAA-paper 2000-2337, 2000. [11] G. Guruswamy, Coupled finite-difference/finite-element approach for wing-body aeroelasticity, AIAA-paper 92-4680-CP, 1992. [12] J. Hurka, J. Ballmann, Elastic panels in transonic flow, AIAA-paper 2001-2722, 2001. [13] B. Koobus, C. Farhat, Second-order time-accurate and geometrically conservative implicit schemes for flow computations on unstructured dynamic meshes, Comput. Methods Appl. Mech. Engrg. 170 (1999) 103–129.
[14] M. Lesoinne, C. Farhat, Geometric conservation laws for flow problems with moving boundaries, and their impact on aeroelastic computations, Comput. Methods Appl. Mech. Engrg. 134 (1996) 71–80. [15] S. Morton, R. Melville, M. Visbal, Accuracy and coupling issues of aeroelastic Navier–Stokes solutions on deforming meshes, J. Aircraft 35 (5) (1998) 798–805. [16] P. LeTallec, J. Mouro, Fluid structure interaction with large structural displacements, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3039–3067. [17] S. Piperno, Explicit/implicit fluid/structure staggered procedures with a structural predictor and fluid subcycling for 2D inviscid aeroelastic simulations, Int. J. Numer. Methods Fluids 25 (1997) 1207–1226. [18] N. Qin, D. Ludlow, S. Shaw, A matrix-free preconditioned Newton/ GMRES method for unsteady Navier–Stokes solutions, Int. J. Numer. Methods Fluids 33 (2000) 223–248. [19] Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Publishing, New York, 1996. [20] K. Thompson, Time dependent boundary conditions for hyperbolic systems, J. Comput. Phys. 68 (1987) 1–24.