基于真实转角与“试探-修正”插值的四边形单元

Quadrilateral Element Based on True Rotation Angle and “Trial-correction” Interpolation

摘要

为解决不同单元间的兼容性问题，提出一种“试探-修正”位移插值方法，用于构造具有旋转自由度的四节点四边形平面单元.基于“试探-修正”插值的四节点四边形平面单元以平动位移和转动角度作为节点参数，采用双线性插值对位移场进行试探，根据转角的偏差，采用双三次插值对位移场进行校正.通过分片测试证明了该方法构造的四节点四边形平面单元的收敛性，并用三个算例进一步验证了其性能.数值结果表明，基于“试探-修正”插值的四节点四边形平面单元不仅收敛速度快，数值精度高，还能与梁单元兼容，方便混合单元建模.此外，“试探-修正”插值方法可参数化，具有良好的可扩展性，为今后研究其他具有真实转角的单元奠定了基础.

Abstract

To address the problem of incompatibility between different elements， the paper introduces a “trial-correction” displacement interpolation scheme. This scheme is utilized to construct a four-node quadrilateral plane element that incorporates a drilling degree of freedom. The “trial-correction” interpolation-based four-node quadrilateral plane element takes translational displacements and the drilling degree of freedom as nodal parameters， and higher-order interpolation functions are employed to approximate the displacement field. Firstly， bi-linear interpolation is used to trail the displacement fields. And then， according to the deviation of the drilling degree of freedom， the displacement fields are corrected with bi-cubic interpolation. The convergence of the “trial-correction” interpolation-based four-node quadrilateral plane element is proved by the patch test， and its performance is further verified by three examples. The numerical results show that the “trial-correction” interpolation-based four-node quadrilateral plane element has a high convergence rate， a high numerical accuracy， and can be compatible with beam elements for convenient mix-element modeling. Moreover， the “trial-correction” interpolation method is parameterized and has good extendibility， which lays a foundation for future research about other elements with true rotation angle.

Graphical abstract

关键词

有限元法 / 插值 / 数值算法收敛性 / 四边形平面单元 / 真实转角

随着计算机技术的发展，有限元法已被广泛应用于解决各种科学和工程问题^［1］.对于同时具有不同结构的复杂机构，在运用有限元法进行分析的过程中，可能会遇到奇异刚度矩阵的问题，这会导致分析无法成功^［2-3］.而带有旋转自由度的平面单元可以通过与板单元结合来构造壳体单元，从而使壳体单元避免奇异刚度矩阵问题.不过必须指出的是， 开发高性能的低阶带有旋转自由度的平面单元仍然是研究人员面临的挑战.

Allman^［4］提出对节点使用二次插值构造出带真实转角的三角形平面单元.在此基础之上， Cook^［5］证明了通过线性应变三角形单元和角节点旋转的三角形单元之间的转换技术可以得到Allman等效元，并在此基础上提出了一种带转动自由度的四边形单元.近十年间，带转角的平面单元有了多方面的完善与扩充，其中Shang等^［6］提出了一种带有真实旋转自由度的非对称膜单元， 包括开发了具有真实转动自由度的非对称三角形单元^［7］和非对称四边形单元^［8］. Rezaiee-Pajand等^［9］使用泰勒级数的方法对应变场进行插值， 提出了基于应变的三角形膜单元. Boutagouga^［10］提出了一种具有真钻井旋转自由度的基于位移的膜单元公式. 面内转角在杂交有限元、混合有限元中也扮演重要角色， 例如Daszkiewicz等^［11］基于Hu-Washizu原理， 提出了含钻井旋转的混合四节点壳单元， 并将该单元用于非线性分析.Wang^［12］提出了一种具有真实转动自由度的平面单元.文颖等^［13］构建了一种含有面内自由度的四节点板壳单元，用于非线性分析中. Sangtarash等^［14］利用Airy应力函数的解析解以及四节点单元的等参位移构造了一个高性能且带有旋转自由度的非对称四边形单元.基于修正偶应力理论，Wu等^［15］采用罚函数法约束引入的旋转量并构造了一个八节点四边形单元，而Long等^［16］通过引入薄板弯曲的假设，修改偶应力弹性的三维控制方程构造了类似于非线性有限元模型性能的Trefftz板单元.文献［17］对平面单元的发展做了比较详细的归纳.

与传统等参单元相比，这些包含真实转角的单元表现出更好的性能.然而，这些单元在应力/应变贴片测试中未能严格满足要求，其收敛性和数值稳定性有待进一步探讨.此外，带转动自由度的单元在公式推导方面过于复杂，不利于在工程实际中应用.

本文将平动位移和转动位移作为节点参数，提出了一种“试探-校正”插值方法，用于构建参数化位移插值. 基于“试探-校正”插值方法，推导了一个具有真实旋转自由度的四节点四边形平面单元（four-node-quadrilateral plane element with trial-correction interpolation，plane-Q4-TC单元）. plane-Q4-TC单元具有三次位移场， 展现出较高的数值精度和较快的收敛速度， 可以与平面梁单元兼容， 便于进行混合单元建模. 同时， 参数化的“试探-校正”插值方案具有良好的扩展性， 可用于构建梁、壳和实体等多种单元， 从而形成一个多维度兼容的有限元系统.

本文内容安排如下：第1章为plane-Q4-TC单元的“试探-校正”插值方案的理论公式推导，并通过分片测试证明所构建单元的收敛性.第2章用数值算例来验证plane-Q4-TC单元的性能.第3章为结论.

1 理论公式

1.1 节点参数

假设一个无限小平面单元的真实转动如图1所示， 图1（a）和图1（b）分别指刚性单元的真实转角以及可变形单元的真实转角， 刚性单元的真实转角{\theta }_u和{\theta }_v相同， 因此可用转角\overline{\theta }表示：

\overline{\theta }={\theta }_u=-\frac{\partial u}{\partial y}={\theta }_v=\frac{\partial v}{\partial x}（1）

