Optimization with nonstationary, nonlinear monolithic fluid-structure interaction

Within this work, we consider optimization settings for nonlinear, nonstationary fluid-structure interaction. The problem is formulated in a monolithic fashion using the arbitrary Lagrangian-Eulerian framework to set-up the fluid-structure forward problem. In the optimization approach, either optimal control or parameter estimation problems are treated. In the latter, the stiffness of the solid is estimated from given reference values. In the numerical solution, the optimization problem is solved with a gradient-based solution algorithm. The nonlinear subproblems of the FSI forward problem are solved with a Newton method including line search. Specifically, we will formally provide the backward-in-time running adjoint state used for gradient computations. Our algorithmic developments are demonstrated with some numerical examples as for instance extensions of the well-known fluid-structure benchmark settings and a flapping membrane test in a channel flow with elastic walls.


Introduction
This paper is devoted to the study of optimal control and parameter estimation problems of nonstationary, nonlinear fluid-structure interaction (FSI). For general overviews on the FSI forward problem, we refer to the books [10,26,29,9,4,7,52,28]. Fluid-structure interaction is still one of the most challenging problem settings within multiphysics applications. The main reason being that the dynamics of both subproblems are exchanged on the interface and accurate discretizations are therefore necessary. Secondly, numerical algorithms are sensitive in terms of stability to the physical parameters; known as added-mass effect [11,27,2,56]. As coupling strategy we choose, in this work, the well-known arbitrary Lagrangian-Eulerian (ALE) technique [17,42,25].
Employing FSI as forward problem within an optimization framework contains the previously mentioned difficulties and yields significant further challenges when dealing with nonstationary problem settings. Historically, this subject falls into the category of PDE-constrained optimization [46]. Studies concentrating on theoretical and computational aspects for stationary FSI optimization are Our algorithmic techniques are substantiated with several numerical tests in Section 6. We summarize our main findings in Section 7.
2 Modeling the FSI forward problem 2

.1 Notation
We denote by Ω := Ω(t) ⊂ R d , d = 2, the domain of the FSI problem. The domain consists of two time-dependent subdomains Ω f (t) and Ω s (t). The FSI-interface between Ω f (t) and Ω s (t) is denoted by Γ i (t) = ∂Ω f (t) ∩ ∂Ω s (t). The initial (or later reference) domains are denoted by Ω, Ω f and Ω s , respectively, with the interface Γ i = ∂ Ω f ∩ ∂ Ω s . Furthermore, we denote the outer boundary by ∂ Ω = Γ = Γ in ∪ Γ D ∪ Γ out where Γ D and Γ in are Dirichlet boundaries (for the velocities and displacements) and Γ out denotes a fluid outflow Neumann boundary, respectively. The displacements are set to zero on Γ out .
As frequently used in the literature, we denote the L 2 scalar product in Ω with (a, b) := (a, b) Ω :=

Spaces
For the function spaces in the (fixed) reference domains Ω, Ω f , Ω s , we define spaces for spatial discretization only. Rather than employing Bochner-spaces [14,61] for space-time functions, the time t is later explicitly accounted for, e.g., [21]  Next, in the fluid domain, we define further: In the solid domain, we use For the FSI problem using variational-monolithic coupling [40,19,20] the velocity spaces are extended from Ω f and Ω s to the entire domain Ω such that we can work with global H 1 functions.
Thus, we define: By this choice, the kinematic and dynamic coupling conditions are automatically satisfied in a variational sense.
Finally, we notice that the spaces on the current domains Ω, Ω f , Ω s are defined correspondingly without 'hat' notation.

2.3
The ALE concept, transformed fluid flow, and solids in Lagrangian coordinates In this section, we recapitulate the ingredients to formulate a coupled problem (i.e., fluid-structure interaction) with the help of the ALE approach. The ALE mapping A : Ω f → Ω f is defined first.

The ALE transformation and ALE time-derivative
First, we define the ALE transformation: The ALE mapping is defined in terms of the vector-valued (artificial) fluid mesh which is specified through the deformation gradient and its determinant Furthermore, function values in Eulerian and Lagrangian coordinates are identified by Here, I denotes the identity matrix. The mesh velocity is defined by w := ∂ t A.
The key quantity to measure the fluid mesh regularity is J. The artificial fluid displacement u f (the mesh motion) is obtained in this work by solving a biharmonic equation [35,58,19,57].
Finally, the transformation between different coordinate systems requires transformation of derivatives. For a vector-valued function u ∈ Ω and u ∈ Ω it holds, e.g., [39]: Finally, the ALE time-derivative is the total derivative of an Eulerian field and is important when working on moving domains:

