梁单元无疑是有限元分析中应用最为广泛的单元之一。为了满足各种工程需求,梁单元常常可能要求具备以下特性或功能：

(1) 剪切变形;

(2) 端头刚域,梁-柱交汇处的区域;

(3) 部分刚接,梁-柱交汇等处的刚度损失;

(4) 端头铰接,梁-柱交汇等处的构造或刚度缺失;

(5) 截面变化,例如桥梁等大跨结构;

(6) 部分偏心,例如与墙柱相连的外立面梁;

(7) 几何刚度,轴向力引起的附加刚度;

(8) 轴线为曲线,例如：拱。

在中外文献中关于梁单元的形成方法层出不穷,这里不一一列举,但大部分文献或只考虑某一特定的单元形式,或形成方案过于繁琐,不宜在工程软件中获得应用。本文从工程实用的角度出发,着重讨论前5类功能在长为$L$的二维平面梁单元中的统一实现,读者不难将其推广至三维情形。

以Timoshenko梁(T $\!-\!$梁)讨论刚度矩阵以及载荷向量的形成方案,Euler梁(E $\!-\!$梁)可以作为T $\!-\!$梁的一个特例。为了处理变截面梁、端头刚域以及端头部分刚接,我们以柔度法为讨论重点。

考虑一个置于$Oxy$平面内长为$L$的梁单元,其主轴与$x$轴重合,左端与原点$O$重合。与标准的有限元教材^[1-2]稍许不同,考虑到单元的通用性,工程软件中梁单元的$4\times4$刚度矩阵(这里只考虑弯曲部分)往往不是按相应的公式直接形成的,而是通过简支梁或悬臂梁的$2\times2$刚度矩阵扩张得到的。为此,对照梁单元与同跨度的简支梁(图1),它们的挠度和转角满足变换

$$ \begin{eqnarray*} &&w(\xi )=w_{\rm s} (\xi )+(1-\xi )w_{i} +\xi w_{j} \\ && \theta (\xi )=\theta_{\rm s} (\xi )+\frac{w_{j} -w_{i} }{L} \end{eqnarray*} $$

其中$w(\xi )$, $\theta (\xi )$, $w_{\rm s} (\xi )$, $\theta_{\rm s} (\xi)$分别以自然坐标$\xi=x/L$表示梁单元和相应简支梁的挠度与转角。 $w(\xi )$, $\theta (\xi)$满足T $\!-\!$梁的控制方程,当且仅当$w_{\rm s} (\xi )$, $\theta_{\rm s} (\xi )$满足该方程。记${\mathbf u}^{(e)}=\{w_{i} ,\theta_{i} ,w_{j} ,\theta_{j} \}$为梁单元两端的挠度与截面转角,${\mathbf u}_{\rm s}^{(e)} =\{\alpha_{i} ,\alpha _{j} \}$为简支梁两端截面转角,则有

其中,$\mathbf S$是两组结点位移之间的几何关系矩阵,其转置是相应的两端外力关系矩阵,即

$$ \begin{eqnarray*} {\mathbf f}^{(e)}={\mathbf S}^{\rm T}{\mathbf f}_{\rm s}^{(e)} \end{eqnarray*} $$

其中,${\mathbf f}^{(e)}=\left\{{Q_{i} }\ \ {M_{i} }\ \ {Q_{j} } \ \ {M_{j} } \right\}$是梁单元两端的作用力,而${\mathbf f}_{\rm s}^{(e)} =\left\{{M_{i} } \ \ {M_{j} } \right\}$是相应简支梁两端的作用力矩,$Q_{i} =-Q_{j} =(M_{i} +M_{j} )/L$反映了梁端作用力的平衡关系。由此可知,通过这两个关系矩阵,任意梁单元的变形和梁端作用力均可由一个相应的简支梁表达出来,当然 我们也可以由简支梁得到普通梁单元的刚度矩阵与等效载荷向量。

