Continuum Damage Model

docs/documentation/case.md

@@ -886,7 +886,20 @@ Note: In 1D/2D/3D simulations, all three velocity components are treated as stat
 
 Note: For relativistic flow, the conservative and primitive densities are different. `rho_wrt` outputs the primitive (rest mass) density.
 
-### 15. Cylindrical Coordinates
+### 15. Elasticity
+
+| Parameter         | Type    | Description                                         |
+| ---:              | :---:   | :---                                                |
+| `hypoelasticity`  | Logical | Enable hypoelasticity simulation                    |
+| `hyperelasticity` | Logical | Enable hyperelasticity simulation                   |
+| `cont_damage`     | Logical | Enable continuum damage model                       |
+| `tau_star`        | Real    | Threshold stress for continuum damage model         |
+| `cont_damage_s`   | Real    | Power `s` for continuum damage model                |
+| `alpha_bar`       | Real    | Damage factor (rate) for continuum damage model     |
+
+- `cont_damage` activates continuum damage model for solid materials. Requires `tau_star`, `cont_damage_s`, and `alpha_bar` to be set (empirically determined) ([Cao et al., 2019](references.md#cao19)).
+
+### 16. Cylindrical Coordinates
 
 When ``cyl_coord = 'T'`` is set in 3D the following constraints must be met:
 

docs/documentation/references.md

@@ -60,4 +60,6 @@
 
 - <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>
 
-- <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>
+- <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>
+
+- <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>

examples/1D_cont_damage/case.py

@@ -0,0 +1,94 @@
+#!/usr/bin/python
+import json
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0e00,
+            "x_domain%end": 1.0e00,
+            "m": 100,
+            "n": 0,
+            "p": 0,
+            "dt": 5e-10,
+            "t_step_start": 0,
+            "t_step_stop": 50000,
+            "t_step_save": 50000,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "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": 1,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            # Turning on Hypoelasticity
+            "hypoelasticity": "T",
+            "fd_order": 4,
+            "cont_damage": "T",
+            "tau_star": 0.0,
+            "cont_damage_s": 2.0,
+            "alpha_bar": 1e-4,
+            # Acoustic Source
+            "acoustic_source": "T",
+            "num_source": 1,
+            "acoustic(1)%support": 1,
+            "acoustic(1)%loc(1)": 0.1,
+            "acoustic(1)%pulse": 1,
+            "acoustic(1)%npulse": 999,
+            "acoustic(1)%dir": 1.0,
+            "acoustic(1)%mag": 1000.0,
+            "acoustic(1)%frequency": 1e4,
+            "acoustic(1)%delay": 0,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            # Patch 1 L
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.25,
+            "patch_icpp(1)%length_x": 0.5,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%pres": 1e5,
+            "patch_icpp(1)%alpha_rho(1)": 1000 * (1.0 - 1e-6),
+            "patch_icpp(1)%alpha_rho(2)": 1000 * 1e-6,
+            "patch_icpp(1)%alpha(1)": 1.0 - 1e-6,
+            "patch_icpp(1)%alpha(2)": 1e-6,
+            "patch_icpp(1)%tau_e(1)": 0.0,
+            # Patch 2 R
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%x_centroid": 0.75,
+            "patch_icpp(2)%length_x": 0.5,
+            "patch_icpp(2)%vel(1)": 0.0,
+            "patch_icpp(2)%pres": 1e5,
+            "patch_icpp(2)%alpha_rho(1)": 1000 * 1e-6,
+            "patch_icpp(2)%alpha_rho(2)": 1000 * (1.0 - 1e-6),
+            "patch_icpp(2)%alpha(1)": 1e-6,
+            "patch_icpp(2)%alpha(2)": 1.0 - 1e-6,
+            "patch_icpp(2)%tau_e(1)": 0.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
+            "fluid_pp(1)%G": 0e0,
+            "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,
+        }
+    )
+)

examples/2D_cont_damage/case.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python3
+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": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "teno": "T",
+            "teno_CT": 1e-8,
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 1,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "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": 1,
+            "precision": 2,
+            "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,
+        }
+    )
+)

src/common/m_variables_conversion.fpp

@@ -1122,6 +1122,7 @@ contains
                         do i = strxb, strxe
                             ! subtracting elastic contribution for pressure calculation
                             if (G_K > verysmall) then
