Axisymmetric SNS-FEM formulation

本文讨论小变形轴对称 THM 的 SNS-FEM 离散，包括光滑域几何、四个虚拟积分点、耦合残差与切线，以及机械体积力和 Darcy 重力项。公式按本文的符号与章节独立说明；验证记录与检查要求集中在第 14 节。

0. 物理范围与符号约定

本文采用以下假设：小应变、准静态、固定参考几何、单一均匀且各向同性的材料、不可压缩固体颗粒 $K_s\to\infty$，以及标量 Biot 系数 $b=1$。孔隙率及水力、热学系数取给定参考值，不随位置或当前未知量更新。机械体力允许是已知的空间和时间函数，但不依赖当前 Newton 未知量；Darcy 定律中的重力加速度采用给定常向量。含一般体力不意味着材料非均匀，也不意味着可把机械体力自动解释成作用于流体的加速度。

体力、重力加速度和混合物参考密度分别写为

其中 $\boldsymbol f$ 是单位总体积的机械体力，量纲为力/体积；$\boldsymbol g$ 是重力加速度，量纲为长度/时间平方。纯重力时 $\boldsymbol f=\rho\boldsymbol g$；$z$ 轴向上时 $\boldsymbol g=[0,-g]^T$。

注意： 若还含有其他机械体力，应将它加到 $\boldsymbol f$ 中，但不得因此把 Darcy 定律中的 $\boldsymbol g$ 替换为 $\boldsymbol f/\rho$。后者会改变模型中流体受到的作用，不能仅由力学平衡方程推出。

本文使用以下物理符号：

在 $b=1$ 的设定下，孔压耦合直接使用几何列向量 $\boldsymbol m$，不再引入独立的张量型 Biot 参数：

材料参数不带光滑域下标 $k$。 例如 $n,K_f,\alpha_s,\alpha_f,\rho_s,\rho_f,c_s,c_f,k_f,\gamma_f,k_T,S,\beta_T,C_T$ 在整个计算域内相同。几何量 $V_k$、局部矩阵 $K_{ab}^k$、节点向量 $\boldsymbol U_k$ 及积分点状态 $\boldsymbol\sigma'_{kg},C_{kg}^{\mathrm{ep}}$ 仍保留位置索引。均匀材料的各积分点可以经历不同加载路径，因此应力、塑性变量和算法切线未必相同。

本构的四分量应力顺序为 $[rr,zz,\theta\theta,rz]$，应变的第 4 项使用工程剪应变 $\gamma_{rz}=2\varepsilon_{rz}$。矩阵分量排列与张量所选基的排列必须明确对应，不能混用剪切应变因子。

水力传导系数 $k_f$ 与导热系数 $k_T$ 均为常标量。若程序接口需要二维矩阵，则分别传入 $k_f I_2$ 与 $k_T I_2$；$I_2$ 是二维单位矩阵，正文公式直接使用标量。以下 Darcy 重力采用给定竖向常量 $\boldsymbol g=[0,-g]^T$，机械一般体力则允许是已知的 $\boldsymbol f(\boldsymbol x,t)$。

1. 已确定的前提与光滑方案

采用小变形轴对称模型，不考虑绕轴扭转。材料采用第 0 节的均匀、各向同性及 $b=1$ 设定。$r,z$ 分别为径向、轴向坐标，$r\geq0$，$\theta$ 表示环向角度。$u,p,T$ 分别表示位移、孔隙水压力和温度；$B$ 是将节点位移换算为应变的矩阵。位移、孔压和温度整体耦合求解：变形与孔压相互影响，温度变化影响变形和孔压。热对流可通过 Keyword.HeatAdvection 开关选择；第 8 节统一给出含重力的通量、残差和切线。 关闭时，温度方程不依赖位移和孔压；开启时，孔压通过 Darcy 速度影响温度，但仍不加入由变形率、孔压变化率直接产生的热量项，具体见第 8.5 节。实现时逐场检查绝对和相对残差，避免零载荷 patch 的零除；必要时将 Newton 更新减半重试，每次本构试算仍使用上一步已收敛状态。

已确定：所有光滑操作均按四边形子域的实际旋转体积加权。 计算按“平均矩阵（第 4 节）→ 梯度矩阵（第 5 节）→ 虚拟点矩阵（第 6 节）”的顺序展开。

下表汇总本文采用的构造。

