🔢 mathematics

A Jacobian-free Newton-Krylov method for high-order cell-centred finite volume solid mechanics

This paper presents a robust Jacobian-free Newton-Krylov method coupled with third- and fourth-order cell-centred finite volume formulations for solid mechanics, which significantly enhances accuracy in resolving complex stress fields while maintaining computational efficiency and being implemented in the open-source solids4foam toolbox.

Original authors: Ivan Batistic, Pablo Castrillo, Philip Cardiff

Published 2026-07-21

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Ivan Batistic, Pablo Castrillo, Philip Cardiff

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict how a giant, complex structure—like a bridge, a car chassis, or even a human bone—will bend, stretch, or snap under pressure. Scientists use powerful computers to simulate these scenarios, but the math is incredibly tricky. They have to break the object down into millions of tiny puzzle pieces (called "cells") and calculate the forces acting on each one. For decades, the standard way of doing this was like using a low-resolution camera: it gave a decent picture, but if you zoomed in too close, the image got blurry and jagged. To get a clearer picture, you'd have to use millions more puzzle pieces, which slowed the computer down to a crawl.

Recently, scientists have been trying to build "high-definition" simulators that can see the fine details without needing a billion puzzle pieces. However, these high-definition methods are notoriously difficult to solve because the math equations become a tangled, massive knot that is hard to untangle. The usual tools for untangling them are either too slow or require so much computer memory that they crash the system. This paper tackles that exact problem: how do we get a crystal-clear, high-definition simulation of solid objects without breaking the computer's brain?

The authors, a team of researchers from Ireland, Croatia, and Uruguay, have developed a new way to solve these high-definition simulations. Think of their method as a clever shortcut. Instead of trying to map out every single connection in the massive knot of math (which is like trying to draw every thread in a spiderweb), they use a "Jacobian-free" approach. This is like navigating a maze by feeling the walls and taking a few steps at a time, rather than trying to memorize the entire map before you start walking. They combine this with a "high-order" technique, which is like using a super-smart interpolation tool to guess the shape of the curve between puzzle pieces with incredible precision, rather than just drawing straight lines.

The team tested this new "smart navigator" on a variety of virtual objects, from simple beams and plates with holes to thick-walled cylinders and even materials that stretch like rubber. They found that their method could achieve much higher accuracy than the old, standard methods. In fact, for some tests, their high-definition approach was so good that it could reproduce the exact mathematical answer almost perfectly, whereas the old methods would always have a little bit of "fuzziness" or error.

Crucially, they showed that this high-precision method didn't have to be slow or memory-hungry. By using a simplified "preconditioner"—a rough, low-resolution sketch of the problem to help guide the solver—they managed to keep the calculations fast and efficient. They also added a special "stabilizer" (called $\alpha$ -stabilisation) to stop the simulation from getting jittery or shaking uncontrollably on messy, irregular shapes, which is a common problem in these types of simulations.

The results, which were run on a suite of 2D and 3D test cases involving both stiff metals and stretchy rubber-like materials, suggest that this approach is a major step forward. It proves that you can have your cake and eat it too: you can get the high accuracy of a high-definition simulation without the massive computational cost that usually comes with it. The researchers have even made their code available to the public, hoping that other scientists will use it to build even better simulations for engineering and science. While the paper doesn't claim to have solved every problem in the world of physics, it demonstrates that combining these specific high-order techniques with this clever "Jacobian-free" solver is a powerful and robust way to simulate how solid things behave in the real world.

Technical Summary: A Jacobian-free Newton-Krylov Method for High-Order Cell-Centred Finite Volume Solid Mechanics

Problem Statement
While high-order finite-volume discretisations (achieving third-order accuracy or higher) have seen significant development in computational fluid dynamics (CFD) over the past two decades, their application to solid mechanics remains underexplored. Existing high-order solid mechanics formulations often rely on segregated solvers with linear convergence or fully implicit Newton–Raphson approaches that require the explicit assembly and storage of complex Jacobian matrices. The latter is computationally expensive and memory-intensive, particularly for nonlinear material models where deriving exact Jacobians is non-trivial. Furthermore, previous high-order solid mechanics implementations often required custom code prototypes or specific boundary treatments that hindered integration with standard open-source frameworks. There is a need for a robust, efficient, and flexible method that combines high-order accuracy with the computational advantages of Jacobian-free techniques on general unstructured meshes.

Methodology
The authors extend the Jacobian-free Newton–Krylov (JFNK) framework to third- and fourth-order cell-centred finite volume formulations for both linear and nonlinear solid mechanics. The methodology is implemented within the solids4foam toolbox for OpenFOAM, ensuring compatibility with arbitrary polyhedral meshes.

