Skip to content

Commit 068da2c

Browse files
ChrisZYJChrisZYJ
authored
Continuum Damage Model (#816)
1 parent a2a8800 commit 068da2c

36 files changed

+1475
-59
lines changed

Diff for: ‎docs/documentation/case.md

+45-32
Large diffs are not rendered by default.

Diff for: ‎docs/documentation/references.md

+3-1
Original file line numberDiff line numberDiff line change
@@ -60,4 +60,6 @@
6060

6161
- <a id="Miyoshi05">Miyoishi, T., and Kusano, K. (2005). A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. Journal of Computational Physics, 208(1), 315-344.</a>
6262

63-
- <a id="Powell94">Powell, K. G. (1994). An approximate Riemann solver for magnetohydrodynamics: (That works in more than one dimension). In Upwind and high-resolution schemes (pp. 570-583). Springer.</a>
63+
- <a id="Powell94">Powell, K. G. (1994). An approximate Riemann solver for magnetohydrodynamics: (That works in more than one dimension). In Upwind and high-resolution schemes (pp. 570-583). Springer.</a>
64+
65+
- <a id="Cao19">Cao, S., Zhang, Y., Liao, D., Zhong, P., and Wang, K. G. (2019). Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid. International Journal of Solids and Structures, 169:55–71. Elsevier.</a>

Diff for: ‎docs/documentation/visualization.md

+2-2
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
# Flow visualization
22

33
A post-processed database in Silo-HDF5 format can be visualized and analyzed using Paraview and VisIt.
4-
After the post-processing of simulation data (see section [Running](running.md#running-1)), a directory named `silo_hdf5` contains a silo-HDF5 database.
4+
After the post-processing of simulation data (see section [Running](running.md)), a directory named `silo_hdf5` contains a silo-HDF5 database.
55
Here, `silo_hdf5/` includes a directory named `root/` that contains index files for flow field data at each saved time step.
66

77
### Visualizing with Paraview
@@ -38,7 +38,7 @@ For many cases, this step will require resizing the render view window.
3838
VisIt is an alternative open-source interactive parallel visualization and graphical analysis tool for viewing scientific data.
3939
Versions of VisIt after 2.6.0 have been confirmed to work with the MFC databases for some parallel environments.
4040
Nevertheless, installation and configuration of VisIt can be environment-dependent and are left to the user.
41-
Further remarks on parallel flow visualization, analysis, and processing of the MFC database using VisIt can also be found in [Coralic (2015)](references.md#Coralic15) and [Meng (2016)](references.md#Meng16).
41+
Further remarks on parallel flow visualization, analysis, and processing of the MFC database using VisIt can also be found in [Coralic (2015)](references.md) and [Meng (2016)](references.md).
4242

4343
The user can launch VisIt and open the index files under `/silo_hdf5/root`.
4444
Once the database is loaded, flow field variables contained in the database can be added to the plot.

Diff for: ‎examples/1D_cont_damage/case.py

+94
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
#!/usr/bin/python
2+
import json
3+
4+
# Configuring case dictionary
5+
print(
6+
json.dumps(
7+
{
8+
# Logistics
9+
"run_time_info": "T",
10+
# Computational Domain Parameters
11+
"x_domain%beg": 0.0e00,
12+
"x_domain%end": 1.0e00,
13+
"m": 100,
14+
"n": 0,
15+
"p": 0,
16+
"dt": 5e-10,
17+
"t_step_start": 0,
18+
"t_step_stop": 50000,
19+
"t_step_save": 50000,
20+
# Simulation Algorithm Parameters
21+
"num_patches": 2,
22+
"model_eqns": 2,
23+
"alt_soundspeed": "F",
24+
"num_fluids": 2,
25+
"mpp_lim": "F",
26+
"mixture_err": "F",
27+
"time_stepper": 3,
28+
"weno_order": 5,
29+
"weno_eps": 1.0e-16,
30+
"weno_Re_flux": "F",
31+
"weno_avg": "F",
32+
"mapped_weno": "T",
33+
"null_weights": "F",
34+
"mp_weno": "F",
35+
"riemann_solver": 1,
36+
"wave_speeds": 1,
37+
"avg_state": 2,
38+
"bc_x%beg": -3,
39+
"bc_x%end": -3,
40+
# Turning on Hypoelasticity
41+
"hypoelasticity": "T",
42+
"fd_order": 4,
43+
"cont_damage": "T",
44+
"tau_star": 0.0,
45+
"cont_damage_s": 2.0,
46+
"alpha_bar": 1e-4,
47+
# Acoustic Source
48+
"acoustic_source": "T",
49+
"num_source": 1,
50+
"acoustic(1)%support": 1,
51+
"acoustic(1)%loc(1)": 0.1,
52+
"acoustic(1)%pulse": 1,
53+
"acoustic(1)%npulse": 999,
54+
"acoustic(1)%dir": 1.0,
55+
"acoustic(1)%mag": 1000.0,
56+
"acoustic(1)%frequency": 1e4,
57+
"acoustic(1)%delay": 0,
58+
# Formatted Database Files Structure Parameters
59+
"format": 1,
60+
"precision": 2,
61+
"prim_vars_wrt": "T",
62+
"parallel_io": "F",
63+
# Patch 1 L
64+
"patch_icpp(1)%geometry": 1,
65+
"patch_icpp(1)%x_centroid": 0.25,
66+
"patch_icpp(1)%length_x": 0.5,
67+
"patch_icpp(1)%vel(1)": 0.0,
68+
"patch_icpp(1)%pres": 1e5,
69+
"patch_icpp(1)%alpha_rho(1)": 1000 * (1.0 - 1e-6),
70+
"patch_icpp(1)%alpha_rho(2)": 1000 * 1e-6,
71+
"patch_icpp(1)%alpha(1)": 1.0 - 1e-6,
72+
"patch_icpp(1)%alpha(2)": 1e-6,
73+
"patch_icpp(1)%tau_e(1)": 0.0,
74+
# Patch 2 R
75+
"patch_icpp(2)%geometry": 1,
76+
"patch_icpp(2)%x_centroid": 0.75,
77+
"patch_icpp(2)%length_x": 0.5,
78+
"patch_icpp(2)%vel(1)": 0.0,
79+
"patch_icpp(2)%pres": 1e5,
80+
"patch_icpp(2)%alpha_rho(1)": 1000 * 1e-6,
81+
"patch_icpp(2)%alpha_rho(2)": 1000 * (1.0 - 1e-6),
82+
"patch_icpp(2)%alpha(1)": 1e-6,
83+
"patch_icpp(2)%alpha(2)": 1.0 - 1e-6,
84+
"patch_icpp(2)%tau_e(1)": 0.0,
85+
# Fluids Physical Parameters
86+
"fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
87+
"fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
88+
"fluid_pp(1)%G": 0e0,
89+
"fluid_pp(2)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
90+
"fluid_pp(2)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
91+
"fluid_pp(2)%G": 1.0e09,
92+
}
93+
)
94+
)

Diff for: ‎examples/2D_cont_damage/case.py

+104
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
#!/usr/bin/env python3
2+
import json
3+
4+
# Configuring case dictionary
5+
print(
6+
json.dumps(
7+
{
8+
# Logistics
9+
"run_time_info": "T",
10+
# Computational Domain Parameters
11+
"x_domain%beg": 0.0,
12+
"x_domain%end": 0.001,
13+
"y_domain%beg": 0.0,
14+
"y_domain%end": 0.0005,
15+
"m": 50,
16+
"n": 25,
17+
"p": 0,
18+
"dt": 2e-12,
19+
"t_step_start": 0,
20+
"t_step_stop": 40000,
21+
"t_step_save": 2000,
22+
# Simulation Algorithm Parameters
23+
"num_patches": 2,
24+
"model_eqns": 2,
25+
"alt_soundspeed": "F",
26+
"num_fluids": 2,
27+
"mpp_lim": "F",
28+
"mixture_err": "F",
29+
"time_stepper": 3,
30+
"weno_order": 5,
31+
"weno_eps": 1.0e-16,
32+
"teno": "T",
33+
"teno_CT": 1e-8,
34+
"null_weights": "F",
35+
"mp_weno": "F",
36+
"riemann_solver": 1,
37+
"wave_speeds": 1,
38+
"avg_state": 2,
39+
"bc_x%beg": -6,
40+
"bc_x%end": -6,
41+
"bc_y%beg": -2,
42+
"bc_y%end": -6,
43+
# Hypoelasticity
44+
"hypoelasticity": "T",
45+
"fd_order": 4,
46+
"cont_damage": "T",
47+
"tau_star": 0.0,
48+
"cont_damage_s": 2.0,
49+
"alpha_bar": 1e-4,
50+
# Formatted Database Files Structure Parameters
51+
"format": 1,
52+
"precision": 2,
53+
"prim_vars_wrt": "T",
54+
"parallel_io": "T",
55+
# Patch 1 Liquid
56+
"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,
61+
"patch_icpp(1)%vel(1)": 0.0,
62+
"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),
65+
"patch_icpp(1)%alpha(1)": 1.0 - 1e-6,
66+
"patch_icpp(1)%alpha_rho(2)": 1100 * 1e-6,
67+
"patch_icpp(1)%alpha(2)": 1e-6,
68+
# Patch 2 Solid
69+
"patch_icpp(2)%alter_patch(1)": "T",
70+
"patch_icpp(2)%geometry": 3,
71+
"patch_icpp(2)%x_centroid": 0.0005,
72+
"patch_icpp(2)%y_centroid": 0.000125,
73+
"patch_icpp(2)%length_x": 0.0005,
74+
"patch_icpp(2)%length_y": 0.00025,
75+
"patch_icpp(2)%vel(1)": 0.0,
76+
"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
96+
"fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
97+
"fluid_pp(1)%pi_inf": 4.4e00 * 5.57e08 / (4.4e00 - 1.0e00),
98+
"fluid_pp(1)%G": 0.0,
99+
"fluid_pp(2)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
100+
"fluid_pp(2)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
101+
"fluid_pp(2)%G": 1.0e09,
102+
}
103+
)
104+
)

Diff for: ‎src/common/m_variables_conversion.fpp

+7
Original file line numberDiff line numberDiff line change
@@ -1122,6 +1122,7 @@ contains
11221122
do i = strxb, strxe
11231123
! subtracting elastic contribution for pressure calculation
11241124
if (G_K > verysmall) then
1125+
if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp)
11251126
qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
11261127
((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
11271128
! Double for shear stresses
@@ -1149,6 +1150,8 @@ contains
11491150
qK_prim_vf(c_idx)%sf(j, k, l) = qK_cons_vf(c_idx)%sf(j, k, l)
11501151
end if
11511152

1153+
if (cont_damage) qK_prim_vf(damage_idx)%sf(j, k, l) = qK_cons_vf(damage_idx)%sf(j, k, l)
1154+
11521155
end do
11531156
end do
11541157
end do
@@ -1390,6 +1393,8 @@ contains
13901393
do i = strxb, strxe
13911394
! adding elastic contribution
13921395
if (G > verysmall) then
1396+
if (cont_damage) G = G*max((1._wp - q_prim_vf(damage_idx)%sf(j, k, l)), 0._wp)
1397+
13931398
q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
13941399
(q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
13951400
! Double for shear stresses
@@ -1413,6 +1418,8 @@ contains
14131418
q_cons_vf(c_idx)%sf(j, k, l) = q_prim_vf(c_idx)%sf(j, k, l)
14141419
end if
14151420

1421+
if (cont_damage) q_cons_vf(damage_idx)%sf(j, k, l) = q_prim_vf(damage_idx)%sf(j, k, l)
1422+
14161423
end do
14171424
end do
14181425
end do

Diff for: ‎src/post_process/m_data_output.fpp

+3
Original file line numberDiff line numberDiff line change
@@ -346,6 +346,9 @@ contains
346346
! Elastic stresses
347347
if (hypoelasticity) dbvars = dbvars + (num_dims*(num_dims + 1))/2
348348
349+
! Damage state variable
350+
if (cont_damage) dbvars = dbvars + 1
351+
349352
! Magnetic field
350353
if (mhd) then
351354
if (n == 0) then

Diff for: ‎src/post_process/m_global_parameters.fpp

+10
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,7 @@ module m_global_parameters
112112
logical :: elasticity !< elasticity modeling, true for hyper or hypo
113113
integer :: b_size !< Number of components in the b tensor
114114
integer :: tensor_size !< Number of components in the nonsymmetric tensor
115+
logical :: cont_damage !< Continuum damage modeling
115116
logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
116117
!> @}
117118

@@ -134,6 +135,7 @@ module m_global_parameters
134135
type(int_bounds_info) :: xi_idx !< Indexes of first and last reference map eqns.
135136
integer :: c_idx !< Index of color function
136137
type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns.
138+
integer :: damage_idx !< Index of damage state variable (D) for continuum damage model
137139
!> @}
138140

139141
! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
@@ -363,6 +365,7 @@ contains
363365
elasticity = .false.
364366
b_size = dflt_int
365367
tensor_size = dflt_int
368+
cont_damage = .false.
366369

367370
bc_x%beg = dflt_int; bc_x%end = dflt_int
368371
bc_y%beg = dflt_int; bc_y%end = dflt_int
@@ -715,6 +718,13 @@ contains
715718
sys_size = c_idx
716719
end if
717720

721+
if (cont_damage) then
722+
damage_idx = sys_size + 1
723+
sys_size = damage_idx
724+
else
725+
damage_idx = dflt_int
726+
end if
727+
718728
end if
719729

720730
if (chemistry) then

Diff for: ‎src/post_process/m_mpi_proxy.fpp

+2-1
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,8 @@ contains
170170
& 'file_per_process', 'relax', 'cf_wrt', &
171171
& 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', &
172172
& 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', &
173-
& 'rkck_adap_dt', 'output_partial_domain', 'relativity']
173+
& 'rkck_adap_dt', 'output_partial_domain', 'relativity', &
174+
& 'cont_damage' ]
174175
call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
175176
#:endfor
176177

Diff for: ‎src/post_process/m_start_up.f90

+9-1
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,7 @@ subroutine s_read_input_file
8484
relax_model, cf_wrt, sigma, adv_n, ib, num_ibs, &
8585
cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
8686
cfl_target, surface_tension, bubbles_lagrange, rkck_adap_dt, &
87-
sim_data, hyperelasticity, Bx0, relativity
87+
sim_data, hyperelasticity, Bx0, relativity, cont_damage
8888

8989
! Inquiring the status of the post_process.inp file
9090
file_loc = 'post_process.inp'
@@ -418,6 +418,14 @@ subroutine s_save_data(t_step, varname, pres, c, H)
418418
end do
419419
end if
420420

421+
if (cont_damage) then
422+
q_sf = q_cons_vf(damage_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
423+
write (varname, '(A)') 'damage_state'
424+
call s_write_variable_to_formatted_database_file(varname, t_step)
425+
426+
varname(:) = ' '
427+
end if
428+
421429
! Adding the pressure to the formatted database file
422430
if (pres_wrt .or. prim_vars_wrt) then
423431
q_sf = q_prim_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)

Diff for: ‎src/pre_process/m_data_output.fpp

+2
Original file line numberDiff line numberDiff line change
@@ -341,6 +341,8 @@ contains
341341
write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)/nbub
342342
else if (i == n_idx .and. adv_n .and. bubbles_euler) then
343343
write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
344+
else if (i == damage_idx) then
345+
write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
344346
end if
345347
end do
346348
close (2)

Diff for: ‎src/pre_process/m_global_parameters.fpp

+8
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,7 @@ module m_global_parameters
9393
integer :: b_size !< Number of components in the b tensor
9494
integer :: tensor_size !< Number of components in the nonsymmetric tensor
9595
logical :: pre_stress !< activate pre_stressed domain
96+
logical :: cont_damage !< continuum damage modeling
9697
logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
9798

9899
! Annotations of the structure, i.e. the organization, of the state vectors
@@ -111,6 +112,7 @@ module m_global_parameters
111112
type(int_bounds_info) :: xi_idx !< Indexes of first and last reference map eqns.
112113
integer :: c_idx !< Index of the color function
113114
type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns.
115+
integer :: damage_idx !< Index of damage state variable (D) for continuum damage model
114116

115117
! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
116118
! Stands for "InDices With BUFFer".
@@ -340,6 +342,7 @@ contains
340342
pre_stress = .false.
341343
b_size = dflt_int
342344
tensor_size = dflt_int
345+
cont_damage = .false.
343346

344347
mhd = .false.
345348
relativity = .false.
@@ -808,6 +811,11 @@ contains
808811
sys_size = c_idx
809812
end if
810813

814+
if (cont_damage) then
815+
damage_idx = sys_size + 1
816+
sys_size = damage_idx
817+
end if
818+
811819
end if
812820

813821
if (chemistry) then

0 commit comments

Comments
 (0)