几何单元的矩阵列式

Notation

在确定单元的几何类型, 并已知单元的节点坐标之后, 理应可以得到关于该单元的任何几何信息, 比如单元内部积分点的全局坐标, 形函数的取值, 形函数偏导数的取值等等. 这些在转化为代码中实际的矩阵计算时, 需要约定矩阵的行列分别对应的指标类型. 接下来将建立公式与实际计算中矩阵列式之间一一对应的关系.

在进行公式推导和变量定义之前, 有必要首先声明文中用到的 Notation:

• 使用 \(i,j,k,\ldots\) 表示空间和参数指标, 一般放在下标位置
• 使用 \(a, b,c,\ldots\) 表示单元节点的指标, 一般放在上标位置
• 如无特别说明, 均默认指标使用 Einstein 求和约定
• 公式字体 \(varible\) 表示一般的公式推导; 字体 \(\mathrm{matrix}\) 表示对应变量的矩阵形式, 并一般使用方括号 \([\cdot]\) 包围. \([\cdot]_{ij}\) 表示矩阵分量
• 按照 MATLAB 惯例, \([1 ~ 2 ~ 3]\) 或 \([1,2,3]\) 表示行向量; \([1;2;3]\) 表示列向量

定义的标量值

变量名	变量说明
nedim	参数坐标维数
nsdim	空间维数
nndof	节点自由度数
nen	单元节点数
nedof	单元自由度数, nedof = nndof * nen

定义的向量与矩阵值

变量名	变量维度	变量说明
x	1 x nsdim	全局坐标
x'	1 x nsdim	单元坐标
N	1 x nen	形函数在给定某一参数坐标下的取值
encrd	nen x nsdim	单元节点在全局坐标系下的坐标
encrd_ecsys	nen x nsdim	单元节点在单元坐标系下的坐标
Q	nsdim x nsdim	坐标变换矩阵
J	nsdim x nedim	全局坐标系下的 Jacobian 矩阵
J'	nsdim x nedim	单元坐标系下的 Jacobian 矩阵
J''	nedim x nedim	去除零行后的 Jacobian 矩阵 (方阵)
dNdxi	nen x nedim	形函数对参数坐标的偏导数
dNdx	nen x nsdim	形函数对全局坐标的偏导数
dNdx'	nen x nedim	形函数对单元坐标的偏导数
epsilon	[3|4|6] x 1	积分点处的应变
d	nedof x 1	单元自由度变量

参数坐标得到全局坐标

给定单元的某一参数坐标 \(\xi \in \mathbb{R}^{\rm nedim}\) 后, 该点的全局坐标 \(x \in \mathbb{R}^{\rm nsdim}\) 为: \[ x(\xi) = x^a N^a(\xi) \tag{1} \] 其变量关系为 \[ \boxed{ \mathrm{ [x] = [N] [encrd]} } \]

Jacobian

Jacobian 可以看作是参数坐标系中的微元 \(\mathrm{d} \xi\) 到全局坐标系中的微元 \(\mathrm{d} x\) 的仿射变换 \[ \mathrm{d} x = J \cdot \mathrm{d} \xi \] 写成分量形式为 \[ \mathrm{d} x_i = \frac{\partial x_i}{\partial \xi_j} \mathrm{d} \xi_j \] 记 Jacobian 的矩阵形式为 \[ [\mathrm{J}]_{ij} = \frac{\partial x_i}{\partial \xi_j} \] 当使用式 (1) 表示的映射关系时, Jacobian 矩阵可以表示为形函数关于参数坐标导数的线性组合: \[ [\mathrm{J}]_{ij} = x_i^a \frac{\partial N^{a}}{\partial \xi_j} \] 将形函数关于参数坐标的导数组成的矩阵记作变量 dNdXi, Jacobian 矩阵就可以表示为 \[ \boxed{ [\mathrm{J}] = \rm [encrd]^{\top} [dNdxi] } \tag{2} \]

\(\partial N / \partial x\) 矩阵