记$\xi=x/L$为梁单元的自然坐标,($*$)$'$ 为($*$)对$\xi$的导数。由于简支梁是静定结构,我们可以通过考虑剪切效应的最小余能原理求简支T $\!-\!$梁的柔度矩阵

其中的弯矩取线性插值

满足梁平衡关系。对于等截面梁,简支T $\!-\!$梁的柔度矩阵可解析给出

其中无量纲参数$\eta ={12EI}/({GA_{\rm s}L^{2}})$,表示弯曲与剪切刚度之比。其刚度矩阵可由式(4)求逆获得

另一方面,我们还可以利用最小势能原理获得梁单元与相应简支梁单元的刚度矩阵变换关系以及简支梁单元的刚度矩阵^[3]。记$H_{2} (\xi )=\xi (\xi -1)^{2}$, $H_{4} (\xi )=\xi^{2}(\xi -1)$为^{[0, 1]}上标准的三次Hermite插值函数中关于导数的两个,则满足无分布力简支T $\!-\!$梁方程的位移插值形函数可以写为

其中

$$ \begin{eqnarray*} &&\mathbf N_{\rm s}=\left[ {{\begin{array}{*{20}c} {N_{2} } \\ {N_{4} } \\ \end{array} }} \right]=\frac{1}{1+\eta }\left\{ {\left[ {{\begin{array}{*{20}c} {H_{2} } \\ {H_{4} } \\ \end{array} }} \right]+\frac{\eta }{2}\left[ {{\begin{array}{*{20}c} {\xi (1-\xi )} \\ {-\xi (1-\xi )} \\ \end{array} }} \right]} \right\}\\ &&\mathbf M_{\rm s}=\left[ {{\begin{array}{*{20}c} {M_{2} } \\ {M_{4} } \\ \end{array} }} \right]=\frac{1}{1+\eta }\left\{ {\left[ {{\begin{array}{*{20}c} {{H}'_{2} } \\ {{H}'_{4} } \\ \end{array} }} \right]+\eta \left[ {{\begin{array}{*{20}c} {1-\xi } \\ \xi \\ \end{array} }} \right]} \right\} \end{eqnarray*} $$

基于最小势能原理的单元刚度矩阵为

由式(7)中的积分,可以获得与式(5)完全相同的刚度矩阵。式(4)和式(5)互为逆,表明两端截面转角、力矩和插值所得的跨内挠度、弯矩之间的准确关系。

当$\eta =0$时,剪切刚度无穷大,T $\!-\!$梁回归到E $\!-\!$梁。作为$\xi $的函数, ${\mathbf N}_{\rm s} $就是Hermite插值$\mathbf H$,${\mathbf M}_{\rm s}$是其微商${\mathbf {H}'}$。此时,式(5)也退化到E $\!-\!$梁的刚度矩阵。

当截面沿中性轴变化时,位移的三次Hermite插值形函数不再满足梁的平衡方程,因此最小势能原理所导出的刚度矩阵不再是简支梁的精确刚度矩阵。但作为静定结构,简支梁的弯矩分布仍然是线性的,插值函数式(3)仍然给出结构力学意义下的精确内力,式(2)是精确的表达式。由此,最小余能原理所导出的柔度矩阵还是简支梁的精确柔度矩阵。在这种意义下,变截面梁单元的刚度矩阵可以从柔度矩阵求逆而准确得到。

式(2)一般可通过较高阶的Gauss积分或Lobatto积分获得。我们的经验是5点Gauss积分的精度已足够工程应用。变截面梁的刚度矩阵如采用位移法得到,则可能隐含较大的误差^[4]。

在结构力学中梁-柱的计算是所谓的中心线计算。在实际工程中,梁在其与柱的交汇区域内几乎没有弯曲变形,这一部分称为刚域(图2),其长度一般可取为柱的中线到柱外缘距离的75%。

一旦计入不变形的刚域,梁的弹性变形部分缩短,刚度增大。如图3,仍记梁柱中心线处简支梁的转角为${\mathbf u}_{\rm s}^{(e)} =\{\alpha_{i} ,\alpha_{j}\}$,又记刚域外缘弹性梁的转角为${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf u}} }_{\rm s}^{(e)} =\{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\alpha }}_{i} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\alpha }}_{j} \}$,那么由线性化的几何关系得到

