Fluid-solid coupling schemes: fluidSolidInterfaces

Fluid-solid coupling schemes are implemented in solids4foam via the fluidSolidInterface base class and derived runtime-selectable fluid-solid interface coupling schemes.

The base class is implemented in fluidSolidInterface/fluidSolidInterface.{H,C}. It manages the fluid and solid models, interface mapping, residual evaluation, outer-correction controls, and shared coupling data.

Available Schemes

The following schemes are available through the fluidSolidInterface entry in constant/fsiProperties.

weakCoupling

Weak Dirichlet-Neumann coupling without outer FSI correctors in each time step. This is the cheapest option, but it is also the least robust for strongly coupled problems.

Implementation: weakCouplingInterface/weakCouplingInterface.{H,C}

fixedRelaxation

Strong Dirichlet-Neumann coupling with a fixed under-relaxation factor. This is a simple strong-coupling scheme and is often used as a baseline or for mildly coupled problems.

Implementation: fixedRelaxationCouplingInterface/fixedRelaxationCouplingInterface.{H,C}

Aitken

Strong Dirichlet-Neumann coupling with Aitken dynamic under-relaxation. This improves on fixed relaxation by adapting the relaxation factor during the outer iterations.

Implementation: AitkenCouplingInterface/AitkenCouplingInterface.{H,C}

IQNILS

Strong Dirichlet-Neumann coupling accelerated with the interface quasi-Newton inverse least-squares (IQN-ILS) method. This scheme reuses secant information from previous coupling iterations and time steps to improve convergence for challenging FSI problems. The implementation includes filtering of repeated or near-linearly-dependent secant modes before the least-squares solve to improve numerical robustness.

Common IQNILS controls in constant/fsiProperties include:

  • relaxationFactor: startup relaxation used before enough secant information is available.
  • couplingReuse: number of previous time steps whose secant information is retained.
  • minSignificant: absolute filtering tolerance for repeated or near-dependent modes.
  • relMinSignificant: relative QR2-style filtering tolerance based on the fraction of new information in each older secant mode. Defaults to 1e-3 when couplingReuse is greater than zero; otherwise it is disabled.
  • normalizeCouplingColumns: optionally normalize each V/W secant column pair before the QR solve to improve the conditioning of the least-squares system.
  • residualSumPreconditioning: optionally scale each V column and the residual by the inverse of the current pointwise residual magnitude in the preciceStyleCouplingQR solve. Unlike preCICE's residual-sum preconditioner, which weights whole data blocks, the weight is per point.
  • residualPreconditioningEpsilon: regularization added to the residual magnitude when computing the above scaling weight.
  • reusePreviousStepModes: optionally reuse the previous time step's cached secants only in the first IQN-ILS solve of the new time step. If the immediately previous time step converged in one iteration and generated no new secants, the implementation falls back to the latest non-empty cached time step instead, which is closer to preCICE's backup LS-system behavior.
  • maxReuseUpdateNormRatio: optional safeguard that caps the norm of the interface displacement update, relative to the current interface residual, in every IQN-ILS solve that uses secants from previous time steps (couplingReuse or reusePreviousStepModes).
  • combinedCouplingSystem: optionally assemble one aligned least-squares system across all coupled interfaces instead of solving one interface at a time. On single-interface cases this should reproduce the per-interface result; it is mainly a structural step toward a more preCICE-like assembly. This can also be combined with preciceStyleCouplingQR.
  • requireAllResidualMeasures: optionally require all normalized residual measures in the shared FSI convergence check to satisfy the tolerance. Legacy behavior uses the minimum of the available residual measures, which can stop earlier when only one measure has decayed.
  • outerCorrAbsoluteTolerance: optional absolute interface residual tolerance in metres, used by the strong coupling schemes (see Absolute residual tolerance).
  • qrSolveTolerance: relative cutoff used to ignore nearly singular directions during the triangular back-substitution.
  • reorthogonalizeCouplingColumns: optionally apply a second Gram-Schmidt-style correction while assembling the QR system.
  • predictSolid: optionally solve the solid once before the outer FSI loop.
  • skipModesOnSolidPredictor (default yes): do not apply the secant modes to the residual of the solid predictor. That residual comes from a solid solve with the previous step's fluid force, not from the fluid-solid coupling the modes describe, and applying reused modes to it can move the interface by orders of magnitude too much (#587). The first iteration of a step then uses the predictor/relaxation update, and the modes are first applied after a complete fluid-solid iteration, as in preCICE and CoCoNuT.

With couplingReuse greater than zero, the secants of a step accepted unconverged (allowUnconvergedCoupling yes) are discarded rather than reused.

