Solids4Foam: High-Order Mechanics in OpenFOAM

• Solids4Foam is an open-source framework that integrates high-order finite-volume methods for modeling complex solid mechanics and fluid–solid interactions.
• It features a modular, object-oriented architecture with PETSc support, enabling runtime selection of solvers, constitutive models, and stabilization strategies.
• Practical applications include efficient nonlinear convergence via JFNK, validated performance on benchmarks, and seamless integration with OpenFOAM’s fluid solvers.

Solids4Foam is an open-source toolbox for the OpenFOAM CFD platform that implements cell-centred finite-volume (FV) methods for solid mechanics and fluid–solid interaction (FSI), including advanced algorithms for nonlinear solid mechanics, high-order FV discretisations, and tightly coupled multiphysics workflows. It leverages object-oriented design and PETSc library support to offer runtime-selectable solvers, constitutive models, and stabilisation strategies, accommodating a wide range of solid mechanics phenomena and making it possible to combine FV-based solid methods directly with OpenFOAM’s native fluid solvers (Cardiff et al., 2018, Cardiff et al., 24 Feb 2025, Batistic et al., 26 Jan 2026).

1. Core Architecture and OpenFOAM Integration

Solids4Foam adopts a modular architecture consistent with OpenFOAM conventions. The main directory structure includes:

• src/: C++ library code (solvers, material models, numerical utilities)
• applications/solvers/: case-ready executables (e.g., highOrderSolidMechanicsJFNK)
• applications/utilities/: mesh and result post-processors
• tutorials/: sample case setups (Batistic et al., 26 Jan 2026).

The primary class hierarchy is organized as follows:

• PhysicsModel (abstract)
  • SolidModel
  • LinearGeometryModel (segregated/coupled)
  • NonlinearGeometryModel (Total/Updated Lagrangian)
  • FluidModel
  • FSIModel (partitioned Dirichlet–Neumann) (Cardiff et al., 2018).

The modular design promotes encapsulation, inheritance, and polymorphism, enabling seamless runtime selection of constitutive laws, algebraic solvers, and multiphysics coupling strategies via dictionary files, with no recompilation necessary.

Interaction with PETSc is tightly woven into nonlinear solver workflows, with JFNK and other Krylov methods enabled through dynamic linking and class interface wrappers (Cardiff et al., 24 Feb 2025, Batistic et al., 26 Jan 2026).

2. Governing Equations and Finite-Volume Discretisation

The FV framework integrates the strong forms of linear momentum (with Lagrangian kinematics) over each control volume:

$$\int_{\Omega_0} \rho_0 \frac{\partial^2 \mathbf{u}}{\partial t^2} \, d\Omega_0 = \oint_{\Gamma_0} \mathbf{n}_0 \cdot \boldsymbol{\sigma} \, d\Gamma_0 + \int_{\Omega_0} \rho_0 \mathbf{b} \, d\Omega_0,$$

with cell-centred storage of primary solution variables. Surface divergence, body force, and inertia are assembled with second- or higher-order quadrature, and boundary conditions (displacement, traction, symmetry) are enforced via FV-compatible treatments—e.g., weak/penalty enforcement in LRE stencils (Batistic et al., 26 Jan 2026, Cardiff et al., 2018).

Notably, the toolbox supports third- and fourth-order discretisations via local least-squares reconstruction and rigorous quadrature. The gradient at a quadrature point $\tilde{x}$ is represented by a local Taylor expansion of order $p$:

$$u_\alpha(\mathbf{x}) \approx \sum_{|\mathbf{m}|\leq p} \frac{1}{\mathbf{m}!} D^\mathbf{m} u_\alpha(\tilde{\mathbf{x}}) (\mathbf{x} - \tilde{\mathbf{x}})^\mathbf{m},$$

with weights found via local regression using neighbouring cells (stencil size scales with $p$) (Batistic et al., 26 Jan 2026).

3. Solution Algorithms and Jacobian-Free Newton-Krylov Methods

Solids4Foam implements both traditional segregated quasi-Newton methods and Jacobian-free Newton-Krylov (JFNK) solvers. The JFNK approach is a fully-coupled Newton solver where the Jacobian-vector product is approximated via finite differences:

$$J(\mathbf{u})\,\mathbf{v} \approx \frac{R(\mathbf{u} + \varepsilon\,\mathbf{v}) - R(\mathbf{u})}{\varepsilon}, \quad \varepsilon = \mathrm{fmDiffEpsilon}\,(1+\|\mathbf{u}\|).$$

The nonlinear update at iteration $k$ is

$$\mathbf{u}^{k+1} = \mathbf{u}^k + \delta \mathbf{u},$$

where $\delta \mathbf{u}$ solves

$$J(\mathbf{u}^k)\,\delta \mathbf{u} = -R(\mathbf{u}^k)$$

using a Krylov method (flexible GMRES, PETSc implementation), right-preconditioned with a compact-stencil Laplacian corresponding to the original semi-implicit FV operator. Preconditioner options include ILU($k$), MUMPS LU, and HYPRE BoomerAMG (Cardiff et al., 24 Feb 2025, Batistic et al., 26 Jan 2026).

The JFNK approach is tightly embedded as JFNKSolidSolver (C++ class), inheriting from the solids4foam solidModel infrastructure and registered for run-time selection. Configuration uses standard OpenFOAM dictionaries (e.g., selection of JFNK type and PETSc SNES options in fvSolution) (Cardiff et al., 24 Feb 2025).

