PARTITIONED SOLUTION OF THE UNSTEADY ADJOINT EQUATIONS FOR THE ONE-DIMENSIONAL FLOW IN A FLEXIBLE TUBE

For gradient-based optimization, the gradient of the objective function needs to be calculated repeatedly. If the number of design variables is high, this gradient can be obtained efficiently from adjoint equations. This research focuses on the gradient calculation for an objective function which involves a fluid-structure interaction (FSI) simulation. The interaction can be calculated in a partitioned way by coupling a flow solver with a structural solver. In this work, quasi-Newton coupling iterations with an approximation for the inverse of the Jacobian from a least-squares model (IQN-ILS) are employed for the state equations as well as for their adjoint equations. The problem at hand is the unsteady, one-dimensional flow of an incompressible, inviscid fluid in an elastic tube. Special attention has been given to the interface variables which are exchanged between the adjoint flow and structural solver, to avoid the communication of system matrices between them.


INTRODUCTION
In this paper, unsteady adjoint equations are derived and solved for an elastic tube through which an incompressible fluid flows.The model for the tube wall is a one-dimensional generalised string model and the stiffness of each tube segment is determined by a separate parameter.The cost function measures the difference between the actual displacement of the tube wall and a prescribed displacement.
Several possibilities exist for the calculation of the cost function's gradient with respect to the stiffness parameters.It can be calculated using finite differences, the direct method and the adjoint method [1,29].The adjoint method is chosen because its computational cost does not significantly depend on the number of parameters, which is advantageous for a high number of parameters.The difference between the discrete and continuous adjoint formulation is that the discretisation occurs respectively before and after the derivation of the adjoint equations.Here, the discrete adjoint formulation is selected to obtain the exact gradient of the discrete equations.
The fluid-structure interaction in the model can be dealt with in two ways.In the monolithic approach, the flow equations and structural equations are solved simultaneously [15,17,12].Conversely, the flow equations and structural equations are solved separately in the partitioned approach.To obtain software modularity, the partitioned approach is selected for this research.In weakly (or explicitly, staggered) coupled partitioned simulations, the flow equations and the structural equations are solved once (or a limited, fixed number of times) per time step [10,8].As a result, the equilibrium of force and velocity on the fluid-structure interface is not guaranteed.In strongly (or implicitly) coupled partitioned simulations, coupling iterations are performed between the flow equations and the structural equations in each time step.Consequently, the equilibrium conditions on the interface are satisfied (up to an error proportional to the convergence tolerance of the coupling iterations).The weakly coupled methods are suitable for compressible aeroelastic simulations, while strong coupling should be used for simulations with incompressible fluids [35].
Several strongly coupled partitioned techniques treat the flow solver and the structural solver as black boxes, i.e. programmes which calculate an output for a given input, but whose internal algorithms can neither be accessed nor modified.These methods include Gauss-Seidel iterations, Gauss-Seidel iterations with Aitken relaxation [25,39,18], interface GMRES [24] and interface quasi-Newton (IQN-ILS) iterations [5,16].However, Gauss-Seidel iterations are often unstable if the ratio of the fluid density to the structure density is high, among other reasons [2,4].In comparisons consisting of several test cases, IQN-ILS requires fewer coupling iterations per time step than Aitken relaxation or Interface GMRES [5], and the computational cost ratio of a partitioned simulation with IQN-ILS to a monolithic simulation ranges from 0.56 to 3.16 [7].Therefore, IQN-ILS is selected as coupling algorithm.
In this work, both the forward and the adjoint fluid-structure interaction equations are solved in a partitioned way using the IQN-ILS coupling algorithm.A gradient for a spatially distributed parameter calculated using adjoint equations has already been combined with optimisation of fluid-structure interaction by several authors.For example, it has been applied for the aerostructural shape optimisation of an aeroplane [20,21] and its wings [22,23,9,11].However, these are steady problems which can be solved using Gauss-Seidel iterations, while Gauss-Seidel iterations are unstable for the unsteady model at hand.In [27,28], the adjoint of the steady Euler equations is calculated and the resulting gradient is inserted in an unsteady linear aeroelastic model with four equations to derive a controller for flutter reduction.
Unsteady adjoint solvers have been developed for the optimisation of fluids and structures alone.For example, an adjoint solver of the unsteady Euler equations has been used for the optimisation of blast mitigation devices [34].Also adjoint solvers of the unsteady Navier-Stokes equations exist [31], even for moving and deforming grids [26].Checkpointing is often applied to reduce the memory requirements of unsteady adjoint simulations [14,32,38], but this is not necessary in this case.
The remainder of this article is organised as follows.The equations and solution pro- cedure are described in Section 2 for the forward problem and in Section 3 for the adjoint problem.Section 4 presents the results, followed by the conclusions in Section 5.