可变形单元的真实转角{\theta }_u和{\theta }_v不相同， 因此真实转角{\theta }_z可定义为：

{\theta }_z=\frac{\mathrm{1}}{\mathrm{2}}\left({\theta }_u+{\theta }_v\right)=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)（2）

包含真实转角的四节点四边形平面单元模型如图2所示，图2（a）和图2（b）分别指物理坐标系和自然坐标系下的平面单元.每个单元节点p_i(i=\mathrm{1}~\mathrm{4})都包含3个参数，分别是平动位移u、v以及绕z轴旋转的真实转角{\theta }_z，位移场u(\xi ,\eta )和v(\xi ,\eta )均由以上3个节点参数插值得到.其中， (\xi ,\eta )为自然坐标， 与物理坐标之间的变换关系为：

\begin{array}{l}x(\xi ,\eta )=\mathit{N}(\xi ,\eta ){\mathit{x}}_{\mathrm{e}}\\ y(\xi ,\eta )=\mathit{N}(\xi ,\eta ){\mathit{y}}_{\mathrm{e}}\end{array}（3）

式中： {\mathit{x}}_{\mathrm{e}}和{\mathit{y}}_{\mathrm{e}}分别是x和y方向节点坐标矢量； \mathit{N}(\xi ,\eta )是形函数矩阵， 可表达为：

\begin{array}{c}\mathit{N}(\xi ,\eta )=\left[\begin{array}{cccc}{n}_{\mathrm{1}}& n_{\mathrm{2}}& n_{\mathrm{3}}& n_{\mathrm{4}}\end{array}\right]\\ n_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{4}}\left(\mathrm{1}-\xi \right)(\mathrm{1}-\eta )\\ n_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{4}}\left(\mathrm{1}+\xi \right)(\mathrm{1}-\eta )\\ n_{\mathrm{3}}=\frac{\mathrm{1}}{\mathrm{4}}\left(\mathrm{1}+\xi \right)(\mathrm{1}+\eta )\\ n_{\mathrm{4}}=\frac{\mathrm{1}}{\mathrm{4}}\left(\mathrm{1}-\xi \right)(\mathrm{1}+\eta )\end{array}（4）

1.2 位移插值

本节采用“试探-修正”法构造位移场， 假设位移场u(\xi ,\eta )和v(\xi ,\eta )均由试探项位移u^{\mathrm{t}}(\xi ,\eta )、v^{\mathrm{t}}(\xi ,\eta )和修正项位移u^{\mathrm{c}}(\xi ,\eta )、v^{\mathrm{c}}(\xi ,\eta )组成， 表达式为：

\begin{array}{l}u(\xi ,\eta )=u^{\mathrm{t}}(\xi ,\eta )+u^{\mathrm{c}}(\xi ,\eta )\\ v(\xi ,\eta )=v^{\mathrm{t}}(\xi ,\eta )+v^{\mathrm{c}}(\xi ,\eta )\end{array}（5）

式中：试探项位移场u^{\mathrm{t}}(\xi ,\eta )、v^{\mathrm{t}}(\xi ,\eta )由双线性插值得到， 表达式为：

\begin{array}{l}{u}^{\mathrm{t}}(\xi ,\eta )=\mathit{N}(\xi ,\eta ){\mathit{u}}_{\mathrm{e}}={\mathit{N}}_u^{\mathrm{t}}(\xi ,\eta ){\mathit{d}}_{\mathrm{e}}\\ v^{\mathrm{t}}(\xi ,\eta )=\mathit{N}(\xi ,\eta ){\mathit{v}}_{\mathrm{e}}={\mathit{N}}_v^{\mathrm{t}}(\xi ,\eta ){\mathit{d}}_{\mathrm{e}}\end{array}（6）

式中： {\mathit{u}}_{\mathrm{e}}和{\mathit{v}}_{\mathrm{e}}分别是沿着x轴和y轴的节点位移；{\mathit{d}}_{\mathrm{e}}指节点参数的位移向量， 表达式为：

\begin{array}{l}{\mathit{d}}_{\mathrm{e}}={\left[\begin{array}{cccc}{\mathit{d}}_{\mathrm{1}}^{\mathrm{T}}& {\mathit{d}}_{\mathrm{2}}^{\mathrm{T}}& {\mathit{d}}_{\mathrm{3}}^{\mathrm{T}}& {\mathit{d}}_{\mathrm{4}}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}\\ {\mathit{d}}_i=\left[\begin{array}{ccc}{u}_i& v_i& {\theta }_{z,i}\end{array}\right]\end{array}（7）

形函数{\mathit{N}}_u^{\mathrm{t}}(\xi ,\eta )和{\mathit{N}}_v^{\mathrm{t}}(\xi ,\eta )的表达式为：

\begin{array}{l}{\mathit{N}}_u^{\mathrm{t}}(\xi ,\eta )=\\ \left[\begin{array}{cccc}{\mathit{N}}_{u,\mathrm{1}}^{\mathrm{t}}(\xi ,\eta )& {\mathit{N}}_{u,\mathrm{2}}^{\mathrm{t}}(\xi ,\eta )& {\mathit{N}}_{v,\mathrm{3}}^{\mathrm{t}}(\xi ,\eta )& {\mathit{N}}_{u,\mathrm{4}}^{\mathrm{t}}(\xi ,\eta )\end{array}\right]\\ {\mathit{N}}_{u,i}^{\mathrm{t}}=\left[\begin{array}{ccc}{n}_i& \mathrm{0}& \mathrm{0}\end{array}\right]\end{array}（8）

\begin{array}{l}{\mathit{N}}_v^{\mathrm{t}}(\xi ,\eta )=\\ \left[\begin{array}{cccc}{\mathit{N}}_{v,\mathrm{1}}^{\mathrm{t}}(\xi ,\eta )& {\mathit{N}}_{v,\mathrm{2}}^{\mathrm{t}}(\xi ,\eta )& {\mathit{N}}_{v,\mathrm{3}}^{\mathrm{t}}(\xi ,\eta )& {\mathit{N}}_{v,\mathrm{4}}^{\mathrm{t}}(\xi ,\eta )\end{array}\right]\\ {\mathit{N}}_{v,i}^{\mathrm{t}}=\left[\begin{array}{ccc}\mathrm{0}& n_i& \mathrm{0}\end{array}\right]\end{array}（9）

由式（2）得知， 试探项构造的真实旋转角可以计算为：

{\theta }_z^{\mathrm{t}}(\xi ,\eta )=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial v^{\mathrm{t}}}{\partial x}-\frac{\partial u^{\mathrm{t}}}{\partial y}\right)={\mathit{B}}_z^{\mathrm{t}}(\xi ,\eta ){\mathit{d}}_{\mathrm{e}}（10）

