图1所示的一端固支均匀杆的二阶控制微分方程为

\[EA\dfrac{d ^2u}{d x^2} + f = 0 \ \ (1)\]式中, \(A\)为截面面积, \(E\)为弹性模量, \(u(x)\)为轴向位移函数, \(f(x)\)为任意次分布载荷. 边界条件为\(u(0) = 0\), \(d u(L) / d x = 0\). 设\(f\)的形式为\[f(x) = \sum_{k = 0}^n q_k x^k , \ \ n = 0, 1, 2, \cdots \ \ (2)\]其中\(q_{k}\)为已知系数. 在第2和第3节讨论梁单元情况时, 所用横向分布载荷与式(2)相同.根据边界条件求解方程(1)可得位移精确解为\[u(x) = \dfrac{1}{EA}\sum_{k = 0}^n {\dfrac{q_k }{k + 1}\left( {L^{k + 1}x - \dfrac{x^{k + 2}}{k + 2}}\right)} \ \ (3)\]

图1

	Figure Option ViewDownloadNew Window
	图1 受分布载荷作用的一端固支均匀杆

二次杆单元位移函数为

\[ u(x) = \phi _1 u_1 + \phi _2 u_2 + \phi _3 u_3\ \ (4)\]

式中, \(u_1 \), \(u_2 \)和\(u_3 \)分别为单元左端、右端和中点的位移, \(\phi_i\) (\(i =1, 2, 3)\)为利用拉格朗日插值多项式得到的节点位移形函数. 二次杆单元的刚度矩阵\( {k}\)和载荷向量\({f}\)分别为

\[{\pmb k} = \dfrac{EA}{3L}\left[ \begin{array}{ccc} 7 & 1 & { - 8} \\ 1 & 7 & { - 8} \\ { - 8} & { - 8} & {16} \end{array}\!\! \right] \ \ (5a)\]

\[{\pmb f} = \left[ \begin{array}{c} {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}(1 - k)}{(k + 1)(k + 2)(k + 3)}}} \\ {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}(k + 1)}{(k + 2)(k + 3)}} }\ \ {\sum_{k = 0}^n {\dfrac{4q_k L^{k + 1}}{(k + 2)(k + 3)}} }\end{array} \!\! \right] \ \ (5b)\]

式(5)对应的节点位移列向量为\({\pmb u}^{\rm T} = [ {u_1 } \ \ {u_2 } \ \ {u_3 } ]\). 由此得到的二次杆单元节点位移是

\[\left.\begin{array}{l} u_2 = \dfrac{1}{EA}\sum_{k = 0}^n {\dfrac{q_k L^{k + 2}}{k + 2}} \ \ u_3 = \dfrac{1}{EA}\sum_{k = 0}^n {\dfrac{q_k L^{k + 2}(2k + 9)}{4(k +2)(k + 3)}} \end{array} \right \} \ \ (6)\]

由式(3)可以得到自由端以及中点位移的精确解分别为

\[\left.\begin{array}{l} u_L = \dfrac{1}{EA}\sum_{k = 0}^n {\dfrac{q_k L^{k + 2}}{k + 2}} \\ u_{\tfrac{L}{2}} = \dfrac{1}{EA}\sum_{k\mbox{ = }0}^n {\dfrac{q_k L^{k+ 2}\left[ {k + 2 - 2^{ - (k + 1)}} \right]}{2(k + 1)(k + 2)}} \end{array} \right \} \ \ (7)\]

把有限元解式(6)与精确解式(7)对比, 可以得到如下结论：

(1)二次杆单元的外部节点或端点位移总是精确的.

(2)仅当\(n=0\)即载荷均布时, 内部节点位移才是精确的.实际上, 均布载荷作用下的精确位移函数是二次的, 因此二次杆单元给出的位移函数场是精确的.

对于结论(1), 也可以通过静力凝聚方法来证明, 下面为具体过程.

首先写出如下分块形式的单元平衡方程\[\left[\!\! \begin{array}{cc} {\pmb k}_{\rm ee} & {\pmb k}_{\rm ei} \\ {\pmb k}_{\rm ie} & {\pmb k}_{\rm ii}\end{array} \right]\left[\!\! \begin{array}{c} {\pmb u}_{\rm e} \\ {\pmb u}_{\rm i}\end{array} \!\! \right] = \left[ \begin{array}{c}{\pmb f}_{\rm e} \\ {\pmb f}_{\rm i} \end{array} \!\! \right] \ \ (8)\]

