Computation of viscoelastic shear shock waves using finite volume schemes with artificial compressibility

TL;DR

This work develops a 3D shock‑capturing finite‑volume framework for incompressible viscoelastic solids using Fung–Simo QLV theory with Yeoh elasticity to model soft tissues. An artificial compressibility (AC) approach replaces strict incompressibility, enabling a robust split scheme that couples hyperbolic flux updates with memory‑variable relaxation, and MUSCL reconstruction enhances accuracy in nonlinear regimes. The method is validated against linear viscoelastic benchmarks and extended to nonlinear elastic and viscoelastic scenarios, where shear shocks, dispersion, and compression–shear coupling are observed in 1D, 2D, and 3D settings. Results demonstrate the method’s ability to capture shock formation and focusing phenomena relevant to brain mechanics, while highlighting numerical costs and parameter choices that balance accuracy and efficiency. The code base is made available, offering a platform for future enhancements such as adaptive meshing and more advanced predictor–corrector strategies.

Abstract

The formation of shear shock waves in the brain has been proposed as one of the plausible explanations for deep intracranial injuries. In fact, such singular solutions emerge naturally in soft viscoelastic tissues under dynamic loading conditions. To improve our understanding of the mechanical processes at hand, the development of dedicated computational models is needed. The present study concerns three-dimensional numerical models of incompressible viscoelastic solids whose motion is analysed by means of shock-capturing finite volume methods. More specifically, we focus on the use of the artificial compressibility method, a technique that has been frequently employed in computational fluid dynamics. The material behaviour is deduced from the Fung--Simo quasi-linear viscoelasiticity theory (QLV) where the elastic response is of Yeoh type. We analyse the accuracy of the method and demonstrate its applicability for the study of nonlinear wave propagation in soft solids. The numerical results cover accuracy tests, shock formation and wave focusing.

Figures (9)

• Figure 1: Dissipation factor $-\mathfrak{Im}(\kappa^2)/\mathfrak{Re}(\kappa^2)$ of infinitesimal harmonic waves deduced from Eq. \ref{['Dispersion2D']}, compared to the case of power-law attenuation tripathi19b. The vertical dashed lines mark the relaxation frequencies $\omega_\ell/(2\pi)$ of Table \ref{['tab:Parameters']}.
• Figure 2: Right-going shear waves in linear incompressible elasticity with sinusoidal forcing. Left: shearing velocity obtained numerically using the artificial compressibility method \ref{['FVDisp1D']} for $\epsilon = 0.9$ and various mesh sizes $\Delta x$. Right: numerical error in $L^2$-norm with mesh-dependent AC parameter $\epsilon = 4\, (\Delta x/L)^{0.3}$ (solid lines) and with constant AC parameter $\epsilon = 0.9$ (dashed lines).
• Figure 3: Evolution of the relative error in kinematic energy norm for several values of the compressibility ratio $\epsilon$.
• Figure 4: Right-going shear waves in linear incompressible elasticity, Cauchy problem \ref{['CauchyProb']}. Left: shearing velocities obtained numerically using the artificial compressibility method \ref{['FVDisp1D']} for $\epsilon = 0.9$ and various mesh sizes $\Delta x$. Right: numerical error in $L^2$-norm with mesh-dependent AC parameter $\epsilon = 4\, (\Delta x/L)^{0.3}$ (solid lines) and with constant AC parameter $\epsilon = 0.9$ (dashed lines).
• Figure 5: Right-going shear wave in linear incompressible viscoelasticity. Shearing velocity obtained numerically using the artificial compressibility method \ref{['FVDisp1D']} with MUSCL reconstruction and $\epsilon = 0.9$ at various mesh sizes $\Delta x$.

Theorems & Definitions (1)

• Remark