Cont damage fix

docs/documentation/case.md

@@ -1149,7 +1149,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.
 

docs/documentation/equations.md

@@ -572,6 +572,26 @@ 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`
+
+A scalar damage field \f$D \in [0,1]\f$ is transported with the damageable-solid partial mass
+\f$m_s = \sum_{i:\,G_i > 0} \alpha_i \rho_i\f$:
+
+\f[\frac{\partial (m_s D)}{\partial t} + \nabla \cdot (m_s D\, \mathbf{u}) = m_s\,\dot{D}\f]
+
+Damage grows when the maximum principal Cauchy stress
+\f$\sigma_1 = \lambda_{\max}(-p\mathbf{I} + \boldsymbol{\tau}^e)\f$ exceeds \f$\tau^*\f$:
+
+\f[\dot{D} = \bigl(\bar{\alpha}\,\max(\sigma_1 - \tau^*,\, 0)\bigr)^{s}\f]
+
+The damaged shear modulus is
+
+\f[G = G_0(1-D),\f]
+
+and is used in the elastic stress evolution and HLL/HLLC wave speeds; elastic energy uses the undamaged modulus \f$G_0\f$.
+
 ## 8. Phase Change (`relax = .true.`) (\cite Wilfong26 Sec. 4.1.3)
 
 **Source:** `src/common/m_phase_change.fpp`

docs/references.bib

@@ -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}
+}

examples/2D_cont_damage/case.py

@@ -1,4 +1,16 @@
 #!/usr/bin/env python3
+# 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
@@ -9,18 +21,18 @@
             "run_time_info": "T",
             # Computational Domain Parameters
             "x_domain%beg": 0.0,
-            "x_domain%end": 0.001,
+            "x_domain%end": 0.01,
             "y_domain%beg": 0.0,
-            "y_domain%end": 0.0005,
-            "m": 50,
-            "n": 25,
+            "y_domain%end": 0.01,
+            "m": 199,
+            "n": 199,
             "p": 0,
-            "dt": 2e-12,
+            "dt": 5.0e-9,
             "t_step_start": 0,
-            "t_step_stop": 40000,
-            "t_step_save": 2000,
+            "t_step_stop": 1200,
+            "t_step_save": 100,
             # Simulation Algorithm Parameters
-            "num_patches": 2,
+            "num_patches": 3,
             "model_eqns": "5eq",
             "alt_soundspeed": "F",
             "num_fluids": 2,
@@ -29,76 +41,86 @@
             "time_stepper": "rk3",
             "weno_order": 5,
             "weno_eps": 1.0e-16,
-            "teno": "T",
-            "teno_CT": 1e-8,
+            "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": -2,
+            "bc_y%beg": -6,
             "bc_y%end": -6,
-            # Hypoelasticity
+            # Hypoelasticity + continuum damage
             "hypoelasticity": "T",
             "fd_order": 4,
             "cont_damage": "T",
-            "tau_star": 0.0,
+            "tau_star": 2.0e7,
             "cont_damage_s": 2.0,
-            "alpha_bar": 1e-4,
+            "alpha_bar": 2.0e-5,
             # Formatted Database Files Structure Parameters
             "format": "silo",
             "precision": "double",
             "prim_vars_wrt": "T",
-            "parallel_io": "T",
-            # Patch 1 Liquid
+            "parallel_io": "F",
+            # Background water (ambient)
             "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)%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": 1e05,
-            "patch_icpp(1)%alpha_rho(1)": 1100 * (1.0 - 1e-6),
+            "patch_icpp(1)%pres": 1.0e5,
+            "patch_icpp(1)%alpha_rho(1)": 1000.0 * (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_rho(2)": 1000.0 * 1e-6,
             "patch_icpp(1)%alpha(2)": 1e-6,
-            # Patch 2 Solid
+            "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)%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)%x_centroid": 0.00175,
+            "patch_icpp(2)%y_centroid": 0.005,
             "patch_icpp(2)%length_x": 0.0005,
-            "patch_icpp(2)%length_y": 0.00025,
+            "patch_icpp(2)%length_y": 0.01,
             "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
+            "patch_icpp(2)%pres": 4.0e8,
+            "patch_icpp(2)%alpha_rho(1)": 1000.0 * (1.0 - 1e-6),
+            "patch_icpp(2)%alpha(1)": 1.0 - 1e-6,
+            "patch_icpp(2)%alpha_rho(2)": 1000.0 * 1e-6,
+            "patch_icpp(2)%alpha(2)": 1e-6,
+            "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)%alter_patch(1)": "T",
+            "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)%vel(1)": 0.0,
+            "patch_icpp(3)%vel(2)": 0.0,
+            "patch_icpp(3)%pres": 1.0e5,
+            "patch_icpp(3)%alpha_rho(1)": 1000.0 * 1e-6,
+            "patch_icpp(3)%alpha(1)": 1e-6,
+            "patch_icpp(3)%alpha_rho(2)": 1000.0 * (1.0 - 1e-6),
+            "patch_icpp(3)%alpha(2)": 1.0 - 1e-6,
+            "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.0e00 / (4.4e00 - 1.0e00),
-            "fluid_pp(1)%pi_inf": 4.4e00 * 5.57e08 / (4.4e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (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,
+            "fluid_pp(2)%G": 1.0e9,
         }
     )
 )