Skip to content

Cont damage fix #1723

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
ChrisZYJ wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

Cont damage fix #1723

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2039d6a
Fix continuum damage: primitive principal-Cauchy driver, solid-mass t…
ChrisZYJ Aug 11, 2026
27e9460
Validate continuum damage inputs and require hypoelasticity
ChrisZYJ Aug 11, 2026
b69d337
Make damage goldens threshold-discriminating and add an axisymmetric …
ChrisZYJ Aug 11, 2026
6bffc9e
Document the continuum damage model and rework the damage examples
ChrisZYJ Aug 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/documentation/case.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1144,7 +1144,7 @@ Note: For relativistic flow, the conservative and primitive densities are differ
| `rburn%%n` | Real | Reactive-burn pressure-drive exponent |
| `rburn%%ta` | Real | Reactive-burn activation temperature [K] (0 = off) |

- `cont_damage` activates continuum damage model for solid materials. Requires `tau_star`, `cont_damage_s`, and `alpha_bar` to be set (empirically determined) (\cite Cao19).
- `cont_damage` activates the continuum damage model for hypoelastic solid materials (requires `hypoelasticity = T`; HLL/HLLC only). Damage is produced by tensile maximum principal Cauchy stress beyond `tau_star` (\f$\geq 0\f$) at rate `(alpha_bar*(sigma_1 - tau_star))**cont_damage_s` and is transported with the damageable-solid partial mass; see @ref equations for the model statement (\cite Cao19; \cite Spratt24). `tau_star`, `cont_damage_s` (\f$> 0\f$), and `alpha_bar` (\f$\geq 0\f$) are empirically determined.

- `reactive_burn` converts a "reactant" fluid into a "product" fluid (`num_fluids = 2`, ``chemistry = 'F'``) via a programmed pressure burn `dlambda/dt = rburn%%k (1 - lambda) ((p - rburn%%pign)/rburn%%pref)^rburn%%n`. The two fluids share the same `gamma`/`pi_inf` and differ only in `qv`, so the conversion releases `qv` through the mixture EOS — a reactive-Euler/ZND detonation model on the diffuse-interface framework. It runs on the 5-equation (`model_eqns = 2`) and 6-equation (`model_eqns = 3`) multi-fluid models. Setting `rburn%%ta > 0` multiplies the rate by an Arrhenius factor `exp(-rburn%%ta/T)`, where `T` is the reactant phasic temperature, giving temperature-driven ignition instead of a pure pressure switch.

Expand Down
42 changes: 42 additions & 0 deletions docs/documentation/equations.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -572,6 +572,48 @@ where \f$\mathbf{l} = \nabla \mathbf{u}\f$ is the velocity gradient and \f$\math

This adds 6 additional transport equations in 3D (symmetric stress tensor: \f$\tau_{xx}^e, \tau_{xy}^e, \tau_{yy}^e, \tau_{xz}^e, \tau_{yz}^e, \tau_{zz}^e\f$).

### 7.2 Continuum Damage (`cont_damage = .true.`) (\cite Cao19; \cite Spratt24 Sec. 4.1.2)

**Source:** `src/simulation/m_hypoelastic.fpp` (damage source), Riemann solvers (transport)

A scalar damage field \f$D \in [0,1]\f$ weakens the shear response of hypoelastic material.
The model is a two-way hybrid of the hypoelastic formulation above with the
Tuler–Butcher-type overstress law of \cite Cao19, in the shear-modulus form derived in
\cite Spratt24: it is not a literal reproduction of Cao's constant-Poisson-ratio
finite-element model (no element erosion, and the volumetric response remains governed
by the equation of state).

**Damage transport.** Damage is carried by the damageable-solid partial mass
\f$m_s = \sum_{i:\,G_i > 0} \alpha_i \rho_i\f$ (a fluid is damageable if its `fluid_pp%%G`
is positive):

\f[\frac{\partial (m_s D)}{\partial t} + \nabla \cdot (m_s D\, \mathbf{u}) = m_s\,\dot{D}\f]

so damage moves with the solid material, pure fluid can neither produce nor carry
damage mass, and \f$D = m_s D / m_s\f$ is recovered clamped to \f$[0,1]\f$.

