Efficiency of fluid-structure interaction simulations with adaptive underrelaxation and multigrid acceleration

In the present paper the efficiency of acceleration techniques for fluidstructure interaction computations are investigated. The solution procedure involves the finite volume flow solver FASTEST, the finiteelement structural solver FEAP, and the coupling interface MpCCI. Within the employed partitioned solution approach, a geometric multigrid solution strategy on moving grids for the fluid domain is introduced. In particular, the order in which the convective fluxes have to be treated within the pressurecorrection smoothing procedure is addressed. For reducing the coupling iteration steps an adaptive underrelaxaation algorithm is employed. Both acceleration techniques are investigated separately and in combination with respect to numerical efficiency. As test configuration a representative three-dimensional ullsteady coupled problem is considered.


INTRODUCTION
Coupled fluid solid problems usually require a high computational effort.Especially in threedimensional cases the computation time for the fluid part often increases dramatically with the number of unknowns when applying simple iterative solution algorithms.
Within the frame of a partitioned approach, on the one hand, both solution algorithms for the fluid and structure problems may be irmproved and, on the other hand, the coupling strategy can be optirnized.
In the present investigation a multigrid procedure for moving meshes together with an adaptive underrelaxation strategy is applied for accelerating the coupled computations.A few investigations concerning the efficiency of adaptive coupling schellles have been presented in earlier works.In [9] a comparison between three different methods, the Aitken, the Tschebyscheff, and the method of deepest descent is presented.In [15] modified Aitken-like methods are investigated including a study of their speed-up behavior.The works concentrate on the investigation of the adaptive coupling procedure.The present paper extends these studies by employing additionally a multigrid method and investigating the combined acceleration.
Concerning multigrid techniques a huge numbers of works have been published.An introduction is given in [2] also describing non-linear cases.In [16] multigrid techniques for the incompressible Navier-Stokes equations in combination with finite volume discretizations are presented.However, detailed descriptions of multigrid techniques for fluid computations on moving meshes as they are needed for coupled fluid-structure problems can hardly be found in the literature.In the present paper a modified multigrid algorithm is outlined especially taking into account the correct treatmeIlt of convective fluxes in case of moving grids.
The investigations are performed on the basis of a three-dimeIlsional unsteady fluidstructure interaction configuration.

GOVERNING EQUATIONS
We consider a problem domain Ω consisting of a fluid part Ω f and a solid part Ω s .For the fluid domain part Ω f we assume a flow of an incompressible Newtonian fluid.In this case the basic conservation equations governing transport of mass and momentlml for a fluid control vollmle V f with surface S f are given by (1) (2) where v is the velocit,y vector with respect to Cartesian coordinates x, t is the time, ρ f iS the fluid density, n is the outward normal vector and f f are external volume forces (e.g., buoyancy forces).v 9 is the velocity with which S f may move (grid velocity) due to displacements of solid parts.The Cauchy stress tensor T f for incompressible Newtonian fluids is defined by (3) with the pressure p,s the dynamic viscosity µ f and the identity tensor I.For the structure we denote a material point in the reference configuration as X whose position in the current configuration is given by (4) The displacements are evaluated by (5) For more details see [14,10].The basic balance equation for momentum for the solid u = x − X .
domain Ω s can he written as (6) where x ¨= ∂ 2 x(X, t)/∂t 2 is the acceleration, S s denotes the second PiolaKirchhoff stress tensor, ρ s is the density of the solid, and f s are external volume forces acting on the solid (e.g., gravitational forces).F s = ∂ x /∂X denotes the deformation gradient.
In the present investigation we consider the Saint Venant-Kirchhoff material law (7) with the Green-Lagrangian strain tensor (8) as kinematic property.λ s and µ s are the two Larmé constants.
The problem formulation has to be closed by prescribing suitable boundary and interface conditions.On solid and fluid boundaries Γ S and Γ f standard conditions as for individual solid and fluid problems can be prescribed.For the velocities and the stresses on a fluid-solid interface Γ i we have the conditions (9) where v b is the velocity of the interface and T s = F s S s F s T s T / det F s is the Cauchy stress tensor.