其中, \({\pmb u}_{\rm e}^{\rm T} = \left[ {u_1 } \ \ {u_2 } \right]\)为外部节点参数列向量, \({u}_{\rm i}^{\rm T} = \left[ {u_3 } \right]\)为内部节点参数列向量; 分块刚度矩阵为\[\left.\begin{array}{cc}{\pmb k}_{\rm ee} = \dfrac{EA}{3L}\left[ \begin{array}{cc} 7 & 1 \\ 1 & 7 \end{array} \right] , \ \ {\pmb k}_{\rm ei} =\dfrac{EA}{3L}\left[ \begin{array}{c} { - 8} \\ { - 8}\end{array} \right] \\ {\pmb k}_{\rm ie} = {\pmb k}_{\rm ei}^{\rm T} , \ \ {\pmb k}_{\rm ii} = \dfrac{EA}{3L}\left[ {16} \right]\end{array} \!\!\right\} \ \ (9)\]

而\({\pmb f}_{\rm e} \)和\({\pmb f}_{\rm i} \)分别是与\({\pmb u}_{\rm e}\)和\({\pmb u}_{\rm i}\)对应的载荷矢量

\[\left.\begin{array}{c} {\pmb f}_{\rm e}^{\rm T} =\Bigg [ {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}(1 - k)}{(k + 1)(k + 2)(k + 3)}} } \\ \qquad {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}(k + 1)}{(k + 2)(k + 3)}} } \Bigg ] \\ {\pmb f}_{\rm i}^{\rm T} = \left[ {\sum_{k = 0}^n {\dfrac{4q_k L^{k + 1}}{(k + 2)(k + 3)}} } \right] \end{array} \right\} \ \ (10)\]

利用凝聚方法从方程(8)可以得到\[\hat {\pmb k}{\pmb u}_{\rm e} = \hat{\pmb f} \ \ (11)\]

其中

$$\left. \begin{array}{l} \hat {\pmb k} = {\pmb k }_{{\rm ee}} - {\pmb k}_{\rm ei} {\pmb k }_{\rm ii}^{ - 1} {\pmb k}_{\rm ie} = \dfrac{EA}{L}\left[\!\! \begin{array}{cc}1 & { - 1} \\ { - 1} & 1\end{array} \!\!\right] \\ \hat {\pmb f} = {\pmb f}_{\rm e} - {\pmb k}_{\rm ei}{\pmb k}_{\rm ii}^{ - 1} {\pmb f}_{\rm i} = \left[ \begin{array}{c} {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}}{(1 + k)(2 + k)}} } \\ {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}}{(2 + k)}} }\end{array} \!\! \right] \end{array}\right\} \ \ (12)$$

式(12)给出的就是两节点线性或等应变杆单元的刚度矩阵和载荷列向量. 内部节点参数向量为

\[ {\pmb u}_{\rm i} = {\pmb k}_{\rm ii}^{ - 1} \left( {{\pmb f}_{\rm i} - {\pmb k}_{\rm ie} {\pmb u}_{\rm e} } \right) = \dfrac{1}{EA}\sum_{k = 0}^n \dfrac{q_k L^{k + 2}(2k + 9)}{4(k + 2)(k + 3)}\ \ (13)\]

因此, 凝聚内部节点参数后, 二次杆单元刚度矩阵和载荷列向量分别退化成为线性杆单元的刚度矩阵和载荷列向量, 而已有文献或利用功的互等定理已经证明线性杆单元在任意次多项式分布载荷作用下其节点位移是精确的, 由此再次证明二次杆单元端点位移的精确性. 上述这种 静力凝聚方法适用于下面要讨论的梁.

将上述结论推广如下：在任意次分布载荷作用下, 任意次杆单元的端点位移总是精确的, 对于二次及以上杆单元情况都可以利用静力凝聚方法来证明; 当单元次数与精确位移函数的次数相同时, 单元内部节点位移才是精确的.该结论适用于轴单元, 对于欧拉梁和铁木辛柯梁单元, 可以得到类似的结论, 参见 第2节和第3节.

图2所示一端固支均匀欧拉梁的域内非齐次四阶欧拉平衡方程为

\[EI\dfrac{d ^4w}{d x^4} = f \ \ (14)\]

式中, \(I\)为截面惯性矩, \(w(x)\)为挠度函数. 边界条件为\(w(0) = 0\), \({d w(0)} /{d x} = 0\), \({d ^2w(L)} / {d x^2} = 0\)和\({d ^3w(L)} / {d x^3} = 0\). 根据边界条件求解方程(14), 可得挠度及其导数的精确解

