OPTIMAL PRESSURE BOUNDARY CONTROL OF STEADY MULTISCALE FLUID-STRUCTURE INTERACTION SHELL MODEL DERIVED FROM KOITER EQUATIONS

. The ﬂuid-structure interaction (FSI) problem has been extensively studied, and many papers and books are available in the literature on the subject. In this work, we consider some optimal FSI pressure boundary control applications by using a membrane model derived from the Koiter shell equations where the thickness of the solid wall can be neglected and the computational cost of the numerical problem reduced. We study the inverse problem with the aim of achieving a certain objective by changing some design parameters (e.g. forces, boundary conditions or geometrical domain shapes) by using an optimal control approach based on Lagrange multipliers and adjoint variables. In particular, a pressure boundary optimal control is presented in this work. The optimality system is derived from the ﬁrst-order optimality condition by taking the Fr´echet derivatives of the Lagrangian with respect to all the variables involved. This system is solved by using a ﬁnite element code with mesh-moving capabilities. In order to support the proposed approach, we perform numerical tests where the pressure on a ﬂuid domain boundary controls the displacement that occurs in a well-deﬁned region of the solid domain.


Introduction
Recently, the numerical simulations of fluid-structure interaction (FSI) problems have become more and more popular, and many papers and books have been published on this topic (see [1,2,3,4,5]). The applications of the numerical modeling of FSI systems are various, raging from wind turbines and aircraft to hemodynamics. In FSI problems the fluid flow changes the tensional state of a solid structure that is left free to move and the solid deformation has an important effect on the fluid flow.
Several techniques have been developed to reduce the computational cost of FSI problems. In this respect, this work is based on the reduction of the dimensionality of the solid, through a model built on the Koiter shell equations [6]. In order to couple the fluid and the structure domains, the Koiter shell equations are embedded into the fluid equations as a Robin boundary condition [7]. The coupling fluid-structure conditions are automatically treated implicitly, so the stability of this numerical scheme is preserved. This model has many applications in cases where a fluid interacts with a thin membrane that deforms mainly in the normal direction.
In the last years optimization problems have gained popularity among the research community by using gradient-based adjoint methods, see [8,9]. For adjoint FSI optimization the interested reader can see [10,11,12,13] and the references therein. In this paper, we solve a stationary displacement matching problem where the control variable is the fluid boundary pressure using the Lagrangian multipliers method to obtain the optimality system. The rest of this paper is organized as follows. In Section 2 we introduce the mathematical model describing our multi-scale FSI problem and in Section 3 we derive the optimality system arising from the minimization of the augmented Lagrangian is presented. In Section 4 some simple two-dimensional numerical results are then reported.