Continuous equations
The unsteady flow of an incompressible, inviscid fluid in a straight, flexible tube is modelled (see Figure 1).The model is one-dimensional in an axisymmetric (r, φ, z) coordinate system.The non-overlapping domains of the fluid and the structure are indicated as Ω f and Ω s , respectively.The common boundary of these domains is the fluid-structure interface Γ, defined as with ∂Ω indicating the boundary of domain Ω.
The governing equations for the flow are the conservation of mass and momentum, formulated as with a the cross-sectional area of the tube, t the time, u the axial velocity and p the pressure.
For the structure, a so-called generalised string model is applied.This model is derived from the linear elasticity theory for a cylindrical tube with small thickness under the assumption of membrane deformations [30,13].It disregards the axial and circumferential displacement of the tube wall.The governing equation for the structure is given by with r the inner radius (a = πr 2 ), ρ s the structural density and h the thickness of the tube wall.The other coefficients are the shear modulus G, the Young modulus E and the Poisson coefficient ν.The viscoelastic term is further omitted (γ = 0), which is a common simplification [13].The Timoshenko shear correction factor κ is calculated from the Poisson coefficient

Discrete equations
The tube with length is discretised in space using m e segments with length ∆z.The pressure and velocity are stored in the cell centres.All terms in the flow equations are discretised using a central scheme, except for first-order upwind discretisation of the convective term in the momentum equation.The time discretisation is first-order backward Euler with time step size ∆t.The superscript n − 1 indicates the previous time level (t = (n − 1)∆t), while the superscript n for the current time level (t = n∆t) is omitted.
To obtain structural equations with only first derivatives in time, the radial velocity of the tube wall is introduced (v = ∂r/∂t).Although relatively uncommon, backward Euler time discretisation is applied for the structure as well, to avoid difficulties due to different time discretisation of the flow equations and the structural equations [36].Moreover, the spatial discretisation is performed using a finite difference scheme instead of the typical finite elements [19].
The elasticity modulus E m of each segment is determined by the corresponding parameter s m which varies from -1 to 1.
This definition of the elasticity modulus ensures that it can vary over a realistic range.The goal is to calculate the sensitivity of the cost function with respect to the value of the parameters s m (m ∈ {1, . . ., m e }) which appear in this equation.The discrete equations are subsequently linearised with respect to the reference values r o , p o , u o and v o .From this point on, r, p, u and v denote perturbations with respect to these reference values.Moreover, p o , u o and v o are set to zero to simplify the resulting equations.It has been demonstrated in previous research [6] that this particular choice of the reference values results in a model with the same numerical behaviour as the model with all nonlinear terms.Also physical properties such as wave propagation are preserved by this linearisation and reference.

Boundary conditions
As the segments of the tube are indicated with subscripts 1 to m e , the inlet (left-hand side) is indicated with a subscript 0 and the outlet (right-hand side) with a subscript m e + 1.
The flow rate at the inlet is prescribed as a function of the time t and the pressure is calculated using linear extrapolation.At the outlet of the tube, the velocity is obtained from a linear extrapolation and a Windkessel model relates this velocity with the outlet pressure [37].The parameter s me+1 determines the stiffness of this Windkessel model.The capacitor c represents the compliance, while the resistors r p and r d model the viscous resistance.The value of c is modified by the parameter s me+1 which also varies from -1 to 1.
For the structure, a zero-curvature boundary condition is applied at both the inlet and the outlet.

Matrix notation
The linearised equations for time step n can be written in the following block-matrix format In this equation, M f and M s are the system matrices of the flow solver and the structural solver, respectively.The interaction between the fluid and the structure is captured by the offdiagonal blocks C f and C s .The matrix C f describes how the residual of the flow equations changes due to a displacement of the fluid-structure interface, whereas the matrix C s converts the pressure on the interface into a contribution (nodal forces) to the residual of the structural equations.The state vector x is divided into a fluid part and a structure part, indicated with the subscripts f and s, respectively.Both x f and x s are subdivided once more using the subscripts Γ and Ω.The subscript Γ refers to variables on the fluid-structure interface, while the subscript Ω refers to variables that only appear in either the flow equations or the structural equations.
The right-hand side consists of time-dependent vectors b f and b s which do not depend on the state x, and a contribution which depends on the state in the previous time step x n−1 .The dimensions of x f Ω , x f Γ , x sΓ and x sΩ are respectively given by m f Ω × 1, m f Γ × 1, m sΓ ×1 and m sΩ ×1.The dimensions of the matrices in Equation 8 follow from the dimension of these vectors.For this specific case, the state vectors are given by with the superscript T indicating a transpose.The inlet and outlet model for the flow as described in Section 2.3 are included in m f Ω .The dimensions mentioned above can thus be written as a function of the number of tube segments m e .
Equation 8 can be abbreviated as with The initial state x 0 is set to zero.For the specific case that is analysed in this work, the block D s in this general equation is filled with zeros.Due to the linearisation, the matrices A and B are independent of x and therefore they remain constant during the simulation.The residual of the governing equations for time step n is then given by with

