Source-linked AI summary
High order ADER schemes for a unified first order hyperbolic formulation of continuum mechanics: viscous heat-conducting fluids and elastic solids
Michael Dumbser, Ilya Peshkov, Evgeniy Romenski, Olindo Zanotti
TL;DR
The paper addresses numerical solution of the unified HPR first-order hyperbolic model, whose equations include stiff sources and non-conservative products. It applies high-order ADER finite volume and discontinuous Galerkin schemes and reports agreement with reference solutions across viscous-fluid tests, while extending the study to elastic solids. Future work includes non-Newtonian fluids and complex visco-plastic solids.
Problem
The full first-order HPR model with heat conduction had not previously been solved numerically in multiple dimensions with all terms, despite its challenging nonlinear hyperbolic structure.
Method
The paper applies one-step high-order ADER-WENO finite volume and ADER discontinuous Galerkin schemes to HPR systems with stiff source terms and non-conservative products.
Results
The schemes produce numerical convergence results and detailed agreement with analytical or numerical reference solutions for viscous-fluid and elastic-solid test problems.
Takeaways & Limitations
The study demonstrates that the ADER schemes can discretize the full HPR model across test problems ranging from viscous compressible fluids to elastic solids.
Takeaways & Limitations
Future applications are needed for non-Newtonian fluids, complex visco-plastic solids, moving meshes, and adaptive-mesh crack problems.
Abstract
from arXiv · showhide
This paper is concerned with the numerical solution of the unified first order hyperbolic formulation of continuum mechanics recently proposed by Peshkov & Romenski, denoted as HPR model. In that framework, the viscous stresses are computed from the so-called distortion tensor A, which is one of the primary state variables. A very important key feature of the model is its ability to describe at the same time the behavior of inviscid and viscous compressible Newtonian and non-Newtonian fluids with heat conduction, as well as the behavior of elastic and visco-plastic solids. This is achieved via a stiff source term that accounts for strain relaxation in the evolution equations of A. Also heat conduction is included via a first order hyperbolic evolution equation of the thermal impulse, from which the heat flux is computed. The governing PDE system is hyperbolic and fully consistent with the principles of thermodynamics. It is also fundamentally different from first order Maxwell-Cattaneo-type relaxation models based on extended irreversible thermodynamics. The connection between the HPR model and the classical hyperbolic-parabolic Navier-Stokes-Fourier theory is established via a formal asymptotic analysis in the stiff relaxation limit. From a numerical point of view, the governing partial differential equations are very challenging, since they form a large nonlinear hyperbolic PDE system that includes stiff source terms and non-conservative products. We apply the successful family of one-step ADER-WENO finite volume and ADER discontinuous Galerkin finite element schemes in the stiff relaxation limit, and compare the numerical results with exact or numerical reference solutions obtained for the Euler and Navier-Stokes equations. To show the universality of the model, the paper is rounded-off with an application to wave propagation in elastic solids.
1. Introduction
The HPR model provides a unified first-order hyperbolic framework for fluids and solids, representing deformation, rearrangement, and transport through additional state variables and finite-speed waves. The paper investigates high-order ADER discretizations for its challenging nonlinear system with stiff sources and non-conservative products.
- Unified model: The HPR model aims to describe inviscid and viscous fluids, heat conduction, and elastic or visco-plastic solids within one first-order hyperbolic formulation.Its scope extends across fluid and solid mechanics when the continuum description applies.
- Material description: The distortion tensor A represents deformation and rotation of finite-sized material elements, while the relaxation time characterizes their rearrangement ability.The formulation introduces internal material-element structure and a continuum analogue of Frenkel’s rearrangement time.
- Wave formulation: Hyperbolicity follows from convexity of the total energy potential and yields longitudinal and shear waves for momentum transfer.The model treats transverse viscous momentum transfer as dissipative shear-wave propagation.
- Relaxation: Dissipative rearrangements are represented by stiff algebraic source terms, keeping characteristic speeds finite as the relaxation time decreases.This contrasts with hyperbolic Maxwell-Cattaneo-type models whose characteristic speeds can diverge in the zero-relaxation limit.
- Numerical objective: The paper applies ADER-WENO finite volume and ADER discontinuous Galerkin methods to the full HPR system, which includes stiff sources and non-conservative products.The study targets a system that had not previously been solved numerically in multiple dimensions with all terms included.
2.1. Formulation of the model
The HPR formulation is closed by a total energy potential depending on density, entropy, velocity, distortion, and thermal impulse. Its derivatives generate constitutive fluxes and relaxation sources for viscous stress, heat flux, and material-element rearrangement.
- Governing variables: The HPR equations include mass, momentum, distortion, thermal-impulse, entropy, and total-energy evolution equations.The thermal-impulse equation is formally analogous to the momentum equation, with temperature playing the role of pressure.
- Constitutive closure: The energy potential generates all constitutive fluxes and dissipative source terms through its state-variable derivatives.Specifying the energy is therefore a central closure step for the model.
- Energy structure: The total energy is decomposed into molecular, mesoscopic, and macroscopic contributions: E(ρ, s, v, A, J) = E1(ρ, s) + E2(A, J) + E3(v).This decomposition assigns energy to the molecular, material-element, and flow scales.
- Dissipation and heat conduction: The distortion relaxation source models shear-strain dissipation, while the thermal-impulse source models heat exchange between material elements.The associated relaxation times are τ1 for strain rearrangement and τ2 for heat conduction.
- Stiff relaxation limit: The model’s formal stiff-limit connection to Navier–Stokes–Fourier theory is obtained as τ1 →0 and τ2 →0.The relaxation parameters are selected to recover the classical viscous and heat-conducting behavior in this limit.
2.2. Discussion
The HPR model is constructed to be thermodynamically compatible and symmetric hyperbolic when its energy potential is convex. Its characteristic structure contains finite-speed longitudinal and shear waves, with frequency-dependent phase behavior analyzed through dispersion relations.
- Thermodynamic compatibility: The HPR system is overdetermined, with 18 PDEs for 17 unknowns, so consistency requires compatibility with total-energy conservation.The paper addresses this through the thermodynamic structure of the equations.
- Thermodynamic compatibility: The total-energy equation follows automatically when the remaining HPR equations and the required constitutive constraints are satisfied.This closes the overdetermined system and establishes its consistency.
- Well-posedness: Convexity of the total energy with respect to the state variables makes the HPR system symmetric hyperbolic and gives a locally well-posed initial-value problem.The equivalent convexity conditions can be expressed in conservative or primitive variables.
- Characteristic waves: The characteristic structure includes transverse shear waves and longitudinal pressure waves, reflecting the model’s wave-based treatment of transport.Shear waves are finite-speed transverse modes associated with momentum transfer across the mean flow.
- Characteristic speeds: HPR perturbations propagate at finite speeds for every frequency, with low-frequency sound waves approaching the fluid sound speed c0.The paper uses dispersion analysis to obtain phase velocities and attenuation factors as functions of angular frequency.
2.3. Formal asymptotic analysis, Newton’s viscous law and Fourier’s law of heat conduction
A Chapman–Enskog analysis connects the HPR model to Euler and Navier–Stokes–Fourier behavior as relaxation times become small. The stiff limit yields vanishing viscous stresses at zeroth order, classical viscosity at first order, and Fourier heat conduction.
- The formal analysis links the HPR model to classical Navier–Stokes–Fourier theory in the stiff relaxation limit τ1 ≪1 and τ2 ≪1.
- 2.3.1. Asymptotic limit of the viscous stress tensor: In the stiff limit, the distortion tensor A tends toward an orthogonal matrix.
- 2.3.1. Asymptotic limit of the viscous stress tensor: At zeroth order, viscous stresses vanish and the HPR model recovers the inviscid compressible Euler equations.The material elements change volume but not shape in this limit.
- 2.3.1. Asymptotic limit of the viscous stress tensor: A first-order expansion of the distortion tensor produces the classical compressible Navier–Stokes stress tensor under Stokes’ hypothesis.The derivation expands G around an isotropic leading term and uses its deviatoric evolution.
- 2.3.1. Asymptotic limit of the viscous stress tensor: The HPR stress is derived from the quadratic energy contribution E2(A, J), whereas classical Navier–Stokes theory postulates the stress as a constitutive relation.
- 2.3.2. Asymptotic limit of the heat flux: A parallel expansion of the thermal impulse J yields the familiar Fourier heat flux and its heat-conduction coefficient.
- 2.3.3. On the experimental measurement of the model parameters: The conventional viscosity and heat-conductivity coefficients do not uniquely determine the HPR parameters, so experiments are needed to measure relaxation and wave-speed quantities.High-frequency sound experiments can estimate parameters such as c_s and dissipation time, while heat-wave experiments are needed for c_h.
3. ADER finite volume and ADER discontinuous Galerkin finite element schemes
The paper solves the nonlinear hyperbolic HPR system with non-conservative products and stiff sources using one-step ADER finite-volume and discontinuous-Galerkin methods. A local space-time predictor and path-conservative corrector handle high-order evolution and interface jumps.
- The HPR equations are a nonlinear hyperbolic system containing conservative fluxes, non-conservative products, and stiff source terms.
- One-step ADER-FV and ADER-DG schemes provide high-order accuracy in space and time without Runge–Kutta sub-stages.
- The solution is represented by piecewise polynomials using modal bases on simplex elements and tensor-product nodal bases on quadrilateral elements.
- The PNPM framework includes discontinuous Galerkin and finite-volume schemes as limiting cases, with N = M giving DG and N = 0 giving WENO finite volume.
- A local space-time Galerkin predictor evolves each element independently before the corrector couples neighboring elements.The predictor replaces the earlier Cauchy–Kowalevski procedure.
- The weak formulation uses space-time basis functions, temporal upwinding, and the governing flux, non-conservative product, and source terms.
- The path-conservative corrector accounts for jumps across element boundaries and evaluates non-conservative products through interface paths.A Rusanov jump term and straight-line path integral are used, with the maximum interface signal speed controlling dissipation.
- The implementation uses the Rusanov local Lax–Friedrichs solver, although other Riemann solvers can also be employed.
4. Numerical results
The numerical tests compare HPR solutions with exact or reference solutions for inviscid Euler and viscous Navier–Stokes flows. Across vortex, Stokes, and boundary-layer problems, the results support accurate behavior in relevant relaxation limits.
- Numerical convergence studies in the stiff inviscid limit: The convected isentropic vortex is evaluated against the exact compressible Euler solution in the inviscid limit on successively refined meshes.The setup uses periodic boundaries on Ω = [0, 10] × [0, 10] and reaches t = 1.0.
- The first problem of Stokes: The first problem of Stokes compares HPR results from an ADER-DG P3P3 scheme with the exact incompressible Navier–Stokes solution at t = 1.The comparison is performed for multiple viscosities, including µ = 10^-2, 10^-3, and 10^-4.
- Laminar boundary layer over a flat plate: The flat-plate boundary-layer test uses third-order ADER-WENO at t = 10 and compares a vertical velocity cut at x = 0.5 with the Blasius solution.The figure also shows the boundary-layer thickness δ0.99, velocity contours, and velocity profiles.
- Laminar boundary layer over a flat plate: The flat-plate results confirm HPR behavior in the stiff relaxation limit as τ1 → 0, where it reproduces known Navier–Stokes results.The same test also displays components of the distortion tensor A.
4.4. Hagen-Poiseuille flow in a duct
The Hagen–Poiseuille test evaluates steady viscous flow in a duct against the Navier–Stokes reference solution. The HPR computation reproduces the laminar flow and its velocity profile under the stated low-Mach-number setup.
- Hagen-Poiseuille flow in a duct: The test models steady Newtonian flow in a rectangular duct driven by a constant pressure gradient, using the known Navier–Stokes parabolic velocity profile.The setup imposes Δp = −4.8, with mean velocity ū = 1 and maximum velocity u_max = 1.5.
- Hagen-Poiseuille flow in a duct: The numerical results very well reproduce the laminar duct flow, with only a moderate increase in velocity from x = 0 to x = 10.The authors conclude that HPR successfully solves this classical steady viscous-flow test.
4.6. Double shear layer
The HPR model is tested on double-shear-layer and cylinder-flow problems, with ADER-WENO solutions compared against Navier–Stokes references and distortion-based flow visualization.
- Double shear layer: The double shear layer uses a fourth-order ADER-WENO scheme with ν = 2·10−4 on a 200 × 200 periodic grid through t = 1.8.The HPR results are compared with an incompressible Navier–Stokes reference solution.
- Double shear layer: Vorticity contours compare HPR solutions at t = 0.8, 1.2, and 1.8 with a staggered semi-implicit space-time DG Navier–Stokes solution.
- Double shear layer: The distortion component A12 is plotted with 41 contour colors over [-1,1] at four times to visualize the double-shear-layer evolution.
- Circular cylinder: The cylinder case at Re = 150 produces a von Kármán vortex street, whose sound signal has Strouhal number St = 0.175.The reported value is in reasonable agreement with reference values St = 0.183 and St = 0.182.
- Compressible mixing layer: For the compressible mixing layer, HPR results show reasonable qualitative agreement with two Navier–Stokes references, while a distortion component reveals vortex structures more clearly than vorticity.
4.9. 2D Taylor-Green vortex
Taylor–Green vortex tests compare HPR simulations with exact Navier–Stokes and DNS references in two and three dimensions, while distortion components expose vortex structures.
- 2D Taylor-Green vortex: The two-dimensional Taylor–Green vortex is simulated to t = 10 with a fourth-order ADER-DG P3P3 scheme on a 50 × 50 grid.The HPR calculation uses μ = 10−2 and periodic boundaries.
- 2D Taylor-Green vortex: The two-dimensional HPR solution shows excellent agreement with the exact incompressible Navier–Stokes solution for velocity and pressure.The component A11 also reveals the vortex structures.
- 3D Taylor-Green vortex: The three-dimensional test uses the full HPR model at Re = 100 and Re = 200 with a third-order ADER-WENO finite volume scheme through t = 10.The simulations use 224^3 elements and include heat-conduction parameters.
- 3D Taylor-Green vortex: At Re = 100, kinetic-energy dissipation agrees well with DNS data, whereas at Re = 200 the scheme is too dissipative for t > 6.The authors propose refined grids and higher polynomial degrees to investigate the deviation.
- 3D Taylor-Green vortex: The component A11 visualizes developing small-scale structures, and all components of A are shown at t = 10.
4.11. Heat conduction in a gas
The heat-conduction test evaluates HPR against classical viscous-flow formulations, while the viscous-shock test compares HPR profiles with Becker’s exact Navier–Stokes solution.
- Heat conduction: The heat-conduction problem initializes a density discontinuity with uniform pressure and zero velocity, then evolves it using an ADER-DG P3P3 scheme.The HPR parameters include μ = 10−2, α = 2, and κ = 10−2.
- Heat conduction: The heat-conduction figure compares temperature distributions and heat fluxes, using q1 = −κTx for Navier–Stokes and q1 = α2J1T for HPR.
- Viscous shock profile: The viscous-shock test considers a supersonic shock with Ms > 1 and uses Becker’s exact traveling-wave solution available at Pr = 0.75.
- Viscous shock profile: For Ms = 2 and Res = 100, the HPR viscous shock profile agrees excellently with the exact Navier–Stokes solution, aside from a small spurious wave.The wave may result from start-up error caused by the non-equilibrium initialization of A and J.
4.13. Viscous double Mach reflection problem
A viscous double Mach reflection problem tests the HPR model on a Mach 10 shock striking a 30° ramp at two shock Reynolds numbers and against an inviscid reference.
- Problem setup: The problem models a Mach 10 shock impacting a 30° ramp with reflecting-slip, inflow, outflow, and imposed oblique-shock boundary conditions.
- Numerical setup: A third-order P0P2 ADER-WENO finite volume scheme computes the HPR solution on a 1400 × 400 grid to t = 0.2.The parameters use Pr = 0.75 and initialize A = 3√ρ I and J = 0.
- Results: The density contours compare viscous cases with μ = 10−1 and μ = 10−2, corresponding to Res = 100 and Res = 1000, against an inviscid Euler reference.
4.14. Application to solid mechanics
The HPR model applies to solid mechanics through a unified PDE formulation, and its Lamb’s problem results agree closely with linear elasticity.
- The HPR model describes fluid and solid mechanics within one PDE system.
- Lamb’s problem tests wave propagation in an elastic solid with a free surface, a setting conventional Navier–Stokes equations cannot solve.
- The solid-mechanics setup uses an ADER-DG P4 scheme through final time t = 1.3 on a 200 × 100 element grid.
- The HPR simulation is compared with classical linear elasticity using the same ADER-DG scheme and grid.
- The HPR and linear-elasticity computations show excellent agreement in wave fields and very good agreement in the recorded velocity signal.
5. Conclusion
The paper demonstrates that high-order ADER methods can solve the full first-order HPR model across fluid and solid mechanics, while identifying several directions for extending the approach.
- The authors report the first numerical application of a method to the full first-order HPR model with heat conduction.
- High-order one-step ADER finite volume and discontinuous Galerkin schemes handle non-conservative products and stiff source terms across tests from viscous fluids to elastic solids.
- The HPR model represents fluid and solid mechanics as limiting cases of one first-order hyperbolic mathematical model.
- Future work includes moving unstructured meshes, non-Newtonian fluids, complex visco-plastic solids, adaptive meshes, crack generation, and crack propagation.
- The paper also identifies extension of the finite-wave-speed HPR formulation to relativistic regimes as future research.
A. Eigenvalues of matrices Ak
This appendix gives the eigenvalue structure of the viscous subsystem matrices and states the convexity assumption used to ensure real eigenvalues.
- The appendix provides formulas for the eigenvalues of the viscous subsystem matrices A_k when heat conduction is ignored.
- The notation lists variables and matrix components used in the viscous subsystem eigenvalue analysis.
- The eigenvalues are expressed as v_k − λ_3, v_k − λ_2, v_k − λ_1, v_k, v_k, v_k + λ_1, v_k + λ_2, and v_k + λ_3.
- Assuming a convex equation of state makes the HPR model symmetric hyperbolic, so the eigenvalues are assumed real and obtained from a cubic characteristic polynomial.