Physical model
In this section we introduce the mathematical model for the FSI problem. We consider the classic Navier-Stokes equations to model the fluid motion, and a shell model to describe the solid behavior. In particular, the structural model is based on the Koiter shell approach that considers the model of an elastic thin membrane. We introduce now some functional spaces defined on the domain Ω used in the rest of the paper: we denote with L 2 (Ω) the space of square integrable functions, and with H s (Ω) the standard Sobolev space with norm · s . Moreover, we denote with H s 0 (Ω) the space of all functions in H s (Ω) that vanish on the boundary of Ω.
The Koiter shell approach relies on the assumptions that the structure displacements are small and normal to the shell surface. The domain of the structure is denoted by Γ s , the displacement and the external surface force vectors by η and f s , respectively [7,14]. The weak form of the considered shell equation results for appropriate test functions ψ belonging to a functional space to be determined on the basis of the imposed boundary conditions. Furthermore, ρ s and h s are the density and the thickness of the shell, and E αβλδ and γ αβ are the elasticity and the change of metric, respectively.
In this work negligible bend, shear stresses, and linear elastic constitutive law with a homogeneous and isotropic material are considered [6]. Under these hypotheses the structure model (1) reduces to a simple scalar equation. The dimension of the structure is then reduced by one. So the following simplified model is obtained where ρ f and u are the density and the velocity vector of the fluid, respectively, and η 3 represents the displacement normal to the reference solid surface. In particular, when one considers cylindrical geometries of radius R, it can be demonstrated that [7]  The fluid is modeled as Newtonian, homogeneous and incompressible, described in ALE form as [2,15] where ρ f and u are the density and the velocity vector of the fluid. Moreover, σ f is the Cauchy stress tensor of the fluid written as σ f = −pI + µ(∇u + ∇u T ), µ and p are the dynamic viscosity and the pressure of the fluid, respectively. The system of equations (4) is completed with appropriate boundary conditions. The fluid domain is Ω f , and w is the ALE velocity that determines step by step the position of the nodes of the fluid domain as x f (t) = x 0 + t 0 wdτ . This shell model allows us to reduce by one the dimension of the solid, so the structural equations can be reduced to a boundary condition on Γ s for the fluid problem. The two sub-systems (4) and (2) can be coupled by imposing where Γ D,f are the boundaries of Ω f where a Dirichlet condition is imposed. In order to satisfy the continuity of the test functions φ · n = ψ over the interface surface Γ s in the coupled system, we introduce the following functional space Now we can derive the weak form of the coupled final system for all (φ, ψ) ∈ W 0 , q ∈ L 2 (Ω f ). A finite element technique is used to obtain the discrete weak formulation of (6). Following the work in [7], we treat explicitly the position of the fluid domain, and consider an implicit discretization of the coupling conditions. With this approach, the structural equation can be incorporated into the fluid equations as a boundary condition (Robin scheme).
However, this work is based on the stationary solution of the presented system. Under this hypothesis, the system (6) becomes 3 Optimality system In this work, we are interested in solving a given shell deformation by controlling the fluid pressure over a boundary. For this purpose, we start introducing the following objective functional where the first term is the distance in norm between the actual displacement and the desired value over the controlled boundary Γ d , and the second term is a standard Tychonov regularization term that limits the L 2 -norm of the fluid boundary pressure p c , i.e. the control variable. The regularization parameter λ weights the importance of the two terms over the cost functional. In general, too much regularization leads to smoother but less effective controls, while a lack of regularization may cause numerical issues since usually the optimal solution lies in distributional spaces.
We now introduce the following augmented Lagrangian functional L, that is obtained by adding to the objective functional J the FSI state equations (7) multiplied by a set of Lagrange multipliers. In the following, we refer to the Lagrange multipliers as adjoint variables.
where we integrate by parts the contributions of the fluid stress tensor σ f . The surface integrals can be rewritten by substituting the definition of τ n . The stationary points of the Lagrangian functional can be found by setting to zero the Fréchet derivatives taken for all the problem variables. When the derivatives are taken with respect to the adjoint variables the weak form of the state system (7) is recovered as well as the boundary conditions. By taking the derivatives in the direction δp we get By considering the volume integral we get the following continuity equation for the adjoint velocity while, with the surface contributions, we recover the control equation over the controlled boundary Γ c and the boundary conditions on Γ D On Γ N we have δp = 0 since we prescribe Neumann boundary conditions with fixed pressure.
For δη we have Recalling that Γ d ⊂ Γ t we obtain the following boundary conditions for the adjoint system We collect δu terms and integrate by parts to obtain The strong form of the adjoint velocity reads with boundary conditions Moreover we have to consider the contribution on L given by the motion of the boundary Γ t along the direction δη where χ represents the shell curvature. Under the hypothesis of small deformation we can safely neglect the terms where χ appears. We also have that DL DΓ δη = 0 since the term with u a is defined on the surface Γ t , and a constant extension of it towards the normal direction to the surface leads to a null normal gradient of this term.
In short the optimality system consists of the state system (7), the control equation (12), the adjoint system (11)-(16) and the boundary conditions (13)-(15)-(18). Since the optimality system doubles the state variables the use of a one-shot method is not appropriate and we use a segregated approach for the solution of the state, adjoint, and gradient equations. An advantage of our monolithic aproach is that we can reuse the same solver for both the solution of the state (7) and adjoint systems (11)-(16) with few modifications.

Numerical results
In this section, we report some numerical results obtained by using the mathematical model shown in the previous sections. We consider a rectangular domain Ω = {(x, y) : x ∈ [0, 0.1], y ∈ [0, 0.3]} as shown in Figure 1 on the left. The fluid has density ρ f = 1000 kg/m 3 and dynamic viscosity µ = 100 P a · s. For the solid, we consider β = 60kP a/m and thickness h s = 0.0075 m. For all the presented simulations, the domain was uniformly divided with a regular rectangular mesh.
We implement a standard steepest descent algorithm in the multigrid finite element code FEMuS, that relies on PETSc libraries for the solution of the multigrid discretized linear solver with MPI paralellization.