**Damage rate.** Production is driven by the maximum principal value of the total
Cauchy stress \f$\sigma = -p\,\mathbf{I} + \boldsymbol{\tau}^e\f$, evaluated in the full
three-dimensional sense in every dimensionality (out-of-plane components follow from
the traceless closure in 1D and plane 2D; the stored hoop stress enters in axisymmetry):

\f[\dot{D} = \bigl(\bar{\alpha}\,\max(\sigma_1 - \tau^*,\, 0)\bigr)^{s}\f]

Only tension beyond the threshold \f$\tau^*\f$ produces damage; production stops once
\f$D = 1\f$. Note that 2D planar simulations are plane-strain by construction, so the
out-of-plane constraint stress participates in \f$\sigma_1\f$.

**Two-way coupling.** The damaged shear modulus \f$G = G_0(1 - D)\f$ enters the elastic
stress evolution and the Riemann wave speeds (tangent-modulus interpretation: damage
weakens subsequent stress growth but does not relax previously accumulated stress).
The elastic energy in \f$E\f$ is booked with the undamaged modulus \f$G_0\f$, consistent
with the tangent interpretation and bounded as \f$D \to 1\f$. The elastic energy released
by softening is not separately tracked, and there is no failure/erosion mechanism: a
fully damaged region behaves as a fluid with residual (frozen) deviatoric stress. This
is a bounded damage-only model.

Supported with the HLL and HLLC Riemann solvers.

## 8. Phase Change (`relax = .true.`) (\cite Wilfong26 Sec. 4.1.3)

**Source:** `src/common/m_phase_change.fpp`
Expand Down
8 changes: 8 additions & 0 deletions docs/references.bib
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -714,3 +714,11 @@ @article{Magnaudet2000
year = {2000},
doi = {10.1146/annurev.fluid.32.1.659}
}

@phdthesis{Spratt24,
author = {J.-S. A. Spratt},
title = {Numerical simulations of cavitating bubbles in elastic and viscoelastic materials for biomedical applications},
school = {California Institute of Technology},
year = {2024},
doi = {10.7907/g34e-6p65}
}
94 changes: 0 additions & 94 deletions examples/1D_cont_damage/case.py
View file
Open in desktop

This file was deleted.

228 changes: 126 additions & 102 deletions examples/2D_cont_damage/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,104 +1,128 @@
#!/usr/bin/env python3
#!/usr/bin/python
# 2D continuum-damage demonstration: solid disk in water hit by a planar pulse.
#
# A planar high-pressure strip in the water launches a pulse that diffracts around
# and transmits through a damageable solid disk. Tension concentrations (poles of
# the disk during passage, interior wave focusing afterwards) exceed tau_star and
# accumulate damage; the surrounding water stays undamaged (damage is carried by
# the solid partial mass, U_D = m_s D).
#
# Parameters are deliberately sub-critical: with the tangent model, stress is
# retained as D grows, and driving D all the way to 1 under sustained tension
# leaves a zero-shear-stiffness cell with residual stress (a documented model
# limitation). tau_star and alpha_bar here keep max D around 0.2.
import json

