Optimal Pressure Boundary Control of Steady Multiscale Fluid-Structure Interaction Shell Model Derived from Koiter Equations

Abstract

Fluid-structure interaction (FSI) problems are of great interest, due to their applicability in science and engineering. However, the coupling between large fluid domains and small moving solid walls presents numerous numerical difficulties and, in some configurations, where the thickness of the solid wall can be neglected, one can consider membrane models, which are derived from the Koiter shell equations with a reduction of the computational cost of the algorithm. With this assumption, the FSI simulation is reduced to the fluid equations on a moving mesh together with a Robin boundary condition that is imposed on the moving solid surface. In this manuscript, we are interested in the study of inverse FSI problems that aim to achieve an objective by changing some design parameters, such as forces, boundary conditions, or geometrical domain shapes. We study the inverse FSI membrane model by using an optimal control approach that is based on Lagrange multipliers and adjoint variables. In particular, we propose a pressure boundary optimal control with the purpose to control the solid deformation by changing the pressure on a fluid boundary. We report the results of some numerical tests for two-dimensional domains to demonstrate the feasibility and robustness of our method.

1. Introduction

Recently, numerical simulations of fluid–structure interaction (FSI) problems have gained popularity in the research community due to a large variety of possible applications, which range from wind turbines and aircraft to hemodynamics. In FSI models, the fluid changes the tensional state of a solid structure that is left free to move and deforms the fluid flow domain. Many approaches have been proposed and implemented to model the behavior of the interaction between a fluid and solid, and a large variety of articles and books is available on this topic [1,2,3,4,5]. A first attempt to classify the FSI algorithms can be made while considering the kind of coupling between fluid and solids. In partitioned approaches, the equations of fluid mechanics, structural mechanics, and mesh motion are solved sequentially. This allows re-using existing fluid and structure solvers, such as commercial software (e.g., Ansys [6], etc.). However, convergence difficulties are sometimes encountered, most commonly when the added mass effect is not negligible or when an incompressible fluid is fully enclosed by the structure. The interested reader in partitioned approaches can see [7,8,9,10]. In monolithic approaches, the equations of fluid, structure, and mesh motion are solved simultaneously, in a fully-coupled fashion, see [11,12,13]. The main advantage is that these solvers are more robust and many of the problems encountered with the staggered approaches are avoided. However, strongly-coupled approaches necessitate writing a fully-integrated FSI solver, which virtually precludes the use of existing fluid and structure solvers.

It can be easily understood the interest among the scientific community towards methods that try to combine the stability and convergence properties of monolithic approaches, with the reduced computational costs that are required by partitioned ones. In this context, this work is based on the reduction of the dimensionality of the solid, through a model derived from the Koiter shell equations [14,15]. This model has many applications in vascular dynamics, where a fluid interacts with a thin membrane that mainly deforms in the normal direction [16,17]. To couple the fluid and structure domains, the Koiter shell equations are embedded into the fluid equations as a Robin boundary condition  [18,19]. The fluid–structure coupling conditions are automatically treated in an implicit way as for monolithic approaches, so that the stability of this numerical scheme is preserved.

In this paper, we investigate a class of FSI steady inverse control problems modeled by the Koiter shell equations. We use a Lagrange multiplier approach and keep a monolithic formulation in velocity-displacement fields that implicitly takes care of the force balance. We introduce an auxiliary velocity field that keeps the optimality system in the state and adjoint variables symmetric. We propose solving a stationary displacement matching problem with a boundary optimal control approach. We consider a Neumann boundary control, where the control variable is the fluid boundary pressure. In the evaluation of the control, great care must be taken in choosing the appropriate functional spaces. In this case, the inlet pressure is considered in the space of square integrable functions. Although the mathematical background of gradient-based adjoint methods is of great interest to investigate the existence of optimal solutions, see [20,21,22], only a few works that deal with adjoint FSI optimization can be found in literature [23,24,25,26]. In [23], the authors study an inverse parameter estimation problem to recover the arterial stiffness from the measured displacement. Pressure boundary control problems are solved in [24] with the Newton method and in [25] through an auxiliary displacement field that permits extending the velocity field to the solid domain obtaining a symmetric adjoint formulation. This extension method is applied in [26] to a larger class of problems, such as distributed control and inverse Young modulus estimation with control inequalities.

The rest of this paper is organized, as follows. In Section 2, we first introduce the equations used to model the interaction between an incompressible, viscous, Newtonian fluid, and a deformable wall described by the elastic Koiter shell equations. Subsequently, we recover the weak formulation of our Koiter fluid–structure model. In Section 3, we introduce the optimal control problem and the optimality system that arises from the minimization of the augmented Lagrangian. In Section 4, we present some numerical results, proving the effectiveness of the method to find optimal control solutions.

2. Physical Model

In this section, we introduce the mathematical model for FSI problem. While the fluid has been modeled with the classic Navier–Stokes equations, a shell model was used for the solid modeling. In particular, the structural model is based on the Koiter shell approach that considers the model of an elastic thin membrane. In the following, we denote with $L^2(\Omega )$ the space of square integrable functions on the domain $\Omega$, and with $H^s(\Omega )$ the standard Sobolev space with norm ${\parallel ·\parallel }_s$. We also denote, with $H_0^s(\Omega )$, the space of all functions in $H^s(\Omega )$ that vanish on the boundary of $\Omega$, and with $H^{-s}(\Omega )$ the dual space of $H_0^s(\Omega )$. The trace space for the functions in $H^1(\Omega )$ is denoted by $H^{1/2}(\Gamma )$.

2.1. The Linear Koiter Shell Model

The Koiter shell approach relies on the assumption that the structure displacements are small and normal to the surface of the shell. Let $\mathit{\theta }\left(\mathit{x}\right)$ be a linear mapping between a reference fixed point and a generic point on the middle surface. We denote, with $\omega$, the points of the reference configuration and with ${\Gamma }_s$ the deformed membrane configuration, as shown in Figure 1.

We now consider the tangential base ${\tilde{\mathit{a}}}_{\alpha }={\displaystyle \frac{\partial \mathit{\theta }}{\partial x_{\alpha }}}$, for $\alpha =1,2$. Thus, ${\tilde{\mathit{a}}}_1$ and ${\tilde{\mathit{a}}}_2$ define the plane tangential to the shell, while ${\tilde{\mathit{a}}}_n=({\tilde{\mathit{a}}}_1\times {\tilde{\mathit{a}}}_2)/|{\tilde{\mathit{a}}}_1\times {\tilde{\mathit{a}}}_2|$ defines the unit vector normal. Now, we define the covariant metric tensor of the middle surface as $A_{\alpha \beta }={\tilde{\mathit{a}}}_{\alpha }·{\tilde{\mathit{a}}}_{\beta }$. Another fundamental tensor considered in this work is

$$B_{\alpha \beta }={\tilde{\mathit{a}}}_n·{\displaystyle \frac{\partial }{\partial x_{\alpha }}}\left({\displaystyle \frac{\partial \mathit{\theta }}{\partial x_{\beta }}}\right)={\tilde{\mathit{a}}}_n·{\displaystyle \frac{\partial {\mathit{a}}_{\alpha }}{\partial x_{\beta }}}.$$

The deformation can now be specified by the difference between ${\tilde{A}}_{\alpha \beta }$ and ${\tilde{B}}_{\alpha \beta }$ in the deformed and reference configuration. Subsequently, we introduce the strain measures

$${\gamma }_{\alpha \beta }=\frac{1}{2}({\tilde{A}}_{\alpha \beta }-A_{\alpha \beta }),{\rho }_{\alpha \beta }={\tilde{B}}_{\alpha \beta }-B_{\alpha \beta },$$

where ${\gamma }_{\alpha \beta }$ is the change of metric tensor, while ${\rho }_{\alpha \beta }$ denote the covariant components of the linearized change of curvature tensor that is associated with a displacement $\mathit{\eta }={\eta }_i{\mathit{a}}^i$. In the rest of the paper, given a tensor ${\mathcal{T}}_{\alpha \beta }$, we define ${\mathcal{T}}_{\beta }^{\lambda }={\left(A_{\alpha \lambda }\right)}^{-1}{\mathcal{T}}_{\alpha \beta }$. The strain measures can be written in terms of the displacement components $\eta$, as

