My fourth tutorial: 3dTube
Tutorial Aims
 Demonstrate the solution of an internal flow fluidsolid interaction problem;
 Demonstrate the fluidsolid interaction coupling approaches available in solids4foam;
 Provide insight into the relative performance of the different coupling approaches for internal flows;
 Examine the effect of timestep on the performance of the coupling approaches.
Case Overview
This case consists of a pressure pulse applied in a thickwalled elastic tube (Figure 1).
Figure 1: Wave propagation in an elastic pipe
The fluid is assumed incompressible, Newtonian and isothermal, with a density of 1000.0 kg/m3 and kinematic viscosity of 3e6 m2/s. A pressure wave, with a peak of 1333.3 Pa, is applied at the tube inlet for a duration of 3e3 s. The outlet pressure is held at 0 Pa throughout.
The tube wall is assumed an isotropic elastic body under the smallstrain regime, modelled with Hooke's law, for consistency with the original publication that proposed this benchmark. The density is 1200.0 kg/m3, Young's modulus is 300 kPa, and Poisson's ratio is 0.3.
The case demonstrates a strong coupling between the fluid and the solid due to the high fluidtosolid density ratio. When using a partitioned solution approach (as in solids4foam), the socalled ‘addedmass operator' makes the problem difficult to solve due to numerical instabilities.
Currently, in solids4foam, there are two classes of approaches for partitioned FSI coupling:

DirichletNeumann coupling, where a Dirichlet condition is applied to the fluid velocity and Neumann conditions to the fluid pressure and solid displacement. This approach does not require modification of the underlying fluid solver; see Tuković Ž, Karač A, Cardiff P, Jasak H, Ivanković A, 2018, OpenFOAM finite volume solver for fluid–solid interaction. Trans FAMENA, 42(3):1–31.10.21278/TOF.42301;

RobinNeumann coupling, where a Dirichlet condition is applied to the fluid velocity, a Robin condition to the fluid pressure, and a Neumann condition to the solid displacement. This approach does require modification of the underlying fluid solver, and hence it cannot be considered a blackbox coupling approach; see Tuković Ž, Bukač M, Cardiff P, Jasak H, Ivanković A, 2018, Added mass partitioned fluid–structure interaction solver based on a robin boundary condition for pressure. In: OpenFOAM selected papers of the 11th workshop. Springer, Berlin, pp 1–23.
For each of these two classes of approach, we can employ different methods to accelerate the FSI iteration loop convergence; within solids4foam, we can use:

Aitken's dynamic relaxation;

the IQNILS algorithm of Joris Degroote, Robby Haelterman, Sebastiaan Annerel, Peter Bruggeman, Jan Vierendeels, Performance of partitioned procedures in fluid–structure interaction, Computers & Structures, 88, 7–8, 2010, 10.1016/j.compstruc.2009.12.006.
As well as using these acceleration algorithms, we can also use a weakly compressible fluid model rather than the standard fully incompressible model; for FSI, weakly compressible fluid models how been shown to improve convergence, for example, see E. Tandis and A. Ashrafizadeh, "A numerical study on the fluid compressibility effects in strongly coupled fluid–solid interaction problems," Engineering with Computers, 2019, doi: 10.1007/s00366019008804..
In this tutorial, we will compare six variants of the approaches above:
 DirichletNeumann formulation with Aitken's acceleration and an incompressible fluid model: We use the
pimpleFluid
fluid model.  DirichletNeumann formulation with IQNILS acceleration and an incompressible fluid model: We use the
pimpleFluid
fluid model.  DirichletNeumann formulation with Aitken's acceleration and a weakly compressible fluid model: We use the
sonicLiquidFluid
fluid model.  DirichletNeumann formulation with IQNILS acceleration and a weakly compressible fluid model: We use the
sonicLiquidFluid
fluid model.  RobinNeumann formulation with an incompressible fluid: Support for this approach is implemented in the
pimpleFluid
model, where we apply special interface boundary conditions:elasticWallPressure
for the fluid pressure field; andelasticWallVelocity
for the fluid velocity field.  DirichletNeumann formulation with Aitken's acceleration and an incompressible fluid model using preCICE: This approach is the same as approach 1, except the preCICE coupling implementation is used. preCICE is an opensource coupling library for partitioned multiphysics simulations. In this case,
solids4Foam
is used as the solid solver, and standardpimpleFoam
is used as the fluid solver (as opposed to solids4foam'spimpleFluid
fluid model). Due to implementation differences, this preCICE IQNILS approach may perform differently than the builtin solids4foam approach.
In all approaches, the solid domain setup is exactly the same, where an incremental small strain formulation is used (the linearGeometry
solid model). Onequarter of the tube's crosssection is considered, although the case could actually be modelled as 2D axisymmetric. The test is run for 0.02 s. A relatively tight FSI loop tolerance of 1e6 is used for all approaches based on the interface motion. For approach 6 (preCICE), the relative displacement tolerance was set to 1e6, and the relative force tolerance was set to 1e3.
In all cases, the firstorder Euler time scheme is used for the solid and fluid.
Running the Case
For approaches 1 to 4, the tutorial case is located at solids4foam/tutorials/fluidSolidInteraction/3dTube
. The case can be run using the included Allrun
scripts: Allrun.pimpleFluid
for approaches 1 and 2, i.e../Allrun.pimpleFluid
, and ./Allrun.sonicLiquidFluid
for approaches 3 and 4. These scripts update the case with links to the correct files to be used by each approach. The acceleration algorithm, Aitken's or IQNILS, is specified in the constant/fsiProperties
file.
For approach 5 (RobinNeumann formulation), you should change the boundary conditions in the 0/fluid/p
and 0/fluid/U
files to the elastic*
versions; currently, the Robin approach is only available in foamextend4.1.
The tutorial case for approach 6, which uses preCICE, is located at solids4foam/tutorials/fluidSolidInteractionpreCICE/3dTube
; once again, the Allrun
script runs the case.
Remember that a tutorial case can be cleaned and reset using the included Allrun
script, i.e. ./Allclean
.
Expected Results
This case has been proposed as a benchmark for FSI problems. The solution for point A's (see Figure 1) axial and radial displacement is available in: A. Lozovskiy, M. A. Olshanskii, and Y. V. Vassilevski, "Analysis and assessment of a monolithic FSI finite element method," Computers and Fluids, vol. 179, pp. 277–288, 2019, doi: 10.1016/j.compfluid.2018.11.004.
Large Time Step
Initially, we compare the solutions using a relatively large time step size of 1e4 s. Figure 2 shows the radial displacement at point A vs time, Figure 3 shows the axial displacement at point A vs time, and Figure 4 shows the iterations per time step. In the figures, the approaches are designated as:
 DNFAitken: DirichletNeumann formulation with Aitken's acceleration and an incompressible fluid model;
 DNFIQNILS: DirichletNeumann formulation with IQNILS acceleration and an incompressible fluid model;
 WCDNFAitken: DirichletNeumann formulation with Aitken's acceleration and a weakly compressible fluid model;
 WCDNFIQNILS: DirichletNeumann formulation with IQNILS acceleration and a weakly compressible fluid model;
 RNF: RobinNeumann formulation with an incompressible fluid;
 preCICEDNIQNILS: DirichletNeumann formulation with Aitken's acceleration and an incompressible fluid model using preCICE.