\[\left.\begin{array}{l} w(x) = \dfrac{1}{EI}\sum_{k = 0}^n {q_k } \left[ {\dfrac{x^{k + 4}}{(k+ 1)(k + 2)(k + 3)(k + 4)} - }\right. \left.{ \dfrac{L^{k + 1}x^3}{6(k + 1)} + \dfrac{L^{k +2}x^2}{2(k + 2)}} \right] \\ \dfrac{d w(x)}{d x} = \dfrac{1}{EI}\sum_{k = 0}^n {q_k }\left[ {\dfrac{x^{k + 3}}{(k + 1)(k + 2)(k + 3)} - }\right. \\ \qquad \left.{ \dfrac{L^{k + 1}x^2}{2(k +1)} + \dfrac{L^{k + 2}x}{(k + 2)}} \right] \end{array} \!\!\right\} \ \ (15)\]

图2

	Figure Option ViewDownloadNew Window
	图2 分布载荷作用下的均匀欧拉梁

欧拉梁单元容许位移函数要保证挠度及其一阶导数的连续性. 利用三点一阶艾米(Hermite)插值多项式得到的五次欧拉梁单元容许位移函数为

\[ w(x) = \phi _1 (x)w_1 + \phi _2 (x)\theta _1 + \phi _3 (x)w_2 + \phi _4 (x)\theta _2 +\phi _5 (x)w_3 + \phi _6 (x)\theta _3 \ \ (16)\]

式中, \(w_1 \), \(w_2 \)和\(w_3 \)分别为左端、右端和内部节点挠度, 参见图3, \(\theta _i = {d w_i } / {d x}\) (\(i =1, 2, 3\))为挠度斜率或截面转角, 规定从\(x\)轴转向\(z\)轴为正.

图3

	Figure Option ViewDownloadNew Window
	图3 五次欧拉梁单元节点参数

根据结构刚度矩阵和载荷列向量得到的节点位移为

\[\left.\!\!\begin{array}{l} w_2 = \dfrac{1}{EI}\sum_{k = 0}^n {q_k L^{k + 4}\dfrac{2k + 9}{6(k +3)(k + 4)}} \\ \theta _2 = \dfrac{1}{EI}\sum_{k = 0}^n {q_k L^{k + 3}\dfrac{1}{2(k +3)}} \\ w_3 = \dfrac{1}{EI}\sum_{k = 0}^n {q_k L^{k + 4}\dfrac{10k^2 + 98k +255}{96(k + 3)(k + 4)(k + 5)}} \\ \theta _3 = \dfrac{1}{EI}\sum_{k = 0}^n {q_k L^{k + 3}\dfrac{6k^3 +94k^2 + 495k + 840}{16(k + 3)(k + 4)(k + 5)(k + 6)}} \end{array} \right \} \ \ (17)\]

对比式(15)和式(17)可以看出, 五次欧拉梁单元在任意次分布载荷作用下的端点位移都是精确的. 当\(n \leqslant\)1时, 内部节点位移才是精确的. 通过静力凝聚消除内部节点参数可以将五次欧拉梁单元退化为三次欧拉梁单元.

承受横向分布载荷作用的均匀铁木辛柯梁如图4所示, 多项式形式的分布弯矩为

\[m(x) = \sum_{j = 0}^n {m_j x^j} \quad (n = 0, 1, 2, \cdots) \ \ (18)\]

图4

	Figure Option ViewDownloadNew Window
	图4 分布载荷作用下的均匀铁木辛柯梁

以挠度\(w\)和剖面转角\(\psi \)作为自变函数, 其域内欧拉平衡方程为

\[ \dfrac{d }{d x}\left[ {\kappa GA\left( {\dfrac{d w}{d x} - \psi } \right)} \right]+ f = 0 \ \ (19a)\]

\[\dfrac{d }{d x}\left( {EI\dfrac{d \psi }{d x}} \right) + \kappa GA\left({\dfrac{d w}{d x} - \psi } \right) + m = 0\ \ (19b)\]

式中, \(\kappa \)为截面剪切修正系数. 边界条件为\(w(0) = 0\), \(\psi (0) = 0\), \({d \psi (L)} / {d x} = 0 \)和\(\left. {\left( {{d w} / {d x} - \psi } \right)} \right|_{x = L} = 0\). 方程(19)的精确解为