对于某一点处形函数对 \(x\) 的偏导数, 可以通过链式法则获得: \[ \frac{\partial N^{a}}{\partial x_i} = \frac{\partial N^{a}}{\partial \xi_j} \frac{\partial \xi_j}{\partial x_i} \] 上式中, 当单元类型确定之后, 形函数关于参数坐标 \(\xi\) 的导数就是已知并有解析表达式的, 但是参数坐标关于全局坐标的函数表示, 并不总是那么方便地得到. 但庆幸地是, 我们只需要知道某一点上的取值, 而这可以通过反函数定理很快地得到: \[ \left.\frac{\partial x_i}{\partial \xi_k}\right|_{\xi} \left.\frac{\partial \xi_k}{\partial x_j}\right|_{x} = \delta_{ij} \] 当 nsdim = nedim 时, 可以通过 Jacobian 逆矩阵得到 \[ \boxed{ \mathrm{[dNdx] = [dNdxi]} [J]^{-1} } \] 当 nsdim > nedim 时 (空间维度总是大于或等于单元维度), 建立单元坐标系 \(\{e_1, e_2, e_3\}\) 为

1. 二维单元: \(e_3\) 垂直于单元所在平面, \(e_1\) 平行于边 \(1-2\) 并指向节点 \(2\), \(e_2 = e_3 \times e_1\) ;
2. 一维单元: \(e_1\) 从节点 \(1\) 指向节点 \(2\), 其它基矢量与 \(e_1\) 满足右手坐标系关系.