2. 先看清光滑域，再区分面积、体积和组装

2.1 一个完整光滑域由哪些部分组成

先选一个网格节点 $k$。在每个与它相连的三角形中，依次连接“节点 $k$—一侧边的中点—三角形的面积形心—另一侧边的中点”，得到一个四边形。

$D_{ke}$ 是三角形 $e$ 内的这一个四边形；$D_k$ 是围绕节点 $k$ 的所有这些四边形合在一起的区域。 下文把三角形 $e$ 称为“母三角形”，只是说明 $D_{ke}$ 属于哪个单元。

图 1。 左图中，所有蓝色四边形合起来才是 $D_k$。右图只放大其中一个三角形：蓝色四边形是 $D_{ke}$，其余部分属于该三角形另外两个顶点的光滑域。橙色轮廓表示同一个母三角形。一个完整光滑域通常跨越多个三角形，一个子域则只占其中一个三角形的一部分。计算光滑域的矩阵时，要用到这些母三角形的全部节点，包括图中位于蓝色区域外的三角形顶点。

2.2 面积、旋转体积与虚拟点距离

$A_{ke}$ 是四边形 $D_{ke}$ 在 r–z 平面内的面积；$V_{ke}$ 是它绕 z 轴旋转一周后的体积。$A_k,V_k$ 则分别表示完整光滑域的面积、旋转体积。因此，

这里 $e\ni k$ 表示“遍历所有包含节点 $k$ 的三角形”。面积和体积都是标量，直接相加即可。

$r$ 是到旋转轴的距离，$A_k$ 是二维面积，$V_k$ 是旋转一周后的体积，$\ell_k$ 是构造四个虚拟点矩阵时使用的偏移距离。式中的 $\mathrm{d}s$ 是平面内的一小段边长，$\mathrm{d}S$ 是这段边旋转后对应的表面积。当前轴对称代码用 virtualPointDistance = sqrt(domain.area/pi) 计算 $\ell_k$，与径向坐标 $r$ 区分。

轴对称内力、刚度等要按体积积分，但不能把代码中每一个面积都换成体积。 例如，$\ell_k=\sqrt{A_k/\pi}$ 表示与光滑域面积相等的圆的半径；如果改为 $\sqrt{V_k/\pi}$，结果的单位就不再是长度。

本文在体积和表面积分中统一保留 $2\pi$：长度用 m、应力用 kPa 时，等效节点力就是 kN。基本应变公式和半径在积分中的位置可参见 FEniCSx 轴对称弹性示例；下面有关 SNS 离散的讨论按本笔记中的几何与组装定义展开。

2.3 组装：按节点编号对齐，再相加

本文统一用 $\operatorname{Assemble}$ 表示组装。它的意思就是：把对应节点自由度的行或列，按节点编号和自由度类型对齐，再相加。 刚度矩阵需要对齐行和列；B、N、E 只需对齐节点列，它们的应变、标量或梯度分量行保持原顺序。位移的径向与轴向分量、孔压、温度分别对应各自的位置。

例如，只看每个节点的一个标量贡献（数值仅用于示意）。三角形 1 的节点依次为 4、8、12，三角形 2 的节点依次为 4、12、19。按节点编号对应后，组装结果是：

这些数字只用于说明组装。三角形 1 的第 3 个节点和三角形 2 的第 2 个节点都是节点 12，所以对应的 7 和 11 要相加。刚度矩阵组装时，行和列都按这个规则对应。

图 1 说明的是“哪些四边形组成一个光滑域”；组装说明的是“这些单元或光滑域的矩阵怎样对应相加”。

2.4 后文共用的记号

下面几种记号会反复出现：

矩阵块 $K_{ab}^{k}$ 的上标 $k$ 表示光滑域 $k$ 的局部贡献；不带这个上标的 $K_{ab}$ 表示组装后的全局矩阵。下标 $a,b$ 标记方程和未知量，具体对应见第 8 节。正文中的向量用粗体表示；B、N、E 等矩阵则保留与代码接近的字母写法。

本文的试函数（test function）统一使用字母 v，用下标区分方程、粗体区分向量与标量。 力学方程的试函数为虚位移向量 $\boldsymbol v_u=[v_r,v_z]^T$，其中 $v_r,v_z$ 分别是径向、轴向分量；孔压方程的试函数为标量 $v_p$；温度方程的试函数为标量 $v_T$。下标 $u,p,T$ 标明所属方程，各方程中的试函数独立选取。每个弱式前都会说明其含义。