\[ w(x) = \dfrac{1}{EI}\sum_{k = 0}^n q_k \left[ {\dfrac{x^{k + 4}}{(k + 1)(k + 2)(k + 3)(k +4)} - }\right. \]

\[\qquad \left. { \dfrac{L^{k + 1}x^3}{6(k + 1)} + \dfrac{L^{k + 2}x^2}{2(k + 2)}} \right] -\]

\[\qquad \dfrac{1}{\kappa GA}\sum_{k = 0}^n {q_k \left[ {\dfrac{x^{k + 2}}{(k + 1)(k + 2)} -\dfrac{L^{k + 1}x}{k + 1}} \right]} -\]

\[\qquad \dfrac{1}{EI}\sum_{j = 0}^n {m_j \left[ {\dfrac{x^{j + 3}}{(j + 1)(j + 2)(j + 3)} -\dfrac{L^{j + 1}x^2}{2(j + 1)}} \right]} \ \ (20a)\]

\[ \psi (x) = \dfrac{1}{EI}\sum_{k = 0}^n q_k \left[ {\dfrac{x^{k + 3}}{(k + 1)(k + 2)(k +3)} - }\right.\]

\[\qquad \left.{ \dfrac{L^{k + 1}x^2}{2(k + 1)} + \dfrac{L^{k + 2}x}{k + 2}} \right] -\]

\[\qquad \dfrac{1}{EI}\sum_{j = 0}^n {m_j \left[ {\dfrac{x^{j + 2}}{(j + 1)(j + 2)} -\dfrac{L^{j + 1}x}{j + 1}} \right]} \ \ (20b)\]

从方程(19)和式(20b)可以看出, 剖面转角与剪切修正系数\(\kappa \)无关.

图5给出了三次铁木辛柯梁单元节点参数配置情况. 有独立插值和非独立插值两种方法来确定\(w\)和\(\psi \).

图5

	Figure Option ViewDownloadNew Window
	图5 三次铁木辛柯梁单元

首先来看独立插值方法. 用两点一阶艾米插值得到的\(w\)为

\[ w(x) = N_1 (x)w_1 + N_2 (x)\dfrac{d w_1 }{d x} + N_3 (x)w_2 + N_4 (x)\dfrac{d w_2 }{d x} \ \ (21)\]

用拉格朗日插值方法得到的\(\psi \)为

\[\psi (x) = N_5 (x)\psi _1 + N_6 (x)\psi _2 + N_7 (x)\psi _3 \ \ (22)\]

其中, 形函数\(N_1 , N_2 , N_3 , N_4 \)为三次艾米插值函数, \(N_5 , N_6 , N_7 \)为二次拉格朗日插值函数.

在式(21)和式(22)中, 端点的\(w_1 \), \(w_2 \), \(\psi _1 \)和\(\psi _2 \)与保证梁单元\(C_0\)连续性要求有关, 故称其为单元的外部节点参数. 配置在端点的\({d w_1 } /{d x}\)和\({d w_2 }/{d x}\)以及\(\psi _3\)与单元间位移协调无关, 称其为内部节点参数.用静力凝聚方法消除内部节点参数后, 得到如下刚度矩阵和载荷列向量^[6]

\[{\pmb k} = \dfrac{EI}{(1 + S)L^3}\left[\!\!\begin{array}{cccc} {12} & {6L} & { - 12} & {6L} \\ & {(4 + S)L^2} & { - 6L} & {(2 - S)L^2} \\ \hbox{对} & & {12} & { - 6L} \\ & \hbox{称} & & {(4 + S)L^2}\end{array}\!\!\right] \ \ (23a)\]

\[{\pmb f} = \dfrac{1}{S + 1}\left[\!\! \begin{array}{c} {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}(12S + 7Sk + Sk^2 + 6k + 12)}{(k + 1)(k + 2)(k + 3)(k + 4)}} - \sum_{j = 0}^n {\dfrac{6m_j L^j}{(j + 2)(j + 3)}} } \\ {\sum_{k = 0}^n {\dfrac{q_k L^{k + 2}(4S + Sk + 4)}{2(k + 2)(k + 3)(k+ 4)}} + \sum_{j = 0}^n {\dfrac{m_j L^{j + 1}(3S + Sj - 2j)}{(j + 1)(j + 2)(j + 3)}} } \\ {\sum_{k = 0}^n {\dfrac{q_k L^{k + 1}(12S + 7Sk + Sk^2 + k^2 + 8k +12)}{(k + 2)(k + 3)(k + 4)}} + \sum_{j = 0}^n {\dfrac{6m_j L^j}{(j + 2)(j + 3)}} }\\ {\sum_{k = 0}^n {\dfrac{ - q_k L^{k + 2}(4S + Sk + 2k + 4)}{2(k +2)(k + 3)(k + 4)}} + \sum_{j = 0}^n {\dfrac{m_j L^{j + 1}(3S + Sj + j)}{(j + 2)(j + 3)}} }\end{array} \!\! \right] \ \ (23b)\]