Equations for fluids and solids
In this section, we briefly state the basic underlying equations first separately. In the following, we first present fluid flow and then the solid part.

Strong forms
The isothermal, incompressible Navier-Stokes equations in an ALE setting read: Given v in , h f , and where the (symmetric) Cauchy stress is given by with the density ρ f and the kinematic viscosity ν f . Later in the FSI problem, the function h f will be given by the solid velocity v s . The normal vector is denoted by n f .
The equations for geometrically non-linear elastodynamics are given as follows: Given h s , u 0 , and v 0 ; find u s : Ω s × I → R d such that The constitutive law is given by the tensor: Here, µ and λ are the Lamé coefficients for the solid. The solid density is denoted by ρ s and the solid deformation gradient is F = I + ∇ u s . Later in FSI, the vector-valued function h s will be given by the normal stress from the fluid problem. Furthermore, n s denotes the normal vector.

Variational forms
The previous Navier-Stokes equations in a variational ALE framework described in a reference domain Ω f are given by: Let v D f a suitable extension of Dirichlet inflow data. Find vector-valued velocities and a scalar-valued pressure f are satisfied, and for almost all times t ∈ I holds: Here, g f := − ρ f ν f F −T ∇ v T f denotes a correction term on the outflow boundary and n f is the outer normal vector. The transformed Cauchy stress tensor reads: The variational formulation for elastodynamics can be formulated as a first-order-in-time system: Formulation 2.3 (First order system in time weak formulation of elasticity including strong damping). Find u s ∈ V 0 s and v s ∈ L s with the initial data u s (0) = u 0 and v s (0) = v 0 such that for almost all times t ∈ I: Dirichlet conditions and the second condition is a (non-homogeneous) Neumann condition.
In variational-monolithic coupling these Dirichlet conditions are built into the corresponding function space by employing a globalized Sobolev space V 0 (see (1)). Neumann type conditions are weakly incorporated through interface integrals, which actually cancel out in the later models because of their weak continuity thanks to working with the space V 0 .
For the fluid problem, continuity of velocities is required (i.e., a kinematic coupling condition): To complete the solid problem, we must enforce the balance of the normal stresses on the interface (i.e., a dynamic coupling condition): For the geometric problem we have from which we obtain immediately ∂ t u s = v s = v f on Γ i by temporal differentiation.

The FSI model using biharmonic mesh motion
Combining the previous equations for fluids and solids and applying biharmonic mesh motion for realizing the ALE mapping, we obtain the following FSI model [19,57,58]: . Let the constitutive laws from before be given and α > 0 be a small parameter. Find a global vector-valued velocity, vector-valued displacements, additional displacements (due to the splitting of the biharmonic mesh motion model into two secondorder equations) and a scalar-valued fluid pressure, i.e., and u s (0) = u 0 s are satisfied, and for almost all times t ∈ I holds: The Neumann coupling conditions on Γ i are fulfilled in a variational way and cancel in monolithic modeling due to the global test space V 0 in which the test functions from both the fluid and the solid subdomains coincide on the interface. Thus, the condition is implicitly contained in the above system.

Discretization
In this section, we discuss temporal and spatial discretization of the forward problem. Our derivation contains many details on all terms of the FSI forward problem. The overall problem can be posed, however, in an abstract fashion, which facilitates the derivation of the backward-in-time adjoint problem in Section 4.