Implementation: IQNILSCouplingInterface/IQNILSCouplingInterface.{H,C}

oneWayCoupling

Pseudo-coupling scheme for one-way FSI. The fluid solution is assumed to be known already, and the fluid fields are read from a pre-run fluid case and applied to the solid.

Implementation: oneWayCouplingInterface/oneWayCouplingInterface.{H,C}

thermal

Strong thermal coupling interface with fixed under-relaxation. This extends the mechanical coupling infrastructure to temperature and heat-flux transfer and can optionally include mechanical coupling as well.

Implementation: thermalCouplingInterface/thermalCouplingInterface.{H,C}

Immersed interfaces

A solid interface patch may drive an immersed body of the fluid (the immersedBoundaryForce finite volume option of src/immersedBoundary, OpenFOAM.com only) instead of a fluid patch: the fluid mesh is then fixed and needs no motion solver or interface mapping. The fluid patch of the interface is none, and the immersedInterfaces dictionary, keyed by the solid patch name, gives the body, whose motion must be fsiDriven, and optionally the closurePatches: other patches of the solid (e.g. a clamped root) such that the union of the interface and closure patches is a closed surface with outward normals. In a two-dimensional fluid mesh the patches of the empty direction are not needed: the surface is extended through the mesh and capped in that direction.

solidPatch      plate;
fluidPatch      none;

immersedInterfaces
{
    plate
    {
        body            flag;          // body of the fluid fvOptions
        closurePatches  (plateFix);    // optional
    }
}

The lists solidPatches/fluidPatches may mix immersed (none) and body-fitted interfaces. In each coupling iteration the immersed surface is moved to the solid interface displacement of the coupling scheme (relaxed or accelerated as for a fluid patch), with the velocity of the displacement increment over the time step, and the traction on the solid patch is the surface traction of the immersed boundary averaged over the quadrature points of each face. The log compares its total with the momentum exchange between the fluid and the body. The immersed interfaces work with fixedRelaxation, Aitken and IQNILS, and not with weakCoupling, oneWayCoupling, thermal or the Robin (elasticWallPressure) interface conditions. See the fluids/immersedBoundary/immersedHronTurekFsi2 tutorial.

Implementation: fluidSolidInterface/fsiImmersedInterface.{H,C} builds the surface and maps the tractions, and fluidSolidInterface/fsiImmersedBoundary.H is the header-only interface implemented by immersedBoundaryForce.

Fluid mesh motion

The fluid mesh is moved by the mesh motion solver in constant/fluid/dynamicMeshDict. At each mesh update, the interface gives a fluid mesh motion coupler the change in the displacement of each fluid interface patch, and the coupler sets it in the form its family of motion solvers expects. The coupler is selected from the motion solver, or set explicitly in constant/fsiProperties:

fluidMeshMotion
{
    type    auto;   // default
}
  • velocity: velocity-based solvers, e.g. velocityLaplacian; the pointMotionU patch value is the displacement increment divided by the time-step size.
  • displacement: displacement-based solvers, e.g. displacementLaplacian and displacementSBRStress; the pointDisplacement patch value is the total displacement relative to the solver reference points, points0.
  • solidModel: solids4foamSolidModel (OpenFOAM.com); the total displacement is set on the solid model pointD and fixedDisplacement D.
  • RBF: RBFMeshMotionSolver; the patch face displacement increments are set.
  • tetFemVelocity: foam-extend tetFEM solvers; the motionU patch value is the displacement increment divided by the time-step size.
  • subsetMesh: foam-extend newSubsetMotionSolverFvMesh; pointMotionU is set on the sub-mesh.

The FSI patches must have a fixedValue condition in the motion field. Component solvers, e.g. velocityComponentLaplacian, are not supported.

Velocity- or displacement-based motion

Velocity-based solvers are path dependent: the interior mesh is advanced from its current state each time step, so under periodic or long-running motion the interior mesh can drift cycle after cycle and periodic quantities of interest may not become stationary. A displacement-based solver, e.g. displacementLaplacian with the same diffusivity, always solves from the reference points, so the mesh for a given interface position does not depend on how the interface got there. It is recommended for periodic or long-running FSI cases; add a 0/fluid/pointDisplacement file, with fixedValue on the FSI patches, in place of pointMotionU.

The displacement coupler sets the interface value from the current geometry, points + increment - points0, rather than accumulating increments, so it is the same whatever the number of outer iterations, it includes the interface points moved directly when interfaceDeformationLimit is used, and it does not depend on the earlier boundary values. A restart still needs a pointDisplacement field to read, and points0 must be the reference points of the run being continued.

Solid model mesh motion

