降阶均质化方法层级分块计算

发表于 2025-01-12 更新于 2025-05-19 分类于多尺度阅读次数：

在有限元求解框架下的增量步与迭代步

降阶均质化方法是在 ABAQUS UMAT 层面实现的. ABAQUS 在每个增量步的迭代步时调用 UMAT, 考虑有限元计算框架中第 $n$ 个增量步, 第 $i+1$ 个迭代步. 在进入 UMAT 时, 已知前 $n-1$ 个增量步的历史信息, 和当前迭代步 $i+1$ 的位移. 为了使接下来方程的增量表达式更加清晰, 首先引入以下增量记号, 它们写在力学变量之前, 表示不同的含义:

记号 $\delta$ 表示变量该迭代步的校正量;
记号 $\Delta_{n}^{i+1}$ 表示在第 $i+1$ 个迭代步时确定的第 $n$ 个增量步的增量. 若在第 $i+1$ 个迭代步收敛, 那么 $\Delta_{n} :=\Delta_{n}^{i+1}$.

我们再引入花括号 $\{\cdot\}$ 表示分块物理量的集合. 例如 $\{ \Delta \varepsilon \}$ 表示分块上的应变 $\Delta \varepsilon^{(\alpha)}$, $\alpha=1,2,\ldots,N$.

降阶多尺度计算

在接下来的记号中, 使用 $[\cdot]$ 表示应用 Voigt 记法矩阵化的张量; $[:]=\mathrm{diag}(1,1,1,2,2,2)$, $[\boldsymbol{\epsilon}]$ 和 $[\boldsymbol{\varepsilon}]$ 分别表示张量应变和工程应变, $[:][\boldsymbol{\epsilon}]=[\boldsymbol{\varepsilon}]$. 计算在第 $n+1$ 步时各分块上的平均应变分量. 待求解的方程为 \[ \boldsymbol{\epsilon}_{n+1}^{(\beta)} = \mathbb{E}^{(\beta)} : \boldsymbol{\epsilon}_{n+1}^{c} + \sum_{\alpha} \mathbb{P}^{(\beta\alpha)} : \boldsymbol{\mu}_{n+1}^{(\alpha)} \tag{1} \] 写成增量形式, 并定义残差函数为 \[ \boldsymbol{r}^{(\beta)}(\{ \Delta \boldsymbol{\epsilon}^{(\alpha)} \}) \triangleq \Delta \boldsymbol{\epsilon}^{(\beta)} - \mathbb{E}^{(\beta)} : \Delta \boldsymbol{\epsilon}^{c} - \sum_{\alpha} \mathbb{P}^{(\beta\alpha)} : \Delta \boldsymbol{\mu}^{(\alpha)} \] 应用 Voigt 记法, 将上述张量方程写成矩阵形式为 \[ \left[ \boldsymbol{r}^{(\beta)} \right] = \left[ \Delta \boldsymbol{\epsilon}^{(\beta)} \right] - \left[ \mathbb{E}^{(\beta)} \right] \left[ \Delta \boldsymbol{\varepsilon}^{c} \right] - \sum_{\alpha} \left[ \mathbb{P}^{(\beta\alpha)} \right] \left[ \Delta \boldsymbol{\varkappa}^{(\alpha)} \right] \] 式中, $\left[ \Delta \boldsymbol{\varkappa}^{(\alpha)} \right]_{i} = [:]\left[ \Delta \boldsymbol{\mu}^{(\alpha)} \right]_{i}$ 上式对平均应变求偏导数, 就得到 Newton 迭代法中 Jacobian 矩阵的表达式: \[ \frac{ \partial \left[ \boldsymbol{r}^{(\beta)} \right]}{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{(\gamma)} \right]} = \delta^{\beta\gamma} [:]^{-1} - \sum_{\alpha} \left[ \mathbb{P}^{(\beta\alpha)} \right] \frac{ \partial \left[ \Delta \boldsymbol{\varkappa}^{(\alpha)} \right]}{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{(\gamma)} \right]} \] 由于假设分块上的本征应变只与该分块上的平均应变相关, 所以 \[ \frac{ \partial \left[ \boldsymbol{r}^{(\beta)} \right]}{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{(\gamma)} \right]} = \delta^{\beta\gamma} [:]^{-1} - \left[ \mathbb{P}^{(\beta\gamma)} \right] \frac{ \partial \left[ \Delta \boldsymbol{\varkappa}^{(\gamma)} \right]}{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{(\gamma)} \right]} \]

多层级分块状态变量