Temporal discretization
Our goal is to apply A-stable finite differences in time. Specifically, time discretization is based on a One-step-θ scheme as presented for the pure FSI problem, Formulation 2.4, in [58].
In more detail, semi-discretization in time yields a sequence of generalized steady-state problems that are completed by appropriate boundary values at every time step. Let Time derivatives are discretized with a backward difference quotient such that where u := u n := u(t n ), v := v n := v(t n ), u n−1 := u(t n−1 ), v n−1 := v(t n−1 ). Furthermore, the mesh Formulation 3.1 (The time-discretized abstract problem). We aim to find where the semi-linear form A(·)(·) is split into Details of this decomposition are provided in Definition 3.2. all 'standard' terms (e.g., stress terms, fluid convection). We then obtain the decomposition: where the fluid stress tensor σ f is further split into σ f,vu , σ f,p : The (nonlinear) time derivative in A T ( U )( Ψ) is approximated by a backward difference quotient. For the time step t n ∈ I, for n = 1, 2, . . . , N ( N ∈ N), we compute v := v n , u i := u n i (i = f, s) via where we introduce the parameter θ ∈ [0, 1]. Furthermore, we use and u n i := u i (t n ), v n := v(t n ), and J := J n := J(t n ). In our computations in Section 6, we always consider J n,0.5 . The former time step is given by v n−1 , etc. for i = f, s. , w n−1 , p n−1 f } and the time step k := k n = t n − t n−1 be given. In order to solve (12), we seek U n = { v n , u n f , u n s , w n , p n f } ∈ X 0 by employing one-step-θ splitting: The concrete scheme depends on the choice of the parameter θ ∈ [0, 1]. For θ = 1 we obtain the strongly A-stable backward Euler scheme (BE). If k < 0.5, for θ = 0.5 + k, we obtain the 2nd order (shown for linear parabolic problems in [50,36]), A-stable, globally stabilized, Crank-Nicolson scheme (CNs).

Spatial discretization
The time discretized formulation is the starting point for the Galerkin discretization in space. To this end, we construct a finite dimensional subspace X 0 h ⊂ X 0 to find an approximate solution to the continuous problem. As previously explained, in the context of our variational-monolithic formulation, the computations are done on the reference configuration Ω. We use two dimensional shape-regular meshes. A mesh consists of quadrilateral cells K. They perform a non-overlapping cover of the computation domain Ω ⊂ R d , d = 2. The corresponding mesh is given by T h = { K}. The discretization parameter in the reference configuration is denoted by h and is a cell-wise constant that is given by the diameter h K of the cell K.
On T h , the conforming finite element space for { v h , u f,h , u s,h , p f,h , w h } is denoted by the space X h ⊂ X. For Navier-Stokes flow, i.e., { v h , p f,h }, we prefer the biquadratic, discontinuous-linear Q c 2 /P dc 1 element. For the specific definitions of the single elements, we refer the reader to [13]. The property of the Q c 2 /P dc 1 element is continuity of the velocity values across different mesh cells [31]. However, the pressure is defined by discontinuous test functions. Therefore, this element preserves local mass conservation, is of low order, gains the inf-sup stability, and is an optimal choice for both fluid problems and fluid-structure interaction problems. The two displacement variables, namely u h , w h are discretized with Q c 2 elements. In total, the discretized forward problem consists of This abstract formulation serves as basis to derive the adjoint state in Section 4.

Gradient computation
We are interested in identifying material parameters, e.g., µ in (6). To this end, we denote by q ∈ R p , with p ≥ 1, the collection of these parameters, and will define suitable cost functionals J (q, U ) to be minimized. To highlight the dependence of the equation on the parameters q, we add an additional q argument to the form A E , e.g., we consider (15). Assuming that (15) admits a unique solution for any given q, we can obtain an optimality system by the standard Lagrange formalism, see, e.g., [46,55,38]. For a rigorous proof of the required differentiability properties, some progress has been made for stationary FSI-problems in [60]. A rigorous derivation of the corresponding adjoints in the context of shape optimization can be found in [34].
The formal Lagrange approach provides an adjoint equation to (15) as Here ∂ U n A denotes the directional derivative of the form A with respect to its U n argument, and the first argument of the second parentheses denotes the respective direction.
With this adjoint, we obtain the total derivative of the cost functional q → J (q) := J (q, U ) in a direction δq as allowing the calculation of the reduced gradient ∇J(q) ∈ R p of the cost functional by cf., e.g., [5].