Coupling iterations
Equation 8 is solved in a partitioned way, which signifies that the flow equations and the structural equations are solved separately.Consequently, coupling iterations need to be performed between the flow equations and the structural equations to obtain the solution of the coupled problem.At convergence of the coupling iterations, the solution is identical to what a monolithic solver would calculate (up to an error proportional to the convergence tolerance of the coupling iterations).In every coupling iteration in time step n, the flow equations and the structural equations are solved with given values of x sΓ and x f Γ , respectively.In the following sections, the superscript k indicates coupling iteration k in time step n.
In the first coupling iteration, an extrapolation (E) based on previous time steps is performed to obtain the input for the flow solver.If no reuse from previous time steps is applied in the least-squares model, then a relaxation with factor ω determines the input for the flow solver in the second coupling iteration.At the beginning of coupling iteration k, the coupling algorithm calculates x k sΓ .The flow equations are then given by . (15) Once x k f Ω and x k f Γ have been calculated, x k f Γ is given to the structural solver.This calculation is further referred to as The structural solver subsequently solves the structural equations for xk sΓ and x k sΩ .The tilde is used to distinguish between the motion of the fluid-structure interface calculated by the coupling algorithm at the beginning of the coupling iteration and that calculated by the structural solver at the end.This calculation is further referred to as xk sΓ = S(x k f Γ ).The most trivial coupling algorithm is Gauss-Seidel iteration, which means that the motion of the fluid-structure interface calculated by the structural solver at the end of the previous coupling iteration is applied by the flow solver at the beginning of the next coupling iteration. x However, it is well-understood that this coupling algorithm is unstable for fluid-structure interaction problems with an incompressible fluid and comparable density of the fluid and the structure [2,6,4].Therefore, the interface quasi-Newton coupling algorithm with an approximation for the inverse of the Jacobian from a least-squares model (IQN-ILS) is applied [5].This coupling algorithm treats both the flow solver and the structural solver as black boxes.Using the interface displacement applied in the flow solver (x k sΓ ) and calculated by the structural solver (x k sΓ ) in all coupling iterations, the coupling algorithm constructs an approximation (indicated with a hat) for the inverse of the Jacobian of the interface residual with respect to the interface displacement x k sΓ .This results in the following update at the beginning of each coupling iteration Because only the product of the approximation for the inverse of the Jacobian with a vector is required, this matrix does not have to be constructed explicitly.In a matrix-free implementation of IQN-ILS, the computational cost and the memory requirements of the coupling algorithm scale linearly with the number of degrees of freedom in the interface displacement (m sΓ ).The least-squares model which predicts ∆x k−1 sΓ as a function of −r k−1 sΓ is further represented by M k .
Coupling iterations are performed until ||r k sΓ || 2 < c ||r 1 sΓ || 2 , in which c denotes the relative convergence criterion of the coupling iterations.The complete procedure for the forward simulation is described in Algorithm 1.

Algorithm 1
The IQN-ILS coupling algorithm for the forward simulation.for k = 1, . . ., k e do 3: else if k = 2 and no reuse then 6: else 8: end if 11: end for 18: end for 3. SENSITIVITY ANALYSIS

Cost function
The dimensionless least-squares cost function j is defined using a normalised difference between a prescribed displacement and a simulation of the model described above.It is given by with the superscript ref referring to the prescribed displacement (or reference).The vector s contains the m e + 1 parameters which are defined in Equation 5 and Equation 6.The vector x is the combination of the state vectors of all time steps The vector xsΓ is identical to x, but all entries which do not correspond to x n sΓ have been set to zero, giving So, the cost function is a sum over all time steps and all tube segments of the squared difference between the radius in the simulation and in the measurement.

Minimisation problem
With the definition of the cost function in Equation 20, a minimisation problem can be formulated as min s,x j(s, x) (23) with the governing equations as constraints The parameter vector s and the state vector x are defined in Equations 5, 6 and 21.Here, r is the combination of the residual of the governing equations in all time steps