式中： {\mathit{B}}_z^{\mathrm{t}}(\xi ,\eta )的表达式可以写成：

\begin{array}{c}{\mathit{B}}_z^{\mathrm{t}}(\xi ,\eta )=\left[\begin{array}{cccc}{\mathit{B}}_{z,\mathrm{1}}^{\mathrm{t}}& {\mathit{B}}_{z,\mathrm{2}}^{\mathrm{t}}& {\mathit{B}}_{z,\mathrm{3}}^{\mathrm{t}}& {\mathit{B}}_{z,\mathrm{4}}^{\mathrm{t}}\end{array}\right]\\ {\mathit{B}}_{z,i}^{\mathrm{t}}=\left[\begin{array}{ccc}-\frac{\mathrm{1}}{\mathrm{2}}\frac{\partial n_i}{\partial y}& \frac{\mathrm{1}}{\mathrm{2}}\frac{\partial n_i}{\partial x}& \mathrm{0}\end{array}\right]\end{array}（11）

结合式（3）， 将式（11）的偏导数写成：

\left[\begin{array}{c}\frac{\partial n_i}{\partial x}\\ \frac{\partial n_i}{\partial y}\end{array}\right]={\mathit{J}}^{-\mathrm{1}}(\xi ,\eta )\left[\begin{array}{c}\frac{\partial n_i}{\partial \xi }\\ \frac{\partial n_i}{\partial \eta }\end{array}\right]（12）

式中：\mathit{J}(\xi ,\eta )是雅可比矩阵， 表达式为：

\mathit{J}(\xi ,\eta )=\left[\begin{array}{c}\frac{\partial \mathit{N}(\xi ,\eta )}{\partial \xi }\\ \frac{\partial \mathit{N}(\xi ,\eta )}{\partial \eta }\end{array}\right]\left[\begin{array}{cc}{\mathit{x}}_{\mathrm{e}}& {\mathit{y}}_{\mathrm{e}}\end{array}\right]（13）

对于平面单元， 节点p_i(i=\mathrm{1}~\mathrm{4})可以满足以下表达式：

\left\{\begin{array}{l}{u}_i^{\mathrm{t}}=u^{\mathrm{t}}(\xi ,\eta )=u_i\\ v_i^{\mathrm{t}}=v^{\mathrm{t}}(\xi ,\eta )=v_i\\ {\theta }_i^{\mathrm{t}}={\theta }_i^{\mathrm{t}}(\xi ,\eta )={\mathit{B}}_z^{\mathrm{t}}({\xi }_i,{\eta }_i){\mathit{d}}_{\mathrm{e}}\end{array}\right.（14）

式中：{\xi }_i和{\eta }_i是节点p_i的自然坐标.修正项u^{\mathrm{c}}(\xi ,\eta )、v^{\mathrm{c}}(\xi ,\eta )和.{\theta }_z^{\mathrm{c}}(\xi ,\eta )需要满足以下情况：

\left\{\begin{array}{l}{u}_i^{\mathrm{c}}=u_i-u_i^{\mathrm{t}}=\mathrm{0}\\ v_i^{\mathrm{c}}=v_i-v_i^{\mathrm{t}}=\mathrm{0}\\ {\theta }_{z,i}^{\mathrm{c}}={\theta }_{z,i}-{\theta }_{z,i}^{\mathrm{t}}={\mathit{B}}_{z,i}^{\mathrm{c}}{\mathit{d}}_{\mathrm{e}}\end{array}\right.（15）

式中：

\begin{array}{c}{\mathit{B}}_{z,i}^{\mathrm{c}}=\left[\begin{array}{cccc}{\mathit{B}}_{z,\mathrm{1}}^{\mathrm{c}}& {\mathit{B}}_{z,\mathrm{2}}^{\mathrm{c}}& {\mathit{B}}_{z,\mathrm{3}}^{\mathrm{c}}& {\mathit{B}}_{z,\mathrm{4}}^{\mathrm{c}}\end{array}\right]\\ {\mathit{B}}_{z,j}^{\mathrm{c}}({\xi }_i,{\eta }_i)=\left[\begin{array}{ccc}\frac{\mathrm{1}}{\mathrm{2}}\frac{\partial n_j}{\partial y}& -\frac{\mathrm{1}}{\mathrm{2}}\frac{\partial n_j}{\partial x}& n_j({\xi }_i,{\eta }_i)\end{array}\right]\end{array}（16）

通过三次插值的方法构造修正项位移， u^{\mathrm{c}}(\xi ,\eta )和v^{\mathrm{c}}(\xi ,\eta )可表达为：

\begin{array}{c}{u}^{\mathrm{c}}(\xi ,\eta )=\mathit{P}(\xi ,\eta )\mathit{\alpha }\\ v^{\mathrm{c}}(\xi ,\eta )=\mathit{P}(\xi ,\eta )\mathit{\beta }\end{array}（17）

式中：\mathit{P}(\xi ,\eta )是由12个插值基函数组成的行向量， 如表1所示； α 和 β 是待定系数， 可表达为：

\begin{array}{c}\mathit{\alpha }={\left[\begin{array}{cccc}{\alpha }_{\mathrm{1}}& {\alpha }_{\mathrm{2}}& \cdots & {\alpha }_{\mathrm{12}}\end{array}\right]}^{\mathrm{T}}\\ \mathit{\beta }={\left[\begin{array}{cccc}{\beta }_{\mathrm{1}}& {\beta }_{\mathrm{2}}& \cdots & {\beta }_{\mathrm{12}}\end{array}\right]}^{\mathrm{T}}\end{array}（18）

为了求解待定系数， 在节点p_i处建立以下方程：

\left[\begin{array}{c}\mathit{P}({\xi }_i,{\eta }_i)\\ {\mathit{P}}_{\xi }({\xi }_i,{\eta }_i)\\ {\mathit{P}}_{\eta }({\xi }_i,{\eta }_i)\end{array}\right]\mathit{\alpha }=\left[\begin{array}{c}{u}^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ u_{{}^{\xi }}^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ u_{{}^{\eta }}^{\mathrm{c}}({\xi }_i,{\eta }_i)\end{array}\right]=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathit{J}({\xi }_i,{\eta }_i)\end{array}\right]\left[\begin{array}{c}{u}_i^{\mathrm{c}}\\ u_x^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ u_y^{\mathrm{c}}({\xi }_i,{\eta }_i)\end{array}\right]（19）

\left[\begin{array}{c}\mathit{P}({\xi }_i,{\eta }_i)\\ {\mathit{P}}_{\xi }({\xi }_i,{\eta }_i)\\ {\mathit{P}}_{\eta }({\xi }_i,{\eta }_i)\end{array}\right]\mathit{\beta }=\left[\begin{array}{c}{v}^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ v_{{}^{\xi }}^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ v_{{}^{\eta }}^{\mathrm{c}}({\xi }_i,{\eta }_i)\end{array}\right]=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathit{J}({\xi }_i,{\eta }_i)\end{array}\right]\left[\begin{array}{c}{v}_i^{\mathrm{c}}\\ v_x^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ v_y^{\mathrm{c}}({\xi }_i,{\eta }_i)\end{array}\right]（20）