现在,长为$L-c-d$的弹性简支梁两端的截面转角${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf u}} }_{\rm s}^{(e)} $与外力矩${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf f}} }_{\rm s}^{(e)} $满足

$$ \begin{eqnarray*} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf u}} }_{\rm s}^{(e)} ={\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf F}} }_{\rm s}^{(e)} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf f}} }_{\rm s}^{(e)} \end{eqnarray*} $$

其中弹性部分的柔度矩阵${ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf F}} }_{\rm s}^{(e)} $按式(2)计算,但长度替换为$L-c-d$。由刚域外缘与中心线的力矩平衡得出

综合起来,中心线位置的简支梁柔度矩阵为

梁端有时与柱等其他构件的连接不紧密,不能完全传递力矩,此时一般称为部分刚接,它也可在弹塑性计算中模拟理想塑性铰。在有限元计算中,部分刚接用梁端的附加扭转弹簧来模拟。

较为方便地获得刚度矩阵的方法仍然是利用最小余能原理,修改不含转角弹簧的柔度矩阵${ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf F}} }_{\rm s}^{(e)} $。例如,用最小余能原理得到图4所示含右端转动弹簧的柔度矩阵其逆即为相应的刚度矩阵。

对于等截面直梁,当$\dfrac1k \to 0$时,${\mathbf K}_{\rm s}^{(e)} \to \dfrac{3EI}{L}\cdot$ $\left[ {{\begin{array}{*{20}c} 1 & \\[-2mm] & 0 \\ \end{array} }} \right]$,表示右端不能传递力矩,退化为完全铰接;当$\dfrac1k\to +\infty $时,${\mathbf K}_{\rm s}^{(e)} \to \dfrac{2EI}{L}\left[ {{\begin{array}{*{20}c} 2 & 1 \\ 1 & 2 \\ \end{array} }} \right]$,退化为完全刚接,与1.2节得到的刚度相同。

由以上讨论可以得到梁单元刚度矩阵形成方式：首先计算弹性简支梁的柔度矩阵${ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf F}} }_{\rm s}^{(e)} $,如果存在刚域或部分刚接的情形,则需进行端头刚域与集中弹簧修正,获得${\mathbf F}_{\rm s}^{(e)} $,最后通过变换${\mathbf K}^{(e)}={\mathbf S}^{\rm T}\left[ {{\mathbf F}_{\rm s}^{(e)} } \right]^{-1}{\mathbf S}$计算梁单元的刚度矩阵。

梁载荷的等效结点力大约是工程有限元软件中最为复杂的部分。由于工程实际需求,必须要考虑多种形式的载荷,例如：集中力、均布载荷、三角形分布载荷、梯形分布载荷等,且载荷作用在纯弹性梁、组合刚域-弹性梁、组合弹簧-弹性梁。按梁的形式分工况形成载荷无疑是一种可能性,但我们在这里再一次从简支梁出发,给出一个简洁的方案。限于篇幅,仅叙述集中力的作用下的等效载荷,其他分布载荷情况可考虑为单位集中力的等效结点力的积分。

假设集中力$p$作用在梁上,距$i$端的距离为$a$,记$b=L-a$,用力的平衡关系可得简支梁跨内集中力作用在两端支座上的等效作用力

用柔度法可得两端的转角

其中

是简支梁在集中力$p$作用下的弯矩。对于等截面简支梁^[5]可得到

$$ \begin{eqnarray*}\%1 \left[ {{\begin{array}{*{20}c} {\alpha_{i} } \\ {\alpha_{j} } \\ \end{array} }} \right]=\frac{p}{6LEI}\left[ {{\begin{array}{*{20}c} {ab(L+b)} \\ {-ab(L+a)} \\ \end{array} }} \right] \end{eqnarray*} $$

