Implemented Blanco-Martin model for rock salt #24

.github/workflows/CI.yml

@@ -34,7 +34,7 @@ jobs:
     - name: Install MOOSE packages
       shell: bash -l {0}
       run: |
-        mamba install moose-dev=2024.10.17=mpich
+        mamba install moose-dev=2024.12.23=mpich
     - name: Compile BEAVER
       shell: bash -l {0}
       run: |
@@ -63,7 +63,7 @@ jobs:
     - name: Install MOOSE packages
       shell: bash -l {0}
       run: |
-        mamba install moose-dev=2024.10.17=mpich lcov
+        mamba install moose-dev=2024.12.23=mpich lcov
     - name: Compile and test BEAVER (with coverage)
       shell: bash -l {0}
       run: |

examples/viscoelasticity/blanco-martin/RTL_S1-1.py

@@ -0,0 +1,31 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+plt.style.use('publication.mplstyle')
+plt.rcParams['text.usetex'] = False  # Ensure LaTeX is not used for rendering text
+
+def read_csv(filename):
+    # Read data
+    time, EpsAx = np.loadtxt(filename, delimiter=',', skiprows=1, usecols=[0, 4], unpack=True)
+    return time, EpsAx
+
+def plot_figure():
+    fig, ax = plt.subplots()
+    ax.set_title('Axial_strain (-) vs time (d)')  # Use Unicode directly
+
+# Reading data
+    time, EpsAx = read_csv("RTL_S1-1.csv")   
+# Plot AE data
+    ax.plot(time, EpsAx, color="k", label="EpsAx")  
+# Additional plot settings
+    ax.set_xlabel('Time (days)')
+    ax.set_ylabel('EpsAx (-)')
+  #  ax.set_yscale('log')  # Set y-axis to log base 10
+    ax.legend()  
+# Axes limit
+#    ax.set_xlim(4, 8)
+# Save figure
+    fig.savefig("RTL_S1-1.pdf", format="pdf", dpi=300, bbox_inches="tight")
+
+if __name__ == '__main__':
+    plot_figure()

examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.i