solids4foamSolidModel moves the fluid mesh with the meshMotionTotalDisplacement solid model, solved on the fluid mesh itself, typically with the diffusionElastic or diffusionHyperElastic law, whose stiffness is scaled by a mesh motion diffusivity so that the mesh is stiffer near the interface:

// constant/fluid/dynamicMeshDict
dynamicFvMesh   dynamicMotionSolverFvMesh;
motionSolver    solids4foamSolidModel;
diffusivity     quadratic inverseDistance (interface);

As with displacementLaplacian, the solid model solves for the total displacement from points0, discretised on the current fluid mesh, so no copy of the mesh is needed; with finiteStrain yes, for the diffusionHyperElastic law, its kinematics are those of the current configuration (see the solid model page). The stiffness scale factor is recalculated on the current mesh before each solution.

The solid model reads solidProperties and mechanicalProperties from constant/fluid, D and pointD from 0/fluid, and its schemes and solver settings from system/fluid/fvSchemes and fvSolution: grad(D) (leastSquaresS4f), laplacian(DD,D), and the D solver and relaxation factor. On the FSI patches, D must be fixedDisplacement with enforcePointPatchField no, and pointD must be fixedValue. The solid model fields whose names the fluid may also use, e.g. U, p, rho and g, are private to it, and are neither read nor written; it writes its D, pointD and sigma fields, and the state of its law, to the fluid time directories. The beamInCrossFlow tutorial has an example: ./Allrun diffusionElastic.

Notes

  • The partitioned schemes share the common fluidSolidInterface machinery for transferring tractions and displacements across the interface.
  • Different schemes have different stability and cost trade-offs. In general, weakCoupling is the cheapest, fixedRelaxation and Aitken are simple strong-coupling options, and IQNILS is the most sophisticated partitioned acceleration scheme in this directory.
  • fixedRelaxation, Aitken and IQNILS start each time step from the solid predictor (predictSolid, default yes): the solid is solved once with the fluid force of the previous time step. With predictor (default yes), the first iteration moves the fluid interface with the full predicted solid displacement; later iterations are relaxed or accelerated as usual. The interface displacement is an increment within the time step, so with predictor no the first iteration moves the fluid interface by only relaxationFactor times the solid's step, and the interface stops, or with backward time differencing reverses, as seen by the fluid. On cases with a strong added-mass effect and a small relaxationFactor, the resulting first-iteration fluid force can be hundreds of times the converged one and make the solid diverge (#489).
  • For IQNILS, increasing couplingReuse can improve convergence, but it is not guaranteed to make a case monotonically more robust. The best reuse level is case dependent and interacts with time-step size and the quality of the inner fluid and solid solves. Conversely, couplingReuse 0 with a small relaxationFactor is not robust on cases with a strong added-mass effect: after the full predictor step of the first iteration, the next iterations have no secants yet and barely move the interface, and the solid can diverge within the first few iterations of a step. Reused secants recover the step from the second iteration.

Absolute residual tolerance

The shared convergence check normalizes the interface residual in two ways: residual1, relative to the largest residual in the time step, and residual2, relative to the largest total interface displacement. The legacy default accepts a step when min(residual1, residual2) < outerCorrTolerance, which can stop after one iteration through residual2 while residual1 is still large in quasi-steady phases, so the coupling error accumulates (#552). The log prints a note when the legacy test is passed by one measure only. requireAllResidualMeasures yes uses the maximum instead, but near a steady state the residual is tiny from the first iteration, so the relative reduction residual1 can be unattainable and every step reaches nOuterCorr.

outerCorrAbsoluteTolerance (default 0, disabled) adds an absolute test: a step is also converged if the maximum point mismatch between the solid and fluid interface, max |residual| over the interface points in metres, is below the value. With the default incrementalResiduals yes this is the mismatch of the interface displacement increments, not the accumulated separation of the interfaces, so small accepted increments can still add up over many steps; with incrementalResiduals no it is the point separation of the interfaces. The test is in addition to the existing test (legacy or requireAllResidualMeasures), so it does not reopen the min() loophole. The maximum is used rather than an RMS because it bounds the local error and cannot hide a localised mismatch. The absolute residual is printed next to the relative residuals, and a line states when it accepts a step. On Robin interfaces the pressure and flux criteria still apply, and the accepted displacement residual enters the Robin convergence history clamped to outerCorrTolerance.

requireAllResidualMeasures yes;
outerCorrAbsoluteTolerance 1e-9; // m

Choose a small fraction of the smallest interface cell size, or of the displacement accuracy you need, e.g. 1e-3 of the cell size or below. With the option unset the behaviour is unchanged.