如果是变截面梁,式(12)的数值积分需要按函数的光滑性分段进行。

用式(13)跨内载荷产生的转角,可以得到跨内载荷在两端的等效作用力矩

同时按照平衡关系,这两个等效力矩也会带来在支座上的附加作用力,应用式(1)中的转换矩阵${\mathbf S}$,我们获得最终等效结点载荷

对于等截面梁,等效结点力的表达式是

如不考虑$\eta $,即剪切刚度无穷大,式(17)退化为等截面E $\!-\!$梁的等效结点力。

上面的推导过程中,我们没有区分是否是等截面梁。如果是变截面梁,封闭形式的等效结点力一般难于得到,在程序中建议用分段数值积分计算。

如果要考虑带刚域梁的等效结点力,则需要分别处理作用在弹性部分与刚域部分的跨内载荷。我们首先用前述的方法计算出作用在弹性部分载荷的等效结点力${ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\mathbf p}} }^{(e)}$,沿用图3的记号,记中心线的等效结点力为${\mathbf p}^{(e)}$,用平衡即可得到它们之间的转换关系

如果集中力$p$作用在左、右端的刚域上,即$a\leqslant c$或$b\leqslant d$,那么等效集中力仅出现在一端,分别为

在实际计算中并不需要显式地采用式(17),而是通过计算简支梁的支反力与转角间接获得等效结点力,然后再实施关于刚域的变换。

特别注意,部分刚接或杆端扭转弹簧对简支梁的反力与转角计算在目前的方案中均无直接影响,它们对等效结点载荷的影响间接地通过刚度矩阵的修改计入。当然我们假定,梁的一端不同时出现矛盾的部分刚接与刚域。

如图5所示,假设我们已经通过求解有限元方程获得了梁两端的位移${\mathbf u}^{(e)}=\{w_{i} ,\theta_{i} ,w_{j} ,\theta_{j} \}$,则两端的真实作用力可按产生结点位移所必需的作用力与跨间载荷的等效结点力叠加计算,

$$ \begin{eqnarray*} {\mathbf f}^{(e)}={\mathbf K}^{(e)}{\mathbf u}^{(e)}-{\mathbf p}^{(e)}=\left[ {{\begin{array}{*{20}c} {Q_{i} } \\ {M_{i} } \\ {Q_{j} } \\ {M_{j} } \\ \end{array} }} \right] \end{eqnarray*} $$

其中${\mathbf p}^{(e)}$是等效结点力^[6]。对于梁单元,这样得到的端头作用力与端头内力之间仅有符号的差异。内力输出时应按工程习惯进行正负号 的调整。无论是否含端头刚域与弹簧,跨内分布横向载荷$q(x)$的梁内力用梁段的平衡关系给出

$$ \begin{eqnarray*} &&Q(x)=Q_{i} +\int_0^x {q(s){\rm d}s}\\ &&M(x)=M_{i} -Q_{i} x-\int_0^x {(x-s)q(s){\rm d}s} \end{eqnarray*} $$

原则上从柔度法可以得到简支梁两个弯曲方向刚度矩阵${\mathbf K}_{{\rm s,1\mbox{-}2}}^{(e)}$, ${\mathbf K}_{{\rm s,1\mbox{-}3}}^{(e)} \in {\bf R}^{2\times 2}$,其中1-2代表主弯曲平面,1-3代表次弯曲平面。除了弯曲之外,梁单元还要考虑沿中性轴向的固定-自由杆拉压刚度${\mathbf K}_{11}^{(e)} \in {\bf R}^{1\times 1}$与绕中性轴的固定-自由杆扭转刚度${\mathbf K}_{44}^{(e)} \in {\bf R}^{1\times 1}$,合在一起就是

$$ \begin{eqnarray*} {\tilde{{\mathbf K}}}^{(e)}=\mbox{diag}({\begin{array}{*{20}c} {{\mathbf K}_{{\rm s,1\mbox{-}2}}^{(e)} } & {{\mathbf K}_{{\rm s,1\mbox{-}3}}^{(e)} } & {{\mathbf K}_{11} } & {{\mathbf K}_{44} } \\ \end{array} })\in {\bf R}^{6\times 6}, \end{eqnarray*} $$

