A Multi-Vector Interface Quasi-Newton Method with Linear Complexity for Partitioned Fluid-Structure Interaction

The response of the structure to external loads is expressed via the displacement field $`\vekt{d}^s(\vekt{x},t)`$, which is governed by the equation of motion, stating a dynamic balance of inner and outer stresses:

\begin{alignat}
{2}
    \rho^s \frac{d^2 \vekt{d}^s}{dt^2} &= \boldsymbol{\nabla} \cdot \CauchyStress^s + \vekt{b}^s \qquad &&\text{in } \Omega_t^s ~~\forall t \in [0,T]~.
\end{alignat}

In this relation, $`\rho^s`$ denotes the material density and $`\vekt{b}^s`$ the resultant of all body forces per unit volume, whereas $`\CauchyStress^s`$ represents the Cauchy stress tensor.

As constitutive relation, the St. Venant-Kirchhoff material model is used: It relates the 2nd Piola-Kirchhoff stresses $`\vekt{S} := \det (\vekt{F}) ~ \vekt{F}^{-1} \, \vekt{T}^s \, \vekt{F}^{-T}`$ to the Green-Lagrange strains $`\vekt{E} := \frac{1}{2} \left( \vekt{F}^T \vekt{F} - \vekt{I} \right)`$ via a linear stress-strain law, reading

\begin{align}
\vekt{S} = \lambda^s \text{tr} \left( \vekt{E} \right) + 2 \mu^s \vekt{E} ~.
\end{align}

Therein, $`\vekt{F}`$ denotes the deformation gradient, while $`\lambda^s`$ and $`\mu^s`$ are the Lamé constants. As the Green-Lagrange strain definition forms a nonlinear kinematic relation, the structural model is geometrically nonlinear. Hence, it is capable of representing large displacements and rotations, but only small strains .

Collecting all this information, the equation of motion can be expressed in the (undeformed) reference configuration $`\Omega^s_0`$ in a total Lagrangian fashion:

\begin{alignat}
{2}
\rho^s \frac{d^2 \vekt{d}^s}{dt^2} &= \boldsymbol{\nabla}_0 \cdot \left( \vekt{S} \vekt{F}^T \right) + \vekt{b}^s \qquad &&\text{in } \Omega_0^s ~~\forall t \in [0,T]~.
\end{alignat}

Again, the problem is closed by defining an initial displacement $`\vekt{d}_0`$, which is typically zero, and a set of boundary conditions on two complementary subsets of $`\Gamma_0^s = \partial \Omega^s_0`$: Prescribing the displacement $`\vekt{g}^s`$ on the Dirichlet part $`\Gamma_{D,0}^s`$ and the tractions $`\vekt{h}^s`$ on the Neumann part $`\Gamma_{N,0}^s`$, with the outer normal in the reference state $`\vekt{n}_0`$, the conditions read

\begin{alignat}{2}
    \vekt{d}^s(\vekt{x},t=0) =& ~\vekt{d}_0^s(\vekt{x} )\qquad &&\text{in } \Omega_0^s ~, \\
    \vekt{d}^s =& ~ \vekt{g}^s       &&\text{on }  \Gamma_{D,\,0}^{s} ~~\forall t \in [0,T] ~, \\
    \vekt{F \,S} \, \vekt{n_0} =& ~ \vekt{h}^s       &&\text{on }  \Gamma_{N,\,0}^{s} ~~\forall t \in [0,T] ~.
    \end{alignat}

We gratefully acknowledge the computing time granted by the Jülich-Aachen Research Alliance JARA-HPC.

In general, fluid-structure interaction considers a fluid domain $`\Omega_t^f \subset \mathbb{R}^{\nsd}`$ and a structure $`\Omega_t^s \subset \mathbb{R}^{\nsd}`$, that are connected via a shared boundary $`\Gamma_t^{fs} = \partial \Omega_t^f \cap \partial \Omega_t^s \subset \mathbb{R}^{\nsd-1}`$, the FSI interface. The subscript $`t`$ refers to the time level, while $`\nsd`$ denotes the number of spatial dimensions. This section introduces the models and equations employed for the two subproblems along with the boundary conditions interlinking their solutions at the shared interface.

This paper presents a novel interface quasi-Newton method improving the partitioned simulation of fluid-structure interaction. Building on the inverse multi-vector Jacobian (IQN-IMVJ) variant, it overcomes the main drawback: As the cost of handling an explicit inverse Jacobian approximation increases quadratically with the problem size, the IQN-IMVJ method significantly slows down large-scale simulations. The new interface quasi-Newton implicit multi-vector least-squares (IQN-IMVLS) method, in contrast, systematically avoids any explicit Jacobian approximation, so that the multi-vector update is realized with linear complexity. As a consequence, its computational cost stays negligible even and particularly for industrial-scale applications.