$${\gamma }_{\alpha \beta }=\frac{1}{2}({\eta }_{\alpha |\beta }+{\eta }_{\beta |\alpha })-B_{\alpha \beta }{\eta }_n,$$

(1)

$${\rho }_{\alpha \beta }={\eta }_{|\alpha \beta }-B_{\alpha }^k{B}_{k\beta }{\eta }_n+B_{\alpha }^k{\eta }_{k|\beta }+B_{\beta }^k{\eta }_{k|\alpha }+B_{\alpha |\beta }^k{\eta }_k,$$

(2)

where the additional subscript preceded by a vertical stroke that is applied to the generic quantity $\chi$, e.g., ${\chi }_{\beta |\alpha }$, means the covariant surface differentiation with respect to the coordinate $x_{\alpha }$.

In this work, we only consider normal displacements, thus there is not any change in length between the reference middle surface and the deformed one. Note that the normals are completely specified by arbitrary variations of $A_{\alpha \beta }$ and $B_{\alpha \beta }$ under the compatibility conditions (1) and (2). We can now introduce the linear Koiter shell equations. Let ${\Gamma }_s$ be the considered shell domain, and let ${\upsilon }_s$ be a measurable subset of $\partial {\Gamma }_s$. In the following, ${\partial }_t$ will denote the outer normal derivative operator along ${\Gamma }_s$. Because ${\gamma }_{\alpha \beta }\left(\mathit{\eta }\right),{\rho }_{\alpha \beta }\left(\mathit{\eta }\right)\in L^2\left({\Gamma }_s\right)$, we can introduce the functional space $\mathit{V}\left({\Gamma }_s\right)$, such that

$$\mathit{V}\left({\Gamma }_s\right)=\{\mathit{\eta }={\eta }_i{\mathit{a}}^i:{\eta }_{\alpha }\in H^1\left({\Gamma }_s\right),{\eta }_n\in H^2\left({\Gamma }_s\right),{\eta }_i={\partial }_t{\eta }_n=0\mathrm{on}{\upsilon }_s\}.$$

For more information regarding the introduced functional spaces and the derivation of the equations used in this section, see [27]. We can define the energy functional of the Koiter shell equation as

$$\begin{array}{cc}\hfill \mathcal{K}\left(\mathit{\eta }\right)& =\frac{1}{2}{\int }_{{\Gamma }_s}\left(\epsilon A^{\alpha \beta \sigma \tau }{\gamma }_{\sigma \tau }\left(\mathit{\eta }\right){\gamma }_{\alpha \beta }\left(\mathit{\eta }\right)+\frac{{\epsilon }^3}{3}{A}^{\alpha \beta \sigma \tau }{\rho }_{\sigma \tau }\left(\mathit{\eta }\right){\rho }_{\alpha \beta }\left(\mathit{\eta }\right)\right)\sqrt{a}dy\hfill \\ \hfill & -{\int }_{{\Gamma }_s}{f}^{i,\epsilon }{\mathit{\eta }}_i\sqrt{a}dy\forall \mathit{\eta }\in \mathit{V}\left({\Gamma }_s\right),\hfill \end{array}$$

(3)

where $a\left(\mathit{\eta }\right)=det\left(A_{\alpha \beta }\left(\mathit{\eta }\right)\right)$ and $\epsilon$ is the thickness of the solid wall. The contravariant components of the shell elasticity tensor $A^{\alpha \beta \sigma \tau }$ are defined as

$$A^{\alpha \beta \sigma \tau }={\displaystyle \frac{2E\nu }{1-{\nu }^2}}{A}^{\alpha \beta }{A}^{\sigma \tau }+{\displaystyle \frac{E}{1+\nu }}(A^{\alpha \sigma }{A}^{\beta \tau }+A^{\alpha \tau }{A}^{\beta \sigma }),$$

where $\nu$ is the Poisson coefficient and E the Young modulus of the solid material. The given functions $f^{i,\epsilon }\in L^2\left({\Gamma }_s\right)$ take the forces applied to the shell into account. We impose the boundary condition ${\eta }_i={\partial }_t{\eta }_n=0$ on ${\upsilon }_s$, therefore the shell is clamped on the boundaries of its middle surface. For all of the admissible displacements ${\mathit{\xi }}^{\epsilon }={\xi }_1^{\epsilon }{a}^i\in \mathit{V}\left({\Gamma }_s\right)$, we can derive the variational equation for the vector field $\mathit{\eta }$, as

$${\int }_{{\Gamma }_s}\left(\epsilon A^{\alpha \beta \sigma \tau }{\gamma }_{\sigma \tau }\left({\mathit{\xi }}^{\epsilon }\right){\gamma }_{\alpha \beta }\left(\mathit{\eta }\right)+\frac{{\epsilon }^3}{3}{A}^{\alpha \beta \sigma \tau }{\rho }_{\sigma \tau }\left({\mathit{\xi }}^{\epsilon }\right){\rho }_{\alpha \beta }\left(\mathit{\eta }\right)\right)\sqrt{a}dy={\int }_{{\Gamma }_s}{f}^{i,\epsilon }{\eta }_i\sqrt{a}dy.$$

(4)

In addition to the hypothesis of normal displacement, we also introduce the following assumptions: small deformations of the solid shell, negligible bending terms, and isotropic and homogeneous material. Therefore, we can neglect the term with ${\epsilon }^3$ in Equation (4). Moreover, we can assume ${\mathit{\xi }}^{\epsilon }$ as a test function for the variational formulation, obtaining

$${\int }_{{\Gamma }_s}\epsilon A^{\alpha \beta \sigma \tau }{\gamma }_{\alpha \beta }\left(\mathit{\eta }\right){\gamma }_{\sigma \tau }\left(\mathit{\psi }\right)d\Gamma ={\int }_{{\Gamma }_s}{\mathit{f}}_s·\mathit{\psi }d\Gamma ,\forall \mathit{\psi }\in \mathit{V}\left({\Gamma }_s\right).$$

(5)

If we restrict the membrane displacements to only normal direction, we can further simplify the variational formulation (5), and reduce it to a simple scalar equation for ${\eta }_n$. In strong form, it reads

$$\begin{array}{cccc}\hfill & \beta {\eta }_n=f_s\hfill & \hfill \mathrm{on}& {\Gamma }_s,\hfill \\ \hfill & \mathrm{with}{\eta }_n{|}_{t=0}={\eta }_0,{\displaystyle \frac{\partial {\eta }_n}{\partial t}}{|}_{t=0}={\eta }_v\hfill & \hfill \mathrm{on}& {\Gamma }_s,\hfill \end{array}$$

(6)

where

$$\beta (x_1,x_2)=\frac{\epsilon E}{1-{\nu }^2}\left((1-\nu )B_{\beta }^k{B}_k^{\beta }+\nu B_{\beta }^{\beta }{B}_k^k\right).$$

(7)

2.2. The Prestressed Term

The presented Koiter model does not account for prestressed loading, along with the shell structure. Let us consider the deformed non-shell configuration ${\Omega }_s$. Note that ${\Omega }_s$ has thickness $\epsilon$, and the shell surface is defined as the middle surface of it. In a weak formulation, the prestress term is [28]

$${\int }_{{\Omega }_s}\nabla \left(\mathit{\eta }\mathit{P}\right):\nabla {\mathit{\xi }}^{\epsilon }d\mathit{x},$$

(8)

where $\mathit{P}$ is the Cauchy stress tensor in the deformed configuration for only tangential stresses in ${\Omega }_s$. In the rest of this section, we will follow the procedure that is presented in [18]. By taking the limit for $\epsilon \to 0$, we lead Equation (8) back to the membrane case. Thus, we can write the deformation field as

$$\mathit{\eta }={\eta }_i(x_1,x_2){\mathit{a}}^i+x_n\left({\displaystyle \frac{\partial {\eta }_n}{\partial x_{\alpha }}}+B_{\alpha }^k{\eta }_k\right){\mathit{a}}^{\alpha }.$$

Because the integration is performed on ${\Omega }_s={\Gamma }_s\times [-\epsilon /2,+\epsilon /2]$, the terms where $x_n$ appears are of higher order in $\epsilon$ (in the order of ${\epsilon }^2$ or more), so we can neglect them. The interested reader can see the full expansion in [18]. Now, we introduce the surface covariant derivative of a vector field (see [29]). The covariant derivative of $\mathit{\eta }$ is defined as