4. Constitutive Models and Stabilisation Strategies

A broad range of mechanical behaviour is supported via the MechanicalLaw base class. Implemented models include:

• Small-strain linear elasticity: $$\boldsymbol{\sigma} = 2\mu \boldsymbol{\epsilon} + \lambda \operatorname{tr}(\boldsymbol{\epsilon}) I$$
• Hyperelasticity: neo-Hookean (large deformation), St Venant–Kirchhoff, Mooney–Rivlin
• Linear viscoelasticity (Prony series)
• J2 plasticity with isotropic hardening
• Orthotropic elasticity, thermoelasticity, poroelasticity (Cardiff et al., 2018, Batistic et al., 26 Jan 2026).

Zero-energy modes and checker-boarding are suppressed using an $\alpha$-stabilisation scheme, which augments the FV traction at each face by:

$$\mathbf{t}_{\alpha, f} = \alpha \bar{K} |d_{PN}\cdot\mathbf{n}_f|\,(\mathbf{u}^*_{f,N} - \mathbf{u}^*_{f,P}),$$

where $\bar{K} = 2\mu + \lambda$ and $\alpha \approx 0.1$ is optimal for robustness on irregular meshes. High-order extrapolations are used for face displacements $\mathbf{u}^*$. The stabilisation is essential for higher-order FV accuracy on general polyhedral meshes (Batistic et al., 26 Jan 2026).

5. Fluid–Solid Interaction Coupling

Solids4Foam implements FSI via partitioned Dirichlet–Neumann coupling:

1. The fluid solver (e.g., OpenFOAM PISO/SIMPLE) advances with boundary conditions from the current solid displacement.
2. Fluid traction ($\mathbf{n}\cdot\boldsymbol{\sigma}_f$) is computed and transferred to the solid domain.
3. Solid displacement is updated under the new interface traction via the solid solver.
4. The process iterates (Aitken under-relaxation) until interface residuals converge to tolerance (Cardiff et al., 2018).

Interface boundary conditions are mapped as strong traction continuity on the FSI interface:

$$\mathbf{n}\cdot\boldsymbol{\sigma}_s = \mathbf{n}\cdot\boldsymbol{\sigma}_f$$

The FSI workflow leverages OpenFOAM’s mesh mapping infrastructure, maintaining generality for unstructured or moving meshes.

6. Practical Usage, Configuration, and Performance

User setup for solids4foam cases closely matches core OpenFOAM practices, with fields (0/), properties (constant/materialProperties, transportProperties), and solver/numerics controls (system/fvSchemes, fvSolution, controlDict). Key configuration options:

• fvSolution: selection of solver (JFNK, segregated), nonlinear and linear tolerances, maximum iterations, preconditioner, line-search.
• fvSchemes: selection and order of gradient, divergence, interpolation, and Laplacian schemes.
• Dirichlet and Neumann BCs are natively supported; weak enforcement is used for higher-order FV (Batistic et al., 26 Jan 2026, Cardiff et al., 2018).

Typical workflow includes mesh generation, dictionary editing, and starting the appropriate executable. For higher-order FV, stencil sizes and $\alpha$-stabilisation must be tuned (see $p$, $\alpha$ parameters).

Performance benchmarks indicate:

• Quadratic convergence of JFNK in elastic and hyperelastic cases.
• Significant speedups (up to $\sim$300$\times$) over segregated methods in large-scale, nonlinear, or complex FSI problems.
• Robustness on irregular meshes is governed by the $\alpha$-stabilisation; over-integration alone is insufficient (Cardiff et al., 24 Feb 2025, Batistic et al., 26 Jan 2026).
• Memory overhead for GMRES and preconditioner is typically $1$–$3\times$ that of segregated CG+ILU(0).
• LU is fastest for $\leq$50k DOF; multigrid (GAMG/AMG) preferred at larger scales or in parallel (Cardiff et al., 24 Feb 2025).

7. Validation, Extensibility, and Benchmark Problems

Comprehensive validation has been performed on a range of canonical and challenging benchmark problems:

Problem	Mesh Sizes	Features/Findings
MMS cube	125–4M	JFNK and segregated both achieve 2nd order convergence
Spherical cavity	1k–1.6M	JFNK 2nd order displ., 1st order stress; speedup $2$–$10\times$
Cook’s membrane	9–147k	JFNK convergent for elastic/hyperelastic; diverges for current elastoplastic preconditioning
Cantilever beam	—	4th-order FV reproduces cubic solution to machine precision
Thick-walled cylinder	—	High-order method achieves mesh-independent accuracy; p=1 exhibits locking
FSI microbeam	—	Accurate transient tip displacement; robust FSI coupling (Cardiff et al., 24 Feb 2025, Batistic et al., 26 Jan 2026, Cardiff et al., 2018)

Extensibility follows from the OOP design: new constitutive models inherit from MechanicalLaw, new solvers from PhysicsModel. Factory registration and dictionary-driven configuration ensure rapid prototyping and development.

Solids4Foam, through its high-order FV discretisation, robust JFNK algorithms, and open-source availability, serves as a reference platform for research in computational solid mechanics, FSI, and methodological advances combining accuracy, robustness, and scalability. Source code, tutorials, and extensive documentation are openly available at https://github.com/solids4foam (Batistic et al., 26 Jan 2026, Cardiff et al., 24 Feb 2025, Cardiff et al., 2018).