式(23)对应的节点位移列向量为\({\pmb w }^{\rm T} = [ {w_1 } \ \ {\psi _1 }\) \ \ \({w_2 } \ \ {\psi _2 } ]\). 缩聚后的单元刚度矩阵能够适用于各种 \(S = \dfrac{12EI}{\kappa GAL^2}\)值. 当\(S = 0\), \(m_{j} =0\)时, 式(23)将退化为三次欧拉梁单元的刚度矩阵和载荷列向量. 利用式(23)求得的外部节点位移为

\[\left. \!\! \begin{array}{l} w_2 = \dfrac{1}{EI}\sum_{k = 0}^n {\dfrac{q_k L^{k + 4}(12S + 7Sk +Sk^2 + 4k^2 + 26k + 36)}{12(k + 2)(k + 3)(k + 4)}} + \dfrac{1}{EI}\sum_{j = 0}^n {\dfrac{m_j L^{j + 3}(j +4)}{2(j + 2)(j + 3)}} \\ \psi _2 = \dfrac{1}{EI}\sum_{k = 0}^n {\dfrac{q_k L^{k + 3}}{2(k +3)}} + \dfrac{1}{EI}\sum_{j = 0}^n {\dfrac{m_j L^{j + 2}}{(j + 2)}} \end{array} \right \} \ \ (24a)\]

内部节点位移为

\[\left.\begin{array}{l} \dfrac{d w_1 }{d x} = \dfrac{1}{EI}\sum_{k = 0}^n{\dfrac{q_i L^{k + 3}S(k^2 + 5k + 24)}{12(k + 2)(k + 3)(k + 4)}} \\ \dfrac{d w_2 }{d x} = \dfrac{1}{EI}\sum_{k = 0}^n{\dfrac{q_k L^{k + 3}(Sk^2 - Sk + 6k^2\mbox{ + }36k\mbox{ + }48)}{12(k + 2)(k + 3)(k + 4)}} +\dfrac{1}{EI}\sum_{j = 0}^n {\dfrac{m_j L^{j + 2}}{(j + 2)}} \\ \psi _3 = \dfrac{1}{EI}\sum_{k = 0}^n {\dfrac{q_k L^{k + 3}(3k\mbox{ +1}4)}{8(k + 3)(k + 4)}} \mbox{ + }\dfrac{1}{EI}\sum_{j = 0}^n {\dfrac{m_j L^{j + 2}(2j + 9)}{4(j + 2)(j +3)}} \end{array} \right \} \ \ (24{\rm b})\]

把式(24)和式(20)的精确解进行比较可以发现, 对于分布载荷作用下的均匀铁木辛柯梁, 用上述独立插值得到的三次铁木辛柯梁单元的外部节点挠度和转角都是精确的.内部节点位移参数\({d w_1 } /{d x}\)和\({d w_2 } /{d x}\)在\(n \leqslant \)1时是精确的; 但\(\psi _3\)总是不精确的, 这是因为该单元转角函数是二次的, 而\(n=\)0和\(n=\)1时转角精确解的次数分别是3和4, 参见式(20b).

用这种独立插值方法还可以构造5次以上的铁木辛柯梁单元, 其外部位移参数总是精确的, 内部位移参数是否精确的分析方法同前.

利用非独立插值得到的三次铁木辛柯梁单元的方法如下^[3]. 令\(w\)和\(\psi \)分别为完备的三次和二次代数多项式, 把它们代入齐次欧拉方程, 可以得到\(w\)和\(\psi \)的系数关系, 最后有

\[\left.\begin{array}{l} w(x) = ax^3 + bx^2 + cx + d \\ \psi (x) = 3ax^2 + 2bx + c + \dfrac{6EI}{\kappa GA}a \end{array} \right \} \ \ (25)\]

其中的\(a\), \(b\), \(c\) 和\(d\)可以用\(w_1 \), \(w_2 \), \(\psi _1 \)和\(\psi _2 \)来替换.利用这种梁单元得到的刚度矩阵和载荷向量与式(23)相同.