$${\mathit{\eta }}_{\alpha |k}^s={\displaystyle \frac{\partial {\eta }_{\alpha }}{\partial x_k}}-{\Gamma }_{\alpha k}^{\beta }{\eta }_{\beta },\mathrm{with}{\Gamma }_{\alpha k}^{\beta }=a^{\beta }·{\displaystyle \frac{\partial a_k}{\partial x_{\alpha }}}.$$

With this notation, we can write the three-dimensional covariant derivatives of $\mathit{\eta }$ as

$$\begin{array}{cccc}\hfill & {\mathit{\eta }}_{\alpha |\beta }={\mathit{\eta }}_{\alpha |\beta }^s-B_{\alpha \beta }{\eta }_n,\hfill & \hfill & {\mathit{\eta }}_{\alpha |n}=-\left({\displaystyle \frac{\partial {\eta }_n}{\partial x_{\alpha }}}+B_{\alpha }^k{\eta }_k\right),\hfill \\ \hfill & {\mathit{\eta }}_{n|\alpha }={\displaystyle \frac{\partial {\eta }_n}{\partial x_{\alpha }}}+B_{\alpha }^k{\eta }_k,\hfill & \hfill & {\mathit{\eta }}_{n|n}=0.\hfill \end{array}$$

Now by integrating (8) and considering ${\Omega }_s={\Gamma }_s\times [-\epsilon /2,+\epsilon /2]$, we have

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_s}{\mathit{\eta }}_{k|\alpha }{P}^{\alpha \beta }{\mathit{\xi }}_{k|\beta }^{\epsilon }d\mathit{x}={\int }_{{\Gamma }_s}{\int }_{-\epsilon /2}^{\epsilon /2}{B}_{k\alpha }{\eta }_n{P}^{\alpha \beta }{B}_{k\beta }bb{\xi }_n^{\epsilon }dxd\Gamma ={\int }_{{\Gamma }_s}\epsilon B_{k\alpha }{P}^{\alpha \beta }{B}_{k\beta }{\eta }_n{\xi }_n^{\epsilon }d\Gamma ,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & {\int }_{{\Omega }_s}{\mathit{\eta }}_{n|\alpha }{P}^{\alpha \beta }{\mathit{\xi }}_{n|\beta }^{\epsilon }d\mathit{x}={\int }_{{\Gamma }_s}{\int }_{-\epsilon /2}^{\epsilon /2}{\displaystyle \frac{\partial {\eta }_n}{\partial x_{\alpha }}}{P}^{\alpha \beta }{\displaystyle \frac{\partial {\xi }_n^{\epsilon }}{\partial x_{\beta }}}dxd\Gamma ={\int }_{{\Gamma }_s}\epsilon P^{\alpha \beta }{\displaystyle \frac{\partial {\eta }_n}{\partial x_{\alpha }}}{\displaystyle \frac{\partial {\xi }_n^{\epsilon }}{\partial x_{\beta }}}d\Gamma .\hfill \end{array}$$

(10)

Equation (9) can be incorporated into the coefficient $\beta$ in the Equation (7). The Equation (10) gives a second derivative in space to be added inside the model (6), obtaining

$$\begin{array}{cccc}\hfill & {\beta }^{*}{\eta }_n-\nabla ·\left(\mathit{P}\nabla {\eta }_n\right)=f_s\hfill & \hfill \mathrm{on}& {\Gamma }_s,\hfill \\ \hfill & \mathrm{with}{\eta }_n{|}_{t=0}={\eta }_0,{\displaystyle \frac{\partial {\eta }_n}{\partial t}}{|}_{t=0}={\eta }_v\hfill & \hfill \mathrm{on}& {\Gamma }_s,\hfill \end{array}$$

(11)

and

$${\beta }^{*}(x_1,x_2)=\epsilon \left(\frac{E}{1-{\nu }^2}\left((1-\nu )B_{\beta }^k{B}_k^{\beta }+\nu B_{\beta }^{\beta }{B}_k^k\right)+B_{k\alpha }{P}^{\alpha \beta }{B}_{k\beta }\right).$$

Model (11) must be completed with the proper boundary conditions, i.e., ${\eta }_n{|}_{\partial {\Gamma }_s}=0$. Note that the presented prestressed model can be used when the deformed configuration is close enough to the reference one, so we can consider an isotropic elastic tensor.

2.3. The Cylindrical Geometry

Let now consider a cylindrical geometry, in order to show how to explicitly calculate all of the introduced tensors in a simple three-dimensional geometry. If we consider a system of cylindrical coordinates and a cylinder of radius R, we have $r(\theta ,z)=\{(x,y,z)\in {\mathbb{R}}^3|x=Rcos\left(\theta \right),y=Rsin\left(\theta \right),z=z\}$. The covariant basis is given by

$${\mathit{a}}_1=\left[\begin{array}{c}-Rsin\theta \\ Rcos\theta \\ 0\end{array}\right],{\mathit{a}}_2=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],{\mathit{a}}_n=\left[\begin{array}{c}cos\theta \\ sin\theta \\ 0\end{array}\right],$$

this implies that

$$A_{\alpha \beta }=\left[\begin{array}{cc}{R}^2& 0\\ 0& 1\end{array}\right],B_{\alpha \beta }=\left[\begin{array}{cc}-R& 0\\ 0& 0\end{array}\right],$$

therefore

$$B_{\beta }^{\alpha }=A^{\alpha k}{B}_{k\beta }=\left[\begin{array}{cc}1/R^2& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}-R& 0\\ 0& 0\end{array}\right]=\left[\begin{array}{cc}1/R& 0\\ 0& 0\end{array}\right].$$

Subsequently, it is easy to show that

$$\beta ={\displaystyle \frac{\epsilon E}{1-{\nu }^2}}\left((1-\nu )B_{\beta }^k{B}_k^{\beta }+\nu B_{\beta }^{\beta }{B}_k^k\right)={\displaystyle \frac{\epsilon E}{1-{\nu }^2}}\left(\frac{1-\nu }{R^2}+\frac{\nu }{R^2}\right)={\displaystyle \frac{\epsilon E}{1-{\nu }^2}}{\displaystyle \frac{1}{R^2}}.$$

(12)

Moreover, because the prestress term acts only in the longitudinal dimension, in the cylindrical case we have

$$P=\left[\begin{array}{cc}0& 0\\ 0& P^{zz}\end{array}\right].$$

In this case, the system (11) can be simplified as a one dimensional equation

$$\beta {\eta }_n-{\mu }_s{\displaystyle \frac{{\partial }^2{\eta }_n}{\partial z^2}}=f_s,$$

(13)

where ${\mu }_s=P^{zz}$. See [30,31] for more information on the applications of this model to hemodynamic.

2.4. The Coupled Fluid-Structure Model

The fluid is modeled as Newtonian, homogeneous, and incompressible. The fluid satisfies the following equations in ALE form  [2,32]

$$\begin{array}{cccc}\hfill & {\rho }_f{\displaystyle \frac{\partial \mathit{u}}{\partial t}}{|}_{\mathcal{A}}+{\rho }_f[(\mathit{u}-\mathit{w})·\nabla ]\mathit{u}-\nabla ·{\sigma }^f=\mathit{0}\hfill & \hfill & \mathrm{on}{\Omega }_f,\hfill \end{array}$$

(14)

$$\begin{array}{cccc}\hfill & \nabla ·\mathit{u}=0\hfill & \hfill & \mathrm{on}{\Omega }_f,\hfill \end{array}$$

(15)

where ${\rho }_f$ and $\mathit{u}$ are the density and the velocity vector of the fluid, respectively. The Cauchy stress tensor of the fluid is written as ${\sigma }^f=-p\mathit{I}+\mu (\nabla \mathit{u}+\nabla {\mathit{u}}^T)$, where p and $\mu$ are the pressure and dynamic viscosity of the fluid, respectively. The system of Equation (14) is completed with appropriate boundary conditions. Moreover, ${\Omega }_f$ is the fluid domain and $\mathit{w}$ is the ALE velocity that determines, step-by-step, the position of the nodes of the fluid domain as ${\mathit{x}}_f\left(t\right)={\mathit{x}}_0+{\int }_0^t\mathit{w}d\tau$.