Plane channel test
In this first test, the fluid flows vertically from the bottom to the top. The region of the boundary Γ 2 represents a solid wall with no-slip boundary condition (u = 0) and Γ 3 is the membrane where we impose the generalized Robin boundary condition (2). In Figure 1, we also report the results obtained simulating the system without control and in steady-state. The displacement of the points in the domain (center) and the pressure (right) are shown. The pressure presents a linearly decreasing trend from the bottom, where p = 6000 P a is imposed, to the top, where we fixed p = 0 P a, and it is shown with lines of iso-magnitude.
x d1 x d2 We solve multiple simulations with different regularization parameter λ. The results are presented in Table 1. Note that the smaller is λ, the closer the displacement of the controlled point x d is to the desired one. This result is expected, since with larger λ the contribution of the regularization term in the minimization of the functional is more relevant. This is introduced in the functional in order to limit the control parameter p to the space of square integrable functions. Therefore, with larger λ we find more regular optimization parameter p, but less precise displacement field η. In Table, we also report the number of iterations needed for the implemented algorithm to find the optimal solution.
We focus now on the controlled inlet pressure field updated through the formula In fact, depending on the regularization parameter, different inlet pressure fields can be obtained. In Figure 2, the controlled pressure field along the boundary Γ 1 is reported for various values of λ. Note that the choice of the regularization parameter strongly affects the controlled pressure field. With weak regularization, the objective term dominates in the functional and the pressure attains large values, thus effectively controlling the membrane displacement. In Figure 2 it is also reported the reference starting pressure. The comparison between the uncontrolled field and all the controlled pressure fields show that the control strongly affects the solution on Γ c = Γ 1 in order to obtain the desired displacement η d . Moreover, in Figure 2

Airbag test
We consider the same geometry of the previous test, togheter with airbag-like boundary conditions: we impose a no-slip condition on Γ 1 ∪ Γ 4 and the Koiter boundary condition on Γ 3 . The controlled pressure is imposed on Γ c = Γ 2 , where we initially impose p = 300 P a. We also consider all the physical properties introduced above. In this framework we want to control the displacement field on the two different points, x d1 and x d2 (see Figure 1, left).
Observing the results obtained in a steady state without control, shown in Figure 3, one can see that the problem is symmetric since we are neglecting the buoyancy forces. The two points x d 1 and x d 2 present a displacement η = 0.001152 m and the pressure in the channel is almost constant everywhere. The goal is to control the displacement of the two points x d 1 and x d 2 , optimizing the pressure of the fluid on the boundary Γ 2 . In order to break the symmetry of the problem, we choose a desired displacement of point x d 1 equal to η d 1 = −0.01 m and a desired displacement of x d 2 equal to η d 2 = 0.01 m. The functional associated with the presented numerical problem reads Many simulations have been done for different regularization parameter λ. In Figure 4, the control parameter p obtained with different values of λ is reported. Note that the trend obtained with λ = 10 −10 shows slightly different values compared to the other two reported cases.
The reason of that is explained in  2.54909 · 10 −5 −0.010032 0.004951 4 η d 2 from a lower value of it. On the contrary, if we consider the results obtained with λ = 10 −10 , we see that we are getting to the objective from a larger displacement. In fact, the functional has a different local minimum and the algorithm in this case moved to a different one. In Figure 4, we also reported the pressure of the fluid in Ω for λ = 10 −10 . Again, note that the value of λ affects the pressure field to be imposed on the control domain Γ 2 . In this case, the resulting velocity implies an inlet of fluid from a region of Γ 2 close to the point x d 2 of which we want to increase the displacement, and an outlet of the fluid in a region of Γ 2 close to the point x d 1 , of which we want to reduce the displacement. In general, this case is an ill-posed control problem, due to various local minima close to each other.

Conclusions
In this work a mathematical and numerical method have been proposed to solve an optimal boundary control problem applied to a fluid-structure interaction model based on the Koiter's equation. Using this formulation, the equations for the solid become boundary conditions for the fluid equations reducing the computational cost of the numerical simulation. Then, we have obtained the optimality system applied to the pressure boundary control problem applied to the introduced Koiter FSI model to find the required