Figure 2: Axial displacement at point A vs time with deltaT = 1e4 s
Figure 3: Radial displacement at point A vs time with deltaT = 1e4 s
Figure 4: Number of FSI iterations per timestep with deltaT = 1e4 s
The predictions from all approaches agree closely. Examining the number of FSI iterations per time step, both implementations (solids4foam and preCICE) of DirichletNeumann coupling with IQNILS acceleration are seen to require the least number of iterations. The weakly compressible approach is the next best performing approach, while the incompressible Aitken'saccelerated DirichletNeumann and RobinNeumann approaches are seen to perform the poorest.
Small Time Step
To observe the effect of the time step size, the cases were rerun with a smaller time step of 2.5e5 s, where Figures 5, 6 and 7 show the radial displacement, axial displacement and the number of iterations.
Figure 5: Axial displacement at point A vs time with deltaT = 2.5e5 s
Figure 6: Radial displacement at point A vs time with deltaT = 2.5e5 s
Figure 7: Number of FSI iterations per timestep with deltaT = 2.5e5 s
Unlike the larger time step, the RobinNeumann approach now requires the least number of iterations per time step (exactly 4 for every time step). The weakly compressible approaches are the next best, where the IQNILSaccelerated compressible model outperforms the Aitken'saccelerated compressible model. The incompressible IQNILS approaches are the next best (both solids4foam and preCICE). Finally, in this case, the Aitken'saccelerated incompressible model shows the poorest performance, requiring an order of magnitude greater number of iterations than the best approach. The impressive performance of the RobinNeumann approach can also be observed for smaller time steps; in general, for cases like this, if the time step is sufficiently small, the RobinNeumann approach requires minimal iterations; however, once the time step is large, then the Robin approach diverges or becomes uncompetitive.
FirstOrder Euler vs SecondOrder Backward Time Schemes
Finally, we demonstrate the effect of the time scheme by comparing the firstorder Euler
and secondorder backward
time schemes, where both the solid and fluid use the same time schemes. For this, we rerun the case with five time step sizes: 8e4, 4e4, 2e4, 1e4 and 5e5 s. Figure 8 shows the firstorder Euler radial displacements results, Figure 9 shows the secondorder backward results, and Figure 10 compares both. The secondorder scheme shows small time discretisation errors for each time step and is seen to converge more quickly as the time step size is reduced, as expected. The results were generated using the weakly compressible fluid model and Aitken's coupling; however, the same behaviour is expected with the other modelling and coupling approaches.
Figure 8: The effect of the time step size when using the firstorder Euler
time scheme
Figure 9: The effect of the time step size when using the secondorder backward
time scheme
Figure 10: Comparing the firstorder Euler
and secondorder backward
schemes as the time step size is reduced
Table 1 gives the total clock time for each model when running these models in serial. Surprisingly, the fastest model is the backward method with the smallest time step (5e5 s). The fastest Euler model uses the largest time step (8e4 s), but next fastest model uses the second smallest time step (1e4 s). The explanation for this is that when using larger time steps, a much greater number of outer fluidsolid interaction iterations are required; whereas for smaller time steps, much fewer iterations are required, possibly resulting in a lower overall clock time. If the time step sizes were reduced further, it would be expected that at some point, the models would start to become slower (e.g., as shown by the Euler method when going from 1e4 to 5e5 s); however, these results highlight the importance of selecting a sufficiently small time step size for fluidsolid interaction simulations.
Time Step Size (in s)  Euler (in s)  backward (in s) 

8e4  1568  2233 
4e4  2618  2710 
2e4  3009  2451 
1e4  2021  1999 
5e5  2718  1648 
Table 1: Clock times for the Euler and backward with different time step sizes
Data Availability
The results and gnuplot scripts used to generate the figures above are available in the solids4foam tutorials benchmark data repository.