In the following, $(·,·)$ will denote the $L^2(\Omega \left(t\right))$ inner product, ${(·,·)}_{{\Gamma }_s}$ will denote the $L^2\left({\Gamma }_s\right)$ inner product, and the bilinear form $a(·,·)$ is defined as

$$a(\mathit{w},\mathit{v})=\mu (\nabla \mathit{w}+{\left(\nabla \mathit{w}\right)}^T,\nabla \mathit{v}).$$

The optimal control that is introduced in this work is based on the stationary solution of the fluid–structure system. Therefore, in the following, we neglect all of the time derivatives. We can write the weak formulation for the fluid problem as

$$\begin{array}{cc}\hfill & a(\mathit{u},\mathit{\varphi })+{\rho }_f(((\mathit{u}-\mathit{w})·\nabla )\mathit{u},\mathit{\varphi })-(p,\nabla ·\mathit{\varphi })\hfill \\ \hfill & ={\int }_{\Gamma }\left({\sigma }_f{\tilde{\mathit{a}}}_n\right)·{\tilde{\mathit{a}}}_n(\mathit{\varphi }·{\tilde{\mathit{a}}}_n)d\Gamma +{\int }_{{\Gamma }_N}\mathit{h}·\mathit{\varphi }d\Gamma ,\hfill \\ \hfill & (q,\nabla \mathit{u})=0,\hfill \end{array}$$

(16)

for all $\mathit{\varphi }\in {\mathit{H}}^1{(\Omega ):\mathit{\varphi }|}_{{\Gamma }_D}=0$ and $(\mathit{\varphi }·{\tilde{\mathit{a}}}_{\alpha }){|}_{\Gamma }=0$ and for all $q\in L^2(\Omega )$.

With a similar approach, we can write the weak formulation for the prestressed shell Equation (11), as

$${({\beta }^{*}{\eta }_n,\psi )}_{{\Gamma }_s}+{(\mathit{P}\nabla {\eta }_n,\nabla \psi )}_{{\Gamma }_s}={(f_s,\psi )}_{{\Gamma }_s},\forall \psi \in H^1\left({\Gamma }_s\right).$$

(17)

By using this shell model, the structure equations can be reduced to a boundary condition on ${\Gamma }_s$ for the fluid problem. Therefore, the two sub-systems (16) and (17) are coupled by imposing ${\sigma }^f·\mathit{n}-f_s=0$ on ${\Gamma }_s$. Now, we define the functional space $V^0=\{\mathit{\varphi }\in H^1\left({\Omega }_f\right):\mathit{\varphi }{|}_{{\Gamma }_{D,f}}=\mathit{0}\}$, where ${\Gamma }_{D,f}$ are the boundaries of ${\Omega }_f$ where Dirichlet conditions are imposed. In order to satisfy the continuity of the test functions $\mathit{\varphi }·\mathit{n}=\psi$ over the interface surface ${\Gamma }_s$ in the coupled system, a new functional space is introduced as

$$W^0=\{(\mathit{\varphi },\psi )\in V^0\times H^1\left({\Gamma }_s\right):\mathit{\varphi }·\mathit{n}=\psi \mathrm{over}{\Gamma }_s\}.$$

(18)

Now, we can derive the weak form of the coupled final system by simplifying the two terms ${\int }_{\Gamma }\left({\sigma }_f{\tilde{\mathit{a}}}_n\right)·{\tilde{\mathit{a}}}_n(\mathit{\varphi }·{\tilde{\mathit{a}}}_n)d\Gamma$ and ${(f_s,\psi )}_{{\Gamma }_s}$. Note that $\psi =(\mathit{\varphi }·{\tilde{\mathit{a}}}_n)\circ \mathit{\theta }\left(\mathit{x}\right)$. Subsequently, the coupled system reads

$$\begin{array}{cc}\hfill & {\rho }_f\left(\left[\right(\mathit{u}-\mathit{w})·\nabla ]\mathit{u},\mathit{\varphi }\right)+a(\mathit{u},\mathit{\varphi })-(p,\nabla \mathit{\varphi })+{({\beta }^{*}{\eta }_n,\psi )}_{{\Gamma }_s}\hfill \\ \hfill & +{(\mathit{P}\nabla {\eta }_n,\nabla \psi )}_{{\Gamma }_s}={\int }_{{\Gamma }_{N,t}}\mathit{h}·\mathit{\varphi }d\Gamma ,\hfill \\ \hfill & (\nabla ·\mathit{u},q)=0.\hfill \end{array}$$

(19)

for all $(\mathit{\varphi },\psi )\in W^0$, $q\in L^2\left({\Omega }_f\right)$. A finite element technique is used to obtain the discrete weak formulation of (19). With this approach, the structural equation can be incorporated in the fluid equations as Robin boundary conditions.

3. Optimality System

In this work, we are interested in solving a displacement matching optimal control problem. We aim to obtain a desired shell deformation profile by controlling the fluid pressure over a portion of its boundary. Therefore, the objective functional that we intend to minimize reads

$$\mathcal{J}(\mathit{\eta },p)=\frac{1}{2}{\int }_{{\Gamma }_d}{\parallel \mathit{\eta }-{\mathit{\eta }}_d\parallel }^{2}d\Gamma +\frac{\lambda }{2}{\int }_{{\Gamma }_c}{p}^{2}d\Gamma .$$

(20)

The first term expresses the distance in norm between the actual displacement and its desired value over the controlled boundary ${\Gamma }_d$, which can either be the whole moving wall or one sub-portion. The second term has been added to limit the $L^2$-norm of the control, the fluid boundary pressure p. The regularization parameter $\lambda$ weighs the importance of the two terms over the cost functional and the choice of its value can be challenging. As a general rule, too much regularization leads to smoother, but less effective, controls, while a lack of regularization may lead to numerical and convergence issues.

In this work, we use a standard Lagrange multiplier approach to derive the first-order necessary conditions. We first introduce the augmented functional $\mathcal{L}$, which is obtained by adding to the objective functional $\mathcal{J}$ the FSI state Equation (19) multiplied by the appropriate Lagrange multipliers, also known as adjoint variables.

$$\begin{array}{cc}\hfill & \mathcal{L}(\mathit{\eta },\mathit{u},{\mathit{u}}_a,p,p_a,\Gamma )=\mathcal{J}(\mathit{\eta },p)-({\rho }^f(\mathit{u}·\nabla )\mathit{u},{\mathit{u}}_a)+(p,\nabla ·{\mathit{u}}_a)\hfill \\ \hfill & -(p_a,\nabla ·\mathit{u})-\mu (\nabla \mathit{u},\nabla {\mathit{u}}_a)+{\int }_{\Gamma }\mu (\nabla \mathit{u}·\mathit{n})·{\mathit{u}}_{a}d\Gamma \hfill \\ \hfill & -{\int }_{\Gamma }\left(p\mathit{n}\right)·{\mathit{u}}_{a}d\Gamma -{\int }_{{\Gamma }_s}{\mathit{u}}_a·(\beta \mathit{\eta }+\mathit{P}{\nabla }^2\mathit{\eta }-{\mathit{f}}_s-{\tau }_n)d\Gamma ,\hfill \end{array}$$

(21)

Since our control variable is the fluid pressure on the boundary, we integrate by parts the contributions of the fluid stress tensor ${\mathit{\sigma }}^f$. The surface integrals can be rewritten by substituting the definition of ${\mathit{\tau }}_n$. The stationary points of the Lagrangian functional can be found by setting to zero the Fréchet derivatives taken with respect to all the problem variables. When considering the variations of the adjoint variables, we recover the weak form of the state system (19) and its boundary conditions. By taking the derivative in the direction $\delta p$, we get

$$\frac{D\mathcal{L}}{Dp}\delta p=(\nabla ·{\mathit{u}}_a,\delta p)-{\int }_{\Gamma -{\Gamma }_s}({\mathit{u}}_a·\mathit{n})\delta pd\Gamma +{\int }_{{\Gamma }_c}\lambda p\delta pd\Gamma =0\forall \delta p\in L^2(\Omega ).$$

(22)

