6. Calculation of the coherent tangent operator#

In pure mechanics or in modeling THM with the Hoek-Brown law used according to its first aspect, the stress tensor represents the effective stress tensor which depends solely on the deformation tensor. The law of mechanical behavior thus only provides the derivative \(\frac{\partial \sigma }{\partial \varepsilon }\) in pure mechanics noted \(\frac{\partial \sigma \text{'}}{\partial \varepsilon }\) in thermo-hydro-mechanical modeling.

The calculation of the quantity \(\frac{\partial \sigma }{\partial \varepsilon }\) \(\frac{\partial \sigma }{\partial {p}_{c}}=\frac{\partial \sigma }{\partial {\sigma }_{p}}\frac{\partial {\sigma }_{p}}{\partial {p}_{c}}\) of the coherent tangent matrix is the same for both aspects of the law. The detail of this derivative is presented in sub-paragraph 6.1.

In modeling THM with the Hoek-Brown law used according to its second aspect, the stress tensor representing the total stress tensor, the variation in the stress tensor depends on the variation in the strain tensor \(\frac{\partial \sigma }{\partial \varepsilon }\), on the variation in the strain tensor, on the variation in the gas pressure \(\frac{\partial \sigma }{\partial {p}_{g}}\), and on the variation in the capillary pressure \(\frac{\partial \sigma }{\partial {p}_{c}}\).

Gas pressure and capillary pressure are not always involved in the calculation but depending on the case of the simulation. In the code, we will therefore calculate \(\frac{\partial \sigma }{\partial {p}_{p}}\) then \(\frac{\partial \sigma }{\partial {p}_{g}}=\frac{\partial \sigma }{\partial {\sigma }_{p}}\frac{\partial {\sigma }_{p}}{\partial {p}_{g}}\) and.

\(\frac{\partial {\sigma }_{p}}{\partial {p}_{c}}\) and \(\frac{\partial {\sigma }_{p}}{\partial {p}_{g}}\) depend on hydraulic simulation, and are independent of mechanical law (11).

6.1. Calculation of \(\frac{\partial \sigma }{\partial \varepsilon }\)#

We want to calculate the coherent matrix: \(\frac{\partial \sigma }{\partial \varepsilon }\text{=}\frac{\partial s}{\partial \varepsilon }\text{+}\frac{1}{3}\frac{\partial {I}_{1}}{\partial \varepsilon }I\) \(\frac{\partial {s}_{i}^{e}}{\partial {\varepsilon }_{\text{pq}}}\). Gold: \(\begin{array}{}\frac{\partial s}{\partial \varepsilon }\text{=}\frac{\partial {s}^{e}}{\partial \varepsilon }(1\text{-}\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}\Delta \lambda )\text{+}\frac{3\mu \Delta \lambda }{({\sigma }_{\text{eq}}^{e}{)}^{2}}{s}^{e}\frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial \varepsilon }\text{-}\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}{s}^{e}\frac{\partial \Delta \lambda }{\partial \varepsilon }\\ \frac{\partial {I}_{1}}{\partial \varepsilon }\text{=}\frac{\partial {I}_{1}^{e}}{\partial \varepsilon }\text{-}\mathrm{9K}(\frac{\partial \eta }{\partial \Delta \lambda }\frac{\partial \Delta \lambda }{\partial \varepsilon }\Delta \lambda \text{+}\eta \frac{\partial \Delta \lambda }{\partial \varepsilon })\end{array}\)

\(\tilde{P}\text{.}{s}^{e}\text{.}P=\stackrel{ˉ}{{s}^{e}}\)

with: \(\begin{array}{}\frac{\partial {s}_{\text{ij}}^{e}}{\partial {\varepsilon }_{\text{pq}}}=2\mu ({\delta }_{\text{ip}}{\delta }_{\text{jq}}-\frac{1}{3}{\delta }_{\text{ij}}{\delta }_{\text{pq}})\\ \frac{\partial {I}_{1}^{e}}{\partial {\varepsilon }_{\text{pq}}}=\mathrm{3K}{\delta }_{\text{pq}}\\ \frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial {\varepsilon }_{\text{pq}}}=\sqrt{\frac{3}{2}}\frac{\partial {s}_{\text{II}}^{e}}{\partial {\varepsilon }_{\text{pq}}}=\frac{3}{2{\sigma }_{\text{eq}}^{e}}\sum _{i,j}{s}_{\text{ij}}^{e}\frac{\partial {s}_{\text{ij}}^{e}}{\partial {\varepsilon }_{\text{pq}}}=\frac{3}{2{\sigma }_{\text{eq}}^{e}}2\mu \sum _{i,j}{s}_{\text{ij}}^{e}({\delta }_{\text{ip}}{\delta }_{\text{jq}}-\frac{1}{3}{\delta }_{\text{ij}}{\delta }_{\text{pq}})=\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}{s}_{\text{pq}}^{e}\\ \frac{\partial \eta }{\partial \Delta \lambda }=\frac{\partial \eta }{\partial \Delta \gamma }\frac{\partial \Delta \gamma }{\partial \Delta \lambda }=(\eta +1)\frac{\partial \eta }{\partial \Delta \gamma }\end{array}\)

