目前，针对梁扭转的计算，一般以薄壁箱梁为主要研究对象。当薄壁截面发生扭转时，截面会产生纵向的翘曲，根据翘曲有无约束^[1]，将扭转分为两种类型：自由度扭转(也称圣维南扭转)和约束扭转。对于标准截面(等厚度的“圆管截面”和“正多边形”截面)，一般认为没有畸变^[2]；但在实际工程中，梁板块的横向弯曲等因素会使截面产生翘曲^[3]，非标准截面又会使翘曲受到约束，每块板进一步在垂直于板平面的横向平面内发生弯曲^[3]，最终导致截面产生畸变。且大量研究表明^{[2, 4]}：翘曲、畸变效应不能忽略，因此合理考虑翘曲、畸变的影响，对于类似于桥梁这样的结构设计、施工、监测与检测有重要意义。

目前关于翘曲畸变的计算方法，主要分为解析理论和数值理论。解析理论：主要包括弹性理论、乌曼斯基的第一理论和第二理论、詹涅里杰和巴诺夫柯理论、符拉索夫理论^[5]等。孙训方等^[6]基于弹性理论和实验现象，认为实心矩形截面承受扭矩时，角点处剪应力几乎为零，最大剪应力发生在长边或短边中点处。强士中等^[7]分析乌曼斯基的第一理论和第二理论等，认为不同的假定条件和处理方式，是导致不同理论的计算结果存在显著差异的主要原因。在数值有限元分析中，杆单元(或梁单元)因其精度高、求解速度快、适用性广，备受关注。但由于传统梁单元^[8](单结点6自由度)通常忽略了翘曲、畸变效应，使得数值结果存在一定的计算误差，于是为了提高计算精度，国内外很多学者持续展开了“构造梁单元”的研究。例如：包世华等^[2]在Bernoulli-Euler梁单元的基础上，增设了第7自由度扭率。杨绿峰等^[9]对扭转角和扭率进行3次Hermite插值，构造了单结点7自由度的梁单元。钟新谷等^[10]对梁单元的弯曲、扭转、翘曲均进行3次插值，构造了单结点7自由度的梁单元。韦成龙等^[11]提出一种28自由度(单结点14自由度)，考虑了翘曲、畸变、剪力滞效应的梁单元。谢旭等^[12]通过引入扭率、畸变角、畸变率3个自由度，提出一种新的27个自由度(单结点9自由度)梁单元。张叔辉^[13]通过引入翘曲、畸变、畸变率、上翼缘剪切位移差4个自由度，提出一种新的20个自由度(单结点10自由度)梁单元。虽然上述梁单元考虑了截面翘曲、畸变的影响，在一定程度上弥补了Bernoulli-Euler梁单元约束扭转的计算缺陷，但在应用时仍存在一定的局限性，例如：

1) 积分精度较低：由于特殊函数(例如对数、指数、分式等^{[1 − 13]})的存在，使得单元的刚度、非结点荷载、质量等矩阵的弱势积分十分复杂，最终导致单元的数值积分精度下降；且由于单一结点自由度数目与其他单元不同(例如Q4薄膜元、R8实体单元等)，共结点处的精度难以保证；

2) 适用范围较低：因为上述理论均以开口或闭口的薄壁箱梁为主要研究对象，使得参数具备一定的特殊性，导致理论难以涵盖所有截面类型，例如实心截面、变截面、非对称截面等；

3) 求解效率较低：由于增设了一个或多个自由度，使得矩阵的维度不断扩大，导致方程组解算的效率呈几何倍数的降低；在同等条件下相较于传统的梁单元(单结点6自由度)，会耗费计算机大量的存储空间。

4) 编程难度较大：由于增设了一个或多个自由度，使得有限元程序中不同单元之间，单一结点的自由度数目不同，增大了不同单元共结点处自由度处理方法的编程难度。

因此，针对上述不足，本文对前期得到的平面两结点双线性插值梁单元(Bm_ZY)^[14]进行了空间拓展，较合理地考虑了多个自由度之间的耦合关系，并通过近似函数简化单元积分，尝试构造了一类考虑翘曲、畸变的单结点6自由度低次梁单元；最后借助于自编的Python有限元程序(ZQFEM)进行计算分析，并和Midas Civil商业软件计算结果对比，验证了所构造的单元位移模式的正确性及可靠性。

1   空间两结点双线性插值梁单元(Bm_ZY)

1.1   基本假定与坐标定义

在空间梁单元上建立局部坐标系，如图1所示。以梁的中轴线为x轴，在梁单元的截面处建立正交坐标系yoz，并引入如下基本假定(本文所有符号的含义见附表1)：

图  1  梁单元坐标系

Figure  1.  Beam element coordinate system

下载: 全尺寸图片 幻灯片