Because this derivative must vanish for every choice of $\delta p$, then the volume term gives the following continuity equation for the adjoint velocity

$$\nabla ·{\mathit{u}}_a=0\mathrm{on}\Omega ,$$

(23)

while, by setting to zero the surface terms, we recover the following control equation over the controlled boundary ${\Gamma }_c$

$${\int }_{\Gamma -{\Gamma }_s}({\mathit{u}}_a·\mathit{n})\delta pd\Gamma -{\int }_{{\Gamma }_c}\lambda p\delta pd\Gamma =0\forall \delta p\in L^2(\Omega )\Rightarrow p=\frac{{\mathit{u}}_a·\mathit{n}}{\lambda }\mathrm{on}{\Gamma }_c,$$

(24)

and the boundary conditions on ${\Gamma }_D$

$${\int }_{\Gamma -{\Gamma }_s-{\Gamma }_c}({\mathit{u}}_a·\mathit{n})\delta pd\Gamma =0\forall \delta p\in L^2(\Omega )\Rightarrow {\mathit{u}}_a·\mathit{n}=0\mathrm{on}{\Gamma }_D.$$

(25)

We remark that, where we impose Neumann boundary conditions with fixed pressure (i.e., ${\Gamma }_N$). we have $\delta p=0$.

Collecting the variations in the direction $\delta \mathit{\eta }$, we have

$$\frac{D\mathcal{L}}{D\mathit{\eta }}\delta \mathit{\eta }=-{\int }_{{\Gamma }_s}{\mathit{u}}_a·(\beta \delta \mathit{\eta }+\mathit{P}{\nabla }^2\delta \mathit{\eta })d\Gamma +{\int }_{{\Gamma }_d}(\mathit{\eta }-{\mathit{\eta }}_d)·\delta \mathit{\eta }d\Gamma =0\forall \delta \mathit{\eta }\in W^0.$$

(26)

We integrate by parts (26), obtaining the equation for the boundary conditions of the adjoint velocity

$$-{\int }_{{\Gamma }_s}{\mathit{u}}_a·\beta \delta \mathit{\eta }d\Gamma +{\int }_{{\Gamma }_s}\nabla {\mathit{u}}_a·\mathit{P}\nabla \delta \mathit{\eta }d\Gamma +{\int }_{{\Gamma }_d}(\mathit{\eta }-{\mathit{\eta }}_d)·\delta \mathit{\eta }d\Gamma =0\forall \delta \mathit{\eta }\in W^0.$$

(27)

We now collect $\delta \mathit{u}$ terms and integrate by parts, obtaining

$$\begin{array}{cc}\hfill & {\rho }^f\left(\left[(\delta \mathit{u}·\nabla )\mathit{u}+(\mathit{u}·\nabla )\delta \mathit{u}\right],{\mathit{u}}_a\right)-(\nabla p_a,\delta \mathit{u})-\mu ({\nabla }^2{\mathit{u}}_a,\delta \mathit{u})+{\int }_{\Gamma }\left(p_a\mathit{n}\right)·\delta \mathit{u}d\Gamma \hfill \\ \hfill & -{\int }_{\Gamma -{\Gamma }_s}\mu (\nabla \delta \mathit{u}·\mathit{n})·{\mathit{u}}_{a}d\Gamma +{\int }_{\Gamma }\mu (\nabla {\mathit{u}}_a·\mathit{n})·\delta \mathit{u}d\Gamma =0\forall \delta \mathit{u}\in W^0,\hfill \end{array}$$

(28)

which gives the equation for the adjoint velocity

$${\rho }^f{\left(\nabla \mathit{u}\right)}^T·{\mathit{u}}_a+{\rho }^f(\mathit{u}·\nabla ){\mathit{u}}_a-\nabla p_a-\mu {\nabla }^2{\mathit{u}}_a=0\mathrm{on}\Omega ,$$

(29)

with the boundary conditions

$${\mathit{u}}_a=0\mathrm{on}{\Gamma }_D,{\tau }_{na}=0\mathrm{on}{\Gamma }_N\cup {\Gamma }_s.$$

(30)

Moreover, we have to consider the change of $\mathcal{L}$ due to the motion of the boundary ${\Gamma }_s$ along the direction $\delta \mathit{\eta }$

$$\frac{D\mathcal{L}}{D\Gamma }\delta \mathit{\eta }={\int }_{{\Gamma }_s}\beta (\nabla {\mathit{u}}_a·\mathit{n}+\chi {\mathit{u}}_a)·\delta \mathit{\eta }d\Gamma =0\forall \delta \mathit{\eta }\in W^0,$$

(31)

where $\chi$ represents the shell curvature. Under the hypothesis of small deformations, we can safely neglect the terms where $\chi$ appears. The term with ${\mathit{u}}_a$ is defined on the surface ${\Gamma }_s$ and a constant extension of this value towards the normal direction to the surface leads to a null normal gradient, so we have that the moving shape terms $\frac{D\mathcal{L}}{D\Gamma }\delta \eta$ vanish.

To summarize, the optimality system is composed by the state system (19), by the control Equation (24) and by the adjoint system (23)–(28) with boundary conditions (25), (27) and (30). Since the optimality system doubles the number of the state variables, we use a segregated approach for the solution of the state, adjoint, and gradient equations. In this case, we can reuse the same solver for both the solution of the state (19) and adjoint systems (23)–(28) with few modifications. In Algorithm 1, we describe the steepest descent method used.

Algorithm 1 Description of the Steepest Descent algorithm.

1. Set a state $({\mathit{u}}^0,p^0,{\eta }^0)$ satisfying (19) ▹ Setup of the state - Reference case 2. Compute the functional ${\mathcal{J}}^0$ in (20) 3. Set $r^0=1$ for$i=1\to i_{max}$do     4. Solve the system (23)–(28) to obtain the adjoint state $({\mathit{u}}_a^i,p_a^i)$     5. Set the control update $\delta p^i=-(p_c^{i-1}+{\mathit{u}}_a^i·\mathit{n}/\lambda )$     6. Set $r^i=r^0$     while ${\mathcal{J}}^i(p_c^{i-1}+r^i\delta p^i)>{\mathcal{J}}^{i-1}\left(p_c^{i-1}\right)$ do ▹ Line search         7. Set $r^i=\rho r^i$         8. Solve (19) for the state $({\mathit{u}}^i,p^i,{\eta }^i)$ with $p_c^i=p_c^{i-1}+r^i\delta p^i$         if $r^i<toll$ then            Line search not successful ▹ End of the algorithm         end if     end while end for

4. Numerical Results

In this section, we report some of the numerical results that were obtained by using the mathematical model shown in the previous sections. We implemented the presented Koiter FSI model in the multigrid finite element code FEMuS [33], which relies on PETSc libraries for the solution of the multigrid discretized linear solver with MPI parallelization. We compare the results that were obtained with the presented Koiter FSI model with a simplified analytic problem. Subseequently, we further validate the implemented model through the comparison with a validated FSI model.

We also implemented the presented steepest descent algorithm in the code FEMuS, and we present different tests, showing the dependence of the optimal solution on the regularization parameter $\lambda$, the convergence with the grid refinement of the implemented algorithm, and the effectiveness of the used method when non-constant desired target fields are requested.

4.1. Validation of the Koiter FSI Model

Because the presented fluid–structure model based on Koiter shell equations is still not widespread, in the literature there are only a few works on a benchmarking of it. In this section, we will refer to the benchmarks that are present in [34]. In particular, we first consider a simple numerical case, where an analytical solution is available. Another benchmark that is based on a comparison between the Koiter model and a fully validated monolithic fluid–structure model is then presented.

4.1.1. Exact Koiter Solution

We consider a fluid flowing through a cylindrical channel, obtained through the rotation of the rectangle $\Omega =\left\{\right(x,y):x\in [0,0.005],y\in [0,0.06\left]\right\}$ around the z-axis, obtaining a cylinder of radius $R=0.005$ m and height $L=0.06$ m. We consider an inlet at the bottom and an outlet at the top. On the outer wall, the Robin boundary condition for Koiter shell is imposed. We consider a solid shell with density ${\rho }_s=1100\frac{\mathrm{kg}}{{\mathrm{m}}^3}$, a Poisson coefficient ${\nu }_s=0.1$, an elastic modulus $E=124$ kPa, and a wall thickness $\epsilon =2\times {10}^{-4}$ m. Moreover, we consider a fluid with density ${\rho }_f=1000\frac{\mathrm{kg}}{{\mathrm{m}}^3}$ and viscosity ${\nu }_f=0.001\frac{{\mathrm{m}}^2}{\mathrm{s}}$.