如果记 \(\{e_1, e_2, e_3\}\) 在旧坐标系下的坐标矩阵为 \[ [Q] \triangleq [e_1 ~ e_2 ~ e_3] \] 那么, 新旧坐标之间的变换关系为 \[ [x] = [x'] [Q]^{\top} , \quad [x'] = [x] [Q] \] 节点在新坐标系下坐标为 \(x'^{a}\), 并且单元内所有点的坐标分量 \(x'_i\) ( 二维单元 \(i=3\), 一维单元 \(i=2,3\) ) 为常数. 如果记变换后的单元的节点坐标变量为 encrd_ecsys, 那么与变量 encrd 之间的关系为 \[ \boxed{ \rm [encrd\_ecsys] = [encrd] [Q] } \tag{3} \] 在新坐标系下, 用形函数表示为参数坐标的映射关系为 \[ x' = x'^a N^a(\xi) \] 式中的形函数形式上与未变换之前的形函数是一致的, 因此微元之间的关系同样为 \[ \mathrm{d} x' = J' \cdot \mathrm{d} \xi \] 式中, \(J'\) 为新坐标系下的 Jacobian, 类似于式 (2) \[ [\mathrm{J}'] = \rm [encrd\_ecsys]^{\top} [dNdxi] \] 再根据坐标变换关系 (3), 就得到不同坐标系下 Jacobian 之间的关系 \[ \boxed{ [\mathrm{J}'] = \rm [Q]^{\top} [encrd]^{\top} [dNdxi] = [Q]^{\top} [J] } \]

矩阵 \([J']\) 大于 nedim 的行元素全部等于 0. 去除零行之后得到 \([J'']\): nedim x nedim, 得到在新坐标系下的形函数偏导数为: \[ \mathrm{[dNdx'] = [dNdxi]} [J'']^{-1} \] 如果希望得到关于旧坐标系的偏导数值, 那么根据坐标变换关系: \[ \frac{\partial ~~}{\partial x_{\beta}} = \frac{\partial ~~}{\partial x'_{i}} \frac{\partial x'_{i}}{\partial x_{\beta}} = Q_{\beta i} \frac{\partial ~~}{\partial x'_{i}} \] 需要强调的是上式中指标的范围. \(\beta\) 是旧坐标系下的指标, 因此范围为 \(1,2,3\). \(i\) 是新坐标系下的指标, 在新坐标系下, \(i\) 的范围是 \(1,2\) (二维单元) 或 \(1\) (一维单元). 写成矩阵的形式为 \[ \boxed{ \rm [dNdx] = [dNdx'] [Q(0,0,nsdim,nedim)]^{\top} } \]

六节点单元, 形函数没有问题, 积分点上每次计算得到的应变增量没有问题, 应力计算应该也是没有问题的, 但是在线性第二次迭代时, 得不到正确结果? 这是为什么?

问题应该出现在 src\PField\FElement\FSolid2D\FElemSld2D.cpp 中的迭代过程. 需要再想想

应该不是迭代方法的问题, 还是第一次迭代步中最后的位移计算出错了

B 矩阵

\(B\) 矩阵的形式与应变向量的排布相关. 在小位移假设下, 应变总是关于位移对 \(x\) 偏导数的线性组合, 也即关于 \(\partial N / \partial x\) 的线性组合: \[ \varepsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) \] 对位移使用形函数进行插值, 得到 \[ u_{i}(x) = u_{i}^{a} N^{a}(x) \tag{4} \] 式中, \(u_{i}^{a}\) 为节点上的位移分量, 为待求解的未知量. 将上式代入到小应变的表达式当中得到 \[ \varepsilon_{ij} = \frac{1}{2}(u_{i}^{a} \frac{\partial N^{a}}{\partial x_j} + u_{j}^{a} \frac{\partial N^{a}}{\partial x_i}) \tag{5} \] 将上式写成矩阵形式, 得到 \[ \boxed{ \rm [epsilon] = [B] [d] } \] 接下来我们将阐述上式中变量具体的排布方式. 对于变量 \(\rm [epsilon]\), 我们使用 Voigt 记法, 将应变存储为向量形式: \[ [\mathrm{epsilon}] = [\varepsilon_{11}; ~ \varepsilon_{22}; ~ 2\varepsilon_{12}] \quad \text{二维情景} \\ [\mathrm{epsilon}] = [\varepsilon_{11}; ~ \varepsilon_{22}; ~ \varepsilon_{33}; ~ 2\varepsilon_{23}; ~ 2\varepsilon_{13}; ~ 2\varepsilon_{12}] \quad \text{三维情景} \] 对于变量 \([d]\), 设节点 \(a\) 上自由度变量排布为 \[ [\mathrm{u}^a] = [u_1^a ; ~ u_2^a ; ~ \cdots u_{\mathrm{nndof}}^a] \] 单元自由度向量通过 \([u^{a}]\) 进行组装 \[ \rm [d] = [u^1 ; u^2 ; \cdots ; u^{nen}] \] 在定义单元自由度向量和应变向量之后, 式 (5) 中将之联系起来的矩阵 \(B\) 就唯一确定了. 我们首先使用定义的单元自由度向量表示式 (4) 中的 \(u(x)\) \[ [\mathrm{u}](x) = N^{a}(x) [\mathrm{u}^a] = [N^{1}(x)[\mathrm{I}] , N^2(x)[\mathrm{I}], \ldots , N^{\mathrm{nen}}(x)[\mathrm{I}]] [\mathrm{d}] \] 式中, \([\mathrm{I}]\) 为维度 nndof 的单位阵. 定义算子矩阵 \[ [\mathrm{D}] = \begin{bmatrix} \frac{\partial }{\partial x_1} & \\ & \frac{\partial }{\partial x_2} \\ \frac{\partial }{\partial x_2} & \frac{\partial }{\partial x_1} \end{bmatrix} (\text{二维情景}), \quad \begin{bmatrix} \frac{\partial }{\partial x_1} & & \\ & \frac{\partial }{\partial x_2} & \\ & & \frac{\partial }{\partial x_3} \\ & \frac{\partial }{\partial x_3} & \frac{\partial }{\partial x_2} \\ \frac{\partial }{\partial x_3} & & \frac{\partial }{\partial x_1} \\ \frac{\partial }{\partial x_2} & \frac{\partial }{\partial x_1} & \end{bmatrix} (\text{三维情景}) \] 由此应变可以表示为 \[ [\mathrm{epsilon}] = [\mathrm{D}] [\mathrm{u}](x) = [\mathrm{D}] [N^{1}(x)[\mathrm{I}] , N^2(x)[\mathrm{I}], \ldots , N^{\mathrm{nen}}(x)[\mathrm{I}]] [\mathrm{d}] \] 如果定义 \[ [\mathrm{D}] N^{a}(x) = [\mathrm{B}^{a}] \] 那么 \[ [\mathrm{B}] = [\mathrm{D}] [N^{1}(x)[I] , N^2(x)[I], \ldots , N^{\mathrm{nen}}(x)[I]] = [\mathrm{B}^1, \mathrm{B}^2, \ldots, \mathrm{B}^{\mathrm{nen}}] \]