1) 梁的变形始终处于线弹性范围内，且变形对梁自身力学性质的影响忽略不计(“弹性小变形”假定)；

2) 梁单元的位移模式和应变模式严格满足《弹性力学》几何方程和本构关系^[15]；

3) 梁纵向纤维的伸缩、截面翘曲，由截面形状、面内外自由度耦合等多种因素产生。

1.2   构造位移模式和应变模式

以u、v、w分别表征单元x、y、z方向的位移。由《弹性力学简明教程》可知^[15]，任意单元上的几何方程为：

式(1)和式(2)中，变量γ和θ分别为剪应变和转动矢量(也称“旋转度”)。当截面自由扭转时，由《材料力学》可得到扭转的几何关系^[6]：

式中，$R = R(y,{\textit{z}}) = \sqrt {{y^2} + {{\textit{z}}^2}}$。

联立式(1)~式(3)可得：

从式(4)可以看出，截面处的位移v、w与扭转角θ_yz、扭率θ_yz,x有关。当截面扭转时，截面各点处会存在不同的横向位移，位移差导致弯曲角θ_zx、θ_xy发生变化，并进一步产生横向弯曲，最后使得截面发生翘曲、畸变现象。为了便于后续的单元积分，现尝试牺牲部分截面特性，将式(4)进行简化，得到：

进一步引入z方向自由度，得到非耦合位移模式^[14]：

引入形函数N矩阵和几何B矩阵，可导出梁单元的位移向量δ和应变分量：

式中：${{\boldsymbol{q}}_i} = {\left[ {{{\boldsymbol{u}}_i}}\;\;\;{{{\boldsymbol{v}}_i}}\;\;\;{{{\boldsymbol{w}}_i}}\;\;\;{{{\boldsymbol{\theta}} _{y{\textit{z}}i}}}\;\;\;{{{\boldsymbol{\theta}} _{{\textit{z}}xi}}}\;\;\;{{{\boldsymbol{\theta}} _{xyi}}} \right]^{\mathrm{T}}}$，i为局部结点号；${\boldsymbol{ \varepsilon }}= {\left[ {{{\boldsymbol{\varepsilon}} _x}}\;\;\;{{{\boldsymbol{\varepsilon}} _y}}\;\;\;{{{\boldsymbol{\varepsilon}} _{\textit{z}}}}\;\;\;{{{\boldsymbol{\gamma}} _{y{\textit{z}}}}}\;\;\;{{{\boldsymbol{\gamma}} _{{\textit{z}}x}}}\;\;\;{{{\boldsymbol{\gamma}} _{xy}}} \right]^{\mathrm{T}}}$。

由此得到形函数N矩阵(位移模式)和几何B矩阵(应变模式)的显示表达：

1.3   单元积分

由单元的本构关系得到D矩阵^{[14 − 17]}：

最后通过弱式积分，并结合积分表(附表2)，得到空间Bm_ZY单元的刚度矩阵。由于公式的复杂性，这里仅给出实心圆形变截面的轴向拉压刚度和扭转刚度(i、j分别表示刚度矩阵的行号、列号)：

在求解刚度矩阵时，也可通过“高斯积分”对式(4)所表征的非简化的位移模式进行计算，但本文不再赘述。

2   算例

为验证空间Bm_ZY单元的可靠性，且便于成果的复现，本文基于自编的ZQFEM有限元程序进行了算例验证。不妨考虑一“悬臂梁模型”，并定义圆形、矩形、正方形等三种不同的截面外轮廓，并保持内外轮廓形状相似、形心重合。定义模型参数：梁长L=1 m，弹性模量E=1 Pa，泊松比µ =0，扭矩M_yz=1 N·m，矩形宽b=0.03 m，矩形高h=0.04 m，圆形半径r=0.01 m，正方形边长a=0.02 m。

如表1所示，为实心圆形截面、矩形截面、正方形截面扭转角的计算结果，其中“Python_初等梁^[14]”为Python内置的Bernoulli-Euler梁单元；“Midas Civil_翘曲”表示使用的梁单元开启了第七自由度；“Midas Civil_非翘曲”表示使用的梁单元关闭第七自由度；“《材料力学》^[6]”表示自由扭转的解析解计算值。表1以计算值与自由扭转计算结果的相对误差，表征单元截面翘曲畸变的影响，同类截面相对误差的绝对值越大，翘曲畸变的影响也越大。保持截面外轮廓不变，设置一系列截面内轮廓的尺寸梯度，不难获得如图2所示不同类型截面的扭转角与壁厚的变化规律。

图  2  空心截面扭转角与壁厚曲线图

Figure  2.  A Curve of torsional Angle and wall thickness of hollow section

下载: 全尺寸图片 幻灯片