式中：下标ξ、η、x和y表示对应的偏导数. 假设修正项u^{\mathrm{c}}(\xi ,\eta )和v^{\mathrm{c}}(\xi ,\eta )呈现刚性节点特征， 即修正项在节点处不引起应变， 表达为：

\left\{\begin{array}{c}{u}_x^{\mathrm{c}}({\xi }_i,{\eta }_i)=v_y^{\mathrm{c}}({\xi }_i,{\eta }_i)=\mathrm{0}\\ -u_y^{\mathrm{c}}({\xi }_i,{\eta }_i)=v_x^{\mathrm{c}}({\xi }_i,{\eta }_i)={\theta }_{z,i}^{\mathrm{c}}\end{array}\right.（21）

将式（21）代入式（15）可得到：

\left[\begin{array}{c}{u}_i^{\mathrm{c}}\\ u_x^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ u_y^{\mathrm{c}}({\xi }_i,{\eta }_i)\end{array}\right]={\mathit{S}}_i^u{\mathit{d}}_{\mathrm{e}},\left[\begin{array}{c}{v}_i^{\mathrm{c}}\\ v_x^{\mathrm{c}}({\xi }_i,{\eta }_i)\\ v_y^{\mathrm{c}}({\xi }_i,{\eta }_i)\end{array}\right]={\mathit{S}}_i^v{\mathit{d}}_{\mathrm{e}}（22）

式中：{\mathit{S}}_i^u和{\mathit{S}}_i^v分别为：

{\mathit{S}}_i^u=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ -{\mathit{B}}_i^{\mathrm{c}}\end{array}\right],{\mathit{S}}_i^v=\left[\begin{array}{c}\mathrm{0}\\ {\mathit{B}}_i^{\mathrm{c}}\\ \mathrm{0}\end{array}\right]（23）

将式（22）代入式（19）和式（20）， 得到：

\left\{\begin{array}{l}{\mathit{A}}_i\mathit{\alpha }={\mathit{J}}_i^{\mathrm{e}}{\mathit{S}}_i^u{\mathit{d}}_{\mathrm{e}}\\ {\mathit{A}}_i\mathit{\beta }={\mathit{J}}_i^{\mathrm{e}}{\mathit{S}}_i^v{\mathit{d}}_{\mathrm{e}}\end{array}\right.（24）

其中，

{\mathit{A}}_i=\left[\begin{array}{c}\mathit{P}({\xi }_i,{\eta }_i)\\ {\mathit{P}}_{\xi }({\xi }_i,{\eta }_i)\\ {\mathit{P}}_{\eta }({\xi }_i,{\eta }_i)\end{array}\right],{\mathit{J}}_i^{\mathrm{e}}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathit{J}({\xi }_i,{\eta }_i)\end{array}\right]（25）

将式（24）写成联立方程的形式， 可表达为：