NUMERICAL FLUID-STRUCTURE COUPLING SCHEME
The discretization of the problem domain is based on a block-structuring technique.Fluid and solid parts are assigned to different blocks.Solid blocks are treated by the finite-element solver FEAP (see [13]).For the fluid blocks, which can be defined as moving or fixed, the parallel multigrid finite-volume flow solver FASTEST is employed (see [5,12]).Both solvers involve second-order spatial discretizations and fully implicit second-order time discretizations.
For the fluid-structure coupling an implicit partitioned approach is employed.In Fig. 1 a schematic view of the iteration process is given.After the initializations the flow field is determined in the actual flow geometry.From this the friction and pressure forces on the interacting walls are computed.These are passed to the structural solver as bolmdary conditions.The structural solver computes the deformations, with which then the fluid mesh is modified.Afterwards the flow solver is started again.The fluid-structure interaction (FSI) iteration loop is repeated until a convergence criterion is reached, which is defined by the change of the mean displacements (10) Where m is the FSI iteration counter, N is the number of interface nodes, and ʈ .ʈ ∞ denotes the maximum norm.Note that an explicit coupling method would be obtained, if only one FSI iteration is performed.
The data transfer between the flow and solid solvers within the partitioned solution procedure is performed via an interface realized by the coupling library MpCCI (see [6]) that controls the data communication and also carries out the interpolations of the data from the fluid and solid grids.More details are given in [11].
Various test computations have shown that the coupling scheme is rather sensitive with respect to the deformations especially in the first FSI iterations.Here, situations that are far away from the physical equilibrium can arise, which may lead to instabilities or even the divergence of the FSI iterations.In order to counteract this effect an adaptive underrelaxation is employed.
By using an relaxation factor α m FSI the actually computed displacements u ~ m are linearly weighted with the values u m-1 from the preceding iteration to give the new displacements u m+1 :  (11) where 0 < α m FSI ≤ 1.Note that the underrelaxation does not change the final converged solution.The basic effects of the underrelaxation has already been shown in [11].
For the adaptive determination of α m FSI different methods are known.We employ here an approach based on the Aitken method for vectorial equations, which is an extrapolation approach frequently applied in the context of Newton-Raphson iterations.The basis of this approach was proposed by Aitken in 1937 [1] and later improved by Irons and Tuck [8].It was identified as very efficient for computations in the field of fluid-structure interaction by Mok [9].
Employing the vahles from two preceding iterations the so-called Aitken factor γ m is extrapolated by: ( 12) The actual underrelaxation factor α m FSI then is defined by (14) As first Aitken factor in each time step γ 0 the last one from the preceding time step can be taken.For the first time step some reasonable value can be chosen, e.g., γ 0 = 0.

MULTIGRID METHODS ON MOVING GRIDS
For a better understanding of the multigrid procedure on moving meshes with all the aspects that have to be considered, first, the iteration algorithm for coupled computations on single grids will be described.Afterwards the multigrid method without mesh moving is introduced, and finally, the extension to multigrid in combiIlatioIl with moving grids is presented.

SINGLE GRID METHOD ON MOVING MESHES
Considering the momentum equation ( 2), an additional term for moving meshes (15) appears.This can be seen as a flux correction term due to the grid movement that has to be subtracted from the usual convective flux F = in case of fixed meshes.As suggested in [3] and [4], involving the space conservation law Thus, if one ensures that the space conservation law is fulfilled at any time the mass conservation equation can be remained unchanged and the only difference between fixed and moving meshes is the flux correction term F corr in the momentum equation.
Let us denote the routine which assembles the iteration matrix for the momentum equation as "BUILD-A".BUILD-A takes the flux F from the previous iteration and assembles the matrix coefficients for the convective part (note that fluxes are considered in all three spatial dimensions and are denoted as a vector from now on).Now, we first modify the flux F by subtracting F corr so that we can use the new corrected flux F * = F -F corr to use the original BUILD_A for fixed meshes.This flux correction plays a very important role for applying the multigrid method to moving meshes.
The following Algolithm 1 shows the solution procedure for the coupled computation on single grids:

MULTIGRID METHODS 0N NON-MOVING MESHES
Employing a simple iterative scheme on a single grid, it turns out, that the high frequencies of the error are reduced very quickly but the low frequencies decrease gradually.The idea of multigrid acceleration is to build an interpolated equation for the error on a coarser grid.Regarding the larger grid spacing the lower error frequencies from the fine grid becomes higher frequencies and can be reduced much faster due to the smaller number of control volumes.For further details of multigrid techniques see [2,16,7].
We consider the non-linear system of equations on grid level k with the coefficient matrix A k and the source term S k , both depending on the solution Φ k = (p,v.F) k .Starting on grid level k with φ 0 k from the previous solution or initialization, after n iterations the variables satisfy the equations only up to a residual R n k (19) Subtracting equation ( 19) From (18) leads to (20) that can be seen as an equation for corrections to reduce the residual term on the right hand side and solved therefore on a coarser grid k -1 with double grid spacing (the nonlinear dependency of A and S now is omitted for simplicity).The assembly of the coarse grid equation begins with transferring the residual R k n to R const k-1 by summing up the corresponding fine grid residuals.The same procedure is applied to the fluxes.With an appropriate interpolation operator I k k-1 the fine grid variables φ n k = (P,υ) n k are restricted to the coarse grid With these restricted properties the coarse grid equation can be written as (21) Α const k-1 and S const k-1 are computed in the same way as on the fine grid using the restricted fine grid variables.While the underlined term in (21), remain unchanged during the iteration process, the other terms change.Only in the very first iteration φ 0 k-1 iS equal to φ const k-1 .After m coarse grid iterations we obtain If the underlined term of equation ( 21) is now joined together with S m k-1 and R m k-1 is subtracted, the resulting equation has the same appearance as (19) and the whole process can be recursively repeated for coarser grid levels k-2, etc.
The corrections are evaluated as (22) and prolongated to the finer grid using an appropriate interpolation operator: (23) ψ k is added to φ n k resulting in an improved solution φ impr k .This can be used either to do some damping iterations necessary due to interpolation errors, or to begin with a new V-cycle when the finest grid is reached.
The following Algorithm 2 summarizes the multigrid solution procedure for fixed meshes (for simplicity only for 2 grid levels):

MULTIGRID METHODS ON MOVING MES~E.S
Now the single grid approach for the flux correction for moving grids as described in section 4.1 is applied to the multigrid solution strategy for fixed meshes as introduced in the previous section.We employ the fact that the only difference between fixed and moving meshes is the flux correction F * = F -F corr , provided that the space conservation law is fulfilled at any time.
As in case of moving single grids, (p, v, F) 0 k is given from the previous iteration or initialization.Notice, that the routine BUILD A which assembles the momentum iteration matrix A is not modified for the extension to moving grids.Therefore, first the correction of fluxes F * = F -F corr has to be performed and afterwards the discrete equation can be built.The remaining solution procedure on the fine grid is unchanged, i.e., solve momentum equation, solve pressure-correction equation, and update variables.After n iterations the preliminary solution (p, v, F) n k is found.These quantities, and additionally R n k , have to be restricted to the coarser grid.The routine BUILD-A is used to compute Rk. but since BUILD A has not been modified for moving grids, once again the flux correction has to be applied to calculate the correct residual on grid k employing BUILD A. Doing so, the flux F n k , which also needs to be restricted, would have been changed.Thus, before flux correction F n k must be saved, then the flux is corrected, R n k is calculated with BUILD A, and the saved flux is transfered back.Now all quantities of the fine grid (p, v, F, R) n k , may be restricted to the coarser grid k-1 in the same way as described in the previous section.
Since only corrections are computed on the coarser grid, flux correction has not to be applied there.The rermaining solution steps, i.e., solve momentuu equation, solve pressure correction, calculate corrections, and prolongate corrections are identical to the case of fixed grids.After the improved solution (p, v, F) imp k is achieved a new V-cycle can be started.
The following Algorithm 3 indicates the whole solution procedure for the multigrid approach on moving meshes for coupled cormputations:

TEST CASE
As test case a three dimensional lid driven cubic cavity is considered (see Figure 2).The side length of the cube is L = 1m.At the bottom a flexible membrane with thickness t mem = 0.1 m is situated.The membrane is fixed at the edges and the lid moves with a time dependent lid velocity v lid defined by v lid (t) = 0.5 (1.5 -cos (2πt/T 0 )) with period T 0 = 5s.The inlet and outlet have 10% of the total height h in = h out =0.1 m and are situated in the upper part on the left and right sides, respectively.The inflow velocity corresponds to the lid velocity.For the investigation three successive refined fluid grids are considered.The multigrid computations are done with three grid levels for each grid set.The finest grid level of the multigrid computation also is taken as discretization for the single grid computation.The settings with the corresponding numbers of control volumes are summarized in Table 1.Further, two different types of structural underrelaxation are analyzed, i.e., the standard one with fixed α m FSI = 0.8 and the adaptive Aitken method as outlined in Section 3 for the medium grid set.
The structure discretization is fixed for all computations.Here, linear solid brick elements are used with a discretization of 20 × 2 × 20 for the x-, y-, and z-directions, respectively.To give an illustration of the problem solution in Fig. 3 a snapshot of the velocity field and in Fig. 4 the structural deformations for one period (in steps of 45 º time phase angles) are shown.In Fig. 5 the temporal displacement evolution of the center of the membrane is depicted.The point denoted as start point of measure indicates the initial point from where on the computing time measurements are taken.In each case a simulation time of 1.5s from a predefined flow field computed.The computer employed for all computations was a Intel Pentium4 PC with a clock rate of 2539 MHz.

RESULTS
The objectivc of this investigation is to quantify the accelerations that can be achieved by the adaptivc underrelaxation and the multigrid technique.For this the averaged computing times, numbers of fine grid iterations, and numbers of coupling steps per time steps are compsred for the different setups.
At first.single grid and multigrid computation with fixed underrelaxation parameter α m FSI = 0.8 are considered.In Table 2 the numbers of fine grid iterations and the computing times for the different grid sets are summarized together with the corresponding acceleration factors.One can observe that even on the coarsest grid the multigrid method already is nearly two times faster as the single grid method.The acceleration factor increases to more than 12 for the fine grid.However, the finer the grid the larger is the portion of time that is spent on the coarser grid levels.In Fig. 6 the accelerating effect of the multigrid approach is depicted graphically.It shows the coursc of the residuals of the mass during one time step for the single grid and the multigrid method for the finest grid.The accelerating effect of the multigrid method clearly call be seen.Next, the standard and adaptive underrelaxations are compared.The average numbers of FSI iterations per time step and the computing times for single grid and multigrid methods start point of measure and the corresponding acceleration factors for the medium grid set are indicated in Table 3.
The average numbers of coupling steps per time step hardly vary between the single grid and multigrid methods.In both cases the number is nearly halved when using the adaptive Aitken method.The reduction of the number of FSI iterations also has a direct effect on the total computation time.By using the method with the Aitken accelerator one can reduce the computation time by nearly a factor of 2 in both cases.Combining both approaches the computation with multigrid and adaptive Imderrelaxation is more than 14 times faster than the single grid computation with fixed underrelaxation.

