Goal-Oriented A-Posteriori Error Estimation and Adaptivity for Fluid-Structure Interaction

(1)

Delft Aerospace Computational Science 0000 0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111

000000000

111111111

ENGINEERING MECHANICS

REPORT

May 2008

Goal-Oriented A-Posteriori Error

Estimation and Adaptivity for

Fluid-Structure Interaction

P.W. Fick, K.G. van der Zee and E.H. van Brummelen P.W.Fick@tudelft.nl

(2)

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

(3)

Abstract

Numerical simulation of fluid-structure interaction generally requires vast computational resources. Para-doxically, the computational work is dominated by the complexity of the subsystem that is of least practical interest, viz. the fluid. The resolution of each of the many small-scale features in the fluid is prohibitively expensive. However, quantitative concern is generally restricted to the structural response. Goal-oriented adaptation strategies provide a paradigm to bypass this paradox. The crucial point is that in order to reliably determine the structural response, it is generally not necessary to fully resolve all the small-scale features in the fluid subsystem. Over-all computational cost can therefore be significantly reduced by resolving only those features that bear a pronounced influence on a specific structural output quantity of interest. In the context of fluid-structure interaction, meaningful output quantities are, for instance, the space-time average of the structural displacement and the nett energy that is transferred from the fluid to the structure. To enable the design of finite element meshes specifically tailored to the efficient computation of one such target quantity, an a posteriori error estimate must be obtained relating the error in the quantity of interest to local errors in the computed solution. This is accomplished by approximately solving an appropriate dual problem. The solution to this problem can be used to identify the sensitivity of the target-quantity error to local discretization residuals. Using this information, it is then possible to construct goal-oriented error indicators, providing the basis for an optimal-adaptive mesh refinement strategy, capable of delivering efficient and reliable control of the error in the target quantity to within a user-defined tolerance.

In this thesis, we apply the existing framework for duality-based a posteriori error estimation and goal-oriented adaptivity to a prototypical fluid-structure-interaction model problem. Specifically, we consider the two-dimensional panel problem concerning the aeroelastic response of a flexible panel with infinite aspect ratio immersed in an inviscid fluid flow. For general linear and nonlinear output quantities of interest, we formulate an appropriate dual problem and derive dual-weighted Type I a posteriori error estimates. The sharpness of these estimates is demonstrated through a series of numerical experiments for physically stable as well as physically unstable test conditions in both the subsonic and supersonic regime. On the basis of the derived a posteriori error estimates, we then design and implement an adap-tive algorithm capable of producing space-time meshes specifically tailored to the efficient computation of a certain target quantity of interest. Numerical results are presented, highlighting the superiority of the proposed duality-based approach over a more traditional mesh refinement algorithm employing a residual-based indicator. Furthermore, comparisons between h- and hp-refinement strategies are made to illustrate the extra increase in efficiency, which can be gained from the use of hp-refinement techniques.

Keywords: Adaptive finite element methods, fluid-structure interaction, discontinuous Galerkin, space-time finite elements, a posteriori error estimation, adjoint consistency, goal-oriented adaptivity

(4)

(5)

Introduction

1.1 Motivation

Fluid-structure-interaction (FSI) phenomena play an important role in a variety of engineering disci-plines including aerospace, mechanical, civil and biomedical engineering. A traditional application area has been in aircraft design, where the flexible aircraft structure interacts with the surrounding airflow and may encounter violent flow-induced instabilities such as flutter and buffet. These instabilities restrict the admissible operating conditions of the aircraft and must therefore be controlled by a proper design. Outside aerospace engineering, examples of fluid-structure-interaction phenomena include wind-induced vibrations of buildings and bridges, dynamic instabilities of wind turbine blades and flow-induced vibra-tions in oil pipes. An important and interesting new application area has recently arisen in the field of biomechanics, where it is believed that many vascular diseases (for example astherosclerosis) and disor-ders (such as aneurysms) are related to fluid-structure-interaction phenomena. Moreover, fluid-structure interaction plays a crucial role in the development and design of biomedical prosthetic devices such as artificial heart valves and stents.

With such a diverse range of engineering applications, there is an increasing demand for accurate and reliable solution techniques to analyze the dynamic behavior of fluid-structure-interaction systems. Until now, this has been achieved primarily using highly simplified linearized models operating in the frequency domain. Because of the underlying assumptions, these tools do not possess a high fidelity and must there-fore be supplemented with expensive experimental tests to guarantee a safe and reliable design. Driven by the necessity to analyze more and more complicated problems, and fueled by the advances in com-putational science and parallel processing, there exists an increasing impetus in industry and academia to replace these expensive experimental tests by fully-coupled simulations in the time domain employ-ing more accurate nonlinear descriptions for flow and structure, see e.g. [14]. As these simulations are computationally very expensive, recent years have seen a rapid increase in the demand and development of robust and efficient computational methods to enable the effective use of such simulation tools for realistic engineering applications.

The development of computational methods for fluid-structure interaction is a challenging endeavor on account of three major difficulties, which are closely linked to the generic features of fluid-structure-interaction phenomena. First of all, fluid-structure fluid-structure-interaction constitutes a multiphysics problem as it requires the simultaneous solution of both a fluid and a structure problem. Hence, a coupling must be established between the fields of computational fluid dynamics (CFD) and computational solid dynam-ics (CSD). In order to accommodate the structural deformations, the fluid mesh must deform in time,

(8)

requiring complicated mesh deformation algorithms. The latter is often achieved by modeling the mesh as a pseudo-structural system using a spring network analogy, see e.g. [50, 42]. This then leads to a so-called three-field formulation for the fluid, the structure and the motion of the fluid mesh. Secondly, because the position of the fluid-structure interface is unknown a priori and needs to evolve from the solution of the coupled problem, fluid-structure interaction problems exhibit a free-boundary character, giving rise to a complicated interconnection between the governing initial-boundary-value problems, and the domains on which these are defined. Finally, fluid-structure interaction also constitutes a multiscale problem, owing to the wide range of length and time scales present in the fluid and structure subsystems. These disparate scales in fluid and structure translate into different resolution requirements for the dis-cretization and typically lead to severely ill-conditioned systems.

Faced by these issues, it may be clear that the numerical solution of fluid-structure interaction is pro-foundly more difficult than the solution of the fluid and structure separately. Accordingly, numerical simulation of FSI generally requires vast computational resources. Paradoxically, the computational work is dominated by the complexity of the subsystem that is generally of least practical interest, viz. the fluid.1 _{The resolution of each of the many small-scale features in the fluid is prohibitively expensive.} However, practical interest is typically restricted to the structural response. Goal-oriented adaptation strategies could provide a paradigm to bypass this paradox. The crucial point is that in order to reliably determine the structural response, it is generally not necessary to fully resolve all the small-scale features in the fluid subsystem. Over-all computational cost can therefore be significantly reduced by resolving only those features that bear a pronounced influence on the structural quantity of interest. Hence, this strongly motivates the use of optimal-adaptive discretization techniques, tailored to the specific goals of the computation.

1.2 Research objectives

Motivated by the considerations outlined above, we wish to investigate optimal-adaptive computational techniques for fluid-structure interaction. These techniques rely on a posteriori error estimates for out-put functionals J (·) as a function of the comout-puted solution uh, which express certain target quantities that are of particular engineering or scientific interest. Examples of target quantities that often arise in traditional aeronautical CFD applications include for instance the lift or drag around an airfoil, the nett mass flow through a turbomachinery cascade or the radiated noise from an engine intake. Turning toward the numerical simulation of fluid-structure interaction phenomena, we may select the output functional in connection with our interest for the structural response. Hence, a relevant target quantity of interest could for example be the (local) stress or strain in the structure or the average structural displacement. Alternatively, since the response of the system is governed by the energy balance between the two subsys-tems, another relevant functional would be an energy-transfer functional that expresses the nett energy transferred from the fluid to the structure.

