Uncertainty Quantification in Fluid-Structure Interaction Simulations Using a Simplex Elements Stochastic Collocation Approach Jeroen A.S. Witteveen∗, Hester Bijl† Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629HS Delft, The Netherlands

An efficient uncertainty quantification method for unsteady problems is presented in order to achieve a constant accuracy in time for a constant number of samples. The approach is applied to the aeroelastic problems of a transonic airfoil flutter system and the AGARD 445.6 wing benchmark with uncertainties in the flow and the structure.

I.

Introduction

Numerical errors in multi-physics simulations start to reach acceptable engineering accuracy levels due to the increasing availability of computational resources. Nowadays, uncertainties in modeling multi-scale phenomena in fluid-structure, fluid-thermal, and aero-acoustic interactions have a larger effect on the accuracy of computational predictions than discretization errors. It is therefore especially important in multi-physics applications to systematically quantify the effect of physical variations on a routine basis. Furthermore, unsteady fluid-structure interaction applications are practical aeronautical examples of dynamical systems which are known to amplify initial variations with time. In these problems, natural irreducible input variability can trigger the earlier onset of unstable flutter behavior, which can lead to unexpected fatigue damage and structural failure. Polynomial Chaos uncertainty quantification methods1, 2, 7, 8, 12, 13, 17, 20, 21, 29 however usually result in a fast increasing number of samples with time to resolve the effect of random parameters in dynamical systems with a constant accuracy.16 Resolving the asymptotic stochastic effect, which is of practical interest in postflutter analysis, can in these long time integration problems lead to thousands of required samples. The increasing number of samples is caused by the increasing nonlinearity of the response surface for increasing integration times. This effect is especially profound in problems with oscillatory solutions in which the frequency of the response is affected by the random parameters. The frequency differences between the realizations lead to increasing phase differences with time, which in turn result in an increasingly oscillatory response surface and more required samples. In order to enable efficient uncertainty quantification in time-dependent simulations, a special uncertainty quantification methodology for unsteady oscillatory problems is developed. The approach based on time-independent parameterization of oscillatory samples achieves a constant uncertainty quantification interpolation accuracy in time with a constant number of samples.22 A parameterization in terms of the time-independent functionals frequency, relative phase, amplitude, a reference value, and the normalized period shape is used for period-1 responses. The extension with a damping factor and an algorithm for identifying higher-period shape functions is also applicable to more complex and non-periodic realizations.23 A second uncertainty quantification formulation for achieving a constant accuracy in time with a constant number of samples is developed based on interpolation of oscillatory samples at constant phase instead of at constant time.24 The scaling of the samples with their phase eliminates the effect of the increasing phase ∗ Postdoctoral † Full

researcher, Member AIAA, +31(0)15 2782046, [email protected], http://www.jeroenwitteveen.com. Professor, Member AIAA, +31(0)15 2785373.

1 of 16 American Institute of Aeronautics and Astronautics

differences in the response, which usually leads to the fast increasing number of samples with time. The resulting formulation is not subject to a parameterization error and it can resolve time-dependent functionals that occur for example in transient behavior. It is also proven that phase-scaled uncertainty quantification results in a bounded error as function of the phase for periodic responses and under certain conditions also in a bounded error in time.26 The unsteady approaches are independent of the employed non-intrusive uncertainty quantification to perform the actual interpolation of the oscillatory samples. Here the Simplex Elements Stochastic Collocation (SESC) method based on Newton-Cotes quadrature in simplex elements is employed.25 Since the piecewise polynomial approximation of the response of SESC satisfies the extrema diminishing robustness criterion extended to probability space, it enables us to resolve also bifurcation phenomena reliably.26 The formulation is also extended to multi-frequency responses of continuous structures by using a wavelet decomposition preprocessing step in order to treat the different frequency components separately.28 After the introduction of the mathematical statement of the uncertainty quantification problem in section II, the efficient uncertainty quantification method for unsteady problems is presented in section III. The developed approach is applied to the unsteady fluid-structure interaction test cases of a two-dimensional airfoil flutter problem and the three-dimensional aeroelastic AGARD 445.6 wing. In section IV the stochastic post-flutter analysis of a two-degree-of-freedom rigid airfoil in Euler flow with nonlinear structural stiffness in pitch and plunge with uncertainty in a combination of randomness in two system parameters is considered. The effect of free stream velocity fluctuation is analyzed in section V for the AGARD aeroelastic wing in the transonic flow regime 5% below the deterministic flutter speed. Additional applications to other aeroelastic systems are reported in.18, 19, 27 Finally, the conclusions are summarized in section VI.

II.

Mathematical formulation of the uncertainty quantification problem

Consider a dynamical system subject to na uncorrelated second-order random input parameters a(ω) = {a1 (ω), .., ana (ω)} ∈ A with parameter space A ∈ Rna , which governs an oscillatory response u(x, t, a) L(x, t, a; u(x, t, a)) = S(x, t, a),

(1)

with operator L and source term S defined on domain D × T × A, and appropriate initial and boundary conditions. The spatial and temporal dimensions are defined as x ∈ D and t ∈ T , respectively, with D ⊂ Rd , d = {1, 2, 3}, and T = [0, tmax ]. A realization of the set of outcomes Ω of the probability space (Ω, F, P ) is denoted by ω ∈ Ω = [0, 1]na , with F ⊂ 2Ω the σ-algebra of events and P a probability measure. Here we consider a non-intrusive uncertainty quantification method l which constructs a weighted approximation w(x, t, a) of response surface u(x, t, a) based on ns deterministic solutions vk (x, t) ≡ u(x, t, ak ) of (1) for different parameter values ak ≡ a(ωk ) for k = 1, .., ns . The samples vk (x, t) can be obtained by solving the deterministic problem L(x, t, ak ; vk (x, t)) = S(x, t, ak ),

(2)

for k = 1, .., ns , using standard spatial discretization methods and time marching schemes. A non-intrusive uncertainty quantification method l is then a combination of a sampling method g and an interpolation s method h. Sampling method g defines the ns sampling points {ak }nk=1 and returns the deterministic samples v(x, t) = {v1 (x, t), .., vns (x, t)}. Interpolation method h constructs an interpolation surface w(x, t, a) through the ns samples v(x, t) as an approximation of u(x, t, a). We are eventually interested in an approximation of the probability distribution and statistical moments µui (x, t) of the output u(x, t, a), which can be obtained by sorting and weighted integration of w(x, t, a) Z w(x, t, a)i fa (a)da. (3) µui (x, t) ≈ µwi (x, t) = A

