diff --git a/doc/content/bib/blackbear.bib b/doc/content/bib/blackbear.bib index 29c911b6..f5fa6ac0 100644 --- a/doc/content/bib/blackbear.bib +++ b/doc/content/bib/blackbear.bib @@ -244,3 +244,50 @@ @article{poyet2009temperature year={2009}, publisher={Elsevier} } + + @article{lee1998plastic, + title={Plastic-damage model for cyclic loading of concrete structures}, + author={Lee, Jeeho and Fenves, Gregory L}, + journal={Journal of engineering mechanics}, + volume={124}, + number={8}, + pages={892--900}, + year={1998}, + publisher={American Society of Civil Engineers} +} + +@book{lee1996theory, + title={Theory and implementation of plastic-damage model for concrete structures under cyclic and dynamic loading}, + author={Lee, Jeeho}, + year={1996}, + publisher={University of California, Berkeley} +} + +@article{wilkins2020method, + title = {A method for smoothing multiple yield functions}, + author = {Andy Wilkins and Benjamin W. Spencer and Amit Jain and Bora Gencturk}, + year = {2020}, +journal = {International Journal for Numerical Methods in Engineering}, + volume = {121}, + number = {3}, + pages = {434--449}, + doi = {10.1002/nme.6215} +} + +@article{lubliner1989plastic, + title={A plastic-damage model for concrete}, + author={Lubliner, Jacob and Oliver, Javier and Oller, Sand and Onate, EJIJos}, + journal={International Journal of solids and structures}, + volume={25}, + number={3}, + pages={299--326}, + year={1989}, + publisher={Elsevier} +} + +@article{krabbenhoft2002basic, + title={Basic computational plasticity}, + author={Krabbenh{\o}ft, KRISTIAN}, + journal={University of Denmark}, + year={2002} +} diff --git a/doc/content/media/Return_mapping_flow_chart.png b/doc/content/media/Return_mapping_flow_chart.png new file mode 100644 index 00000000..aebf3738 Binary files /dev/null and b/doc/content/media/Return_mapping_flow_chart.png differ diff --git a/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md b/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md new file mode 100644 index 00000000..eeaceb64 --- /dev/null +++ b/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md @@ -0,0 +1,21 @@ +# Compute Multiple Inelastic Damage Stress + +!syntax description /Materials/ComputeMultipleInelasticDamageStress + +## Description + +`ComputeMultipleInelasticDamageStress` computes the damage stress. + +## Example Input Files + +The input settings for the inelastic material model is as follows + +!listing test/tests/damage_plasticity_model/uniaxial_test.i block=Materials/stress + +!syntax parameters /Materials/ComputeMultipleInelasticDamageStress + +!syntax inputs /Materials/ComputeMultipleInelasticDamageStress + +!syntax children /Materials/ComputeMultipleInelasticDamageStress + +!bibtex bibliography diff --git a/doc/content/source/materials/DamagePlasticityStressUpdate.md b/doc/content/source/materials/DamagePlasticityStressUpdate.md new file mode 100644 index 00000000..689cd492 --- /dev/null +++ b/doc/content/source/materials/DamagePlasticityStressUpdate.md @@ -0,0 +1,346 @@ +# Damage Plasticity Model + +The [!cite](lee1996theory) model accounts for the independent damage in tension and compression. It also accounts for degradation of the elastic modulus of the concrete as the loading goes beyond yielding in either tension or compression. The model uses the incremental theory of plasticity and decomposes the total strain, $\boldsymbol{\varepsilon}$, into elastic strain, $\boldsymbol{\varepsilon}^{e}$, and plastic strain, $\boldsymbol{\varepsilon}^{p}$, as follows +\begin{equation} + \boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^{e} + \boldsymbol{\varepsilon}^{p} \label{eps_def} +\end{equation} +where bold symbol represents a vectoral or tensorial quantity. The relation between elastic strain and the stress, $\boldsymbol{\sigma}$, is given by +\begin{equation} + \boldsymbol{\varepsilon}^{e} = \boldsymbol{\mathfrak{E}}^{-1}:\boldsymbol{\sigma} \label{eps_e_def} +\end{equation} +where $\boldsymbol{\mathfrak{E}}$ is the elasticity tensor. Using Eqs. \eqref{eps_def}-\eqref{eps_e_def}, the relation between $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}^{p}$ is expressed as +\begin{equation} + \boldsymbol{\sigma} = \boldsymbol{\mathfrak{E}}:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{e}\right) +\end{equation} +Since the model considers the effect of damage in elastic stiffness, an effective stress, +$\boldsymbol{\sigma}^{e}$, is defined, where the stress for a given strain always corresponds to the +undamaged elastic stiffness of the material, $\boldsymbol{\mathfrak{E}}_{0}$ The relation between +$\boldsymbol{\sigma}^{e}$, $\boldsymbol{\varepsilon}$, and $\boldsymbol{\varepsilon}^{p}$ is given by +\begin{equation} + \boldsymbol{\sigma}^e = \boldsymbol{\mathfrak{E}}_0:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{e}\right) +\end{equation} +To consider the degradation of reinforced-concrete structures, an isotropic damage was +considered in concrete material. Hence, the relation between $\boldsymbol{\sigma}^e$ and $\boldsymbol{\sigma}$ can be established by +the isotropic scalar degradation damage variable, D, as follows +\begin{equation} + \boldsymbol{\sigma} = \left(1-D\right)\boldsymbol{\sigma}^e \label{sigma_def} +\end{equation} +\begin{equation} + \boldsymbol{\sigma} = \left(1-D\right)\boldsymbol{\mathfrak{E}}_0:\left(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{e}\right)\label{sigma_def2} +\end{equation} +The Damage Plasticity Model has various attributes to define the mechanical behavior of concrete +in tension and compression such as the yield function, plastic potential, strength of material +in tension and compression, and hardening and softening of the yield surface. These attributes +are discussed in detail in the following sections. A method for the implementation of the +Damage Plasticity Model and for the estimation of crack width are also presented in the upcoming +sections. + +## Yield Function + +The yield function, $\mathfrak{F}$ is a function of $\boldsymbol{\sigma}$, the strength of the material in uniaxial tension, $f_t$, and the strength of the material in uniaxial compression, $f_c$. It was used to describe the admissible stress space. For this implementation, the yield function in stress space is defined as follows +\begin{equation} \label{yf} +\begin{gathered} + \mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right) = \frac{1}{1-\alpha} \\ + \left(\alpha I_1 + \sqrt{3J_2} + \beta\left(\boldsymbol{\kappa}\right)<{\hat{\boldsymbol{\sigma}}_{max}}>\right) - f_c\left(\boldsymbol{\kappa}\right) +\end{gathered} +\end{equation} +where $I_1$ and $J_2$ is first invariant of stress and second invariant of the deviatoric +component of the stress, respectively, $ =\frac{x+|x|}{2}$ is the Macaulay bracket function, ${\hat{\boldsymbol{\sigma}}_{max}}$ is algebraically maximum principal +stress, $\alpha = \frac{f_{b0}-f_{c0}}{2 f_{b0}-f_{c0}}$ is a parameter that relates +uniaxial, $f_{c0}$, and biaxial, $f_{b0}$, yield strength of concrete in compression, +$\beta\left(\boldsymbol{\kappa}\right)=\frac{f_c\left(\boldsymbol{\kappa}\right)}{f_t\left(\boldsymbol{\kappa}\right)}\left(\alpha-1\right)-\left(1+\alpha\right)$ is a parameter that +relates tensile, $f_t\left(\boldsymbol{\kappa}\right)$, and compressive, $f_c\left(\boldsymbol{\kappa}\right)$, yield strength which are +function of a vector of damage variable, $\boldsymbol{\kappa} = \{\kappa_t, \kappa_c\}$ and $\kappa_t$ +and $\kappa_c$ are the damage variables in tension and compression, respectively. + +The implementation first solves the given problem in the effective stress space and then transform the effective stress to stress space using Eq. \eqref{sigma_def2}. Thus, the yield strength of the concrete under uniaxial loading is expressed as effective yield strength as follows +\begin{equation} + f_t\left(\boldsymbol{\kappa}\right) = \left(1-D_t \left(\kappa_t\right)\right)f_{t}^{e}\left(\kappa_t\right) \label{ft} +\end{equation} +\begin{equation} + f_c\left(\boldsymbol{\kappa}\right) = \left(1-D_c \left(\kappa_c\right)\right)f_{c}^{e}\left(\kappa_c\right) \label{fc} +\end{equation} +where $f_{t}^{e}$ and $f_{c}^{e}$ are the yield strength of the concrete in tension and +compression, respectively and $D_t$ and $D_c$ are the degradation damage variables in +tension and compression, respectively such that $0\leq D_t$\textless 1 and $0\leq D_c$\textless 1. +The scalar degradation damage variable is expressed in terms of $D_t$ and $D_c$ as follows +\begin{equation} + D\left(\boldsymbol{\kappa}\right) = 1-\left(1-D_t\left(\kappa_t\right)\right)\left(1-D_c\left(\kappa_c\right)\right) \label{D} +\end{equation} +Hence, for uniaxial tension, $D=D_t$, while for uniaxial compression, $D=D_c$.The yield strength for multi-axial loading, i.e., Eqs. \eqref{ft}-\eqref{fc}, can be rewritten as +\begin{equation} + f_t\left(\boldsymbol{\kappa}\right) = \left(1-D\left(\boldsymbol{\kappa}\right)\right)f_{t}^{e}\left(\kappa_t\right) \label{ft_new} +\end{equation} +\begin{equation} + f_c\left(\boldsymbol{\kappa}\right) = \left(1-D\left(\boldsymbol{\kappa}\right)\right)f_{c}^{e}\left(\kappa_c\right) \label{fc_new} +\end{equation} +Similarly, the first invariant of $\boldsymbol{\sigma}^e$, $I_1^e$, and second invariant of the deviatoric component of $\boldsymbol{\sigma}^e$, $J_2^e$, can be rewritten in terms of $I_1$ and $J_2$ as follows +\begin{equation} + I_1^e = \left(1-D\left(\boldsymbol{\kappa}\right)\right)I_1 \label{I1e} +\end{equation} +\begin{equation} + J_2^e = \left(1-D\left(\boldsymbol{\kappa}\right)\right)^2J_2 \label{J2e} +\end{equation} +Since $D$ \textless 1, the maximum principal effective stress ${\hat{\boldsymbol{\sigma}}_{max}}^e$ is expressed in the terms of ${\hat{\boldsymbol{\sigma}}_{max}}$ as follows +\begin{equation} + {\hat{\boldsymbol{\sigma}}_{max}}^e = \left(1-D\left(\boldsymbol{\kappa}\right)\right){\hat{\boldsymbol{\sigma}}_{max}} \label{sig_max_e} +\end{equation} +Consequently, yield function $\left(\mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right)\right)$ is a homogenous +function, i.e., $x \mathfrak{F}\left(\boldsymbol{\sigma},f_t,f_c\right) = \mathfrak{F}\left(x \boldsymbol{\sigma},x f_t,x f_c\right)$ Hence, using Eqs. \eqref{ft_new}-\eqref{sig_max_e}, the yield function in the effective stress space was obtained by multiplying by a factor $\left(1-D\right)$ of both sides of Eq. \eqref{yf}, as follows +\begin{equation}\label{yf_e} +\begin{gathered} + \mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) = \frac{1}{1-\alpha} \\ + \left(\alpha I_1^e + \sqrt{3J_2^e} + \beta\left(\boldsymbol{\kappa}\right)<{\hat{\boldsymbol{\sigma}}_{max}}^e>\right) - f_c^e\left(\boldsymbol{\kappa}\right) +\end{gathered} +\end{equation} + +## Plastic Potential + +It was found that for concrete, the Drucker-Prager flow rule describes the experimentally observed volumetric expansion of the material as opposed to the Von-Mises flow rule, which results in no volumetric expansion ([!cite](krabbenhoft2002basic)). Since all the equations are solved in the effective stress space, the plastic flow potential is also defined in the effective stress space ([!cite](lee1996theory)) as follows +\begin{equation} + \dot{\boldsymbol{\varepsilon}^p} = \dot{\gamma} \nabla_{\boldsymbol{\sigma}^e}\Phi\left(\boldsymbol{\sigma}^e\right) +\end{equation} +\begin{equation} \label{flowRule} + \Phi\left(\boldsymbol{\sigma}^e\right)=\alpha_p I_1^e+\|{s^e}\| +\end{equation} +where $\Phi$ is the plastic potential function, $s^e$ is the deviatoric component of the +$\boldsymbol{\sigma}^e$, and $\|\cdot\|$ is $L_2$ norm of $\alpha_p$ is a parameter that controls the +dilatancy of concrete, and $\dot{\gamma}$ is the plastic consistency parameter. + +## Strength Function + +Since the concrete shows strain-softening in tension and strain hardening and softening in compression, the concrete strength is expressed as a combination of two exponential functions as follows +\begin{equation} + f_N = f_{N0} \left(\left(1+a_N\right) e^{-b_N \varepsilon^p}- a_N e^{-2b_N \varepsilon^p}\right) \label{fN} +\end{equation} +where $f_{N0}$ is the initial yield stress of the material, $N = t$, for the uniaxial +tension, $N = c$, for uniaxial compression, $a_N$ and $b_N$, are the material constants +that describe the softening and hardening behavior of the concrete. Similarly, the +degradation of the elastic modulus is also expressed as another exponential function as +follows +\begin{equation} + D_N = 1 - e^{-d_N \varepsilon^p} \label{DN} +\end{equation} +where $d_N$ is a constant that determine the rate of degradation of $\boldsymbol{\mathfrak{E}}$ with the +increase in plastic strain. The strength of the material in the effective stress space was +obtained using Eqs. \eqref{ft_new}-\eqref{fc_new}, and \eqref{fN}-\eqref{DN}, as follows +\begin{equation}\label{fNe} + f_N^e = f_{N0} \left(\left(1+a_N\right) \left(e^{-b_N \varepsilon^p}\right)^{1-\frac{d_N}{b_N}}- + a_N \left(e^{-b_N \varepsilon^p}\right)^{2-\frac{d_N}{b_N}}\right) +\end{equation} +The damage variable, $\kappa_N$ is defined as +\begin{equation} + \kappa_N = \frac{1}{g_N}\int_0^{\varepsilon^p} {f_N\left(\varepsilon^p\right)d\varepsilon^p} \label{kN_def} +\end{equation} +where $g_N$ $\left(=\int_0^\infty {f_N\left(\varepsilon^p\right)d\varepsilon^p}=\frac {f_{N0}}{b_N}\left(1+\frac{a_N}{2}\right)\right)$ +is the fracture energy density during the process of cracking, which is derived from the +fracture energy, $G_N$, which is a material property. The relation between $G_N$ and $g_N$ +is expressed as follows +\begin{equation} + g_N = \frac{G_N}{l_N} \label{GN_def} +\end{equation} +where $l_N$ is characteristic length or the size of the deformation localization zone. +Thus, the plastic strain can be presented in terms of damage variable as follows +\begin{equation} + \varepsilon^p = \frac{1}{b_N} \log{\frac{\sqrt{\Phi_N}}{a_N}} \label{eps_p} +\end{equation} +where $\Phi_N = 1 + a_N \left(2+a_N \right)\kappa_N$. Using Eqs. \eqref{fN} and \eqref{eps_p}, the +strength of the concrete can be expressed in terms of the damage variable as follows +\begin{equation} + f_N = f_{N0} \frac{1+a_N-\sqrt{\Phi_N\left(\kappa_N\right)}}{a_N}\sqrt{\Phi_N\left(\kappa_N\right)} \label{fN_new} +\end{equation} +Thus, the strength of the material and degradation damage variable in the effective stress space can be written as + +\begin{equation} + f_N^e = f_{N0} \left(\frac{1+a_N-\sqrt{\Phi_N\left(\kappa_N\right)}}{a_N}\right)^{1-\frac{d_N}{b_N}} \sqrt{\Phi_N\left(\kappa_N\right)} \label{fNe_new} +\end{equation} +\begin{equation} + D_N = 1- \left(\frac{1+a_N-\sqrt{\Phi_N\left(\kappa_N\right)}}{a_N}\right)^{\frac{d_N}{b_N}} \label{DN_new} +\end{equation} +where $a_N$, $b_N$,and $d_N$ are the modeling parameters, which are evaluated from the +material properties. Since the maximum compressive strength of concrete, $f_{cm}$, was used +as a material property, $f_{cm}$ was obtained in terms of $a_c$ by finding maximum value of +compressive strength in Eq. \eqref{fNe} as follows +\begin{equation} + f_{cm} = \frac{f_{c0}\left(1+a_c\right)^2}{4a_c} \label{fcm} +\end{equation} +Thus, $a_c$ can be expressed as follows +\begin{equation} + a_c = 2\frac{f_{cm}}{f_{c0}}-1+2\sqrt{\left(\frac{f_{cm}}{f_{c0}}\right)^2-\frac{f_{cm}}{f_{c0}}} \label{ac} +\end{equation} +Similarly, if $G_c$ and $l_c$ are known then $b_c$ can be expressed in term of known quantities as follows +\begin{equation} + b_c = \frac{f_{c0}}{\frac{G_c}{l_c}}\left(1+\frac{a_c}{2}\right) \label{bc} +\end{equation} +A relationship between $a_t$ and $b_t$ is written as follows +\begin{equation} + b_t = \frac{f_{t0}}{\frac{G_t}{l_t}}\left(1+\frac{a_t}{2}\right) \label{bt} +\end{equation} +[!cite](lubliner1989plastic) suggested that if the slope of $\sigma$ versus $\varepsilon^p$ curve is +known at $\varepsilon^p=0$, then another relationship between $a_t$ and $b_t$ will be obtained +as follows +\begin{equation} + \left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0} = f_{t0}b_t\left(a_t-1\right) \label{slope} +\end{equation} +Thus, $a_t$ was obtained using Eqs. \eqref{bt}-\eqref{slope} as follows +\begin{equation} + a_t = \sqrt{\frac{9}{4}+\frac{2\frac{G_t}{l_t} \left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0}}{f_{t0}^2}}\label{at} +\end{equation} +The minimum slope of the $\sigma$ versus $\varepsilon^p$ curve is +$\left(\left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0}\right)_{min}= +-\frac{9}{8}\frac{f_{t0}^2}{\frac{G_t}{l_t}}$, which is a function of the characteristic length in tension. +Therefore, a mesh independent slope parameter $\omega\in\left(0,1\right)$, is defined such that +\begin{equation} + \left(\frac{d\boldsymbol{\sigma}}{d\varepsilon^p}\right)_{\varepsilon^p=0} = \omega \left(\left(\frac{d\sigma}{d\varepsilon^p}\right)_{\varepsilon^p=0}\right)_{min} \label{slope_new} +\end{equation} +Using Eqs. \eqref{at}-\eqref{slope_new}, $a_t$ is rewritten as follows +\begin{equation} + a_t = \frac{3}{2}\sqrt{1-\omega}-\frac{1}{2}\label{at_new} +\end{equation} +The ratio of $\frac{d_c}{b_c}$ was obtained by specifying degradation values for uniaxial +compression case from experiments. If the degradation in the elastic modulus is known, +denoted as $\widetilde{D}_c$, when the concrete is unloaded from $\sigma =f_{cm}$, then $\frac{d_c}{b_c}$ will be obtained using the following relation +\begin{equation} + \widetilde{D}_c = 1 - \left(\frac{1+a_c}{2a_c}\right)^{\frac{d_c}{b_c}} \label{Dc_fcm} +\end{equation} +\begin{equation} + \frac{d_c}{b_c} = \frac{\log\left(1-\widetilde{D}_c\right)}{\log\left(\frac{1+a_c}{2a_c}\right)} \label{dcbc_fcm} +\end{equation} +Similarly, if degradation in the elastic modulus is known, denoted as $\widetilde{D}_t$, when the material is unloaded from $\sigma=\frac {f_{t0}}{2}$, on softening branch, then $\frac{d_t}{b_t}$ will be obtained using the following relation +\begin{equation} + \widetilde{D}_t = 1 - \left(\frac{1+a_t-\sqrt{1+a_t^2}}{2a_t}\right)^{\frac{d_t}{b_t}} \label{Dt_ft0_2} +\end{equation} +\begin{equation} + \frac{d_t}{b_t} = \frac{\log\left(1-\widetilde{D}_t\right)}{\log\left(\frac{1+a_t-\sqrt{1+a_t^2}}{2a_t}\right)} \label{Dt_ft0} +\end{equation} +Thus, material modeling parameters $a_N$,$b_N$, and $d_N$ were obtained, which were used in +defining the strength of concrete in both tension and compression as given in Eq. +\eqref{fNe_new}. These parameters are also used to define the degradation damage variable in +both tension and compression as indicated in Eq. \eqref{DN_new}. + +## Hardening Potential + +The vector of two damage variables, $\boldsymbol{\kappa}=\{\kappa_t, \kappa_c\}$, was used in the implementation as the state variable to store the state of damage in tension and compression, separately. The evolution of these damage variables is defined in terms of the hardening potential, $H$, as +\begin{equation} + \dot{\boldsymbol{\kappa}} = \dot{\gamma} H\left(\boldsymbol{\sigma}^e, \boldsymbol{\kappa}\right)\label{kappa} +\end{equation} +The evolution of the damage variable is expressed in terms of the evolution of $\boldsymbol{\varepsilon}^p$ as follows +\begin{equation} + \dot{\boldsymbol{\kappa}} = \frac{1}{g_N}f_N^e\left(\kappa_N\right)\dot{\boldsymbol{\varepsilon}^p} \label{kappa_ep} +\end{equation} +where $g_N$ is dissipated energy density during the process of cracking. The scalar $\dot{\boldsymbol{\varepsilon}^p}$, is extended to multi-dimensional case as follows +\begin{equation} + \dot{\boldsymbol{\varepsilon}^p} = \delta_{tN} r\left(\hat{\boldsymbol{\sigma}^e}\right)\dot{\varepsilon}^{p}_{max} + \delta_{cN} \left(1-r\left(\hat{\boldsymbol{\sigma}^e}\right)\right)\dot{\varepsilon}^{p}_{min} \label{ep_dot} +\end{equation} +where $\delta_{ij}$ is the Dirac delta function and $\hat{\boldsymbol{\sigma}^e}$ are eigenvalues of the $\boldsymbol{\sigma}^e$, +\begin{equation}\label{r_sige} + r\left(\hat{\boldsymbol{\sigma}^e}\right) = + \begin{cases} + 0,& \text{if } \boldsymbol{\sigma}^e = \boldsymbol{0}\\ + \frac{\sum_{i=1}^3<\sigma^e_i>}{\sum_{i=1}^3|\sigma^e_i|}, & \text{otherwise} + \end{cases} +\end{equation} +$\dot{\varepsilon}^{p}_{max}$ and $\dot{\varepsilon}^{p}_{min}$ +are the maximum and minimum principal plastic strain, respectively. From Eqs. \eqref{kappa_ep} - \eqref{r_sige}, the evolution of $\boldsymbol{\kappa}$ was obtained as +\begin{equation} + \dot{\boldsymbol{\kappa}} = \boldsymbol{h}\left(\hat{\boldsymbol{\sigma}^e}\right):\dot{\boldsymbol{\varepsilon}}^{\hat{p}} \label{kappa_h_ep} +\end{equation} +where +\begin{equation}\label{h} + \boldsymbol{h}\left(\hat{\boldsymbol{\sigma}^e}\right)= + \begin{bmatrix} + \frac{r\left(\hat{\boldsymbol{\sigma}^e}\right)}{g_t}f_t^e\left(\kappa_t\right)&0&0\\ + 0&1&0\\ + 0&0&\frac{1-r\left(\hat{\boldsymbol{\sigma}^e}\right)}{g_c}f_c^e\left(\kappa_c\right)\\ + \end{bmatrix} +\end{equation} +and ‘:’ represents products of two matrices. Hence, $H\left(\boldsymbol{\sigma}^e,\boldsymbol{\kappa}\right)$ in Eq. \eqref{kappa} was obtained as follows +\begin{equation} + H\left(\boldsymbol{\sigma}^e, \boldsymbol{\kappa}\right) = \boldsymbol{h}\cdot \nabla_{\hat{\boldsymbol{\sigma}^e}}\Phi\left(\hat{\boldsymbol{\sigma}^e}\right) \label{H_def} +\end{equation} +where ‘$\cdot$’ represents the dot product of a matrix and a vector, and $\nabla_{\hat{\boldsymbol{\sigma}^e}}$ is the gradient with respect to principal effective stress components, $\hat{\boldsymbol{\sigma}^e}$. Thus, the hardening potential that governs the evolution of damage variables is expressed in terms of effective stress space. + +## Return Mapping Algorithm + +Return mapping algorithm is summarized in a flowchart [flowChart]. For given strain increment and previous state of stress, trial effective stress is obtained by elastic increment in the elastic predictor step according to +\begin{equation} + \boldsymbol{\sigma}_{n+1}^{e^{tr}} = \boldsymbol{\mathfrak{E}}_0:\left(\varepsilon_{n+1}-\varepsilon_n^p\right) \label{predictor_step} +\end{equation} +If the state of trial effective stress lies inside the admissible domain, i.e, $\mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) < 0$, the trial effective stress is considered as an admissible effective stress and old damage variables are supplied for the next step. +If the state of trial effective stress is outside admissible domain, i.e, $\mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) > 0$, the effective stress is obtained using the plastic corrector step + +!media media/Return_mapping_flow_chart.png + id=flowChart + style=width:30%; + caption=Flow chart for return mapping algorithm. + +\begin{equation} +\boldsymbol{\sigma}_{n+1}^e = \boldsymbol{\sigma}_{n+1}^{e^{tr}}-\boldsymbol{\mathfrak{E}}_0:\varepsilon_n^p \label{plasticCorrector} +\end{equation} +Once admissible effective stress is obtained, degradation corrector step is utilized to account for the stiffness degradation on the state of stress according to +\begin{equation} +\boldsymbol{\sigma}_{n+1} = \left(1-D_{n}\right)\boldsymbol{\sigma}_{n+1}^{e} \label{degradation_corrector} +\end{equation} + + + +During the plastic corrector step, the returned effective stress should satisfy the Kuhn-Tucker conditions on the $\mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right)$ and $\dot{\gamma}$, according to +\begin{equation} + \begin{split} \label{khunTuckerConditions} + \dot{\gamma} > 0 \\ + \dot{\gamma}\mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) = 0 \\ + \mathfrak{F}\left(\boldsymbol{\sigma}^e,f_t^e,f_c^e\right) = 0 + \end{split} +\end{equation} +As per flow rule in Eq. \eqref{flowRule}, the plastic corrector step, i.e., Eq. \eqref{plasticCorrector} can be rewritten as +\begin{equation} +\boldsymbol{\sigma^e}_{n+1} = \boldsymbol{\sigma}_{n+1}^{e^{tr}}-\dot{\gamma}\left(2G\frac{\boldsymbol{s}_{n+1}^e}{\|\boldsymbol{s}_{n+1}^e\|} + 3K\alpha_p\boldsymbol{I}\right) \label{returnMap1} +\end{equation} +where $G$ is shear modulus and $K$ is bulk modulus. After separating the volumetric and deviatoric components from Eq. \eqref{returnMap1} following relations can be obtained +\begin{equation} + I_{1|n+1} = I_{1|n+1}^{e^{tr}} - 9K\alpha \alpha_p \dot{\gamma} \label{stressRelation1} +\end{equation} +\begin{equation}\label{stressRelation2} + \begin{gathered} + \frac{\boldsymbol{s}_{n+1}^e}{\|\boldsymbol{s}_{n+1}^e\|} = \frac{\boldsymbol{s}_{{n+1}}^{e^{tr}}}{\|\boldsymbol{s}_{n+1}^{e^{tr}}\|} \\ + {\|\boldsymbol{s}^{e}_{n+1}\|} = {\|\boldsymbol{s}_{n+1}^{e^{tr}}\|} - 2G\dot{\gamma} + \end{gathered} +\end{equation} +Using Eqs. \eqref{stressRelation1} and \eqref{stressRelation2}, Eq. \eqref{returnMap1} can be written as +\begin{equation} + \boldsymbol{\sigma}_{n+1}^e = \boldsymbol{\sigma}_{n+1}^{e^{tr}}-\dot{\gamma}\left(2G\frac{\boldsymbol{s}^{e^{tr}}_{n+1}}{\|\boldsymbol{s}_{{n+1}}^{e^{tr}} \|}+ 3K\alpha_p\boldsymbol{I}\right) \label{returnMap2} +\end{equation} +In case of plastic deformation, the returned state of stress should lie on the yield surface as per Kuhn-Tucker conditions (Eq. \eqref{khunTuckerConditions}, therefore $\mathfrak{F}\left(\boldsymbol{\sigma}_{n+1}^e,f_t^e,f_c^e\right) = 0$, i.e., +\begin{equation} \label{yfnext} + \begin{gathered} + \alpha I_{1|n+1}^e + \sqrt{3J_{2|n+1}^e} + \beta\left(\boldsymbol{\kappa}\right)<\hat{\boldsymbol{\sigma}}^e_{n+1|max}> \\ - \left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right) = 0 + \end{gathered} +\end{equation} +Using Eq. \eqref{stressRelation1}, \eqref{stressRelation2}, and \eqref{returnMap2}, Eq. \eqref{yfnext} can be written as +\begin{equation} \label{yfzero} + \begin{gathered} + \alpha\left(I_{1|n+1}^{e^{tr}} - 9K\alpha \alpha_p \dot{\gamma}\right) + + \left(\sqrt{\frac{3}{2}}\|\boldsymbol{s}_{{n+1}}^{e^{tr}}\| - \sqrt{6}G\dot{\gamma}\right)\\+ + \beta\left(\boldsymbol{\kappa}\right)<\hat{\boldsymbol{\sigma}}^e_{n+1|max}> - \left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right) = 0 + \end{gathered} +\end{equation} +Thus, the plastic multiplier can be by solving Eq. \eqref{yfzero} as +\begin{equation}\label{gammaDef} + \dot{\gamma} = + \begin{cases} + \frac{\alpha I_{1|n+1}^{e^{tr}}+\sqrt{\frac{3}{2}}\|\boldsymbol{s}_{{n+1}}^{e^{tr}}\|-\left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right)} + {9K \alpha_p + \sqrt{6}G}, & \text{if $\sigma_{m|n+1}^e < 0$}\\ + \frac{\alpha I_{1|n+1}^{e^{tr}}+\sqrt{\frac{3}{2}}\|\boldsymbol{s}_{{n+1}}^{e^{tr}}\|+\beta\left(\boldsymbol{\kappa}\right) \sigma_{m|n+1}^{e^{tr}}-\left(1-\alpha\right)f_c^e\left(\boldsymbol{\kappa}\right)} + {9K \alpha_p + \sqrt{6}G + \beta\left(\boldsymbol{\kappa}\right)\left(2G\frac{s^{e^{tr}}_{m|n+1}}{\|\boldsymbol{s}_{{n+1}}^{e^{tr}} \|}+ 3K\alpha_p\right)}, & \text{otherwise}. + \end{cases} +\end{equation} +where $\sigma_{m|n+1}^e$, $\sigma_{m|n+1}^{e^{tr}}$, and $s^{e^{tr}}_{m|n+1}$ are the $m^{th}$ component of the $\hat{\boldsymbol{\sigma}}_{n+1}^e$, $\boldsymbol{\sigma}_{n+1}^{e^{tr}}$, and $\boldsymbol{s}^{e^{tr}}_{n+1}$, respectively, which corresponds to maximum principal effective stress in $\left(n+1\right)^{th}$ step. Eq. \eqref{gammaDef} is solved iteratively. + + +!syntax parameters /Materials/DamagePlasticityStressUpdate + +!syntax inputs /Materials/DamagePlasticityStressUpdate + +!syntax children /Materials/DamagePlasticityStressUpdate + +!bibtex bibliography diff --git a/include/materials/ComputeMultipleInelasticDamageStress.h b/include/materials/ComputeMultipleInelasticDamageStress.h new file mode 100644 index 00000000..a569e8b7 --- /dev/null +++ b/include/materials/ComputeMultipleInelasticDamageStress.h @@ -0,0 +1,42 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* BlackBear */ +/* */ +/* (c) 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +#pragma once + +#include "ComputeMultipleInelasticStress.h" + +class ComputeMultipleInelasticDamageStress : public ComputeMultipleInelasticStress +{ +public: + static InputParameters validParams(); + ComputeMultipleInelasticDamageStress(const InputParameters & parameters); + +protected: + /// damage parameter for DamagePlasticityStressUpdate model + const MaterialProperty & _D; + const MaterialProperty & _D_old; + const MaterialProperty & _D_older; + + virtual void computeQpJacobianMult() override; + + virtual void computeAdmissibleState(unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & inelastic_strain_increment, + RankFourTensor & consistent_tangent_operator) override; + + virtual void + updateQpStateSingleModel(unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & combined_inelastic_strain_increment) override; +}; diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h new file mode 100644 index 00000000..51b958a9 --- /dev/null +++ b/include/materials/DamagePlasticityStressUpdate.h @@ -0,0 +1,184 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* BlackBear */ +/* */ +/* (c) 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +#pragma once + +#include "MultiParameterPlasticityStressUpdate.h" + +class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate +{ +public: + static InputParameters validParams(); + DamagePlasticityStressUpdate(const InputParameters & parameters); + + /** + * Does the model require the elasticity tensor to be isotropic? + */ + bool requiresIsotropicTensor() override { return true; } + +protected: + virtual void initQpStatefulProperties() override; + virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment) override; + +private: + const Real _f_tol; + const Real _alfa; + const Real _alfa_p; + const Real _s0; + + const Real _Chi; + const Real _Dt; + const Real _ft; + const Real _FEt; + + const Real _fyc; + const Real _Dc; + const Real _fc; + const Real _FEc; + + const Real _at; + const Real _ac; + const Real _zt; + const Real _zc; + const Real _dPhit; + const Real _dPhic; + const Real _sqrtPhit_max; + const Real _sqrtPhic_max; + const Real _dt_bt; + const Real _dc_bc; + const Real _ft0; + const Real _fc0; + const Real _small_smoother2; + + const Real _sqrt3; + + /// Whether to provide an estimate of the returned stress, based on perfect plasticity + const bool _perfect_guess; + + /// Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to stress space + RankTwoTensor _eigvecs; + + MaterialProperty & _max_principal; + MaterialProperty & _min_principal; + MaterialProperty & _intnl0; + MaterialProperty & _intnl1; + MaterialProperty & _ele_len; + MaterialProperty & _gt; + MaterialProperty & _gc; + + MaterialProperty & _tD; + MaterialProperty & _cD; + MaterialProperty & _D; + MaterialProperty & _min_ep; + MaterialProperty & _mid_ep; + MaterialProperty & _max_ep; + MaterialProperty & _sigma0; + MaterialProperty & _sigma1; + MaterialProperty & _sigma2; + + Real ft(const std::vector & intnl) const; /// tensile strength + Real dft(const std::vector & intnl) const; /// d(tensile strength)/d(intnl) + Real fc(const std::vector & intnl) const; /// compressive strength + Real dfc(const std::vector & intnl) const; /// d(compressive strength)/d(intnl) + Real beta(const std::vector & intnl) const; + Real dbeta0(const std::vector & intnl) const; + Real dbeta1(const std::vector & intnl) const; + void weighfac(const std::vector & stress_params, Real & wf) const; /// weight factor + void dweighfac(const std::vector & stress_params, + Real & wf, + std::vector & dwf) const; /// d(weight factor)/d(stress) + Real damageVar(const std::vector & stress_params, const std::vector & intnl) const; + + void computeStressParams(const RankTwoTensor & stress, + std::vector & stress_params) const override; + + std::vector dstress_param_dstress(const RankTwoTensor & stress) const override; + + std::vector d2stress_param_dstress(const RankTwoTensor & stress) const override; + + void setEffectiveElasticity(const RankFourTensor & Eijkl) override; + + virtual void preReturnMapV(const std::vector & trial_stress_params, + const RankTwoTensor & stress_trial, + const std::vector & intnl_old, + const std::vector & yf, + const RankFourTensor & Eijkl) override; + + virtual void setStressAfterReturnV(const RankTwoTensor & stress_trial, + const std::vector & stress_params, + Real gaE, + const std::vector & intnl, + const yieldAndFlow & smoothed_q, + const RankFourTensor & Eijkl, + RankTwoTensor & stress) const override; + + void yieldFunctionValuesV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & yf) const override; + + void computeAllQV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & all_q) const override; + + virtual void flowPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & r) const; + + virtual void dflowPotential_dstress(const std::vector & stress_params, + const std::vector & intnl, + std::vector> & dr_dstress) const; + + virtual void dflowPotential_dintnl(const std::vector & stress_params, + const std::vector & intnl, + std::vector> & dr_dintnl) const; + + virtual void hardPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & h) const; + + virtual void dhardPotential_dstress(const std::vector & stress_params, + const std::vector & intnl, + std::vector> & dh_dsig) const; + + virtual void dhardPotential_dintnl(const std::vector & stress_params, + const std::vector & intnl, + std::vector> & dh_dintnl) const; + + void initialiseVarsV(const std::vector & trial_stress_params, + const std::vector & intnl_old, + std::vector & stress_params, + Real & gaE, + std::vector & intnl) const; + + void setIntnlValuesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl_old, + std::vector & intnl) const override; + + void setIntnlDerivativesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl, + std::vector> & dintnl) const override; + + virtual void consistentTangentOperatorV(const RankTwoTensor & stress_trial, + const std::vector & trial_stress_params, + const RankTwoTensor & stress, + const std::vector & stress_params, + Real gaE, + const yieldAndFlow & smoothed_q, + const RankFourTensor & Eijkl, + bool compute_full_tangent_operator, + const std::vector> & dvar_dtrial, + RankFourTensor & cto) override; +}; diff --git a/src/materials/ComputeMultipleInelasticDamageStress.C b/src/materials/ComputeMultipleInelasticDamageStress.C new file mode 100644 index 00000000..4f92615e --- /dev/null +++ b/src/materials/ComputeMultipleInelasticDamageStress.C @@ -0,0 +1,73 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* BlackBear */ +/* */ +/* (c) 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +#include "ComputeMultipleInelasticDamageStress.h" +#include "StressUpdateBase.h" + +registerMooseObject("BlackBearApp", ComputeMultipleInelasticDamageStress); + +InputParameters +ComputeMultipleInelasticDamageStress::validParams() +{ + InputParameters params = ComputeMultipleInelasticStress::validParams(); + return params; +} + +ComputeMultipleInelasticDamageStress::ComputeMultipleInelasticDamageStress( + const InputParameters & parameters) + : ComputeMultipleInelasticStress(parameters), + _D(getMaterialProperty("elemental_damage_variable")), + _D_old(getMaterialPropertyOld("elemental_damage_variable")), + _D_older(getMaterialPropertyOlder("elemental_damage_variable")) +{ +} + +void +ComputeMultipleInelasticDamageStress::computeQpJacobianMult() +{ + ComputeMultipleInelasticStress::computeQpJacobianMult(); + _Jacobian_mult[_qp] = (1.0 - _D_older[_qp]) * _Jacobian_mult[_qp]; + // _Jacobian_mult[_qp] = (1.0 - _D[_qp]) * _Jacobian_mult[_qp]; +} + +void +ComputeMultipleInelasticDamageStress::updateQpStateSingleModel( + unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & combined_inelastic_strain_increment) +{ + ComputeMultipleInelasticStress::updateQpStateSingleModel( + model_number, elastic_strain_increment, combined_inelastic_strain_increment); + _Jacobian_mult[_qp] = (1.0 - _D_older[_qp]) * _Jacobian_mult[_qp]; +} + +void +ComputeMultipleInelasticDamageStress::computeAdmissibleState( + unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & inelastic_strain_increment, + RankFourTensor & consistent_tangent_operator) +{ + _models[model_number]->updateState(elastic_strain_increment, + inelastic_strain_increment, + _rotation_increment[_qp], + _stress[_qp], + _stress_old[_qp] / (1.0 - _D_older[_qp]), + // _stress_old[_qp] / (1.0 - _D[_qp]), + _elasticity_tensor[_qp], + _elastic_strain_old[_qp], + _tangent_operator_type == TangentOperatorEnum::nonlinear, + consistent_tangent_operator); + _stress[_qp] *= (1.0 - _D_older[_qp]); +} diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C new file mode 100644 index 00000000..2615573d --- /dev/null +++ b/src/materials/DamagePlasticityStressUpdate.C @@ -0,0 +1,678 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* BlackBear */ +/* */ +/* (c) 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ +#include "DamagePlasticityStressUpdate.h" +#include "libmesh/utility.h" + +registerMooseObject("BlackBearApp", DamagePlasticityStressUpdate); + +InputParameters +DamagePlasticityStressUpdate::validParams() +{ + InputParameters params = MultiParameterPlasticityStressUpdate::validParams(); + params.addParam( + "yield_function_tolerance", + "If the yield function is less than this amount, the (stress, internal parameters) are " + "deemed admissible. A std::vector of tolerances must be entered for the multi-surface case"); + params.addRangeCheckedParam("factor_relating_biaxial_unixial_cmp_str", + 0.1, + "factor_relating_biaxial_unixial_cmp_str < 0.5 & " + "factor_relating_biaxial_unixial_cmp_str >= 0", + "Material parameter that relate biaxial and uniaxial " + "compressive strength, i.e., \alfa = (fb0-fc0)/(2*fb0-fc0)"); + params.addRequiredParam("factor_controlling_dilatancy", + "controls the dilation of concrete"); + params.addRangeCheckedParam("stiff_recovery_factor", + 0., + "stiff_recovery_factor <= 1. & stiff_recovery_factor >= 0", + "stiffness recovery factor"); + params.addRangeCheckedParam( + "ft_ep_slope_factor_at_zero_ep", + "ft_ep_slope_factor_at_zero_ep <= 1 & ft_ep_slope_factor_at_zero_ep >= 0", + "slope of ft vs plastic strain curve at zero plastic strain"); + params.addRequiredParam( + "tensile_damage_at_half_tensile_strength", + "fraction of the elastic recovery slope in tension at 0.5*ft0 after yielding"); + params.addRangeCheckedParam("yield_strength_in_tension", + "yield_strength_in_tension >= 0", + "Tensile yield strength of concrete"); + params.addRangeCheckedParam("fracture_energy_in_tension", + "fracture_energy_in_tension >= 0", + "Fracture energy of concrete in uniaxial tension"); + + params.addRangeCheckedParam("yield_strength_in_compression", + "yield_strength_in_compression >= 0", + "Absolute yield compressice strength"); + params.addRequiredParam("compressive_damage_at_max_compressive_strength", + "damage at maximum compressive strength"); + params.addRequiredParam("maximum_strength_in_compression", + "Absolute maximum compressive strength"); + params.addRangeCheckedParam("fracture_energy_in_compression", + "fracture_energy_in_compression >= 0", + "Fracture energy of concrete in uniaxial compression"); + + params.addRequiredRangeCheckedParam( + "tip_smoother", + "tip_smoother>=0", + "Smoothing parameter: the cone vertex at mean = cohesion*cot(friction_angle), will be " + "smoothed by the given amount. Typical value is 0.1*cohesion"); + params.addParam("perfect_guess", + true, + "Provide a guess to the Newton-Raphson proceedure " + "that is the result from perfect plasticity. With " + "severe hardening/softening this may be " + "suboptimal."); + params.addClassDescription("Damage Plasticity Model for concrete"); + return params; +} + +DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters & parameters) + : MultiParameterPlasticityStressUpdate(parameters, 3, 1, 2), + _f_tol(getParam("yield_function_tol")), + + _alfa(getParam("factor_relating_biaxial_unixial_cmp_str")), + _alfa_p(getParam("factor_controlling_dilatancy")), + _s0(getParam("stiff_recovery_factor")), + + _Chi(getParam("ft_ep_slope_factor_at_zero_ep")), + _Dt(getParam("tensile_damage_at_half_tensile_strength")), + _ft(getParam("yield_strength_in_tension")), + _FEt(getParam("fracture_energy_in_tension")), + + _fyc(getParam("yield_strength_in_compression")), + _Dc(getParam("compressive_damage_at_max_compressive_strength")), + _fc(getParam("maximum_strength_in_compression")), + _FEc(getParam("fracture_energy_in_compression")), + + _at(1.5 * std::sqrt(1 - _Chi) - 0.5), + _ac((2. * (_fc / _fyc) - 1. + 2. * std::sqrt(std::pow((_fc / _fyc), 2.) - _fc / _fyc))), + + _zt((1. + _at) / _at), + _zc((1. + _ac) / _ac), + _dPhit(_at * (2. + _at)), + _dPhic(_ac * (2. + _ac)), + _sqrtPhit_max((1. + _at + sqrt(1. + _at * _at)) / 2.), + _sqrtPhic_max((1. + _ac) / 2.), + _dt_bt(log(1. - _Dt) / log((1. + _at - sqrt(1. + _at * _at)) / (2. * _at))), + _dc_bc(log(1. - _Dc) / log((1. + _ac) / (2. * _ac))), + _ft0(0.5 * _ft / + ((1. - _Dt) * pow((_zt - _sqrtPhit_max / _at), (1. - _dt_bt)) * _sqrtPhit_max)), + _fc0(_fc / ((1. - _Dc) * pow((_zc - _sqrtPhic_max / _ac), (1. - _dc_bc)) * _sqrtPhic_max)), + _small_smoother2(std::pow(getParam("tip_smoother"), 2)), + + _sqrt3(sqrt(3.)), + _perfect_guess(getParam("perfect_guess")), + _eigvecs(RankTwoTensor()), + _max_principal(declareProperty("max_principal_stress")), + _min_principal(declareProperty("min_principal_stress")), + _intnl0(declareProperty("damage_state_in_tension")), + _intnl1(declareProperty("damage_state_in_compression")), + _ele_len(declareProperty("element_length")), + _gt(declareProperty("elemental_fracture_energy_in_tension")), + _gc(declareProperty("elemental_fracture_energy_in_compression")), + _tD(declareProperty("elemental_tensile_damage")), + _cD(declareProperty("elemental_compression_damage")), + _D(declareProperty("elemental_damage_variable")), + _min_ep(declareProperty("min_ep")), + _mid_ep(declareProperty("mid_ep")), + _max_ep(declareProperty("max_ep")), + _sigma0(declareProperty("damaged_min_principal_stress")), + _sigma1(declareProperty("damaged_mid_principal_stress")), + _sigma2(declareProperty("damaged_max_principal_stress")) +{ +} + +void +DamagePlasticityStressUpdate::initQpStatefulProperties() +{ + // if (_current_elem->n_vertices() < 3) + // _ele_len[_qp] = _current_elem->length(0, 1); + // else if (_current_elem->n_vertices() < 5) + // _ele_len[_qp] = (_current_elem->length(0, 1) + _current_elem->length(1, 2)) / 2.; + // else + // _ele_len[_qp] = + // (_current_elem->length(0, 1) + _current_elem->length(1, 2) + _current_elem->length(0, 4)) + // / 3.; + _ele_len[_qp] = std::cbrt(_current_elem->volume()); + + _gt[_qp] = _FEt / _ele_len[_qp]; + _gc[_qp] = _FEc / _ele_len[_qp]; + + _min_ep[_qp] = 0.; + _mid_ep[_qp] = 0.; + _max_ep[_qp] = 0.; + _sigma0[_qp] = 0.; + _sigma1[_qp] = 0.; + _sigma2[_qp] = 0.; + _intnl0[_qp] = 0.; + _intnl1[_qp] = 0.; + _tD[_qp] = 0.; + _cD[_qp] = 0.; + _D[_qp] = 0.; + MultiParameterPlasticityStressUpdate::initQpStatefulProperties(); +} + +void +DamagePlasticityStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/) +{ + std::vector eigstrain; + _plastic_strain[_qp].symmetricEigenvalues(eigstrain); + _min_ep[_qp] = eigstrain[0]; + _mid_ep[_qp] = eigstrain[1]; + _max_ep[_qp] = eigstrain[2]; +} + +void +DamagePlasticityStressUpdate::computeStressParams(const RankTwoTensor & stress, + std::vector & stress_params) const +{ + stress.symmetricEigenvalues(stress_params); +} + +std::vector +DamagePlasticityStressUpdate::dstress_param_dstress(const RankTwoTensor & stress) const +{ + std::vector sp; + std::vector dsp; + stress.dsymmetricEigenvalues(sp, dsp); + return dsp; +} + +std::vector +DamagePlasticityStressUpdate::d2stress_param_dstress(const RankTwoTensor & stress) const +{ + std::vector d2; + stress.d2symmetricEigenvalues(d2); + return d2; +} + +void +DamagePlasticityStressUpdate::setEffectiveElasticity(const RankFourTensor & Eijkl) +{ + // Eijkl is required to be isotropic, so we can use the + // frame where stress is diagonal + for (unsigned a = 0; a < _num_sp; ++a) + for (unsigned b = 0; b < _num_sp; ++b) + _Eij[a][b] = Eijkl(a, a, b, b); + _En = _Eij[2][2]; + const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]); + for (unsigned a = 0; a < _num_sp; ++a) + { + _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom; + for (unsigned b = 0; b < a; ++b) + _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom; + } +} + +void +DamagePlasticityStressUpdate::preReturnMapV(const std::vector & /*trial_stress_params*/, + const RankTwoTensor & stress_trial, + const std::vector & /*intnl_old*/, + const std::vector & /*yf*/, + const RankFourTensor & /*Eijkl*/) +{ + std::vector eigvals; + stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs); +} + +void +DamagePlasticityStressUpdate::setStressAfterReturnV(const RankTwoTensor & /*stress_trial*/, + const std::vector & stress_params, + Real /*gaE*/, + const std::vector & intnl, + const yieldAndFlow & /*smoothed_q*/, + const RankFourTensor & /*Eijkl*/, + RankTwoTensor & stress) const +{ + // form the diagonal stress + stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0); + // rotate to the original frame + stress = _eigvecs * stress * (_eigvecs.transpose()); + // _dir[_qp] = _eigvecs; + Real D = damageVar(stress_params, intnl); + _sigma0[_qp] = (1. - D) * stress_params[0]; + _sigma1[_qp] = (1. - D) * stress_params[1]; + _sigma2[_qp] = (1. - D) * stress_params[2]; + _intnl0[_qp] = intnl[0]; + _intnl1[_qp] = intnl[1]; + _D[_qp] = D; +} + +void +DamagePlasticityStressUpdate::yieldFunctionValuesV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & yf) const +{ + Real I1 = stress_params[0] + stress_params[1] + stress_params[2]; + Real J2 = + (pow(stress_params[0] - stress_params[1], 2.) + pow(stress_params[1] - stress_params[2], 2.) + + pow(stress_params[2] - stress_params[0], 2.)) / + 6. + + _small_smoother2; + Real sqrtJ2 = sqrt(J2); + yf[0] = 1. / (1. - _alfa) * + (_alfa * I1 + _sqrt3 * sqrtJ2 + + beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - + fc(intnl); +} + +void +DamagePlasticityStressUpdate::computeAllQV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & all_q) const +{ + Real I1 = stress_params[0] + stress_params[1] + stress_params[2]; + Real J2 = + (pow(stress_params[0] - stress_params[1], 2.) + pow(stress_params[1] - stress_params[2], 2.) + + pow(stress_params[2] - stress_params[0], 2.)) / + 6. + + _small_smoother2; + Real sqrtJ2 = sqrt(J2); + std::vector DevSt(3); // vector of principal deviatoric stress + for (unsigned i = 0; i < 3; ++i) + DevSt[i] = stress_params[i] - I1 / 3.; + + // yieldFunctionValuesV(stress_params, intnl, all_q[0].f); + all_q[0].f = 1. / (1. - _alfa) * + (_alfa * I1 + _sqrt3 * sqrtJ2 + + beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - + fc(intnl); + + for (unsigned i = 0; i < _num_sp; ++i) + all_q[0].df[i] = 1. / (1. - _alfa) * + (_alfa + _sqrt3 * DevSt[i] / (2. * sqrtJ2) + + beta(intnl) * (stress_params[2] < 0. ? 0. : (i == 2))); + all_q[0].df_di[0] = + 1. / (1. - _alfa) * (dbeta0(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])); + all_q[0].df_di[1] = + 1. / (1. - _alfa) * (dbeta1(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) - + dfc(intnl); + + flowPotential(stress_params, intnl, all_q[0].dg); + dflowPotential_dstress(stress_params, intnl, all_q[0].d2g); + dflowPotential_dintnl(stress_params, intnl, all_q[0].d2g_di); +} + +void +DamagePlasticityStressUpdate::flowPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & r) const +{ + Real J2 = + (pow(stress_params[0] - stress_params[1], 2.) + pow(stress_params[1] - stress_params[2], 2.) + + pow(stress_params[2] - stress_params[0], 2.)) / + 6. + + _small_smoother2; + Real invsqrt2J2 = 1. / sqrt(2. * J2); + std::vector DevSt(3); + DevSt[0] = (2. * stress_params[0] - stress_params[1] - stress_params[2]) / 3.; // dJ2/dsig0 + DevSt[1] = (2. * stress_params[1] - stress_params[2] - stress_params[0]) / 3.; // dJ2/dsig1 + DevSt[2] = (2. * stress_params[2] - stress_params[0] - stress_params[1]) / 3.; // dJ2/dsig2 + + Real D = damageVar(stress_params, intnl); + + for (unsigned int i = 0; i < _num_sp; ++i) + r[i] = (_alfa_p + (J2 < _f_tol ? 0. : DevSt[i] * invsqrt2J2)) * pow((1. - D), 1); +} + +void +DamagePlasticityStressUpdate::dflowPotential_dstress( + const std::vector & stress_params, + const std::vector & intnl, + std::vector> & dr_dstress) const +{ + Real J2 = + (pow(stress_params[0] - stress_params[1], 2.) + pow(stress_params[1] - stress_params[2], 2.) + + pow(stress_params[2] - stress_params[0], 2.)) / + 6. + + _small_smoother2; + Real invsqrt2J2 = 1. / sqrt(2. * J2); + std::vector DevSt(3); + DevSt[0] = (2. * stress_params[0] - stress_params[1] - stress_params[2]) / 3.; // dJ2/dsig0 + DevSt[1] = (2. * stress_params[1] - stress_params[2] - stress_params[0]) / 3.; // dJ2/dsig1 + DevSt[2] = (2. * stress_params[2] - stress_params[0] - stress_params[1]) / 3.; // dJ2/dsig2 + + Real D = damageVar(stress_params, intnl); + + for (unsigned i = 0; i < _num_sp; ++i) + for (unsigned j = 0; j < (i + 1); ++j) + { + if (i != j) + { + dr_dstress[i][j] = J2 < _f_tol ? 0. + : invsqrt2J2 * (-1. / 3. - DevSt[i] * DevSt[j] / (2. * J2)) * + pow((1. - D), 2); + dr_dstress[j][i] = dr_dstress[i][j]; + } + else + dr_dstress[i][i] = J2 < _f_tol ? 0. + : invsqrt2J2 * (2. / 3. - DevSt[i] * DevSt[j] / (2. * J2)) * + pow((1. - D), 2); + } +} + +void +DamagePlasticityStressUpdate::dflowPotential_dintnl( + const std::vector & /* stress_params */, + const std::vector & /* intnl */, + std::vector> & dr_dintnl) const +{ + for (unsigned i = 0; i < _num_sp; ++i) + for (unsigned j = 0; j < _num_intnl; ++j) + dr_dintnl[i][j] = 0.; +} + +void +DamagePlasticityStressUpdate::hardPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & h) const +{ + Real wf; + weighfac(stress_params, wf); + std::vector r(3); + flowPotential(stress_params, intnl, r); + h[0] = wf * ft(intnl) / _gt[_qp] * r[2]; + h[1] = -(1. - wf) * fc(intnl) / _gc[_qp] * r[0]; +} + +void +DamagePlasticityStressUpdate::dhardPotential_dstress(const std::vector & stress_params, + const std::vector & intnl, + std::vector> & dh_dsig) const +{ + Real wf; + std::vector dwf(3); + dweighfac(stress_params, wf, dwf); + + std::vector r(3); + flowPotential(stress_params, intnl, r); + std::vector> dr_dsig(3, std::vector(3)); + dflowPotential_dstress(stress_params, intnl, dr_dsig); + + for (unsigned i = 0; i < _num_sp; ++i) + { + dh_dsig[0][i] = (wf * dr_dsig[2][i] + dwf[i] * r[2]) * ft(intnl) / _gt[_qp]; + dh_dsig[1][i] = -((1. - wf) * dr_dsig[0][i] - dwf[i] * r[0]) * fc(intnl) / _gc[_qp]; + } +} + +void +DamagePlasticityStressUpdate::dhardPotential_dintnl( + const std::vector & stress_params, + const std::vector & intnl, + std::vector> & dh_dintnl) const +{ + Real wf; + weighfac(stress_params, wf); + std::vector r(3); + flowPotential(stress_params, intnl, r); + + dh_dintnl[0][0] = wf * dft(intnl) / _gt[_qp] * r[2]; + dh_dintnl[0][1] = 0.; + dh_dintnl[1][0] = 0.; + dh_dintnl[1][1] = -(1 - wf) * dfc(intnl) / _gc[_qp] * r[0]; +} + +void +DamagePlasticityStressUpdate::initialiseVarsV(const std::vector & trial_stress_params, + const std::vector & intnl_old, + std::vector & stress_params, + Real & /* gaE */, + std::vector & intnl) const +{ + setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl); +} + +void +DamagePlasticityStressUpdate::setIntnlValuesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl_old, + std::vector & intnl) const +{ + Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2]; + Real J2_trial = (pow(trial_stress_params[0] - trial_stress_params[1], 2.) + + pow(trial_stress_params[1] - trial_stress_params[2], 2.) + + pow(trial_stress_params[2] - trial_stress_params[0], 2.)) / + 6. + + _small_smoother2; + Real invsqrt2J2_trial = 1. / sqrt(2. * J2_trial); + Real G = 0.5 * (_Eij[0][0] - _Eij[0][1]); // Lame's mu + Real K = _Eij[0][1] + 2. * G / 3.; // Bulk modulus + Real C1 = (2. * G * invsqrt2J2_trial); + Real C2 = -(I1_trial / 3. * G * invsqrt2J2_trial - 3. * K * _alfa_p); + Real C3 = 3. * K * _alfa_p; + + RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0], + trial_stress_params[1] - current_stress_params[1], + trial_stress_params[2] - current_stress_params[2], + 0., + 0., + 0.); + RankTwoTensor fac = J2_trial < _f_tol ? C3 * RankTwoTensor(1., 1., 1., 0., 0., 0.) + : RankTwoTensor(C1 * trial_stress_params[0] - C2, + C1 * trial_stress_params[1] - C2, + C1 * trial_stress_params[2] - C2, + 0., + 0., + 0.); + + Real lam = dsig.L2norm() / fac.L2norm(); + std::vector h(2); + hardPotential(current_stress_params, intnl_old, h); + + intnl[0] = intnl_old[0] + lam * h[0]; + intnl[1] = intnl_old[1] + lam * h[1]; +} + +void +DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl, + std::vector> & dintnl) const +{ + Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2]; + Real J2_trial = (pow(trial_stress_params[0] - trial_stress_params[1], 2.) + + pow(trial_stress_params[1] - trial_stress_params[2], 2.) + + pow(trial_stress_params[2] - trial_stress_params[0], 2.)) / + 6.; + Real invsqrt2J2_trial = 1. / sqrt(2. * J2_trial); + Real G = 0.5 * (_Eij[0][0] - _Eij[0][1]); // Lame's mu + Real K = _Eij[0][1] + 2. * G / 3.; // Bulk modulus + Real C1 = (2. * G * invsqrt2J2_trial); + Real C2 = -(I1_trial / 3. * G * invsqrt2J2_trial - 3. * K * _alfa_p); + Real C3 = 3. * K * _alfa_p; + + RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0], + trial_stress_params[1] - current_stress_params[1], + trial_stress_params[2] - current_stress_params[2], + 0., + 0., + 0.); + RankTwoTensor fac = J2_trial < _f_tol ? C3 * RankTwoTensor(1., 1., 1., 0., 0., 0.) + : RankTwoTensor(C1 * trial_stress_params[0] - C2, + C1 * trial_stress_params[1] - C2, + C1 * trial_stress_params[2] - C2, + 0., + 0., + 0.); + + Real lam = dsig.L2norm() / fac.L2norm(); + + std::vector dlam_dsig(3); + for (unsigned i = 0; i < _num_sp; ++i) + dlam_dsig[i] = dsig.L2norm() == 0. ? 0. + : -(trial_stress_params[i] - current_stress_params[i]) / + (dsig.L2norm() * fac.L2norm()); + + std::vector h(2); + hardPotential(current_stress_params, intnl, h); + std::vector> dh_dsig(2, std::vector(3)); + dhardPotential_dstress(current_stress_params, intnl, dh_dsig); + std::vector> dh_dintnl(2, std::vector(2)); + dhardPotential_dintnl(current_stress_params, intnl, dh_dintnl); + + for (unsigned i = 0; i < _num_intnl; ++i) + for (unsigned j = 0; j < _num_sp; ++j) + dintnl[i][j] = dlam_dsig[j] * h[i] + lam * dh_dsig[i][j]; +} + +Real +DamagePlasticityStressUpdate::ft(const std::vector & intnl) const +{ + Real sqrtPhi_t = sqrt(1. + _at * (2. + _at) * intnl[0]); + if (_zt > sqrtPhi_t / _at) + return _ft0 * pow(_zt - sqrtPhi_t / _at, (1. - _dt_bt)) * sqrtPhi_t; + else + return _ft0 * 1.E-6; +} + +Real +DamagePlasticityStressUpdate::dft(const std::vector & intnl) const +{ + Real sqrtPhi_t = sqrt(1. + _at * (2. + _at) * intnl[0]); + if (_zt > sqrtPhi_t / _at) + return _ft0 * _dPhit / (2 * sqrtPhi_t) * pow(_zt - sqrtPhi_t / _at, -_dt_bt) * + (_zt - (2. - _dt_bt) * sqrtPhi_t / _at); + else + return 0.; +} + +Real +DamagePlasticityStressUpdate::fc(const std::vector & intnl) const +{ + Real sqrtPhi_c; + if (intnl[1] < 1.0) + sqrtPhi_c = sqrt(1. + _ac * (2. + _ac) * intnl[1]); + else + sqrtPhi_c = sqrt(1. + _ac * (2. + _ac) * 0.99); + return _fc0 * pow((_zc - sqrtPhi_c / _ac), (1. - _dc_bc)) * sqrtPhi_c; +} + +Real +DamagePlasticityStressUpdate::dfc(const std::vector & intnl) const +{ + if (intnl[1] < 1.0) + { + Real sqrtPhi_c = sqrt(1. + _ac * (2. + _ac) * intnl[1]); + return _fc0 * _dPhic / (2. * sqrtPhi_c) * pow(_zc - sqrtPhi_c / _ac, -_dc_bc) * + (_zc - (2. - _dc_bc) * sqrtPhi_c / _ac); + } + else + return 0.; +} + +Real +DamagePlasticityStressUpdate::beta(const std::vector & intnl) const +{ + return (1. - _alfa) * fc(intnl) / ft(intnl) - (1. + _alfa); +} + +Real +DamagePlasticityStressUpdate::dbeta0(const std::vector & intnl) const +{ + return -(1. - _alfa) * fc(intnl) * dft(intnl) / pow(ft(intnl), 2.); +} + +Real +DamagePlasticityStressUpdate::dbeta1(const std::vector & intnl) const +{ + return dfc(intnl) / ft(intnl) * (1. - _alfa); +} + +void +DamagePlasticityStressUpdate::weighfac(const std::vector & stress_params, Real & wf) const +{ + Real Dr = 0.; + Real Nr = 0.; + for (unsigned i = 0; i < _num_sp; ++i) + { + if (stress_params[i] > 0.) + { + Nr += stress_params[i]; + Dr += stress_params[i]; + } + else + Dr += -stress_params[i]; + } + wf = Nr / Dr; +} + +void +DamagePlasticityStressUpdate::dweighfac(const std::vector & stress_params, + Real & wf, + std::vector & dwf) const +{ + std::vector dNr(3, 0.), dDr(3, 0.); + Real Dr = 0.; + Real Nr = 0.; + for (unsigned i = 0; i < _num_sp; ++i) + { + if (stress_params[i] > 0.) + { + Nr += stress_params[i]; + dNr[i] = 1.; + Dr += stress_params[i]; + dDr[i] = 1.; + } + else + { + Dr += -stress_params[i]; + dDr[i] = -1.; + } + } + wf = Nr / Dr; + + for (unsigned i = 0; i < _num_sp; ++i) + dwf[i] = (dNr[i] - wf * dDr[i]) / Dr; +} + +Real +DamagePlasticityStressUpdate::damageVar(const std::vector & stress_params, + const std::vector & intnl) const +{ + Real sqrtPhi_t = sqrt(1. + _at * (2. + _at) * intnl[0]); + if (_zt > sqrtPhi_t / _at) + _tD[_qp] = 1. - pow(_zt - sqrtPhi_t / _at, _dt_bt); + else + _tD[_qp] = 1. - 1.E-6; + + Real wf; + weighfac(stress_params, wf); + Real s = _s0 + (1. - _s0) * wf; + + Real sqrtPhi_c; + if (intnl[1] < 1.0) + sqrtPhi_c = sqrt(1. + _ac * (2. + _ac) * intnl[1]); + else + sqrtPhi_c = sqrt(1. + _ac * (2. + _ac) * 0.99); + + _cD[_qp] = 1. - pow((_zc - sqrtPhi_c / _ac), _dc_bc); + return 1. - (1. - s * _tD[_qp]) * (1. - _cD[_qp]); +} + +void +DamagePlasticityStressUpdate::consistentTangentOperatorV( + const RankTwoTensor & /* stress_trial */, + const std::vector & /* trial_stress_params */, + const RankTwoTensor & /*stress*/, + const std::vector & /* stress_params */, + Real /*gaE*/, + const yieldAndFlow & /*smoothed_q*/, + const RankFourTensor & elasticity_tensor, + bool /* compute_full_tangent_operator */, + const std::vector> & /* dvar_dtrial */, + RankFourTensor & cto) +{ + cto = elasticity_tensor; + return; +} diff --git a/test/tests/damage_plasticity_model/gold/uniaxial_compression_out.csv b/test/tests/damage_plasticity_model/gold/uniaxial_compression_out.csv new file mode 100644 index 00000000..fc1cc3f2 --- /dev/null +++ b/test/tests/damage_plasticity_model/gold/uniaxial_compression_out.csv @@ -0,0 +1,52 @@ +time,D,displacement_x,e_xx,e_xy,e_xz,e_yy,e_yz,e_zz,intnl0,intnl1,max_ep,mid_ep,min_ep,react_x,react_y,react_z,s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,sigma0,sigma1,sigma2 +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 +10,0,-0.0001,-0.0001,-1.99E-21,-8.34E-22,1.80E-05,-2.26E-21,1.80E-05,0,0,0,0,0,3.17,5.33E-17,2.24E-17,-3.17,-5.33E-17,-2.24E-17,3.95E-13,-6.06E-17,3.95E-13,0,0,0 +20,0,-0.0002,-0.0002,-2.20E-21,-3.44E-22,3.60E-05,-3.27E-21,3.60E-05,0,0,0,0,0,6.34,5.90E-17,9.24E-18,-6.34,-5.90E-17,-9.24E-18,-4.11E-13,-8.78E-17,-4.11E-13,0,0,0 +30,0,-0.0003,-0.0003,4.76E-22,6.35E-22,5.40E-05,2.67E-21,5.40E-05,0,0,0,0,0,9.51,-1.28E-17,-1.71E-17,-9.51,1.28E-17,1.71E-17,5.05E-15,7.18E-17,5.02E-15,0,0,0 +40,0,-0.0004,-0.0004,-1.27E-21,-1.11E-21,7.20E-05,-9.26E-22,7.20E-05,0,0,0,0,0,12.68,3.41E-17,2.99E-17,-12.68,-3.41E-17,-2.99E-17,1.54E-12,-2.49E-17,1.54E-12,0,0,0 +50,0,-0.0005,-0.0005,-1.06E-21,-8.58E-21,9.00E-05,-3.71E-22,9.00E-05,0,0,0,0,0,15.85,2.84E-17,2.30E-16,-15.85,-2.84E-17,-2.30E-16,7.03E-13,-9.96E-18,7.04E-13,0,0,0 +60,0.00350959,-0.0006,-0.000588195,-2.87E-20,-8.49E-21,0.000105875,-4.68E-21,0.000105875,4.16E-17,0.001954762,1.28E-05,1.28E-05,-1.18E-05,18.64578867,3.07E-16,2.24E-16,-18.64578867,-3.07E-16,-2.24E-16,1.72E-13,-1.63E-16,1.71E-13,-18.5803496,1.70E-13,1.72E-13 +70,0.019535827,-0.0007,-0.000637024,-1.30E-20,-2.46E-21,0.000114664,-1.86E-20,0.000114664,1.54E-13,0.011340486,6.85E-05,6.85E-05,-6.30E-05,20.19367563,1.83E-15,1.19E-16,-20.19367563,-1.83E-15,-1.19E-16,1.45E-10,-4.26E-16,1.45E-10,-19.79917548,1.42E-10,1.42E-10 +80,0.035679756,-0.0008,-0.000684583,4.46E-20,1.27E-20,0.000123225,2.98E-20,0.000123225,2.56E-13,0.0215388,0.000125601,0.000125601,-0.000115417,21.62512169,-2.02E-15,-1.82E-16,-21.62512169,2.02E-15,1.82E-16,1.03E-10,7.64E-16,1.03E-10,-20.92698777,9.97E-11,9.97E-11 +90,0.051941449,-0.0009,-0.000730836,2.89E-21,-2.63E-21,0.00013155,1.00E-20,0.00013155,3.26E-13,0.032540569,0.000184091,0.000184091,-0.000169164,22.71489264,-8.84E-16,2.51E-17,-22.71489264,8.84E-16,-2.51E-17,7.41E-11,3.13E-16,7.41E-11,-21.96413576,7.17E-11,7.17E-11 +100,0.068320452,-0.001,-0.000775745,-1.04E-20,2.71E-22,0.000139634,4.89E-20,0.000139634,3.75E-13,0.044335289,0.000244043,0.000244043,-0.000224255,23.71370544,-1.16E-15,3.01E-17,-23.71370544,1.16E-15,-3.01E-17,5.42E-11,1.23E-15,5.42E-11,-22.91103448,5.24E-11,5.24E-11 +110,0.084815783,-0.0011,-0.000819272,4.20E-20,1.50E-20,0.000147469,-5.31E-20,0.000147469,4.11E-13,0.056910948,0.000305499,0.000305499,-0.000280728,24.62195323,-2.21E-15,-3.81E-16,-24.62195323,2.21E-15,3.81E-16,4.01E-11,-1.26E-15,4.01E-11,-23.7681765,3.87E-11,3.87E-11 +120,0.101425895,-0.0012,-0.000861377,-5.31E-21,6.75E-21,0.000155048,2.28E-19,0.000155048,4.36E-13,0.070253874,0.000368503,0.000368503,-0.000338623,25.44011227,-1.08E-16,5.35E-17,-25.44011227,1.08E-16,-5.35E-17,2.99E-11,5.52E-15,2.99E-11,-24.53614674,2.89E-11,2.89E-11 +130,0.118148636,-0.0013,-0.000902018,-3.78E-21,1.05E-20,0.000162363,-6.76E-19,0.000162363,4.55E-13,0.084348583,0.000433099,0.000433099,-0.000397982,26.16875604,-2.58E-17,-2.87E-16,-26.16875604,2.58E-17,2.87E-16,2.25E-11,-1.66E-14,2.25E-11,-25.21563723,2.17E-11,2.17E-11 +140,0.134981205,-0.0014,-0.000941154,9.53E-21,-8.16E-21,0.000169408,1.98E-18,0.000169408,4.69E-13,0.099177623,0.000499334,0.000499334,-0.000458846,26.80857009,-4.99E-17,2.69E-16,-26.80857009,4.99E-17,-2.69E-16,1.70E-11,4.79E-14,1.70E-11,-25.80746194,1.63E-11,1.64E-11 +150,0.151920106,-0.0015,-0.00097874,3.32E-21,-2.78E-20,0.000176173,-7.51E-18,0.000176173,4.79E-13,0.114721426,0.000567256,0.000567256,-0.00052126,27.3603668,-7.12E-17,6.78E-16,-27.3603668,7.12E-17,-6.78E-16,1.28E-11,-1.78E-13,1.28E-11,-26.31257139,1.21E-11,1.25E-11 +160,0.168961097,-0.0016,-0.001014733,3.92E-20,6.80E-20,0.000182652,2.61E-17,0.000182652,4.87E-13,0.130958154,0.00063691,0.00063691,-0.000585267,27.82510011,-9.94E-16,-1.58E-15,-27.82510011,9.94E-16,1.58E-15,9.63E-12,6.06E-13,9.63E-12,-26.73206731,8.63E-12,9.88E-12 +170,0.18609914,-0.0017,-0.00104909,6.96E-20,3.68E-20,0.000188836,-1.24E-16,0.000188836,4.93E-13,0.147863558,0.000708345,0.000708345,-0.00065091,28.20387957,-1.59E-15,-8.38E-16,-28.20387957,1.59E-15,8.38E-16,7.68E-12,-2.83E-12,7.68E-12,-27.06721621,4.44E-12,1.03E-11 +180,0.203328288,-0.0018,-0.001081767,2.40E-19,2.09E-19,0.000194718,5.31E-16,0.000194718,5.00E-13,0.165410787,0.000781609,0.000781609,-0.000718233,28.49798676,-5.35E-15,-4.67E-15,-28.49798676,5.35E-15,4.67E-15,5.57E-12,1.19E-11,5.57E-12,-27.31946704,-6.91E-12,1.76E-11 +190,0.220641758,-0.0019,-0.001112719,6.18E-20,9.94E-20,0.000200289,5.99E-16,0.000200289,5.58E-13,0.183570384,0.000856749,0.000856749,-0.000787281,28.70888792,-1.35E-15,-2.17E-15,-28.70888792,1.35E-15,2.17E-15,7.44E-11,1.31E-11,7.44E-11,-27.49045923,5.77E-11,8.47E-11 +200,0.238031743,-0.002,-0.001141906,1.44E-20,-2.71E-20,0.000205543,7.65E-16,0.000205543,6.30E-13,0.202310046,0.000933811,0.000933811,-0.000858094,28.83824946,-3.09E-16,6.05E-16,-28.83824946,3.09E-16,-6.05E-16,9.16E-11,1.64E-11,9.16E-11,-27.58203956,7.07E-11,1.04E-10 +210,0.255489392,-0.0021,-0.001169285,-1.53E-20,1.20E-20,0.000210471,-4.34E-18,0.000210471,6.30E-13,0.221594565,0.00101284,0.00101284,-0.000930715,28.88794706,1.51E-16,-1.25E-16,-28.88794706,-1.51E-16,1.25E-16,-1.25E-08,-9.09E-14,-1.25E-08,-27.59627328,-1.20E-08,-1.20E-08 +220,0.273004767,-0.0022,-0.001194817,1.02E-20,-6.65E-20,0.000215067,-4.47E-18,0.000215067,6.30E-13,0.241385762,0.001093879,0.001093879,-0.001005183,28.86008204,-2.34E-16,1.31E-15,-28.86008204,2.34E-16,-1.31E-15,-1.32E-08,-9.14E-14,-1.32E-08,-27.53545421,-1.26E-08,-1.26E-08 +230,0.290566778,-0.0023,-0.001218466,-2.08E-20,-1.17E-20,0.000219324,-5.96E-19,0.000219324,6.30E-13,0.261642433,0.001176968,0.001176968,-0.001081534,28.75699116,7.63E-16,2.43E-16,-28.75699116,-7.63E-16,-2.43E-16,-1.35E-08,-1.19E-14,-1.35E-08,-27.40211448,-1.29E-08,-1.29E-08 +240,0.308163145,-0.0024,-0.001240196,-2.69E-20,-3.94E-20,0.000223235,-3.38E-19,0.000223235,6.30E-13,0.282320351,0.001262143,0.001262143,-0.001159804,28.58125514,2.02E-16,6.41E-16,-28.58125514,-2.02E-16,-6.41E-16,-1.34E-08,-6.53E-15,-1.34E-08,-27.1990307,-1.27E-08,-1.27E-08 +250,0.325780369,-0.0025,-0.001259979,2.00E-20,-5.67E-20,0.000226796,-5.79E-18,0.000226796,6.30E-13,0.303372292,0.001349439,0.001349439,-0.001240021,28.33570579,-4.90E-16,9.49E-16,-28.33570579,4.90E-16,-9.49E-16,-1.27E-08,-1.10E-13,-1.27E-08,-26.92922817,-1.21E-08,-1.21E-08 +260,0.343403702,-0.0026,-0.001277787,-6.07E-21,1.67E-20,0.000230002,1.29E-18,0.000230002,6.30E-13,0.324748129,0.001438884,0.001438884,-0.001322213,28.02342987,-6.59E-17,-3.27E-16,-28.02342987,6.59E-17,3.27E-16,-1.15E-08,2.40E-14,-1.15E-08,-26.59598167,-1.09E-08,-1.09E-08 +270,0.361017155,-0.0027,-0.001293598,1.68E-20,-1.86E-20,0.000232848,-6.59E-18,0.000232848,6.30E-13,0.346394971,0.001530501,0.001530501,-0.001406402,27.64776983,-3.04E-16,1.77E-16,-27.64776983,3.04E-16,-1.77E-16,-1.17E-08,-1.19E-13,-1.17E-08,-26.20281256,-1.11E-08,-1.11E-08 +280,0.378603504,-0.0028,-0.001307398,-2.43E-20,2.99E-20,0.000235332,3.12E-18,0.000235332,6.30E-13,0.368257371,0.001624307,0.001624307,-0.001492602,27.21232051,4.04E-16,-4.78E-16,-27.21232051,-4.04E-16,4.78E-16,-1.29E-08,5.51E-14,-1.29E-08,-25.7534815,-1.23E-08,-1.23E-08 +290,0.396144337,-0.0029,-0.001319177,3.85E-20,5.83E-20,0.000237452,2.99E-18,0.000237452,6.30E-13,0.390277596,0.001720313,0.001720313,-0.001580823,26.72092154,-7.96E-16,-8.02E-16,-26.72092154,7.96E-16,8.02E-16,-1.25E-08,5.12E-14,-1.25E-08,-25.25197648,-1.19E-08,-1.19E-08 +300,0.413620116,-0.003,-0.001328931,-1.17E-19,-1.24E-19,0.000239208,5.46E-17,0.000239208,6.30E-13,0.412395977,0.001818522,0.001818522,-0.001671069,26.17764477,1.89E-15,2.16E-15,-26.17764477,-1.89E-15,-2.16E-15,-1.16E-08,9.11E-13,-1.16E-08,-24.70249578,-1.09E-08,-1.09E-08 +310,0.431010268,-0.0031,-0.001336667,-5.02E-19,-9.53E-19,0.0002406,1.23E-15,0.0002406,6.30E-13,0.434551324,0.001918927,0.001918927,-0.001763333,25.58677634,8.14E-15,1.55E-14,-25.58677634,-8.14E-15,-1.55E-14,-5.09E-09,1.99E-11,-5.09E-09,-24.10942533,-4.82E-09,-4.78E-09 +320,0.448293315,-0.0032,-0.001342397,9.28E-20,-6.40E-19,0.000241631,-3.06E-18,0.000241631,6.30E-13,0.456681414,0.002021516,0.002021516,-0.001857603,24.95279339,-1.20E-15,1.01E-14,-24.95279339,1.20E-15,-1.01E-14,-5.47E-10,-4.82E-14,-5.47E-10,-23.47731101,-5.15E-10,-5.15E-10 +330,0.465447025,-0.0033,-0.001346142,1.47E-19,-8.67E-20,0.000242306,3.89E-18,0.000242306,6.30E-13,0.478723545,0.002126264,0.002126264,-0.001953858,24.28033472,-2.13E-15,1.33E-15,-24.28033472,2.13E-15,-1.33E-15,-6.39E-10,5.95E-14,-6.39E-10,-22.81082494,-6.01E-10,-6.00E-10 +340,0.4824486,-0.0034,-0.001347935,8.00E-19,7.36E-19,0.000242628,-5.71E-18,0.000242628,6.30E-13,0.500615153,0.002233137,0.002233137,-0.002052065,23.57416623,-1.17E-14,-1.09E-14,-23.57416623,1.17E-14,1.09E-14,-1.66E-13,-8.47E-14,-1.82E-13,-22.11472704,-2.49E-13,-7.75E-14 +350,0.499274891,-0.0035,-0.001347813,-1.90E-20,-4.42E-20,0.000242606,-6.33E-17,0.000242606,6.30E-13,0.522294465,0.002342093,0.002342093,-0.002152187,22.83914111,3.34E-16,5.87E-16,-22.83914111,-3.34E-16,-5.87E-16,-8.27E-13,-9.09E-13,-8.25E-13,-21.39382242,-1.69E-12,1.44E-13 +360,0.515902638,-0.0036,-0.001345827,3.20E-20,-3.72E-20,0.000242249,-6.33E-18,0.000242249,6.30E-13,0.543701184,0.002453078,0.002453078,-0.002254173,22.08015565,-5.12E-16,5.46E-16,-22.08015565,5.12E-16,-5.46E-16,-9.64E-09,-8.80E-14,-9.64E-09,-20.65291506,-9.02E-09,-9.02E-09 +370,0.532308734,-0.0037,-0.001342035,3.10E-20,6.09E-20,0.000241566,-3.25E-18,0.000241566,6.30E-13,0.564777199,0.002566029,0.002566029,-0.002357965,21.30210089,-2.15E-16,-8.51E-16,-21.30210089,2.15E-16,8.51E-16,-9.92E-09,-4.39E-14,-9.92E-09,-19.89675844,-9.27E-09,-9.27E-09 +380,0.548470497,-0.0038,-0.001336502,5.93E-20,5.73E-20,0.00024057,-3.53E-18,0.00024057,6.30E-13,0.585467262,0.002680874,0.002680874,-0.002463498,20.50981323,-8.17E-16,-6.48E-16,-20.50981323,8.17E-16,6.48E-16,-9.10E-09,-4.58E-14,-9.10E-09,-19.13000669,-8.49E-09,-8.49E-09 +390,0.564365951,-0.0039,-0.001329305,-4.39E-21,1.19E-19,0.000239275,-5.88E-18,0.000239275,6.30E-13,0.605719647,0.00279753,0.00279753,-0.002570695,19.70802252,4.65E-17,-1.44E-15,-19.70802252,-4.65E-17,1.44E-15,-8.94E-09,-7.41E-14,-8.94E-09,-18.35716479,-8.32E-09,-8.32E-09 +400,0.579974111,-0.004,-0.001320524,1.16E-19,1.46E-19,0.000237694,-3.73E-18,0.000237694,6.30E-13,0.625486754,0.002915909,0.002915909,-0.002679476,18.90130194,-1.16E-15,-1.72E-15,-18.90130194,1.16E-15,1.72E-15,-8.21E-09,-4.53E-14,-8.21E-09,-17.5825413,-7.64E-09,-7.64E-09 +410,0.595275243,-0.0041,-0.001310249,2.93E-20,-4.09E-20,0.000235845,-4.68E-18,0.000235845,6.30E-13,0.64472563,0.003035915,0.003035915,-0.002789751,18.09401946,-3.13E-16,3.95E-16,-18.09401946,3.13E-16,-3.95E-16,-8.32E-09,-5.48E-14,-8.32E-09,-16.8102049,-7.73E-09,-7.73E-09 +420,0.610251115,-0.0042,-0.001298575,-2.26E-19,-1.19E-19,0.000233744,-4.75E-18,0.000233744,6.30E-13,0.663398396,0.003157443,0.003157443,-0.002901425,17.29029338,2.62E-15,1.35E-15,-17.29029338,-2.62E-15,-1.35E-15,-8.07E-09,-5.35E-14,-8.07E-09,-16.04394575,-7.49E-09,-7.49E-09 +430,0.62488521,-0.0043,-0.0012856,-8.23E-21,-6.81E-20,0.000231408,-6.51E-17,0.000231408,4.83E-12,0.681472568,0.003280387,0.003280387,-0.0030144,16.49395305,7.20E-17,7.40E-16,-16.49395305,-7.20E-17,-7.40E-16,3.72E-09,-7.07E-13,3.72E-09,-15.28724308,3.45E-09,3.45E-09 +440,0.639162904,-0.0044,-0.001271425,2.15E-19,4.07E-20,0.000228857,6.48E-17,0.000228857,1.07E-11,0.698921251,0.003404636,0.003404636,-0.003128575,15.70850686,-2.26E-15,-4.26E-16,-15.70850686,2.26E-15,4.26E-16,5.07E-09,6.78E-13,5.07E-09,-14.54324106,4.69E-09,4.69E-09 +450,0.653071611,-0.0045,-0.001256155,-6.02E-20,-1.16E-20,0.000226108,-2.15E-16,0.000226108,1.85E-11,0.715723226,0.003530077,0.003530077,-0.003243845,14.93711621,5.94E-16,1.17E-16,-14.93711621,-5.94E-16,-1.17E-16,6.62E-09,-2.17E-12,6.62E-09,-13.81473032,6.12E-09,6.12E-09 +460,0.666600864,-0.0046,-0.001239895,2.65E-18,3.01E-18,0.000223181,-5.20E-18,0.000223181,1.85E-11,0.731862896,0.003656597,0.003656597,-0.003360105,14.18257945,-2.63E-14,-2.91E-14,-14.18257945,2.63E-14,2.91E-14,-2.09E-13,-5.04E-14,-2.13E-13,-13.10413978,-2.46E-13,-1.45E-13 +470,0.67974237,-0.0047,-0.001222747,-2.69E-20,5.31E-19,0.000220095,-3.51E-18,0.000220095,1.85E-11,0.747330141,0.003784081,0.003784081,-0.003477253,13.44732218,4.39E-17,-4.95E-15,-13.44732218,-4.39E-17,4.95E-15,-2.71E-14,-3.27E-14,-2.98E-14,-12.4135345,-5.90E-14,6.41E-15 +480,0.692490003,-0.0048,-0.001204816,-3.34E-18,-2.18E-18,0.000216867,-3.26E-18,0.000216867,1.85E-11,0.76212005,0.003912419,0.003912419,-0.003595184,12.73339745,2.98E-14,1.94E-14,-12.73339745,-2.98E-14,-1.94E-14,-3.60E-14,-2.92E-14,-4.25E-14,-11.74462255,-7.62E-14,3.85E-15 +490,0.70483976,-0.0049,-0.001186199,-6.08E-18,-5.61E-18,0.000213516,-3.84E-18,0.000213516,1.85E-11,0.776232588,0.004041502,0.004041502,-0.003713801,12.04249192,5.16E-14,4.85E-14,-12.04249192,-5.16E-14,-4.85E-14,-6.72E-14,-3.28E-14,-6.48E-14,-11.09876697,-1.29E-13,7.02E-15 +500,0.716789685,-0.005,-0.001166995,-1.19E-17,-1.16E-17,0.000210059,-4.30E-18,0.000210059,1.85E-11,0.789672185,0.004171225,0.004171225,-0.003833005,11.37593971,9.75E-14,9.62E-14,-11.37593971,-9.75E-14,-9.62E-14,-1.06E-13,-3.55E-14,-1.00E-13,-10.47700401,-2.06E-13,1.55E-14 \ No newline at end of file diff --git a/test/tests/damage_plasticity_model/gold/uniaxial_tension_out.csv b/test/tests/damage_plasticity_model/gold/uniaxial_tension_out.csv new file mode 100644 index 00000000..1f1b54ba --- /dev/null +++ b/test/tests/damage_plasticity_model/gold/uniaxial_tension_out.csv @@ -0,0 +1,32 @@ +time,D,displacement_x,e_xx,e_xy,e_xz,e_yy,e_yz,e_zz,intnl0,intnl1,max_ep,mid_ep,min_ep,react_x,react_y,react_z,s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,sigma0,sigma1,sigma2 +0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 +5,0,-5.00E-05,-5.00E-05,-8.01E-22,-5.62E-22,9.00E-06,5.36E-22,9.00E-06,0,0,0,0,0,1.585,2.15E-17,1.51E-17,-1.585,-2.15E-17,-1.51E-17,1.15E-13,1.44E-17,1.15E-13,0,0,0 +10,0,-0.0001,-0.0001,-2.38E-22,1.11E-21,1.80E-05,-5.36E-22,1.80E-05,0,0,0,0,0,3.17,6.40E-18,-2.99E-17,-3.17,-6.40E-18,2.99E-17,2.04E-14,-1.44E-17,2.03E-14,0,0,0 +15,0,-0.00015,-0.00015,-1.24E-21,-1.27E-21,2.70E-05,7.81E-22,2.70E-05,0,0,0,0,0,4.755,3.34E-17,3.41E-17,-4.755,-3.34E-17,-3.41E-17,1.51E-15,2.10E-17,1.44E-15,0,0,0 +20,0,-0.0002,-0.0002,-2.67E-21,-6.35E-22,3.60E-05,8.74E-22,3.60E-05,0,0,0,0,0,6.34,7.18E-17,1.71E-17,-6.34,-7.18E-17,-1.71E-17,6.40E-14,2.35E-17,6.40E-14,0,0,0 +25,0,-0.00025,-0.00025,-5.29E-22,-4.55E-21,4.50E-05,2.09E-21,4.50E-05,0,0,0,0,0,7.925,1.42E-17,1.22E-16,-7.925,-1.42E-17,-1.22E-16,-2.18E-13,5.62E-17,-2.18E-13,0,0,0 +30,0,-0.0003,-0.0003,-1.16E-21,-2.54E-21,5.40E-05,-9.26E-22,5.40E-05,0,0,0,0,0,9.51,3.13E-17,6.83E-17,-9.51,-3.13E-17,-6.83E-17,-4.37E-13,-2.49E-17,-4.37E-13,0,0,0 +35,0,-0.00035,-0.00035,-1.69E-21,-1.06E-21,6.30E-05,3.31E-21,6.30E-05,0,0,0,0,0,11.095,4.55E-17,2.84E-17,-11.095,-4.55E-17,-2.84E-17,-1.52E-13,8.89E-17,-1.52E-13,0,0,0 +40,0,-0.0004,-0.0004,-1.27E-21,6.35E-22,7.20E-05,-2.41E-21,7.20E-05,0,0,0,0,0,12.68,3.41E-17,-1.71E-17,-12.68,-3.41E-17,1.71E-17,-3.44E-13,-6.47E-17,-3.44E-13,0,0,0 +45,0,-0.00045,-0.00045,-2.12E-21,4.24E-22,8.10E-05,2.57E-21,8.10E-05,0,0,0,0,0,14.265,5.69E-17,-1.14E-17,-14.265,-5.69E-17,1.14E-17,3.18E-13,6.90E-17,3.18E-13,0,0,0 +50,0,-0.0005,-0.0005,-3.49E-21,-9.00E-21,9.00E-05,4.82E-21,9.00E-05,0,0,0,0,0,15.85,9.39E-17,2.42E-16,-15.85,-9.39E-17,-2.42E-16,-4.66E-13,1.29E-16,-4.66E-13,0,0,0 +55,0,-0.00055,-0.00055,-2.12E-21,6.35E-22,9.90E-05,1.11E-21,9.90E-05,0,0,0,0,0,17.435,5.69E-17,-1.71E-17,-17.435,-5.69E-17,1.71E-17,-1.62E-13,2.99E-17,-1.62E-13,0,0,0 +60,0.00350959,-0.0006,-0.000588195,2.43E-20,9.93E-21,0.000105875,4.41E-20,0.000105875,1.06E-15,0.001954762,1.28E-05,1.28E-05,-1.18E-05,18.64578867,-3.60E-17,-2.54E-16,-18.64578867,3.60E-17,2.54E-16,4.36E-12,1.21E-15,4.36E-12,-18.5803496,4.34E-12,4.35E-12 +65,0.011355885,-0.00065,-0.000612303,1.79E-20,1.61E-20,0.000110215,-6.99E-20,0.000110215,5.88E-15,0.006456462,4.10E-05,4.10E-05,-3.77E-05,19.41000639,-1.06E-15,-4.23E-16,-19.41000639,1.06E-15,4.23E-16,9.14E-12,-1.84E-15,9.14E-12,-19.18958859,9.03E-12,9.03E-12 +70,0.019228393,-0.0007,-0.000636103,2.35E-20,2.10E-20,0.000114499,2.57E-21,0.000114499,9.85E-15,0.01115344,6.95E-05,6.95E-05,-6.39E-05,20.09369807,5.20E-16,-4.42E-16,-20.09369807,-5.20E-16,4.42E-16,7.76E-12,6.25E-17,7.76E-12,-19.77673677,7.64E-12,7.64E-12 +75,0.027127075,-0.00075,-0.000659591,2.69E-21,4.63E-21,0.000118726,2.52E-20,0.000118726,1.31E-14,0.016044534,9.84E-05,9.84E-05,-9.04E-05,20.67160689,-8.20E-16,-1.80E-16,-20.67160689,8.20E-16,1.80E-16,6.60E-12,6.52E-16,6.60E-12,-20.34184634,6.49E-12,6.49E-12 +80,0.03505187,-0.0008,-0.000682764,3.34E-20,3.00E-20,0.000122898,-8.66E-20,0.000122898,1.59E-14,0.021128514,0.000127581,0.000127581,-0.000117236,21.22744789,1.64E-16,-7.37E-16,-21.22744789,-1.64E-16,7.37E-16,5.63E-12,-2.21E-15,5.63E-12,-20.88497057,5.54E-12,5.54E-12 +85,0.043002679,-0.00085,-0.000705617,-2.74E-20,-3.22E-21,0.000127011,1.83E-19,0.000127011,1.82E-14,0.026404062,0.000157123,0.000157123,-0.000144383,21.76127228,1.26E-15,1.69E-16,-21.76127228,-1.26E-15,-1.69E-16,4.82E-12,4.76E-15,4.82E-12,-21.40616595,4.74E-12,4.75E-12 +90,0.050979365,-0.0009,-0.000728145,-2.36E-21,1.96E-20,0.000131066,-4.91E-19,0.000131066,2.01E-14,0.031869774,0.000187019,0.000187019,-0.000171855,22.27313447,1.16E-16,-4.25E-16,-22.27313447,-1.16E-16,4.25E-16,4.14E-12,-1.27E-14,4.14E-12,-21.9054927,4.06E-12,4.09E-12 +95,0.058981754,-0.00095,-0.000750346,2.63E-21,-9.75E-21,0.000135062,1.14E-18,0.000135062,2.17E-14,0.037524154,0.000217272,0.000217272,-0.000199654,22.76309295,-1.73E-17,2.80E-16,-22.76309295,1.73E-17,-2.80E-16,3.57E-12,2.93E-14,3.57E-12,-22.38301543,3.48E-12,3.54E-12 +100,0.067009666,-0.001,-0.000772212,-2.86E-20,-2.91E-20,0.000138998,-3.17E-18,0.000138998,2.49E-14,0.043365632,0.000247887,0.000247887,-0.000227788,23.23120768,7.01E-16,7.20E-16,-23.23120768,-7.01E-16,-7.20E-16,6.99E-12,-8.09E-14,6.99E-12,-22.83879971,6.79E-12,6.96E-12 +105,0.075062766,-0.00105,-0.000793742,2.01E-20,6.23E-21,0.000142874,7.80E-18,0.000142874,2.73E-14,0.049392452,0.00027887,0.00027887,-0.000256258,23.67755258,-4.89E-16,-9.81E-17,-23.67755258,4.89E-16,9.81E-17,5.65E-12,1.97E-13,5.65E-12,-23.27292811,5.35E-12,5.76E-12 +110,0.083140825,-0.0011,-0.00081493,1.01E-21,-6.10E-21,0.000146687,-2.56E-17,0.000146687,2.93E-14,0.055602858,0.000310224,0.000310224,-0.00028507,24.1021984,-9.81E-18,1.53E-16,-24.1021984,9.81E-18,-1.53E-16,4.58E-12,-6.42E-13,4.58E-12,-23.68547769,3.82E-12,5.18E-12 +115,0.091243516,-0.00115,-0.000835771,1.45E-20,6.40E-22,0.000150439,7.19E-17,0.000150439,3.09E-14,0.061994947,0.000341956,0.000341956,-0.000314229,24.5052274,-3.64E-16,-1.59E-17,-24.5052274,3.64E-16,1.59E-17,3.73E-12,1.79E-12,3.73E-12,-24.07653564,1.78E-12,5.56E-12 +120,0.099370476,-0.0012,-0.00085626,4.28E-20,2.18E-20,0.000154127,-2.64E-16,0.000154127,3.33E-14,0.068566699,0.000374071,0.000374071,-0.00034374,24.88672472,-1.06E-15,-5.36E-16,-24.88672472,1.06E-15,5.36E-16,3.05E-12,-6.51E-12,3.05E-12,-24.44619594,-3.90E-12,9.89E-12 +125,0.107521297,-0.00125,-0.000876394,8.37E-21,6.97E-21,0.000157751,3.21E-18,0.000157751,4.10E-14,0.075315983,0.000406573,0.000406573,-0.000373606,25.24678399,-1.95E-16,-1.59E-16,-25.24678399,1.95E-16,1.59E-16,1.85E-11,7.86E-14,1.85E-11,-24.79455984,1.81E-11,1.82E-11 +130,0.115695536,-0.0013,-0.000896166,-2.44E-21,-5.34E-20,0.00016131,1.26E-18,0.00016131,4.18E-14,0.082240547,0.000439468,0.000439468,-0.000403834,25.58550629,-2.69E-16,1.29E-15,-25.58550629,2.69E-16,-1.29E-15,2.06E-12,3.05E-14,2.06E-12,-25.12173631,1.99E-12,2.05E-12 +135,0.123892703,-0.00135,-0.000915573,1.72E-20,6.56E-20,0.000164803,7.54E-18,0.000164803,4.25E-14,0.089338013,0.00047276,0.00047276,-0.000434427,25.90300056,-4.23E-16,-1.57E-15,-25.90300056,4.23E-16,1.57E-15,1.70E-12,1.81E-13,1.70E-12,-25.4278424,1.48E-12,1.86E-12 +140,0.132112269,-0.0014,-0.000934609,9.96E-20,7.31E-20,0.00016823,-1.99E-17,0.00016823,4.31E-14,0.096605879,0.000506456,0.000506456,-0.000465391,26.19938403,-2.52E-15,-1.74E-15,-26.19938403,2.52E-15,1.74E-15,1.41E-12,-4.73E-13,1.41E-12,-25.71300371,8.82E-13,1.88E-12 +145,0.140353658,-0.00145,-0.00095327,2.23E-19,2.42E-19,0.000171589,9.13E-17,0.000171589,4.38E-14,0.104041511,0.000540561,0.000540561,-0.00049673,26.47478261,-5.26E-15,-5.70E-15,-26.47478261,5.26E-15,5.70E-15,1.17E-12,2.15E-12,1.17E-12,-25.97735474,-1.12E-12,3.42E-12 +150,0.14861625,-0.0015,-0.00097155,8.06E-20,1.01E-19,0.000174879,8.65E-17,0.000174879,4.85E-14,0.111642141,0.000575079,0.000575079,-0.00052845,26.72933131,-1.86E-15,-2.35E-15,-26.72933131,1.86E-15,2.35E-15,1.17E-11,2.02E-12,1.17E-11,-26.22103932,9.35E-12,1.36E-11 \ No newline at end of file diff --git a/test/tests/damage_plasticity_model/tests b/test/tests/damage_plasticity_model/tests new file mode 100644 index 00000000..cb8d4b43 --- /dev/null +++ b/test/tests/damage_plasticity_model/tests @@ -0,0 +1,26 @@ +[Tests] + issues = '#184' + design = 'DamagePlasticityStressUpdate.md ComputeMultipleInelasticDamageStress.md' + [Damage_Plasticity_Model] + requirement = 'Blackbear shall demonstrate stress-strain behavior of concrete ' + [uniaxial_tension] + type = 'CSVDiff' + input = 'uniaxial_test.i' + csvdiff = 'uniaxial_tension_out.csv' + allow_test_objects = True + cli_args = 'Outputs/file_base=uniaxial_tension_out' + detail = ' in uniaxial tensile loading.' + [] + [uniaxial_compression] + type = 'CSVDiff' + input = 'uniaxial_test.i' + csvdiff = 'uniaxial_compression_out.csv' + allow_test_objects = True + cli_args = 'BCs/x_r/function=-2E-5*x*t + Executioner/end_time=500 + Executioner/dt=10 + Outputs/file_base=uniaxial_compression_out' + detail = ' in uniaxial compressive loading.' + [] + [] +[] diff --git a/test/tests/damage_plasticity_model/uniaxial_test.i b/test/tests/damage_plasticity_model/uniaxial_test.i new file mode 100644 index 00000000..845cf53b --- /dev/null +++ b/test/tests/damage_plasticity_model/uniaxial_test.i @@ -0,0 +1,470 @@ +[Mesh] + type = GeneratedMesh + dim = 3 + nx = 1 + ny = 1 + nz = 1 + xmin = -0.5 + xmax = 0.5 + ymin = -0.5 + ymax = 0.5 + zmin = -0.5 + zmax = 0.5 +[] + +[Variables] + [disp_x] + [] + [disp_y] + [] + [disp_z] + [] +[] + +[Kernels] + [TensorMechanics] + displacements = 'disp_x disp_y disp_z' + save_in = 'resid_x resid_y resid_z' + [] +[] + +[AuxVariables] + [resid_x] + [] + [resid_y] + [] + [resid_z] + [] + [stress_xx] + order = CONSTANT + family = MONOMIAL + [] + [stress_xy] + order = CONSTANT + family = MONOMIAL + [] + [stress_xz] + order = CONSTANT + family = MONOMIAL + [] + [stress_yy] + order = CONSTANT + family = MONOMIAL + [] + [stress_yz] + order = CONSTANT + family = MONOMIAL + [] + [stress_zz] + order = CONSTANT + family = MONOMIAL + [] + [strain_xx] + order = CONSTANT + family = MONOMIAL + [] + [strain_xy] + order = CONSTANT + family = MONOMIAL + [] + [strain_xz] + order = CONSTANT + family = MONOMIAL + [] + [strain_yy] + order = CONSTANT + family = MONOMIAL + [] + [strain_yz] + order = CONSTANT + family = MONOMIAL + [] + [strain_zz] + order = CONSTANT + family = MONOMIAL + [] + [min_ep] + order = CONSTANT + family = MONOMIAL + [] + [mid_ep] + order = CONSTANT + family = MONOMIAL + [] + [max_ep] + order = CONSTANT + family = MONOMIAL + [] + [sigma0] + order = CONSTANT + family = MONOMIAL + [] + [sigma1] + order = CONSTANT + family = MONOMIAL + [] + [sigma2] + order = CONSTANT + family = MONOMIAL + [] + [f0] + order = CONSTANT + family = MONOMIAL + [] + [D] + order = CONSTANT + family = MONOMIAL + [] + [intnl0] + order = CONSTANT + family = MONOMIAL + [] + [intnl1] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + [stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [] + [stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [] + [stress_xz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xz + index_i = 0 + index_j = 2 + [] + [stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [] + [stress_yz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yz + index_i = 1 + index_j = 2 + [] + [stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [] + [strain_xx] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [] + [strain_xy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [] + [strain_xz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xz + index_i = 0 + index_j = 2 + [] + [strain_yy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [] + [strain_yz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yz + index_i = 1 + index_j = 2 + [] + [strain_zz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_zz + index_i = 2 + index_j = 2 + [] + [f0_auxk] + type = MaterialStdVectorAux + property = plastic_yield_function + index = 0 + variable = f0 + [] + [D_auxk] + type = MaterialRealAux + property = elemental_damage_variable + variable = D + [] + [min_ep] + type = MaterialRealAux + property = min_ep + variable = min_ep + [] + [mid_ep] + type = MaterialRealAux + property = mid_ep + variable = mid_ep + [] + [max_ep] + type = MaterialRealAux + property = max_ep + variable = max_ep + [] + [sigma0] + type = MaterialRealAux + property = damaged_min_principal_stress + variable = sigma0 + [] + [sigma1] + type = MaterialRealAux + property = damaged_mid_principal_stress + variable = sigma1 + [] + [sigma2] + type = MaterialRealAux + property = damaged_max_principal_stress + variable = sigma2 + [] + [intnl0_auxk] + type = MaterialRealAux + property = damage_state_in_tension + variable = intnl0 + [] + [intnl1_auxk] + type = MaterialRealAux + property = damage_state_in_compression + variable = intnl1 + [] +[] + +[Postprocessors] + [react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [] + [react_y] + type = NodalSum + variable = resid_y + boundary = 'left' + [] + [react_z] + type = NodalSum + variable = resid_z + boundary = 'left' + [] + [displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [] + [s_xx] + type = PointValue + point = '0 0 0' + variable = stress_xx + [] + [s_xy] + type = PointValue + point = '0 0 0' + variable = stress_xy + [] + [s_xz] + type = PointValue + point = '0 0 0' + variable = stress_xz + [] + [s_yy] + type = PointValue + point = '0 0 0' + variable = stress_yy + [] + [s_yz] + type = PointValue + point = '0 0 0' + variable = stress_yz + [] + [s_zz] + type = PointValue + point = '0 0 0' + variable = stress_zz + [] + [e_xx] + type = PointValue + point = '0 0 0' + variable = strain_xx + [] + [e_xy] + type = PointValue + point = '0 0 0' + variable = strain_xy + [] + [e_xz] + type = PointValue + point = '0 0 0' + variable = strain_xz + [] + [e_yy] + type = PointValue + point = '0 0 0' + variable = strain_yy + [] + [e_yz] + type = PointValue + point = '0 0 0' + variable = strain_yz + [] + [e_zz] + type = PointValue + point = '0 0 0' + variable = strain_zz + [] + [max_ep] + type = PointValue + point = '0 0 0' + variable = max_ep + [] + [mid_ep] + type = PointValue + point = '0 0 0' + variable = mid_ep + [] + [min_ep] + type = PointValue + point = '0 0 0' + variable = min_ep + [] + [sigma0] + type = PointValue + point = '0 0 0' + variable = sigma0 + [] + [sigma1] + type = PointValue + point = '0 0 0' + variable = sigma1 + [] + [sigma2] + type = PointValue + point = '0 0 0' + variable = sigma2 + [] + [D] + type = PointValue + point = '0 0 0' + variable = D + [] + [intnl0] + type = PointValue + point = '0 0 0' + variable = intnl0 + [] + [intnl1] + type = PointValue + point = '0 0 0' + variable = intnl1 + [] +[] + +[BCs] + [x_r] + type = FunctionDirichletBC + variable = disp_x + boundary = 'right' + function = '-2E-5*x*t' + [] + [x_l] + type = FunctionDirichletBC + variable = disp_x + boundary = 'left' + function = '0' + [] + [y_l] + type = FunctionDirichletBC + variable = disp_y + boundary = 'bottom' + function = '0' + [] + [z_l] + type = FunctionDirichletBC + variable = disp_z + boundary = 'back' + function = '0' + [] +[] + +[Materials] + [elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [] + [strain] + type = ComputeIncrementalSmallStrain + displacements = 'disp_x disp_y disp_z' + [] + [pdm] + type = DamagePlasticityStressUpdate + factor_relating_biaxial_unixial_cmp_str = 0.109 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_strength_in_tension = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.70 + tensile_damage_at_half_tensile_strength = 0.51 + fracture_energy_in_tension = 12.3E-3 + + yield_strength_in_compression = 18.30 + maximum_strength_in_compression = 27.60 + compressive_damage_at_max_compressive_strength = 0.25 + fracture_energy_in_compression = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [] + [stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + [] +[] + +[Executioner] + end_time = 150 + dt = 5 + type = Transient +[] + +[Outputs] + file_base = ./test/tests/damage_plasticity_model/uniaxial_tension_out + exodus = false + [csv] + type = CSV + [] +[]