Having established a proper target quantity, the aim of the optimal-adaptive strategy will then be to construct a mesh Th and corresponding finite element space Uh,p that optimally control the error |J (u) − J (uh)|, where “optimal” in this case can either mean “most economical for achieving a pre-scribed level of accuracy TOL” or “most accurate for a given computational cost”. To this end, we do not mind that errors in the small-scale features of the flow field are locally large when their contribution to the structural quantity of interest is only small. Accordingly, it shall be our goal to derive dual-weighted

1_{As demonstrated by Piperno and Farhat [14], for a realistic application (in this case the aeroelastic response of an F16}

aircraft) typically more than 98% of the total simulation time is consumed by the fluid subsystem, i.e. inside the fluid solver and the mesh update algorithm.

(9)

a posteriori error estimates of the form

J (u) − J (uh) ≈ hRp(uh) , ωh(z)i .

Here, ωh(z) acts as a weight factor that measures the sensitivity of the error in the output quantity of interest, J (uh), to the primal equation residual, Rp(uh). These weights are obtained by approximately solving an auxiliary dual (or adjoint) problem; see for instance the review papers [3, 20] and the refer-ences therein. Upon localization of the a posteriori error estimate, dual-weighted error indicators can be constructed providing the basis for our optimal-adaptive discretization strategy.

The aforementioned framework for a posteriori error estimation and goal-oriented adaptivity has seen a steady development in the past decade. Its performance has been demonstrated for several generic linear and nonlinear boundary-value problems in both fluid and solid mechanics. However, as pointed out by Becker and Rannacher [3], the realization of the developed concepts for time-dependent problems, though theoretically straightforward, is still in a relatively immature state. Moreover, the extension to multiphysics problems - such as fluid-structure interaction - is in its infancy and constitutes a major challenge for the further development of this framework. In particular, for fluid-structure interaction, a fundamental complication arises due to the free-boundary character of the fluid-structure interface, as this leaves the underlying fluid domain unknown a priori. Consequently, the formulation of an appropriate dual problem becomes a highly nontrivial task. Very recently, some methods have been introduced to accommodate this issue [9, 54, 55]. In this work, however, we choose to focus on the multiphysics aspect of fluid-structure interaction. Accordingly, to simplify the setup and analysis of the coupled problem, we will bypass the implications related to the free-boundary character of the fluid-structure interface by introducing a so-called geometric linearization using surface transpiration techniques.

As a prototypical model problem, we consider the two-dimensional panel problem, cf. [48]. This problem constitutes the aeroelastic response of a flexible panel with infinite aspect ratio immersed in an inviscid fluid flow. We discretize this problem using a space-time finite element formulation for both the fluid and structure, with discontinuous approximation spaces in time. The latter is advantageous, particu-larly because of the fact that the spatial mesh can then be conveniently changed from one time step to the next, which greatly facilitates the use of adaptive mesh refinement techniques to track and cap-ture local solution feacap-tures that propagate through the spatial domain. Subsequently, we formulate an appropriate dual problem and complete an adjoint consistency analysis. By construction, the discrete dual problem involves the transpose of the primal bilinear form. It will be shown that this leads to an adjoint-inconsistent formulation. To recover adjoint consistency, we propose a consistent modification of the kinematic coupling condition at the fluid-structure interface. Based on the solution of the so-obtained modified primal and dual problems, we then derive dual-weighted a posteriori error estimates for general linear and nonlinear output quantities of interest. The sharpness of these error estimates is demonstrated through a series of numerical experiments for physically stable as well as physically unstable solutions in both the subsonic and supersonic regime. On the basis of the derived a posteriori error estimates, we finally design and implement an adaptive algorithm that optimally controls the error in a specific output quantity of interest to within a user-defined tolerance. Comparisons with residual-based error indicators are used to highlight the superiority of the proposed duality-based approach over more traditional ad hoc strategies. Furthermore, we investigate the extra increase in efficiency, which can be gained from the use of hp-refinement techniques. In the latter case, the decision between h- and p-refinement is based on estimating the local smoothness of the unknown exact solution using analyticity and Sobolev regularity estimation techniques [32, 10, 29]. If the underlying solution is predicted to be locally smooth, then poly-nomial enrichment is used; else, if the solution is non-smooth, h-refinement is performed. Accordingly, regions in the computational domain where the solution is locally smooth are isolated from regions with limited regularity, thereby further improving the efficiency of the adaptive algorithm.

(10)

1.3 Thesis outline

The contents of this thesis are organized as follows. Following the general introduction in Chapter 1, Chapter 2 is concerned with the definition of the model problem and the corresponding mathematical formulation. Governing equations are introduced for the fluid, the structure and the kinematic and dynamic interface conditions. The equations are presented both in strong form and in weak form, and are finally rephrased into a single variational formulation in space-time for the aggregated fluid-structure-interaction system. Subsequently, a discretization of the model problem is presented in Chapter 3. Chapter 4 is concerned with duality-based a posteriori error analysis for general output quantities of interest. In particular, we will derive dual-weighted a posteriori error estimates, providing the basis for an optimal-adaptive mesh refinement algorithm. The performance of this algorithm, and the quality of the a posteriori error estimates are then studied in Chapter 5 through a series of numerical experiments. Finally, in Chapter 6 we provide a summary of our work and sum up the most important conclusions.

(11)

Chapter 2

Problem Statement

In this work we apply the framework of a posteriori error estimation and goal-oriented adaptivity to a prototypical fluid-structure-interaction model problem. This chapter is concerned with the definition of this model problem and the corresponding mathematical formulation.

The contents of this chapter are organized as follows. First, in section 2.1 we introduce the model problem, viz. the two-dimensional panel problem. In sections 2.2–2.4 we present the governing equations for the fluid, the structure and the interface conditions, both in strong form and in weak form. In section 2.5 these are then recast into a single variational formulation for the aggregated fluid-structure-interaction problem in space-time. Finally, some relevant target quantities of interest for our model problem are established in section 2.6.

2.1 The two-dimensional panel problem

As a model problem we consider the two-dimensional panel problem, cf. [48]. This problem constitutes the aeroealastic response of a flexible panel with infinite aspect ratio immersed in an inviscid fluid flow. The upper side of the panel is exposed to an airstream with uniform inflow conditions; the lower side is connected to a cavity with air at rest. For an illustration, see Figure 2.1. The freestream Mach number, speed of sound, pressure and density are denoted by M∞, a∞, p∞and ρ∞, respectively. The pressure in the cavity underneath the panel is denoted by β, which in this case is taken to be equal to the freestream pressure p∞.

The dynamic motion of the panel is driven by the pressure force exerted by the fluid. Depending on the system parameters, the dynamic behavior can be either stable or unstable. In the subsonic regime the panel may exhibit steady-state divergence or self-induced buckling; in the supersonic regime the critical instability mode is governed by flutter. For a more detailed elaboration of these phenomena we refer to [8]. The panel problem is a classical problem in the field of aeroelasticity and has been extensively investigated in the past, both numerically and experimentally. In particular, we refer to the paper by Piperno and Farhat (2001) [48], where the two-dimensional panel problem is introduced as a benchmark problem in FSI. We use this paper as a comparison to validate our method in the supersonic regime. In the low-speed subsonic regime, we make use of exact analytic expressions derived by Kornecki et. al. [39].

(12)

y x p|_Γ_s π (x, t) β α (x, t) w (x, t) Γs Γff Ωf

Figure 2.1: Illustration of the two-dimensional panel problem with an expanded view of the interface region.

2.2 Governing equations for the fluid

2.2.1 The linearized small-disturbance potential flow equation