\left\{\begin{array}{c}\mathit{A}\mathit{\alpha }={\mathit{J}}_{\mathrm{e}}{\mathit{S}}_u{\mathit{d}}_{\mathrm{e}}\\ \mathit{A}\mathit{\beta }={\mathit{J}}_{\mathrm{e}}{\mathit{S}}_v{\mathit{d}}_{\mathrm{e}}\end{array}\right.（26）

其中，

\mathit{A}=\left[\begin{array}{c}{\mathit{A}}_{\mathrm{1}}\\ {\mathit{A}}_{\mathrm{2}}\\ {\mathit{A}}_{\mathrm{3}}\\ {\mathit{A}}_{\mathrm{4}}\end{array}\right],{\mathit{J}}_{\mathrm{e}}=\left[\begin{array}{cccc}{\mathit{J}}_{\mathrm{1}}^{\mathrm{e}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{J}}_{\mathrm{2}}^{\mathrm{e}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathit{J}}_{\mathrm{3}}^{\mathrm{e}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{J}}_{\mathrm{4}}^{\mathrm{e}}\end{array}\right]（27）

{\mathit{S}}_u=\left[\begin{array}{c}{\mathit{S}}_{\mathrm{1}}^u\\ {\mathit{S}}_{\mathrm{2}}^u\\ {\mathit{S}}_{\mathrm{3}}^u\\ {\mathit{S}}_{\mathrm{4}}^u\end{array}\right],{\mathit{S}}_v=\left[\begin{array}{c}{\mathit{S}}_{\mathrm{1}}^v\\ {\mathit{S}}_{\mathrm{2}}^v\\ {\mathit{S}}_{\mathrm{3}}^v\\ {\mathit{S}}_{\mathrm{4}}^v\end{array}\right]（28）

求解方程（26）， 求出待定系数：

\left\{\begin{array}{l}\mathit{\alpha }={\mathit{A}}^{-\mathrm{1}}{\mathit{J}}_{\mathrm{e}}{\mathit{S}}_u{\mathit{d}}_{\mathrm{e}}\\ \mathit{\beta }={\mathit{A}}^{-\mathrm{1}}{\mathit{J}}_{\mathrm{e}}{\mathit{S}}_v{\mathit{d}}_{\mathrm{e}}\end{array}\right.（29）

将式（29）代入方程（17），得到：

u^{\mathrm{c}}(\xi ,\eta )=\mathit{P}(\xi ,\eta ){\mathit{A}}^{-\mathrm{1}}{\mathit{J}}_{\mathrm{e}}{\mathit{S}}_u{\mathit{d}}_{\mathrm{e}}

v^{\mathrm{c}}(\xi ,\eta )=\mathit{P}(\xi ,\eta ){\mathit{A}}^{-\mathrm{1}}{\mathit{J}}_{\mathrm{e}}{\mathit{S}}_v{\mathit{d}}_{\mathrm{e}}（30）

因此面内位移的形函数为：

\begin{array}{l}{\mathit{N}}_u^c=\mathit{P}(\xi ,\eta ){\mathit{A}}^{-\mathrm{1}}{\mathit{J}}_{\mathrm{e}}{\mathit{S}}_u\\ {\mathit{N}}_v^c=\mathit{P}(\xi ,\eta ){\mathit{A}}^{-\mathrm{1}}{\mathit{J}}_{\mathrm{e}}{\mathit{S}}_v\end{array}（31）

将式（6）、式（30）、式（31）代入式（5）可得到：

u(\xi ,\eta )=\left({\mathit{N}}_u^{\mathrm{t}}(\xi ,\eta )+{\mathit{N}}_u^{\mathrm{c}}(\xi ,\eta )\right){\mathit{d}}_{\mathrm{e}}

v(\xi ,\eta )=\left({\mathit{N}}_v^{\mathrm{t}}(\xi ,\eta )+{\mathit{N}}_v^{\mathrm{c}}(\xi ,\eta )\right){\mathit{d}}_{\mathrm{e}}（32）

1.3 单元刚度矩阵

根据式（32）， 膜应变可计算为：

\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right]=\left[\begin{array}{c}\frac{\partial u}{\partial x}\\ \frac{\partial v}{\partial y}\\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\end{array}\right]=\left({\mathit{B}}_{\mathrm{m}}^{\mathrm{t}}+{\mathit{B}}_{\mathrm{m}}^{\mathrm{c}}\right){\mathit{d}}_{\mathrm{e}}（33）

式中：{\epsilon }_x和{\epsilon }_y代指膜单元的正应变；{\gamma }_{xy}代指剪应变；{\mathit{B}}_{\mathrm{m}}^{\mathrm{t}}和{\mathit{B}}_{\mathrm{m}}^{\mathrm{c}}分别对应试探项和修正项的应变矩阵， 可定义为：

{\mathit{B}}_{\mathrm{m}}^k=\left[\begin{array}{c}\frac{\partial {\mathit{N}}_u^k(\xi ,\eta )}{\partial x}\\ \frac{\partial {\mathit{N}}_v^k(\xi ,\eta )}{\partial y}\\ \frac{\partial {\mathit{N}}_u^k(\xi ,\eta )}{\partial y}+\frac{\partial {\mathit{N}}_v^k(\xi ,\eta )}{\partial x}\end{array}\right],(k=\mathrm{t},\mathrm{c})（34）

其中，关于x和y的偏导数可以计算为：

\left[\begin{array}{c}\frac{\partial {\mathit{N}}_l^k(\xi ,\eta )}{\partial x}\\ \frac{\partial {\mathit{N}}_l^k(\xi ,\eta )}{\partial y}\end{array}\right]={\mathit{J}}^{-\mathrm{1}}(\xi ,\eta )\left[\begin{array}{c}\frac{\partial {\mathit{N}}_l^k(\xi ,\eta )}{\partial \xi }\\ \frac{\partial {\mathit{N}}_l^k(\xi ,\eta )}{\partial \eta }\end{array}\right]=

{\mathit{J}}^{-\mathrm{1}}(\xi ,\eta )\left[\begin{array}{c}\frac{\partial {\mathit{P}}_l^k(\xi ,\eta )}{\partial \xi }\\ \frac{\partial {\mathit{P}}_l^k(\xi ,\eta )}{\partial \eta }\end{array}\right]{\mathit{A}}^{-\mathrm{1}}{\mathit{J}}_{\mathrm{e}}{\mathit{S}}_l,(l=u,v)（35）

令：

{\mathit{B}}_{\mathrm{m}}^{\mathrm{s}}={\mathit{B}}_{\mathrm{m}}^{\mathrm{t}}+{\mathit{B}}_{\mathrm{m}}^{\mathrm{c}}（36）

式中：{\mathit{B}}_{\mathrm{m}}^{\mathrm{s}}指膜单元的应变矩阵，膜刚度矩阵可表达为：

\begin{array}{l}{\mathit{k}}_{\mathrm{m}}^{\mathrm{e}}=h{\int }_{\Omega }{\mathit{B}}_{\mathrm{m}}^{\mathrm{T}}\mathit{C}{\mathit{B}}_{\mathrm{m}}\mathrm{d}x\mathrm{d}y=\\ h{\int }_{-\mathrm{1}}^{\mathrm{1}}{\int }_{-\mathrm{1}}^{\mathrm{1}}{\mathit{B}}_{\mathrm{m}}^{\mathrm{T}}\mathit{C}{\mathit{B}}_{\mathrm{m}}\mathit{J}(\xi ,\eta )\mathrm{d}\xi \mathrm{d}\eta \end{array}（37）

式中：h和Ω分别为单元的厚度和积分域； C 为本构矩阵. 对于平面应力问题和平面应变问题， C 可分别表示为：

\mathit{C}=\frac{E}{\mathrm{1}-{\nu }^{\mathrm{2}}}\left[\begin{array}{ccc}\mathrm{1}& \nu & \mathrm{0}\\ \nu & \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}-\nu }{\mathrm{2}}\end{array}\right]（38）