Since we are considering a cylindrical geometry, we consider the simplified model (13). By neglecting the prestress term, the problem turns out to be linear, and can be exactly solved under certain conditions. By imposing a constant inlet pressure $p_{in}$, in fact, it is possible to obtain the analytical solution of the pressure and the displacement fields for stationary solutions.

Under the presented parameter setting we obtain $\beta =\frac{\epsilon E}{(1-{\nu }_s)R^2}=1010$ kPa/m. Then, by setting a time step of $0.0005$ s, an inlet pressure of $p_{in}=25$ Pa and an outlet pressure of $p_{out}=0$ Pa, after a time $t=0.25$ s the steady state is reached.

As mentioned in [34], the fluid pressure is linear within the channel, as

$$p_e(r,z)=p_e\left(z\right)=\frac{p_{out}z+p_{in}(L-z)}{L}z\in [0,L],\forall r\in [0,R].$$

(32)

Subsequently, the exact radial displacement of the structure is simply given by

$${\eta }_e\left(z\right)=\frac{p_e\left(z\right)}{\beta }.$$

(33)

The comparison between the displacement field that is simulated with the implemented algorithm and the exact displacement (33) is reported in Figure 2 on the left, showing good agreement between the expected and simulated values of the considered field. However, some discrepancies between the exact and the simulated solutions can be found at the extremes of the cylinders (at the inlet and outlet). This is due to boundary conditions on the displacement field that should be improved in future works. At the same time, in Figure 2 on the right the comparison of the exact (32) and the simulated pressure fields along the cylinder axis are reported. It can be seen that there is total agreement between the simulated pressure and the exact one.

4.1.2. Comparison with a Validated FSI Model

The Koiter FSI model is now compared with the result of a monolithic fluid–structure model that is implemented in the FEMuS code and validated through the Turek FSI benchmarks [1]. The validation and the introduction of the mathematical and numerical model can be found in [35].

Differently from the benchmark that was proposed in [34], we do not consider a multi-layer case, and the comparison is developed on a single solid layer. Because the monolithic model satisfies the Turek benchmark, we use it as a reference for the testing of the new algorithm. The test has been carried out using the same parameters of the previous test, with the only exception of the elastic modulus ($E=15625$ Pa) and the wall thickness ($\epsilon =1.2\times {10}^{-3}$ m). From (12), we obtain $\beta \simeq 757$ kPa/m.

Because the benchmarks available in the literature usually involve transient problems, we consider the time-dependent Koiter model (see [18] for more information) by reformulating (13) as

$${\rho }_s\epsilon {\displaystyle \frac{{\partial }^2{\eta }_n}{\partial t^2}}+\beta {\eta }_n-{\mu }_s{\displaystyle \frac{{\partial }^2{\eta }_n}{\partial z^2}}=f_s.$$

Now, we consider a variable inlet pressure value as