Adjoint equations
As the state vector depends on the parameters, the gradient of the cost function j(s, x) = j(s, x(s)) requires application of the chain rule.The total derivatives of j with respect to the parameters are The partial derivatives in this equation can be calculated quickly and easily from Equation 20, as the state x remains constant for ∂j/∂s and the parameters s remain constant for ∂j/∂x.By contrast, the total derivative dx/ds requires the solution of the unsteady fluid-structure interaction problem and it is thus time-consuming to calculate.
As the governing equations of the fluid-structure interaction problem always need to be satisfied, r should always be equal to 0. Consequently, also the total derivative dr/ds should be equal to 0. By applying the chain rule, this total derivative is given by Rewriting the previous equation as ∂r ∂x dx ds = − ∂r ∂s (28) results in a system that can be solved for the total derivative dx/ds.The result is subsequently substituted in Equation 26, yielding ) −1 ∂r ∂s .(29) In this equation, the matrix inversion is a symbolic notation for the solution of a large system, which combines all equations of all time steps of the unsteady fluid-structure interaction problem.Obviously, this matrix is neither constructed explicitly nor inverted or factorised due to excessive memory requirements.As explained above, the time steps are calculated consecutively with a partitioned solution technique.Depending on how the factors are grouped in the last term of Equation 29, two methods can be distinguished.The direct method solves Equation 28 for dx/ds and substitutes the result in Equation 29.By contrast, the adjoint method first calculates the product of the first two factors in the last term of Equation 29.

∂j ∂x ( ∂r ∂x
) The vector a is the so-called adjoint (or dual) state, which is the solution of the adjoint equation ( ∂r ∂x The sought-after gradient dj/ds is finally calculated as The computational cost of solving Equation 31 is independent of the number of parameters, so the adjoint method is suitable for a large number of parameters.However, there are as many right-hand sides as objective functions.As for the direct method, the matrix in Equation 31is too large to be factorised and stored in the case of unsteady fluid-structure interaction.So, the adjoint equations need to be solved for every objective function.Alternatively, the adjoint method can be derived by introducing Lagrange multipliers, for example as in [33].From the explanation above, it is clear that the choice between the direct and the adjoint method depends on the number of parameters and objective functions.The direct method requires the solution of a problem comparable to the governing equations for every parameter, while its cost is almost independent of the number of objective functions.The opposite is true for the adjoint method.If central (respectively forward) finite differences are used for the calculation of the gradient, the governing equations need to be solved twice (respectively once) for every parameter and for every objective function.In this specific case, there is only one objective function and there are many parameters, so the adjoint method is selected.
In this case, the cost function (Equation 20) does not explicitly depend on the parameters s, so ∂j ∂s = 0.
Conversely, the residual r (Equation 14) depends on the parameters s, giving The only non-zero contributions to ∂r/∂s are due to the terms that contain the elasticity modulus E m (Equation 5) and the compliance of the Windkessel model c (Equation 6).
In the adjoint equation (Equation 31), both (∂j/∂x) T and (∂r/∂x) T are required.The former is calculated analytically, giving The calculation of ∂r/∂x and its transpose are more involved.This matrix is not constructed or stored, but the solution to Equation 31 is calculated using time steps.Using Equation 14, ∂r/∂x is given by Because ∂r/∂x is lower bidiagonal, the forward simulation consists of forward time steps.By contrast, its transpose is upper bidiagonal, which is the reason for the backward time steps in the adjoint simulation.

Matrix notation
Considering Equation 38, the block-matrix format of a time step for the solution of the adjoint equation (Equation 31) is similar to Equation 8.
This equation is written in a general form.The structure is self-adjoint which allows for a simplification (M T s = M s ), but this is not used here to keep the formulation general.In abbreviated form, the previous equation yields which is similar to Equation 11, except for the backward time steps (n ∈ {n e , . . ., 1}).The initial adjoint state a ne+1 is again set to zero.

Coupling iterations
Equation 39 is solved in a partitioned way as well, by performing coupling iterations between the adjoint flow equations and the adjoint structural equations until the solution of the coupled adjoint problem has been found.The superscript k again indicates coupling iteration k in time step n.
The adjoint flow equations in coupling iteration k are given by (41) which are solved for a k f Ω and a k f Γ .This equation contains C s and D s , which belong to the structural solver.The exchange of matrices between the flow solver and the structural solver is highly unwanted in the partitioned approach.Therefore, C T s a k sΓ and D T s a n+1 sΓ are given to the flow solver, instead of a k sΓ and a n+1 sΓ .While C T s a k sΓ is provided to the flow solver in every coupling iteration, D T s a n+1 sΓ is only exchanged once at the beginning of every time step.In the special case that C s and D s are identical, giving D T s a n+1 sΓ to the flow solver is not required because the flow solver can store the information from t n+1 .
The structural solver subsequently solves the adjoint structural equations ) for ãk sΓ and a k sΩ .The tilde is used to distinguish between the adjoint state calculated by the coupling algorithm at the beginning of the coupling iteration and that calculated by the structural solver at the end.This equation contains C f and D f , which belong to the flow solver.To avoid the exchange of matrices, C T f a k f Γ and D T f a n+1 f Γ are provided to the structural solver.While the former is given in every coupling iteration, the latter is only given once at the beginning of the time step.
In conclusion, the flow solver calculates a k f for a given C T s a k sΓ and gives C T f a k f Γ to the structural solver, whereas the structural solver calculates a k s for a given C T f a k f Γ and returns C T s ãk sΓ .These calculations are further referred to as 41) and structural equations (Equation 42) are coupled using the IQN-ILS algorithm, similarly to the forward equations.The partitioned adjoint simulation is described in Algorithm 2. The differences with the partitioned forward simulation are as follows.First, the vectors D T f a n+1 f Γ and D T s a n+1 sΓ are exchanged between the flow solver and the structural solver at the beginning of every time step (lines 2-3).Then, the extrapolation (E) is adapted in a straightforward way to the backward time steps (line 6).Moreover, the interface residual is defined as for k = 1, . . ., k e do 5:

The adjoint flow equations (Equation
else if k = 2 and no reuse then 8: Equation 41 end for 20: end for

RESULTS
For the results, the tube is discretised in m e = 100 segments.The parameters of the fluid-structure interaction model and the Windkessel model are listed in Table 1.
Figure 2 depicts the radius at the middle of the tube during one period for the minimal (s = −1), nominal (s = 0) and maximal (s = 1) values of all parameters.The differences in the other variables and at other locations are of the same order of magnitude.The cost function's gradient for s = 0 is shown in this figure, with the reference in Equation 20 calculated using a sinusoidal variation of the tube's stiffness.
The adjoint calculation of the cost function's gradient is verified by comparing it to central finite differences ( dj ds m Table 1.The parameters of the fluid-structure interaction model and the Windkessel model [37].with s m and ∆s m respectively the value and the perturbation of element m in vector s.This comparison is performed for different values of m and s, with ∆s m = 10 −3 and 10 −4 .Table 2 lists the values from this comparison between the adjoint calculation and the finite difference calculation with ∆s m = 10 −4 .All parameters in the reference calculation (required for the cost function calculation in Equation 20) have been set to 1.For the forward and the adjoint calculation, the elements of the vector s have consecutively been set to -1, 0 and 1. Elements m ∈ {1, 10, 101} of the gradient are analysed.The difference between both gradients is at least 5 orders of magnitude smaller than the absolute value of the gradient.

CONCLUSIONS
In this work, the cost function is calculated with an unsteady fluid-structure interaction simulation and its gradient is obtained from an unsteady adjoint simulation.Both the forward and the adjoint simulation are partitioned with a quasi-Newton coupling algorithm (IQN-ILS).The exchange of operators between the flow solver and the structural solver is avoided by a particular choice of the variables that are exchanged between the solvers.Table 2.The verification of the gradient calculation by means of a comparison with finite differences using step size ∆s m = 10

Figure 1 .
Figure 1.The model for the flow in a straight, flexible tube with details of the cross-section and a control volume used in the discretisation of the governing equations.Also the prescribed velocity at the inlet and the Windkessel model at the outlet are depicted.

Algorithm 2
) in the adjoint simulation.Finally, xi sΓ is replaced by C T s ãi sΓ in the least-squares model (M on line 10).The IQN-ILS coupling algorithm for the adjoint simulation.1: for n = n e , . . ., 1 do 2: send D T f a n+1 f Γ from flow solver to structural solver 3: send D T s a n+1 sΓ from structural solver to flow solver 4: c o 6.35•10 −10 m 3 /Pa r d 1.768•10 9 Pa•s/m 3 r p 2.834•10 8 Pa•s/m 3

Figure 2 .
Figure 2. (left) The radius linearised with respect to r o at the middle of the tube (z = /2) as a function of time for minimal (s = −1), nominal (s = 0) and maximal values (s = 1) of all parameters.(right) The cost function's gradient with respect to the stiffness of the tube segments for s = 0 if the reference is calculated with a sinusoidal variation of the stiffness.