或者

\mathit{C}=\frac{E(\mathrm{1}-\nu )}{(\mathrm{1}+\nu )(\mathrm{1}-\mathrm{2}\nu )}\left[\begin{array}{ccc}\mathrm{1}& \frac{\nu }{\mathrm{1}-\nu }& \mathrm{0}\\ \frac{\nu }{\mathrm{1}-\nu }& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}-\mathrm{2}\nu }{\mathrm{2}(\mathrm{1}-\nu )}\end{array}\right]（39）

式中：E为弹性模量；\nu为泊松比.

1.4 分片试验

位移场的设置如表2所示， 分别产生恒定应变{\epsilon }_x,{\epsilon }_y和γ. 以{\epsilon }_x为例， 将节点坐标p_i代入位移函数， 可得到节点位移：

d_i=\left[\begin{array}{c}{u}_i\\ v_i\\ {\theta }_i\end{array}\right]=\left[\begin{array}{c}{a}_{\mathrm{0}}+a_{\mathrm{1}}{x}_i\\ \mathrm{0}\\ \mathrm{0}\end{array}\right]（40）

由式（40）可得：

\begin{array}{l}{\mathit{S}}_u{\mathit{d}}_{\mathrm{e}}=\mathrm{0}\\ {\mathit{S}}_v{\mathit{d}}_{\mathrm{e}}=\mathrm{0}\end{array}（41）

因此可推导出位移场：

\begin{array}{l}u(\xi ,\eta )=u^{\mathrm{t}}(\xi ,\eta )\\ v(\xi ,\eta )=v^{\mathrm{t}}(\xi ,\eta )\end{array}（42）

即在应变{\epsilon }_x为常数的情况下， plane-Q4-TC单元等价于经典等参元， 在{\epsilon }_x或{\epsilon }_y为常数的情况下也可得到类似结论. 因此， plane-Q4-TC单元能够通过分片试验.

2 数值算例

本章主要通过4个不同的算例：悬臂梁、MacNeal梁、Cook问题以及plane-Q4-TC单元在工程应用中的可行性分析， 来验证“试探-修正”法构造plane-Q4-TC单元的精度以及收敛性.

2.1 悬臂梁

悬臂梁的长宽比为10， 梁的左端受到固定约束， 在右端A点处施加垂直向下的拉力F_y， 具体的受力情况以及结构参数如图3所示. 该悬臂梁使用不同的有限单元网格（网格尺寸l_{\mathrm{e}}=1， 2.5， 5， 10， 20， 25， 50， 100 mm）进行离散化， 在A点处沿y方向的位移结果如表3所示.

作为比较， 本节使用等参元和ABAQUS S4R单元对悬臂梁案例进行静力学分析， 分析结果如表3所示， 以此来验证本研究构建的单元收敛性. 3种单元的位移结果如图4所示. 此外， 将ABAQUS的高阶四边形单元作为参考解绘制在图4中来直观体现各种单元的收敛趋势.

由图4可知，在相同的单元网格划分条件下， 使用“试探-修正”法构造的plane-Q4-TC单元的数值结果与参考值更加接近，比等参元和ABAQUS S4R的收敛速度更快.此外，从表3的数值结果可以得知，plane-Q4-TC单元使用10 mm网格离散时的计算结果优于其他两种单元使用2.5 mm网格划分的结果.plane-Q4-TC单元使用2.5 mm网格离散时求解的位移值与另外两种单元使用1 mm网格划分得到的位移值基本相等.

综合上述分析，使用“试探-修正”位移插值法构造的plane-Q4-TC单元具有很高的数值精度，在相同的单元尺寸下，该单元更加接近参考值，而且比等参元和ABAQUS S4R单元具有更快的收敛速度.

2.2 MacNeal梁

图5所示为MacNeal提出的细长悬臂梁模型， 其通常用来检验单元剪切自锁现象以及非常规网格对计算精度的影响，本次算例应用了两种梯形网格［图5（a）和图5（b）］来做本次测试.悬臂梁长度L= 6 mm，高度h=0.2 mm，厚度t=0.1 mm，杨氏模量E=107 MPa， 泊松比\nu=0.3.在载荷M作用下右端挠度的参考值为-0.005 4 mm；在载荷P作用下右端挠度的参考值为0.108 1 mm^［6］.将计算结果归一化并记录在表4中.

根据表4可以看出等参元的精度在处理MacNeal梁时出现剪切自锁现象，主要表现为精度欠佳， 在不同测试环境下最优精度仅为3.1%.为解决该问题，众多学者在理论上进行了不断完善， Shang等^［6］提出的US-Q4θ单元在该问题上甚至无限趋于解析解.