应用 ROH 作为材料本构进行计算时, 状态变量可分为两类:

各分块非弹性材料演化时的历史状态变量, 称为材料状态变量;
多尺度方程引入的状态变量, 比如分块上的平均应变和平均应力等, 称为多尺度状态变量.

新引入分块层级状态变量 MESFE_P5001226_NSV_NCLUSTER, 并在结构体 ROH_ptn6 中加入新的数据 cidx, 标识分块所属分块组序号. 分块组所对应的分块编号通过外部文件读取. 根据概念大小, 有如下包含关系: \[ \text{Partition} \subseteq \text{PTNGroup} \subset \text{Phase} \]

PTNGRP 状态变量更新

设在当前增量步时, PTNGRP 的数量为为 NPTNGRP(int). 首先更新分块所属分块组的序号, 然后组装分块组, 分块组的数据类型也是是 ROH_ptn6, 接下来降阶多尺度方程的求解在分块组层级上执行. 分块组状态变量应等于分块组所包含的任意单胞分块上的状态变量:

\[ \left[ \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n} := \left[ \boldsymbol{\sigma}^{(\alpha_{i})} \right]^{n}, \quad \left[ \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n} := \left[ \boldsymbol{\varepsilon}^{(\alpha_{i})} \right]^{n}, \quad \left[ \mathbb{E}_{\mathrm{tan}}^{\langle \alpha \rangle} \right]^{n} \triangleq \frac{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n} } := \frac{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{(\alpha_{i})} \right]^{n} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n} } \] 注意, 对于某一分块组所拥有的所有子分块, 它们的材料和多尺度状态变量应该是完全相同的: \[ \left[ \boldsymbol{\sigma}^{(\alpha_{i})} \right]^{n} = \left[ \boldsymbol{\sigma}^{(\alpha_{j})} \right]^{n}, \quad \left[ \boldsymbol{\varepsilon}^{(\alpha_{i})} \right]^{n} = \left[ \boldsymbol{\varepsilon}^{(\alpha_{j})} \right]^{n}, \quad \left[ \mathbb{E}_{\mathrm{tan}}^{(\alpha_{i})} \right]^{n} = \left[ \mathbb{E}_{\mathrm{tan}}^{(\alpha_{j})} \right]^{n} \\ \alpha_{i},\alpha_{j} \in \{ i \mid \text{Partition}(i) \subseteq \text{PTNGRP} \langle \alpha \rangle \} \] 式中, $\langle \alpha \rangle$ 表示分块组的编号, $(\alpha_{i})$ 表示分块的编号. 分块组的材料状态变量初始化为上一次增量步最终收敛得到的值: \[ \left[ \mathrm{statv}^{\langle \alpha \rangle} \right]_{0} := \left[ \mathrm{statv}^{(\alpha_{i})} \right]^{n} \] 应用上一增量步计算得到的瞬时应变影响张量, 根据当前宏观应变增量 $\left[ \Delta \boldsymbol{\varepsilon}^{c} \right]$, 更新各个分块组的初始平均应变增量: \[ \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{1} := \left[ \mathbb{E}_{\mathrm{tan}}^{\langle \alpha \rangle} \right]^{n} \left[ \Delta \boldsymbol{\varepsilon}^{c} \right] \] 接下来应用 Newton 迭代法迭代求解各个分块组的平均应变. 设在第 $i$ 次迭代步时分块组的平均应变增量为 $\left[ \Delta \varepsilon^{\langle \alpha \rangle} \right]_{i}$,代入到各相材料的 UMAT 中进行计算, 得到相应的应力, 一致切线模量和状态变量:

输入: $\left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i}$ $\left[ \mathrm{statv}^{\langle \alpha \rangle} \right]^{n}$
输出: $\left[\boldsymbol{\sigma}^{\langle \alpha \rangle} \right]_{i} := \left[\boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n} + \left[ \Delta \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]_{i}$ $\frac{ \partial \left[ \Delta \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]_{i} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i} }$ $\left[ \mathrm{statv}^{\langle \alpha \rangle} \right]_{i}$

更新本征应变: \[ \left[ \boldsymbol{\varkappa}^{\langle \alpha \rangle} \right]_{i} := \left[ \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n} + \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i} - [\mathbb{L}^{\langle \alpha \rangle}]^{-1} : \left[ \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]_{i} \] 以及本征应变关于分块应变的偏导数 \[ \frac{ \partial \left[ \Delta \boldsymbol{\varkappa}^{\langle \alpha \rangle} \right]_{i} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i} } = [\boldsymbol{I}] - [\mathbb{L}^{\langle \alpha \rangle}]^{-1} \frac{ \partial \left[ \Delta \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]_{i} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i} } \]