High-Order Discretisation: The method employs a deconvolution finite volume approach where solution nodes are located at cell centroids. High-order accuracy is achieved through:
- Gradient Reconstruction: A local least-squares reconstruction (Local Regression Estimator, LRE) is used to evaluate displacement gradients and higher-order derivatives at cell faces and quadrature points.
- Flux Integration: Gaussian quadrature is applied at cell faces, with polygonal faces subdivided into triangular subfaces in 3D to facilitate accurate integration.
- Stabilisation: An $\alpha$ -stabilisation scheme is incorporated to dampen high-frequency error modes (zero-energy modes) on irregular meshes. Unlike second-order Rhie–Chow stabilisation, this scheme preserves the formal high-order accuracy of the discretisation.
JFNK Solution Strategy: The nonlinear system of equations is solved using a JFNK algorithm.
- Jacobian-Free: The Jacobian matrix is not explicitly formed. Instead, matrix-vector products are approximated via finite differencing of the residual vector ( $Jv \approx [R(u+\epsilon v) - R(u)]/\epsilon$ ).
- Preconditioning: To ensure convergence and robustness, a compact-stencil approximate Jacobian is used as a preconditioner. This preconditioner is derived from a lower-order (second-order) spatial discretisation, providing a "physics-based" approximation that maintains the efficiency of the JFNK approach without the memory overhead of the full high-order Jacobian.
- Implementation: The solver utilizes the PETSc library for the Krylov subspace method (GMRES) and is integrated into the OpenFOAM ecosystem, handling general unstructured meshes without the need for additional boundary unknowns required in previous custom implementations.

Key Contributions

First Demonstration of JFNK for High-Order Solids: This work presents the first successful application of a JFNK solver to high-order (beyond second-order) finite-volume solid mechanics formulations.
Integration with Open-Source Frameworks: Unlike previous custom-code prototypes, this implementation is built within the solids4foam toolbox for OpenFOAM. It supports arbitrary polyhedral cells and streamlines boundary condition treatment, removing the need for additional boundary unknowns.
Robust Preconditioning Strategy: The study demonstrates that a compact-stencil Jacobian derived from a lower-order scheme remains an effective preconditioner for high-order Krylov solvers, preserving the efficiency and fast convergence of the JFNK approach.
$\alpha$ -Stabilisation for High-Order Schemes: The paper introduces and validates an $\alpha$ -stabilisation scheme specifically adapted for high-order solid mechanics, showing its effectiveness in suppressing oscillations on unstructured grids while scaling consistently with the order of accuracy.

Results
The proposed methodology was benchmarked against a suite of two- and three-dimensional test problems involving linear elastic, hyperelastic (Mooney–Rivlin, neo-Hookean), and large-strain regimes. Key findings include:

Accuracy: The higher-order formulations (3rd and 4th order) deliver substantial accuracy improvements over standard second-order schemes. In the "patch test" (Cantilever Beam), the fourth-order scheme reproduced the exact analytical solution (machine precision) across all mesh types.
Convergence Rates: The solver achieved the expected theoretical orders of convergence for stress and displacement fields on various mesh topologies (hexahedral, tetrahedral, polyhedral, regular, and irregular). Notably, the presence of curved boundaries (e.g., in the hole-in-plate and pressurised cylinder cases) did not degrade convergence rates, even without specialized isoparametric boundary treatments.
Performance: The JFNK approach achieved strong performance with minimal modifications to existing segregated frameworks. The use of a lower-order preconditioner proved sufficient to maintain robustness and efficiency, avoiding the computational cost of assembling full high-order Jacobians.
Robustness: The method successfully handled near-incompressible materials (Poisson's ratio $\approx 0.5$ ) and complex nonlinear material laws without the pressure oscillations typically seen in lower-order displacement-based formulations.

Significance
The paper concludes that combining higher-order finite-volume methods with JFNK solvers represents a significant advancement in computational solid mechanics. This synergy offers a pathway to solve complex stress and deformation fields with high accuracy while retaining the flexibility of finite-volume methods on unstructured meshes. By eliminating the need for explicit Jacobian assembly and demonstrating robustness on irregular grids, the proposed method addresses key bottlenecks in large-scale simulations. The authors emphasize that these findings underscore the potential of this approach to advance the state of the art, particularly for multiphysics simulations where a unified finite-volume framework is advantageous. The implementation is openly released to facilitate further exploration and adoption by the community.

More like this