各方向定义在图6中给出。

为了获得整体坐标系下的梁(弯-拉-扭)的$12\times12$单元刚度矩阵,需要经过两次变换,第一次是将整体的位移$\{u_{x}^{(i)}$ $u_{y}^{(i)}$ $u_{z}^{(i)}$ $\theta_{x}^{(i)}$ $\theta_{y}^{(i)}$ $\theta _{z}^{(i)}$ $u_{x}^{(j)}$ $u_{y}^{(j)}$ $u_{z}^{(j)}$ $\theta_{x}^{(j)}$ $\theta_{y}^{(j)}$ $\theta_{z}^{(j)} \}^{\rm T}$通过坐标变换转换到梁局部坐标系中$\{u_{1}^{(i)}$ $u_{2}^{(i)}$ $u_{3}^{(i)}$ $\theta_{1}^{(i)}$ $\theta_{2}^{(i)}$ $\theta_{3}^{(i)}$ $u_{1}^{(j)}$ $u_{2}^{(j)}$ $u_{3}^{(j)}$ $\theta_{1}^{(j)}$ $\theta _{2}^{(j)}$ $\theta_{3}^{(j)} {\}}^{\rm T}$,其中$\{i;{\mathbf n}_{1}$, ${\mathbf n}_{2}$, ${\mathbf n}_{3} \}$是梁的右手局部坐标系,${\mathbf n}_{1} $是从$i$到$j$的方向,${\mathbf n}_{2} $是次弯曲轴方向,${\mathbf n}_{3} $是主弯曲轴方向。第二次是将局部坐标系的位移变换到"简支"坐标系的位移,它的形式是

$$ \begin{eqnarray*} \left[ {{\begin{array}{*{20}c} {\alpha_{1-2}^{(i)} } \\ {\alpha_{1-2}^{(j)} } \\ {\alpha_{1-3}^{(i)} } \\ {\alpha_{1-3}^{(j)} } \\ {u^{(j-i)}} \\ {\theta^{(j-i)}} \\ \end{array} }} \right]=\left[ {{\begin{array}{*{20}c} & 1/L & & & & 1 & & -1/L & & & & \\ & 1/L & & & & & & -1/L & & & & 1 \\ & & -1/L & & 1 & & & & 1/L & & & \\ & & -1/L & & & & & & 1/L & & 1 & \\ {-1} & & & & & & 1 & & & & & \\ & & & {-1} & & & & & & 1 & & \end{array} }} \right]\cdot\left[ {\begin{array}{l} u_{1}^{(i)} \\ u_{2}^{(i)} \\ u_{3}^{(i)} \\ \theta_{1}^{(i)} \\ \theta_{2}^{(i)} \\ \theta_{3}^{(i)} \\ u_{1}^{(j)} \\ u_{2}^{(j)} \\ u_{3}^{(j)} \\ \theta_{1}^{(j)} \\ \theta_{2}^{(j)} \\ \theta_{3}^{(j)} \end{array}} \right] \end{eqnarray*} $$

梁的刚度矩阵与等效结点载荷的形成方案要考虑多方面的因素,如图7所示。

梁单元刚度矩阵的形成方案按以下顺序进行:

(1) 由式(2)求简支梁的柔度矩阵,如有刚域或部分刚接情形,需作修正;

(2) 求逆获得简支梁的刚度矩阵;

(3) 按式(7)求梁的刚度矩阵。

前两步针对变截面梁(无显式),对于有显式表达的等截面梁,可直接写出${\mathbf K}_{\rm s} $,再变换即可。

梁的刚度矩阵与等效结点载荷在有限元软件中应用非常广泛,本文从有限元软件实现的角度给出了一个通用梁单元的框架,便于工程软件开发人员参考借鉴。以上仅是我们对有限元梁单元的一点思考,有不当与谬误之处,请同行批评指正。