Continuum Damage Model (#816)

Author: ChrisZYJ

docs/documentation/case.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>

docs/documentation/references.md

@@ -1,7 +1,7 @@
 # Flow visualization
 
 A post-processed database in Silo-HDF5 format can be visualized and analyzed using Paraview and VisIt.
-After the post-processing of simulation data (see section [Running](running.md#running-1)), a directory named `silo_hdf5` contains a silo-HDF5 database.
+After the post-processing of simulation data (see section [Running](running.md)), a directory named `silo_hdf5` contains a silo-HDF5 database.
 Here, `silo_hdf5/` includes a directory named `root/` that contains index files for flow field data at each saved time step.
 
 ### Visualizing with Paraview
@@ -38,7 +38,7 @@ For many cases, this step will require resizing the render view window.
 VisIt is an alternative open-source interactive parallel visualization and graphical analysis tool for viewing scientific data.
 Versions of VisIt after 2.6.0 have been confirmed to work with the MFC databases for some parallel environments.
 Nevertheless, installation and configuration of VisIt can be environment-dependent and are left to the user.
-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).
+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).
 
 The user can launch VisIt and open the index files under `/silo_hdf5/root`.
 Once the database is loaded, flow field variables contained in the database can be added to the plot.

docs/documentation/visualization.md

@@ -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/1D_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,
+        }
+    )
+)

examples/2D_cont_damage/case.py

@@ -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/common/m_variables_conversion.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_data_output.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_global_parameters.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_mpi_proxy.fpp

@@ -84,7 +84,7 @@ subroutine s_read_input_file
             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 @@ subroutine s_save_data(t_step, varname, pres, c, H)
             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)