# Configuring case dictionary
print(
json.dumps(
{
# Logistics
"run_time_info": "T",
# Computational Domain Parameters
"x_domain%beg": 0.0,
"x_domain%end": 0.001,
"y_domain%beg": 0.0,
"y_domain%end": 0.0005,
"m": 50,
"n": 25,
"p": 0,
"dt": 2e-12,
"t_step_start": 0,
"t_step_stop": 40000,
"t_step_save": 2000,
# Simulation Algorithm Parameters
"num_patches": 2,
"model_eqns": "5eq",
"alt_soundspeed": "F",
"num_fluids": 2,
"mpp_lim": "F",
"mixture_err": "F",
"time_stepper": "rk3",
"weno_order": 5,
"weno_eps": 1.0e-16,
"teno": "T",
"teno_CT": 1e-8,
"null_weights": "F",
"mp_weno": "F",
"riemann_solver": "hll",
"wave_speeds": "direct",
"avg_state": "arithmetic",
"bc_x%beg": -6,
"bc_x%end": -6,
"bc_y%beg": -2,
"bc_y%end": -6,
# Hypoelasticity
"hypoelasticity": "T",
"fd_order": 4,
"cont_damage": "T",
"tau_star": 0.0,
"cont_damage_s": 2.0,
"alpha_bar": 1e-4,
# Formatted Database Files Structure Parameters
"format": "silo",
"precision": "double",
"prim_vars_wrt": "T",
"parallel_io": "T",
# Patch 1 Liquid
"patch_icpp(1)%geometry": 3,
"patch_icpp(1)%x_centroid": 0.0005,
"patch_icpp(1)%y_centroid": 0.00025,
"patch_icpp(1)%length_x": 0.001,
"patch_icpp(1)%length_y": 0.0005,
"patch_icpp(1)%vel(1)": 0.0,
"patch_icpp(1)%vel(2)": 0.0,
"patch_icpp(1)%pres": 1e05,
"patch_icpp(1)%alpha_rho(1)": 1100 * (1.0 - 1e-6),
"patch_icpp(1)%alpha(1)": 1.0 - 1e-6,
"patch_icpp(1)%alpha_rho(2)": 1100 * 1e-6,
"patch_icpp(1)%alpha(2)": 1e-6,
# Patch 2 Solid
"patch_icpp(2)%alter_patch(1)": "T",
"patch_icpp(2)%geometry": 3,
"patch_icpp(2)%x_centroid": 0.0005,
"patch_icpp(2)%y_centroid": 0.000125,
"patch_icpp(2)%length_x": 0.0005,
"patch_icpp(2)%length_y": 0.00025,
"patch_icpp(2)%vel(1)": 0.0,
"patch_icpp(2)%vel(2)": 0.0,
"patch_icpp(2)%pres": 1e05,
"patch_icpp(2)%alpha_rho(1)": 1100 * 1e-6,
"patch_icpp(2)%alpha(1)": 1e-6,
"patch_icpp(2)%alpha_rho(2)": 1100 * (1.0 - 1e-6),
"patch_icpp(2)%alpha(2)": 1.0 - 1e-6,
# Acoustic source
"acoustic_source": "T",
"num_source": 1,
"acoustic(1)%support": 5,
"acoustic(1)%loc(1)": 0.00005,
"acoustic(1)%loc(2)": 0.0,
"acoustic(1)%pulse": 1,
"acoustic(1)%npulse": 999,
"acoustic(1)%mag": 100.0,
"acoustic(1)%wavelength": 0.0001,
"acoustic(1)%foc_length": 0.00045,
"acoustic(1)%aperture": 0.0008,
"acoustic(1)%delay": 0.0,
# Fluids Physical Parameters
"fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
"fluid_pp(1)%pi_inf": 4.4e00 * 5.57e08 / (4.4e00 - 1.0e00),
"fluid_pp(1)%G": 0.0,
"fluid_pp(2)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
"fluid_pp(2)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
"fluid_pp(2)%G": 1.0e09,
}
)
)
GAMMA = 4.4
PI_B = 6.0e8
RHO = 1000.0
EPS = 1e-6