偏导数与标签分开写。 为避免歧义，本文将空间偏导直接展开，例如 $\frac{\partial N_I}{\partial r}$、$\frac{\partial p}{\partial r}$，不用下标逗号表示求导。下列已命名离散量中保留的逗号只分隔标签或索引，不表示再求一次导数：

上标中的 $V,\mathrm c$、$V,\mathrm{incomplete}$、$\mathrm{adv,grav}$，以及 $K_{33\mathrm{Adv,grav}}$ 中的逗号，也只分隔方案或贡献标签。例如 $n_{\mathrm{out},r}$ 表示外法向的径向分量，其含义与偏导运算不同。函数自变量、向量分量和并列公式之间的普通逗号仍按通常标点理解。

3. 四边形光滑子域的体积如何计算

3.1 三角形的体积可以用三角形中心计算

母三角形顶点为 $i,j,l$，面积为 $A_e$；$r_i,r_j,r_l$ 是三个顶点的半径，$r_{ec}$ 是三角形面积形心的半径，$V_e$ 是该三角形旋转一周的体积。因为 $r$ 在直边三角形上是线性函数，所以

是精确几何公式。

3.2 子域体积一般不是三角形体积的三分之一

本节用 $\boldsymbol x_{ie}^{A}=[r_{ie}^{A},z_{ie}^{A}]^T$ 表示四边形子域 $D_{ie}$ 的面积形心坐标向量，也就是把整个四边形看成均匀薄片时的重心坐标。下标 $i$ 是子域所对应的节点编号，$e$ 是母三角形编号；上标 $A$ 表示按面积求平均。计算子域实际旋转体积时，使用的是面积形心的径向坐标 $r_{ie}^{A}$。面积形心一般不同于四个顶点坐标的算术平均；后者在第 3.3 节记为 $\boldsymbol x_{ie}^{\mathrm{vertices}}$。

节点 $i$ 对应的四边形依次连接“节点 $i$—边 $ij$ 中点—三角形中心—边 $il$ 中点”。它的面积确实是 $A_e/3$，面积形心的坐标为：

这里 $\boldsymbol x_{ie}^{A}$ 是子域 $D_{ie}$ 的面积形心坐标，$A_{ie}=A_e/3$ 是子域面积，$V_{ie}$ 是其真实旋转体积。因此

将四边形沿“节点—三角形中心”分成两个等面积小三角形后分别积分，即可得到这个结果。它满足

但一般不满足 $V_{ie}=V_e/3$。对子域 $i$，等式 $V_{ie}=V_e/3$ 成立的条件是 $2r_i=r_j+r_l$。

例如，取三角形的三个顶点为 $(r,z)=(0,0),(1,0),(0,1)$，按此顺序将顶点局部编号为 $1,2,3$。三个子域面积相等，但旋转体积分配为

图 2。 左图按实际坐标画出三个子域，右图用相同颜色比较面积占比和体积占比。这里的分数只对应这个三角形，不是所有网格的固定分配比例。

远离旋转轴的子域，转一圈后扫过的体积更大。另外，轴线节点虽然 $r_i=0$，但它的光滑域仍向轴线外侧延伸，旋转后的体积不为零。

旧手稿中的子域体积计算错误： 旧手稿采用 $2\pi r_{ec}A_{ke}$，相当于把母三角形的旋转体积平均分成三份，一般不等于各四边形子域的实际体积。本文按子域面积形心的半径计算实际旋转体积，公式见上文。

3.3 四个顶点的平均坐标不等于面积形心（注意）

上述四边形的顶点平均坐标记为 $\boldsymbol x_{ie}^{\mathrm{vertices}}$，其中上标 $\mathrm{vertices}$ 只标明“顶点平均”：

四个顶点的平均坐标 $\boldsymbol x_{ie}^{\mathrm{vertices}}$ 一般不等于第 3.2 节定义的面积形心坐标 $\boldsymbol x_{ie}^{A}=[r_{ie}^{A},z_{ie}^{A}]^T$。这里 $i$ 是子域所对应的节点编号，$e$ 是母三角形编号，上标 $A$ 表示按面积求平均。$\boldsymbol x_{ie}^{A}$ 是把整个四边形子域 $D_{ie}$ 看成均匀薄片时的重心坐标，其中 $r_{ie}^{A}$ 是计算该子域实际旋转体积所用的半径。