Update
In the first step of the previous algorithm for solving the nonlinear primal problem, we employ the following Newton method. At each time point the following problem is given: Algorithm 2 (Residual-based Newton's method). We omit h and n for the convenience of the reader.
This choice is justified since the spatial numbers of degrees of freedom is moderate in our numerical examples. Moreover, in order to save computational cost, we adopt simplified Newton steps; i.e., the matrix A ( U k )(δ U k , Ψ) is only rebuild when λ l k < 1 (defined below) or A( U k ) ∈ [0.001, 1] A( U k−1 ) .

The criterion for convergence is the contraction of the residuals:
3. If (20) (20) is fulfilled for a l * < l M or l M is reached. In the latter case, no convergence is obtained and the program aborts.
4. In case of l * < l M we check next the (relative) stopping criterion: If this is criterion is fulfilled, set U n := U k+1 . Otherwise, we increment k → k + 1 and goto Step 1.

Numerical tests
We conduct several numerical tests in this section. These are implemented in the open-source package DOpElib [18,32] using the finite elements of deal.II [1].

Material parameters
The material parameters for the FSI-1 and FSI-3 tests are chosen as proposed in [41,9] and listed in Table 1. The parameters for the flapping example are a mixture of [30,59,24].

Example 1: Parameter estimation within the FSI 1 benchmark
In this first numerical example, we consider a quasi-stationary setting based on the FSI 1 benchmark [41,9]. The optimization problem reads: identify the Lamé parameter µ s from measurements of the beam-tip displacement at A := (0.6, 0.2). The exact values u d are taken from a reference solution with µ = 0.5 · 10 6 kg m s 2 . The forward problem is solved with the backward Euler scheme, i.e., θ = 1, since the configuration is stationary and we only use a time-dependent method to find the stationary limit.

Configuration
The geometry of the FSI-1 and FSI-3 settings are displayed in Figure 1. An elastic beam is attached to a cylinder and is surrounded by an incompressible fluid. The initial geometry is once uniformly refined in space. On the cylinder and outer boundary Γ w we enforce zero Dirichlet boundary conditions for v and u.
On the outflow boundary Γ out we prescribe the do-nothing outflow condition [37]. The inflow profile on Γ in is given by: v(0, y) := 1.5 y (0.41 − y) 4 0.41 2 v mean (t). The mean inflow v mean (t) is 0.2 m/s for Example 1 (FSI 1) and 2.0 m/s in Example 2 (FSI 3). In the FSI 1 test case, we compute n = 25 time steps using k = 1 s and in the FSI 3 example, we work with k = 0.001 s with T = 0.6 s corresponding to n = 6000 time steps.

Discussion of the FSI 1 findings
Our results for three different configurations are displayed in the Tables 2, 3, and 4. In the first run with α = 0.001, the algorithm converges slowly in order to estimate q k and to reduce the cost functional J (q k ). The main reason is due to the low regularization, which is confirmed by two further runs with α = 0.1 and 1.
In Table 3, the value of α is enlarged to 0.1. Here, in 155 gradient iterations, the cost functional is reduced by a order to 10 14 from an initial control q 0 = 5000 to q 155 = 10 6 .
Increasing further α to 1 (Table 3) yields a reduction in J (q k ) from about 10 11 to 10 −6 . The gradient algorithm converges in 5 iterations.

Example 2: Optimal control within the FSI 3 benchmark
In this second numerical test, we employ the same geometry as in Example 1. The material parameters and boundary data can be found in Table 1 and Section 6.2.2. We now consider an optimal control problem in which µ s is detected in such a way to match the displacement value at the beam tip at (0.6, 0.4) obtained by the FSI 1 simulation in Example 1. Since this numerical test is nonstationary with periodic solutions in the original forward run, we use the shifted Crank-Nicolson time-stepping scheme with minimal numerical dissipation.

Cost functional
The cost functional is given by:  Table 3: Optimization results for the FSI 1 example with α = 0.1 and q d = 10 6 . The initial Residual in q 0 = 5 000 is |∇J (q 0 )| = 9.875 · 10 4 with q d = 2.27007 · 10 −5 , i.e., the functional is similar to Example 1 and the reference value comes from FSI 1, but not FSI 3.

Discussion of the FSI 3 findings
Graphical plots of the solution are provided in Figure 2. Our quantitative results are shown in Table 5.

Example 3: Two-dimensional flapping membranes
In this third example, we consider two-dimensional flap dynamics. This test is a challenge because of the thin flaps and the mesh regularity. The original setups for forward simulations were inspired by [30]. Our configuration here is a further extension, towards FSI-optimization, of [59] and [24].

Cost functional
The cost functional is given by: where T is the end time value as in the other examples and F (·) is the drag functional defined as where n is the unit normal vector pointing outward of the domain Ω s and e 1 the first unit vector in R 2 . The boundary part, where the drag is evaluated is Moreover, we notice that we only control µ in the valves, while in the rest of the solid, the value is as in Table 1.

Conclusions
In this work, we developed settings for FSI-based optimization. Therein, the FSI problem is nonlinear