观察Midas Civil的计算结果。从表1中可以发现无论是否“开启第七自由度(翘曲自由度)”，实心矩形和实心正方形截面的差异值均随着截面的不同而不断变化。说明Midas Civil中梁的单元刚度矩阵已经考虑了扭转自由度与其他自由度的耦合关系，且耦合关系以刚度修正系数的形式内置于梁单元的刚度矩阵中(类比剪切修正系数k_s^[8])，开启第七自由度的作用仅为“是否让翘曲自由度参与方程组的解算”。从图2中可以发现，当开启第七自由度时，截面的扭转刚度显著增大；且随着壁厚的不断增大，截面的扭转刚度也不断增大。当关闭第七自由度时，空心正方形截面的扭转角呈现明显的不连续性，且矩形截面与正方形截面的相对误差，均在壁厚0.09处出现不连续点。开启第七自由度前后，正方形或矩形截面的差异值，呈相反的变化趋势。

观察Bm_ZY单元的计算结果。对于圆形截面和矩形截面，相对误差与截面的壁厚无关，与截面类型有关，分别为−51.35%、−57.14%，与文献[2]的结论相吻合；对于矩形截面，随着壁厚的增大，扭转刚度也不断的增大，差异值从−59.42%逐渐变化到−60.22%，翘曲畸变的影响也逐渐增大。因此，相较于Midas Civil，本文所采用的空间Bm_ZY单元，具有运算效率较高(自由度少)、翘曲畸变的影响较合理等优点。

表  1  实心截面扭转角计算结果

Table  1.  Calculation results of torsional angle of solid belly section

截面形状	Midas Civil_翘曲	Midas Civil_非翘曲	Python_初等梁[14]	《材料力学》[6]	Bm_ZY
实心圆形 (相对误差)	127 323 954.474 (+0.00%)	127 323 954.474 (+0.00%)	127 323 954.474 (+0.00%)	127 323 954.474 (+0.00%)	61 941 383.257 (−51.35%)
实心矩形 (相对误差)	9 121 250.417 (+14.02%)	10 288 856.640 (+28.61%)	8 000 000.000 (+0.00%)	8 000 000.000 (+0.00%)	3 182 634.608 (−60.22%)
实心正方形 (相对误差)	79 038 701.805 (+5.38%)	88 888 888.889 (+18.51%)	75 000 000.000 (+0.00%)	75 000 000.000 (+0.00%)	32 142 857.143 (−57.14%)

下载: 导出CSV 

| 显示表格

3   结论

本文基于平面Bm_ZY单元，构造了考虑翘曲畸变的6自由度的空间Bm_ZY单元，得到以下结论：

(1) 相较于其他梁单元，Bm_ZY单元具有积分简单、矩阵维数低、计算结果趋势较为可靠、适用任何截面等优点。

(2) 截面翘曲畸变效应，会显著增大截面的扭转刚度，且同类截面随着壁厚的增大，扭转刚度也会不断增大。

(3) 圆形截面和正方形截面无翘曲畸变效应；矩形截面翘曲畸变的影响，随着截面壁厚呈的增大而逐渐增大。

从表1和图2可以看出，本文的Bm_ZY单元扭转效应的计算结果，与商业软件Midas Civil差距较大；文献[17]虽然指出“截面方向离散程度不够，是造成Bm_ZY单元误差较大的主要原因”，但从式(14)和式(15)可以看出，截面的泊松效应和弹性模量，也会显著影响单元的计算结果。因此，进行梁的高精度足尺实验，并探究泊松效应、弹性模量两种常数的独立性，是作者接下来的主要工作。

附录1 ：名词解释

  1  名词解释

  1.  Explanation of nouns

符号	含义		符号	含义		符号	含义		符号	含义
γyz	局部坐标系yoz剪应变		θyz	局部坐标系yoz面转角		“,”(下角标)	偏导简化符号		r	截面半径(外轮廓)
γzx	局部坐标系xoz剪应变		θzx	局部坐标系xoz面转角		A	梁单元截面积		h	截面高(外轮廓)
γxy	局部坐标系xoy剪应变		θxy	局部坐标系xoy面转角		L	梁单元长度		b	截面宽(外轮廓)
εx	局部坐标系x轴向应变		u	局部坐标系x轴向位移		µ	泊松比		E	弹性模量
εy	局部坐标系y轴向应变		v	局部坐标系y切向位移		U2	弹性常数[15]		U1	弹性常数[15]
εz	局部坐标系z轴向应变		w	局部坐标系z切向位移		D	本构矩阵		B	几何矩阵
ks	剪切修正系数(剪切系数)		kA	变截面参数		q	位移列阵		N	形函数矩阵
ε	应变矩阵		−	−		−	−		−	−

下载: 导出CSV 

| 显示表格