原 markdown 文件（本讨论整理前的技术笔记）记载，参考函数 Smooth_T3_axisym.m 使用 xc_sub=(x1+xm1+xc+xm3)/4 计算子域中心坐标，并使用 N_ino/center_SD_ino 计算 B 的环向应变行。按这条四点算术平均表达式，代码变量 xc_sub 对应本文的 $\boldsymbol x_{ie}^{\mathrm{vertices}}$，即四个顶点坐标的平均，一般不等于面积形心 $\boldsymbol x_{ie}^{A}$。下面分析这两条表达式对环向应变的影响。

按照原 markdown 文件对 N_ino/center_SD_ino 的描述，若分子取型函数的精确面积平均，分母取四个顶点的平均半径，则两者不对应。对第 3.2 节的单三角形、节点 1 和 $u_r=r$，由此得到 $\varepsilon_{\theta\theta}=14/15$，正确值应为 $1$。这一结果说明，混用型函数的面积平均与四顶点平均半径会导致环向应变误差。旧手稿中的“子域形心”采用的就是四顶点平均坐标；将其当作四边形子域的面积形心是错误的。

4. 平均 B 矩阵的三种构造

方法一、二都先用母三角形中心矩阵 $B_e$ 近似单元内的 B，再分别按面积、旋转体积平均。方法三直接对原始 B 做精确旋转体积平均，本文采用方法三。

4.1 方法一：中心矩阵的面积平均（不采用）

对母单元中心矩阵做面积平均，会产生轴对称不平衡内力，因此不采用。 下面给出构造及反例。

母三角形 $e$ 的中心半径为 $r_{ec}$，中心处的型函数值为 $1/3$。节点 $I$ 的径向、轴向位移对应的两列为

$B_e$ 由三个节点的列块组成，应变分量顺序为 $[rr,zz,\theta\theta,rz]$，第 4 项为工程剪应变。对各四边形子域按面积加权，得到

反例： 四个角点为 $(1,0),(3,0),(3,2),(1,2)$，内部点为 $(1.7,0.8)$；将内部点与四角连接，得到四个三角形。长度单位为 m，孔压、体力均为零，施加均匀应力

外边界施加与该应力对应的面力，内部节点没有外荷载，其合内力应为零。孔压为零时，有效应力 $\boldsymbol\sigma'$ 等于总应力；各光滑域的骨架内力为

组装后，内部节点的径向内力为 $3.4867281154\ \mathrm{kN}$，正确值应为零。这就是该面积平均方法不满足力平衡的反例。

4.2 方法二：中心矩阵的旋转体积平均

仍用同一套 $B_e$，改用各子域的实际旋转体积加权：

上标 $\mathrm c$ 表示使用母单元中心值。该方法保持组装后的 B 体积积分：

对每个母单元，中心环向值满足

其余三行的 T3 型函数导数为常数，因此组装后的积分与原始 B 的体积积分相同。应力处处相同时，该方法与原始轴对称 T3 得到相同的节点内力。 在上述反例中，内部节点径向内力为 $-1.8\times10^{-15}\ \mathrm{kN}$，即舍入误差范围内的零。

这里保持的是组装后的积分；单个四边形子域内的环向值仍采用母单元中心近似。

4.3 方法三：原始 B 的精确旋转体积平均（本文采用）

直接积分原始位置相关矩阵 $B(r,z)$：

B 的第 3 行计算环向应变 $\varepsilon_{\theta\theta}=u_r/r$。 其径向位移列的原始系数为 $N_I/r$，精确平均后记为 $\overline h_{kI}^V$：

这里用 $\mathrm{d}V=2\pi r\,\mathrm{d}A$ 抵消了 $1/r$。右端虽写成面积积分，除数仍为真实体积 $V_k$，所以仍是旋转体积平均。节点 $I$ 不属于某个母单元时，该单元对该系数的贡献为零。

对本文的“节点—边中点—单元中心—另一边中点”四边形，节点 $k,j,l$ 按母三角形顺序排列，有

平均矩阵第 3 行的径向位移列取 $\overline h_{kI}^V$，轴向位移列取零。第 1、2、4 行的型函数导数为常数，其精确体积平均与方法二相同。