$$p_{in}=\left\{\begin{array}{cc}\frac{p_{max}}{2}\left(1-cos\left(\frac{2\pi t}{t_{max}}\right)\right)\hfill & \mathrm{if}t<t_{max},\hfill \\ 0\hfill & \mathrm{if}t\ge t_{max},\hfill \end{array}\right.$$

where $p_{max}=1.333$ kPa and $t_{max}=0.15$ s. With all of the presented hypotheses, the two different tests have been carried out on the meshes that are presented in Figure 3, with a time step $\Delta t=5\times {10}^{-4}$ s. As boundary conditions, we impose an inlet of the fluid from the bottom and an outlet from the top. On the outer wall (external surface of the cylinders), structure conditions are imposed.

In Figure 4 and Figure 5, the comparison of the displacement field between the Koiter and monolithic model is reported. The displacement of the structure is calculated along the line between the two points $(0.005,0,0)$ and $(0.005,0,0.6)$. It can be noted that the two models have similar behavior in time for all of the reported time steps. In particular, the two models differ from each other for $t=0.02$ s and for $t=0.04$ s. After the initial transition, the two solutions seem to converge to a similar solution, as can be seen for $t=0.06$ s and $t=0.08$ s.

Despite the slight differences between the Koiter and monolithic model noted above, from the presented preliminary results we can assert that the implement Koiter model is consistent with the solution that is obtained with the monolithic model, which has been fully tested in previous works.

Note that the reduction of the computational cost that i sobtained through the Koiter model depends on the solid mesh: the finer is the solid grid, the higher the reduction of the computational cost.

4.2. Regularization Term Tests

We now consider the numerical results of the implemented steepest descent algorithm. We show a simple numerical test to show the dependence of the implemented control algorithm for different values of the regularization term $\lambda$. We consider the simplified case, where the domain ${\Omega }_d$ is reduced to a single point ${\mathit{x}}_d$. We consider a rectangular domain $\Omega =\left\{\right(x,y):x\in [0,0.1],y\in [0,0.3\left]\right\}$, as shown in Figure 6 on the left. The fluid has density ${\rho }^f=1000\mathrm{kg}/{\mathrm{m}}^3$ and dynamic viscosity $\mu =100\mathrm{Pa}·\mathrm{s}$. For the solid, we consider $\beta =60$ kPa/m and thickness $h_s=0.0075$ m. The domain is uniformly divided in a regular rectangular mesh.

We only present one simple case, where the fluid flows vertically from the bottom to the top. The region of the boundary ${\Gamma }_2$ represents a solid wall with no-slip boundary condition ($\mathit{u}=\mathit{0}$) and ${\Gamma }_3$ is the membrane where we impose the generalized Robin boundary condition (6). As an initial condition, we impose a linear pressure field, decreasing from the bottom, where $p=6000$ Pa is imposed, to the top, where we fixed $p=0$ Pa. In Figure 6, we also report the steady solutions of the displacement field (center) and the pressure field (right) when the control algorithm is not considered.

The simulations aim to control the displacement of the point ${\mathit{x}}_d$ of the membrane, in $\widehat{ı}$ direction, optimizing the pressure of the fluid on ${\Gamma }_1$. If we do not apply the control algorithm, the point ${\mathit{x}}_{\mathit{d}}$ shows a displacement $\eta =0.015824$ m. The goal of this test is to reduce it to ${\eta }_d=0.005$ m, thus changing the pressure of the fluid on ${\Gamma }_1$. Thus, the objective functional of the problem reads

$$\mathcal{J}(\eta ,p)=\frac{1}{2}{(\eta {|}_{{\mathit{x}}_d}-{\eta }_d)}^2+\frac{\lambda }{2}{\int }_{{\Gamma }_1}{p}^{2}d\Gamma .$$

(34)

Table 1. Regularization term test. Objective functional and displacement values obtained for different $\lambda$ values. The total number of iteration of the algorithm is also reported.

$\mathit{\lambda }$	$\mathcal{J}$ ($\mathit{\eta },{\mathit{p}}_{\mathit{c}}$)	${\mathit{\eta }}_{\mathit{opt}}$ (m)	Iterations
∞	$5.85839\times {10}^{-05}$	0.015824	-
${10}^{-08}$	$2.19246\times {10}^{-06}$	0.002906	4
${10}^{-09}$	$5.54438\times {10}^{-09}$	0.004895	8
${10}^{-10}$	$2.18941\times {10}^{-10}$	0.004979	10
${10}^{-11}$	$6.10506\times {10}^{-12}$	0.004997	12
${10}^{-12}$	$3.86734\times {10}^{-15}$	0.005000	26

presents the results of multiple simulations, with variations of the regularization parameter $\lambda$, which has been introduced to regularize the control parameter p and force it to the functional space $L^2(\Omega )$. Note that the solution is closer to the objective for small values of $\lambda$. This result was expected, since, with larger $\lambda$, the contribution of the regularization term in the minimization of the functional is more relevant. Thus, with larger $\lambda$, we find smooth pressure profiles, but less precise displacement field $\eta$. In Table, we also report the number of iterations needed to the Algorithm 1 to find the optimal solution.

The controlled pressure $p_c$ is updated via the formula

$$p_c^i=p_c^{i-1}-r^{i,j}\left(p_c^{i-1}-\frac{{\mathit{u}}_a^i·\mathit{n}}{\lambda }\right).$$

(35)

It clearly depends on the regularization parameter. In Figure 7, the controlled pressure field along the boundary ${\Gamma }_1$ is reported for various values of $\lambda$. Note that the choice of the regularization parameter strongly affects the controlled pressure field, as expected. With less regularization, the objective term dominates in the functional and the pressure can have larger values, thus effectively controlling the membrane displacement. The comparison between the imposed reference initial condition and controlled pressure fields in all of the studied cases show that the control algorithm changes the solution on ${\Gamma }_c={\Gamma }_1$ in order to obtain the desired displacement ${\eta }_d$. Moreover, in Figure 7 on the right, the velocity field for $\lambda ={10}^{-15}$ is reported. Note that the pressure field that is imposed by the control algorithm forces an inverse flow at the inlet.

4.3. Grid Convergence Tests

Now, we show some results on the grid convergence of the proposed method. In particular, we require the reduction of the integral of the difference between $\eta$ and ${\eta }_d$ on ${\Omega }_d$ when the grid is refined. Thus, we report different tests, depending on the requested displacement on ${\Omega }_d$.

In all of the considered cases, we consider the same physical values used in the last section. Moreover, we will consider $\lambda ={10}^{-10}$, unless stated otherwise. We consider a $2\times 2$ mesh of the domain that is introduced in Figure 8 on the left, and we refine it with a multi-grid approach. See [36] and reference therein for more information regarding the multi-grid method.

4.3.1. Displacement Reduction Test

We first compute a displacement reduction concerning the reference case, in which we request a desired displacement ${\eta }_d=0.005$ m on ${\Omega }_d$. We consider the geometry in Figure 8 (left) with ${\Omega }_d=\{(x,y):x\in [0.075,0.1],y\in [0.075,0.225]\}$. Moreover, we consider the same boundary conditions that were introduced in the last numerical test. We control the pressure p field over ${\Gamma }_c={\Gamma }_1$. The average displacement on ${\Omega }_d$ in the steady state without control case is ${\overline{\eta }}_d=0.00751$ m. The functional that is associated with the presented control problem in all of the tests reported in the following reads

$$\mathcal{J}(\eta ,p)=\frac{1}{2}{\int }_{{\Gamma }_d}{(\eta -{\eta }_d)}^{2}d\Gamma +\frac{\lambda }{2}{\int }_{{\Gamma }_1}{p}^{2}d\Gamma .$$

(36)

In Figure 8 on the right, the comparison between the controlled pressure field for different grid refinements is reported along the controlled boundary ${\Gamma }_1$. In particular, all of the presented solutions seem to converge with the grid to a certain pressure field. All of the fields differ from the reference pressure field and, in particular, it can be noted that the pressure is reduced by the control algorithm at the inlet. This result is expected, since we are requesting a reduction of the displacement on ${\Omega }_d$.

In

Table 2. Displacement reduction test. Distance from the objective ${\eta }_d$ as a function of the number of refinement levels.

Levels (ndof)	$\mathit{\lambda }$	${\int }_{{\mathbf{\Omega }}_{\mathit{d}}}$ ($\mathit{\eta }-{\mathit{\eta }}_{\mathit{d}}{)}^2\mathit{dx}$	Iterations	R
$2\left(81\right)$	∞	$3.45\times {10}^{-8}$	-	-
$2\left(81\right)$	${10}^{-10}$	$4.01\times {10}^{-9}$	10	$1.16\times {10}^{-1}$
$3\left(289\right)$	${10}^{-10}$	$3.41\times {10}^{-9}$	10	$9.88\times {10}^{-2}$
$4\left(1089\right)$	${10}^{-10}$	$3.26\times {10}^{-9}$	12	$9.45\times {10}^{-2}$
$5\left(4225\right)$	${10}^{-10}$	$3.22\times {10}^{-9}$	10	$9.32\times {10}^{-2}$

, the values of the distance between the desired solution and the solution found with the optimal control algorithm are reported. In particular, we compute the objective distance as ${\int }_{{\Omega }_d}{(\eta -{\eta }_d)}^{2}dx$, where the integral over ${\Omega }_d$ is calculated on the same grid with the same quadrature rule for all of the considered grid refinements.

Note that the distance from the objective decreases when the refinement of the grid is increased. In Table 2 , we also report the number of iterations of the algorithm to find the optimal solution. All of these values are compared with the reference simulation (no control, $\lambda =\infty$). We also report the reduction rate, defined as

$$R={\displaystyle \frac{{\int }_{{\Omega }_d}{(\eta -{\eta }_d)}^{2}dx}{{\int }_{{\Omega }_d}{({\overline{\eta }}_d-{\eta }_d)}^{2}dx}}.$$

(37)

Note that all of the reported values of R are always $R<1$. Therefore, the solution is always improved with respect to the reference one, which indicates that the control algorithm always finds a solution better than the reference solution.

4.3.2. Displacement Increase Test

Now, we consider a displacement increase with respect to the reference configuration. In particular, we request a desired displacement of ${\eta }_d=0.02$ m on ${\Omega }_d=\{(x,y):x\in [0.075,0.1],y\in [0.075,0.15]\}$, which is different when compared to the last considered case. We consider the domain that is reported in Figure 9 on the left. All of the physical quantities are considered equal to the previous case. The average displacement field in ${\Omega }_d$ in the uncontrolled case is ${\overline{\eta }}_d=0.00833$ m.

In Figure 9, on the right, we report the comparison between the control pressure fields over ${\Gamma }_1$ for different grid refinements. In particular, in contrast to the tests that were previously reported, a meaningful increase of the controlled pressure can be noted in this case. This is a direct consequence of the request for an increase of the desired displacement on ${\Omega }_d$. Note that the pressure fields for 4 and 5 levels are coincident, as expected.

In

Table 3. Displacement increase test. Distance from the objective ${\eta }_d$ as a function of the number of refinement levels.

Levels (ndof)	$\mathit{\lambda }$	${\int }_{{\mathbf{\Omega }}_{\mathit{d}}}$ ($\mathit{\eta }-{\mathit{\eta }}_{\mathit{d}}{)}^2\mathit{dx}$	Iterations	R
$2\left(81\right)$	∞	$2.91\times {10}^{-7}$	-	-
$2\left(81\right)$	${10}^{-10}$	$2.39\times {10}^{-8}$	12	$8.21\times {10}^{-2}$
$3\left(289\right)$	${10}^{-10}$	$1.76\times {10}^{-8}$	7	$6.05\times {10}^{-2}$
$4\left(1089\right)$	${10}^{-10}$	$1.58\times {10}^{-8}$	12	$5.43\times {10}^{-2}$
$5\left(4225\right)$	${10}^{-10}$	$1.53\times {10}^{-8}$	11	$5.26\times {10}^{-2}$

, the comparison between all of the tested cases is reported. In particular, all of the cases are compared with the reference one without control (i.e., $\lambda =\infty$). As in the previous case, we report the distance value from the objective calculated through the integral ${\int }_{{\Omega }_d}{(\eta -{\eta }_d)}^{2}dx$. All of the controlled simulations show a minor distance from the objective when compared to the reference case. The reduction rate, as introduced in (37), is $R<1$ for all of the tested grids. This means that the solution always improves when compared to the reference one, and, for higher refinement levels, reduces the distance from the objective.

4.4. Non-Constant Desired Field Tests

We now consider a desired displacement field that depends on the position $\mathit{x}$. This test can be useful for practical applications where a non-constant objective is required. We present two different problems, with a sinusoidal and a step desired field. The geometry of the domain is reported in Figure 8 (left), with ${\Omega }_d=\{(x,y):x\in [0.075,0.1],y\in [0.075,0.225]\}$.

We consider a fluid density ${\rho }^f=1000\mathrm{kg}/{\mathrm{m}}^3$, fluid dynamic viscosity $\mu =1\mathrm{Pa}·\mathrm{s}$, and, for the approximation of the solid to mono-dimensional membrane, we consider $\beta =60$ kPa/m and thickness $h_s=0.0075$ m. Moreover, we will consider $\lambda ={10}^{-6}$, unless stated otherwise. We consider a $2\times 2$ mesh of the domain, refined 3 times with a multigrid technique. We consider an airbag-like case, therefore we impose a no-slip condition on ${\Gamma }_1\cup {\Gamma }_4$, an inlet pressure on ${\Gamma }_2$ and the Koiter boundary condition on ${\Gamma }_3$. We also have ${\Gamma }_c={\Gamma }_2$. We impose an inlet pressure of $p_{in}=600Pa$ and consider $\lambda ={10}^{-6}$.

4.4.1. Sinusoidal Desired Displacement Test

We consider a simple sinusoidal case, where ${\eta }_d$ is defined as

$${\eta }_d=0.01+0.0025sin\left(\frac{2\pi (y-0.075)}{0.3}\right)\left(m\right).$$

The shape of this function has been chosen to be easily reproduced while using the control strategy that is presented in this work.

In

Table 4. Sinusoidal desired displacement test. Distance from the objective and reduction rate R of the objective functional as a function of control algorithm iterations.

It	$\mathit{\lambda }$	${\int }_{{\mathbf{\Omega }}_{\mathit{d}}}$$(\mathit{\eta }-{\mathit{\eta }}_{\mathit{d}}{)}^2\mathit{d}\mathit{x}$	R
-	∞	$4.15\times {10}^{-7}$	-
1	${10}^{-6}$	$1.90\times {10}^{-7}$	$4.59\times {10}^{-1}$
3	${10}^{-6}$	$9.90\times {10}^{-8}$	$2.39\times {10}^{-1}$
5	${10}^{-6}$	$1.37\times {10}^{-8}$	$3.30\times {10}^{-2}$
7	${10}^{-6}$	$9.45\times {10}^{-9}$	$2.28\times {10}^{-2}$
9	${10}^{-6}$	$7.46\times {10}^{-9}$	$1.80\times {10}^{-2}$

, we report the distance from the objective, depending on the considered iteration of the implemented Steepest Descent algorithm. We also report the reduction rate, as defined in (37). For simplicity, we only report the odd iterations. Note that the distance from the objective decreases in all the iteration of the optimization algorithm, and reach a final value of $7.46\times {10}^{-9}$. The final reduction rate of $1.80\times {10}^{-2}$ suggests that the final solution is closer to the objective when compared to the initial solution.

In Figure 10, the displacement field along the line ℓ between the points $(0.1,0.075)$ and $(0.1,0.225)$ (such that $\ell \in {\Omega }_d$) is reported. Note that the optimal and reference solutions have a similar behavior; therefore, we can conclude that the optimization algorithm found a good solution, in agreement with the desired displacement.

4.4.2. Step Function Desired Displacement Test

To test the presented algorithm on sharp desired functions, we now consider a step function ${\eta }_d$, defined as

$${\eta }_d=\left\{\begin{array}{cc}0.0075m\hfill & \mathrm{if}y\in [0.075,0.15],\hfill \\ 0.015m\hfill & \mathrm{if}y\in [0.15,0.225].\hfill \end{array}\right.$$

All of the physical parameters are the same used in the last cases. Again, we consider an airbag case, where ${\Gamma }_c={\Gamma }_2$. The presented case is not a straightforward control problem, since it is not easy to represent a step function displacement by controlling the pressure on ${\Gamma }_2$. We split the distance from the objective into two terms ${\epsilon }_1$ and ${\epsilon }_2$, such that

$${\epsilon }_i={\int }_{{\Omega }_{d,i}}{(\eta -{\eta }_d)}^{2}d\mathit{x},i=1,2,$$

where ${\Omega }_{d,1}=\{(x,y):x\in [0.075,0.1],y\in [0.075,0.15\left[\right\}$ and ${\Omega }_{d,2}=\{(x,y):x\in [0.075,0.15],y\in [0.075,0.225]\}$ with ${\int }_{{\Omega }_d}{(\eta -{\eta }_d)}^{2}d\mathit{x}={\epsilon }_1+{\epsilon }_2$.

In

Table 5. Step function desired displacement test. Distance from the objective for both the domains ${\Omega }_{d,1}$ and ${\Omega }_{d,2}$, and overall distance from the objective $\epsilon$ as a function of the algorithm iteration. The reduction rate R is also reported.

It	$\mathit{\lambda }$	${\mathit{\epsilon }}_1$	${\mathit{\epsilon }}_2$	$\mathit{\epsilon }$	R
-	∞	$3.64\times {10}^{-7}$	$7.63\times {10}^{-8}$	$4.40\times {10}^{-7}$	-
1	${10}^{-6}$	$2.23\times {10}^{-8}$	$2.17\times {10}^{-7}$	$2.39\times {10}^{-7}$	5.43 $\times {10}^{-1}$
3	${10}^{-6}$	$1.05\times {10}^{-8}$	$7.97\times {10}^{-8}$	$9.02\times {10}^{-8}$	$2.05\times {10}^{-1}$
5	${10}^{-6}$	$7.74\times {10}^{-8}$	$5.49\times {10}^{-9}$	$8.28\times {10}^{-8}$	$1.88\times {10}^{-1}$
7	${10}^{-6}$	$2.31\times {10}^{-8}$	$1.33\times {10}^{-8}$	$3.46\times {10}^{-8}$	$7.87\times {10}^{-2}$
9	${10}^{-6}$	$1.28\times {10}^{-8}$	$1.49\times {10}^{-8}$	$2.77\times {10}^{-8}$	$6.29\times {10}^{-2}$

, we report the values of ${\epsilon }_1$ and ${\epsilon }_2$ for various iterations of the algorithm. Note that the values of ${\epsilon }_1$ and ${\epsilon }_2$ are not uniformly decreasing with the algorithm iterations, but their sum does. In fact, the algorithm is designed to reduce the value of $\epsilon$: such value and, consequently, the reduction rate R, decrease to the values $2.77\times {10}^{-8}$ and $6.29\times {10}^{-2}$, respectively, after 9 iterations. The overall reduction rate indicates that the algorithm finds a solution that is consistently better than the initial one.

In Figure 11, the displacement field along ℓ is reported. It is very difficult to match perfectly the desired displacement ${\eta }_d$ due to the nature of the desired function. The optimization algorithm works to minimize the distance between the optimal and desired solutions, and the solution that is found after 9 iterations of the optimization algorithm is reported in Figure. We remark that, even if the proposed optimal solution does not represent a perfect matching with the desired one, it is an upgrade of the initial solution by $1/R\approx 16$ times.

5. Conclusions

In this work, we have studied a fluid–structure interaction optimization problem based on a multi-scale model, where the thickness of the solid wall can be neglected and, therefore, greatly reduces rhe computational cost of the algorithm. We have stated the problem in monolithic form, where the membrane equation has been taken into account as a Robin boundary condition imposed on the moving boundary. We have obtained the optimality system of a pressure boundary control problem, to find a desired displacement of the membrane. The optimal control problem has been solved by using the Lagrange multiplier method and the gradient of functional has been determined by solving the adjoint problem.

The Koiter shell model has been validated through some numerical benchmarks, which involve the comparison with an analytical result on a simple case, and the comparison with a fully-validated monolithic FSI model.

Subsequently, a Steepest Descent algorithm has been introduced to solve the presented control problem, and it has been implemented in a finite element code. Some two-dimensional numerical tests have been introduced to validate the implemented algorithm. In particular, the dependence of the algorithm on the regularization parameter and its convergence with the grid refinement have been reported. The algorithm shows good convergence properties, and the dependence on the regularization term is consistent with the proposed mathematical model. Finally, the algorithm has been tested with non-constant objectives and, in particular, we considered a sinusoidal and a step-function desired displacement field. The algorithm found a sinusoidal solution in agreement with the desired one, and an improved, although not perfect, step-function solution.

The implemented algorithm has been proven to find a solution closer to the desired one in all of the tested cases in comparison to the uncontrolled solution. Moreover, better solutions are found as the regularization term tends to zero and the solutions converge with the grid refinement, as expected. Therefore, the presented algorithm aims to be an efficient tool for optimal boundary control problems that are applied to fluid–structure interaction simulations.