This information can be used for reducing design safety factors and robust design optimization, in contrast to reliability analysis in which the probability of failure is determined.14

III.

An efficient uncertainty quantification method for unsteady problems

The efficient uncertainty quantification formulation for oscillatory responses based on interpolation of scaled samples at constant phase is developed in section B. Scaling the samples with their phase is proven 2 of 16 American Institute of Aeronautics and Astronautics

to result in a bounded error for non-periodic responses in section C. The robust extrema diminishing uncertainty quantification method based on Newton-Cotes quadrature in simplex elements employed in the unsteady approach is first presented in the next section. A.

Robust extrema diminishing uncertainty quantification

A multi-element uncertainty quantification method l evaluates integral (3) by dividing a bounded parameter space A into ne non-overlapping simplex elements Aj ⊂ A µwi (x, t) =

ne Z X j=1

w(x, t, a)i fa (a)da.

(4)

Aj

Here we consider a multi-element Polynomial Chaos method based on Newton-Cotes quadrature points and simplex elements.25 A piecewise polynomial approximation w(x, t, a) is then constructed based on ns deterministic solutions vj,k (x, t) = u(x, t, aj,k ) for the values of the random parameters aj,k that correspond to the n ˜ s Newton-Cotes quadrature points of degree d in the elements Aj µwi (x, t) =

ne X n ˜s X

cj,k vj,k (x, t)i ,

(5)

j=1 k=1

where cj,k is the weighted integral of the Lagrange interpolation polynomial Lj,k (a) through Newton-Cotes quadrature point k in element Aj Z Lj,k (a)fa (a)da, (6) cj,k = Aj

for j = 1, .., ne and k = 1, .., n ˜ s . Here, second degree Newton-Cotes quadrature with d = 2 is considered in combination with adaptive mesh refinement in probability space, since low order approximations are more effective for approximating response response surfaces with singularities. The initial discretization of parameter space A for the adaptive scheme consists of the minimum of neini = na ! simplex elements and nsini = 3na samples, see Figure 1. The example of Figure 1 for two random input parameters can geometrically be extended to higher dimensional probability spaces. The elements Aj are adaptively refined using a refinement measure ρj based on the largest absolute eigenvalue of the Hessian Hj , as measure of the curvature of the response surface approximation in the elements, weighted by the probability fj contained by the elements Z fa (a)da, (7) fj = Aj

Pne

with j=1 fj = 1. The stochastic grid refinement is terminated when convergence measure δne is smaller than a threshold value δne < δ¯ where   k µw⌊ne /2⌋ (x, t) − µwne (x, t) k∞ k σw⌊ne /2⌋ (x, t) − σwne (x, t) k∞ , (8) , δne = max k µwne (x, t) k∞ k σwne (x, t) k∞ with µw (x, t) and σw (x, t) the mean and standard deviation of w(x, t, ω), or when a maximum number of samples n ¯ s is reached. Convergence measure δne can be extended to include also higher statistical moments of the output. In elements where the quadratic second degree interpolation results in an extremum other than in a quadrature point, the element is subdivided into n ˜ e = 2na subelements with a linear first degree NewtonCotes approximation of the response without performing additional deterministic solves. It is proven in26 that the resulting approach satisfies the extrema diminishing (ED) robustness concept in probability space min(w(a)) ≥ min(u(a)) ∧ max(w(a)) ≤ max(u(a)) A

A

A

A

∀u(a),

(9)

where the arguments x and t are omitted for simplicity of the notation. The ED property leads to the advantage that no non-zero probabilities of unphysical realizations can be predicted due to overshoots or undershoots at discontinuities in the response. Due to the location of the Newton-Cotes quadrature points the deterministic samples are also reused in successive refinements and the samples are used in approximating the response in multiple elements. 3 of 16 American Institute of Aeronautics and Astronautics

(a) Element

(b) Initial grid

(c) Adapted grid

Figure 1. Discretization of two-dimensional parameter space A using 2-simplex elements and second-degree Newton-Cotes quadrature points given by the dots.

B.

Efficient uncertainty quantification interpolation at constant phase

Polynomial Chaos methods usually require a fast increasing number of samples with time to maintain a constant accuracy. Performing the uncertainty quantification interpolation of oscillatory samples at constant phase instead of at constant time results, however, in a constant accuracy with a constant number of samples. Assume, therefore, that solving equation (2) for realizations of the random parameters ak results in oscillatory samples vk (t) = u(ak ), of which the phase vφk (t) = φ(t, ak ) is a well-defined monotonically increasing function of time t for k = 1, .., ns . In order to interpolate the samples v(t) = {v1 (t), .., vna (t)} at constant phase, first, their phase as function of time vφ (t) = {vφ1 (t), . . . , vφna (t)} is extracted from the deterministic solves v(t). Second, the time series ˆ (vφ (t)) according for the phase vφ (t) are used to transform the samples v(t) into functions of their phase v to vˆk (vφk (t)) = vk (t), (10) for k = 1, .., ns , see Figure 2. And, third, the sampled phases vφ (t) are interpolated to the function wφ (t, a) wφ (t, a) = h(vφ (t)),

(11)

ˆ (vφ (t)) are interpolated at a constant phase as approximation of φ(t, a). Finally, the transformed samples v ϕ ∈ wφ (t, a) to w(ϕ, ˆ a) = h(ˆ v(ϕ)). (12) Repeating the latter interpolation for all phases ϕ ∈ wφ (t, a) results in the function w(w ˆ φ (t, a), a). The interpolation w(w ˆ φ (t, a), a) is then transformed back to an approximation in the time domain w(t, a) as follows w(t, a) = w(w ˆ φ (t, a), a). (13) The resulting function w(t, a) is an approximation of the unknown response surface u(t, a) as function of time t and the random parameters a(ω). The actual sampling g and interpolation h is performed using the extrema diminishing uncertainty quantification method l based on Newton-Cotes quadrature in simplex elements described in the previous section. The phases vφ (t) are extracted from the samples based on the local extrema of the time series v(t). A trial and error procedure identifies a cycle of oscillation based on two or more successive local maxima. The selected cycle is accepted if the maximal error of its extrapolation in time with respect to the actual sample is smaller than a threshold value ε¯k for at least one additional cycle length. The functions for the phases vφ (t) in the whole time domain T are constructed by identifying all successive cycles of v(t) and linear extrapolation to t = 0 and t = tmax before and after the first and last complete cycle, respectively. The phase is normalized to zero at the start of the first cycle and a user defined parameter determines whether the sample is assumed to attain a local extremum at t = 0. The interpolation at constant phase is restricted to the time domain that corresponds to the range of phases that is reached by all samples in each of the elements. If the phase vφk (t) cannot be extracted from one of the samples vk (t) for k = 1, .., ns , then uncertainty quantification interpolation h is directly applied to the time-dependent samples v(t).