\({\stackrel{ˉ}{s}}^{e}=\text{diag}({s}_{1}^{e},{s}_{2}^{e},{s}_{3}^{e})\)

from where: \(\frac{\partial \sigma }{\partial \varepsilon }\text{=}\frac{\partial {s}^{e}}{\partial \varepsilon }(1\text{-}\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}\Delta \lambda )\text{+}\frac{9{\mu }^{2}\Delta \lambda }{({\sigma }_{\text{eq}}^{e}{)}^{3}}{s}^{e}\text{.}{s}^{e}\text{-}\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}{s}^{e}\text{.}\frac{\partial \Delta \lambda }{\partial \varepsilon }\text{+}\frac{1}{3}\frac{\partial {I}_{1}^{e}}{\partial \varepsilon }\text{-}\mathrm{3K}(\frac{\partial \eta }{\partial \Delta \lambda }\frac{\partial \Delta \lambda }{\partial \varepsilon }\Delta \lambda \text{+}\eta \frac{\partial \Delta \lambda }{\partial \varepsilon })\)

\(P\)

To calculate \(\frac{\partial \Delta \lambda }{\partial {\varepsilon }_{\text{pq}}}\), we use the equation \(\dot{F}(\Delta \lambda )=0\) \(\frac{\partial {s}^{e}}{\partial {\varepsilon }_{\text{pq}}}=P\text{.}\frac{\partial {\stackrel{ˉ}{s}}^{e}}{\partial {\varepsilon }_{\text{pq}}}\text{.}\tilde{P}\). We get:

\(\begin{array}{c}\frac{\partial \Delta \lambda }{\partial {\varepsilon }_{\mathrm{pq}}}\left[\begin{array}{c}\text{-}({s}_{3}^{e}\text{-}{s}_{1}^{e})\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}}\text{-}\frac{\partial b}{\Delta \lambda }\left[1\text{+}\frac{1}{{\sigma }_{3}^{b\text{-}d}}\left[{s}_{3}^{e}\text{+}\frac{{I}_{1}^{e}}{3}\text{-}(3K\eta \text{+}{s}_{3}^{e}\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}})\Delta \lambda \right]\right]\\ \text{+}\frac{b}{{\sigma }_{3}^{b\text{-}d}}(\frac{\partial \eta }{\partial \Delta \gamma }3K(\eta \text{+}1)\Delta \lambda \text{+}(3K\eta \text{+}{s}_{3}^{e}\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}}))\\ \text{-}\frac{1}{2}(\begin{array}{c}\frac{\partial (S{\sigma }_{c}^{2})}{\partial \Delta \lambda }\text{-}\frac{\partial (m{\sigma }_{c})}{\partial \Delta \lambda }({s}_{3}^{e}\text{+}\frac{{I}_{1}^{e}}{3}\text{-}(3K\eta \text{+}{s}_{3}^{e}\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}})\Delta \lambda )\\ \text{+}{\sigma }_{c}m(\frac{\partial \eta }{\partial \Delta \gamma }3K(\eta \text{+}1)\Delta \lambda \text{+}(3K\eta \text{+}{s}_{3}^{e}\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}}))\end{array})\times \\ (S{\sigma }_{c}^{2}\text{-}m{\sigma }_{c}{({s}_{3}^{e}\text{+}\frac{{I}_{1}^{e}}{3}\text{-}(3K\eta \text{=}{s}_{3}^{e}\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}}\Delta \lambda ))}^{\text{-}\frac{1}{2}})\\ \text{=}\text{-}(\frac{\partial {s}_{3}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}\text{-}\frac{\partial {s}_{1}^{e}}{{s}_{3}^{e}})\left[1\text{-}\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}}\Delta \lambda \right]\text{-}({s}_{3}^{e}\text{-}{s}_{1}^{e})\frac{3\mu }{{({\sigma }_{\mathrm{eq}}^{e})}^{2}}\frac{\partial {\sigma }_{\mathrm{eq}}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}\Delta \lambda \\ \text{+}(\frac{\partial {s}_{3}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}\text{+}\frac{1}{3}\frac{\partial {I}_{1}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}\text{-}(\text{-}\frac{1}{{\sigma }_{\mathrm{eq}}^{e}}\frac{\partial {\sigma }_{\mathrm{eq}}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}{s}_{3}^{e}\text{+}\frac{\partial {s}_{3}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}})\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}})\times \\ (\frac{b}{{\sigma }_{3}^{b\text{-}d}}\text{-}\frac{1}{2}{\sigma }_{c}m{(S{\sigma }_{c}^{2}\text{-}m{\sigma }_{c}({s}_{3}^{e}\text{+}\frac{{I}_{1}^{e}}{3}\text{-}(3K\eta \text{+}{s}_{3}^{e}\frac{3\mu }{{\sigma }_{\mathrm{eq}}^{e}})\Delta \lambda ))}^{\text{-}\frac{1}{2}})\end{array}\right]\end{array}\)

