Delft Aerospace Computational Science ENGINEERING MECHANICS
REPORT
May 2006Numerical investigation of the membrane
fluid-structure-interaction problem
R. in ’t Groen and E.H. van Brummelen rob@intgroen.net
All rights reserved. This work may not be translated or copied in whole or in part without the written permission of the author.
TUD/LR/EM
Kluyverweg 1, 2629HS Delft, The Netherlands
P.O. Box 5058, 2600 GB Delft, The Netherlands
Phone: +31 15 278 5460
Fax: +31 15 261 1465
Numerical investigation of the membrane
fluid-structure-interaction problem
Rob in ’t Groen
Delft University of Technology, Faculty of Aerospace Engineering P.O. Box 5058, 2600 GB Delft, The Netherlands
for the smile that inspired me to finalize . . .
ABSTRACT
3
Table of Contents
1 Introduction 5 2 Model problem 7 2.1 Problem statement . . . 7 2.1.1 Fluid problem . . . 7 2.1.2 Structure problem . . . 9 2.1.3 Interface conditions . . . 9 2.2 Non-dimensionalized equations . . . 92.2.1 Non-dimensional fluid problem . . . 10
2.2.2 Non-dimensional structure problem . . . 11
2.2.3 Non-dimensional interface conditions . . . 11
2.2.4 Similarity parameters . . . 12
3 Discretization 13 3.1 Fluid discretization . . . 13
3.2 Structure discretization . . . 15
3.2.1 First-order-system formulation . . . 15
3.2.2 A space-time finite element domain . . . 16
3.2.3 Weak formulation . . . 16
3.2.4 Riemann solution on inter element faces . . . 17
3.2.5 Riemann solution on domain boundaries . . . 19
3.2.6 Displacement . . . 19
3.2.7 Aggregated structural model . . . 20
3.2.8 Non-linear model extension . . . 21
4 Numerical results 25 4.1 Structural model convergence . . . 25
4.1.1 Norm definitions . . . 26
4.1.2 Model error . . . 27
4.2 Computational template and reference . . . 27
4.2.2 Reference case . . . 30 4.3 Assessment of string behavior . . . 31
5 Conclusions and recommendations 35
I Subiteration method 37
5
Chapter 1
Introduction
Numerical solution methods for fluid-structure-interaction problems are of great importance in a multitude of physical and engineering disciplines, e.g., aerospace engineering, offshore engineering, bio-mechanics and civil engineering. When using numerical methods to solve non linear systems of equations many different types of instability can occur. The underlying mathematical problem can have an unstable solution, i.e. a solution that grows unbounded in time, or the numerical solution procedure can be unstable. These instabilities may cause the performing computer program to crash, or even worse, result in misleading solutions.
In the present work a membrane interacting along its surface with an aerodynamic flow is investigated. This particular problem is of interest for research focused on the behavior of flexible aerospace struc-tures. To get a full understanding of the fluid-structure-interaction behavior we start with a simplified problem, such that solutions are not dominated by complex physical effects. This should give a clear view on the stability behavior, both physically and numerically, and allow us to determine the main influence factors.
The next chapter elaborates the mathematical formulation of the fluid-structure-interaction model considered in this thesis. It sheds light on the 1D vibrating membrane equation, i.e. the string equation, the Euler flow equations driving the string and the conditions at the interface between both models. Secondly, a non-dimensionalization will be presented showing the similarity parameters for the membrane behavior.
Chapter 4 presents the numerical results, giving insight in the convergence and stability behavior of the vibrating string in a fluid-structure interaction. For this purpose, 2 different norm types are employed, viz. the L2-norm and the energy norm. We present convergence plots of the structural
model in both norms. In addition, we present the typical solutions of the vibrating string driven by an Euler flow.
7
Chapter 2
Model problem
2.1. Problem statement
The fluid-structure-interaction problem we consider in this thesis comprises the 1D vibrating mem-brane equation, i.e. the string equation, subject to a pressure force induced by an Euler flow. 2.1.1 Fluid problem
The Euler equations describe inviscid compressible fluid flows. In our problem we consider an Euler flow on an open bounded space-time domain with spatial coordinates x, y and temporal coordinate t: Qf = {(x, y, t) |0 < t < T ; −∞ < x < ∞; 0 < y < ∞}, where T defines the end of the period of
interest. Let us remark that in the fluid-structure-interaction problem, the lower boundary at y = 0 in fact sustains a deformation. However, in this work we apply a so-called geometric linearization. This implies that the fluid domain remains undeformed.
PSfrag replacements x y z α L
The boundary of this domain comprises the fixed boundary ∂Qf and the interface Γi = {(x, t) |0 <
t < T ; 0 < x < L}, with L the length of the string.
The unsteady Euler equations supplemented with the continuity and energy equation in conservative form are given by:
∂qf ∂t + ∂f (qf) ∂x + ∂g(qf) ∂y = 0, (2.1a) with: qf = ρ ρu ρv E , f(qf) = ρu ρu2+ p ρuv (E + p) u , g(qf) = ρv ρuv ρv2+ p (E + p) v . (2.1b)
In equation (2.1), qf(x, y, t) denotes the fluid state vector in which ρ, u, v and E denote the fluid
density, horizontal velocity, vertical velocity and total energy. Additionally, p denotes the fluid pressure which is related to the fluid state variables through the equation of state:
p = (γ − 1) E −1 2ρ u
2
+ v2 , (2.2)
with γ the specific heats ratio. In our problem γ = 1.4, which is the typical γ-value for a flow of air. In case of a supersonic flow, due to the characteristic directions, 4 boundary conditions need to be prescribed at the inflow boundary and none at the outflow boundary. More about boundary conditions for the 2D Euler flow can be read in ref. [7]. In case of a supersonic flow problem we prescribe:
qf(−∞, y, t) = ρ∞ ρ∞u∞ 0 E∞
T
. (2.3)
In contrast, when a subsonic flow is considered, we need to enforce 3 boundary conditions at the inflow and 1 at the outflow boundary. In that case we provide:
qf(−∞, y, t) = ρ∞ ρ∞u∞ 0 E T
, (2.4a)
qf(∞, y, t) = ρ ρu ρv E∞ T
. (2.4b)
Let us remark that by taking v∞ = 0 at y = 0 in equation (2.3) and (2.4a) the inflow is specified parallel to the stationary wall. The fluid problem within our fluid-structure-interaction problem can be described as a fluid subject to a moving boundary. We label the boundary offset as α : Γi → R.
Along Γi and along the stationary boundary at y = 0 we define tangent wall conditions i.e. the fluid
flow is tangent to the wall. In addition, we define the total energy on the wall boundary equal to the total energy from the inflow conditions i.e. Ewall= E∞.
To complete the description of the fluid problem the equations are supplemented with the initial conditions:
qf(x, y, 0) = qf 0(x, y), (2.5)
2.2. Non-dimensionalized equations 9
2.1.2 Structure problem
The equation describing the vibrating string is based on a simple equilibrium of forces in the transversal direction. It is assumed that the the string does not move in the longitudinal direction. The structure is defined on a space-time domain Qs= {(x, t)|0 < t < T ; 0 < x < L}. We solve for the 1D vibrating
membrane equation i.e. a vibrating string, which is driven by the pressure difference over the interface. One side of the string gets its pressure from the Euler flow and the other side is subject to a reference pressure p0. The 1D vibrating membrane equation is then given by:
m∂
2z
∂t2 − σ
∂2z
∂x2 = π, (2.6)
with z the string deflection, m the structure density, σ the in plane tension and π = p0−p the pressure
difference (net pressure force) that drives the string. At the boundary we keep the string fixed to a zero deflection:
z(0, t) = z(L, t) = 0. (2.7)
To complete the description we supplement the problem with appropriate initial conditions:
z(x, 0) = z0(x), (2.8a)
˙z(x, 0) = ˙z0(x). (2.8b)
2.1.3 Interface conditions
On the interface we have a set of conditions which determine the coupling of the fluid equations and the structure equation. These conditions can be split into kinematics and dynamics. The kinematic conditions are given by:
∂α ∂xu +
∂α
∂t − v = 0, (2.9a)
α − z = 0. (2.9b)
The first condition (2.9a) states that the vertical component of the fluid velocity has to match the interface displacement velocity, see figure 2.2. We remark that when the fluid domain is not geomet-rically linearized, the second condition (2.9b) implies that the displacement of the membrane equals the displacement of the fluid boundary. The dynamic condition is given by:
π − (p0− p) = 0. (2.10)
This condition simply states that the force π on the structure is due to the pressure difference between the fluid and the reference pressure p0. Throughout, we select p0equal to the pressure of the in flowing
fluid, such that an equilibrium is achieved with a uniform flow and a zero string deflection.
2.2. Non-dimensionalized equations
PSfrag replacements v u∂α ∂x + u ∂α ∂t ∆α ∆x
Figure 2.2: Illustration of the interface condition.
2.2.1 Non-dimensional fluid problem
We introduce the scaling parameters: ˆx= x/L, ˆτf = tc∞/L, ˆρ = ρ/ρ∞ and ˆE = E/E∞, where c∞
denotes the speed of sound. Additionally, we introduce the function ˆqf(ˆx, ˆτf) = L−1qf Lˆx, Lc−1∞ˆτf.
We use this scaling to transform equation (2.1a) into its non-dimensional form: ρ0c∞ L ∂ˆq f ∂ ˆτf +∂f (ˆqf) ∂ ˆx + ∂g(ˆqf) ∂ ˆy = 0, (2.11) and thus: ∂ˆqf ∂ ˆτf +∂f (ˆqf) ∂ ˆx + ∂g(ˆqf) ∂ ˆy = 0. (2.12)
One may note that the non-dimensional fluid equations have the same form as the original dimensioned equations.
This set of non-dimensional equations is supplemented with the appropriate non-dimensional initial conditions:
ˆ
qf(ˆx, 0) = 1 u∞/c∞ 0 1 T
, (2.13)
and boundary conditions:
ˆ
qf(ˆxbdy, ˆτf) = ˆqbdy(ˆx, ˆτf) . (2.14)
Moreover we specify for the solid wall boundaries: ˆ
qf(ˆxwall, ˆτf) = ρˆ ˆu 0 1
T
. (2.15)
2.2. Non-dimensionalized equations 11
for the interface boundary condition, and the in-/outflow conditions. The latter depend on whether a subsonic or supersonic flow case is considered.
2.2.2 Non-dimensional structure problem
For the structural non-dimensionalization we scale with the parameters: ˆx = x/L and ˆτs = tcs/L.
Herein, the structural characteristic velocity cs is chosen as the velocity corresponding with the
fre-quency of the first eigenmode of the string:
cs=
r σ
m. (2.16)
In addition, we introduce the function ˆz (ˆx, ˆτs) = L−1z Lˆx, Lc−1s ˆτs. This scaling transforms equation
(2.6) into its non-dimensional form: ∂2zˆ ∂ ˆτ2 s − ∂2zˆ ∂ ˆx2 = ρ∞L m c2 ∞ c2 s ˆ π. (2.17)
It is to be noted that in equation (2.17) there are no parameters in the left-hand side. Further, we notice that the source term contains ratios of fluid parameters over structure parameters. The right-hand side presents a mass ratio ρ∞L/m and a characteristic velocities ratio cs/c∞.
To complete the non-dimensional structure problem, we provide the non-dimensional initial and boundary conditions:
ˆ
z (ˆx, 0) = ˆz0(ˆx) , ∂ ˆz (ˆx, 0) /∂ ˆτs= ∂ ˆz0(ˆx) /∂ ˆτs, z (0, ˆˆ τs) = ˆz (1, ˆτs) = 0. (2.18)
2.2.3 Non-dimensional interface conditions
Using the previously introduced scaling parameters, the kinematic condition from equation (2.9b) is cast into:
ˆ
α (ˆx, ˆτf) − ˆz(ˆx, ˆτs) = ˆα (ˆx, ˆτf) − ˆz ˆx, csc−1∞τˆf = 0. (2.19)
We remark that the non-dimensionalization of this kinematic condition is only visible in the scaling of the arguments. Using the same scaling parameters, equation (2.9a) is cast into:
∂ ˆα ∂ ˆxu +ˆ
∂ ˆα
∂ ˆτf − ˆv = 0,
(2.20)
which is in the same form as the original dimensioned equation. However due to the argument scaling from equation (2.19), it holds that:
∂ ˆα ∂ ˆx = ∂ ˆz ∂ ˆx and ∂ ˆα ∂ ˆτf = cs c∞ ∂ ˆz ∂ ˆτs . (2.21)
Thus, we find again the ratio of characteristic velocities cs/c∞. Non-dimensionalization of our dynamic
condition (2.10) results in:
ˆ
2.2.4 Similarity parameters
Within the non-dimensionalized equations 3 similarity parameters characterizing the fluid-structure-interaction system appear:
ν = u∞/c∞ (Mach number),
κ = cs/c∞ (characteristic velocities ratio),
µ = ρ∞L/m (mass ratio),
(2.23)
13
Chapter 3
Discretization
We use a spectral finite-element method to discretize the fluid and structure equations. Anticipating the extension to deforming fluid domains, we opt for a space/time discretization, as such a dis-cretization allows for a straightforward implementation of moving domains, see ref. [5]. The fluid discretization is a standard stabilized continuous Galerkin discretization with weakly enforced initial conditions. We provide a concise specification in section 3.1. Further details can be found in, e.g., ref. [1]. The discontinuous Galerkin discretization for the structure and the corresponding treatment of the inter-element discontinuities are elaborated in section 3.2.
3.1. Fluid discretization
Since in this thesis the focus is on the discretization and behavior of the string we will not fully elaborate the discretization of the fluid equations. However, for completeness we provide a brief
PSfrag replacements x x y −0.25 1.25 0.00 1.00
specification of the fluid discretization.
We define a space-time domain extended with elements in front and aft, w.r.t. the flow direction, of the perturbed interface boundary. Figure 3.1 gives an illustration of the space domain. It shows the fluid-structure-interaction domains perpendicular to the time axis and visualizes clearly the element extension. These additional elements avoid conflicting definitions between boundary conditions at the moving interface and boundary conditions at the inflow.
To discretize the fluid problem with a Galerkin finite-element method, we cast equation (2.12) into a space-time variational statement as described in ref. [3]:
find qf ∈ U such that:
Z Q φ· ∂q f ∂t + ∂f (qf) ∂x + ∂g(qf) ∂y dQ = 0 ∀φ ∈ U, (3.1)
in which the product space U = U1× U2× U3× U4 accommodates the trial functions for qf and the
test functions φ. Integrating by parts and adding a stabilizing term yields:
Z ∂Q φ· qfnt+ f (qf) nx+ g(qf) ny dS − Z Q ∂φ ∂t · qf+ ∂φ ∂x · f (qf) + ∂φ ∂y · g(qf) dQ +X e Z Qe (Lφ) · τ rf qˆf dQe= 0. (3.2)
In equation (3.2), nt, nx and ny denote the t, x, y components of the space/time outward normal
vector along the boundary ∂Q. Furthermore the last integral in this equation is the least-squares operator as described in ref. [8]. Herein, is L the differential operator for the Euler equations:
L = ∂t∂ +∂f qf ∂qf ∂ ∂x+ ∂g qf ∂qf ∂ ∂y, (3.3) is rf qˆf
the residual of the Euler equations and τ a matrix that enhances the stability of the numerical solution.
We consider a Galerkin method that is continuous in space and discontinuous in time. Accordingly, the initial conditions are enforced weakly, and the boundary conditions are enforced strongly. Weak enforcement of the initial conditions is accomplished by integrating equation (3.2) by parts and re-placing the terms at t = 0 by their specification according to the initial conditions. To enforce the boundary conditions strongly, the trial functions are set in accordance with the boundary conditions and the corresponding test functions are constrained to zero. Furthermore, since we cannot compute on an infinite domain, we restrict the extend of the domain in y-direction. This introduces a boundary at y = ymax for which we provide a solid wall boundary condition.
We expand the test and trial functions with a finite number of shape functions: φ = Pˆ
φiN˜i and
qf =P ˆqf jN˜j, where ˜Ni = Ni1(x)Ni2(y)Ni3(t) combined with a renumbering scheme from i1, i2, i3
to i. In these expansions ˆφi, ˆqf j denote the weights of the shapes and Np(ξ) denotes a shape function
3.2. Structure discretization 15 X i X j ˆ φi Z ∂Q n ˜NiN˜jˆq f j+ ˜Nif( ˜Njqˆf j) + ˜Nig( ˜Njˆqf j) o dS −X i X j ˆ φi Z Q ( ∂ ˜Ni ∂t N˜jˆqf j+ ∂ ˜Ni ∂x f( ˜Njˆqf j) + ∂ ˜Ni ∂y g( ˜Njqˆf j) ) dQ = 0 (3.4)
We consider equation (3.4) in the form of a residual equation: rf ˆqf = 0. In order to solve this
residual equation with a Newton-Raphson method we require the Jacobian of our residual: Jf(ˆqf) = X i X j Z ∂Q ˜ NiN˜j+ ˜Ni ∂f (qf) ∂qf N˜j+ ˜Ni ∂g(qf) ∂qf N˜j dS −X i X j Z Q ( ∂ ˜Ni ∂t N˜j+ ∂ ˜Ni ∂x ∂f (qf) ∂qf ˜ Nj+ ∂ ˜Ni ∂y ∂g(qf) ∂qf ˜ Nj ) dQ (3.5)
The Newton-Raphson iterative update procedure is given by: ˆ
qn+1f = ˆqnf+ δ, (3.6)
in which the update δ is determined from:
Jf qˆf δ = −rf ˆqf . (3.7)
Herein, is Jf the Jacobian of the residual as in equation (3.5) and rf is the residual itself i.e. the
left-hand side of equation (3.2).
3.2. Structure discretization
This section elaborates the (non-standard) discretization of the structural equation. A discontinuous Galerkin (DG) method is very suitable for hyperbolic problems, see ref. [2], [6]. Since the 1D vibrating string equation in space-time is genuinely hyperbolic the choice for a DG method is comprehensible. To introduce Riemann solutions on the inter-element boundaries we rewrite the string equation in a system of first-order wave equations in section 3.2.1. In the next section we introduce the space-time domain for our structural equations. Then, in section 3.2.3, we present a weak formulation for the first-order system. In sections 3.2.4 and 3.2.5 we provide the Riemann solutions on the inter-element edges. In section 3.2.6 we determine the displacement from the wave system solution. In section 3.2.7 we condense the discretized problems into one spectral DG method for the structure. Finally, we provide some expressions to incorporate a nonlinear elastic effect in the vibrating string.
3.2.1 First-order-system formulation
To formulate the DG method, one requires additional information on the boundaries relating the discontinuous elements. We notice that equation (2.17) is a second order wave equation with a source. Rewriting this equation as a system of first order wave equations provides a very natural setting to introduce Riemann solutions for the inter element fluxes and to provide the required boundary relations. We repeat equation (2.17), dropping hats for convenience:
PSfrag replacements x t h ξ τ ζ
Figure 3.2: Illustration of the discretized space-time domain coordinates.
Defining us = ∂z/∂t and vs = ∂z/∂x equation (3.8) translates into a system of first-order wave
equations: ∂qs ∂t + Aij ∂qs ∂x = b, (3.9a) with qs= us vs , Aij= 0 −1 −1 0 , and b = (µ/κ2)π 0 . (3.9b)
3.2.2 A space-time finite element domain
The structure space-time domain Qs as defined in section 2.1.2 is discretized using finite elements.
Since there is no need for difficult element shapes we choose these space-time elements to be rectangular as depicted in figure 3.2. Each element has lengths h and τ in space and time, respectively. We define local coordinates ξ and ζ on each element by the following coordinate transformations:
ξ = x − xc
h/2 , (3.10a)
ζ = t − tc
τ /2 , (3.10b)
where xc and tc represent the center of the element in global space-time coordinates. In section 3.2.7
the local coordinates ξ, ζ are used to define the basis of shape functions. 3.2.3 Weak formulation
To discretize the system of wave equations by means of the discontinuous Galerkin method, we define the following variational problem:
3.2. Structure discretization 17
The product space Q = Q1× Q2 accommodates the test functions ψ and the trial functions for qs.
For the integration over the space-time domain Q we sum the integrations over the element domains Qe. In addition, to introduce a weak enforcement of continuity and boundary conditions we apply the
divergence theorem in space/time. This results in the wave system weak formulation:
X e Z Ωe (ψ · qs) ntdΩ + X γ Z Γi ψ(x− γ, t) − ψ(x+γ, t) · AijqF Ws nxdΓ +X b Z Γb (ψ · Aijqs) nxdΓ − X e Z Qe ∂ψ ∂t · qs+ ∂ψ ∂x · Aijqs dQ =X e Z Qe ψ· bdQ (3.12)
Herein, the sums over e, γ and b are, respectively, the sums over the elements, the inter-element boundaries and the element domain boundaries. Additionally, nt and nx denote, respectively, the
time and space component of the outward normal vector of the element boundary, and ψ(x− γ, t),
ψ(x+
γ, t) denote the test functions respectively on the left and right side of the inter-element boundary
Γi. To enable a weak enforcement of the initial conditions qs0we replace the argument of the integral
on the lower time boundary t = 0 by the initial conditions. Furthermore, we remark that to employ the weak formulation (3.12) we require qF W
s on the inter-element boundaries. To this end, we introduce
in section 3.2.4 a Riemann solution at these inter-element discontinuities.
We summarize the discretization of the first-order system in a variational format:
find qs∈ Q such that: a(ψ, qs) = b(ψ) ∀ψ ∈ Q, (3.13)
with: a(ψ, qs) =X e Z h 0 {ψ(x, τ) · qs(x, τ )} dx +X γ Z τ 0 ψ(x− γ, t) − ψ(x+γ, t) · AijqF Ws (xγ, t) dt +X b Z τ 0 {ψ(L, t) · A ijqs(L, t) − ψ(0, t) · Aijqs(0, t)} dt −X e Z τ 0 Z h 0 ∂ψ ∂t · qs+ ∂ψ ∂x · Aijqs dxdt (3.14) and b(ψ) =X e Z τ 0 Z h 0 ψ· b dxdt +X e Z h 0 {ψ(x, 0) · qs0(x)} dx. (3.15)
3.2.4 Riemann solution on inter element faces
Our discretization allows discontinuities on the element boundaries. This implies that our solution is not uniquely defined on those boundaries. Therefore, we are interested in expressions which determine qs at those locations.
18 Chapter 3. Discretization x x x t t t q− s q− s q+s q+ s qF Ws λ1 λ2 n− n+
Figure 3.3: Illustration of the Riemann problem on inter element faces.
running wave emanate with velocities λ1 and λ2 respectively, see figure 3.3. The Riemann solutions
for hyperbolic problems are elaborated in more detail in ref. [7]. For such a Riemann problem it is known that the state vector changes across waves in space-time and that these states can be related by means of Riemann invariants. The Riemann invariants will thus provide us the required relations on our inter-element boundaries. The Riemann invariants ψk can be found from:
(∇ψk, Rk) = 0. (3.16)
Supplemented with the right eigenvectors Rk expression (3.16) provides two equations for the two
Riemann invariants.
Matrix Aijof the first-order-system (3.9) has two distinct eigenvalues with corresponding eigenvectors:
h λ1= −1, R1= 1 1 Ti ,hλ2= 1, R2= −1 1 Ti . (3.17)
Substituting these eigenvectors into equation (3.16) results in: ∂ψ1 ∂u + ∂ψ1 ∂v =0 ⇒ ψ1= u + v, (3.18a) −∂ψ∂u2+∂ψ2 ∂v =0 ⇒ ψ2=−u + v. (3.18b)
Using the property that for a linear system a Riemann invariant remains constant across its corre-sponding wave we can relate the inter element boundary variables uF W
s , vsF W with the left and right
state variables: uF W s + vsF W = u−s + v−s, (3.19a) −uF W s + vsF W =−u + s + v + s. (3.19b)
With these relations and introducing the element normal vectors as in figure 3.3 we find for the state vector on the inter element faces:
qF Ws =1 2 1 n− x n− x 1 q− s +12 1 n+ x n+ x 1 q+s, (3.20) where n−
x and n+x denote the components in x-direction of the normal vectors n− and n+. Equation
(3.20) gives us an expression for our wave system state qs on a inter-element discontinuity only
3.2. Structure discretization 19
3.2.5 Riemann solution on domain boundaries
The discontinuities on the domain boundaries can be treated similar to those on the inter element faces. However, since we only have one incoming characteristic on these locations versus two unknowns, we need to provide the system with a boundary condition. In section 2.1.2 we introduced the boundary condition z(0, t) = z(L, t) = 0. Since this condition does not define the fluxes we substitute this condition, with the following boundary condition:
∂z(0, t)
∂t =
∂z(L, t)
∂t = 0. (3.21)
Provided with the initial condition z0(0) = z0(L) = 0 this yields an equivalent enforcement of the
boundary conditions.
We first consider the right domain boundary. For this boundary we deduct one relation from the Riemann invariants and we prescribe us= 0:
uRWs + vsRW = u−s + v−s, (3.22a)
uRWs = 0. (3.22b)
Using the introduced normal vectors we can now write an expression for qs at the right boundary
depending only on left hand side information: qRWs = 0 0 n− x 1 q− s. (3.23)
Using a similar procedure for the left domain boundary we find: qLWs = 0 0 n+ x 1 q+s. (3.24)
Equations (3.20), (3.23) and (3.24) provide us expressions needed to weakly incorporate the boundary conditions into our wave system weak formulation (3.12).
3.2.6 Displacement
In formulating our first order system we defined us= ∂z/∂t, which is in itself a differential equation
for the displacement. We can solve this equation on the same finite elements as we used for solving the first order wave system. We define the variational problem as:
find z ∈ Z such that: Z Q φ∂z ∂tdQ = Z Q φusdQ ∀φ ∈ Z (3.25)
Where Z accommodates the test function φ and the trial function z. Using integration by parts in time and summing the element integrations for the entire space/time domain integration, we get a weak formulation very similar to the weak formulation of the first-order system:
X e Z Ωe φ z ntdS − X e Z Qe ∂φ ∂t z dQ = X e Z Qe φ usdQ (3.26)
exerted indirectly in a weak sense through the first order system. We summarize the displacement discretization similar to the first-order summary in variational format:
find z ∈ Z such that: ¯a(φ, z) = ¯b(φ) ∀φ ∈ Z, (3.27)
with: ¯ a(φ, z) =X e Z h 0 φ(x, τ ) z(x, τ ) dx −X e Z τ 0 Z h 0 ∂φ ∂t z dxdt (3.28) and ¯b(φ) = X e Z τ 0 Z h 0 φ usdxdt + X e Z h 0 φ(x, 0) z0(x) dx. (3.29)
3.2.7 Aggregated structural model
We now combine the wave system equations and the displacement equation. Accumulating the equa-tions in one system allows us to solve the equaequa-tions for the structure at once. Solving the aggregated system is slightly more expensive. However, within our fluid-structure-interaction the fluid computa-tion is much more expensive, and this higher expense for the structure is negligible. For this purpose we define: q∗= q s z T and ψ∗= ψ φ T , (3.30)
such that we can condense the two variational formulations (3.13) and (3.27) of respectively the first-order system and the displacement into one:
find q∗∈ Q × Z such that: a(ψ, q
s) + ¯a(φ, z) = b(ψ) + ¯b(φ) ∀ ψ∗∈ Q × Z. (3.31)
To approximate the test and trial functions, we use a basis of shape functions derived from the Jacobi polynomials. The Jacobi polynomials, denoted by Pα,β
p (x), possess the important property
of orthogonality. The Jacobi polynomials can be constructed using a recursion relationship. More details about these polynomials and the modal expansion basis we use can be found in ref. [4]. The recursion relationship for Jacobi polynomials is given by:
3.2. Structure discretization 21
We remark that when α = β = 0 the Jacobi polynomials reduce to Legendre polynomials and when α = β = −1 to Chebychev polynomials.
We expand the test and trial functions with a finite number of shape functions: ψ = Pˆ ψiN˜i,
qs=P ˆqs jN˜j, φ =PφˆkN˜kand z =P ˆzlN˜l, where ˜Ni= Ni1(x)Ni2(t) combined with a renumbering
scheme from i1, i2 to i. In these expansions ˆψi, ˆqs j, ˆφk, ˆzldenote the weights of the shapes and Np(ξ)
denotes a shape function of order p from our modal basis. This shape is constructed on a reference element Qst= {(ξ)| − 1 < ξ < 1} as: Np(ξ) = 1−ξ 2 p = 0 1−ξ 2 1+ξ 2 Pp−1α,β(ξ) 0 < p < P 1+ξ 2 p = P (3.34)
With the introduced shape functions we rewrite the weak formulation of the aggregated structure system into: X i X j a ˆψiN˜i, ˆqs jN˜j +X k X l ¯ a ˆφkN˜k, ˆzlN˜l =X i b ˆψiN˜i +X k ¯b ˆφkN˜k ∀ ˆψi, ˆφk (3.35)
Equation (3.35) gives us a representation of our problem, which we can cast in the form Ax = b, for wich many standard solution procedures are available.
3.2.8 Non-linear model extension
For a more realistic representation of the membrane behavior, the system of equations can be enhanced with a nonlinear elasticity. In this section we propose an energy-conservative non-linearity. For this purpose we introduce: ∂2z ∂t2 − k (zx) ∂2z ∂x2 = 0, with k (zx) = k0+ k1 ∂z ∂x 2 . (3.36)
To cast equation (3.36) into the format that fits the discretization introduced previously, we present the non-linear wave system and the Riemann solution on the discontinuous boundaries.
First, we ensure that the nonlinearity k (zx) leads to an energy conservative formulation. To this end
we multiply equation (3.36) by ∂z/∂t and integrate over the entire space-time domain: Z L 0 Z T 0 ∂z ∂t ∂2z ∂t2 − k0 ∂2z ∂x2 − k1 ∂z ∂x 2∂2z ∂x2 ! dtdx = 0 (3.37)
Since we have defined stationary boundary conditions for the entire space-time domain the second part vanishes, which leaves us with an expression telling that our non-linear structural model conserves the quantity: Es= Z L 0 ( 1 2 ∂z ∂t 2 +k0 2 ∂z ∂x 2 +k1 12 ∂z ∂x 4) . (3.39)
The quantity Esaccording to equation (3.39) can be identified as the energy of the non-linear structure.
We cast the non-linear structure equation into a system of first order equations using the definitions introduced in this work for the linear structure case: us= ∂z/∂t and vs= ∂z/∂x. We rewrite equation
(3.36) in: us vs t + 0 −k0− k1vs2 −1 0 us vs x = 0. (3.40)
For the matrix in equation (3.40), the eigenvalues and corresponding eigenvectors are given by: h λ1= −pk0+ k1vs2, R1= pk0+ k1v2s 1 Ti , h λ2=pk0+ k1v2s, R2= −pk0+ k1v2s 1 Ti . (3.41)
Using the same approach to find relations for the discontinuous boundaries as in the linear case we recall definition (3.16) from which we deduct again 2 differential equations to find the Riemann invariants: ∂ψ1 ∂us pk0+ k1v2s+ ∂ψ1 ∂vs = 0, (3.42a) −∂ψ2 ∂us pk0+ k1vs2+ ∂ψ2 ∂vs = 0. (3.42b)
Using the ansatz:
ψ1,2= ±us+ f (vs), (3.43)
we find for the Riemann invariants:
ψ1,2= ±us+ Z vs v∗ p k0+ k1s2ds (3.44a) = ±us+ vs 2pk0+ k1v 2 s+ k0 2√k1 ln pk1vs+pk0+ k1v 2 s + C. (3.44b)
For the integration constant we choose C = − k0
2√k1 ln √ k0
such that if we take k0= 1 and the limit
for k1to zero we retrieve the Riemann invariants of the linear system:
lim
3.2. Structure discretization 23
We remark that the Riemann approach used in the linear case is not exact for the nonlinear case and can only be used for a limited type of waves. Albeit, assuming that a Riemann invariant is almost constant across its corresponding wave we find two equations that determine the structural state variables along our discrete boundary:
uw=12(u++ u−) + Z v+ v− p k0+ k1s2ds, (3.46a) (u+− u−) + Z v+ vw p k0+ k1s2ds + Z v− vw p k0+ k1s2ds = 0. (3.46b)
With equation (3.46a) we find an explicit expression for uw. On the other hand vwis only implicitly
determined by equation (3.46b). To obtain vw, Newton’s method can be used. Newton’s method is
defined by the iterative process:
vn+1
w = vnw+ δn, (3.47)
where n denotes the iteration number and δn the sequential update, which is determined from:
r0(vn
w) δn = −r (vnw) . (3.48)
Here r (vn
w) is the residual of equation (3.46b) and r0(vwn) is its derivative. To initiate the update
routine we provide an initial guess based on the value in the linear case, from equation (3.20): vw0 = 1 2 u−w+ v−w + 1 2 −u + w+ v + w (3.49)
25
Chapter 4
Numerical results
In this chapter we present numerical results obtained with the previously introduced discretization. These results serve two purposes. First, we verify the stability of the introduced discretization for the structural model. Second, we are interested in the typical interacting-string solutions. For the vibrating string problem there are no systematic studies available in the literature. Furthermore, the behavior is characterized by 3 similarity parameters and thus the solution depends on a substantial parameter space. This renders the qualification of the behavior complex. We resolved the qualifica-tion by performing numerous computaqualifica-tions with different parameter settings. For this purpose we used the finite-element package MPF, see appendix II. Here we present the typical findings of these computations.
Before we present the results of the string problem we show that the non-conventional discretization of the structure equations leads to a stable and convergent numerical scheme. In section 4.1 we provide two norm definitions for the string deflection and we show that the numerical error converges in these norms at the appropriate rate. Thereafter, in section 4.2, we set up a template for our fluid-structure-interaction computations. We quantify the discretization parameters, initial conditions and boundary conditions, but we leave the similarity parameters as unspecified blanks to be filled for the different computations. Finally, in section 4.3 we show the typical behavior of the string interacting with an Euler flow.
4.1. Structural model convergence
To investigate the convergence of the structural model we employ two different norm types. Firstly, the L2-norm over the integrated space-time domain is considered. Secondly, we consider an energy
4.1.1 Norm definitions
We denote the error between the exact and approximate string deflection as e(x, t) = ˜z(x, t) − z(x, t), where ˜z(x, t) is the approximate solution of z(x, t). The L2-norm of the error over the structural
domain is then determined by:
kek2 L2 = Z Qs |e(x, t)|2dQ s. (4.1)
To determine the error over the discretized space-time domain we compute: kek2L2 = X e Z Qe |˜z(x, t) − z(x, t)|2dxdt, (4.2)
where Qerepresents an element domain.
For the derivation of the energy norm we multiply equation (2.17) by ∂z/∂t, ignore the source term and integrate over the space-time domain:
Z L 0 Z T 0 ∂z ∂t ∂2z ∂t2 − ∂2z ∂x2 dtdx = 0. (4.3)
Integrating by parts we rewrite this to: " Z L 0 ( 1 2 ∂z ∂t 2 +1 2 ∂z ∂x 2) dx #T 0 − " Z T 0 ∂z ∂t ∂z ∂x dt #L 0 = 0. (4.4)
Since we have defined stationary boundary conditions for the entire space-time domain the second part vanishes, which leaves us with an expression showing that our model conserves the quantity:
Es= Z L 0 ( 1 2 ∂z ∂t 2 +1 2 ∂z ∂x 2) dx (4.5)
In expression (4.5) the first component of the integral exhibits kinetic energy and the second potential energy. From this expression we may extract a sensible norm definition on energy level:
kzk2 E= Z Qs ( ∂z ∂t 2 + ∂z ∂x 2) dQs. (4.6)
The error e in this norm over the structural space-time domain can be computed with:
kek2 E= X e Z Qe ( ∂ ˜z ∂t − ∂z ∂t 2 + ∂ ˜z ∂x− ∂z ∂x 2) dxdt (4.7)
Recalling our wave system variables us and vswe can also write:
kek2
E∗= k˜us− usk
2
4.2. Computational template and reference 27
Here ˜us and ˜vs are the numerical approximations of usand vs. The superscript ∗ indicates that the
energy norm based on primitive variables is considered. 4.1.2 Model error
To investigate the convergence of the structural model we perform error calculations on a space time domain Qs = {(x, t)|0 < t < 1; 0 < x < 1}. On this domain we consider square space-time elements
i.e. h = τ , see figure 3.2. The numerical solution is compared against an exact solution of equation (3.8). For this purpose we define the initial condition:
z0(x) = sin (2πx) . (4.9)
For the structure separately, i.e. in the absence of the source term, this initial condition yields the solution:
z(x, t) = sin (2πx) cos (2πt) . (4.10)
In the hp-convergence study we present the error reduction by h-refinement for different polynomial orders p. Here h = τ denotes the mesh width and p denotes the polynomial expansion order of the shape functions accommodated on the discrete elements. We expect the error in the numerical solution to converge at a rate defined by, see ref. [4]:
k˜ekL2 ≤ chp+1, (4.11a)
k˜ekE ≤ chp. (4.11b)
We remark that a norm based on derivatives, like the considered energy norm, converges with a rate one order lower than the L2-norm. However, if we consider the energy norm as in equation (4.8),
computed from the first order system variables, we can expect the same rate of convergence as in the L2-norm:
k˜ekE∗ ≤ chp+1. (4.12)
In figures 4.1-4.3 the error in respectively the L2-norm and the previously introduced energy norm is
plotted on a log-log scale against the number of elements. From the inclination of the curves as h → 0 we conclude that the convergence rate of the string model shows agreement with the expected rates of convergence stated in equations (4.11) and (4.12).
4.2. Computational template and reference
In this section we provide a template for further computations. We intend to specify all discretization parameters and leave only the similarity parameters variable. Additionally, we provide a reference case to compare the solutions presented later on in this chapter.
4.2.1 Computation parameters
1e-10 1e-08 1e-06 1e-04 0.01 1 1 10 100 k z − ˆzkL 2 h−1 1 1 1 1 2 3 4 5
Figure 4.1: h-type extension error convergence in the L2
-norm for polynomial expanded shape functions up to orders: p = 1, p = 2, p = 3 and p = 4
1e-10 1e-08 1e-06 1e-04 0.01 1 1 10 100 k z − ˆzkE h−1 1 1 1 1 1 2 3 4
4.2. Computational template and reference 29 1e-10 1e-08 1e-06 1e-04 0.01 1 1 10 100 k z − ˆzk E ∗ h−1 1 1 1 1 2 3 4 5
Figure 4.3: h-type extension error convergence in the energy norm, based on primitive variables, for polynomial expanded shape functions up to orders: p = 1, p = 2, p = 3 and p = 4
we fit 96 hexahedron elements per time slab with ∆x = ∆y = h = 0.125 and ∆t = 0.05. On these elements we use second order polynomial expansions: pfluid= 2.
Our structure domain is defined Qs = {(x, t) |0 < t < T ; 0 < x < 1} and comprises 8 equally sized
elements, such that the element vertex points coincide with the vertices of the fluid elements on the interface, i.e. ∆x = h = 0.125 and ∆t = 0.05. On the structure elements we use fourth order shape functions, i.e. pstruct.= 4.
We remark that the polynomial expansion order of the structure and fluid shape functions cannot be chosen arbitrarily. To illustrate this, figures 4.4 and 4.5 display the deflection z for the reference case in section 4.2.2 with pfluid= pstruct.= 3. The figures clearly show non-decaying discontinuities at the
boundaries of the rightmost element. Such discontinuities form an invalid solution for the string. We did not investigate to find the exact cause of this disparity, but we observed that it vanishes when sufficiently high-order shape functions are used for the structural discretization. Conjecturally, this artificial discontinuity results from the weak enforcement of the boundary conditions. During the time evolution, large gradients can develop at the right boundary. Due to the weak enforcement of the boundary conditions and the large gradient us= ∂z/∂t can deviate significantly from zero. Since z
is obtained from usby integration in time, this implies that z will be nonzero at the boundary. Even
if the large gradients dissappear in time and usvanishes at the right boundary, the discrepancy in z
will persist.
ap-−10−4
z +10−4
0.7 x 1.0 0.7 x 1.0 0.7 x 1.0
Figure 4.4: non-physical discontinuity in reference case deflection z for pfluid = pstruct. = 3 at
t= 1.0 (left), t = 2.0 (middle) and t = 5.0 (right)
-4e-05 -2e-05 0 2e-05 4e-05 6e-05 8e-05 0.0001 x t 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5
Figure 4.5: reference case deflection z in space-time for pfluid= pstruct.= 3
proach avoiding many complications ref. [5]. In appendix I, we elucidate the subiteration method applied to our discretized fluid-structure-interaction problem. For our subiteration loop the number of iterations is restricted to 50 loops. Within the Newton iterations we update the left hand matrix when the update tolerance exceeds a limit of 10−7 accuracy.
We complete the computation template with the initial conditions for the string. We use the first free vibration mode of the string, with a small amplitude:
qs0= 10−4 0 π cos (πx) sin (πx) T
(4.13)
Here, the first two components of qs0are the first order system variables usand vs corresponding to
the initial string deflection; the third component of qs0.
4.2.2 Reference case
4.3. Assessment of string behavior 31
comp. Nelem, fluid Nelem, struct. h ∆t pfluid pstruct. Niter update tol.
param. 12 × 8 = 96 8 0.125 0.05 2 4 50 10−7
Table 4.1: Summary of settings within the computational template.
-2e-05 0 2e-05 4e-05 6e-05 8e-05 0.0001 x t 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5
Figure 4.6: reference case: ν = 1.0, µ = 1.0, κ = 1.0
characteristic parameters, see table 4.2
case ν µ κ
ref. 1.0 1.0 1.0
Table 4.2: Similarity parameter settings for the reference case.
The structure solution is presented in a space-time plot, see figure 4.6, where the grayness indicates the deflection of the string. We observe the effects of the flow on the structure motion and we notice that the string oscillation is damped. Secondly we notice that the waves in the structure move in longitudinal direction due to the air flowing past. To contrast this result to the free vibrating case, we note that our structural equation, which is in effect a second order wave equation, yields undamped oscillations when no force is applied.
4.3. Assessment of string behavior
-0.0001 -8e-05 -6e-05 -4e-05 -2e-05 0 2e-05 4e-05 6e-05 8e-05 1e-04 x t 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10
Figure 4.7: stable string solution; damped: ν = 1.2, µ = 0.01, κ = 1.0
instability is called a divergent solution and an oscillatory unstable solution is also referred to as flutter.
We present the three typical cases for the string problem. Figure 4.7 shows a slightly damped solution. In this case the mass ratio µ 1, i.e., the string is heavy with respect to the fluid flow. We notice that a heavier string renders the string less sensitive to variations in the fluid state. In figure 4.8, we present a divergent solution for a subsonic flow case. Conceivably, for a specific range of subsonic Mach numbers this divergence may occur as by Bernouilli’s theorem the pressure force is larger than the retraction force by the longitudinal tension in the string. Lastly, in figure 4.9, we show the occurrence of flutter for the vibrating string. With a sufficient low velocity ratio κ and the other parameters near the reference case, a growing, unbounded oscillatory instability appears.
4.3. Assessment of string behavior 33 -0.0003 -0.0002 -0.0001 0 1e-04 0.0002 0.0003 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 t x
Figure 4.8: unstable, divergent string solution: ν = 0.5, µ = 1.0, κ = 1.0
-0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 x t 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
-0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 x t 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
35
Chapter 5
Conclusions and recommendations
We investigated the fluid-structure-interaction problem of a vibrating string driven by an Euler flow. To assess the vibrating string behavior we non-dimensionalized the system of equations, which resulted in 3 similarity parameters: the Mach number ν, the mass ratio µ and the characteristic velocities ratio κ. We solved the system of equations with a spectral space-time finite-element method. Moreover we introduced a non-conventional discretization technique for the structure equations. We used a discontinuous Galerkin method and to evaluate the solution on the discontinuous boundaries we solved a Riemann problem on those locations. In order to set up the Riemann problem we translated the string equation into a system of first order wave equations. From this system we determined the string fluxes in space and time and we used the time flux to compute the required string deflection. Numerical results for the convergence of the non-conventional string discretization were presented. We showed that the rate of convergence is in agreement with expectations from theory in two considered norms, the L2-norm as well as the energy norm.
A limitation on the choice of the order of shape functions was found. When the numerical string model is formed with insufficiently high-order shape functions, non-physical discontinuities arise, which do not decay in time. We did not investigate the exact cause of this disparity, but we remarked that with a weak enforcement of the boundary conditions it is conceivable that these discontinuities can exist in the numerical solution. Moreover we observed that this disparity vanishes when sufficient high order shape functions are chosen.
At last we investigated the typical behavior of the vibrating string. We presented examples of the 3 specific solutions, viz., damping, divergence and flutter.
37
Appendix I
Subiteration method
The subiteration method is a member of the strongly-coupled partitioned methods. It is sometimes referred to as Picard iteration or successive approximation. The fluid and the structure are solved subsequently on a space-time slab. After a sufficient number of subiterations, when the system is converged and a fully-implicit coupling between the aerodynamics and the structural model is achieved, we proceed to the next time slab. For the considered fluid-structure-interaction model there exist 4 iteration steps i.e. the kinematics, the fluid problem, the dynamics and the structure problem. The connectivity between these different iteration steps is illustrated in table I.1.
Step 1: Kinematics
Within the kinematics step of the subiteration procedure the fluid boundary displacement is deter-mined from the initial structural displacement or the structural displacement received from the former subiteration loop. In order to solve the kinematics step with a Galerkin method the problem for this step can be given in variational format as:
find αj ∈ A such that K(µ, αj, zj−1) = 0 ∀µ ∈ A. (I.1)
Here A accommodates the test and trial functions and K(µ, αj, zj−1) is the kinematics operator. The
subscript j denotes the number of the current iteration loop and the subscript j − 1 refers to the
α qf π qs
K × 0 0 ×
F × × 0 0
D × × × 0
S 0 0 × ×
previous iteration loop or the provided initial approximation if j − 1 = 0. For the fluid-structure-interaction in this report the kinematics operator can be determined by multiplying equation (2.19) with a test function µ and integrating over the space-time domain.
Step 2: Fluid problem
For the fluid problem a fluid flow subject to a particular set of boundary conditions including the conditions along the perturbed interface is solved. The boundary conditions on the interface are defined by the solution of the kinematics step. The fluid problem in variational format is given by:
find qf j∈ U such that F (φ, qf j, αj) = 0 ∀φ ∈ U (I.2)
Here U accommodates the test and trial functions for the fluid problem. Furthermore the operator F (φ, qf j, αj) is given by the weak formulation (3.2), which is elucidated in chapter 3.
At the perturbed interface the fluid is next to the tangent wall condition subject to a vertical velocity due to the oscillating structure. In order to determine this velocity v at the interface we solve for the following variational formulation:
find v ∈ V such that B(ψ, v) = b(ψ, u, αx, αt) ∀ψ ∈ V. (I.3)
Herein V accommodates the test and trial functions and αx and αt denote the derivatives of the
interface in respectively space and time directions. Furthermore for the considered fluid-structure-interaction the weak formulation B(ψ, v) = b(ψ, u, αx, αt) is determined by multiplying equation (2.20)
with test function ψ and integrating over the entire domain. Step 3: Dynamics
In the dynamics step the pressure force acting on the structure is determined. The pressure force depends on the structure location and the fluid state with respect to the reference pressure:
find πj ∈ P such that D(ν, uj, αj, πj) = 0 ∀ν ∈ P (I.4)
Here P accommodates the test and trial functions for the dynamic condition. In addition, the dynamics operator D(ν, uj, αj, πj) = 0 is derived from equation (2.22) by multiplying with test function ν and
integrating over the space-time domain. Step 4: Structure problem
The last step of the subiteration loop solves for the structure subject to a source term. The source term is received from the dynamics step. In variational format this problem is given by:
find zj∈ Z such that S(ψ, zj, πj) = s(ψ) ∀ψ ∈ Z (I.5)
39
Appendix II
MPF programming environment
All solution methods are performed within MPF. The department of engineering mechanics of the faculty of aerospace engineering of Delft University of Technology developed this framework for multi-physics finite-element discretizations. The framework is based on two C++ libraries, JEM and JIVE from Habanera numerical software. JEM is a C++ toolkit providing a set of classes and functions that facilitate the development of robust, modular software. And JIVE is a programming toolkit aimed at solving partial differential equations. It provides a set of classes, functions and data structures for transforming a partial differential equation into a system of equations; for solving such a system; and for computing quantities derived from the solution.
References
[1] Thomas J.R. Hughes Arif Masud. A space-time galerkin/least-squares finite element formulation of the navier-stokes equations for moving domain problems. Computer Methods in Applied Mechanics and Engineering, 146(91-126), 1997.
[2] Chi-Wang Shu Bernardo Cockburn, George E. Karniadakis. The development of discontinuous galerkin methods. Discontinuous Galerkin methods, theory computation and applications, Lecture notes in computational science and engineering, 11(3-50), 2000.
[3] T.J.R. Hughes F. Shakib and Z. Johan. A new finite element formulation for computational fluid dynamics: X. the compressible euler and navier-stokes equations. Computer Methods in Applied Mechanics and Engineering, 89(141-219), 1991.
[4] Spencer J. Sherwin George Em Karniadakis. Spectral/hp element methods for CFD. Oxford University Press, Inc., 1999.
[5] Christian Michler. Efficient numerical methods for fluid-structure interaction. PrintPartners Ip-skamp, Enschede, 2005.
[6] S.J. Sherwin-J. Peiro O.C. Zienkiewicz, R.L. Taylor. On discontinuous galerkin methods. Inter-national Journal for Numerical Methods in Engineering, 58(1119-1148), 2003.
[7] S.P. Spekreijse. Multigrid solution fo the steady Euler equations. Stichting Mathematisch Centrum, Amsterdam, 1988.