@@ -0,0 +1,296 @@
+# Modified Lemaitre creep model
+# See Blanco-Martin et al. (2023)
+# Parameters
+# Units: stress in MPa, time in days, strain in m / m
+E = 12000
+nu = 0.3
+alpha = 0.575
+kr1 = 1.302
+beta1 = 3.053
+kr2 = 0.091
+beta2 = 1.053
+A = 100
+n = 9
+A1 = 0.034
+n1 = 1.499
+P = 5.0
+Q = 5.0
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  # nx = 5
+  # ny = 10
+  # nz = 5
+  nx = 1
+  ny = 1
+  nz = 1
+  xmin = 0
+  xmax = 65e-03
+  ymin = 0
+  ymax = 130e-03
+  zmin = 0
+  zmax = 65e-03
+[]
+
+[Variables]
+  [disp_x]
+    order = FIRST
+    family = LAGRANGE
+  []
+  [disp_y]
+    order = FIRST
+    family = LAGRANGE
+  []
+  [disp_z]
+    order = FIRST
+    family = LAGRANGE
+  []
+[]
+
+[Kernels]
+  [stress_x]
+    type = BVStressDivergence
+    component = x
+    variable = disp_x
+  []
+  [stress_y]
+    type = BVStressDivergence
+    component = y
+    variable = disp_y
+  []
+  [stress_z]
+    type = BVStressDivergence
+    component = z
+    variable = disp_z
+  []
+[]
+
+[AuxVariables]
+  [eqv_stress]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [eqv_strain]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [eqv_strain_rate]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [eqv_creep_strain_L]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+  [stress_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+[]
+
+[AuxKernels]
+  [eqv_stress_aux]
+    type = BVMisesStressAux
+    variable = eqv_stress
+    execute_on = 'TIMESTEP_END'
+  []
+  [eqv_strain_aux]
+    type = BVEqvStrainAux
+    variable = eqv_strain
+    execute_on = 'TIMESTEP_END'
+  []
+  [eqv_strain_rate_aux]
+    type = BVEqvStrainRateAux
+    variable = eqv_strain_rate
+    execute_on = 'TIMESTEP_END'
+  []
+  [eqv_creep_strain_L_aux]
+    type = ADMaterialRealAux
+    variable = eqv_creep_strain_L
+    property = eqv_creep_strain_L
+    execute_on = 'TIMESTEP_END'
+  []
+  [strain_zz_aux]
+    type = BVStrainComponentAux
+    variable = strain_yy
+    index_i = y
+    index_j = y
+    execute_on = 'TIMESTEP_END'
+  []
+  [stress_zz_aux]
+    type = BVStressComponentAux
+    variable = stress_yy
+    index_i = y
+    index_j = y
+    execute_on = 'TIMESTEP_END'
+  []
+[]
+
+[Functions]
+  [loading2]
+    type = ParsedFunction 
+    expression = 'if(t<=20,5.5,if(t<=55,6,if(t<=76,7,if(t<=97,10,if(t<=118,15,if(t<=139,20,25))))))' 
+  []
+[]
+
+[BCs]
+  [no_x]
+    type = DirichletBC
+    variable = disp_x
+    boundary = 'left'
+    value = 0.0
+  []
+  [no_y]
+    type = DirichletBC
+    variable = disp_y
+    boundary = 'bottom'
+    value = 0.0
+  []
+  [no_z]
+    type = DirichletBC
+    variable = disp_z
+    boundary = 'back'
+    value = 0.0
+  []
+  [BVPressure]
+    [pressure_right]
+      boundary = 'right'
+      displacement_vars = 'disp_x disp_y disp_z'
+      value = ${P}
+    []
+    [pressure_front]
+      boundary = 'front'
+      displacement_vars = 'disp_x disp_y disp_z'
+      value = ${P}
+    []
+    [pressure_top]
+      boundary = 'top'
+      displacement_vars = 'disp_x disp_y disp_z'
+      function = loading2
+    []
+  []
+[]
+
+[Materials]
+  [elasticity]
+    type = BVMechanicalMaterial
+    displacements = 'disp_x disp_y disp_z'
+    young_modulus = ${E}
+    poisson_ratio = ${nu}
+    initial_stress = '-${P} -${Q} -${P}'
+    inelastic_models = 'viscoelastic'
+  []
+  [viscoelastic]
+    type = BVBlancoMartinModelUpdate 
+    alpha = ${alpha}
+    kr1 = ${kr1}
+    beta1 = ${beta1} 
+    kr2 = ${kr2}
+    beta2 = ${beta2}
+    A1 = ${A1}
+    n1 = ${n1}
+    A = ${A}
+    n = ${n}
+    B = 0.0
+    m = ${n}
+  []
+[]
+
+[Postprocessors]
+   [eqv_strain_rate]
+    type = ElementAverageValue
+    variable = eqv_strain_rate
+    outputs = csv
+  []
+  [eqv_strain]
+    type = ElementAverageValue
+    variable = eqv_strain
+    outputs = csv
+  []
+  [eqv_strain_L]
+    type = ElementAverageValue
+    variable = eqv_creep_strain_L
+    outputs = csv
+  []
+  [strain_zz]
+    type = ElementAverageValue
+    variable = strain_yy
+    outputs = csv
+  []
+  [stress_zz]
+    type = ElementAverageValue
+    variable = stress_yy
+    outputs = csv
+  []
+  [q]
+    type = ElementAverageValue
+    variable = eqv_stress
+    outputs = csv
+  []
+[]
+
+[Preconditioning]
+  active = 'hypre'
+  [hypre]
+    type = SMP
+    full = true
+    petsc_options = '-snes_ksp_ew'
+    petsc_options_iname = '-pc_type -pc_hypre_type
+                           -snes_atol -snes_rtol -snes_stol -snes_max_it
+                           -snes_linesearch_type'
+    petsc_options_value = 'hypre boomeramg
+                           1.0e-10 1.0e-12 0 20
+                           basic'
+  []
+  [superlu]
+    type = SMP
+    full = true
+    petsc_options = '-snes_ksp_ew -snes_converged_reason -ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_diagonal_scale -ksp_diagonal_scale_fix'
+    petsc_options_iname = '-snes_type
+                           -snes_atol -snes_rtol -snes_max_it
+                           -pc_type -pc_factor_mat_solver_package
+                           -snes_linesearch_type'
+    petsc_options_value = 'newtonls
+                           1e-10 1e-12 50
+                           lu superlu_dist
+                           l2'
+  []
+  [asm]
+    type = SMP
+    petsc_options = '-snes_ksp_ew'
+    petsc_options_iname = '-ksp_type
+                           -pc_type
+                           -sub_pc_type
+                           -snes_type -snes_atol -snes_rtol -snes_max_it -snes_linesearch_type
+                           -ksp_gmres_restart'
+    petsc_options_value = 'fgmres
+                           asm
+                           ilu
+                           newtonls 1e-10 1e-10 120 basic
+                           201'
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'NEWTON'
+  automatic_scaling = true
+  start_time = 0.0
+  end_time = 200 # 200 days
+  dt = 0.02
+  timestep_tolerance = 1.0e-10
+[]
+
+[Outputs]
+  perf_graph = true
+  exodus = true
+  [csv]
+  type = CSV
+  execute_on = 'timestep_end'
+  []
+[]

examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.py

@@ -0,0 +1,41 @@
+import os, csv
+import numpy as np
+import matplotlib.pyplot as plt
+plt.style.use('../../publication.mplstyle')
+
+def plot_experimental_data(ax):
+  # Read data from file
+  t, stress, strain = np.loadtxt("S1-1.txt", skiprows=7, usecols=[0, 1, 3], unpack=True)
+
+  # ax.plot(t, stress, color="k")
+  ax[0].plot(t, strain / 1.0e+04, color="k")
+  ax[1].plot(t, stress, color="k", label=r"Experimental data")
+
+def plot_simulation_data(ax):
+  # Read data from file
+  t, stress, strain = np.loadtxt("blanco-martin-RTL-2024_csv.csv", delimiter=',', skiprows=1, usecols=[0, 6, 5], unpack=True)
+
+  # ax.plot(t, -stress, color="xkcd:red")
+  ax[0].plot(t, -strain * 1.0e+02, color="xkcd:red")
+  ax[1].plot(t, -stress, color="xkcd:red", label=r"Numerical fit")
+
+if __name__ == "__main__":
+
+  # Figure stress
+  fig, axes = plt.subplots(2, 1, figsize=(3.25, 5), sharex=True)
+  # plt.subplots_adjust(hspace=0.12)
+
+  plot_experimental_data(axes)
+  plot_simulation_data(axes)
+
+  axes[1].set_xlim(0, 200)
+  axes[0].set_ylim(0, 8)
+  axes[1].set_ylim(0, 30)
+  axes[1].set_xlabel(r"Time, (days)")
+  axes[0].set_ylabel(r"Axial strain, $\varepsilon_{\text{ax}}$ (\%)")
+  axes[1].set_ylabel(r"Axial stress, $\sigma_{\text{ax}}$ (MPa)")
+  axes[0].set_title(r"Modified RTL2020 model (Blanco-Martìn et al., 2024)")
+  axes[1].legend(loc="lower right")
+
+  # plt.show()
+  fig.savefig("blanco-martin-RTL-2024.png", format="PNG", dpi=300, bbox_inches="tight")