The effectiveness of the IQN-IMVLS method is confirmed by the numerical test cases in Section 5. Although the number of required coupling iterations stays virtually unchanged compared to the IQN-IMVJ variant, the results show a substantial reduction of computational cost. In accordance to the discussion above, this difference becomes more and more distinct for finer resolutions of the structural discretization. One major strength of the multi-vector concept is that its implicit incorporation of past data does not rely on problem-dependent parameters – and yet keeps up with an optimal parameter choice of the interface quasi-Newton implicit least-squares (IQN-ILS) method, regarding the required coupling iterations. While the IQN-IMVJ approach had to pay a high price in computational cost for this capability, the IQN-IMVLS method is comparable to the IQN-ILS variant concerning the runtime.

Although the IQN-IMVLS variant reintroduces the number of reused time steps $`q`$, this parameter is not problematic at all as the quality of the quasi-Newton update in general benefits from including more steps.
In conclusion, the new IQN-IMVLS method combines the advantages of the IQN-IMVJ and the IQN-ILS variant without their individual drawbacks. Coping without critical parameters, it succeeds in realizing the favorable implicit reutilization of past data with a computational complexity that grows linearly with the problem size.

In recent years, interface quasi-Newton methods have gained growing attention in the fluid-structure interaction community by significantly improving partitioned solution schemes: They not only help to control the inherent added-mass instability, but also prove to substantially speed up the coupling’s convergence. In this work, we present a novel variant: The key idea is to build on the multi-vector Jacobian update scheme first presented by Bogaers et al. and avoid any explicit representation of the (inverse) Jacobian approximation, since it slows down the solution for large systems. Instead, all terms involving a quadratic complexity have been systematically eliminated. The result is a new multi-vector interface quasi-Newton variant whose computational cost scales linearly with the problem size.

The second example is a cylindrical deformable tank partially filled with an incompressible fluid. While Figure [fig:SloshinTank] illustrates the basic configuration, the test case parameters are listed in Table [tab:SloshingTank].

The simulation is started from the steady state where the fluid is at rest and the deformation has adjusted to the hydrostatic pressure. During the simulation, the system is excited in horizontal direction by a periodic movement of the tank bottom with the prescribed displacement $`\Delta x(t) = A \, \sin (2 \pi f \,t)`$. Due to inertia, this excitation results in complex and highly intertwined motions of the deforming tank structure and the liquid sloshing inside.

The evolving fluid domain is discretized by a space-time mesh of $`16243`$ elements and $`4024`$ nodes per time level; since it is adapted to the free-surface via interface tracking , the tank wall and its bottom involve slip-conditions for both the fluid and the mesh. The cylindrical tank wall is modeled by $`32 \times 30 =960`$ nonlinear isogeometric shell elements of degree $`2`$ defined by a NURBS with $`1184`$ control points. The simulation runs for $`1600`$ time steps of size $`\Delta t=0.01`$s; the absolute and relative convergence criteria for the interface deformation are $`\varepsilon_{abs}=10^{-5}`$ and $`\varepsilon_{rel}=10^{-4}`$.

Results

To give an impression of the simulation results, Figure 1 depicts three different snapshots of the sloshing tank. Both the free surface and the tank structure exhibit large displacements.

The efficiency of different update techniques is compared in Table 1, both in terms of coupling iterations and their computational cost. All in all, the results are in very good agreement with those of the elastic tube scenario: While neither a constant nor Aitken’s dynamic relaxation converge fast enough to be a reasonable option, all interface quasi-Newton variants drastically speed up the coupling scheme. Again, incorporating information from past time steps is extremely beneficial. The IQN-ILS method performs very well for this test case, showing a quite flat optimum of choosing the parameter $`q`$; the cost of solving the least-squares problem, however, is clearly increasing with $`q`$.

Although being slightly less efficient for this test case, the implicit reutilization of the multi-vector methods almost keeps up with the IQN-ILS approach in terms of the required coupling iterations. While the IQN-IMVJ variant pays for this fast convergence with the expensive explicit Jacobian, the computational cost of both the IQN-IMVLS and the IQN-MVJ-RS-SVD method is comparable to that of the IQN-ILS method, thanks to the linear complexity. In direct comparison, the new IQN-IMVLS method outruns the restart-based variant for this test case; however, this statement is again best put into the perspective of the specific implementation and the chosen parameters.

Increasing the value $`q`$ brings the convergence speed closer to that of the explicit multi-vector approach. Although a rather high value is required here to reach this speed, the implicit Jacobian product proves to be very cheap as the computational cost stays non-critical even for the maximum of reusing all time steps, i.e., $`q=1600`$. By including the data pairs of the most recent time step explicitly in the least-squares problem, for this test case the convergence of the IQN-IMVLS method is accelerated by about one coupling iteration per time step at moderate cost.

Altogether, the test case proves that the IQN-IMVLS method competes with the IQN-ILS approach not only in terms of the convergence speed, but also regarding the computational cost.