图 3。 第 1 步采用方法三，第 4 步使用每点 $V_k/4$ 的体积权重。有效应力由本构更新得到，输入的位移、温度增量均从上一步已收敛状态算起。四个虚拟点的骨架内力分别计算后组装，孔压内力另行组装。

5. 光滑梯度的构造

第 4 节已给出平均矩阵 $\overline B_k^V$；本节构造径向、轴向两个梯度矩阵。它们是分别计算的离散量，最后在第 6.1 节组合成虚拟点矩阵。下面依次说明梯度的记号、输入矩阵和计算公式。

5.1 梯度的记号与节点列排列

径向、轴向梯度矩阵分别记为 $G_{B,r,k}^V$、$G_{B,z,k}^V$，上标 $V$ 表示体积加权。下标依次标明矩阵种类 B、方向 r 或 z、光滑域编号 $k$；逗号只分隔标签，不表示求导。后面的 $G_{E,r,k}^V,G_{E,z,k}^V$ 也采用这一顺序。

$\widehat B$ 表示梯度积分所使用的母单元矩阵，其取值在第 5.2 节给出。积分前，各母单元矩阵的列都按光滑域 $k$ 的节点顺序排列；该单元不含的节点对应列补零。

$\partial D_k$ 是光滑域边界，$\boldsymbol n_{\mathrm{out}}=[n_{\mathrm{out},r},n_{\mathrm{out},z}]^T$ 是其单位外法向，两个分量分别对应径向、轴向。

5.2 梯度计算所用的母单元中心矩阵

本文在每个母三角形内取 $\widehat B=B_e$，其中 $B_e$ 是第 4.1 节的母单元中心矩阵，并按第 5.1 节对齐节点列。 对属于该三角形的节点 $I$，中心点 $(r_{ec},z_{ec})$ 处有

因此，$\widehat B$ 的第 3 行径向位移列取 $1/(3r_{ec})$，轴向位移列取零；其余行使用该 T3 单元的常数型函数导数。

每个母三角形内部的 $\widehat B$ 是常数，相邻三角形的取值一般不同。第 5.3 节的公式计入这些跨单元的差别。不能仅因各母单元内的矩阵为常数，就把光滑梯度设为零。

5.3 体积加权梯度与径向修正

先用一般矩阵场 $B$ 推导，暂不指定其离散取值。在下列积分存在时，径向光滑梯度定义为 $\partial B/\partial r$ 的体积平均。

乘积法则：

将它代入体积平均，再对每个矩阵分量使用二维散度定理，得到

这就是径向修正项的来源。轴向推导相同，但 $r$ 不随 $z$ 变化，因此没有对应的域内修正项。对跨单元分片不连续的 $B$，这里的导数按包含单元间跳跃的弱导数理解。

注意： 若直接代入原始环向项，域内积分含有 $\int N_I/r\,\mathrm{d}A$。当子域沿轴线有一段边界，且 $N_I$ 在该段非零时，这个积分对数发散，没有有限值。

离散计算时，将边界项和域内项中的 $B$ 同时替换为第 5.2 节的母单元中心矩阵 $\widehat B$。其中的 $1/(3r_{ec})$ 有限，可避免上述奇异性，得到

注意： 边界积分和域内积分必须使用第 5.2 节给定的同一套 $\widehat B$。域内积分按 $\operatorname*{Assemble}_{e\ni k}\left[2\pi A_{ke}B_e\right]$ 组装，不能替换成 $2\pi A_k\overline B_k^V$。

均匀应变再现检查： 给节点施加位移场 $u_r=\varepsilon_{r0}r$、$u_z=\varepsilon_{z0}z+u_{z0}$，其中 $\varepsilon_{r0},\varepsilon_{z0},u_{z0}$ 均为常数，应得到恒定应变

此时 $\widehat B\boldsymbol U_k=\boldsymbol\varepsilon_0$ 在光滑域内处处相同，径向、轴向应变梯度都应为零。以非零的 $\boldsymbol\varepsilon_0$ 为例，省略径向公式的第二项却会得到

这里先代入 $\widehat B\boldsymbol U_k=\boldsymbol\varepsilon_0$，将常量应变提出积分；再用二维散度定理将边界积分转为面积积分，并利用 $\partial r/\partial r=1$。上标 $\mathrm{incomplete}$ 表示漏掉径向第二项的计算结果。也就是说，本来没有变化的应变会被算成有径向变化，四个虚拟点的应力和刚度因此受到影响。加上第二项后，结果才为零。

