diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e8f1e0b..5fcc64a 100644 --- a/.github/workflows/CI.yml +++ b/.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: | diff --git a/examples/viscoelasticity/blanco-martin/RTL_S1-1.py b/examples/viscoelasticity/blanco-martin/RTL_S1-1.py new file mode 100644 index 0000000..c374f24 --- /dev/null +++ b/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() diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.i b/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.i new file mode 100644 index 0000000..501cece --- /dev/null +++ b/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' + [] +[] \ No newline at end of file diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.py b/examples/viscoelasticity/blanco-martin/blanco-martin-RTL-2024.py new file mode 100644 index 0000000..bf4569d --- /dev/null +++ b/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") diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.i b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.i new file mode 100644 index 0000000..ef9cefa --- /dev/null +++ b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.i @@ -0,0 +1,300 @@ +# Modified Lemaitre creep model +# See Blanco-Martin et al. (2024) +# Parameters +# Units: stress in MPa, time in days, strain in m / m +E = 12000 +nu = 0.3 +alpha = 0.326 +kr1 = 0.7 +beta1 = 2.922 +kr2 = 0.009 +beta2 = 0.867 +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 + 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} + # [] + [viscoelastic] + type = BVModifiedLemaitreModelUpdate + alpha = ${alpha} + kr1 = ${kr1} + beta1 = ${beta1} + kr2 = ${kr2} + beta2 = ${beta2} + [] +[] + +[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' + [] +[] \ No newline at end of file diff --git a/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.py b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-2024.py new file mode 100644 index 0000000..8d447f8 --- /dev/null +++ b/examples/viscoelasticity/blanco-martin/blanco-martin-lemaitre-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-lemaitre-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 Lemaitre model (Blanco-Martìn et al., 2024)") + axes[1].legend(loc="lower right") + + # plt.show() + fig.savefig("blanco-martin-lemaitre-2024.png", format="PNG", dpi=300, bbox_inches="tight") diff --git a/include/materials/BVBlancoMartinModelUpdate.h b/include/materials/BVBlancoMartinModelUpdate.h new file mode 100644 index 0000000..6d5b8f4 --- /dev/null +++ b/include/materials/BVBlancoMartinModelUpdate.h @@ -0,0 +1,64 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#pragma once + +#include "BVTwoCreepUpdateBase.h" + +class BVBlancoMartinModelUpdate : public BVTwoCreepUpdateBase +{ +public: + static InputParameters validParams(); + BVBlancoMartinModelUpdate(const InputParameters & parameters); + +protected: + virtual void initQpStatefulProperties() override; + virtual ADReal creepRate(const std::vector & eqv_strain_incr, + const unsigned int i) override; + virtual ADReal creepRateR(const std::vector & eqv_strain_incr); + virtual ADReal creepRateLemaitre(const std::vector & eqv_strain_incr); + virtual ADReal creepRateMunsonDawson(const std::vector & eqv_strain_incr); + virtual ADReal creepRateDerivative(const std::vector & eqv_strain_incr, + const unsigned int i, + const unsigned int j) override; + virtual ADReal creepRateRDerivative(const std::vector & eqv_strain_incr); + virtual ADReal creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, + const unsigned int /*j*/); + virtual ADReal creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, + const unsigned int j); + virtual ADReal lemaitreCreepStrain(const std::vector & eqv_strain_incr); + virtual ADReal munsondawsonCreepStrain(const std::vector & eqv_strain_incr); + virtual void preReturnMap() override; + virtual void postReturnMap(const std::vector & eqv_strain_incr) override; + + // Lemaitre creep strain rate parameters + const Real _alpha; + const Real _kr1; + const Real _kr2; + const Real _beta1; + const Real _beta2; + + // Munson-Dawson creep strain rate parameters + const Real _A1; + const Real _n1; + const Real _A; + const Real _B; + const Real _m; + const Real _n; + + // Internal variable for Lemaitre and Munson-Dawson creep strain + ADMaterialProperty & _eqv_creep_strain_L; + const MaterialProperty & _eqv_creep_strain_L_old; + ADMaterialProperty & _eqv_creep_strain_R; + const MaterialProperty & _eqv_creep_strain_R_old; +}; diff --git a/include/materials/BVModifiedLemaitreModelUpdate.h b/include/materials/BVModifiedLemaitreModelUpdate.h new file mode 100644 index 0000000..148d340 --- /dev/null +++ b/include/materials/BVModifiedLemaitreModelUpdate.h @@ -0,0 +1,44 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#pragma once + +#include "BVCreepUpdateBase.h" + +class BVModifiedLemaitreModelUpdate : public BVCreepUpdateBase +{ +public: + static InputParameters validParams(); + BVModifiedLemaitreModelUpdate(const InputParameters & parameters); + +protected: + virtual void initQpStatefulProperties() override; + virtual ADReal creepRateR(const ADReal & eqv_strain_incr); + virtual ADReal creepRate(const ADReal & eqv_strain_incr) override; + virtual ADReal creepRateRDerivative(const ADReal & eqv_strain_incr); + virtual ADReal creepRateDerivative(const ADReal & eqv_strain_incr) override; + virtual ADReal lemaitreCreepStrain(const ADReal & eqv_strain_incr); + virtual void preReturnMap() override; + virtual void postReturnMap(const ADReal & eqv_strain_incr) override; + + // Lemaitre creep strain rate parameters + const Real _alpha; + const Real _kr1; + const Real _kr2; + const Real _beta1; + const Real _beta2; + + // Internal variable for creep strain + ADMaterialProperty & _eqv_creep_strain; + const MaterialProperty & _eqv_creep_strain_old; +}; diff --git a/include/materials/BVRTL2020ModelUpdate.h b/include/materials/BVRTL2020ModelUpdate.h index fba7bbb..93f134d 100644 --- a/include/materials/BVRTL2020ModelUpdate.h +++ b/include/materials/BVRTL2020ModelUpdate.h @@ -25,11 +25,13 @@ class BVRTL2020ModelUpdate : public BVTwoCreepUpdateBase virtual void initQpStatefulProperties() override; virtual ADReal creepRate(const std::vector & eqv_strain_incr, const unsigned int i) override; + virtual ADReal creepRateR(const std::vector & eqv_strain_incr); virtual ADReal creepRateLemaitre(const std::vector & eqv_strain_incr); virtual ADReal creepRateMunsonDawson(const std::vector & eqv_strain_incr); virtual ADReal creepRateDerivative(const std::vector & eqv_strain_incr, const unsigned int i, const unsigned int j) override; + virtual ADReal creepRateRDerivative(const std::vector & eqv_strain_incr); virtual ADReal creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, const unsigned int /*j*/); virtual ADReal creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, diff --git a/src/materials/BVBlancoMartinModelUpdate.C b/src/materials/BVBlancoMartinModelUpdate.C new file mode 100644 index 0000000..6fe41e7 --- /dev/null +++ b/src/materials/BVBlancoMartinModelUpdate.C @@ -0,0 +1,229 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#include "BVBlancoMartinModelUpdate.h" + +registerMooseObject("BeaverApp", BVBlancoMartinModelUpdate); + +InputParameters +BVBlancoMartinModelUpdate::validParams() +{ + InputParameters params = BVTwoCreepUpdateBase::validParams(); + params.addClassDescription( + "Material for computing a RTL2020 creep update. See Azabou et al. (2021), Rock salt " + "behavior: From laboratory experiments to pertinent long-term predictions."); + // Modified Lemaitre creep strain rate parameters + params.addRequiredRangeCheckedParam("alpha", "0.0 < alpha & alpha < 1.0", "The alpha parameter."); + params.addRequiredRangeCheckedParam("kr1", "kr1 > 0.0", "The kr1 parameter."); + params.addRequiredRangeCheckedParam("kr2", "kr2 > 0.0", "The kr2 parameter."); + params.addRequiredRangeCheckedParam("beta1", "beta1 > 0.0", "The beta1 parameter."); + params.addRequiredRangeCheckedParam("beta2", "beta2 > 0.0", "The beta2 parameter."); + // Munson-Dawson creep strain rate parameters + params.addRequiredRangeCheckedParam("A1", "A1 > 0.0", "The A1 parameter."); + params.addRequiredRangeCheckedParam("n1", "n1 > 0.0", "The n1 parameter."); + params.addRequiredRangeCheckedParam("A", "A >= 0.0", "The A parameter."); + params.addRequiredRangeCheckedParam("B", "B >= 0.0", "The B parameter."); + params.addRequiredRangeCheckedParam("m", "m > 1.0", "The m parameter."); + params.addRequiredRangeCheckedParam("n", "n > 1.0", "The n parameter."); + return params; +} + +BVBlancoMartinModelUpdate::BVBlancoMartinModelUpdate(const InputParameters & parameters) + : BVTwoCreepUpdateBase(parameters), + // Modified Lemaitre creep strain rate parameters + _alpha(getParam("alpha")), + _kr1(getParam("kr1")), + _kr2(getParam("kr2")), + _beta1(getParam("beta1")), + _beta2(getParam("beta2")), + // Munson-Dawson creep strain rate parameters + _A1(getParam("A1")), + _n1(getParam("n1")), + _A(getParam("A")), + _B(getParam("B")), + _m(getParam("m")), + _n(getParam("n")), + // Internal variable for Lemaitre and Munson-Dawson creep strain + _eqv_creep_strain_L(declareADProperty(_base_name + "eqv_creep_strain_L")), + _eqv_creep_strain_L_old(getMaterialPropertyOld(_base_name + "eqv_creep_strain_L")), + _eqv_creep_strain_R(declareADProperty(_base_name + "eqv_creep_strain_R")), + _eqv_creep_strain_R_old(getMaterialPropertyOld(_base_name + "eqv_creep_strain_R")) +{ +} + +void +BVBlancoMartinModelUpdate::initQpStatefulProperties() +{ + _eqv_creep_strain_L[_qp] = 0.0; + _eqv_creep_strain_R[_qp] = 0.0; +} + +ADReal +BVBlancoMartinModelUpdate::creepRate(const std::vector & eqv_strain_incr, const unsigned int i) +{ + if (i == 0) // Lemaitre + return creepRateLemaitre(eqv_strain_incr); + else if (i == 1) // Munson-Dawson + return creepRateMunsonDawson(eqv_strain_incr); + else + throw MooseException("BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRate`!"); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateR(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 0.0; + else + return 1.0e-06 * std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateLemaitre(const std::vector & eqv_strain_incr) +{ + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateR(eqv_strain_incr); + else + return _alpha * creepRateR(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateMunsonDawson(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; + + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); + + if (gamma_ms < saturation_strain) + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + creepRateR(eqv_strain_incr); + else + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateR(eqv_strain_incr); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateDerivative(const std::vector & eqv_strain_incr, + const unsigned int i, + const unsigned int j) +{ + if (i == 0) // Lemaitre + return creepRateLemaitreDerivative(eqv_strain_incr, j); + else if (i == 1) // Munson-Dawson + return creepRateMunsonDawsonDerivative(eqv_strain_incr, j); + else + throw MooseException( + "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateRDerivative(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 1.0; + else + return -1.0e-06 * 3.0 * _G / _alpha * + std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha - 1.0) * + (_beta1 / _kr1 * std::pow(q / _kr1, _beta1 - 1.0) + + _beta2 / _kr2 * std::pow(q / _kr2, _beta2 - 1.0)); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, + const unsigned int j) +{ + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (j == 0) // Lemaitre wrt Lemaitre + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); + else + return std::pow(gamma_l, -1.0 / _alpha) * + (_alpha * gamma_l * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * (_alpha - 1.0) * creepRateR(eqv_strain_incr)); + + else if (j == 1) // Lemaitre wrt Munson-Dawson + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); + else + return _alpha * creepRateRDerivative(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); + + else + throw MooseException( + "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); +} + +ADReal +BVBlancoMartinModelUpdate::creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, + const unsigned int j) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; + + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); + + if (j == 0) // Munson-Dawson wrt Lemaitre + if (gamma_ms < saturation_strain) + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + creepRateRDerivative(eqv_strain_incr); + else + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateRDerivative(eqv_strain_incr); + + else if (j == 1) // Munson-Dawson wrt Munson-Dawson + if (gamma_ms < saturation_strain) + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * + ((1.0 - gamma_ms / saturation_strain) * creepRateRDerivative(eqv_strain_incr) - + 1.0e+06 * _n / saturation_strain * creepRateR(eqv_strain_incr)); + else + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m - 1.0) * + ((gamma_ms / saturation_strain - 1.0) * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * _m / saturation_strain * creepRateR(eqv_strain_incr)); + + else + throw MooseException( + "BVBlancoMartinModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); +} + +ADReal +BVBlancoMartinModelUpdate::lemaitreCreepStrain(const std::vector & eqv_strain_incr) +{ + return _eqv_creep_strain_L_old[_qp] + eqv_strain_incr[0]; +} + +ADReal +BVBlancoMartinModelUpdate::munsondawsonCreepStrain(const std::vector & eqv_strain_incr) +{ + return _eqv_creep_strain_R_old[_qp] + eqv_strain_incr[1]; +} + +void +BVBlancoMartinModelUpdate::preReturnMap() +{ + _eqv_creep_strain_L[_qp] = _eqv_creep_strain_L_old[_qp]; + _eqv_creep_strain_R[_qp] = _eqv_creep_strain_R_old[_qp]; +} + +void +BVBlancoMartinModelUpdate::postReturnMap(const std::vector & eqv_strain_incr) +{ + _eqv_creep_strain_L[_qp] = lemaitreCreepStrain(eqv_strain_incr); + _eqv_creep_strain_R[_qp] = munsondawsonCreepStrain(eqv_strain_incr); +} \ No newline at end of file diff --git a/src/materials/BVModifiedLemaitreModelUpdate.C b/src/materials/BVModifiedLemaitreModelUpdate.C new file mode 100644 index 0000000..56562f3 --- /dev/null +++ b/src/materials/BVModifiedLemaitreModelUpdate.C @@ -0,0 +1,117 @@ +/******************************************************************************/ +/* This file is part of */ +/* BEAVER, a MOOSE-based application */ +/* Multiphase Flow Poromechanics for Induced Seismicity Problems */ +/* */ +/* Copyright (C) 2024 by Antoine B. Jacquey */ +/* Polytechnique Montréal */ +/* */ +/* Licensed under GNU Lesser General Public License v2.1 */ +/* please see LICENSE for details */ +/* or http://www.gnu.org/licenses/lgpl.html */ +/******************************************************************************/ + +#include "BVModifiedLemaitreModelUpdate.h" + +registerMooseObject("BeaverApp", BVModifiedLemaitreModelUpdate); + +InputParameters +BVModifiedLemaitreModelUpdate::validParams() +{ + InputParameters params = BVCreepUpdateBase::validParams(); + params.addClassDescription( + "Material for computing a modified Lemaitre creep update (see Blanco-Martin 2024)."); + params.addRequiredRangeCheckedParam("alpha", "0.0 < alpha & alpha < 1.0", "The alpha parameter."); + params.addRequiredRangeCheckedParam("kr1", "kr1 > 0.0", "The kr1 parameter."); + params.addRequiredRangeCheckedParam("kr2", "kr2 > 0.0", "The kr2 parameter."); + params.addRequiredRangeCheckedParam("beta1", "beta1 > 0.0", "The beta1 parameter."); + params.addRequiredRangeCheckedParam("beta2", "beta2 > 0.0", "The beta2 parameter."); + return params; +} + +BVModifiedLemaitreModelUpdate::BVModifiedLemaitreModelUpdate(const InputParameters & parameters) + : BVCreepUpdateBase(parameters), + // Modified Lemaitre creep strain rate parameters + _alpha(getParam("alpha")), + _kr1(getParam("kr1")), + _kr2(getParam("kr2")), + _beta1(getParam("beta1")), + _beta2(getParam("beta2")), + // Internal variable for creep strain + _eqv_creep_strain(declareADProperty(_base_name + "eqv_creep_strain")), + _eqv_creep_strain_old(getMaterialPropertyOld(_base_name + "eqv_creep_strain")) +{ +} + +void +BVModifiedLemaitreModelUpdate::initQpStatefulProperties() +{ + _eqv_creep_strain[_qp] = 0.0; +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRateR(const ADReal & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * eqv_strain_incr; + + if (q == 0.0) + return 0.0; + else + return 1.0e-06 * std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha); +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRate(const ADReal & eqv_strain_incr) +{ + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateR(eqv_strain_incr); + else + return _alpha * creepRateR(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRateRDerivative(const ADReal & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * eqv_strain_incr; + + if (q == 0.0) + return 1.0; + else + return -1.0e-06 * 3.0 * _G / _alpha * + std::pow(std::pow(q / _kr1, _beta1) + std::pow(q / _kr2, _beta2), 1.0 / _alpha - 1.0) * + (_beta1 / _kr1 * std::pow(q / _kr1, _beta1 - 1.0) + + _beta2 / _kr2 * std::pow(q / _kr2, _beta2 - 1.0)); +} + +ADReal +BVModifiedLemaitreModelUpdate::creepRateDerivative(const ADReal & eqv_strain_incr) +{ + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); + else + return std::pow(gamma_l, -1.0 / _alpha) * + (_alpha * gamma_l * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * (_alpha - 1.0) * creepRateR(eqv_strain_incr)); +} + +ADReal +BVModifiedLemaitreModelUpdate::lemaitreCreepStrain(const ADReal & eqv_strain_incr) +{ + return _eqv_creep_strain_old[_qp] + eqv_strain_incr; +} + +void +BVModifiedLemaitreModelUpdate::preReturnMap() +{ + _eqv_creep_strain[_qp] = _eqv_creep_strain_old[_qp]; +} + +void +BVModifiedLemaitreModelUpdate::postReturnMap(const ADReal & eqv_strain_incr) +{ + _eqv_creep_strain[_qp] = lemaitreCreepStrain(eqv_strain_incr); +} \ No newline at end of file diff --git a/src/materials/BVRTL2020ModelUpdate.C b/src/materials/BVRTL2020ModelUpdate.C index a0bcca4..970efd4 100644 --- a/src/materials/BVRTL2020ModelUpdate.C +++ b/src/materials/BVRTL2020ModelUpdate.C @@ -22,9 +22,11 @@ BVRTL2020ModelUpdate::validParams() params.addClassDescription( "Material for computing a RTL2020 creep update. See Azabou et al. (2021), Rock salt " "behavior: From laboratory experiments to pertinent long-term predictions."); + // Lemaitre creep strain rate parameters params.addRequiredRangeCheckedParam("alpha", "0.0 < alpha & alpha < 1.0", "The alpha parameter."); params.addRequiredRangeCheckedParam("A2", "A2 > 0.0", "The A2 parameter."); params.addRequiredRangeCheckedParam("n2", "n2 > 0.0", "The n2 parameter."); + // Munson-Dawson creep strain rate parameters params.addRequiredRangeCheckedParam("A1", "A1 > 0.0", "The A1 parameter."); params.addRequiredRangeCheckedParam("n1", "n1 > 0.0", "The n1 parameter."); params.addRequiredRangeCheckedParam("A", "A >= 0.0", "The A parameter."); @@ -73,36 +75,42 @@ BVRTL2020ModelUpdate::creepRate(const std::vector & eqv_strain_incr, con throw MooseException("BVRTL2020ModelUpdate: error, unknow creep model called in `creepRate`!"); } +ADReal +BVRTL2020ModelUpdate::creepRateR(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 0.0; + else + return 1.0e-06 * std::pow(q / _A2, _n2); +} + ADReal BVRTL2020ModelUpdate::creepRateLemaitre(const std::vector & eqv_strain_incr) { - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2); + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + + if (gamma_l == 0.0) + return _alpha * creepRateR(eqv_strain_incr); else - return _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); + return _alpha * creepRateR(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); } ADReal BVRTL2020ModelUpdate::creepRateMunsonDawson(const std::vector & eqv_strain_incr) { - ADReal saturation_strain = - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; - ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); if (gamma_ms < saturation_strain) return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2); + creepRateR(eqv_strain_incr); else - return -_B * std::pow(gamma_ms / saturation_strain, _m) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2); + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateR(eqv_strain_incr); } ADReal @@ -119,33 +127,37 @@ BVRTL2020ModelUpdate::creepRateDerivative(const std::vector & eqv_strain "BVRTL2020ModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); } +ADReal +BVRTL2020ModelUpdate::creepRateRDerivative(const std::vector & eqv_strain_incr) +{ + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + + if (q == 0.0) + return 1.0; + else + return -1.0e-06 * 3.0 * _G * _n2 / _A2 * std::pow(q / _A2, _n2 - 1.0); +} + ADReal BVRTL2020ModelUpdate::creepRateLemaitreDerivative(const std::vector & eqv_strain_incr, const unsigned int j) { + ADReal gamma_l = 1.0e+06 * lemaitreCreepStrain(eqv_strain_incr); + if (j == 0) // Lemaitre wrt Lemaitre - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return -3.0 * _G / _A2 * _n2 * _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); else - return _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), -1.0 / _alpha) * - (-3.0 * _G / _A2 * _n2 * lemaitreCreepStrain(eqv_strain_incr) + - (1.0 - 1.0 / _alpha) * - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2); + return std::pow(gamma_l, -1.0 / _alpha) * + (_alpha * gamma_l * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * (_alpha - 1.0) * creepRateR(eqv_strain_incr)); + else if (j == 1) // Lemaitre wrt Munson-Dawson - if (lemaitreCreepStrain(eqv_strain_incr) == 0.0) - return -3.0 * _G / _A2 * _n2 * _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + if (gamma_l == 0.0) + return _alpha * creepRateRDerivative(eqv_strain_incr); else - return -3.0 * _G / _A2 * _n2 * _alpha * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - std::pow(lemaitreCreepStrain(eqv_strain_incr), 1.0 - 1.0 / _alpha); + return _alpha * creepRateRDerivative(eqv_strain_incr) * std::pow(gamma_l, 1.0 - 1.0 / _alpha); + else throw MooseException( "BVRTL2020ModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); @@ -155,34 +167,29 @@ ADReal BVRTL2020ModelUpdate::creepRateMunsonDawsonDerivative(const std::vector & eqv_strain_incr, const unsigned int j) { - ADReal saturation_strain = - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A1, _n1); + ADReal q = _eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1]); + ADReal saturation_strain = (q != 0.0) ? std::pow(q / _A1, _n1) : 1.0e+06; - ADReal gamma_ms = munsondawsonCreepStrain(eqv_strain_incr); + ADReal gamma_ms = 1.0e+06 * munsondawsonCreepStrain(eqv_strain_incr); if (j == 0) // Munson-Dawson wrt Lemaitre if (gamma_ms < saturation_strain) - return -3.0 * _G / _A2 * _n2 * _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n) * + creepRateRDerivative(eqv_strain_incr); else - return 3.0 * _G / _A2 * _n2 * _B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0); + return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m) * + creepRateRDerivative(eqv_strain_incr); + else if (j == 1) // Munson-Dawson wrt Munson-Dawson if (gamma_ms < saturation_strain) - return -_A * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - (_n / saturation_strain * - (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2 + - 3.0 * _G / _A2 * _n2 * (1.0 - gamma_ms / saturation_strain)); + return _A * std::pow(1.0 - gamma_ms / saturation_strain, _n - 1.0) * + ((1.0 - gamma_ms / saturation_strain) * creepRateRDerivative(eqv_strain_incr) - + 1.0e+06 * _n / saturation_strain * creepRateR(eqv_strain_incr)); else return -_B * std::pow(gamma_ms / saturation_strain - 1.0, _m - 1.0) * - std::pow((_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2, - _n2 - 1.0) * - (_m * (_eqv_stress_tr - 3.0 * _G * (eqv_strain_incr[0] + eqv_strain_incr[1])) / _A2 - - 3.0 * _G / _A2 * _n2 * (gamma_ms / saturation_strain - 1.0)); + ((gamma_ms / saturation_strain - 1.0) * creepRateRDerivative(eqv_strain_incr) + + 1.0e+06 * _m / saturation_strain * creepRateR(eqv_strain_incr)); + else throw MooseException( "BVRTL2020ModelUpdate: error, unknow creep model called in `creepRateDerivative`!"); diff --git a/test/tests/viscoelasticity/blanco-martin-lemaitre.i b/test/tests/viscoelasticity/blanco-martin-lemaitre.i new file mode 100644 index 0000000..52ce1cd --- /dev/null +++ b/test/tests/viscoelasticity/blanco-martin-lemaitre.i @@ -0,0 +1,257 @@ +# Modified Lemaitre creep model +# See Blanco-Martin et al. (2024) +# Parameters +# Units: stress in MPa, time in days, strain in m / m +E = 12000 +nu = 0.3 +alpha = 0.326 +kr1 = 0.7 +beta1 = 2.922 +kr2 = 0.009 +beta2 = 0.867 +P = 5.0 +Q = 5.0 +Q1 = 5.5 + +[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 + 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' + [] +[] + +[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 = ${Q1} + [] + [] +[] + +[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} + # [] + [viscoelastic] + type = BVModifiedLemaitreModelUpdate + alpha = ${alpha} + kr1 = ${kr1} + beta1 = ${beta1} + kr2 = ${kr2} + beta2 = ${beta2} + [] +[] + +[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 = 1 + dt = 0.02 + timestep_tolerance = 1.0e-10 +[] + +[Outputs] + perf_graph = true + exodus = true +[] \ No newline at end of file diff --git a/test/tests/viscoelasticity/blanco-martin-rtl.i b/test/tests/viscoelasticity/blanco-martin-rtl.i new file mode 100644 index 0000000..0fc483f --- /dev/null +++ b/test/tests/viscoelasticity/blanco-martin-rtl.i @@ -0,0 +1,253 @@ +# 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 +Q1 = 5.5 + +[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' + [] +[] + +[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 = ${Q1} + [] + [] +[] + +[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} + [] +[] + +[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 = 1 # + dt = 0.02 + timestep_tolerance = 1.0e-10 +[] + +[Outputs] + perf_graph = true + exodus = true +[] \ No newline at end of file diff --git a/test/tests/viscoelasticity/gold/blanco-martin-lemaitre_out.e b/test/tests/viscoelasticity/gold/blanco-martin-lemaitre_out.e new file mode 100644 index 0000000..922e286 Binary files /dev/null and b/test/tests/viscoelasticity/gold/blanco-martin-lemaitre_out.e differ diff --git a/test/tests/viscoelasticity/gold/blanco-martin-rtl_out.e b/test/tests/viscoelasticity/gold/blanco-martin-rtl_out.e new file mode 100644 index 0000000..629c636 Binary files /dev/null and b/test/tests/viscoelasticity/gold/blanco-martin-rtl_out.e differ diff --git a/test/tests/viscoelasticity/tests b/test/tests/viscoelasticity/tests index af8e954..726fbd4 100644 --- a/test/tests/viscoelasticity/tests +++ b/test/tests/viscoelasticity/tests @@ -19,4 +19,14 @@ input = 'lubby2.i' exodiff = 'lubby2_out.e' [] + [blanco-martin-lemaitre] + type = 'Exodiff' + input = 'blanco-martin-lemaitre.i' + exodiff = 'blanco-martin-lemaitre_out.e' + [] + [blanco-martin-rtl] + type = 'Exodiff' + input = 'blanco-martin-rtl.i' + exodiff = 'blanco-martin-rtl_out.e' + [] [] \ No newline at end of file