Skip to content

Commit 488175e

Browse files
ChrisZYJsbryngelson
ChrisZYJ
and
sbryngelson
authored
Cont damage fix (#1723)
Co-authored-by: Spencer Bryngelson <sbryngelson@gmail.com>
1 parent 7f49b7b commit 488175e

30 files changed

Lines changed: 1065 additions & 861 deletions

‎docs/documentation/case.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1152,7 +1152,7 @@ Note: For relativistic flow, the conservative and primitive densities are differ
11521152
| `rburn%%n` | Real | Reactive-burn pressure-drive exponent |
11531153
| `rburn%%ta` | Real | Reactive-burn activation temperature [K] (0 = off) |
11541154

1155-
- `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).
1155+
- `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.
11561156

11571157
- `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.
11581158

‎docs/documentation/equations.md‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -572,6 +572,26 @@ where \f$\mathbf{l} = \nabla \mathbf{u}\f$ is the velocity gradient and \f$\math
572572

573573
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$).
574574

575+
### 7.2 Continuum Damage (`cont_damage = .true.`) (\cite Cao19; \cite Spratt24 Sec. 4.1.2)
576+
577+
**Source:** `src/simulation/m_hypoelastic.fpp`
578+
579+
A scalar damage field \f$D \in [0,1]\f$ is transported with the damageable-solid partial mass
580+
\f$m_s = \sum_{i:\,G_i > 0} \alpha_i \rho_i\f$:
581+
582+
\f[\frac{\partial (m_s D)}{\partial t} + \nabla \cdot (m_s D\, \mathbf{u}) = m_s\,\dot{D}\f]
583+
584+
Damage grows when the maximum principal Cauchy stress
585+
\f$\sigma_1 = \lambda_{\max}(-p\mathbf{I} + \boldsymbol{\tau}^e)\f$ exceeds \f$\tau^*\f$:
586+
587+
\f[\dot{D} = \bigl(\bar{\alpha}\,\max(\sigma_1 - \tau^*,\, 0)\bigr)^{s}\f]
588+
589+
The damaged shear modulus is
590+
591+
\f[G = G_0(1-D),\f]
592+
593+
and is used in the elastic stress evolution and HLL/HLLC wave speeds; elastic energy uses the undamaged modulus \f$G_0\f$.
594+
575595
## 8. Phase Change (`relax = .true.`) (\cite Wilfong26 Sec. 4.1.3)
576596

577597
**Source:** `src/common/m_phase_change.fpp`

‎docs/references.bib‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -714,3 +714,11 @@ @article{Magnaudet2000
714714
year = {2000},
715715
doi = {10.1146/annurev.fluid.32.1.659}
716716
}
717+
718+
@phdthesis{Spratt24,
719+
author = {J.-S. A. Spratt},
720+
title = {Numerical simulations of cavitating bubbles in elastic and viscoelastic materials for biomedical applications},
721+
school = {California Institute of Technology},
722+
year = {2024},
723+
doi = {10.7907/g34e-6p65}
724+
}

‎examples/1D_cont_damage/case.py‎

Lines changed: 0 additions & 96 deletions
This file was deleted.

‎examples/2D_cont_damage/case.py‎

Lines changed: 70 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,16 @@
11
#!/usr/bin/env python3
2+
# 2D continuum-damage demonstration: solid disk in water hit by a planar pulse.
3+
#
4+
# A planar high-pressure strip in the water launches a pulse that diffracts around
5+
# and transmits through a damageable solid disk. Tension concentrations (poles of
6+
# the disk during passage, interior wave focusing afterwards) exceed tau_star and
7+
# accumulate damage; the surrounding water stays undamaged (damage is carried by
8+
# the solid partial mass, U_D = m_s D).
9+
#
10+
# Parameters are deliberately sub-critical: with the tangent model, stress is
11+
# retained as D grows, and driving D all the way to 1 under sustained tension
12+
# leaves a zero-shear-stiffness cell with residual stress (a documented model
13+
# limitation). tau_star and alpha_bar here keep max D around 0.2.
214
import json
315

416
# Configuring case dictionary
@@ -9,18 +21,18 @@
921
"run_time_info": "T",
1022
# Computational Domain Parameters
1123
"x_domain%beg": 0.0,
12-
"x_domain%end": 0.001,
24+
"x_domain%end": 0.01,
1325
"y_domain%beg": 0.0,
14-
"y_domain%end": 0.0005,
15-
"m": 50,
16-
"n": 25,
26+
"y_domain%end": 0.01,
27+
"m": 199,
28+
"n": 199,
1729
"p": 0,
18-
"dt": 2e-12,
30+
"dt": 5.0e-9,
1931
"t_step_start": 0,
20-
"t_step_stop": 40000,
21-
"t_step_save": 2000,
32+
"t_step_stop": 1200,
33+
"t_step_save": 100,
2234
# Simulation Algorithm Parameters
23-
"num_patches": 2,
35+
"num_patches": 3,
2436
"model_eqns": "5eq",
2537
"alt_soundspeed": "F",
2638
"num_fluids": 2,
@@ -29,78 +41,88 @@
2941
"time_stepper": "rk3",
3042
"weno_order": 5,
3143
"weno_eps": 1.0e-16,
32-
"teno": "T",
33-
"teno_CT": 1e-8,
44+
"weno_Re_flux": "F",
45+
"weno_avg": "F",
46+
"mapped_weno": "T",
3447
"null_weights": "F",
3548
"mp_weno": "F",
3649
"riemann_solver": "hll",
3750
"wave_speeds": "direct",
3851
"avg_state": "arithmetic",
3952
"bc_x%beg": -6,
4053
"bc_x%end": -6,
41-
"bc_y%beg": -2,
54+
"bc_y%beg": -6,
4255
"bc_y%end": -6,
43-
# Hypoelasticity
56+
# Hypoelasticity + continuum damage
4457
"hypoelasticity": "T",
4558
"fd_order": 4,
4659
"cont_damage": "T",
47-
"tau_star": 0.0,
60+
"tau_star": 2.0e7,
4861
"cont_damage_s": 2.0,
49-
"alpha_bar": 1e-4,
62+
"alpha_bar": 2.0e-5,
5063
# Formatted Database Files Structure Parameters
5164
"format": "silo",
5265
"precision": "double",
5366
"prim_vars_wrt": "T",
54-
"parallel_io": "T",
55-
# Patch 1 Liquid
67+
"parallel_io": "F",
68+
# Background water (ambient)
5669
"patch_icpp(1)%geometry": 3,
57-
"patch_icpp(1)%x_centroid": 0.0005,
58-
"patch_icpp(1)%y_centroid": 0.00025,
59-
"patch_icpp(1)%length_x": 0.001,
60-
"patch_icpp(1)%length_y": 0.0005,
70+
"patch_icpp(1)%x_centroid": 0.005,
71+
"patch_icpp(1)%y_centroid": 0.005,
72+
"patch_icpp(1)%length_x": 0.01,
73+
"patch_icpp(1)%length_y": 0.01,
6174
"patch_icpp(1)%vel(1)": 0.0,
6275
"patch_icpp(1)%vel(2)": 0.0,
63-
"patch_icpp(1)%pres": 1e05,
64-
"patch_icpp(1)%alpha_rho(1)": 1100 * (1.0 - 1e-6),
76+
"patch_icpp(1)%pres": 1.0e5,
77+
"patch_icpp(1)%alpha_rho(1)": 1000.0 * (1.0 - 1e-6),
6578
"patch_icpp(1)%alpha(1)": 1.0 - 1e-6,
66-
"patch_icpp(1)%alpha_rho(2)": 1100 * 1e-6,
79+
"patch_icpp(1)%alpha_rho(2)": 1000.0 * 1e-6,
6780
"patch_icpp(1)%alpha(2)": 1e-6,
68-
# Patch 2 Solid
81+
"patch_icpp(1)%tau_e(1)": 0.0,
82+
"patch_icpp(1)%tau_e(2)": 0.0,
83+
"patch_icpp(1)%tau_e(3)": 0.0,
84+
# Driver strip (interior so the open boundaries never act as a reservoir)
6985
"patch_icpp(2)%alter_patch(1)": "T",
7086
"patch_icpp(2)%geometry": 3,
71-
"patch_icpp(2)%x_centroid": 0.0005,
72-
"patch_icpp(2)%y_centroid": 0.000125,
87+
"patch_icpp(2)%x_centroid": 0.00175,
88+
"patch_icpp(2)%y_centroid": 0.005,
7389
"patch_icpp(2)%length_x": 0.0005,
74-
"patch_icpp(2)%length_y": 0.00025,
90+
"patch_icpp(2)%length_y": 0.01,
7591
"patch_icpp(2)%vel(1)": 0.0,
7692
"patch_icpp(2)%vel(2)": 0.0,
77-
"patch_icpp(2)%pres": 1e05,
78-
"patch_icpp(2)%alpha_rho(1)": 1100 * 1e-6,
79-
"patch_icpp(2)%alpha(1)": 1e-6,
80-
"patch_icpp(2)%alpha_rho(2)": 1100 * (1.0 - 1e-6),
81-
"patch_icpp(2)%alpha(2)": 1.0 - 1e-6,
82-
# Acoustic source
83-
"acoustic_source": "T",
84-
"num_source": 1,
85-
"acoustic(1)%support": 5,
86-
"acoustic(1)%loc(1)": 0.00005,
87-
"acoustic(1)%loc(2)": 0.0,
88-
"acoustic(1)%pulse": 1,
89-
"acoustic(1)%npulse": 999,
90-
"acoustic(1)%mag": 100.0,
91-
"acoustic(1)%wavelength": 0.0001,
92-
"acoustic(1)%foc_length": 0.00045,
93-
"acoustic(1)%aperture": 0.0008,
94-
"acoustic(1)%delay": 0.0,
95-
# Fluids Physical Parameters
93+
"patch_icpp(2)%pres": 4.0e8,
94+
"patch_icpp(2)%alpha_rho(1)": 1000.0 * (1.0 - 1e-6),
95+
"patch_icpp(2)%alpha(1)": 1.0 - 1e-6,
96+
"patch_icpp(2)%alpha_rho(2)": 1000.0 * 1e-6,
97+
"patch_icpp(2)%alpha(2)": 1e-6,
98+
"patch_icpp(2)%tau_e(1)": 0.0,
99+
"patch_icpp(2)%tau_e(2)": 0.0,
100+
"patch_icpp(2)%tau_e(3)": 0.0,
101+
# Damageable solid disk
102+
"patch_icpp(3)%alter_patch(1)": "T",
103+
"patch_icpp(3)%geometry": 2,
104+
"patch_icpp(3)%x_centroid": 0.005,
105+
"patch_icpp(3)%y_centroid": 0.005,
106+
"patch_icpp(3)%radius": 0.0015,
107+
"patch_icpp(3)%vel(1)": 0.0,
108+
"patch_icpp(3)%vel(2)": 0.0,
109+
"patch_icpp(3)%pres": 1.0e5,
110+
"patch_icpp(3)%alpha_rho(1)": 1000.0 * 1e-6,
111+
"patch_icpp(3)%alpha(1)": 1e-6,
112+
"patch_icpp(3)%alpha_rho(2)": 1000.0 * (1.0 - 1e-6),
113+
"patch_icpp(3)%alpha(2)": 1.0 - 1e-6,
114+
"patch_icpp(3)%tau_e(1)": 0.0,
115+
"patch_icpp(3)%tau_e(2)": 0.0,
116+
"patch_icpp(3)%tau_e(3)": 0.0,
117+
# Fluids: 1 = water (no shear stiffness), 2 = damageable solid
96118
"fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
97119
"fluid_pp(1)%eos": "stiffened_gas",
98-
"fluid_pp(1)%pi_inf": 4.4e00 * 5.57e08 / (4.4e00 - 1.0e00),
120+
"fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
99121
"fluid_pp(1)%G": 0.0,
100122
"fluid_pp(2)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
101123
"fluid_pp(2)%eos": "stiffened_gas",
102124
"fluid_pp(2)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
103-
"fluid_pp(2)%G": 1.0e09,
125+
"fluid_pp(2)%G": 1.0e9,
104126
}
105127
)
106128
)

0 commit comments

Comments
 (0)