case = {
# Logistics
"run_time_info": "T",
# Computational Domain Parameters
"x_domain%beg": 0.0,
"x_domain%end": 0.01,
"y_domain%beg": 0.0,
"y_domain%end": 0.01,
"m": 199,
"n": 199,
"p": 0,
"dt": 5.0e-9,
"t_step_start": 0,
"t_step_stop": 1200,
"t_step_save": 100,
# Simulation Algorithm Parameters
"num_patches": 3,
"model_eqns": "5eq",
"alt_soundspeed": "F",
"num_fluids": 2,
"mpp_lim": "F",
"mixture_err": "F",
"time_stepper": "rk3",
"weno_order": 5,
"weno_eps": 1.0e-16,
"weno_Re_flux": "F",
"weno_avg": "F",
"mapped_weno": "T",
"null_weights": "F",
"mp_weno": "F",
"riemann_solver": "hll",
"wave_speeds": "direct",
"avg_state": "arithmetic",
"bc_x%beg": -6,
"bc_x%end": -6,
"bc_y%beg": -6,
"bc_y%end": -6,
# Hypoelasticity + continuum damage
"hypoelasticity": "T",
"fd_order": 4,
"cont_damage": "T",
"tau_star": 2.0e7,
"cont_damage_s": 2.0,
"alpha_bar": 2.0e-5,
# Formatted Database Files Structure Parameters
"format": "silo",
"precision": "double",
"prim_vars_wrt": "T",
"parallel_io": "F",
# Background water (ambient)
"patch_icpp(1)%geometry": 3,
"patch_icpp(1)%x_centroid": 0.005,
"patch_icpp(1)%y_centroid": 0.005,
"patch_icpp(1)%length_x": 0.01,
"patch_icpp(1)%length_y": 0.01,
"patch_icpp(1)%vel(1)": 0.0,
"patch_icpp(1)%vel(2)": 0.0,
"patch_icpp(1)%pres": 1.0e5,
"patch_icpp(1)%alpha_rho(1)": RHO * (1 - EPS),
"patch_icpp(1)%alpha_rho(2)": RHO * EPS,
"patch_icpp(1)%alpha(1)": 1 - EPS,
"patch_icpp(1)%alpha(2)": EPS,
"patch_icpp(1)%tau_e(1)": 0.0,
"patch_icpp(1)%tau_e(2)": 0.0,
"patch_icpp(1)%tau_e(3)": 0.0,
# Driver strip (interior so the open boundaries never act as a reservoir)
"patch_icpp(2)%geometry": 3,
"patch_icpp(2)%x_centroid": 0.00175,
"patch_icpp(2)%y_centroid": 0.005,
"patch_icpp(2)%length_x": 0.0005,
"patch_icpp(2)%length_y": 0.01,
"patch_icpp(2)%alter_patch(1)": "T",
"patch_icpp(2)%vel(1)": 0.0,
"patch_icpp(2)%vel(2)": 0.0,
"patch_icpp(2)%pres": 4.0e8,
"patch_icpp(2)%alpha_rho(1)": RHO * (1 - EPS),
"patch_icpp(2)%alpha_rho(2)": RHO * EPS,
"patch_icpp(2)%alpha(1)": 1 - EPS,
"patch_icpp(2)%alpha(2)": EPS,
"patch_icpp(2)%tau_e(1)": 0.0,
"patch_icpp(2)%tau_e(2)": 0.0,
"patch_icpp(2)%tau_e(3)": 0.0,
# Damageable solid disk
"patch_icpp(3)%geometry": 2,
"patch_icpp(3)%x_centroid": 0.005,
"patch_icpp(3)%y_centroid": 0.005,
"patch_icpp(3)%radius": 0.0015,
"patch_icpp(3)%alter_patch(1)": "T",
"patch_icpp(3)%vel(1)": 0.0,
"patch_icpp(3)%vel(2)": 0.0,
"patch_icpp(3)%pres": 1.0e5,
"patch_icpp(3)%alpha_rho(1)": RHO * EPS,
"patch_icpp(3)%alpha_rho(2)": RHO * (1 - EPS),
"patch_icpp(3)%alpha(1)": EPS,
"patch_icpp(3)%alpha(2)": 1 - EPS,
"patch_icpp(3)%tau_e(1)": 0.0,
"patch_icpp(3)%tau_e(2)": 0.0,
"patch_icpp(3)%tau_e(3)": 0.0,
# Fluids: 1 = water (no shear stiffness), 2 = damageable solid
"fluid_pp(1)%gamma": 1.0 / (GAMMA - 1.0),
"fluid_pp(1)%pi_inf": GAMMA * PI_B / (GAMMA - 1.0),
"fluid_pp(1)%G": 0.0,
"fluid_pp(2)%gamma": 1.0 / (GAMMA - 1.0),
"fluid_pp(2)%pi_inf": GAMMA * PI_B / (GAMMA - 1.0),
"fluid_pp(2)%G": 1.0e9,
}

print(json.dumps(case))
Loading
Loading