An optimal control method for time-dependent fluid-structure interaction problems

In this article, we derive an adjoint fluid-structure interaction (FSI) system in an arbitrary Lagrangian-Eulerian (ALE) framework, based upon a one-field finite element method. A key feature of this approach is that the interface condition is automatically satisfied and the problem size is reduced since we only solve for one velocity field for both the primary and adjoint system. A velocity (and/or displacement)-matching optimisation problem is considered by controlling a distributed force. The optimisation problem is solved using a gradient descent method, and a stabilised Barzilai-Borwein method is adopted to accelerate the convergence, which does not need additional evaluations of the objective functional. The proposed control method is validated and assessed against a series of static and dynamic benchmark FSI problems, before being applied successfully to solve a highly challenging FSI control problem.


Introduction
Fluid-structure interactions (FSI) are ubiquitous throughout natural and industrial flow systems. In the human body, for example, almost all fluid-conveying vessels are flexible with even small changes in internal fluid pressure inducing strong FSI which can determine a vessel's biological function or dysfunction (Grotberg and Jensen 2004). Important physiological examples include pulse-wave propagation in arteries, wheezing during exhalation and flow-induced deformation and ultimate rupture of arterial and cerebral aneurysms (Heil and Hazel 2011). Industrial examples are equally common, and the motion of fluid through a domain with School of Engineering, University of Edinburgh, Edinburgh, UK elastic boundaries is extremely important in, for example, the adhesion of elastic elements and in hydraulic fracturing for shale-gas production (Box et al. 2018). Aeroelastic phenomena are also extremely important in aviation, and in particular FSI determine whirl flutter instabilities associated with rotorcraft and propeller systems, buffet responses and dynamic responses of aero-structures in gusts (Beran et al. 2017).
Computational FSI has developed rapidly and reached a significant level of maturity, with broad applications including aerodynamics (Bazilevs et al. 2013;Tezduyar and Sathe 2007;Davidson et al. 2012), biomechanics (Moireau et al. 2012;Bazilevs et al. 2010a, b) and ocean mechanics (Bai and Taylor 2009;Finnegan and Goggins 2012;Calderer et al. 2014). A general FSI system has a solid surrounded by a fluid (such as aeroplane) or a solid surrounding a fluid (such as for blood vessels). The interface between the fluid and solid of a FSI system is an unknown of the coupled system, which can be solved by imposing continuity conditions through the interface, such as continuity of velocity and normal stress. Numerical methods for FSI problem have been intensively studied during the past decades. Based upon mesh types, the FSI numerical methods can be broadly categorised into fitted-mesh methods (Heil 2004;Hecht and Pironneau 2017;Tanaka and Kashiyama 2006) and unfittedmesh methods (Peskin 2002;Zhang et al. 2004;Baaijens 2001;Boffi and Gastaldi 2016;Kreissl and Maute 2012); based upon solving strategies, these numerical methods include partitioned/segregated methods (Küttler and Wall 2008;Degroote et al. 2009Degroote et al. , 2013Bazilevs et al. 2013) and monolithic/fully coupled methods (Heil 2004(Heil , 2008Muddle et al. 2012;Wang et al. 2017Wang et al. , 2020; and based upon solving variables, there are one-field (one velocity) methods (Wang et al. 2017(Wang et al. , 2019aHecht and Pironneau 2017) and multi-field (velocity, displacement and Lagrange multiplier) methods (Muddle et al. 2012;Boffi et al. 2015Boffi et al. , 2016. In a recent study , we proposed an energystable one-field monolithic method based on an ALE fitted mesh, which is adopted to derive the optimality condition for the FSI control formulation in this article. Optimal control is a classical theory which seeks solutions for a stationary or dynamical system such that an objective function can be optimised (Tröltzsch 2010). Adjoint-based methods have solid mathematical foundations which can efficiently compute the gradient of a function, and have been widely adopted to solve optimal control problems. Optimal flow control has been fully developed in the past decades, including distributed control (Gunzburger 2003;Abergel and Temam 1990;Hinze et al. 2012): solution existence of a static flow control by a distributed force is analysed in Abergel and Temam (1990) and a piecewise optimal control method for dynamical case has been presented at the same time; a distributed flow control method is implemented by the multigrid method in Hinze et al. (2012). Optimal flow control also includes shape optimisations from pioneering work (Pironneau 1973;1974;Glowinski and Pironneau 1975;Mohammadi and Pironneau 2010) to recent developments in Montenegro-Johnson and Lauga (2015), Henrot and Privat (2010), Dapogny et al. (2018), and Jenkins and Maute (2016): the first-order optimality condition is derived for the Stokes flow (Pironneau 1973) and Navier-Stokes flow (Pironneau 1974) by a variational approach, and a numerical algorithm is presented in Glowinski and Pironneau (1975); a different minimum-drag profile of Stokes flow has been derived in Montenegro-Johnson and Lauga (2015) given a fixed surface, which has a slender shape compared with the classical results in Pironneau (1973Pironneau ( , 1974 given a fixed volume; the classical optimisation of a solid surrounded by a fluid has been extended to the case of a solid surrounding a fluid in Henrot and Privat (2010), in which the authors prove that a cylindrical pipe is not the optimal shape for minimising the viscous dissipation of its inside flow, and a useful implementation of this method using FreeFEM++ is presented in Dapogny et al. (2018). However, few studies of optimal FSI control, in particular time-dependent FSI problems, have appeared in the literature to date and there is no recognised benchmark solution for comparison. Although there has been some progress reported in Heners et al. (2018), Chirco and Manservisi (2020), Failer and Richter (2020), and Wick and Wollner (2020), these are still very challenging both analytically and numerically. The problem falls into the category of inverse FSI problems of moving shape control (Moubachir and Zolesio 2006). The main challenges are summarised below.
A particular challenge in achieving optimal FSI control is to formulate the coupling conditions on the fluid-structure interface into the optimality system (Chirco et al. 2017;Failer et al. 2016). For example, Failer et al. (2016) enforce the coupling condition weakly to analyse optimal control for a linear FSI problem, whereas Chirco and Manservisi (2020) and Chierici et al. (2019) introduce an auxiliary mesh displacement in the solid domain to enforce the coupling condition. In Failer and Richter (2020) and Wick and Wollner (2020), the authors solve both solid velocity and displacement, together with fluid velocity and pressure using a monolithic Newton solver. In our previous work, we formulated the FSI problem using a one-field FEM scheme (Wang et al. 2017(Wang et al. , 2019a and solved it in a fully coupled system, so that the interface conditions are satisfied automatically in the primal FSI equations and are therefore not present in the adjoint FSI equations. As a result, both the primal and adjoint FSI system may be solved in the same manner as a fluid-only problem.  Turek and Hron (2006) Another difficulty in the optimal dynamic control, of either fluid or FSI problems, is that the adjoint equations are coupled with the primal equations in both time and space (Hinze et al. 2012;Abergel and Temam 1990). A major challenge in the time coupling arises because the adjoint problem is propagated backwards in time, which uses the solution of the primal problem as initial conditions. The space coupling challenge is due to the non-linear terms in the primal equations, which introduce the solution of the primal equations to the adjoint equations. In Degroote et al. (2013), a multigrid is used to solve a single spacetime discretisation of the primal problem. Alternatively, one may apply the piecewise-in-time optimal control approach, i.e. approximating the dynamic control by a sequence of local steady-state problems (Abergel and Temam 1990). The time-coupled primal and adjoint problems can also be decoupled and solved in an iterative manner as done in Heners et al. (2018), which also solves the FSI subproblems using a partitioned/decoupled method. In this article, we decouple the primal and adjoint equations in time; however, we solve both the FSI system and the adjoint FSI system using a one-field monolithic approach .
Finally, efficiently selecting a step size for a minimisation algorithm, such as the gradient descent (GD) method, is usually costly. The Armijo rule is widely used in the context of the GD algorithm (Mohammadi and Pironneau 2010;Gerdes et al. 2014); however, this is computationally expensive for dynamic control problems because one needs several evaluations of the objective function, consequently several computations of the state equations. A constant step size (Heners et al. 2018) does not need evaluation of the objective functional; however, this is typically quite inefficient since it may need many iterations in order to reduce the objective sufficiently. In this article, we adopt a stabilised Barzilai-Borwein (BB) method to accelerate the iteration (Burdakov et al. 2019), which saves significant numerical effort and guarantees the quality of step size at the same time.
The original contributions of the paper are summarised as follows: (1) we consider general boundary and initial conditions, consider a complete coupling between the adjoint FSI velocity and the velocity of configuration, and derive the optimality conditions for optimal control of dynamic FSI problems in an arbitrary Lagrangian-Eulerian (ALE) formulation; (2) we formulate the whole control system using one velocity field, and consequently reduce the number of degrees of freedom of the FSI problem. This one-field formulation is very similar to a pure fluid-control formulation, which allows us to do similar analysis as a fluid control method (Abergel and Temam 1990); (3) we adopt a stabilised BB method to select the step size, which does not need additional evaluations of the objective, and has the same cost as using a constant step size but converges faster. Because of the above features, we demonstrate that it is computationally feasible to solve extremely challenging optimal FSI problems, by solving tens (at most hundreds) of CFD simulations.
The paper is organised as follows. In Section 2, the control equations for the FSI problem are introduced in an ALE framework. In Section 3, the optimisation problem is introduced, followed by gradient descent method in Section 4. The main derivation of the optimality system using the Lagrange multiplier method is presented in Section 5. Discretisation and implementation are considered in Section 6. Numerical results are reported in Section 7, with conclusions drawn in Section 8.

PDE for the FSI system in an ALE formulation
We introduce the FSI system using a benchmark test, sketched in Fig. 1, which was first proposed in Turek and Hron (2006)  is a fixed domain with an outer boundary Γ in + Γ w 1 + Γ N and inner boundary Γ w 0 (the circle, as shown in Fig. 1). For notational convenience, denote by Γ w = Γ w 0 + Γ w 1 the wall boundaries, Γ D = Γ in + Γ w on which the Dirichlet boundary condition (BC) is imposed, and Γ = Γ D + Γ N all the boundaries with Γ N being the Neumann boundary on which the zero-normal stress is enforced. In this case, an ALE frame of reference is convenient to describe the FSI system, because it can track the fluid-solid interface Γ t and move arbitrarily elsewhere.
In this article, we consider both an incompressible fluid and an incompressible hyperelastic isotropic solid. We shall only solve for one velocity field in the whole domain, and the conservation of momentum and conservation of mass take the same form in the fluid and solid (just differing in the specific expressions of the stress tensor). Therefore, it is convenient to introduce an indicator function 1 ω (x) = 1 if density, velocity vector and stress tensor respectively. The control partial differential equations, with initial and boundary conditions, for the FSI problem can then be expressed as follows.
Initial condition: u| t=0 = u 0 , ALE initial condition: w| t=0 = w 0 , Continuity of velocity: Continuity of ALE velocity: Continuity of normal stress: In the above, w is the velocity vector of the ALE frame, f is a distributed control variable, μ and λ are artificial Lamé constants for the ALE equation (Richter and Wick 2010), and n and τ are the normal and tangential direction of the outer boundary, respectively, as shown in Fig. 1. The stress tensor of an incompressible Newtonian flow is expressed as: with D(·) = ∇(·) + ∇ T (·), and μ f being the viscosity parameter (which is unrelated to the artificial Lamé constant μ used in the ALE (3)). The stress tensor of an incompressible neo-Hookean solid is expressed as Hecht and Pironneau (2017): with d being the solid displacement, and c 1 being the elasticity parameter. Notice that although the solid stress tensor is expressed as a function of displacement d, we shall not solve for d as an independent variable. Instead we view it as a function of velocity, and solve the whole FSI problem based upon a one-field-velocity method, which uses an ALE description for both the fluid and solid equations. When solving the mesh equation, the mesh velocity follows the fluid velocity only at the interface. Afterwards, both the fluid and solid meshes are updated based on the mesh velocity, and this can improve the mesh quality for the fluid as well as the solid ).

The optimisation problem
Let L 2 (ω) be the square integrable functions in domain ω with inner product (u, v) , which has the boundary data u on Γ . We also denote by L 2 0 (Ω) the subspace of L 2 (Ω) whose functions have zero mean values.
We consider the following optimisation problem: reducing the discrepancy between the state velocity u (and/or displacement d) and the goal velocity u g (and/or displacement d g ) profile, in a control region Ω c ⊆ Ω (and/or Ω s t ), by controlling a distributed body force f. Problem 1 Given an objective velocity profile u g and a time interval subject to (1)- (12), where f are the control variables, β 1 and β 2 are the weights of controlling velocity and displacement, respectively, and γ 1 and γ 2 are the weights in order to control the final velocity profile in Ω c and Ω s t respectively. The first term J (·) in (15) is the real objective to be minimised, and the second term is a regularisation term with a regularisation parameter θ. If θ is too large then the real objective is not achieved accurately, whereas if θ is too small, this may cause convergence issues for the numerical scheme.

The gradient descent algorithm
A general method of iteratively solving the optimisation Problem 1 is to use, from an initial point f 0 , the Taylor expansion to expand J (u (f)) around f k : with δĴ u f k [δf] being the Gâteaux variation (see Rall (2014) or A) with respect to f (NOT u) at point f k along the direction δf: The gradient descent with a line search algorithm seeks a direction δf k = 1 such that: 3. And it is the steepest: Then, f is updated as: The exact line search for the step size α k is costly. In this paper, we use a stabilised Barzilai-Borwein (BB) method proposed in Burdakov et al. (2019) to optimise analytical functions (k > 0): where and In the above, (21) is the original BB formula (Fletcher 2005;Dai et al. 2006). Formula (22) is an improvement for nonconvex functions if α k 1 < 0 (Dai et al. 2015), and formulae (20) and (23)  Notice that formula (20) does not specify how to compute α 0 which is used to start the iteration. In this paper, we manually choose a relatively small α 0 so thatĴ u f 0 + α 0 δf 0 <Ĵ u f 0 to start off -the magnitude of α 0 can be determined by f 0 + α 0 δf 0 ∼ ρ s ∂u ∂t . In order to compute δĴ (·)[δf], we introduce the Lagrange multiplier method in the following section.

The Lagrange multiplier method
The constraints for the optimisation Problem 1 can be eliminated by introducing the Lagrange multipliersû,p,ŵ as follows: Notice that the condition of velocity continuity (10) is not included in the above Lagrangian functional, because it is automatically satisfied by treating the FSI system as a onefield velocity problem (Wang et al. 2019a). The following Karush-Kuhn-Tucker (KKT) conditions are the first-order necessary conditions to minimise (24): These equations will be solved in order to further compute δL(·) [δf] ≡ δĴ (·) [δf] due to the arbitrariness ofû,p,ŵ in the Lagrangian functional (24). Notice that the fluid domain Ω f t , solid domain Ω s t and the fluid-solid interface Γ t are all functions of the ALE velocity w, so these shape variations should be considered when taking variation with respect to w. We will compute the shape variation in the Hadamard form (Schmidt and Schulz 2010;Mohammadi and Pironneau 2010) where ∇ n s is the normal gradient operator, and ∇ Γ t is the tangential gradient operator on Γ t , with the tangential divergence being defined as: for a given function v. Note that replacing n s by n f = −n s in (28) does not change the sign of δΓ t (·).
Notice that if we apply the domain variation (27) to the momentum (1), which is a term in (24), its shape variation is zero.

State equation
Taking Gâteaux variation of (24) with respect to the Lagrange multipliersû,p andŵ gives the state equations in a weak formulation (see Appendix 2 for more details): Problem 2 Given u 0 in (4) and w 0 in (7), for t ∈ (0, T ], where (8) and (5) being satisfied respectively.

Adjoint equation
Taking Gâteaux variation of (24) with respect to the variables u, p and w, we have In the above, we adopt the integration (in time) by parts to obtain and similarly for Θ×Ω ∂ ∂t δw ·ŵ. We also integrate the stress-tensor term by parts (in space) and eliminate the normal-stress term on Γ t using the interface condition (12). In addition, notice that the shape variation Θ×δΓ t (w − u) · w = 0 due to the interface condition (10) using (28).
If we choose the initial condition forû at time t = T (notice that the adjoint equation is solved backwards in time) aŝ (35) and the boundary condition rearrange and integrate the convection terms by parts: and if we also choose the initial condition forŵ at time t = T aŝ and boundary condition then the adjoint-ALE-FSI equation for (û,p,ŵ) is given by: and Remark 1 We give the Dirichlet boundary condition in (36) for Problem 3 in the weak form, where the Neumann boundary condition has been included. If the stress term is integrated by parts in (40), we would have the Neumann boundary condition σ û,p n + (n · u)û = 0 on Γ N for the corresponding problem in the PDE (partial differential equation) form.

Gradient descent direction
Taking Gâteaux variation of (24) with respect to f, we have From which it can be seen the gradient descent direction is δf1 Ω s t =û − θ f.

Discretisation for the FSI control system
We discretise the time interval [0, T ] as t 0 = 0, t 1 , t 2 , . . . , t M = T with t n+1 − t n = t (n = 0, 1, . . . , M), and then discretise the displacement using the backward Euler scheme: where A t n ,t n+1 is the ALE mapping from state Ω t n to Ω t n+1 : x n → x n + tu n+1 . At the same time, we have δd = tδu. The state equations (30) and (31) can then be discretised forward in time, and the adjoint equations (40) and (41) can be discretised backward in time below in Problems 4 and 5 respectively.
Remark 2 Coupling between u n+1 and w n+1 : the contribution of u n+1 to w n+1 is through Γ t n+1 in the boxed term in (45). In practice we may replace u n+1 by u n in this boxed term, decouple (44) and (45), solve them one by one.
and ρ Ωŵ Remark 3 Coupling betweenû n andŵ n : the only contribution ofŵ n toû n is through the boxed term in (46). In practice we may replaceŵ n byŵ n+1 in this boxed term, decouple equations (46) and (47) and solve them one by one. We have numerically investigated the boxed term in (46), and plotted its values for tests 7.3, 7.4 and 7.6. We find that this term is very small and can be neglected, which means it is not necessary to solve (47).
The optimality condition (42) gives the update of the control variable: with α k being computed using the modified Barzilai-Borwein method (20). We use a standard Taylor-Hood finite element (Q 2 /Q 1 ) for the velocity-pressure pair to discretise (44) and (46), and Q 2 element for the mesh velocity to discretise (45) and (47). Finally the whole algorithm is summarised in Algorithm 1.

Remark 4
For time-dependent control problems, the classic piecewise static control method cannot guarantee a global optimisation solution (Abergel and Temam 1990). In the proposed algorithm, if we just run the algorithm in one time step, the algorithm would be one step of the piecewise static control.

Numerical experiments
Considering the complicated features of the control problems, we discuss in this article, namely time-dependent, fluid-structure interaction problems with finite solid deformations, there are limited cases considered in literature. Therefore, in Section 7.1, we first validate the proposed optimal control method against an existing fluid-control result, obtained using a different approach to ours. We then test the proposed method using three dynamic FSI problems in the rest of this section: test 7.3 is a modification of the fluid problem 7.1, which shows that the control variable, the body force f, is feasible; test 7.4 is a static FSI control problem from the literature, for which we use our dynamic control method and show that the FSI control is also tractable. Finally, test 7.6 is a dynamic FSI problem involving complicated solid oscillations, which does not have a steady-state solution, and we show that it is controllable using the proposed methodology.

Lid-driven cavity flow -forcing a Navier-Stokes flow to a Stokes flow
In this example, we reproduce published results for control of a dynamic cavity pure fluid flow, which has been studied in Failer and Richter (2020) using a space-time multigrid method. We show that the proposed scheme, requiring only a relatively small number of CFD simulations, can achieve equivalent results to those reported in Hinze et al. (2012). To achieve this in our algorithm, we turn off the solid part of the proposed scheme, remove the integration in the solid domain, set the ALE velocity w = 0 and apply the control force f in the whole domain (fluid only in this case). Specifically speaking, we only solve (44) and (46) without terms of integration on Ω s t n , Ω s t n+1 and Γ t n . The cavity flow is defined in domain Ω = [0, 1] × [0, 1], with velocity being prescribed as u = (1, 0) at the top of the boundary and u = (0, 0) at all the other boundaries. The fluid's density and viscosity are ρ f = 1.0 and μ f = 0.01 respectively. We use 30 × 30 Q 2 /Q 1 elements, and a fully developed Navier-Stokes (NS) and Stokes flow are shown in Fig. 2a and b respectively. In this flow control test, we choose the fully developed NS flow as the initial flow and determine the body force distribution needed to create the fully developed Stokes flow as the target flow (we run the simulation up to t = 10 using Δt = 0.01), and solve Problem 1 with Ω c = Ω, T = 1, β 1 = 1, β 2 = 0, γ 1 = 1 and γ 2 = 0.
We first test the influence of the regularisation parameter θ as shown in Fig. 3a, from which it can be seen that larger θ stops the objective decreasing at a earlier stage. We then compare the BB method against gradient descent using a constant step size, and the convergence of the objective functional J (·) is shown in Fig. 3b, from which it can be seen that the BB method is not sensitive to the initial step size and converges slightly faster. It may also be seen that a constant step size α = 4 does not converge, while starting from a different initial step size (α 0 =1, 2 or 4) presents similar convergence for the BB method. Overall, the BB method reduces the objective by over 80% after 10 iterations. We also report the L 2 norm of the control force and the gradient descent direction in Fig. 4. Notice that this problem has a steady-state solution and we have a control for the final velocity profile (γ 1 = 1); this is why the control force tends to a δ−function at t = 0 and t = T (see Remark 5). Figure 5 depicts the controlled flow and the control force at different times, from which we can see the two main vortices quickly "push" the NS flow to the Stokes flow. These results are very similar to Figure 2 in Hinze et al. (2012).
Remark 5 If the control problem has a steady-state solution, then the inertia tends to 0 as time involves. In this case, we know that (1) becomes a stationary PDE equation, which naturally holds without any initial conditions or control force. Therefore, if we want u| t=0 = u 0 (initial condition (4)), we expect a control force f that tends to a δ−function at t = 0 in order to balance this initial condition. If we also  (16)), we also expect the control force f tends to a δ−function at t = T .

Cavity flow with an initial pulse -forcing a Navier-Stokes flow to a predefined time-dependent fluid field
In this example, we consider again a fluid-control problem but with a different objective: steering the velocity to be a complicated predefined velocity profile with vortices, which is taken from Hou and Yan (1997). The computational domain and mesh are exactly the same as the test in Section 7.1. We use a wall boundary condition for all four sides of the cavity, and the fluid with ρ f = 1 and μ f = 0.1 is initially driven by: u 0 (x, y) = ((cos(2πx) − 1) sin(2πy), sin(2πx) (1 − cos(2πy))) .
The goal velocity is derived from the following stream function: We use a time step of Δt = 0.01, run the simulation from t = 0 to 1, and compute the control force f(x, y, t) in the whole fluid domain. We compare the goal velocity and the controlled velocity in Fig. 6, from which it can be seen the four vortices have been clearly captured, and the L 2 error has the same magnitude reported in Figure 10 of Hou and Yan (1997). We test three different regularisation parameter θ , and its influence on the objective reduction may be observed from Fig. 7, from which we see again that the BB method performs better than a constant step method (Fig. 8).

Lid-driven cavity flow with an elastic solid wall -reducing the solid deformation
The next problem increases in complexity to consider the case of lid-driven cavity flow with a deformable solid as considered by Zhang et al. (2012), which has the same geometry as the above test in Section 7.1. However, there is a rectangular solid at the bottom of this square as shown in Fig. 9 (l = 1 and h = 0.25). In this case, there is no interior boundary Γ w 0 as shown in Fig. 1. However, this would not change (4) to (47) that we solve. The fluid and solid properties are ρ f = ρ s = 1, μ f = 0.01 and c 1 = 0.2. The purpose of this test is to demonstrate the feasibility of using the proposed algorithm to reduce the solid deformation/displacement. Using the one-field Fig. 11 Convergence of the objective Fig. 12 Steady-state displacement at the FSI interface monolithic algorithm, there is little difference (compared with the NS flow in Section 7.1) from the computational point of view. We can use the same boundary conditions, the same mesh and time step Δt = 0.01 as used for the NS flow in Section 7.1. We first run the forward simulation up to t = 10 in order to get a steady-state solution: the vertical velocity and the deformed solid mesh are shown in Fig. 10. We then compute a distributed force, and enforce it on the solid to reduce the solid displacement. This is to say we solve Problem 1 by setting d g = 0, β 1 = 0, β 2 = 1, and γ 1 = γ 2 = 0.
We compare again gradient descent using a constant step size and using the BB method to compute a step size. The BB method still performs better than a constant step size as shown in Fig. 11, while it is as cheap as the constant step  (46) size: using formula (20) to compute the step size. It can be seen from Fig. 11 that the objective is reduced by more than 70% at 30 iterations, although the convergence becomes slow afterwards. Figure 11 also shows the influence of regularisation parameter θ, and we can see that using a smaller θ allows a greater reduction of the objective to be achieved. However, we have observed numerical instability using θ = 10 −5 for this test. The fluid-solid interface is plotted in Fig. 12, from which it can be seen that the additional iterations from 31 to 120 mainly contribute to changing to another mode, with only slight reductions in the magnitude of the displacement. Notice that the objective (15) has an integration through the time and space domain, so it is not a point-wise reduction of the displacement for this example. We also test this example using a finer mesh (40 × 40), and it can be seen from Figs. 11 and 12 that both the objective the FSI interface show very similar performance when using these two different meshes. For this example, the problem has a steady-state solution and we have no control for the final velocity profile (γ 1 = γ 2 = 0), so we observe that the control force tends to a δ−function at t = 0 as shown in Fig. 13. This is consistent with Remark 5. Finally, we investigate the boxed term in (46) in Fig. 14, from which it can be seen that this term is very small. We have tested the case of neglecting this term and found that all the results were identical.

Channel flow interacting with two flexible beams -minimising the velocity discrepancy in a specified region
This example is taken from Cerroni et al. (2016) which solves a quasi-static problem. Here we solve the full Fig. 15 Computational domain for the oscillating beams dynamic FSI problem and show that it also converges to the same stationary solution as reported in Cerroni et al. (2016). The computational geometry is shown in Fig. 15 which is symmetric about the x-axis. For this example, there is no interior boundary Γ w 0 . However, this would not change (4) to (47) that we solve. A parabolic velocity profile is prescribed at the inlet Γ in : The fluid and solid parameters are ρ f = ρ s = 1, μ f = 0.01, and c 1 = 66.67. We use 1210 Q 2 /Q 1 elements with 5019 nodes as shown in Fig. 16, and a converged time step size of Δt = 0.1. We run the forward FSI simulation up to t = 50 until a steady-state flow is obtained, with the flow field and the solid displacement being displayed in Fig. 16, which is consistent with Fig. 3 in Cerroni et al. (2016). We then use this solution for the control problem: increasing the velocity by 20% for all t = 0 to T in a control region Ω c (see Fig. 15) by controlling a distributed force on the solid. This is to say we solve Problem 1 with β 1 = 1, β 2 = 0 = γ 1 = γ 2 = 0 and u g = 1.2u x (t), u y (t) in Ω c , where (u x (t), u y (t)) is the solution (from t = 0 to t = 50) of this forward FSI problem without control (notice that the solution con-verges to a steady state at t = 50). It can be seen from Fig. 17 that the BB method still performs better than the constant step method, and reducing the regularisation parameter θ from 0.01 to 0.001 allows us to reduce the objective more (Fig. 18). It can also be seen from Fig. 19 that the control force converged to a complicated time-dependent distributed force for this FSI control problem (Fig. 20). The final converged velocity is reported in Fig. 21, and the final controlled velocity matches the target velocity well as plotted in Fig. 18. Finally, we plot the boxed term in (46) in Fig. 20, from which it can be seen that this term is negligible.

Oscillating leaflet oriented across the flow direction -forcing the solid to match a time-dependent displacement
We consider an oscillating leaflet in a fluid channel, which has been widely studied in the FSI literature (Baaijens 2001;Heil 2004;Wang et al. 2017;Chierici et al. 2019). The computational domain is a L × H channel with a h × w leaflet located across it as sketched in Fig. 22, and L = 4.0m, H = 1.0m, w = 0.1m and h = 0.8m in this example. A periodic flow condition is prescribed on the inlet and outlet boundaries, given byū x = 15y (2 − y) sin (2πt).  Simulation of the forward FSI problem can be achieved by solving the state equations (Problem 4) without control. We use a 80 × 20 uniform mesh and time step Δt = 0.01 and solve the forward problem from t = 0 to t = 1 β 1 = 0 and β 2 = 10 5 . The whole control is tractable, which can be seen from the convergence of the objective function as depicted in Fig. 25 for three different regularisation parameters θ. It can be seen from Fig. 23b that the leaflet's deformation has be dramatically increased, and the deflection of the leaflet is very close to the target as shown in Fig. 24. We also investigate the control force, and it can be seen from Fig. 26 that the magnitude of the control force is a periodic function which converges to a fixed function of time. The norm of the gradient descent direction δf and the increment of the control force (αδf) are plotted in Fig. 27, from which it can be seen that the BB method cam successfully compute a varying iteration step α so that the increment of the control force converges, and consequently the control force converges as well.

Oscillating flag oriented along the flow direction -reducing the solid deflection
In this section, we consider the FSI control problem of an oscillating flag attached to a cylinder, where the goal is to minimise the solid deflection through the controlled application of a force. The computational domain is a rectangle (L × H ) with a cut hole of radius r and centre (c, c) as shown in Fig. 28. A leaflet of size l × h is attached to the boundary of the hole (the mesh of the leaflet is fitted to the boundary of the hole, see the solid mesh in Fig. 29). The geometry parameters are L = 2.5, H = 0.41, l = 0.35, h = 0.02, c = 0.2 and r = 0.05. The fluid and solid parameters are as follows: ρ f = ρ s = 10 3 , μ f = 1 and c 1 = 2.0 × 10 6 . The inlet flow is prescribed as: We use a mesh of 10,054 nodes and 2448 biquadratic elements as shown in Fig. 29, and a converged time step of Δt = 10 −3 . The simulation results are first validated against the data provided in Turek and Hron (2006) through the oscillation period and amplitude at the tip of the flag as shown in Fig. 30, with the period and amplitude being around 0.526 and 0.035 respectively. These figures have a good agreement with the reference values given in Turek and Hron (2006) with a period and amplitude being 0.530 and 0.034 respectively. We note here that although our neo-Hookean solid model is different from the Saint Venant-Kirchhoff model in Turek and Hron (2006), however, these two models are equivalent when applying to solving this FSI benchmark problem (Kadapa et al. 2018;Hecht and Pironneau 2017;Wang et al. 2020), in the sense that they  present the same numerical results as first presented in Turek and Hron (2006).
We now extend the analysis to consider the control of this FSI system, using the results of this fully converged FSI system as initial conditions. We use the results of t = 6 as initial values and run the simulation from t = 0 to 0.05. The aim of the control problem is to reduce the oscillation amplitude by solving Problem 1 with d g = 0, β 1 = 0, β 2 = 10 6 , and γ 1 = γ 2 = 0. The BB method converges rapidly and reduces the objective by 60% as shown in Fig. 31, which is also faster than using a constant step size. It can also be seen from Fig. 31 that θ = 10 −8 and 10 −9 presents similar convergence for the objective functional J , but θ = 10 −7 stops the reduction of J at an earlier stage, with the magnitude of the control force being plotted  Fig. 32. For a long-term control, such as from t = 0 to t = 0.40 (about one oscillation period), we insert 8 control points at t = 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35 and 0.40, and continuously solve Problem 1: the following control results use the results of the previous time as initial values for the next time (Fig. 33). We then patch together all these piecewise-in-time results as shown in Fig. 34 from which it can be seen that the oscillation magnitude has been dramatically reduced. Notice that our objective is to reduce the solid deflection by a time-dependent body force, and the solution (body force) of this problem is generally not unique or periodic. Therefore, after applying the control force, we would not expect (as observed in Fig. 34) an oscillation of fixed period and amplitude of the flag. For this problem, we also test a case of applying the control force only at the tip (x > 0.57) of the flag, and we find that the whole control     is still tractable as shown in Fig. 34. We patch together the control force in Fig. 35 after gradually solving eight control problems by using the previous results as initial conditions for the following one. Although the continuity of the control force cannot be guaranteed, the displacement is still continuous in the whole time domain. The control force gradually decreases to zero in one phase of control, and then flips its direction (see Fig. 36) in the following phase and gradually decreases to zero again. It can be seen from Fig. 35 that the magnitude of the overall control force decreases as the solid deflection is reduced. Notice again that the control force is not unique, different piecewise control results would lead to different distributions of the control force, although they all may be able to successfully reduce the solid deformation. Finally, we investigate the boxed term in (46) in Fig. 33, and we find that it is negligible again.

Conclusion and future work
Time-dependent FSI control problems with large solid deformations are very challenging to solve and very few examples have appeared in the literature. This paper has made a number of original contributions to this area, including (1) derivation of the optimality condition for optimal control of time-dependent FSI problems in an ALE formulation; (2) formulation of the whole control system using one velocity field, thereby reducing the size of the FSI problem; and (3) adopting a stabilised Barzilai-Borwein (BB) method to select the gradient descent step size, which does not need additional evaluations of the objective, and has the same cost as using a constant step size but converges faster.
Gradient descent methods are widely used in the context of the adjoint-optimal control, with either the Armijo rule (Mohammadi and Pironneau 2010; Gerdes et al.  (Heners et al. 2018) step size typically being adopted; however, the former is costly and the latter is inefficient when applied to time-dependent control problems. In this paper, we use a stabilised BB method to accelerate the iteration which has the same cost as a constant step size and does not need additional evaluations of the objective function.
The proposed optimal control method is validated and assessed by six numerical tests: two pure fluid-control problems, two FSI problems which have steady-state solutions, and two dynamic FSI problems with complicated solid oscillation. It is shown that the complex FSI control systems can be solved using the proposed numerical method. More generally, the above features mean that it is now computationally feasible to solve extremely challenging optimal FSI problems, by solving as few as tens (at most hundreds) of CFD simulations.
In the future work, we shall consider minimising other objectives, such as the drag force, and apply this FSI control method to dynamic shape optimisation problems, such as morphing structure in aerospace engineering.
Noticing the finite element spaces we use in Problem 2, in which both the trial and test functions satisfy boundary condition (8), we have (31), based upon the optimality condition (25) and the initial condition w(0) = w 0 stated in Problem 2 as well.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.

Replication of results
A piece of Fortran code and script are shared at https://github.com/yongxingwang in order to reproduce results presented in the paper, including simulation for the FSI problems and implementation of the gradient descent method.