4 of 16 American Institute of Aeronautics and Astronautics

v^ k

vk

phase φ

time t (a) samples vk (t)

(b) samples vˆk (φ)

Figure 2. Oscillatory samples as function of time and phase.

C.

Error bound for non-periodic responses

It is shown below that the uncertainty quantification interpolation of non-periodic samples at constant phase results in a bounded relative error. Let u(t, a) be a non-periodic response as function of time t for t ∈ R given by u(t, a) = A(t, a)up (t, a), (14) with amplitude A(t, a) > 0 ∀t ∈ R and up (t, a) a periodic function in time up (t + zT (a), a) = up (t, a),

for all z ∈ Z and a ∈ A,

(15)

with T (a) = 1/f (a) > 0 the period length and f (a) the frequency affected by the random input a(ω). The phase φ(t, a) of the periodic part up (t, a) is given by φ(t, a) = φ0 (a) +

t , T (a)

(16)

with φ0 (a) = φ(0, a). Consider non-periodic response u(t, a) that can then be written as function of φ(t, a) as follows ˆ u(t, a) = u ˆ(φ(t, a), a) = A(φ(t, a))ˆ up (φ(t, a), a), (17) ˆ with A(φ(t, a)) = A(t, a) and u ˆp (φ(t, a), a) = up (t, a). Consider uncertainty quantification method l which results in an approximation w(t, a) of u(t, a) based on applying interpolation method h at constant phase to ns samples v(t) = {v1 (t), . . . , vns } for parameter values ak for k = 1, . . . , ns resulted from sampling method g. Let h be a linear interpolation method in the sense that ch(v) = h(cv),

(18)

with c a constant, which holds for virtually all interpolation techniques including the Simplex Elements Stochastic Collocation approach. Theorem 1. Relative error εˆ(ϕ, a) in approximation w(ϕ, ˆ a) with respect to non-periodic response surface u ˆ(ϕ, a) given by |w(ϕ, ˆ a) − u ˆ(ϕ, a)| εˆ(ϕ, a) = , (19) ˆ A(ϕ) ˆ with A(φ(t, a)) = A(t, a), as results from uncertainty quantification method l applied at constant phase ϕ is bounded for all ϕ ∈ R and a ∈ A by δ for which holds εˆ(ϕ, a) < δ,

for all ϕ ∈ [0, 1] and a ∈ A.

(20)

Proof. For the periodic part up (t, a) of the non-periodic response u(t, a) = A(t, a)up (t, a) holds up (t, a) = u ˆp (φ(t, a), a), 5 of 16 American Institute of Aeronautics and Astronautics

(21)

and   t + zT (a) u ˆp (φ + z, a) = uˆp φ0 (a) + , a = up (t + zT (a), a) T (a) = up (t, a) = uˆp (φ, a), for all z ∈ Z and a ∈ A.

(22)

As function of the phase φ also holds ˆ up (φ, a), u ˆ(φ, a) = A(φ)ˆ

(23)

ˆ with A(φ) = A(t, a). The samples vk (t) resulting from sampling method g vk (t) = gk (u(t, a)) = u(t, ak ),

(24)

for k = 1, . . . , ns can also be written in terms of a periodic part vpk (t) vk (t) = A(t, ak )vpk (t),

(25)

with vpk (t) = gk (up (t, a)) = up (t, ak ). The signals vpk (t) are periodic with period length vTk = T (ak ), since using (15) vpk (t + zvTk ) = up (t + zT (ak ), ak ) = up (t, ak ) = vpk (t) for all z ∈ Z, (26) for k = 1, . . . , ns . The phase vφk (t) = φ(t, ak ) of the signals vpk (t) is then in correspondence with (16) given by t , (27) vφk (t) = vφ0k + vTk for k = 1, . . . , ns , with vφ0k = vφk (0). Scaling the functions vpk (t) with their phase vφk (t) results in vpk (t) = vˆpk (vφk (t)), for k = 1, . . . , ns . Periodicity of vpk (t) gives   t + zvTk vˆpk (vφk (t) + z) = vˆpk vφ0k + vTk = vpk (t + zvTk ) = vpk (t) = vˆpk (vφk (t)),

(28)

for all z ∈ Z,

(29)