使用完整公式后，两个应变梯度均为零；第 4.3 节的平均矩阵也给出 $\boldsymbol\varepsilon_0$。因此，第 6.1 节由平均矩阵加减梯度项构造的四个虚拟点，都得到同一个正确的均匀应变。

压力、温度所用的 E 矩阵及其体积加权梯度见第 7 节。

6. 四个虚拟点、V/4 与真实体积积分

6.1 四个虚拟点及其体积权重

按已确定的体积加权方案，使用第 4.3 节的平均矩阵和第 5.3 节的近似梯度，四个虚拟点的 B 矩阵写为：

$W_{kg}$ 是第 $g$ 个虚拟点的体积积分权重，单位是体积。点 1、2、3、4 依次对应径向负、径向正、轴向负、轴向正的偏移。四个点的权重都为正，四个 B 相加除以四，正好得到平均的 NS 矩阵。下文的 $r_k$ 表示光滑域 $k$ 的中心节点到旋转轴的距离，$r_g$ 表示示意点的径向坐标（省略所属光滑域 $k$）。每个虚拟点分别保存一套塑性状态变量。

用了 $V_k/4$，就不能再乘 $2\pi r_g$。 因为计算体积时已经乘过半径，再乘就是重复计算。本文的权重也不能用 $2\pi r_gA_k/4$ 替代：若四个点围绕节点对称布置，四个权重之和是 $2\pi r_kA_k$，通常不等于光滑域的实际体积 $V_k$，轴线附近还可能出现零权重或负权重。

6.2 “四点乘体积”并不自动等于真实域上的精确积分

如果希望四个点精确计算随位置变化的 B 所产生的积分，仅让四个权重之和等于 $V_k$ 还不够。例如，计算刚度时，B 的分量会两两相乘，会出现位置偏移量的平方，因此还要检查这些平方项的积分是否正确。即使只看面积为 $A_k$、半径为 $\ell_k$ 的均匀二维圆盘，四个等权点放在圆周上时，径向偏移量的平方平均值也是（$\delta r_g$ 是第 $g$ 点相对圆盘中心的径向坐标差，$\delta r$ 是盘内一般位置的同类坐标差）：

真实圆盘区域在积分下标中记为 $\mathrm{disk}$，它满足

四点规则给出的二阶矩是上述圆盘精确值的两倍，因此它不精确积分该圆盘的二阶矩。$\ell_k=\sqrt{A_k/\pi}$ 配合四个等权点，规定了 SNS 梯度项对刚度的贡献大小。保留这一布点与权重，也就保留了原 SNS 的梯度稳定化构造。

本文保留 $\ell_k=\sqrt{A_k/\pi}$ 的点距和每点 $V_k/4$ 的权重。

6.3 靠轴线的“负半径虚拟点”如何理解

本文采用位于 $r\geq0$ 区域的非退化 T3 网格，因此母单元面积 $A_e>0$、中心半径 $r_{ec}>0$，光滑域实际旋转体积 $V_k>0$。预计算 B 及其梯度时，分母来自这些真实几何量。虚拟点的 B 直接按第 6.1 节的四个公式组合，无须在虚拟点坐标处重新求型函数。即使示意图上出现 $r_k-\ell_k<0$，这个示意半径也不进入 $N/r$ 求值或积分权重；内力计算使用组合后的矩阵和正权重 $V_k/4$。

注意： 虚拟点矩阵按第 6.1 节直接组合。不要在负半径处重新求型函数，也不能将示意半径改成 $\lvert r_g\rvert$ 后用于求值或积分。

7. 温度应变，以及压力、温度梯度如何处理

7.1 保留 SNS 的四点 E 矩阵

本文用 SNS 虚拟点矩阵计算渗流和热传导。为区分不同的 E，本文记原始 T3 型函数梯度为 $E_e=\nabla N_e$，常数 NS 矩阵为 $\overline E_k$，SNS 虚拟点矩阵为 $\widetilde E_{kg}$；其中 $k,g$ 分别标明光滑域和虚拟点。四点矩阵由平均 E 加减相应方向的梯度项得到。