尽管本文所提出的plane-Q4-TC单元在精度上不及US-Q4θ单元、HS-A7单元等，但其精度依然超过了86%.

2.3 Cook问题

为了验证plane-Q4-TC单元的性能， 本节对Cook问题进行求解， 即求解非规则几何形状的悬臂梁在端部受剪问题. Cook悬臂梁示意图如图6所示， 梁的左端被固定， 而右端则受到面内剪力的作用. 通过3项指标对plane-Q4-TC单元的性能进行了评估：C点挠度的绝对值（v_C ）、A点的最大主应力（{\sigma }_{A\mathrm{m}\mathrm{a}\mathrm{x}}）以及B点的最小主应力（{\sigma }_{B\mathrm{m}\mathrm{i}\mathrm{n}}）.随后，将所得结果与其他学者的研究成果^［19-20］进行了对比分析，相关数据如表5所示. 因为plane-Q4-TC单元的应力和应变是不连续的， 所以表5中平面单元的每种情况均有3个应变和应力结果：S_{\mathrm{1}}指左侧单元的结果， S_{\mathrm{2}}指右侧单元的结果，S_{\mathrm{m}}指平均值.为了直观地比较不同方法的收敛性， 分别绘制了相应指标的收敛曲线，分别在图7~图9中展示.

图7为C点处的位移收敛曲线，可以观察到，plane-Q4-TC单元具有良好的精度，比Allman^［4］、Q4S^［21］和Pimpinelli^［22］提出的单元具有更好的性能，但略逊于GQ12M8^［16］、Zhang等^［19］和Choi等^［23］提出的单元.对于16×16的网格，只有3个对比数据， plane-Q4-TC平面单元的性能比Q4S^［21］好，但不及Choi等^［23］提出的单元.

图8为A点处的最大主应力的收敛曲线， 可以观察到S_{\mathrm{1}}的收敛曲线呈自上而下的趋势， 而S_{\mathrm{2}}的收敛曲线则与其他学者提出的单元趋势一致， 均为自下而上收敛. 此外，S_{\mathrm{2}}的精度高于S_{\mathrm{1}}. 对于较粗的网格， plane-Q4-TC单元的性能大致与Allman单元和Q4S单元相同， 但不如Zhang等提出的单元. 对于16×16的网格， plane-Q4-TC单元的性能与参考值一致， 优于Q4S单元和Choi等提出的单元.

图9为B点处的最小主应力的收敛曲线，就最小主应力而言，总体上S_{\mathrm{1}}展现出比S_{\mathrm{2}}更优的性能.对于包含2×2和4×4的网格规模，除了GQ12M8^［16］单元外， S_{\mathrm{1}}比其他单元更接近参考值.对于细网格， S_{\mathrm{1}}比Choi等^［23］提出的单元和Q4S单元^［21］收敛性更好， 但稍逊于Zhang等^［19］和Allman^［4］提出的单元.

上述分析结果充分表明， 即使与诸如混合有限元这些先进的单元理论相比， plane-Q4-TC单元也展现出了有竞争力的收敛速度.

2.4 工程应用

刚柔耦合定位平台物理模型如图10所示，主要由3部分组成：刚性框架、柔性铰链和工作平台.柔性铰链的两端分别与工作平台和刚性框架连接.该平台借助柔性铰链的弹性变形来进行微小运动，从而补偿摩擦死区实现纳米级定位，能够满足半导体行业的精密定位需求.结构设计对于平台的精度具有显著影响，而力学分析的作用不可或缺.因此为了提高刚柔耦合定位平台在结构设计上的准确性以提升分析的效率，本节以刚柔耦合定位平台为应用对象，验证平面plane-Q4-TC单元在工程应用中的实用性.

根据刚柔耦合定位平台的结构特点，提取其主视图截面作为本节的简单分析模型，其二维简化模型如图11所示.刚性框架的外侧施加固定约束， 工作台中心施加驱动力F_y .

将刚性框架和工作台使用不同单元尺寸（1，2，3，5，8 mm）的四边形网格离散，然后采用plane-Q4-TC单元建立有限元模型.柔性铰链使用梁单元建模，plane-Q4-TC单元能够与梁单元直接耦合，无须额外添加约束.作为比较，用经典等参元和ABAQUS S4R单元代替plane-Q4-TC单元进行分析.为了避免出现兼容性问题，经典等参元和ABAQUS S4R单元需要在连接节点处添加一些梁单元作为额外约束，具体的添加位置如图11所示.

为了评估3种单元的性能，用其求解平台的载荷点竖直方向位移，结果如表6所示.此外，图12描绘了随着网格的细化，3种单元求解的收敛趋势， 以便更加直观地比较不同单元的性能.同时，为确保分析结果的准确性，将整个二维模型（包括柔性铰链）用0.5 mm的四节点四边形网格进行离散，并使用ABAQUS S4R单元进行建模，旨在避免平面单元和梁单元之间的不兼容问题.然后将分析的结果作为图12的参考解加以对比验证.

从图12中可以看出，随着网格尺寸的减小， plane-Q4-TC单元可以稳定向参考解趋近. 与此相反， 另外两种单元的计算结果对单元尺寸变化的敏感度较低， 其变化趋势呈现出不规则性，表明含转角的平面单元具备更好的性能. 此外， 网格大小为 1 mm时， plane-Q4-TC单元与参考解之间的误差为0.54%， 而ABAQUS S4R单元和等参元的计算误差达到25.08%和24.57%， 表明plane-Q4-TC单元具有更高的计算精度. 由于ABAQUS S4R单元和等参元需要额外添加梁单元约束才能避免兼容问题， 这种处理实质上是一种粗糙的近似， 因此这种近似不可避免地对分析结果造成一定的误差.

3 结 论