for k = 1, . . . , ns . In terms of the phase vφk (t), (25) becomes ˆ φ (t))ˆ vˆk (vφk (t)) = A(v vpk (vφk (t)), k

(30)

ˆ φ (t)). Uncertainty quantification method l results in approxiwith vk (t) = vˆk (vφk (t)) and A(t, ak ) = A(v k ˆ (ϕ) at a constant phase ϕ mation w(ϕ, ˆ a) by applying interpolation method h of the samples v w(ϕ, ˆ a) = h(ˆ v(ϕ)) = h(ˆ v1 (ϕ), . . . , vˆns (ϕ)).

(31)

Introduce also the following notations

and

w ˆp (ϕ, a) = h(ˆ vp1 (ϕ), . . . , vˆpns (ϕ)),

(32)

ˆ wˆp (ϕ, a). w(ϕ, ˆ a) = A(ϕ)

(33)

The relative error εˆ(ϕ, a) as function of phase ϕ in approximation w(ϕ, ˆ a) with respect to u ˆ(ϕ, a) is defined by (19) as |w(ϕ, ˆ a) − u ˆ(ϕ, a)| εˆ(ϕ, a) = . (34) ˆ A(ϕ)

6 of 16 American Institute of Aeronautics and Astronautics

Then holds for εˆ(ϕ + z, a) the following |w(ϕ ˆ + z, a) − u ˆ(ϕ + z, a)| ˆ A(ϕ + z) h(ˆ ˆ + z)ˆ up (ϕ + z, a) v1 (ϕ + z), . . . , vˆns (ϕ + z)) A(ϕ − = ˆ + z) ˆ + z) A(ϕ A(ϕ ! ˆ + z)ˆ ˆ + z)ˆ A(ϕ vpns (ϕ + z) A(ϕ vp1 (ϕ + z) ,..., − uˆp (ϕ + z, a) = h ˆ + z) ˆ + z) A(ϕ A(ϕ

εˆ(ϕ + z, a) =

ˆp (ϕ, a)| = |h(ˆ vp1 (ϕ), . . . , vˆpns (ϕ)) − u ˆ wˆp (ϕ, a) − A(ϕ)ˆ ˆ up (ϕ, a)| |A(ϕ) = ˆ A(ϕ) =

|w(ϕ, ˆ a) − u ˆ(ϕ, a)| = εˆ(ϕ, a), ˆ A(ϕ)

for all z ∈ Z and ϕ ∈ R and a ∈ A.

(35)

Error εˆ(ϕ, a) is, therefore, a periodic function of ϕ. Define δ for which holds (20) εˆ(ϕ, a) < δ,

for all ϕ ∈ [0, 1] and a ∈ A,

then holds εˆ(ϕ + z, a) = εˆ(ϕ, a) < δ,

for all z ∈ Z and ϕ ∈ [0, 1] and a ∈ A,

(36)

for all ϕ ∈ R and a ∈ A,

(37)

and εˆ(ϕ, a) < δ,

Relative error εˆ(ϕ, a) in approximation w(ϕ, ˆ a) is, therefore, bounded by δ for all ϕ ∈ R and a ∈ A. Notice that the proof of Theorem 1 is independent of uncertainty quantification method l, sampling method g, and interpolation method h, as long as the interpolation method satisfies (18), ch(v) = h(cv). For the special case of periodic samples with A(t, a) a constant, the error bound εˆ(ϕ, a) < δ holds for the absolute error εˆ(ϕ, a) = |w(ϕ, ˆ a) − u ˆ(ϕ, a)|.26 ˆ For non-periodic responses which cannot be written as u(t, a) = A(t, a)up (t, a), with A(t, a) = A(φ(t, a)), the relative error εˆ(ϕ, a) is not strictly bounded by δ. However, scaling the resulting samples vk (t) with their phase still eliminates the effect of the increasing phase differences on the increase of the number of required samples. Also for these responses, performing the uncertainty quantification interpolation at constant phase results practically in a constant accuracy with a constant number of samples. The bounded error εˆ(ϕ, a) as function of phase ϕ also results in a bounded error ε(t, a) in time for the Simplex Elements Stochastic Collocation method l, if initial phase φ0 (a) and frequency f (a) of periodic function up (t, a) depend linearly on a, with u(t, a) = A(t, a)up (t, a). Let initial phase φ0 (a), therefore, depend linearly on the random parameters a φ0 (a) = cφ0 ,0 + cφ0 ,1 · a.

(38)

where · denotes the vector inner product, with cφ0 ,0 constant and cφ0 ,1 a vector containing na constants. And let frequency f (a) also depend linearly on a(ω) f (a) = cf,0 + cf,1 · a,

(39)

with cf,0 constant and cf,1 an na -dimensional constant vector. Further is used that Simplex Elements Stochastic Collocation exactly interpolates linear functions. Theorem 2. Relative error ε(t, a) in approximation w(t, a) with respect to non-periodic response surface u(t, a) given by |w(t, a) − u(t, a)| , (40) ε(t, a) = A(t, a) ˆ with A(t, a) = A(φ(t, a)), as resulted from Simplex Elements Stochastic Collocation l applied at constant phase ϕ is bounded for all t ∈ R and a ∈ A by δ for which holds (20) εˆ(ϕ, a) < δ,

for all ϕ ∈ [0, 1] and a ∈ A, 7 of 16

American Institute of Aeronautics and Astronautics

if initial phase φ0 (a) and frequency f (a) of periodic function up (t, a) depend linearly on a, with u(t, a) = A(t, a)up (t, a). Proof. The phase φ(t, a) of the periodic part up (t, a) of non-periodic response u(t, a) is given by (16) φ(t, a) = φ0 (a) +

t = φ0 (a) + f (a)t. T (a)

The linear dependence of φ0 (a) and f (a) on a given by (38) and (39) results in φ(t, a) = cφ0 ,0 + cf,0 t + (cφ0 ,1 + cf,1 t) · a

(41)

The ns sampled phases vφ (t) = {vφ1 (t), . . . , vφns (t)} resulting from sampling method g are, therefore, vφk (t) = cφ0 ,0 + cf,0 t + (cφ0 ,1 + cf,1 t) · ak ,

(42)

for k = 1, . . . , ns . The resulting wφ (t, a) of Simplex Elements Stochastic Collocation interpolation h of the samples vφ (t) then exactly reconstructs the function φ(t, a) wφ (t, a) = h(vφ (t)) = cφ0 ,0 + cf,0 t + (cφ0 ,1 + cf,1 t) · a = φ(t, a).

(43)

Therefore, error εˆ(φ(t, a), a) in the approximation w(w ˆ φ (t, a), a) of response u ˆ(φ(t, a), a) becomes εˆ(φ(t, a), a)

=

|w(w ˆ φ (t, a), a) − u ˆ(φ(t, a), a)| ˆ A(φ(t, a))

=

|w(φ(t, ˆ a), a) − u ˆ(φ(t, a), a)| ˆ A(φ(t, a))

|w(t, a) − u(t, a)| A(t, a) = ε(t, a) < δ for all φ ∈ R and a ∈ A,

=

(44)

according to Theorem 1. Using (16) gives ε(t, a) < δ

IV.

for all t ∈ R and a ∈ A.

(45)

Transonic airfoil flutter

The combined effect of independent randomness in the ratio of natural frequencies ω ¯ (ω) and the free stream velocity U∞ (ω) on the post-flutter behavior of an elastically mounted airfoil is analyzed. The structural model of the pitch-plunge airfoil with cubic nonlinear spring stiffness is given by:6, 11  ω ¯ 2 1 (ξ + βξ ξ 3 ) = − Cl (τ ), U∗ πµ

(46)

1 2 xα ′′ ξ + α′′ + ∗2 (α + βα α3 ) = Cm (τ ), rα2 U πµrα2

(47)

ξ ′′ + xα α′′ +

where βξ = 0m−2 and βα = 300rad−2 are the cubic spring parameters, ξ(τ ) = h/b is the non-dimensional plunge displacement of the elastic axis, see Figure 3, α(τ ) is the pitch angle, and (′ ) denotes differentiation with respect to non-dimensional time τ = U t/b, with half-chord length b = c/2 = 0.5m. The radius of gyration around the elastic axis is rα b = 0.25m, bifurcation parameter U ∗ is defined as U ∗ = U/(bωα ), and the airfoil-air mass ratio is µ = m/πρ∞ b2 = 100, with m the airfoil mass. The elastic axis is located at a distance ah b = −0.25m from the mid-chord position and the mass center is located at a distance xα b = 0.125m from the elastic axis. The ratio of natural frequencies is defined as ω(ω) = ωξ /ωα , with ωξ and ωα the natural frequencies of the airfoil in pitch and plunge, respectively. The randomness in ω ¯ (ω) is described by a uniform distribution around mean value µω¯ = 0.25 with a coefficient of variation of 10%. The 8 of 16 American Institute of Aeronautics and Astronautics

c

b reference position

h

α mid−chord elastic axis centre of mass

ah b x αb

Figure 3. The elastically mounted pitch-plunge airfoil model.

free stream velocity U∞ (ω) is subject to a symmetric unimodal beta distribution with β1 = β2 = 2 with a coefficient of variation of 1% around mean µU∞ = 276.27m/s, which corresponds to M∞ = 0.8. The non-dimensional aerodynamic lift and moment coefficients, Cl (τ ) and Cm (τ ), are determined by solving the unsteady Euler equations. The two-dimensional flow domain is discretized by an unstructured hexahedral mesh of 12 · 103 cells, which was selected based on a grid convergence study. The governing equations are discretized using a second order central finite volume discretization stabilized with artificial dissipation. An Arbitrary Lagrangian-Eulerian formulation is employed to couple the fluid mesh with the movement of the structure. The fluid mesh is deformed using radial basis function interpolation of the boundary displacements.4 Time integration is performed using the second order BDF-2 method until t = 3 with time step ∆t = 0.002, which was established after a time step refinement study. Initially the airfoil is at rest at a deflection of α(0) = 0.1deg and ξ(0) = 0 from its equilibrium position. In order to study the postbifurcation behavior, the bifurcation parameter U ∗ is fixed at 130% of the deterministic linear bifurcation point for the mean values of the random parameters. The stochastic behavior of the angle of attack α(t, ω) is resolved as indicator for the post-flutter airfoil behavior. The Unsteady Adaptive Stochastic Finite Elements response surface approximation of the angle of attack α(t, ω) as function of the random parameters ω ¯ (ω) and U∞ (ω) at t = {0.5; 1.5; 2.5} given in Figure 4 shows an increasingly oscillatory response surface with time. The 10% variation in ω ¯ (ω) has a larger effect on the frequency of the response than U∞ (ω) with 1% variation. Both parameters have a small effect on the amplitude of the oscillation of α(t, ω) of approximately 3o . At t = 0.5 the airfoil exhibits transient behavior from its initial perturbation of α(0) = 0.1o , which is indicated by the smaller amplitude of the response surface variations of approximately 2o . These results are obtained using the time-independent grid in probability space shown in Figure 4d with ns = 9 samples, ne = 2 elements, and nesub = 4096 post-processing subelements. The resulting UASFE approximation of the mean µα (t) and standard deviation σα (t) of the angle of attack α(t, ω) in Figure 5 shows two frequency signals due to the effect of the two random parameters on the frequency of the response. The mean µα (t) exhibits initially an increasing oscillation caused by the deterministic transient of the samples, after which it develops a decaying oscillation due to the effect of the random parameters on the frequency of the response. The large effect of the random parameters on the dynamical system is illustrated by the fast initial increase of the standard deviation σα (t) from its deterministic initial condition. Although the deterministic post-flutter behavior is highly unsteady, the stochastic response reaches a steady asymptotic behavior with a standard deviation of σα = 1.6o , which is a factor 16 larger than the initial angle of attack α(0) = 0.1o . The discretizations with ns = {9, 13, 25} samples and ne = {2, 4, 8} uniformly refined elements, respectively, indicate that the results are uniformly converged in time. The approximation with ns = 25 is converged up to δne = 6.2·10−3 , where δne is defined by (8). The local convergence for µα (t) and σα (t) at t = {0.5; 1.0; 1.5; 2.0; 2.5} given in Tables 1 and 2 for ns = {13, 25} shows no clear increase of convergence measure δ with time. This illustrates that the convergence and the accuracy of the UASFE approximation are in practice constant in time.

9 of 16 American Institute of Aeronautics and Astronautics

(a) t = 0.5

(b) t = 1.5

(c) t = 2.5

(d) Stochastic grid

Figure 4. Response surface of angle of attack α(ω) as function of random natural frequency ratio ω ¯ (ω) and free stream velocity U∞ (ω) for transonic airfoil flutter.

1.8

n =9 s

0.6

n =13

1.6

s

n =25 s

standard deviation σα [deg]

α

mean µ [deg]

0.4 0.2 0 −0.2 −0.4

1.2 1 0.8 0.6 ns=9

0.4

n =13 s

0.2

−0.6 0

1.4

n =25 s

0.5

1

1.5 time t [s]

2

2.5

0 0

3

(a) Mean

0.5

1

1.5 time t [s]

2

2.5

3

(b) Standard deviation

Figure 5. Mean and standard deviation of angle of attack α(ω) for transonic airfoil flutter with random natural frequency ratio ω ¯ (ω) and free stream velocity U∞ (ω).

10 of 16 American Institute of Aeronautics and Astronautics

Table 1. Convergence measure δne for mean angle of attack α(t, ω) for transonic airfoil flutter with random natural frequency ratio ω ¯ (ω) and free stream velocity U∞ (ω).

ns 13 25

ne 4 8

t = 0.5 0.640 · 10−3 0.268 · 10−3

t = 1.0 3.712 · 10−3 2.455 · 10−3

t = 1.5 4.426 · 10−3 3.138 · 10−3

t = 2.0 7.207 · 10−3 4.422 · 10−3

t = 2.5 4.331 · 10−3 2.684 · 10−3

Table 2. Convergence measure δne for the standard deviation of angle of attack α(t, ω) for transonic airfoil flutter with random natural frequency ratio ω ¯ (ω) and free stream velocity U∞ (ω).

ns 13 25

ne 4 8

t = 0.5 2.943 · 10−3 0.973 · 10−3

V.

t = 1.0 3.275 · 10−3 4.388 · 10−3

t = 1.5 7.500 · 10−3 0.859 · 10−3

t = 2.0 5.378 · 10−3 2.194 · 10−3

t = 2.5 3.896 · 10−3 3.344 · 10−3

Three-dimensional transonic wing

The transonic AGARD 445.6 wing30 is a standard benchmark case for the fluid-structure interaction of a three-dimensional continuous structure. The discretization of the aeroelastic configuration is described in section A. In section B randomness is introduced in the free stream velocity. The stochastic response of the system and the flutter probability are determined. A.

AGARD 445.6 wing benchmark problem

The AGARD aeroelastic wing configuration number 330 known as the weakened model is considered here with a NACA 65A004 symmetric airfoil, taper ratio of 0.66, 45o quarter-chord sweep angle, and a 2.5-foot semi-span subject to an inviscid flow. The structure is described by a nodal discretization using an undamped linear finite element model in the Matlab finite element toolbox OpenFEM .15 The discretization contains in the chordal and spanwise direction 6 × 6 brick-elements with 20 nodes and 60 degrees-of-freedom, and at the leading and trailing edge 2 × 6 pentahedral elements with 15 nodes and 45 degrees-of-freedom as in.33 The orthotropic material properties are obtained from3 and the fiber orientation is taken parallel to the quarter-chord line. The Euler equations for inviscid flow5 are solved using a second-order central finite volume discretization on a 60 × 15 × 30m domain using an unstructured hexahedral mesh. The free stream conditions for the density ρ∞ = 0.099468kg/m3 and pressure p∞ = 7704.05Pa are taken from.30 Time integration of the samples is performed using a third-order implicit Runge-Kutta scheme9 until t = 1.25s to determine the stochastic solution until t = 1s. The first bending mode with a vertical tip displacement of ytip = 0.01m is used as initial condition for the structure, see Figure 6. The coupled fluid-structure interaction system is solved using a partitioned IMEX scheme31, 32 with explicit treatment of the coupling terms without sub-iterations. An Arbitrary Lagrangian-Eulerian formulation is employed to couple the fluid mesh with the movement of the structure. The flow forces and the structural displacements are imposed on the structure and the flow using nearest neighbor and radial basis function interpolation,33 respectively. The fluid mesh is also deformed using radial basis function interpolation of the boundary displacements.4 A convergence study has been performed to determine a suitable flow mesh discretization and time step size. Deterministic results for the selected flow mesh with 3.1 · 104 volumes and time step of ∆t = 2.5 · 10−3s agree well with experimental and computational results in the literature.10, 30, 33 The deterministic flutter velocity is found to be Uflut = 313m/s, which corresponds to a Mach number of M∞ = 0.951. B.

Randomness causes non-zero flutter probability

In the following, the effect of randomness in the free stream velocity U∞ (ω) is studied. The mean free stream velocity is chosen 5% below the actual deterministic flutter velocity, µU∞ = 0.95Uflut, to assess the effectiveness of a realistic design safety factor. The coefficient of variation of the assumed unimodal

11 of 16 American Institute of Aeronautics and Astronautics

Figure 6. Initial condition and grid for the AGARD 445.6 wing for mean free stream velocity µU∞ .

beta distribution is set to cvU∞ = 3.5%. The outputs of interest are the lift L(t, ω) and the vertical tip displacement of the tip-node ytip (t, ω). The first Ns = 3 sampled time series of the lift Li (t, ω) of the UASFE discretization with Ne = 1 element show in Figure 7a that the first bending mode is the dominant mode in the system response. A second mode which is initially present in the response, damps out quickly, such that a wavelet decomposition pre-processing step is in this case not necessary to obtain the stochastic solution using UASFE. The samples illustrate that the free stream velocity has a significant effect on the frequency and the damping of the system response, which results in a diverging oscillation for i = 3, and decaying oscillations for i = 1 and mean value µU∞ at i = 2. The same conclusions can be drawn from Figure 7b in which the response surface approximation of the lift L(t, ω) at t = 1 is given for Ne = 5 elements and Ns = 11 samples. The response surface has an oscillatory character due to the effect of the random U∞ (ω) on the frequency of the lift oscillation and consequently on the phase differences in L(t, ω) at t = 1. The adaptive UASFE grid refinement results automatically in a gradually finer mesh in the region of large lift amplitudes at large values of U∞ (ω). 250

250 200

UASFE (N =1,N =3) e

t=1

N =11 samples 200

150 100

s

150

lift L [N]

50

lift [N]

UASFE (N =5) e

s

0

100

−50 50

−100 −150

0

−200 −250 0

i=1

i=2 i=3 0.2

0.4

time t

0.6

0.8

−50 0.88

1

0.9

0.92

0.94 0.96 U /U ∞

(a) lift samples Li (t)

0.98

1

1.02

flut

(b) lift response surface L(t, ω) at t = 1

Figure 7. Results for the AGARD 445.6 wing with random free stream velocity U∞ (ω).

Results for the time evolution of the mean µL (t) and the standard deviation σL (t) of the lift are given in Figure 8 for Ne = 4 and Ne = 5 elements. The two approximations are converged with respect to each other up to 5 · 10−3 . The time history for the mean lift µL (t) shows a decaying oscillation up to t = 0.4s

12 of 16 American Institute of Aeronautics and Astronautics

from the initial value of µL = −23.9N. This behavior can be explained by the decaying lift oscillation for a large range of U∞ (ω) values and the effect of U∞ (ω) on the increasing phase differences with time. For t > 0.4 the decay is approximately balanced by the exponentially increasing amplitude of the unstable part of the U∞ (ω) parameter domain. In contrast, the standard deviation shows an oscillatory increase from the initial σL = 2.46N up to a local maximum of σL = 18.3N at t = 0.31s due to the increasing phase differences with time. For t > 0.31 the standard deviation slightly decreases due to the decreasing lift amplitude in part of the parameter domain. Eventually, the unstable realizations result in an increasing standard deviation which reaches at t = 1 values between σL = 14 and σL = 19, which corresponds to an amplification of the initial standard deviation with a factor 6 to 8. 40

25

UASFE (N =5,N =11) e

s

UASFE (N =4,N =9)

30

s

20

standard deviation lift σ [N]

e

10

L

mean lift µ [N]

L

20

0 −10 −20

15

10

5 UASFE (N =5,N =11)

−30

e

s

UASFE (N =4,N =9) e

−40 0

0.2

0.4

time t [s]

0.6

0.8

0 0

1

0.2

(a) mean µL

0.4

time t [s]

0.6

s

0.8

1

(b) standard deviation σL

Figure 8. Results for the AGARD 445.6 wing with random free stream velocity U∞ (ω).

The nodal description of the structure directly returns the vertical tip-node displacement ytip (t, ω). The approximations of the mean µytip (t) and standard deviation σytip (t) of ytip (t, ω) show in Figure 9 a qualitatively similar behavior as the lift L(t, ω). The standard deviation σytip (t) vanishes, however, initially due to the deterministic initial condition for the structure in contrast with the non-zero σL (t) at t = 0. The standard deviation reaches values between σytip = 4.2 · 10−3 m and σytip = 5.6 · 10−3 m at t = 1, which corresponds to a standard deviation equal to 42% and 56% of the deterministic initial vertical tip deflection.

−3

0.015

7

e

s

y,tip

e

standard deviation tip displacement σ

s

mean tip displacement µ

y,tip

[m]

UASFE (N =4,N =9) 0.01

0.005

0

−0.005

−0.01

−0.015 0

0.2

0.4

time t [s]

0.6

0.8

x 10

[m]

UASFE (N =5,N =11)

6 5 4 3 2

(a) mean µytip

e

s

UASFE (N =4,N =9) e

0 0

1

UASFE (N =5,N =11)

1

0.2

0.4

time t [s]

0.6

0.8

(b) standard deviation σytip

Figure 9. Results for the AGARD 445.6 wing with random free stream velocity U∞ (ω).

13 of 16 American Institute of Aeronautics and Astronautics

s

1

The probability of flutter can be determined by constructing the probability distribution of the damping factor of the system given in Figure 10. The damping factor is here extracted from the last period of oscillation of the sampled vertical tip node displacements. Positive and negative damping factors denote unstable and damped oscillatory responses, respectively. Even though the mean free stream velocity µU∞ is fixed at a safety margin of 5% below the deterministic flutter velocity Uflut , the non-zero probability of positive damping indicates a non-zero flutter probability. The 3.5% variation in U∞ (ω) results actually in a probability of flutter of 6.19%. Taking physical uncertainties into account in numerical predictions is, therefore, a more reliable approach than using safety margins in combination with deterministic simulation results. 1 0.9

probability distribution

0.8 0.7 0.6 0.5 0.4 0.3 0.2

UASFE (N =5,N =11) e

0.1

e

0 −4

s

UASFE (N =4,N =9) −3

−2

s

−1 damping factor

0

1

2

Figure 10. Results for the AGARD 445.6 wing with random free stream velocity U∞ (ω).

VI.

Conclusions

An uncertainty quantification formulation for achieving a constant accuracy in time with a constant number of samples in multi-physics aeroelastic problems is developed based on interpolation of oscillatory samples at constant phase instead of at constant time. The scaling of the samples with their phase eliminates the effect of the increasing phase differences in the response, which usually leads to the fast increasing number of samples with time. The resulting formulation is not subject to a parameterization error and it can resolve time-dependent functionals that occur for example in transient behavior. Phase-scaled uncertainty quantification results in a bounded error as function of the phase for periodic responses and under certain conditions also in a bounded error in time. The developed approach is applied to the unsteady fluid-structure interaction test cases of a twodimensional airfoil flutter problem and the three-dimensional aeroelastic AGARD 445.6 wing. In the stochastic post-flutter analysis of a two-degree-of-freedom rigid airfoil in Euler flow with nonlinear structural stiffness in pitch and plunge, a combination of randomness in two system parameters is considered. The uncertainty in the ratio of natural frequencies of the structure is described by uniform distribution around mean value of 0.25 with a coefficient of variation of 10%. The free stream velocity is subject to a symmetric unimodal beta distribution with a coefficient of variation of 1% around a mean of 276.27m/s, which corresponds to a Mach number of 0.8. In order to study the post-bifurcation behavior, the bifurcation parameter is fixed at 130% of the deterministic linear bifurcation point for the mean values of the random parameters. The behavior of the mean, standard deviation, and response surface of the angle of attack is resolved as indicator for the post-flutter airfoil behavior. The results show that the frequency of the increasingly oscillatory response surface with time is affected more by the frequency ratio than the free stream velocity. Both parameters have a small effect on the amplitude of the oscillation. The asymptotic stochastic behavior of the time-dependent problem is steady with a standard deviation of 1.6 degrees, which is a factor 16 larger than the deterministic initial angle of attack. Convergence results for discretizations with increasing number of samples indicate that the applied uncertainty quantification methodology results in practice in a time-independent accuracy with a constant number of samples.

14 of 16 American Institute of Aeronautics and Astronautics

The AGARD aeroelastic wing configuration number 3 known as the weakened model is also considered with a NACA 65A004 symmetric airfoil, taper ratio of 0.66, 45 degrees quarter-chord sweep angle, and a 2.5-foot semi-span subject to an inviscid flow. The free stream conditions for the density and pressure are 0.099468 kg/m3 and 7704.05 Pa, respectively. The first bending mode with a vertical tip displacement of 0.01 m is used as initial condition for the structure. The deterministic flutter velocity for this configuration is found to be 313 m/s, which corresponds to a Mach number of 0.951. The effect of randomness in the free stream velocity is studied. The mean free stream velocity is chosen 5% below the actual deterministic flutter velocity to assess the effectiveness of a realistic design safety factor. The coefficient of variation of the assumed unimodal beta distribution is set to 3.5%. The outputs of interest are the mean and standard deviation of the lift and the vertical tip displacement of the tip-node, the probability distribution of the damping coefficient, and the probability of flutter. Even though the mean free stream velocity is fixed at a safety margin of 5% below the deterministic flutter velocity, a resulting non-zero probability of positive damping indicates a non-zero flutter probability. The 3.5% variation in free stream velocity results actually in a probability of flutter of 6.19%. Taking physical uncertainties into account in numerical predictions is, therefore, a more reliable approach than using safety margins in combination with deterministic simulation results.

References 1 I. Babuˇ ska, R. Tempone, G.E. Zouraris, Galerkin finite element approximations of stochastic elliptic partial differential equations, SIAM J. Numer. Anal. 42 (2004) 800–825. 2 I. Babuˇ ska, F. Nobile, and R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (2007) 1005–1034. 3 R.Beaubien, F. Nitzsche, D. Feszty, Time and frequency domain solutions for the AGARD 445 wing, in: Proceedings of the International Forum on Aeroelasticity and Structural Dynamics (IFASD), Munich, Germany, 2005. 4 A. de Boer, M.S. van der Schoot, H. Bijl, Mesh deformation based on radial basis function interpolation, Comput. Struct. 85 (2007) 784–795. 5 A.J. Chorin, J.E. Marsden, A mathematical introduction to fluid mechanics, Springer-Verlag, New York, 1979. 6 Y. Fung, An introduction to aeroelasticity, Dover Publications, New York, 1969. 7 R.G. Ghanem, P. Spanos, Stochastic finite elements: a spectral approach, Springer-Verlag, New York, 1991. 8 S. Hosder, R.W. Walters, R. Perez, A non-intrusive polynomial chaos method for uncertainty propagation in CFD simulations, AIAA-2006-891, 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2006. 9 C. Kennedy, M. Carpenter, Additive Runge-Kutta schemes for convection-diffusion-reaction equations, Appl. Numer. Math. 44 (2003) 139-181. 10 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. 11 B.H.K. Lee, L.Y. Jiang, Y.S. Wong, Flutter of an airfoil with a cubic nonlinear restoring force, AIAA-1998-1725, 39th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Long Beach, California, 1998. 12 G.J.A. Loeven, H. Bijl, Probabilistic collocation used in a two-step approach for efficient uncertainty quantification in computational fluid dynamics, CMES 36 (2008) 193–212. 13 L. Mathelin, M.Y. Hussaini, Th.A. Zang, Stochastic approaches to uncertainty quantification in CFD simulations, Num. Alg. 38 (2005) 209–236. 14 R.E. Melchers, Structural reliability: analysis and prediction, Wiley, New York, 1987. 15 Openfem - A finite element toolbox for Matlab and Scilab. Available on: http://www-rocq.inria.fr/OpenFEM/, release 2006a, 2006. 16 C.L. Pettit, P.S. Beran, Spectral and multiresolution Wiener expansions of oscillatory stochastic processes, J. Sound Vib. 294 (2006) 752–779. 17 M.T. Reagan, H.N. Najm, R.G. Ghanem, O.M. Knio, Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection, Combust. Flame 132 (2003) 545–555. 18 S. Sarkar, J.A.S. Witteveen, G.J.A. Loeven, H. Bijl, Effect of uncertainty on the bifurcation behavior of pitching airfoil stall flutter, J. Fluid Struct. 25 (2009) 304–320. 19 J.A.S. Witteveen, S. Sarkar, H. Bijl, Modeling physical uncertainties in dynamic stall induced fluid-structure interaction of turbine blades using arbitrary polynomial chaos, Comput. Struct. 85 (2007) 866–878. 20 J.A.S. Witteveen, H. Bijl, A monomial chaos approach for efficient uncertainty quantification in nonlinear problems, SIAM J. Sci. Comput. 30 (2008) 1296–1317. 21 J.A.S. Witteveen, H. Bijl, Efficient quantification of the effect of uncertainties in advection-diffusion problems using polynomial chaos, Numer. Heat Tr. B-Fund. 53 (2008) 437–465. 22 J.A.S. Witteveen, G.J.A. Loeven, S. Sarkar, H. Bijl, Probabilistic collocation for period-1 limit cycle oscillations, J. Sound Vib. 311 (2008) 421–439. 23 J.A.S. Witteveen, H. Bijl, An unsteady adaptive stochastic finite elements formulation for rigid-body fluid-structure interaction, Comput. Struct. 86 (2008) 2123–2140.

15 of 16 American Institute of Aeronautics and Astronautics

24 J.A.S. Witteveen, H. Bijl, An alternative unsteady adaptive stochastic finite elements formulation based on interpolation at constant phase, Comput. Method Appl. M. 198 (2008) 578–591. 25 J.A.S. Witteveen, G.J.A. Loeven, H. Bijl, An adaptive stochastic finite elements approach based on Newton-Cotes quadrature in simplex elements, Comput. Fluids 38 (2009) 1270–1288. 26 J.A.S. Witteveen, H. Bijl, A TVD uncertainty quantification method with bounded error applied to transonic airfoil flutter, Commun. Comput. Phys. 6 (2009) 406–432. 27 J.A.S. Witteveen, H. Bijl, Higher period stochastic bifurcation of nonlinear airfoil fluid-structure interaction, Math. Probl. Eng. (2009) in press. 28 J.A.S. Witteveen, H. Bijl, Effect of randomness on multi-frequency aeroelastic responses resolved by unsteady adaptive stochastic finite elements, J. Comput. Phys. (2009) submitted. 29 D. Xiu, G.E. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2002) 619–644. 30 E. Yates Jr., AGARD standard aeroelastic configurations for dynamic response. Candidate configuration I.-Wing 445.6, Technical Memorandum 100492, NASA, 1987. 31 A.H. van Zuijlen, H. Bijl, Implicit and explicit higher-order time integration schemes for structural dynamics and fluidstructure interaction computations, Comput. Struct. 83(2-3) (2005) 93-105. 32 A.H. van Zuijlen, H. Bijl, Implicit and explicit higher-order time integration schemes for fluid-structure interaction computations, Int. J. Multiscale Comput. Eng. 4(2) (2006) 255-263. 33 A.H. van Zuijlen, A. de Boer, H. Bijl, Higher-order time integration through smooth mesh deformation for 3D fluidstructure interaction simulations, J. Comput. Phys. 224 (2007) 414-430.

16 of 16 American Institute of Aeronautics and Astronautics