So all we have to do is calculate \(\frac{\partial {s}_{i}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}\). Using the notations from paragraph 4.1, we have: \(\tilde{P}\text{.}{s}^{e}\text{.}P\text{=}{\stackrel{ˉ}{s}}^{e}\) where \({\stackrel{ˉ}{s}}^{e}\text{=}\mathrm{diag}({s}_{1}^{e},{s}_{2}^{e},{s}_{3}^{e})\) and \(P\) is the matrix of the associated eigenvectors. As a result, \(\frac{\partial {s}^{e}}{\partial {\varepsilon }_{\text{pq}}}\text{=}\tilde{P}\text{.}\frac{\partial {\stackrel{ˉ}{s}}^{e}}{\partial {\varepsilon }_{\text{pq}}}\text{.}P\text{+}\frac{\partial \tilde{P}}{\partial {\varepsilon }_{\text{pq}}}\text{.}{\stackrel{ˉ}{s}}^{e}\text{.}P\text{+}\tilde{P}\text{.}{\stackrel{ˉ}{s}}^{e}\text{.}\frac{\partial P}{\partial {\varepsilon }_{\text{pq}}}\).

Using the same reasoning as in Appendix 1, we can show that \(\frac{\partial {s}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}\text{=}P\text{.}\frac{\partial {\stackrel{ˉ}{s}}^{e}}{\partial {\varepsilon }_{\mathrm{pq}}}\text{.}\tilde{P}\) and finally:. \(\tilde{P}\text{.}\frac{\partial {s}^{e}}{\partial {\varepsilon }_{\text{pq}}}\text{.}P=\frac{\partial {\stackrel{ˉ}{s}}^{e}}{\partial {\varepsilon }_{\text{pq}}}\)

6.2. Special case of multiple eigenvalues for \({s}^{e}\)#

In the case where the \({s}^{e}\) matrix has multiple eigenvalues, the remarks made in paragraph 5.2.1 apply. In case \({s}_{3}^{e}\text{=}{s}_{2}^{e}\) for example, the vector \(\frac{\partial {s}_{3}^{e}}{\partial {\mathrm{\varepsilon }}_{\text{pq}}}\) therefore only takes into account one of the two directional derivatives of \({s}_{3}^{e}\text{=}{s}_{2}^{e}\). Hence the idea of writing \({s}_{2}^{e}\text{=}{s}_{3}^{e}\text{=}\frac{1}{2}({s}_{2}^{e}\text{+}{s}_{3}^{e})\) and therefore \(\frac{\partial {s}_{2}^{e}}{\partial {\varepsilon }_{\text{ij}}}\text{=}\frac{\partial {s}_{3}^{e}}{\partial {\varepsilon }_{\text{ij}}}\text{=}\frac{1}{2}(\frac{\partial {s}_{2}^{e}}{\partial {\varepsilon }_{\text{ij}}}\text{+}\frac{\partial {s}_{3}^{e}}{\partial {\varepsilon }_{\text{ij}}})\).

6.3. Calculation of \(\frac{\partial \mathrm{\sigma }}{\partial {\mathrm{\sigma }}_{p}}\)#

  1. Elasticity:

We have: \(\sigma \text{=}{\sigma }^{e}\text{=}H\varepsilon \text{+}{\sigma }_{p}I\) and then \(\frac{\partial \sigma }{\partial {\sigma }_{p}}\text{=}I\).

  1. Plasticity:

We have: \(\frac{\partial \sigma }{\partial {\sigma }_{p}}=\frac{\partial s}{\partial {\sigma }_{p}}+\frac{1}{3}\frac{\partial {I}_{1}}{\partial {\sigma }_{p}}I\) with:

\(\begin{array}{}\frac{\partial s}{\partial {\sigma }_{p}}=\frac{\partial {s}^{e}}{\partial {\sigma }_{p}}(1-\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}\Delta \lambda )+\frac{3\mu \Delta \lambda }{({\sigma }_{\text{eq}}^{e}{)}^{2}}{s}^{e}\frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial {\sigma }_{p}}-\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}{s}^{e}\frac{\partial \Delta \lambda }{\partial {\sigma }_{p}}\\ \frac{\partial {I}_{1}}{\partial {\sigma }_{p}}=\frac{\partial {I}_{1}^{e}}{\partial {\sigma }_{p}}-\mathrm{9K}(\frac{\partial \eta }{\partial \Delta \lambda }\frac{\partial \Delta \lambda }{\partial {\sigma }_{p}}\Delta \lambda +\eta \frac{\partial \Delta \lambda }{\partial {\sigma }_{p}})\end{array}\)

Gold:

\(\begin{array}{}\frac{\partial {s}^{e}}{\partial {\sigma }_{p}}=2\mu \frac{\partial \Delta e}{\partial {\sigma }_{p}}=0\\ \frac{\partial {I}_{1}^{e}}{\partial {\sigma }_{p}}=\mathrm{3K}\frac{\partial \Delta {\varepsilon }_{v}}{\partial {\sigma }_{p}}+3=3\\ \frac{\partial {\sigma }_{\text{eq}}^{e}}{\partial {\varepsilon }_{\text{pq}}}=\sqrt{\frac{3}{2}}\frac{\partial {s}_{\text{II}}^{e}}{\partial {\varepsilon }_{\text{pq}}}=0\\ \frac{\partial \eta }{\partial \Delta \lambda }=\frac{\partial \eta }{\partial \Delta \gamma }\frac{\partial \Delta \gamma }{\partial \Delta \lambda }=(\eta +1)\frac{\partial \eta }{\partial \Delta \gamma }\end{array}\)

from where: \(\frac{\partial \sigma }{\partial {\sigma }_{p}}\text{=}\text{-}\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}{s}^{e}\text{.}\frac{\partial \Delta \lambda }{\partial {\sigma }_{p}}\text{+}I\text{-}\mathrm{3K}(\frac{\partial \eta }{\partial \Delta \lambda }\frac{\partial \Delta \lambda }{\partial {\sigma }_{p}}\Delta \lambda +\eta \frac{\partial \Delta \lambda }{\partial {\sigma }_{p}})I\)

To calculate \(\frac{\partial \mathrm{\Delta \lambda }}{\partial {\mathrm{\sigma }}_{p}}\), we use equation \(\dot{F}(\mathrm{\Delta \lambda })=0\). We get: \(\begin{array}{c}\frac{\partial \Delta \lambda }{\partial {\sigma }_{p}}\left[\begin{array}{c}\frac{3\mu }{{\sigma }_{\text{eq}}^{e}}({s}_{3}^{e}\text{-}{s}_{1}^{e})\text{+}\frac{\partial B}{\partial \Delta \lambda }(1\text{-}\frac{{\sigma }_{3}}{{\sigma }_{3}^{b\text{-}d}})\text{+}3\frac{B}{{\sigma }_{3}^{b\text{-}d}}(\frac{\mu {s}_{3}^{e}}{{\sigma }_{\text{eq}}^{e}}\text{+}K\eta \text{+}K\frac{\partial \eta }{\partial \Delta \lambda }\Delta \lambda )\\ \text{+}(\frac{\partial S{\sigma }_{c}^{2}}{\partial \Delta \lambda }\text{-}\frac{\partial m{\sigma }_{c}}{\partial \Delta \lambda }{\sigma }_{3}\text{+}\mathrm{3m}{\sigma }_{c}\left[\frac{\mu {s}_{3}^{e}}{{\sigma }_{\text{eq}}^{e}}\text{+}K\eta \text{+}K\frac{\partial \eta }{\partial \Delta \lambda }\Delta \lambda \right])\times {(S{\sigma }_{c}^{2}\text{-}{\sigma }_{3}\text{.}m{\sigma }_{c})}^{\text{-}1/2}\end{array}\right]\\ \text{=}\frac{B}{{\sigma }_{3}^{b\text{-}d}}\text{+}\frac{m{\sigma }_{c}}{2\sqrt{S{\sigma }_{c}^{2}\text{-}{\sigma }_{3}\text{.}m{\sigma }_{c}}}\end{array}\)

  1. In the tangent operator it is the derivative \(\frac{\partial \sigma \text{'}}{\partial {\sigma }_{p}}\text{=}\frac{\partial \sigma }{\partial {\sigma }_{p}}\text{-}I\) that should be returned.