本文基于“试探-修正”位移插值方案推导一种plane-Q4-TC单元.通过数值算例表明，该单元具有较快的收敛速度以及较高的数值精度.该单元还能够与平面梁单元直接耦合建模，无需额外的约束条件来确保兼容性，适用于工程实际.此外， “试探-修正”插值法可参数化，且具有良好的可扩展性，可作为推导包含真实转角的梁单元、壳单元和实体单元的理论基础，为相关研究提供参考.

最后，本文单元只涉及位移插值的改进，理论上边界条件和非线性材料本构模型不影响位移插值的有效性.此外，本文单元所使用的真实转角是基于小变形假设的位移梯度张量定义，因此不适用于处理几何非线性问题.今后的工作主要将小变形单元与多体系统动力学的浮动坐标法结合，也可对材料非线性和复杂问题边界进行扩展.

参考文献

[1] HUANG Y Q， YANG Y F， WANG J Z，et al. Unsymmetric extensions of Wilson’s incompatible four-node quadrilateral and eight-node hexahedral elements［J］. International Journal for Numerical Methods in Engineering，2022，123（1）：101-127．

[2] KARPIK A， COSCO F， MUNDO D. Higher-order hexahedral finite elements for structural dynamics： a comparative review［J］. Machines，2023，11（3）： 326．

[3] HUANG G X， LI H X， LU Y H，et al. Hexahedral solid element with rotational degrees of freedom based on a novel trail-correction displacement interpolation scheme［J］. Iranian Journal of Science and Technology，Transactions of Mechanical Engineering， 2024， 48（4）：1717-1730．

[4] ALLMAN D J. A quadrilateral finite element including vertex rotations for plane elasticity analysis［J］．International Journal for Numerical Methods in Engineering，1988，26（3）： 717-730．

[5] COOK R D. On the allman triangle and a related quadrilateral element［J］. Computers & Structures，1986，22（6）：1065-1067．

[6] SHANG Y， OUYANG W G ．4-node unsymmetric quadrilateral membrane element with drilling DOFs insensitive to severe mesh-distortion［J］．International Journal for Numerical Methods in Engineering，2018，113（10）：1589-1606．

[7] SHANG Y， CEN S， QIAN Z H，et al. High-performance unsymmetric 3-node triangular membrane element with drilling DOFs can correctly undertake in-plane moments［J］. Engineering Computations，2018， 35（7）： 2543-2556．

[8] SHANG Y， QIAN Z H， CEN S，et al. A simple unsymmetric 4-node 12-DOF membrane element for the modified couple stress theory［J］. International Journal for Numerical Methods in Engineering，2019，119（9）： 807-825．

[9] REZAIEE-PAJAND M， GHARAEI-MOGHADDAM N， RAME ZANI M. Two triangular membrane elements based on strain［J］. International Journal of Applied Mechanics，2019，11（1）：1950010．

[10] BOUTAGOUGA D. A formulation of membrane finite elements with true drilling rotation［J］. Engineering Computations， 2019， 37（1）： 203-236．

[11] DASZKIEWICZ K， WITKOWSKI W， BURZYŃSKI S，et al ．Robust four-node elements based on Hu-Washizu principle for nonlinear analysis of Cosserat shells［J］. Continuum Mechanics and Thermodynamics， 2019， 31（6）： 1757-1784．

[12] WANG A P. A quadrilateral membrane hybrid stress element with drilling degrees of freedom［J］. Acta Mechanica Sinica，2012， 28（5）： 1367-1373．

[13] 文颖， 戴公连， 曾庆元. 基于带面内转角自由度四节点平板壳单元的板壳非线性分析［J］. 中南大学学报（自然科学版），2013， 44（4）： 1525-1531．

[14] WEN Y， DAI G L， ZENG Q Y. 4-node flat shell element with drilling degrees of freedom for nonlinear analysis of plates and shells［J］. Journal of Central South University （Science and Technology）， 2013， 44（4）： 1525-1531．（in Chinese）

[15] SANGTARASH H， ARAB H G， SOHRABI M R，et al. A high-performance four-node flat shell element with drilling degrees of freedom［J］．Engineering with Computers，2021，37（4）：2837-2852．

[16] WU H P， SHANG Y， CEN S，et al ．Penalty C0 8-node quadrilateral and 20-node hexahedral elements for consistent couple stress elasticity based on the unsymmetric finite element method［J］. Engineering Analysis with Boundary Elements，2023，147： 302-319．

[17] LONG Y Q， XU Y. Generalized conforming quadrilateral membrane element with vertex rigid rotational freedom［J］. Computers & Structures， 1994， 52（4）： 749-755．

[18] BOUTAGOUGA D. A review on membrane finite elements with drilling degree of freedom［J］. Archives of Computational Methods in Engineering， 2021， 28（4）： 3049-3065．

[19] REZAIEE-PAJAND M， KARKON M ．An effective membrane element based on analytical solution［J］. European Journal of Mechanics - A/Solids，2013，39：268-279.

[20] ZHANG H X， KUANG J S ．Eight-node membrane element with drilling degrees of freedom for analysis of in-plane stiffness of thick floor plates［J］. International Journal for Numerical Methods in Engineering， 2008， 76（13）： 2117-2136．

[21] BERGAN P G， FELIPPA C A. A triangular membrane element with rotational degrees of freedom［J］. Computer Methods in Applied Mechanics and Engineering， 1985， 50（1）： 25-69．

[22] MACNEAL R H， HARDER R L. A refined four-noded membrane element with rotational degrees of freedom［J］. Computers & Structures， 1988， 28（1）： 75-84.

[23] PIMPINELLI G. An assumed strain quadrilateral element with drilling degrees of freedom［J］. Finite Elements in Analysis and Design， 2004， 41（3）： 267-283.

[24] CHOI N， CHOO Y S， LEE B C. A hybrid Trefftz plane elasticity element with drilling degrees of freedom［J］. Computer Methods in Applied Mechanics and Engineering， 2006， 195（33-36）： 4095-4105.

基金资助

国家重点研发计划项目（2022YFB4701001）