应力回归算法-Forward Euler method

在塑性力学中，由于材料应力状态的复杂性，单一屈服面无法完整体现材料的应力特性，因此多个屈服面组成的本构模型常出现在岩土或者各向异性材料的模拟中。

其中最经典的莫过于摩尔库伦本构，但摩尔库伦本构在deviatoric平面内各个屈服面之间的夹角是确定的，且各个屈服面是平直的，因此可以采用各种特殊的手段简化应力回归。
对于考虑了材料受压破坏的本构模型，控制受压破坏的屈服面和控制拉-剪破坏的屈服面之间的夹角是任意的，而且各个屈服面的表面是光滑(凸)弯曲的，应力回归的过程会比较复杂。

前欧拉法

针对受压屈服面和受剪屈服面,可以给出屈服面方程：

$$\begin{aligned} &F _c(\sigma, k _c) = 0\\ &F _s(\sigma, k _s) = 0 \end{aligned}$$

式中为应力状态，和为各屈服面对应的内变量(等效应变、等效功等)。

两个屈服面可以将应力空间划分为4个部分(弹性空间E、受压破坏区C、受剪破坏区S、压剪破坏重合区A)。
理论上，C、S、A区都是不存在的(试应力点位于C区应力就需要回归到屈服面，试应力点位于S区应力就需要回归到屈服面，试应力点位于A区应力就需要回归到和交点)。

在和交点处，塑性应变可以按下式计算：

$$d\epsilon^p = d\lambda_s \frac{\partial Q_s}{\partial \sigma} + d\lambda_c \frac{\partial Q_c}{\partial \sigma}$$

应力增量可按下式计算：

$$d\sigma = [K](d\epsilon - d\epsilon^p)$$

对屈服面，将其对(伪)时间求导(consistency condition)：

$$\begin{aligned} dF_s &= \frac{\partial F_s}{\partial \sigma} d\sigma + \frac{\partial F_s}{\partial k_s} d k_s\\ &= \frac{\partial F_s}{\partial \sigma} d\sigma + \frac{\partial F_s}{\partial k_s} \frac{\partial k_s}{\partial \epsilon^p} d\epsilon^p\\ &= \frac{\partial F_s}{\partial \sigma} d\sigma + \frac{\partial F_s}{\partial k_s} \frac{\partial k_s}{\partial \epsilon^p} (d\lambda_s \frac{\partial Q_s}{\partial \sigma} + d\lambda_c \frac{\partial Q_c}{\partial \sigma})\\ &= \frac{\partial F_s}{\partial \sigma} d\sigma + d\lambda_s \frac{\partial F_s}{\partial k_s} \frac{\partial k_s}{\partial \epsilon^p} \frac{\partial Q_s}{\partial \sigma} + d\lambda_c \frac{\partial F_s}{\partial k_s} \frac{\partial k_s}{\partial \epsilon^p} \frac{\partial Q_c}{\partial \sigma} = 0 \end{aligned}$$

令

$$\begin{aligned} &H_{ss} = - \frac{\partial F_s}{\partial k_s} \frac{\partial k_s}{\partial \epsilon^p} \frac{\partial Q_s}{\partial \sigma}\\ &H_{sc} = - \frac{\partial F_s}{\partial k_s} \frac{\partial k_s}{\partial \epsilon^p} \frac{\partial Q_c}{\partial \sigma} \end{aligned}$$

则有

$$\frac{\partial F_s}{\partial \sigma} d\sigma - d\lambda_s H_{ss} - d\lambda_c H_{sc} = \frac{\partial F_s}{\partial \sigma} [K](d\epsilon - d\epsilon^p) - d\lambda_s H_{ss} - d\lambda_c H_{sc} = 0$$

即

$$\begin{aligned} &\frac{\partial F_s}{\partial \sigma} [K] d\epsilon =\frac{\partial F_s}{\partial \sigma} [K] d\epsilon^p + d\lambda_s H_{ss} + d\lambda_c H_{sc}\\ &\frac{\partial F_s}{\partial \sigma} [K] d\epsilon =\frac{\partial F_s}{\partial \sigma} [K] d\lambda_s \frac{\partial Q_s}{\partial \sigma} + \frac{\partial F_s}{\partial \sigma} [K] d\lambda_c \frac{\partial Q_c}{\partial \sigma} + d\lambda_s H_{ss} + d\lambda_c H_{sc} \end{aligned}$$

对屈服面，将其对(伪)时间求导(consistency condition)：