附录2 ：积分表

  1  矩形截面积分示意图

  1.  Integral diagram of rectangular section

下载: 全尺寸图片 幻灯片

将附图1所示矩形区域内的重积分$\iint {f(y,{\textit{z}})}{\mathrm{d}}y{\mathrm{d}}{\textit{z}}$，转化为极坐标系的积分$\iint {{r^i}{{\sin }^j}\theta {{\cos }^k}\theta }{\mathrm{d}}r{\mathrm{d}}\theta = {h^{i + 1}}{\eta _{i {\text{-}} j {\text{-}} k}}$。附表2为附图1中一象限的积分结果。

  2  积分表

  2.  Table of integration

i-j-k	ηi-j-k
1-0-0	$\dfrac{{{b^2}{k_1} + {h^2}{k_2}}}{8}\dfrac{1}{{{h^2}}}$
2-0-0	$\dfrac{{{b^3}[\ln (k{k_1} + {k_1}) + {k_1}k{k_1}] + {h^3}[\ln (k{k_2} + {k_2}) + {k_2}k{k_2}]}}{{48}}\dfrac{1}{{{h^3}}}$
2-2-0	$\dfrac{{{b^3}[ - \ln (k{k_1} + {k_1}) + {k_1}k{k_1}] + 2{h^3}\ln (k{k_2} + {k_2})}}{{48}}\dfrac{1}{{{h^3}}}$
2-0-2	$\dfrac{{2{b^3}\ln (k{k_1} + {k_1}) + {h^3}[ - \ln (k{k_2} + {k_2}) + {k_2}k{k_2}]}}{{48}}\dfrac{1}{{{h^3}}}$
3-0-0	$\dfrac{{{b^4}(k{k_1^3}/3 + {k_1}) + {h^4}(k{k_2^3}/3 + {k_2})}}{{64}}\dfrac{1}{{{h^4}}}$
3-2-0	$\dfrac{{{b^4}k{k_1^3}/3 + {h^4}{k_2}}}{{64}}\dfrac{1}{{{h^4}}}$
3-0-2	$\dfrac{{{b^4}{k_1} + {h^4}k{k_2^3}/3}}{{64}}\dfrac{1}{{{h^4}}}$
3-2-2	$\dfrac{{{b^4}({k_1} - \arctan {k_1}) + {h^4}({k_2} + \arctan {k_1})}}{{64}}\dfrac{1}{{{h^4}}}$
3-4-0	$\dfrac{{{b^4}({k_1^3}/3 - {k_1} + \arctan {k_1}) + {h^4}(\pi /2 - \arctan {k_1})}}{{64}}\dfrac{1}{{{h^4}}}$
3-0-4	$\dfrac{{{b^4}\arctan {k_1} + {h^4}({k_2^3}/3 - {k_2} - \arctan {k_1})}}{{64}}\dfrac{1}{{{h^4}}}$
4-2-0	$\dfrac{{{b^5}[ - 5\ln (k{k_1} + {k_1}) + 3{k_1}k{k_1} + 2k{k_1^3}{k_1}] + 4{h^5}[\ln (k{k_2} + {k_2}) + {k_2}k{k_2}]}}{{1024}}\dfrac{1}{{{h^5}}}$
4-0-2	$\dfrac{{4{b^5}[\ln (k{k_1} + {k_1}) + {k_1}k{k_1}] + {h^5}[ - 5\ln (k{k_2} + {k_2}) + 3{k_2}k{k_2} + 2k{k_2^3}{k_2}]}}{{1024}}\dfrac{1}{{{h^5}}}$
4-4-0	$\dfrac{{{b^5}[3\ln (k{k_1} + {k_1}) + 5{k_1^3}k{k_1} - 3{k_1}k{k_1^3}] + 8{h^5}\ln (k{k_2} + {k_2})}}{{1024}}\dfrac{1}{{{h^5}}}$
4-0-4	$\dfrac{{8{b^5}\ln (k{k_1} + {k_1}) + {h^5}[3\ln (k{k_2} + {k_2}) + 5{k_2^3}k{k_2} - 3{k_2}k{k_2^3}]}}{{1024}}\dfrac{1}{{{h^5}}}$
5-2-0	$\dfrac{{{b^6}(k{k_1^5}/5 + {k_1^3}/3) + {h^6}(k{k_2^3}/3 + {k_2})}}{{320}}\dfrac{1}{{{h^6}}}$
5-0-2	$\dfrac{{{b^6}(k{k_1^3}/3 + {k_1}) + {h^6}(k{k_2^5}/5 + {k_2^3}/3)}}{{320}}\dfrac{1}{{{h^6}}}$

下载: 导出CSV 

| 显示表格

其中，参数k₁、kk₁、k₂、kk₂为Python代码中的变量，满足下式：