由此得到残差和 Jacobian 子矩阵分别为 $$ {i} = {i} - - {} {i}, \

[]{i}^{} { {i}} = ^{} [:]^{-1} - { _{i}} $$

组装多尺度系统的 Jacobian 矩阵 $[\mathbf{A}]_{i}$, 得到分块组的平均应变迭代增量为 \[ \left[ \delta \boldsymbol{\varepsilon} \right]_{i+1} = -[\mathbf{A}]_{i}^{-1} \left[ \boldsymbol{r} \right]_{i} \] 并更新 $i+1$ 步的分块组平均应变:

\[ \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i+1} = \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i} + \left[ \delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]_{i+1} \]

判断迭代收敛之后, 跳出 Newton 迭代过程, 继续判断是否需要分裂分块组. 如果需要分裂分块组, 那么更新分块组的数量 \[ \mathrm{NCLUSTER} := \mathrm{NCLUSTER} + 1 \] 继续重复这一步的过程. 如果不需要分裂分块组的数量, 那么确定最后一次迭代量为当前增量步的增量 $\left[ \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n+1}$ $\frac{ \partial \left[ \Delta \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n+1} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n+1} }$ $\left[ \mathrm{statv}^{\langle \alpha \rangle} \right]^{n+1}$ 并以此计算宏观应力和宏观一致切线模量: \[ \left[ \boldsymbol{\sigma}^{c} \right]^{n+1} = \sum_{\alpha} c^{\langle \alpha \rangle} \left[ \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n+1}, \\ \frac{ \partial \left[ \Delta \boldsymbol{\sigma}^{c} \right]^{n+1} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n+1} } = \sum_{\alpha} c^{\langle \alpha \rangle} \frac{ \partial \left[ \Delta \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n+1} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n+1} } \frac{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n+1} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n+1} } = \sum_{\alpha} c^{\langle \alpha \rangle} \frac{ \partial \left[ \Delta \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n+1} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n+1} } \left[ \mathbb{E}_{\mathrm{tan}}^{\langle \alpha \rangle} \right]^{n+1} \] 式中, $\left[ \mathbb{E}_{\mathrm{tan}}^{\langle \alpha \rangle} \right]^{n+1}$ 为瞬时应变影响张量, 对式 (1) 关于宏观应变增量求偏导数得到 \[ [:]^{-1}\frac{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \beta \rangle} \right]^{n+1} } { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n+1} } = \left[ \mathbb{E}^{\langle \beta \rangle} \right] + \sum_{\gamma} \left[ \mathbb{P}^{\langle \beta\gamma \rangle} \right] \frac{ \partial \left[ \Delta \boldsymbol{\varkappa}^{\langle \gamma\rangle} \right]^{n+1}} { \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \gamma\rangle} \right]^{n+1}} \frac{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \gamma\rangle} \right]^{n+1}} { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n+1}} \] 移项之后, 得到 \[ \sum_{\gamma} [\mathbf{A}]_{i}^{\langle \beta \gamma \rangle} \frac{ \partial \left[ \Delta \boldsymbol{\varepsilon}^{\langle \gamma\rangle} \right]^{n+1}} { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n+1}} = \left[ \mathbb{E}^{\langle \beta \rangle} \right] \] 所以 \[ \frac{ \partial \left[ \Delta \boldsymbol{\varepsilon} \right]^{n+1}} { \partial \left[ \Delta \boldsymbol{\varepsilon}^{c} \right]^{n+1}} = [\mathbf{A}]_{i}^{-1} \left[ \mathbb{E} \right] \] 最后再根据分块组更新后的状态变量, 更新所属分块的状态变量: \[ \left[ \boldsymbol{\sigma}^{(\alpha_{i})} \right]^{n+1} := \left[ \boldsymbol{\sigma}^{\langle \alpha \rangle} \right]^{n+1}, \quad \left[ \boldsymbol{\varepsilon}^{(\alpha_{i})} \right]^{n+1} := \left[ \boldsymbol{\varepsilon}^{\langle \alpha \rangle} \right]^{n+1}, \quad \left[ \mathbb{E}_{\mathrm{tan}}^{(\alpha_{i})} \right]^{n+1} := \left[ \mathbb{E}_{\mathrm{tan}}^{\langle \alpha \rangle} \right]^{n+1} \]