CONCLUSION
An adaptive underrelaxation approach and a multigrid algorithm on moving meshes have been studied with respect to their capabilities to reduce the cormputing tirmes for fully coupled fluid-structure interaction computations.In the nonlinear multigrid method special care has to be taken to the consistent treatment of the additional fluxes due to the grid movement.If this is done properly, the method shows the typical behavior with acceleration factors increasing with the grid size also known from uncoupled simulations.Acceleration factors up to 12 for the finest grid could be achieved for a typical fully coupled three-dimensional FSI simulation.
An improvement also can be achieved by employing an adaptive underrelaxa,tion based on the Aitken method.The acceleration effect is similar for the single grid and multigrid methods.For both an acceleration factor of nearly 2 could be achieved for the considered test case.Since both acceleration techniques hardly effect each other, they can be combined to yield the best performance.Already for the medium sized grid an acceleration factor larger than 14 could be obtained this way.

Figure 1
Figure 1 Flow cllart of coupled solution procedure.
(16) ∂ ∂t the flux correction is discretized by (17) ensuring mass conservativity.The quantities δV c n are called swept volurmes and can easily be calculated out of the change of volurmes due to the rmesh movement.The summation index c runs over the faces of the control volume, the index n denotes the time level tn, and ∆tn is the time step size.

Algorithm 2 ( 1 )
Pressure p 0 k , velocities υ 0 k and fluxes F 0 k are given from the previous iteration or initialization.(2) Do n iterations (2a)-(2c) on grid k: (a) Build discrete momentum equations (BUILD A) on grid k and solve for υ`k employing SIMPLE-Algorithm.(b) Solve pressure-correction equation and update (p, υ, F), k (c) Go to (2a).(3) Perform restriction of (p, υ, F, R), k n : (a) Compute R k n on grid k using routine BUILD-A.(b) R 0

Algorithm 3 ( 1 )
Pressure p 0 k , velocities v 0 k and fluxes F 0 k are given from the previous iteration or initialization.(2) Do n iterations (2a)-(2d) on grid k: (a) Correction of fluxes: F k * = F k -F corr k .(b) Build discrete momentum equations (BUILD-A) on grid k and solve for v ´k employing SIMPLE-Algorithm.(c) Solve pressure-correction equation and update (p, v, F) k .(d) Go to (2a).(3) Perform restriction of (p, v, F, R) n k : (a) Save fluxes: ) Assemble coarse grid equations with restricted values.(5) Do m iterations (5a)-(5c) on grid k-1: (a) Solve momentum equation on grid k-1.(b) Solve pressure-correction equation on grid k-1 and update (p, v, F) k-1 .(c) Go to (5a) .(6) Compute corrections: ψ k-l ≈ (p, v) m k-1 -(p,v) 0 k-1 .(7) Prolongate corrections to fine grid ψ k = I k-1 k ψ k-l and update (p, v, F) imp k .. (8) With the improved fluid field properties on grid k a new V-cycle can be started (go to 1) until convergence is reached for the fluid field properties.(9) Calculate wall forces and send to the structural solver.(10) Receive displacements from the structural solver.(11) Compute new grid and swept volumes δV n c .(12) Go to 1 until FSI-convergence is achieved.(13) Start new time step.

Figure 2
Figure 2 Test case configuration

Figure 3
Figure 3 Snapshot of velocity field for test case.

Figure 4
Figure 4  Structural deformations for one period in steps of 45°time phase angles.

Figure 6
Figure 6  Comparison of residuals of mass conservation for the same whole time step (finest grid levdl).

Table I
Computational gridb and multigrid levelr.

Table 2
Figure 5 Displacement of the center of the membrane in y-direction over time, and start point of measure.

Table 3
Average number of FSI Iterations and computing times for single grid (SG) and multigrid (MG) and acceleration with fixed and adaptive underrelaxation for medium grid set.