The elasticity is assumed linear and isotropic, i.e. given by the Hooke law:
\[ \begin{aligned} \underline{\sigma}&=\underline{\underline{\mathbf{D}}}\,\colon\,\underline{\varepsilon}^{\mathrm{el}}\\ &=\lambda\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\,\underline{I}+2\,\mu\,\underline{\varepsilon}^{\mathrm{el}}\\ &= K^{\star}\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\,\underline{I}+2\,\mu\,{\left(\underline{\varepsilon}^{\mathrm{el}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\,\underline{I}\right)} \\ &= K^{\star}\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\,\underline{I}+2\,\mu\,\underline{s}^{\mathrm{el}} \end{aligned} \] where:
For crushed salt, the buldmodulus is a function of the porosity and can be defined by: \[ K^{\star} = K \cdot e^{-c_k \cdot \eta \cdot \left(\dfrac{1-\eta_0}{1-\eta}\right)} \] where:
The Poisson’s ratio \(\nu^{\star}\) of crushed salt and rock salt \(\nu\) are assumed equal, thus: \[ \nu^{\star} = \nu \]
The Young’ modulus \(E^{\star}\) of crushed salt can thus be derived as: \[ \begin{aligned} E^{\star} &= 3 (1-2 \cdot \nu^{\star}) \cdot K^{\star} \\ &= 3 (1-2 \cdot \nu) \cdot K \cdot e^{-c_k \cdot \eta \cdot \left(\dfrac{1-\eta_0}{1-\eta}\right)} \end{aligned} \]
The inelastic strain \(\underline{\varepsilon}^{\mathrm{in}}\) is split as the sum of two contributions \(\underline{\varepsilon}^{\mathrm{vp}}\) and \(\underline{\varepsilon}^{\mathrm{g}}\) which respectively describe the viscoplastic deformation of single salt grains and the relative displacement between grains, as follows: \[ \underline{\dot{\varepsilon}}^{\mathrm{in}}= \underline{\dot{\varepsilon}}^{\mathrm{vp}}+ \underline{\dot{\varepsilon}}^{\mathrm{g}} \qquad{(1)}\]
In repository conditions, the grain displacement phenomenon can be neglected, but it may a role in certain conditions.
The grain deformation strain rate tensor \(\underline{\dot{\varepsilon}}^{\mathrm{vp}}\) follows an associated Norton-Hoff behaviour based on a Green stress criterion. This criterion is expressed as: \[ \sigma_{\mathrm{eq}}= \sqrt{h_1{\left(\eta\right)} \, p^{2}+h_2{\left(\eta\right)} \, q^{2}} \] where:
Various models can be derived from the choice of the functions \(h_1{\left(\eta\right)}\) and \(h_2{\left(\eta\right)}\). In Korthaus’s work, the following expressions are used: \[ \begin{aligned} h_1{\left(\eta\right)} &= {{\displaystyle \frac{\displaystyle a}{\displaystyle \left( \eta^{-c}-\eta_0^{-c}\right)^m}}} \\ h_2{\left(\eta\right)} &= b_1 + b_2 \, h_1 \end{aligned} \] with:
A potential numerical issue
The expression of \(h_{1}\) diverges as the porosity \(\eta\) tends to the reference porosity \(\eta_{0}\). To avoid numerical issues, a numerical parameter \(\Delta\) is introduced in the implementation.
The normal \(\underline{n}\) to the Green criterion is given by: \[ \begin{aligned} \underline{n}&={\displaystyle \frac{\displaystyle \partial \sigma_{\mathrm{eq}}}{\displaystyle \partial \underline{\sigma}}}={{\displaystyle \frac{\displaystyle 1}{\displaystyle 2\,\sigma_{\mathrm{eq}}}}}\,{\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial \underline{\sigma}}}{\left(h_1{\left(\eta\right)} \, p^{2}+h_2{\left(\eta\right)} \, q^{2}\right)}\\ &={{\displaystyle \frac{\displaystyle 1}{\displaystyle 2\,\sigma_{\mathrm{eq}}}}}{\left(2\,p\,{\displaystyle \frac{\displaystyle \partial p}{\displaystyle \partial \underline{\sigma}}}+2\,q\,{\displaystyle \frac{\displaystyle \partial q}{\displaystyle \partial \underline{\sigma}}}\right)}\\ &={{\displaystyle \frac{\displaystyle 1}{\displaystyle \sigma_{\mathrm{eq}}}}}{\left({{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}} \, h_1{\left(\eta\right)} \, p \, \underline{I} + h_2{\left(\eta\right)} \, \underline{s}\right)} \end{aligned} \]
Finally, the viscoplastic strain rate can be expressed as follows:
\[ \underline{\dot{\varepsilon}}^{\mathrm{vp}}= A_{\mathrm{vp}} \, \exp{\left(-{{\displaystyle \frac{\displaystyle Q_c}{\displaystyle R_m \, T}}}\right)} \, \sigma_{\mathrm{eq}}^{n_{\mathrm{vp}}}\,\underline{n} \qquad{(2)}\]
with:
The expression may we rewritten as follows:
\[ \underline{\dot{\varepsilon}}^{\mathrm{vp}}= \dot{\varepsilon}_{0} \, \exp{\left(-{{\displaystyle \frac{\displaystyle T_{a}}{\displaystyle T}}}\right)} \, {\left({{\displaystyle \frac{\displaystyle \sigma_{\mathrm{eq}}}{\displaystyle \sigma_{0}}}}\right)}^{n_{\mathrm{vp}}}\,\underline{n} \]
where:
The porosity evolution is given by: \[ \dot{\eta} = \left(1-\eta\right)\,{\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{to}}\right)}} \]
The porosity is however bounded between \(0\) and the reference porosity \(\eta_{0}\).
This ordinary differential equation can be integrated exactly between \(t\) and \(t+\theta\,\Delta\,t\) by separation of variables: \[ \int_{t}^{t+\theta\,\Delta\,t}{{\displaystyle \frac{\displaystyle \dot{\eta}}{\displaystyle 1-\eta}}}\,\mathrm{d}t = \int_{t}^{t+\theta\,\Delta\,t}{\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{to}}\right)}}\,\mathrm{d}t \Leftrightarrow -\log{\left({{\displaystyle \frac{\displaystyle 1-{\left.\eta\right|_{t+\theta\,\Delta\,t}}}{\displaystyle 1-{\left.\eta\right|_{t}}}}}\right)} = \theta\,{\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}} \]
Finally, the porosity at the middle of the time step and the porosity at the end of time step are given by: \[ {\left.\eta\right|_{t+\theta\,\Delta\,t}} = 1 - {\left(1-{\left.\eta\right|_{t}}\right)}\,\exp{\left(-\theta\,{\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}}\right)} \qquad{(3)}\] and \[ {\left.\eta\right|_{t+\Delta\,t}} = 1 - {\left(1-{\left.\eta\right|_{t}}\right)}\,\exp{\left(-{\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}}\right)} \qquad{(4)}\]
For the computation of the consistent tangent operator, the expression of the derivatives of those quantities with respect to the total strain increment \(\Delta\underline{\varepsilon}^{\mathrm{to}}\) will be required: \[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} &= {\left(1-{\left.\eta\right|_{t}}\right)}\,\theta\,\exp{\left(-\theta\,{\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}}\right)}\,\underline{I} = \theta\,{\left(1-{\left.\eta\right|_{t+\theta\,\Delta\,t}}\right)}\,\underline{I}\\ {\displaystyle \frac{\displaystyle \partial {\left.\eta\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} &= {\left(1-{\left.\eta\right|_{t}}\right)}\,\exp{\left({\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}}\right)}\,\underline{I} = {\left(1-{\left.\eta\right|_{t+\Delta\,t}}\right)}\,\underline{I} \\ \end{aligned} \right. \qquad{(5)}\]
The elastic strain is the only integration variable required.
The residual associated with the elastic strain can be expressed as follows:
\[ \begin{aligned} f_{\underline{\varepsilon}^{\mathrm{el}}} &= \Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,\underline{\varepsilon}^{\mathrm{vp}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}\\ &= \Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,t\,\dot{\varepsilon}_{0} \, \exp{\left(-{{\displaystyle \frac{\displaystyle T_{a}}{\displaystyle {\left.T\right|_{t+\theta\,\Delta\,t}}}}}\right)} \, {\left({{\displaystyle \frac{\displaystyle {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \sigma_{0}}}}\right)}^{n_{\mathrm{vp}}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}\\ &= \Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,t\,f_{\mathrm{vp}}{\left({\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}\right)}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}} \end{aligned} \] where \(f_{\mathrm{vp}}{\left({\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}\right)}=\dot{\varepsilon}_{0} \, \exp{\left(-{{\displaystyle \frac{\displaystyle T_{a}}{\displaystyle {\left.T\right|_{t+\theta\,\Delta\,t}}}}}\right)} \, {\left({{\displaystyle \frac{\displaystyle {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \sigma_{0}}}}\right)}^{n_{\mathrm{vp}}}\)
The jacobian of the implicit system is given by: \[ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} =\underline{\underline{\mathbf{I}}}+ \Delta\,t\,\underline{n}\,\otimes\,{\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}+ \Delta\,t\,f_{\mathrm{vp}}{\left({\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}\right)}\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} \]
where:
\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &={\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}\,{\displaystyle \frac{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\\ &=\theta\,{\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}}\\ &=\theta\,{\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \]
\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &={\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\\ &=\theta\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}}\\ &=\theta\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \]
While the implementation of the behaviour integration per se is trivial, the computation of the consistent tangent operator is particularly involved.
\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} &= {\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\left[K^{\star}{\left({\left.\eta\right|_{t+\Delta\,t}}\right)}\,{\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\Delta\,t}}\right)}}\,\underline{I}+2\,\mu{\left({\left.\eta\right|_{t+\Delta\,t}}\right)}\,{\left.\underline{s}^{\mathrm{el}}\right|_{t+\Delta\,t}}\,\underline{I}\right]\\ &= {\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\Delta\,t}}\right)}}\,\underline{I}\,\otimes\,{\displaystyle \frac{\displaystyle \partial K^{\star}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}+ 2\,{\left.\underline{s}^{\mathrm{el}}\right|_{t+\Delta\,t}}\,\otimes\,{\displaystyle \frac{\displaystyle \partial \mu^{\star}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}+ {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\\ &= {\left({\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\Delta\,t}}\right)}}\,{\displaystyle \frac{\displaystyle \partial K^{\star}}{\displaystyle \partial {\left.f\right|_{t+\Delta\,t}}}}\,\underline{I}+ 2\,{\displaystyle \frac{\displaystyle \partial \mu^{\star}}{\displaystyle \partial {\left.f\right|_{t+\Delta\,t}}}}\,{\left.\underline{s}^{\mathrm{el}}\right|_{t+\Delta\,t}}\right)}\,\otimes\,{\displaystyle \frac{\displaystyle \partial {\left.\eta\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}+ {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \end{aligned} \]
The derivative \({\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) can be obtained by the implicit function theorem as follows:
\[ \begin{aligned} &&f_{\underline{\varepsilon}^{\mathrm{el}}}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{el}}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)},\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)} = 0\\ &\Rightarrow&\quad{\displaystyle \frac{\displaystyle \mathrm{d}f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}} = 0\\ &\Rightarrow& {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}=-{\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\\ &\Rightarrow& {\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}=-{\left({\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\right)}^{-1}\,\colon\,{\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \end{aligned} \]
with:
\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} &=-\underline{\underline{\mathbf{I}}}+\Delta\,t\,{\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\left[f_{\mathrm{vp}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right]\\ &=-\underline{\underline{\mathbf{I}}}+\Delta\,t\,\left[{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}+f_{\mathrm{vp}}\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\right]\\ &=-\underline{\underline{\mathbf{I}}}+\Delta\,t\,\left[ {\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}+f_{\mathrm{vp}}\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}\right]\,\otimes\,{\displaystyle \frac{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\\ \end{aligned} \]
The derivative of \(f_{\mathrm{vp}}\) with respect to the porosity \({\left.\eta\right|_{t+\theta\,\Delta\,t}}\) is: \[ {\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}={\displaystyle \frac{\displaystyle \partial f_{\mathrm{vp}}}{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}\,{\displaystyle \frac{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}} \]
The partial derivative of the equivalent stress \({\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}\) with respect to the porosity \({\left.\eta\right|_{t+\theta\,\Delta\,t}}\) step is given by:
\[ {\displaystyle \frac{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle 2\,\sigma_{\mathrm{eq}}}}}\left[p^{2}\,{\displaystyle \frac{\displaystyle \partial h_{1}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}+q^{2}\,{\displaystyle \frac{\displaystyle \partial h_{2}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}\right] \]
The partial derivative of the normal \({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\) with respect to the porosity \({\left.\eta\right|_{t+\theta\,\Delta\,t}}\) is given by: \[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}} &=&{\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}\left[{{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left({{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}} \, h_1{\left({\left.\eta\right|_{t+\theta\,\Delta\,t}}\right)} \, p \, \underline{I} + h_2{\left({\left.\eta\right|_{t+\theta\,\Delta\,t}}\right)} \, \underline{s}\right)}\right]\\ &=&{{\displaystyle \frac{\displaystyle -1}{\displaystyle {\left.\sigma_{\mathrm{eq}}^{2}\right|_{t+\theta\,\Delta\,t}}}}}\,\,{\displaystyle \frac{\displaystyle \partial \sigma_{\mathrm{eq}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}\,{\left({{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}} \, h_1{\left({\left.\eta\right|_{t+\theta\,\Delta\,t}}\right)} \, p \, \underline{I} + h_2{\left({\left.\eta\right|_{t+\theta\,\Delta\,t}}\right)} \, \underline{s}\right)}\\ &&+{{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left({{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}} \, {\displaystyle \frac{\displaystyle \partial h_1}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}} \, p \, \underline{I} + {\displaystyle \frac{\displaystyle \partial h_2}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}} \, \underline{s}\right)}\\ &=&{{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,\left[{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}} \, {\displaystyle \frac{\displaystyle \partial h_1}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}} \, p \, \underline{I} + {\displaystyle \frac{\displaystyle \partial h_2}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}} \, \underline{s}-{\displaystyle \frac{\displaystyle \partial \sigma_{\mathrm{eq}}}{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}}\,{\left.n\right|_{t+\theta\,\Delta\,t}}\right]\\ \end{aligned} \]