热传导使用四点 $\widetilde E_{kg}$，热容量使用原始 T3 型函数。关闭热对流时，温度方程不依赖 $u,p$；开启后，仍用四点 $\widetilde E_{kg}$ 计算 Darcy 速度和温度梯度，新增项见第 8.5 节。温度方程不额外加入 PPP 项，E 中的 SNS 梯度项保留。节点光滑在瞬态导热中也可能出现积分不足造成的数值问题，保留这些梯度项有依据，见 Cui 等关于稳态与瞬态导热的研究。

7.2 体积平均的温度与其梯度

采用前面的体积方案，统一定义

T3 的 $\nabla N_e$ 在母单元内为常数，后一式可以直接用子域体积加权。注意： $\overline E_k^V$ 是“先求导、再平均”，不应当对已经存为常数的 $\overline N_k^V$ 直接求导而得到零。

这里 $\nabla=(\partial_r,\partial_z)^T$，$\partial_r,\partial_z$ 分别表示对 $r,z$ 求偏导。$\overline E_{r,k}^V$、$\overline E_{z,k}^V$ 分别是平均梯度矩阵 $\overline E_k^V$ 的第一行、第二行。继续定义

$\alpha_s$ 是固体的线膨胀系数；$\boldsymbol m=[1,1,1,0]^T$ 将标量热膨胀量放入三个正应变分量，并使剪应变为零。$\Delta\boldsymbol\varepsilon_{kg}^{\mathrm{th}}$ 是第 $g$ 个虚拟点的热应变增量，上标 $\mathrm{th}$ 表示热应变。这里的 $\widetilde N_{kg}$ 将节点温度增量换算成该点用于计算热应变的温度增量。它由平均温度矩阵加减梯度项得到，不需要判断虚拟点落在哪个三角形里，再去求那个三角形的型函数。

均匀温升时必须有 $\overline N_k^V\boldsymbol 1=1$、$\overline E_k^V\boldsymbol 1=\boldsymbol 0$，其中 $\boldsymbol 1$ 是与节点数量对应的全 1 列向量，所以四个虚拟点获得相同热应变。在本文均匀材料假设下，无初始应力、无机械体力、无额外机械约束，且孔压保持 $p=0$ 时，自由膨胀位移为 $u_r=\alpha_s\Delta T r,\ u_z=\alpha_s\Delta T z$（$\Delta T$ 为均匀温升）。此时总应变减去热应变得到的机械应变为零，包括环向分量。注意： 不排水升温可能产生孔压，不能直接套用这一自由膨胀检查。

7.3 压力、温度的梯度没有 N/r 行

压力和温度的梯度仍然只有径向和轴向两个分量。轴对称给扩散方程带来的半径因子，已经包含在体积积分中。以热传导项为例：$D$ 表示整个 r–z 计算域，$k_T$ 是常标量导热系数，$v_T$ 是温度方程弱式中的标量试函数，在规定温度的边界上为零。它不是温度 $T$，也不是时间步温度增量 $\Delta T$：

上式已经考虑圆柱坐标的影响，不能再给热传导或 Darcy 的梯度矩阵添加一个类似机械环向应变的 $N/r$ 行，也不需要在刚度里重复增加 $1/r$ 项。

E 的体积加权梯度同样需要径向修正。 记 $G_{E,r,k}^{V},G_{E,z,k}^{V}$ 为计算四点 E 时使用的径向、轴向变化矩阵。对 E 应用第 5.3 节的乘积法则与散度定理，得到

这里 $\widehat E$ 在每个母三角形内取该单元的 $E_e$；各矩阵的节点列同样先按光滑域 $k$ 的节点顺序排列，缺少的列补零。它不含 $N/r$，而且在直边 T3 内是常数，所以新增的面积积分可以按各子域精确求和。没有 $N/r$，不代表可以省略径向修正：这个修正来自体积权重中的 r。 若积分时的 $\widehat E$ 处处相同，径向公式中的两项恰好抵消，结果为零。

例如前两个虚拟点使用

第 3、4 点分别减去、加上轴向项 $\ell_kG_{E,z,k}^{V}$。这里算的是相邻单元 E 值的空间变化，不是对已经存成常数的 $\overline E_k^V$ 直接求导。热应变中用于构造四点温度的 $\overline E_k^V$ 仍按第 7.2 节使用。

8. 内力、残差与刚度应该怎样一起修改