$$\begin{aligned} dF_c &= \frac{\partial F_c}{\partial \sigma} d\sigma + \frac{\partial F_c}{\partial k_c} d k_c\\ &= \frac{\partial F_c}{\partial \sigma} d\sigma + \frac{\partial F_c}{\partial k_c} \frac{\partial k_c}{\partial \epsilon^p} d\epsilon^p\\ &= \frac{\partial F_c}{\partial \sigma} d\sigma + \frac{\partial F_c}{\partial k_c} \frac{\partial k_c}{\partial \epsilon^p} (d\lambda_s \frac{\partial Q_s}{\partial \sigma} + d\lambda_c \frac{\partial Q_c}{\partial \sigma})\\ &= \frac{\partial F_c}{\partial \sigma} d\sigma + d\lambda_s \frac{\partial F_c}{\partial k_c} \frac{\partial k_c}{\partial \epsilon^p} \frac{\partial Q_s}{\partial \sigma} + d\lambda_c \frac{\partial F_c}{\partial k_c} \frac{\partial k_c}{\partial \epsilon^p} \frac{\partial Q_c}{\partial \sigma} = 0 \end{aligned}$$

令

$$\begin{aligned} &H_{cs} = - \frac{\partial F_c}{\partial k_c} \frac{\partial k_c}{\partial \epsilon^p} \frac{\partial Q_s}{\partial \sigma}\\ &H_{cc} = - \frac{\partial F_c}{\partial k_c} \frac{\partial k_c}{\partial \epsilon^p} \frac{\partial Q_c}{\partial \sigma} \end{aligned}$$

则有

$$\frac{\partial F_c}{\partial \sigma} d\sigma - d\lambda_s H_{cs} - d\lambda_c H_{cc} = \frac{\partial F_c}{\partial \sigma} [K](d\epsilon - d\epsilon^p) - d\lambda_s H_{cs} - d\lambda_c H_{cc} = 0$$

即

$$\begin{aligned} &\frac{\partial F_c}{\partial \sigma} [K] d\epsilon =\frac{\partial F_c}{\partial \sigma} [K] d\epsilon^p + d\lambda_s H_{cs} + d\lambda_c H_{cc}\\ &\frac{\partial F_c}{\partial \sigma} [K] d\epsilon =\frac{\partial F_c}{\partial \sigma} [K] d\lambda_s \frac{\partial Q_s}{\partial \sigma} + \frac{\partial F_c}{\partial \sigma} [K] d\lambda_c \frac{\partial Q_c}{\partial \sigma} + d\lambda_s H_{cs} + d\lambda_c H_{cc} \end{aligned}$$

综上，可得矩阵形式

$$\begin{bmatrix} \frac{\partial F_s}{\partial \sigma} [K] d\epsilon\\\\ \frac{\partial F_c}{\partial \sigma} [K] d\epsilon \end{bmatrix} = \begin{bmatrix} \frac{\partial F_s}{\partial \sigma} [K] \frac{\partial Q_s}{\partial \sigma} + H_{ss} & \frac{\partial F_s}{\partial \sigma} [K] \frac{\partial Q_c}{\partial \sigma} + H_{sc}\\\\ \frac{\partial F_c}{\partial \sigma} [K] \frac{\partial Q_s}{\partial \sigma} + H_{cs} & \frac{\partial F_c}{\partial \sigma} [K] \frac{\partial Q_c}{\partial \sigma} + H_{cc} \end{bmatrix} \begin{bmatrix} d \lambda_s\\\\ d \lambda_c \end{bmatrix}$$

令

$$[A] = \begin{bmatrix} \frac{\partial F_s}{\partial \sigma} [K] \frac{\partial Q_s}{\partial \sigma} + H_{ss} & \frac{\partial F_s}{\partial \sigma} [K] \frac{\partial Q_c}{\partial \sigma} + H_{sc}\\\\ \frac{\partial F_c}{\partial \sigma} [K] \frac{\partial Q_s}{\partial \sigma} + H_{cs} & \frac{\partial F_c}{\partial \sigma} [K] \frac{\partial Q_c}{\partial \sigma} + H_{cc} \end{bmatrix}$$

则屈服面塑性增量因子可通过上式求逆计算

$$\begin{bmatrix} d \lambda_s\\\\ d \lambda_c \end{bmatrix} = [A]^{-1} \begin{bmatrix} \frac{\partial F_s}{\partial \sigma} [K] d\epsilon\\\\ \frac{\partial F_c}{\partial \sigma} [K] d\epsilon \end{bmatrix}$$

具体的算法思路如下：