We consider a fluid on the open space-time domain Qf := Ωf(t) × I, consisting of the spatial domain Ωf(t) := {(x, y) : −∞ < x < ∞ ; α (x, t) < y < ∞}, protruded in time over the interval I := ]0, T [. The boundary of Qf is denoted by ∂Qf and consists of the lower time boundary Ωf(0) := Ωf× {t = 0}, the up-per time boundary Ωf(T ) := Ωf× {t = T } and the spatial domain boundaries ∂Ωf(t) × I. Here, ∂Ωf(t) is composed of the fixed wall boundary Γw, the farfield boundary Γff and the moving fluid-structure interface Γs(t) := {(x, y) : 0 < x < L ; y = α (x, t)}. We remark that the definitions of Γs, Ωf and Qf all depend on the position of the fluid-structure interface, parameterized by α (x, t). Since, α (x, t) needs to evolve from the solution to the coupled problem, the problem domain Qf is unknown a priori. Consequently, the considered problem constitutes a free-boundary problem.

To bypass the implications related to the free-boundary character of the fluid-structure interface, we intro-duce a geometric linearization. This is accomplished by setting α (x, t) = α0(x) = 0, so that Γsbecomes the fluid-structure interface in undisplaced position. Accordingly, the fluid domain Ωf(t) does no longer deform in time, reducing the problem from a free-boundary problem to a much simpler fixed-boundary problem.1 _{The movement of the panel will now be modeled by so-called transpiration conditions. We} defer a further elaboration of these aspects to section 2.4.

On Qf, let us adopt a relatively simple description for the fluid. In particular, let us assume inviscid, irrotational, isentropic flow, allowing us to replace the velocity field by a gradient field of a scalar potential

1_{Apart from the fact that this greatly simplifies the solution of the coupled problem, it also places the problem in a more}

suitable framework for duality-based a posteriori error estimation, which has sofar been developed primarily for generic boundary-value problems.

(13)

function Φ, i.e. (u, v) = ∇Φ. Moreover, let us assume small perturbations to this potential of the form Φ = Φ∞+ φ and assume that we are well outside the transsonic regime, i.e. 1 − M2

∞ is large. Under these assumptions, the governing relation for the perturbation potential φ is given by the linearized small-disturbance potential flow equation (see for example [47])

∂2_φ ∂t2 + 2U∞ ∂2_φ ∂x∂t− a 2 ∞ 1 − M∞2 ∂2_φ ∂x2 − a 2 ∞ ∂2_φ ∂y2 = 0 (x, y, t) ∈ Qf , (2.1) subject to suitable initial and boundary conditions, dependent on the specific settings of the prob-lem, see section 2.2.2. The velocity field can be recovered from the solution φ in the following way: (u, v) = (U∞+ φx, φy). Furthermore, the pressure p is given by

p = p∞− ρ∞ ∂φ ∂t + U∞ ∂φ ∂x . (2.2)

In the sequel, we will consider these equations in nondimensional form. To this end, let us introduce the nondimensional variables ˆx = x/L, ˆy = y/L, ˆt = a∞t/L, ˆφ = φ/ (a∞L) and ˆp = p/ ρ∞a2

∞

. The nondimensional counterparts of (2.1)–(2.2) are then as follows:

∂2_φ_ˆ ∂ˆt2 + 2M∞ ∂2_φ_ˆ ∂ ˆx∂ˆt − 1 − M 2 ∞ ∂2_φ_ˆ ∂ ˆx2 − ∂2_φ_ˆ ∂ ˆy2 = 0 ˆx, ˆy, ˆt ∈ ˆQf , (2.3) ˆ p = ˆp∞− ∂ ˆφ ∂ˆt + M∞ ∂ ˆφ ∂ ˆx ! . (2.4)

Equation (2.3) constitutes a second-order hyperbolic equation. The discontinuous Galerkin (DG) method provides a suitable technique for discretizing this class of partial differential equations. Accordingly, in order to enable a DG approach (see section 3.1), it is convenient to reformulate the problem as a system of first-order hyperbolic equations, cf. [13]. For that purpose, we introduce the new set of variables q:= (∂ ˆφ/∂ˆt, ∂ ˆφ/∂ ˆx, ∂ ˆφ/∂ ˆy). In terms of these new variables, equation (2.3) can be rewritten as

∂q1 ∂ˆt + 2M∞ ∂q1 ∂ ˆx − 1 − M 2 ∞ ∂q2 ∂ ˆx − ∂q3 ∂ ˆy = 0 x, ˆˆ y, ˆt ∈ ˆQf . (2.5a)

This equation must be supplemented with two additional relations, for which we may take (a combination of) identities of the form

∂qi ∂ ˆxj −

∂qj

∂ ˆxi = 0 ; i 6= j .

We remark that the choice of constraint relations affects the well-posedness of the problem. In fact, it has been observed that the most straightforward choice

∂q2 ∂ˆt − ∂q1 ∂ ˆx = 0 , ∂q3 ∂ˆt − ∂q1 ∂ ˆy = 0

leads to an ill-posed problem in the supersonic regime (M∞> 1), presumably due to the fact that one of the characteristic directions is then rendered to be coincident with the orientation of the spatial domain boundaries. In other words, one of the characteristic speeds is equal to zero, preventing information from the corresponding characteristic boundary condition to propagate into the domain. This problem can be

(14)

avoided by choosing a slightly different set of constraint relations. Without further rigorous mathematical foundation, in the sequel we consider the choice

∂q2 ∂ˆt − ∂q1 ∂ ˆx = 0 , (2.5b) ∂q3 ∂ˆt − ∂q1 ∂ ˆy + M∞ ∂q3 ∂ ˆx − ∂q2 ∂ ˆy = 0 . (2.5c)

Supplementing the governing equation (2.5a) with the constraint relations (2.5b)-(2.5c) yields the first-order system ∂q ∂t + Ax ∂q ∂x + Ay ∂q ∂y = 0 , (2.6) with Ax=   2M∞ M 2 ∞− 1 0 −1 0 0 0 0 M∞   , Ay =   00 00 −10 −1 −M∞ 0   ,

where, for the purpose of notational convenience, we suppressed the hat superscripts. Regarding the wellposedness of (2.6), we close with the following remarks.

Remark 2.1 (Analogy with Friedrichs systems). Wellposedness results for systems of the form (2.6) are currently only existent under the restriction that the matrices Axand Ayare symmetric (with further restrictions on the boundary conditions), in which case (2.6) reduces to a so-called Friedrichs system, cf. [16]. Unfortunately, this only holds for the special case M∞ = 0, for which (2.6) reduces to the wave equation (in first order-form); for M∞ > 0 the matrices Ax and Ay become nonsymmetric and, accordingly, wellposedness cannot be guaranteed.

Remark 2.2 (Symmetrization). Using the coordinate transformation   xy t   =   10 0 M∞1 0 0 0 1     ξη τ   ⇔   ξη τ   =   1 00 1 −M∞0 0 0 1     xy t   (2.7)

equation (2.3) can be conveniently mapped into the second-order wave equation, for which wellposedness follows directly from the existing literature, see e.g. [44]. However, this property is not automatically inherited by the first-order counterpart of (2.3), i.e. applying the same coordinate transformation to (2.6) does not automatically yield the wave equation in first-order form. The outcome of the coordinate transfor-mation obviously depends on the matrices Axand Ay, and therefore on the choice of constraint relations. In fact, it can easily be shown that the only choice for which the coordinate transformation (2.7) yields a system that is symmetric and equivalent to the wave equation in first-order form, is the choice (2.5b)-(2.5c). In this case the Jacobian of the transformation matrix can be regarded as a right-symmetrizer, comparable to the left-symmetrizers considered by Friedrichs and Lax to analyze the well-posedness of nonlinear conservation laws [17]. A wellposedness analysis for (2.6) should probably follow along these lines. However, this is outside the scope of the present work.

2.2.2 Initial and boundary conditions

To complete the description of the initial-boundary-value problem for the fluid, equation (2.6) must be supplemented with suitable initial and boundary conditions. A necessary condition for well-posedness is that the number of boundary conditions is equal to the number of inbound characteristics, see e.g. [40]. To allow a clear and precise treatment of the boundary conditions, let us first introduce some necessary notation.

(15)

Preliminaries

For x ∈ ∂Qf, let n∂Qf(x) := (nx, ny, nt) denote the outward unit normal vector to ∂Qfat x and consider the flux Jacobian

An(n) := ntI + nxAx+ nyAy ,

where I denotes the 3 × 3 identity matrix. Assuming that (2.6) is hyperbolic, matrix An can be diagonalized, such that we obtain

An= RΛR−1,

where R is a matrix whose columns contain the right eigenvectors of An and Λ is the corresponding diagonal matrix of real eigenvalues λi. These eigenvalues express the characteristic speeds of the Riemann invariants and can be shown to be

λ1= nt+ nxM∞−qn2

x+ n2y , λ2= nt+ nxM∞ , λ3= nt+ nxM∞+ q

n2

x+ n2y . (2.8) The matrix Λ can be decomposed in a negative semidefinite part and a positive semidefinite part, i.e.

Λ = Λ−+ Λ+ , where Λ− = diag λ−i

and Λ+ = diag λ+i

with λ−i = min (0, λi) and λ+i = max (0, λi). With this notation, we now define the matrices

A−n = RΛ−R−1, A+n= RΛ+R −1_, such that An= A+n+ A − n , and accordingly Anq= A+nq −_{+ A}− nq + _, _(2.9)

where q− _{and q}+_{denote the inner and outer trace on ∂Qf, i.e. q}±_{= limǫ→0}_{(q ± ǫn∂Q}

f). We have thus obtained a decomposition of the flux vector Anq on the domain boundary ∂Qf, consisting of an interior contribution from outbound characteristics A+nq− and an exterior contribution from inbound character-istics A−nq+. Now recall that we need to specify a boundary condition for each inbound characteristic. Essentially, this means that we need to specify the exterior contribution A−nq+ by introducing suitable boundary conditions for the Riemann characteristics with negative λi. These boundary conditions can be expressed as

BqB= e (x, y, t) ∈ ∂Qf , (2.10)

where the dimensions of the boundary condition operator B depend on the number of inbound charac-teristics and may therefore be different on different parts of the boundary. This boundary condition will hereafter be referred to as the physical boundary condition.

In agreement with the notational conventions for Friedrichs’ systems (cf. [16]), it is common practice to rewrite the boundary condition (2.10) as a contribution to the flux Anq, leading to a corresponding characteristic boundary condition (this will proof to be useful later on in the definition of the variational problem). To this end, we need to solve an incomplete Riemann problem at the boundary. It is to be noted that the Riemann invariants are given by qc= R−1q= Sq = (S+/S−) q, where (S+/S−) represents the vertical concatenation of S+ and S−, with S+ and S− being the parts of S = R−1corresponding to the Riemann invariants of the outbound and inbound characteristics, respectively. Retaining the m Riemann invariants q+

(16)

3 − m physical boundary conditions of the form (2.10) for the inbound characteristics with negative λi, then leads to the following modified Riemann problem for the boundary state q:

S+ B q= S+ 0 q−+ 0 B qB .

When this problem is linear, i.e. when S+does not depend on q, and for a suitable choice of the physical boundary condition operator B, we can directly solve this expression for q. Accordingly, the boundary condition (2.10) can be weakly imposed through the flux Anqby using

Anq = An S+ B −1 S+ 0 q−+ 0 B qB = An S+ B −1 1 2 S+ B + S+ −B q−+1 2 S+ B − S+ −B qB = 1 2(An+ N ) q −₊1 2(An− N ) q B _, _(2.11)

with the boundary operator

N = An S+ B −1 S+ −B .

Consequently, an alternative expression for the physical boundary condition (2.10) is

(An− N ) q − qB= 0 (x, y, t) ∈ ∂Qf , (2.12) which will hereafter be referred to as the characteristic boundary condition. Having introduced the necessary notation, let us now proceed with the specification of the physical and characteristic boundary conditions on the different parts of the boundary ∂Qf.

Wall boundary conditions

On the fixed wall boundary, Γw×]0, T [, and the fluid-structure interface, Γs×]0, T [, we have n = (0, −1, 0). In accordance with (2.8), there is only one inbound characteristic, for which we need to specify a boundary condition. A natural choice would be to impose flow tangency, giving rise to the physical boundary conditions:

0 0 1 q= 0 on Γw× ]0, T [ , (2.13)

0 0 1 q= ǫ (x, t) on Γs× ]0, T [ . (2.14) Here, ǫ constitutes a so-called transpiration velocity which is induced by the movement of the fluid-structure interface; see section 2.4 for further details. The corresponding characteristic boundary condi-tions are of the form (2.12) with

An= −Ay=   00 00 10 1 M∞ 0   , Nw= Ns=   00 00 −10 1 M∞ 2   and qB=   00 0   on Γw× ]0, T [ , qB=   00 ǫ   on Γs× ]0, T [ .

(17)

Farfield boundary conditions

At the farfield boundary, Γff×]0, T [, it is generally far less trivial what boundary condition to specify. For a numerical treatment it is necessary to truncate the spatial domain Ωf by introducing so-called artificial boundaries in the farfield. These artificial boundaries are only a computational necessity; they have no physical significance. This makes it difficult to choose an appropriate boundary condition. Of course, one could decide to impose asymptotic farfield conditions at the truncated domain boundaries. However, this will produce spurious (i.e. non-physical) reflections of outgoing waves, giving rise to large errors in the computed solution. Instead, a much better approach would be to employ an absorbing boundary condition that limits these spurious reflections. The existing literature in this field provides a variety of methods to construct such an absorbing boundary condition, see e.g. the review papers [22, 33] and the references therein. In the present work we confine our interest to local low-order non-reflecting boundary conditions (NRBCs).

The derivation of such local NRBCs is based on Fourier analysis. The general idea is to assume wave-like solutions at the boundary of the form

q(nx, ny, t) = ˆqei(k(nxx+nyy)+l(−nyx+nxy)−ωt)_.

Here, k and l represent the wavenumbers in normal and transversal direction, respectively, and ω denotes the angular frequency. Substituting this expression into the governing equation (2.6) yields the following dispersion relation

(−ωI + k (nxAx+ nyAy) + l (−nyAx+ nxAy)) ˆq= 0 .

Using the roots of this dispersion relation, we can separate the inbound waves from the outbound waves. Setting the amplitudes of the Fourier modes corresponding to the inbound waves equal to zero, and subsequently making a local backward Fourier transform using rational approximations, a sequence of non-reflecting boundary conditions can be derived of increasing accuracy, for details see the classical pa-per by Engquist and Majda [11]. In the following, we consider the most simple case and derive an NRBC of order p = 1. We proceed along the lines of [19], where the general approach outlined above is used to derive local non-reflecting boundary conditions for first-order multidimensional hyperbolic systems. Let us assume sinusoidal plane waves traveling in the direction normal to the boundary. Hence, we neglect perturbations in the transversal direction by setting l = 0. The Fourier decomposition of q and the dispersion relation then reduce to

q(nx, ny, t) = ˆqei(k(nxx+nyy)−ωt) _, (−ωI + k (nxAx+ nyAy)) ˆq= 0 .

Let rR_{and r}L _{denote the right and left null-vectors of (−ωI + k (nxAx}_{+ nyAy)). Remark that r}R_and rL are also the right and left eigenvectors of (nxAx+ nyAy) corresponding to the eigenvalue ω/kn with orthogonality property rL

n(ω/kn) · rRm(ω/km) = 0 ∀ n 6= m. At the boundary, the solution q can now be decomposed into a sum of Fourier modes with different values for the angular frequency ω. Considering just one particular choice for ω (thus assuming a constant-frequency wave), we obtain

q(nx, ny, t) = " ₃ X n=1 anrRneikn(nxx+nyy) # e−iωt .

(18)

an inbound wave. Note that because of orthogonality we have rL_nq = rL n " ₃ X m=1 amrR meikm(nxx+nyy) # e−iωt = an rLn· rRn eikn(nxx+nyy)_e−iωt _. Hence, an equivalent expression for the non-reflecting boundary condition is

rL_nq= 0 ∀n ∈ {[1, 3] : λn< 0} . (2.15) The accuracy of this NRBC can be expressed by a reflection coefficient, relating the amplitudes of the spuriously reflected inbound waves to the amplitudes of the outbound waves. It can be shown, see [19, section 2.9], that the reflection coefficient is O(l/ω). An interesting observation is that the spurious reflections are theoretically equal to zero for l = 0, i.e. for sinusoidal plane waves traveling in the direc-tion (nx, ny). Hence, an opportunity exists to make the boundary outward normal - as perceived by the NRBC - a function of the evolving solution by adapting it to the local direction of the solution gradient, giving rise to a so-called direction-adaptive NRBC, in the fashion of [45]. In the sequel, however, we do not consider this option any further.

If the time component of the outward boundary normal is equal to zero, then the left and right eigenvectors rL_{and r}R_{are simply equivalent to the rows and columns of R}−1_{and R, respectively. The NRBC in (2.15)} can then be conveniently re-expressed in terms of the local Riemann invariants, yielding:

q−_c = S−q= 0 on Γff× ]0, T [ . (2.16)

The corresponding characteristic boundary condition is of the form (2.12) with Nff = RΓ (S+\ −S−) and qB_{= 0, which can be reduced to}

A−nq= 0 on Γff × ]0, T [ . (2.17)

Remark 2.3 (Higher-order NRBCs). We have restricted our interest here to a local low-order NRBC of order p = 1. Higher-order NRBCs can be devised in a similar fashion, but their implementation becomes problematic beyond p = 2 due to the presence of higher-order derivatives. Recently, methods have been devised to bypass these complications by replacing the higher-order derivatives with specially defined auxiliary variables, see e.g. [23, 24]. In view of the current trend towards higher-order adaptive discretizations, we would like to emphasize that the use of such high-order NRBCs is becoming more and more relevant. After all, an exclusive focus on the reduction of discretization errors is not entirely justified for realistic engineering applications; other sources of errors, such as modeling errors and boundary truncation errors, are equally important and need to be addressed as well.

Remark 2.4 (Grid stretching). The NRBC presented above can be perfectly combined with a grid stretching technique, where the computational mesh is coarsened to create an absorbing buffer zone between the domain of interest close to fluid-structure interface and the artificial farfield. Any outgoing waves and spuriously reflected incoming waves are then be automatically annihilated by this buffer zone due to the numerical dissipation introduced by the coarse mesh. Within an adaptive framework, this can be easily incorporated by starting from a very coarse mesh and refining only in the physical domain of interest close to the fluid-structure interface.

Initial conditions

Finally, we must also specify a boundary condition at the lower time boundary Ωf(0). Here, we have three inbound characteristics, such that A+n = 0 and A−n = I. Accordingly, both the physical and characteristic boundary condition reduce to:

(19)

2.2.3 Functional setting

A functional setting for the initial-boundary-value problem of the fluid subsystem is provided by the graph space of the partial differential operator Lf:= I∂t+ Ax∂x+ Ay∂y, i.e. the space

H (Lf, Qf) := {q ∈L2(Qf)3: Lfq∈L2(Qf)3} , equipped with the norm

kqk2Lf := kqk 2 [L2_(Q f)]3+ kLfqk 2 [L2_(Q f)]3 ;

cf. for example [31]. In summary, the fluid subproblem can now be expressed as follows: find q ∈ H (Lf, Qf) such that (2.6) is satisfied, subject to the initial condition (2.18) and the boundary condi-tions (2.13), (2.14) on Γwand Γs, and (2.16) on Γff.

2.2.4 Variational space-time formulation

To cast the fluid subproblem into a variational formulation in space-time, we multiply (2.6) with appro-priate test functions and integrate the resulting expression over the space-time domain. For this purpose, we define the trial space Q := H (Lf, Qf) and the test space V := H (Lf, Qf). Moroever, we introduce the space E on Γs× ]0, T [ to accomodate the transpiration condition ǫ (x, t). Upon integrating-by-parts and introducing weak enforcements of the initial and boundary conditions, we then obtain the following variational statement for the fluid subproblem: given ǫ ∈ E, find q ∈ Q such that

F (q, v) = f (v, ǫ) ∀ v ∈ V , (2.19a)

with F : Q × V 7→ R and f : V 7→ R defined as F (q, v) := Z T 0 Z Ωf q· L∗fvdΩ dt + Z Ωf(T ) q· v dΩ + Z T 0 Z Γw 1 2(An+ Nw) q · v dΓ dt + Z T 0 Z Γs 1 2(An+ Ns) q · v dΓ dt + Z T 0 Z Γff A+nq· v dΓ dt , (2.19b) f (v, ǫ) := Z Ωf(0) q0· v dΩ − Z T 0 Z Γs 1 2(An− Ns) (0, 0, ǫ) T · v dΓ dt , (2.19c) where L∗

f is the formal adjoint of Lf, defined by L∗f := −I∂t− A T

x∂x− ATy∂y.

2.3 Governing equations for the structure

2.3.1 The Euler-Bernoulli beam equation

For the mathematical formulation of the structure subproblem we confine ourselves to a geometrically linear description of the kinematics. We consider the open bounded space-time domain Qs:= Γs× ]0, T [. The boundary ∂Qs consists of the lower time boundary Γs(0) := Γs× {t = 0}, the upper time boundary Γs(T ) := Γs× {t = T } and the spatial boundary ∂Γs:= {(x, t) ∈ {0, L} × ]0, T [ }. On Qs, the motion of the structure is governed by the Euler-Bernoulli beam equation

M∂ 2_w ∂t2 + ∂2 ∂x2 EI∂ 2_w ∂x2 = P + (β − π) (x, t) ∈ Qs. (2.20) Here, w designates the displacement of the beam from its equilibrium position, and M , E and I denote the beam’s mass, modulus of elasticity and second moment of area, respectively. The right-hand side

(20)

comprises a forcing term, which is composed of the pressure π exerted by the fluid on the structure, the constant pressure β in the cavity underneath the panel, and an external force P , which can for instance be imposed to apply an initial perturbation to the aeroelastic system. In the sequel, we assume that the bending stiffness EI is constant and take β = p∞.

We consider the governing equation for the structure in nondimensional form. To this end, let us introduce the nondimensional variables ˆx = x/L, ˆw = w/L, ˆt = a∞t/L, ˆP = P/ ρ∞a2

∞ , ˆβ = β/ ρ∞a2 ∞ and ˆ π = π/ ρ∞a2 ∞

. The nondimensional counterpart of (2.20) then becomes ∂2_w_ˆ ∂ˆt2 + λ 2∂4wˆ ∂ ˆx4 = µ ˆ P + ˆβ − ˆπ x, ˆˆ t∈ ˆQs, (2.21) with λ = La∞ M1/2_L2_(EI)−1/2 , µ = ρ∞L M .

Here, the parameter λ can be identified as the ratio of characteristic time scales of the fluid and the structure. Furthermore, the parameter µ expresses the fluid-structure mass ratio. Together with the freestream Mach number, M∞, these parameters characterize the behavior of the coupled system.

2.3.2 Initial and boundary conditions

To complete the description of the initial-boundary-value problem for the structure subsystem, equa-tion (2.20) must be supplemented with suitable initial and boundary condiequa-tions. We consider a panel that is clamped at both ends, so that the boundary conditions are given by

w (0, t) = w (L, t) = 0 , ∂w

∂x (0, t) = ∂w

∂x (L, t) = 0 t ∈ ]0, T [ , (2.22) where, for the purpose of notational convenience, we again suppressed the hat superscripts. As initial conditions, we specify an initial displacement w0 _{and an initial velocity w}1_{, i.e.}

w (x, 0) = w0(x) , ∂w

∂t (x, 0) = w

1_(x) _{x ∈ ]0, L[ .} _(2.23)

2.3.3 Functional setting

To provide a functional setting for the problem (2.21)–(2.23), we introduce the Sobolev space Hm_(Γs) of functions on Γs, for which the function itself, as well as the weak derivatives up to and including order |α| = m reside in L2_{. In particular, we consider the special Sobolev space H}m

0 (Γs) of functions in Hm_{(Γs), for which both the function and the weak derivatives of order |α| ≤ m − 1 vanish on the} boundary ∂Γs, i.e.

Hm

0 (Γs) := {w ∈ Hm(Γs) : w|∂Γs = D α_w|∂Γ

s = 0 , |α| ≤ m − 1} .

Evolution problems – such as the one in (2.21)–(2.23) – are naturally set in Bochner spaces. Denoting by X a Hilbert space with norm k·kX, the Bochner space L2_{(0, T ; X) comprises all functions u : ]0, T [ 7→ X} such that kukL2_{(0,T ;X)}< ∞ with

kukL2(0,T ;X)k2:= Z T

0

kuk2Xdt .

The natural solution space for problem (2.21)–(2.23) can now be defined as the Bochner space L2 _{0, T ; H}2

0(Γs)

, with further conditions on the velocity ˙w := ∂w/∂t and the acceleration ¨w := ∂2_w/∂t2_; see the proposition below.

(21)

Proposition 2.1 (Wellposedness). Letting f := P + (β − π) ∈ L2 _{0, T ; L}2_(Γs)_{, w}0 _{∈ H}2

0(Γs) and w1 _{∈ L}2_{(Γs), then there exists a unique solution w ∈ L}2 _{0, T ; H}2

0(Γs)

with ˙w ∈ L2 _{0, T ; L}2_(Γs) satisfying (2.21)–(2.23). In addition, it holds that ¨w ∈ L2 _{0, T ; H}−2

0 (Γs)

. Proof. The proof follows directly from theorems 8.1 and 8.2 in Ref. [44].

In summary, the structure subproblem can now be expressed as follows: find w ∈ L2 _{0, T ; H}2 0(Γs)

with ˙

w ∈ L2 _{0, T ; L}2_(Γs)_{such that (2.21) and (2.23) are satisfied.}

2.3.4 Variational space-time formulation

To cast the structure subproblem into a variational space-time formulation, we define the trial space W:= {w ∈ L2 0, T ; H02(Γs)

: ˙w ∈ L2 0, T ; L2(Γs), ¨w ∈ L2 0, T ; H0−2(Γs)

} and the test space H := {θ ∈ L2 _{0, T ; H}2

0(Γs)

. Moreover, we introduce the space P to accommodate the forcing term π. Upon integrating-by-parts, we then obtain the following variational statement for the structure subproblem: given π ∈ P, find w ∈ W with w = w0 _{and ˙}_{w = w}1 _{on Γs}_{(0), such that}

S (w, θ) = s (θ, π) ∀ θ ∈ H (2.24a)

with S : W × H 7→ R and s : H 7→ R defined as S (w, θ) := Z T 0 Z Γs ¨ wθ + λ2∆w∆θdx dt , (2.24b) s (θ, π) := Z T 0 Z Γs µ (P + p∞− π) θ dx dt . (2.24c)

We remark that the boundary conditions given by (2.22) are automatically embedded in the definition of the approximation space W. The initial conditions (2.23) are weakly enforced.

2.4 Interface conditions

The fluid and structure are interconnected at the fluid-structure interface by kinematic and dynamic interface conditions. In a formal setting, these are provided by

α (x, t) = w (x, t) (x, t) ∈ Γs(t) × ]0, T [ , (2.25a) v|_Γ_s = u|_Γ_s ∂w ∂x + ∂w ∂t (x, t) ∈ Γs(t) × ]0, T [ , (2.25b) π (x, t) = p|_Γ_s (x, t) ∈ Γs(t) × ]0, T [ . (2.25c) Condition (2.25a) relates the displacement of the structure, w (x, t), to the position of the interface, α (x, t). Flow tangency is imposed by condition (2.25b). This condition constitutes a so-called slip boundary condition for the fluid subproblem and translates into the impermeability of the interface; see Figure 2.2. Finally, equilibrium of forces at the fluid-structure interface is prescribed by condition (2.25c). Note that in this formal setting the interface conditions are imposed on the moving fluid-structure inter-face Γs, parameterized by α (x, t). However, α (x, t) (and therefore the position of Γs) is unknown a priori and needs to evolve from the solution of the coupled problem. Consequently, this gives rise to a so-called shape-dependent variational formulation, which is inherently noncanonical. This places the problem in a difficult framework with respect to duality-based a posteriori error estimation, as the formulation of an

(22)

∂w ∂x u|_Γ_s u|_Γ_s ∂w ∂x ∂w ∂t

Figure 2.2: Illustration of the definition of the interface impermeability condition.

appropriate dual problem becomes highly nontrivial; see [54, 55].

Therefore, to bypass these implications, we remove the free-boundary character of the fluid-structure interface by introducing a geometric linearization. This is accomplished by setting α (x, t) = 0 in (2.25a). Accordingly, the fluid domain does no longer deform, which greatly simplifies the coupled problem. This will, however, lead to a violation of the kinematic condition (2.25a), which is justified as long as the structural displacements are small. In fact, this is a very reasonable assumption, since it would be pointless to account for nonlinear effects of large interface displacements when considering a fluid description based on linearized small disturbance potential theory. Hence, we drop condition (2.25a). The two remaining interface conditions, (2.25b) and (2.25c), are applied at the undisplaced fluid-structure interface, yielding:

v|_Γ_s = u|_Γ_s ∂w ∂x +

∂w

∂t (x, t) ∈ Γs× ]0, T [ , (2.26a) π (x, t) = p|_Γ_s (x, t) ∈ Γs× ]0, T [ . (2.26b) In terms of the perturbation potential function φ and the interface transpiration velocity ǫ, these condi-tions (in nondimensional form) translate into

ǫ (x, t) = M∞∂w ∂x + ∂w ∂t (x, t) ∈ Γs× ]0, T [ , (2.27a) π (x, t) = p∞− ∂φ ∂t Γs + M∞ ∂φ ∂x Γs ! (x, t) ∈ Γs× ]0, T [ , (2.27b)

where we made the assumption ∂φ/∂x ≪ M∞. The movement of the interface is now modeled by surface transpiration, whereby condition (2.27a) acts as a so-called transpiration boundary condition for the fluid subproblem, with qB _{= (0, 0, ǫ) at Γs}_{× I. The transpiration method – inspired on the} Lighthill analogy for modeling the boundary layer displacement thickness in panel methods [43] – is a commonly used technique in computational fluid dynamics for simulating surface displacements without actually deforming the computational domain. Applications of this method can be found in the field of aerodynamic shape optimization [35, 56] and unsteady aeroelastic simulation [15]. Numerical and experimental validation has shown that for small to moderate deformations (in the appropriate norm) the surface transpiration technique produces accurate and reliable results [51].

(23)

2.5 Aggregated variational space-time formulation

Having introduced the governing equations for the fluid, the structure and the interface conditions, we are now ready to cast these equations into a single variational space-time formulation for the aggregated fluid-structure-interaction problem. To this end, let us enforce the kinematic impermeability condition by replacing ǫ in the right-hand side of (2.19) by (2.27a). Analogously, we enforce the dynamic interface con-dition by replacing π in (2.24) by (2.27b). Furthermore, let us write ˜π (q) := π (q)− p∞= − (q1+ M∞q2) to denote the nett traction exerted by the fluid on the structure, and let us define the aggregated set of variables x := (q, w) in the product trial space X := Q × W and y := (v, θ) in the product test space Y := V × H. The aggregated variational space-time formulation for the coupled problem can now be expressed as follows: find x ∈ X with w = w0 _{and ˙}_{w = w}1 _{on Γs}_{(0), such that}

P (x, y) = p (y) ∀ y ∈ Y , (2.28a)

with P : X × Y 7→ R and p : Y 7→ R defined as P (x, y) := F (q, v) + Z T 0 Z Γs 1 2(An− Ns) (0, 0, ǫ (w)) T · v dΓ dt + S (w, θ) + Z T 0 Z Γs µ˜π (q) θ dx dt , (2.28b) p (y) := Z Ωf(0) q0· v dΩ + Z T 0 Z Γs µP θ dx dt . (2.28c)

Here, the bilinear forms F (q, v) and S (w, θ) are defined as in (2.19b) and (2.24b), respectively. Moreover, note that the terms involving the kinematic and dynamic interface conditions are no longer part of the data and have been moved from the right-hand side of the segregated variational expressions for the fluid and the structure to the left-hand side of the aggregated variational formulation. The variational formulation presented above will serve as a basis for the subsequent chapters.

2.6 Target quantities of interest

For the model problem outlined above, our main interest will be directed toward the evaluation of specific (bounded and differentiable) target quantities, J (·), as a function of the solution, x := (q, w). In the context of fluid-structure interaction, we note that quantitative concern is generally restricted to the structural response. Accordingly, a relevant target quantity for our model problem can for example be the space-time average of the structural displacement, w, over the time interval ]0, T [ , expressed by

J (q, w) = Z T 0 Z Γs w dΓ dt . (2.29)

Alternatively, as a measure for the maximum stress or strain in the panel, we can consider the space-time average of the second-order spatial derivative of w, giving rise to

J (q, w) = Z T 0 Z Γs ∆w dΓ dt . (2.30)

In most cases, however, overall concern is directed towards the accurate prediction of transient and/or asymptotic structural instabilities, such as divergence and flutter. These instabilities are governed by the overall energy balance across the fluid-structure interface Γs, i.e. the nett energy that is transferred

(24)

from the fluid to the structure. In this context, the fluid can be regarded as an open energy reservoir that is interacting with the structure. Clearly, when the nett energy extracted from this reservoir by the structure is positive, structural excitations are amplified leading to an unstable response (note that there is no structural dissipation for the present model problem). To accurately predict the occurrence of such instabilities, we are therefore interested in the accurate representation of the energy balance at the fluid-structure interface. This interest can be reflected through an energy-transfer functional expressing the nett energy transferred from the fluid to the structure over the time interval ]0, T [ , i.e.

J (q, w) = Z T 0 Z Γs (p∞− π (q))∂w ∂t dΓ dt = Z T 0 Z Γs (q1+ M∞q2)∂w ∂t dΓ dt . (2.31)

Equivalently, based on (2.27a), one can also consider J (q, w) = Z T 0 Z Γs (q1+ M∞q2) q3− M∞∂w ∂x dΓ dt . (2.32)

In practice, the value of the target quantity J (q, w) cannot be evaluated exactly and must be estimated, based on some numerical approximation (qh, wh) ∈ Q × W to the problem (2.28). A measure of the quality of the approximation is provided by the error

E_J _{:= J (q, w) − J (q} h, wh) .

(25)

Chapter 3

Space-Time Finite Element

Discretization

In this chapter we present a finite element formulation for the previously introduced fluid-structure-interaction model problem. We consider higher-order space-time discretizations for both fluid and struc-ture employing discontinuous approximation spaces in time, as introduced by Johnson et al. [36, 37] and further developed by Hughes and his co-workers, see e.g. [34, 49, 46]. Accordingly, the variational expres-sions are formulated directly over the space-time manifold Q, organized in space-time slabs Qn _{= Ω × In,} where Ω is the underlying spatial domain and In = ]tn, tn+1[ are time intervals corresponding to the se-quence {tn} of increasing discrete time levels tn. The space-time approach is advantageous for a number of reasons. First of all, the simultaneous treatment of both space and time eliminates the need for special time integration methods, as the discretization in time is automatically incorporated in the variational expression. The resulting time stepping scheme is fully implicit and can be formulated up to arbitrarily high order. In addition, since we allow the solution to be discontinuous in time, we are at liberty to change the mesh from one slab to the next, which greatly facilitates the use of adaptive mesh refine-ment techniques to track and capture local solution features that propagate through the spatial domain. Finally, the space-time approach is particularly appealing for the discretization of fluid-structure inter-action problems, because it possesses a natural ability to handle time-dependent geometries, rendering our discretization readily extendable to more complicated fluid-structure-interaction problems involving moving and deforming domains [52, 53].

The contents of this chapter are organized as follows. In section 3.1 we present a discretization of the fluid subproblem, which is discontinuous in both space and time. The finite element formulation for the structure subproblem is elaborated in section 3.2. This formulation is based on a spatially conforming discretization, which is C1_{-continuous in space and again discontinuous in time. Finally, we combine the} two discretizations for fluid and structure in section 3.3 to form an aggregated space-time finite element discretizaton for the coupled fluid-structure-interaction problem.

3.1 Discretization of the fluid subproblem

For the fluid subproblem we consider a space-time discontinuous Galerkin finite element method (space-time DGFEM) using discontinuous approximation spaces in both space and (space-time; cf. [38]. Before we present the variational formulation, we first introduce some necessary notation. We shall adhere mostly to the notational conventions and terminology used in [31, 28].

(26)

3.1.1 Preliminaries

Let Ωf be our open, bounded, truncated fluid domain in R2 _{with boundary ∂Ωf} _{:= Γw}_{∪ Γs}_{∪ Γff}_{. The} space-time domain for the fluid subsystem is then defined as Qf := Ωf× I, where I = ]0, T [. Given an ordered series of time levels 0 = t0 < t1 < ... < tN = T , let us partition the time interval I into N subintervals In = ]tn, tn+1[ such that I =SN −1_n=0 In. Accordingly, the space-time domain Qf is split into a sequence of space-time slabs Qn

f := Ωf× In. The boundary of Qnf consists of the lower time boundary Ωn

f := Ωf×{t = tn}, the upper time boundary Ω n+1

f := Ωf×{t = tn+1} and the lateral (spatial) boundary ∂Ωf× In := Pn

w∪ Psn∪ Pffn, with Pwn := Γw× In, Psn:= Γs× In and Pffn:= Γff× In. On the space-time slab Qn

f, let us introduce the space-time tesselation Th,fn, forming a regular or 1-irregular subdivision of the region Qn

f into disjoint open element domains {K} such that Q n

f =

S

K∈Tn

h K.

The restriction to regular or 1-irregular meshes Tn

h,f implies a maximum of one hanging node per face. We note that due to the discontinuous approximation space across time slab interfaces, the tesselation Tn

h,f can be chosen completely independent from one slab to the next. We shall suppose that the family of subdivisions Tn

h,f is shape-regular and that each element K ∈ Th,fn is related to a master element bK through the bijective mapping TK: bK 7→ K, where bK is either the open unit prism bKP or the open unit hexahedron bKH. The spatial and temporal diameter of the element K are denoted by, respectively, hx,K and hτ,K, and the mesh sizes are given by hx= maxK∈Th,fhx,Kand hτ = maxK∈Th,fhτ,K.

To define a basis for the master element bK, let Pm

p represent a polynomial expansion basis of degree p for the m-dimensional unit simplex bKSm. Furthermore, let us assign to each element a pair of nonnegative integers pK= (px, pτ), denoting the local polynomial degree in respectively spatial and temporal direction. A tensorial basis of order pK = (px, pτ) for bK can then be constructed as QpK( bK) := P

2

px ⊗ P

1 pτ if T−1

K (K) = bKP (prismatic reference element) or QpK( bK) := _P₁ px 2 ⊗ P1 pτ if T −1 K (K) = bKH (hexahedral reference element). Collecting the element orders pK and the element mappings TK for slab Qnf in the arrays pn

Q= {pK: K ∈ Th,fn} and T n

f = {TK: K ∈ Th,fn}, we are now ready to introduce the finite element trial and test spaces

Qn_h,p Qnf, Th,fn, pnQ, T n f := qh∈ L2(Qnf) 3 : qh|K◦ TK∈ h Q_p K( bK) i3 ∀ K ∈ Th,fn , Vn_h,p Qn f, Th,fn, pnQ, T n f := vh∈L2_(Qn f) 3 : vh|_K◦ TK∈hQ_p K( bK) i3 ∀ K ∈ Tn h,f . These spaces span a finite-dimensional subspace of the broken anisotropic Sobolev space [Hs_(Qn

f, Th,fn)]3 of composite index s = {(sx,K, sτ,K) : K ∈ Tn

h,f}, which is formally defined as: Hs _Qn f, Th,fn 3 :=nv∈L2(Qnf) 3 : v|_K∈ [HsK (K)]3 ∀ K ∈ Th,fn o .

Moreover, to accommodate the discrete representation of the transpiration velocity, ǫh, we introduce the space En

h on Psn := Γs× In.

For each element K in the subdivision Tn

h,f, we shall refer to ∂K as the union of the faces of K. Furthermore, for x ∈ ∂K we define by nK(x) the outward unit normal vector to ∂K at x. For each element K ∈ Tn

h,f and any v ∈ [H1_(Qn

f, Th,fn)]3we use v−K to denote the inner trace of v on ∂K, i.e. v −

K= limǫ→0(v − ǫnK) (the trace taken from within K). For x ∈ ∂K \ ∂Qn

f we also define the outer trace v +

K of v on ∂K \ ∂Qnf, i.e. v+_K= limǫ→0(v + ǫnK) relative to K, which at the same time is also the inner trace v−_K′ of v relative to the elements K′ _{that have a nonempty two-dimensional intersection with ∂K \ ∂Q}n

(27)

average of vK at x ∈ ∂K \ ∂Qn

f are then defined by [[vK]] := v + K− v − K and {{vK}} := v + K+ v − K /2. Furthermore, we use vK(t±

n) to denote trace values at tn as approached from above and below, i.e. vK(t±

n) = limǫ→0vK(tn± ǫ), and define the jump [[vK]]tn = vK(t +

n) − vK(t−n). In the following it will always be clear from the context which element K in the subdivision Th,fn the quantities nK, v+K and v−K refer to. For the purpose of notational convenience, we shall therefore suppress the subscript (·)_K and simply write n, v+ _{and v}− _instead.

3.1.2 Space-time DGFEM formulation

The space-time DGFEM approximation of (2.19) can now be formulated as follows: given ǫh ∈ Eh := QN −1 n=0 Enh, find qh∈ Qh,p:= QN −1 n=0 Qnh,p such that FDG(qh, vh) = f (vh, ǫh) ∀ vh∈ Vh,p:= N −1_Y n=0 Vn_h,p, (3.1a)

with FDG: Qh,p× Vh,pand f : Vh,p7→ R defined as FDG(qh, vh) := N −1_X n=0 X K∈Tn h,f "Z K qh· L∗fvh dK + Z ∂K\∂Qn f H _q+ h, q − h, n · v−_h dS + Z ∂K∩Pn w 1 2(An+ Nw) q − h · v − hdS + Z ∂K∩Pn s 1 2(An+ Ns) q − h · v − h dS + Z ∂K∩Pn ff A+_nq−_h · v−_h dS + Z ∂K∩Ωn+1_f q−_h t−n+1 · v−_h(tn+1) dS # − N −1_X n=1 X K∈Tn h,f "Z ∂K∩Ωn f q+_h t−nQn−1 f · v − h t + n dS # , (3.1b) f (vh, ǫh) := − N −1_X n=0 X K∈Tn h,f Z ∂K∩Pn s 1 2(An− Ns) (0, 0, ǫh) T · v−_h dS + X K∈T0 h,f Z ∂K∩Ω0 f q0· v−_h t+₀ dS , (3.1c) where, as before, L∗

f represents the formal adjoint of the differential operator Lf, defined by L∗f := −I∂t− ATx∂x− ATy∂y.

We observe that continuity at inter-element boundaries is weakly enforced using numerical fluxes. Firstly, at space-time slab interfaces a simple upwind flux is applied, which arises upon element-wise integration-by-parts and adding the jump term

Z Ωn f v−_h t+n · [[q]]_t_n dS = X K∈Tn h,f Z ∂K∩Ωn f v−_h t+n ·q−_h t+n − qh t−nQn−1 f dS ;

see the last two integrals in (3.1b). This term provides a mechanism by which information is propagated from one space-time slab to the next. Accordingly, the solution qh(t−n) at tn on space-time slab Qn−1f is weakly imposed as an initial condition on the lower time boundary Ωn

f of space-time slab Qnf (note that at n = 0 we simply take qh(t−n) = q0; see the first integral in (3.1c)). As a result of this, the system

(28)

of equations (3.1) decouples and can be solved sequentially in a slab-by-slab fashion. We note that the integral over Ωn

f is evaluated at the common refinement of the interface between T n−1

h,f and Th,fn. Secondly, at interior faces ∂K \ ∂Qn

f inter-element continuity is weakly enforced using the numerical flux function H (·, ·, ·). For DG methods it is common practice to consider monotone flux functions of the form

H _q+_{, q}−_{, n}_{= An}_{(n) {{qh}} −}1

2M[[qh]] . (3.2)

It is well known, see e.g. [6], that the jump [[qh]] has a stability enhancing effect, which is controlled by the upwind matrix M . For M = 0 we simply retrieve a central discretization with zero numerical dissipation. This discretization, however, suffers from oscillations around discontinuities and is generally unstable. The DG literature provides a variety of choices for the flux function H (·, ·, ·). In principle, any two-point monotone Lipschitz function, which is consistent, i.e. H (q, q, n)|_∂K∈FI

n = Anq|∂K∈F I n ∀K ∈ Tn

h,f, and conservative, i.e. H (q+, q−, n)|∂K∈FI

n = − H (q

−_{, q}+_{, −n)|} ∂K∈FI

n, can be used as the numerical flux. A flux function that is often considered in the context of DG methods is for example the (local) Lax-Friedrichs flux, defined by

H_LF _q+_{, q}−_{, n}_{= An}_{(n) {{qh}} −}1

2ρ (An(n)) [[q]] , (3.3) where ρ (An(n)) represents the spectral radius of the Jacobian matrix An(n). Another commonly encountered flux function is the Vijayasundaram flux:

H_V _q+_{, q}−_{, n} _{= A}+ n(n) q −_{+ A}− n(n) q + = 1 2 A + n(n) q −_{+ A}− n(n) q +₋1 2|An(n)| [[q]] ; see also equation (2.9). Here, A+n(n) and A

−

n(n) are defined as in section 2.2.2, and M = |An(n)| = A+n(n)−A

−

n(n). Note that in the linear case, the Vijayasundaram flux – as opposed to the Lax-Friedrichs flux – provides an exact solution to the one-dimensional Riemann problem, cf. Laney [41, section 5.2]. It is therefore equivalent to any of the other well-known approximate Riemann solvers, such as the Godunov, Roe or Engquist-Osher flux. In the sequel, we shall therefore consider the Vijayasundaram flux unless otherwise specified.

At spatial domain boundaries, the numerical flux function H (·, ·, ·) is replaced by appropriate boundary operators derived from the incomplete Riemann problem at the boundary; see section 2.2.2. This results in a weak enforcement of the boundary conditions. Here, Nwand Nsrepresent the boundary operators at, respectively, the fixed wall boundary, Pn

w := Γw×In, and the fluid-structure interface Psn := Γs×In, as defined previously in section 2.2.2. The transpiration condition ǫh∈ Ehat Pn

s is a function of the discrete structural displacement wh and follows from (2.27a). Finally, at the farfield boundary Pn

ff = Γff × In first-order non-reflecting boundary conditions are imposed, cf. (2.17). We remark that the corresponding boundary operator automatically accounts for different conditions at subsonic/supersonic inflow/outflow boundaries.

3.2 Discretization of the structure subproblem

For the discretization of the structure subproblem we consider a space-C1_{-continuous, time-discontinuous} Galerkin finite element method (space-time C1G-DGFEM). Before we present the variational formulation, let us again first introduce and redefine some of the necessary notation.

Goal-Oriented A-Posteriori Error Estimation and Adaptivity for Fluid-Structure Interaction

000000000

000000000

000000000

000000000

000000000

111111111

111111111

111111111

111111111

111111111

REPORT

Goal-Oriented A-Posteriori Error

Estimation and Adaptivity for

Fluid-Structure Interaction

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

Abstract

Contents

Chapter 1

Introduction

1.1

Motivation

1.2

Research objectives

1.3

Thesis outline

Chapter 2

Problem Statement

2.1

The two-dimensional panel problem

2.2

Governing equations for the fluid

2.2.1

The linearized small-disturbance potential flow equation

2.2.2

Initial and boundary conditions

2.2.3

Functional setting

2.2.4

Variational space-time formulation

2.3

Governing equations for the structure

2.3.1

The Euler-Bernoulli beam equation

2.3.2

Initial and boundary conditions

2.3.3

Functional setting

2.3.4

Variational space-time formulation

2.4

Interface conditions

2.5

Aggregated variational space-time formulation

2.6

Target quantities of interest

Chapter 3

Space-Time Finite Element

Discretization

3.1

Discretization of the fluid subproblem

3.1.1

Preliminaries

3.1.2

Space-time DGFEM formulation

3.2

Discretization of the structure subproblem