下面逐块说明单一均匀、各向同性且 $b=1$ 时的公式。矩阵块下标中的 $1,2,3$ 依次对应位移、孔压、温度：第一个数字表示方程，第二个数字表示求导的未知量；关闭对流时各块的实际正负号见第 8.4 节，开启后的新增项见第 8.5 节。$\boldsymbol m=[1,1,1,0]^T$ 是几何列向量，用来取三个法向应变之和，或将同一标量放入三个法向分量；孔压耦合使用它并不引入额外的材料参数。

所有光滑域使用同一组材料参数，但各虚拟点的应力、塑性变量和算法切线 $C_{kg}^{\mathrm{ep}}$ 可以不同。$C_{kg}^{\mathrm{ep}}$ 是本构返回的 $4\times4$ 算法切线矩阵，即有效应力对机械应变试算增量的导数；上标 $\mathrm{ep}$ 表示弹塑性，处于弹性时它就是弹性矩阵。因此这里保留状态索引 $kg$，不能因材料均匀就把所有积分点的切线当作同一常矩阵。局部矩阵仍按节点编号组装。

8.1 骨架内力、K11 与 K13

$\Delta\boldsymbol\varepsilon_{kg}^{\mathrm{mech}}$ 是扣除热应变后的机械应变增量，上标 $\mathrm{mech}$ 表示机械部分；$\boldsymbol f_{u,k}^{\mathrm{skel}}$ 是光滑域 $k$ 的骨架内力贡献，下标 $u$ 表示力学方程。$K_{11}^k$ 是该域的力学切线贡献，$K_{13}^k$ 是温度经热应变影响力学方程的耦合贡献。由于总应变要扣除热应变，温度项带负号，因此全局矩阵中的对应块是 $-K_{13}$。内力与切线必须使用同一组 $\widetilde B_{kg},\widetilde N_{kg},W_{kg}$。

J2、DP 等本构均使用这一组装形式。每轮 Newton 迭代都从上一个已收敛时间步的状态出发，由 UMAT 计算本轮试算的应力和算法切线。只有当前时间步收敛后，才保存四个虚拟点的本构状态，供下一步使用。MEX 使用预先算好的几何矩阵、组装所需的材料参数和最新算法切线，不负责本构更新。

8.2 K12 与 K21：使用平均矩阵

定义型函数的体积积分行向量 $I_{N,k}^V$；这里的 $I$ 只是积分结果的名称，不是单位矩阵，也不是节点编号：

则

此处使用第 4.3 节的平均矩阵 $\overline B_k^V$；它在光滑域内为常数，可以移到积分外面，剩下只对原始 T3 的 N 积分。这里不需要计算四个虚拟点处的原始 T3 型函数。按第 6 节的加减梯度写法，$\overline B_k^V$ 也正好等于四个 $\widetilde B_{kg}$ 的平均。组装后仍有 $K_{21}=K_{12}^T$。

由第 7.2 节的 $\overline N_k^V=I_{N,k}^{V}/V_k$，$K_{12}^k$ 等价于 $V_k(\overline B_k^V)^T\boldsymbol m\overline N_k^V$。

注意： 将面积平均 $\overline N_k^A$ 乘以 $V_k$，通常不等于所需的体积积分。

机械方程中的孔压项保持 $-K_{12}\boldsymbol P$，流体方程中的变形项保持 $K_{21}\Delta\boldsymbol U$。即使孔压处处相同，总应力中的 $-p\boldsymbol m$ 及边界上施加的力仍需共同满足平衡。规定孔压仍不会自动施加机械壁面压力。

8.3 储水、热容量、温度引起的孔压变化与扩散

令

这里 $M_k^V$ 就是把原始 T3 型函数的乘积按体积积分后得到的矩阵（但注意积分区域不是一整个三角形）。以下材料系数在整个域内相同，均不带 $k$；下标 $s,f$ 分别标记固体颗粒、孔隙流体：

则

$K_{22S}$、$K_{22P}$ 分别为储水、Darcy 流动矩阵，$K_{23}$ 为温度引起流体含量变化的耦合矩阵，$K_{33X}$、$K_{33G}$ 分别为热容量、热传导矩阵。下标字母 $S,P,X,G$ 沿用代码中的矩阵名称，不是额外的物理变量。轴对称模型仍有三个法向应变和应力分量。这里 $\boldsymbol m^T\boldsymbol m=3$，所以体应变及热膨胀中的三个法向分量包含环向分量，不能因计算域是 r–z 平面就把 3 改为 2。