+                                if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp)
                                 qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
                                                                 ((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
                                 ! Double for shear stresses
@@ -1149,6 +1150,8 @@ contains
                         qK_prim_vf(c_idx)%sf(j, k, l) = qK_cons_vf(c_idx)%sf(j, k, l)
                     end if
 
+                    if (cont_damage) qK_prim_vf(damage_idx)%sf(j, k, l) = qK_cons_vf(damage_idx)%sf(j, k, l)
+
                 end do
             end do
         end do
@@ -1390,6 +1393,8 @@ contains
                         do i = strxb, strxe
                             ! adding elastic contribution
                             if (G > verysmall) then
+                                if (cont_damage) G = G*max((1._wp - q_prim_vf(damage_idx)%sf(j, k, l)), 0._wp)
+
                                 q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
                                                                (q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
                                 ! Double for shear stresses
@@ -1413,6 +1418,8 @@ contains
                         q_cons_vf(c_idx)%sf(j, k, l) = q_prim_vf(c_idx)%sf(j, k, l)
                     end if
 
+                    if (cont_damage) q_cons_vf(damage_idx)%sf(j, k, l) = q_prim_vf(damage_idx)%sf(j, k, l)
+
                 end do
             end do
         end do

src/post_process/m_data_output.fpp

@@ -346,6 +346,9 @@ contains
             ! Elastic stresses
             if (hypoelasticity) dbvars = dbvars + (num_dims*(num_dims + 1))/2
 
+            ! Damage state variable
+            if (cont_damage) dbvars = dbvars + 1
+
             ! Magnetic field
             if (mhd) then
                 if (n == 0) then

src/post_process/m_global_parameters.fpp

@@ -112,6 +112,7 @@ module m_global_parameters
     logical :: elasticity      !< elasticity modeling, true for hyper or hypo
     integer :: b_size          !< Number of components in the b tensor
     integer :: tensor_size     !< Number of components in the nonsymmetric tensor
+    logical :: cont_damage     !< Continuum damage modeling
     logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
     !> @}
 
@@ -134,6 +135,7 @@ module m_global_parameters
     type(int_bounds_info) :: xi_idx                !< Indexes of first and last reference map eqns.
     integer :: c_idx                               !< Index of color function
     type(int_bounds_info) :: species_idx           !< Indexes of first & last concentration eqns.
+    integer :: damage_idx                          !< Index of damage state variable (D) for continuum damage model
     !> @}
 
     ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
@@ -363,6 +365,7 @@ contains
         elasticity = .false.
         b_size = dflt_int
         tensor_size = dflt_int
+        cont_damage = .false.
 
         bc_x%beg = dflt_int; bc_x%end = dflt_int
         bc_y%beg = dflt_int; bc_y%end = dflt_int
@@ -715,6 +718,13 @@ contains
                 sys_size = c_idx
             end if
 
+            if (cont_damage) then
+                damage_idx = sys_size + 1
+                sys_size = damage_idx
+            else
+                damage_idx = dflt_int
+            end if
+
         end if
 
         if (chemistry) then

src/post_process/m_mpi_proxy.fpp

@@ -170,7 +170,8 @@ contains
             & 'file_per_process', 'relax', 'cf_wrt',                           &
             & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt',          &
             & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange',        &
-            & 'rkck_adap_dt', 'output_partial_domain', 'relativity']
+            & 'rkck_adap_dt', 'output_partial_domain', 'relativity',           &
+            & 'cont_damage' ]
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 

src/post_process/m_start_up.f90

@@ -84,7 +84,7 @@
             relax_model, cf_wrt, sigma, adv_n, ib, num_ibs, &
             cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
             cfl_target, surface_tension, bubbles_lagrange, rkck_adap_dt, &
-            sim_data, hyperelasticity, Bx0, relativity
+            sim_data, hyperelasticity, Bx0, relativity, cont_damage
 
         ! Inquiring the status of the post_process.inp file
         file_loc = 'post_process.inp'
@@ -418,6 +418,14 @@
             end do
         end if
 
+        if (cont_damage) then
+            q_sf = q_cons_vf(damage_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            write (varname, '(A)') 'damage_state'
+            call s_write_variable_to_formatted_database_file(varname, t_step)
+
+            varname(:) = ' '
+        end if
+
         ! Adding the pressure to the formatted database file
         if (pres_wrt .or. prim_vars_wrt) then
             q_sf = q_prim_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)

src/pre_process/m_data_output.fpp

@@ -341,6 +341,8 @@
                             write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)/nbub
                         else if (i == n_idx .and. adv_n .and. bubbles_euler) then
                             write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
+                        else if (i == damage_idx) then
+                            write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
                         end if
                     end do
                     close (2)

src/pre_process/m_global_parameters.fpp

@@ -93,6 +93,7 @@ module m_global_parameters
     integer :: b_size                !< Number of components in the b tensor
     integer :: tensor_size           !< Number of components in the nonsymmetric tensor
     logical :: pre_stress            !< activate pre_stressed domain
+    logical :: cont_damage           !< continuum damage modeling
     logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
 
     ! Annotations of the structure, i.e. the organization, of the state vectors
@@ -111,6 +112,7 @@ module m_global_parameters
     type(int_bounds_info) :: xi_idx                !< Indexes of first and last reference map eqns.
     integer :: c_idx                               !< Index of the color function
     type(int_bounds_info) :: species_idx           !< Indexes of first & last concentration eqns.
+    integer :: damage_idx                          !< Index of damage state variable (D) for continuum damage model
 
     ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
     ! Stands for "InDices With BUFFer".
@@ -340,6 +342,7 @@ contains
         pre_stress = .false.
         b_size = dflt_int
         tensor_size = dflt_int
+        cont_damage = .false.
 
         mhd = .false.
         relativity = .false.
@@ -808,6 +811,11 @@ contains
                 sys_size = c_idx
             end if
 
+            if (cont_damage) then
+                damage_idx = sys_size + 1
+                sys_size = damage_idx
+            end if
+
         end if
 
         if (chemistry) then