基于荷载等效理论的柱体结构剪切变形解析

赵立财

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

为了提高具有铁木辛科梁性质的柱体结构剪切变形的计算效率和精度，基于荷载等效分布理论、节点解析法以及具有正交内插特性的铁木辛柯梁理论，提出了一种高效求解梁柱节点剪切变形位移的等效分布荷载⁃有限元算法，并通过研究2种工况算例验证了算法的可靠性。结果表明：该算法利用非常少的单元（1个或2个单元）便可以获得位移、旋转角、剪切力和弯矩的高精度近似解；其他基于位移的方法，如简化积分法，大约需要40个单元才能获得类似的结果。所提方法可以在不同的连接条件下直接获得铁木辛柯梁柱的稳定函数和弯曲荷载，有效提高了柱体结构剪切变形有限元分析解的精确性。

关键词

有限元法; 精确节点解; 等效分布荷载; 铁木辛柯梁柱; 弯曲荷载

引言

在基于位移的有限元法中，一些学者使用同质微分方程的解作为有限元位移场的形函数，获得了节点精确位移解。文献［

1‑2］提出了一种微分方程的非自伴随算子求解方法，在非对角阵中使用加权函数的条件下，得到了方程伴随算子的齐次微分方程的解。文献［3‑4］研究了基于节点非线性板单元的精度解。在结构领域中，文献［5‑6］将其应用于铁木辛柯梁模型，使用齐次微分方程组解的形函数证明了单元不受阻塞现象的影响。Balobanov等^{［参考文献 7

百度学术}7］采用精确节点解的有限元法，分析铁木辛柯梁模型的自由波传播与频率响应规律。Rajasekaran等^{［参考文献 8

百度学术}8］将精确节点解的有限元法和等效分布荷载概念应用于伯努利⁃欧拉梁模型中，以模拟具有单个或多个裂纹梁的双向功能梯度刚度变化。文献［9⁃11］将精确节点解的有限元法和等效分布荷载概念都应用到铁木辛柯梁模型中。文献［12‑13］将等效分布荷载概念应用于伯努利⁃欧拉模型中柱结构的屈曲分析，揭示了柱结构在轴向拉拔荷载作用下的破坏行为。

笔者认为铁木辛柯柱是伯努利⁃欧拉柱的一种特殊情况，基于等效分布荷载⁃有限元法提出了铁木辛柯柱的精确节点解的有限元理论，该方法可以直接获得2种类型支柱的弯曲荷载和稳定性函数。本研究基于等效分布荷载与正交内插的有限元方法，在每个有限元内对力、位移及旋转角进行近似求解，并通过2个算例说明了该方法的优点。这些有限元法在铁木辛柯模型中的应用非常有效，主要体现在具有剪切变形的梁柱中^［

14⁃16］。

1 精确节点解的有限元计算理论

铁木辛柯柱的微分方程为 $L (U) = F$ ^［

17］，则有

\{\begin{array}{l} \begin{array}{l} L (U) = [\begin{matrix} - (k A G (w^{'} {- ψ))}^{'} + (P w^{'})^{'} \\ - (E I ψ^{'})^{'} - k A G (w^{'} - ψ) \end{matrix}] \\ U = {[\begin{matrix} w & ψ \end{matrix}]}^{T} \end{array} \\ F = {[\begin{matrix} f & 0 \end{matrix}]}^{T} \end{array}

（1）

其中： $w$ 为位移； $ψ$ 为旋转角； $P$ 为荷载； $k$ ， $A$ ， $G$ 分别为剪切力、截面面积和剪切刚度模量的校正因子。

假设 $H$ 为刚性 $E I$ ， $K$ 为 $k A G$ 的积，式（1）可表示为

\begin{array}{l} - (K (w^{'} {- ψ))}^{'} + (P w^{'})^{'} = \\ f - (H ψ^{'})^{'} - K (w^{'} - ψ) = 0 \end{array}

（2）

通过对原系统 $V = {[\begin{matrix} v & θ \end{matrix}]}^{T}$ 进行加权，并在一般子区间 $[α, β]$ 上进行相应的分段积分，得到了类似的变分公式为

b (U, V) = l (V) + \sum_{i = 1}^{4} l_{i} (V) {\hat{l}}_{i} (U)

（3）

双线性和线性形式^［

18］分别为

\{\begin{array}{l} b (U, V) = a (U, V) - \int_{α}^{β} P w^{'} v^{'} d x \\ a (U, V) = \int_{α}^{β} \{H ψ^{'} θ^{'} + K (w^{'} - ψ) (v^{'} - θ)\} d x \\ l (V) = \int_{α}^{β} v f d x \end{array}

（4）

\{\begin{array}{l} l_{1} (V) = v (α) \\ l_{2} (V) = θ (α) \\ l_{3} (V) = v (β) \\ l_{4} (V) = θ (β) \end{array}

（5）

\{\begin{array}{l} {\hat{l}}_{1} (U) = {[- K (w^{'} - ψ)]}_{x = α} + P w^{'} (α) \\ {\hat{l}}_{2} (U) = {[- H ψ^{'}]}_{x = α} \\ {\hat{l}}_{3} (U) = {[- K (w^{'} - ψ)]}_{x = β} - P w^{'} (β) \\ {\hat{l}}_{4} (U) = {[H ψ^{'}]}_{x = β} \end{array}

（6）

例如，在 $x = a$ 端和 $x = b$ 端自由嵌套的零件，在自由端受到恒定轴向力 $P$ 和荷载 $F_{H}$ 作用，此时梁柱弯矩 $M$ 的变分问题可表示为

(w, ψ) \in H_{0}^{1} (Ω) H_{0}^{1} (Ω)

（7）

即

\{\begin{array}{l} \int_{a}^{b} \{H ψ^{'} θ^{'} + K (w^{'} - ψ) (v^{'} - θ) - P w^{'} v^{'}\} d x = \\ \int_{a}^{b} f v d x + F_{H} v (b) + M θ (b) \\ \forall (v, θ) \in H_{0}^{1} (Ω) H_{0}^{1} (Ω) \end{array}

（8）

其中： $H_{0}^{1} (Ω)$ 为函数 $g$ 的索伯列夫空间， $g^{'} \in L^{2} (Ω)$ ， $g (a) = 0$ 。

该问题的唯一解是使能量泛函静止的解，即

I (w, ψ) = \int_{a}^{b} \{\frac{H}{2} ψ^{2} + \frac{K}{2} {(w^{'} - ψ)}^{2} - \frac{P}{2} {w^{'}}^{2}\} d x - (\int_{a}^{b} f w d x + F_{H} w (b) + M ψ (b))

（9）

双线性形式的对称性允许在任何子区间或总区间上表示，即

b (U, V) = \int_{α}^{β} V^{T} L (U) d x + \sum_{i = 1}^{4} l_{i} (V) {\hat{l}}_{i} (U) = \int_{α}^{β} U^{T} L (V) d x + \sum_{i = 1}^{4} l_{i} (U) {\hat{l}}_{i} (V)

（10）

式（10）取加权函数，使其为铁木辛柯柱齐次方程的解，得到

b (U, V) = \sum_{i = 1}^{4} l_{i} (U) {\hat{l}}_{i} (V)

（11）

式（11）不依赖于 $U$ 在子区间内的值的双线性形式，单元的变分方程为

b (U, V) = l (V) + \sum_{i = 1}^{4} l_{i} (V) {\hat{l}}_{i} (U)

（12）

根据式（11）与式（12）可以得出

\sum_{i = 1}^{4} l_{i} (U) {\hat{l}}_{i} (V) = l (V) + \sum_{i = 1}^{4} l_{i} (V) {\hat{l}}_{i} (U)

（13）

因此，有限元平衡方程的构造，即在局部方程和全局方程2种情况下，均得到了方程的解 $N_{i} = {[\begin{matrix} N_{1 i} & N_{2 i} \end{matrix}]}^{T}$ ， $i = 1 ~ 4$ 。这些函数构成了插值的拉格朗日基函数，即 $l_{i} (N_{j}) = δ_{i j}$ ， $i, j = 1 ~ 4$ 。在这种情况下，微分方程组为

\{\begin{array}{l} H_{1} ψ^{‴} + P ψ^{'} = f \\ w^{'} = ψ - m ψ^{″} \end{array}

（14）

其中： $H_{1} = H (1 - P / K)$ ； $m = H / K$ 。

$[α, β]$ 中的平衡方程为

K^{e} u^{e} = f^{e} + q^{e}

（15）

其中： $u^{e}$ 为节点位移； $q^{e}$ 为等效节点荷载； $f^{e}$ 为为平衡节点。

$u^{e}$ ， $q^{e}$ 和 $f^{e}$ 分别表示为

u^{e} = {[\begin{matrix} u_{1}^{e} & u_{2}^{e} & u_{3}^{e} & u_{4}^{e} \end{matrix}]}^{T} = {[\begin{matrix} w (α) & ψ (α) & w (β) & ψ (β) \end{matrix}]}^{T}

（16）

f^{e} = {[\begin{matrix} f_{1}^{e} & f_{2}^{e} & f_{3}^{e} & f_{4}^{e} \end{matrix}]}^{T}

（17）

f_{i}^{e} = \int_{α}^{β} f N_{1 i} d x (i = 1 ~ 4)

（18）

\{\begin{array}{l} q^{e} = {[\begin{matrix} q_{1}^{e} & q_{2}^{e} & q_{3}^{e} & q_{4}^{e} \end{matrix}]}^{T} \\ q_{i}^{e} = {\hat{l}}_{i} (U) (i = 1 ~ 4) \end{array}

（19）

局部刚度矩阵 $K^{e} = [k_{i j}^{e}]$ 的单元为 $k_{i j}^{e} = b (N_{i}, N_{j})$ ，可通过积分形式计算，即

k_{i j}^{e} = b (N_{i}, N_{j}) = \int_{α}^{β} H N_{2 i}^{'} N_{2 j}^{'} d x + \int_{α}^{β} K (N_{1 i}^{'} - N_{2 i}) (N_{1 j}^{'} - N_{2 j}) d x - \int_{α}^{β} P N_{1 i}^{'} N_{1 j}^{'} d x

（20）

也可以通过式（12）和式（20）推导出形函数^［

8，11，18］，即

k_{i j}^{e} = b (N_{i}, N_{j}) = \sum_{n = 1}^{4} l_{n} (N_{i}) {\hat{l}}_{n} (N_{j}) = \sum_{n = 1}^{4} δ_{n i} {\hat{l}}_{n} (N_{j}) = {\hat{l}}_{i} (N_{j})

（21）

铁木辛柯柱 $K_{P T}^{e}$ 和伯努利⁃欧拉柱 $K_{B E}^{e}$ 的精确单元刚度矩阵分别为

K_{P T}^{e} = \frac{H_{1}}{\bar{Δ}} \times [\begin{matrix} - r_{1}^{3} s_{1} & ρ r_{1}^{2} (c_{1} - 1) & r_{1}^{3} s_{1} & ρ r_{1}^{2} (c_{1} - 1) \\ ρ r_{1} (h r_{1} c_{1} - ρ s_{1}) & ρ r_{1}^{2} (1 - c_{1}) & ρ r_{1} (ρ s_{1} - h r_{1}) \\ - r_{1}^{3} s_{1} & ρ r_{1}^{2} (1 - c_{1}) \\ S I M & ρ r_{1} (h r_{1} c_{1} - ρ s_{1}) \end{matrix}]

（22）

其中： $H_{1} = H (1 - P / K)$ ； $m = H / K$ ； $r_{1} = \sqrt[]{P / H_{1}}$ ； $h = β - α$ ； $c_{1} = c o s (r_{1} h)$ ； $s_{1} = s i n (r_{1} h)$ ； $ρ = 1 +$ $r_{1}^{2} m$ ； $\bar{Δ} = 2 ρ (c_{1} - 1) + h r_{1} s_{1}$ ；SIM为局部坐标系中的刚度矩阵。

K_{B E}^{e} = \frac{H}{Δ} [\begin{matrix} - r^{3} s & r^{2} (c - 1) & r^{3} s & r^{2} (c - 1) \\ r (h r c - s) & r^{2} (1 - c) & r (s - h r) \\ - r^{3} s & r^{2} (1 - c) \\ S I M & r (h r c - s) \end{matrix}]

（m=0）

（23）

其中： $r = \sqrt[]{P / H}$ ； $h = β - α$ ； $c = c o s (r h)$ ； $s = s i n (r h)$ ； $Δ = h r s + 2 (c - 1)$ 。

式（22）所示单元内只考虑剪切变形和扭转变形，忽略了轴向变形。因此， $K_{P T}^{e}$ 和 $K_{B E}^{e}$ 的精确单元刚度矩阵分别为4×4矩阵，其中：第1、第3列对应扭转自由度；第2、第4行对应切向自由度。

2 弯曲荷载和稳定性特性

当考虑到剪切变形时，弯曲荷载的值会变小。该荷载是 $P$ 的最小正值，其抵消了还原矩阵 $K_{P T}^{e}$ 的行列式，即相应特征值问题的基本边界条件。

对于长度为L的横梁柱，两端（节点1和节点2）铰接，单元刚度矩阵 $K_{P T}^{e}$ 被简化为式（22）的第2、第4行和第2、第4列中的矩阵（4个单元），取消相应的行列式并用L替代h，令剩下的行列式等于0，可求得稳定临界平衡方程为

(\frac{H_{1}}{\bar{Δ}}) |\begin{matrix} ρ r_{1} (L r_{1} c_{1} - ρ s_{1}) & ρ r_{1} (ρ s_{1} - L r_{1}) \\ ρ r_{1} (ρ s_{1} - L r_{1}) & ρ r_{1} (L r_{1} c_{1} - ρ s_{1}) \end{matrix}| = 0

（24）

将式（24）作相应运算后得到

(H_{1}^{2} / \bar{Δ}) L ρ r_{1}^{3} s_{1} = 0

（25）

假设 $s_{1} = s i n (r_{1} L) = 0$ ，取最小的正根 $r_{1} L = π$ ，得到弯曲荷载为

P = r_{1}^{2} H / (1 + m r_{1}^{2})

（26）

由 $r_{1} = \sqrt[]{P / H_{1}}$ 推导出的两端铰接临界弯曲荷载为

P_{c r T} = (\frac{π^{2} H}{L^{2}}) / (1 + \frac{π^{2} H}{L^{2} K})

（27）

对于长度为L的横梁柱，在一端嵌套（节点1）、另一端自由（节点2）的工况下， $K_{P T}^{e}$ 简化为式（22）中的第3、第4行和第3、第4列中的矩阵（4个单元），取消相应的行列式并用L替代h后，令剩下的行列式等于0，可求得稳定临界平衡方程为

(\frac{H_{1}}{\bar{Δ}}) |\begin{matrix} - r_{1}^{3} s_{1} & ρ r_{1}^{2} (1 - c_{1}) \\ ρ r_{1}^{2} (1 - c_{1}) & ρ r_{1} (L r_{1} c_{1} - ρ s_{1}) \end{matrix}| = 0

（28）

将式（28）作相应运算后得到

(H_{1}^{2} / \bar{Δ}) c_{1} ρ r_{1}^{4} = 0

（29）

假设 $c_{1} = c o s (r_{1} L) = 0$ ，取最小正根 $r_{1} = π / (2 L)$ ，代入式（26）中，则一端自由、另一端固定的临界弯曲荷载为

P_{c r T} = (\frac{π^{2} H}{4 L^{2}}) / (1 + \frac{π^{2} H}{4 L^{2} K})

（30）

对于两端嵌套的工况， $K_{P T}^{e}$ 为式（22）所示的4×4矩阵（10个单元），得到稳定临界平衡方程为

\frac{H_{1}}{\bar{Δ}} \times [\begin{matrix} - r_{1}^{3} s_{1} & ρ r_{1}^{2} (c_{1} - 1) & r_{1}^{3} s_{1} & ρ r_{1}^{2} (c_{1} - 1) \\ ρ r_{1} (h r_{1} c_{1} - ρ s_{1}) & ρ r_{1}^{2} (1 - c_{1}) & ρ r_{1} (ρ s_{1} - h r_{1}) \\ - r_{1}^{3} s_{1} & ρ r_{1}^{2} (1 - c_{1}) \\ S I M & ρ r_{1} (h r_{1} c_{1} - ρ s_{1}) \end{matrix}] = 0

（31）

类似地，求得两端嵌套的临界弯曲荷载为

P_{c r T} = (\frac{2 π^{2} H}{3 L^{2}}) / (1 + \frac{2 π^{2} H}{3 L^{2} K})

（32）

按照文献［

18］的方法，取消了矩阵

K_{P T}^{e}

的行列式，得到了精确的嵌入式加工柱和嵌入式铰接柱的弯曲荷载值。根据伯努利⁃欧拉柱的弯曲荷载，当柱的构型是嵌套⁃自由式、铰接⁃铰接式和嵌套⁃嵌套式时，铁木辛柯柱的临界弯曲荷载可设为

P_{c r T} = P_{c r E} / (1 + P_{c r E} / K)

（33）

式（33）也适用于两端具有相同刚度的弹性嵌套柱。

当柱的构型是嵌套⁃铰接式时，有近似表达式为

P_{c r T} \approx P_{c r E} / (1 + 1.1 P_{c r E} / K)

（34）

由式（22）可以直接得到稳定性函数。对于长度为L、两端铰接的横梁柱，其力矩与弯矩的关系为

[\begin{matrix} M_{1} \\ M_{1} \end{matrix}] = \frac{H_{1}}{\bar{Δ}} [\begin{matrix} ρ r_{1} (L r_{1} c_{1} - ρ s_{1}) & ρ r_{1} (ρ s_{1} - L r_{1}) \\ ρ r_{1} (ρ s_{1} - L r_{1}) & ρ r_{1} (L r_{1} c_{1} - ρ s_{1}) \end{matrix}] [\begin{matrix} ψ_{1} \\ ψ_{2} \end{matrix}]

（35）

在式（35）中代入 $x = r_{1} L (无量纲)$ ，得到

[\begin{matrix} M_{1} \\ M_{1} \end{matrix}] = \frac{H_{1}}{L} [\begin{matrix} \frac{x (x c o s x - ρ s i n x)}{\bar{Δ}} & \frac{x (ρ s i n x - x)}{\bar{Δ}} \\ \frac{x (ρ s i n x - x)}{\bar{Δ}} & \frac{x (x c o s x - ρ s i n x)}{\bar{Δ}} \end{matrix}] [\begin{matrix} ψ_{1} \\ ψ_{2} \end{matrix}]

（36）

假设 $ρ$ 和 $\bar{Δ}$ 为 $x$ 的函数，则有

\{\begin{array}{l} ρ = 1 + {(x / L)}^{2} m \\ \bar{Δ} = 2 [1 + {(x / L)}^{2} m] (c o s x - 1) + x s i n x \end{array}

（37）

同样，将式（26）代入 $x = r_{1} L (无量纲)$ ，得到

P = {(x / L)}^{2} H / (1 + {(x / L)}^{2} m)

（38）

荷载 $P$ 和 $x$ 的二元关系如图1所示。

图1 荷载 $P$ 和 $x$ 的二元关系

Fig.1 Binary relation of load $P$ and $x$

式（36）可以写成经典形式，即

[\begin{matrix} M_{1} \\ M_{1} \end{matrix}] = \frac{E I}{L} [\begin{matrix} \bar{s} & \bar{s} \bar{c} \\ \bar{s} \bar{c} & \bar{s} \end{matrix}] [\begin{matrix} ψ_{1} \\ ψ_{2} \end{matrix}]

（39）

其中： $\bar{s} = \frac{x (x c o s x - ρ s i n x)}{\bar{Δ}}$ ； $\bar{s} \bar{c} = \frac{x (ρ s i n x - x)}{\bar{Δ}}$ ； $\bar{c} = \frac{ρ s i n x - x}{x c o s x - ρ s i n x}$ 。

当L=5 m，m=0.26时，矩形截面的稳定性函数如图2所示。

图2 矩形截面的稳定性函数

Fig.2 Stability function of rectangular section

这些稳定性函数通常用于伯努利⁃欧拉梁柱和铁木辛柯柱以及寻找框架的临界荷载值^［

14］，本研究可以直接从刚度矩阵中构建2个混凝土支柱模型的稳定性函数。图2除了表示稳定性函数

\bar{s}

和

\bar{c}

外，还表示了稳定性函数

\bar{s} \bar{c}

和

{\bar{s}}^{2} (1 - {\bar{c}}^{2})

。后者是式（36）矩阵的行列式，根据式（27）和式（38），该表达式在

x = π

处对应的值被消去，则对应于

P_{c r T}

的值。

3 等效分布荷载

复数域分解形式为 $[a, b] = \cup_{i = 1}^{n - 1} [x_{i}, x_{i + 1}]$ ，考虑了等效分布荷载的概念。在 $[α, β]$ 中，验证了代表等效节点标量积相等的情况，即

\int_{α}^{β} f N_{1 i} d x = \int_{α}^{β} \bar{f} N_{1 i} d x (i = 1 ~ 4)

（40）

其中： $f$ ， $\bar{f}$ 为该分解的等效荷载。

当 $\bar{f}$ 由 $N_{1 i}$ （ $i = 1 ~ 4$ ）分量的线性组合给出时，解是齐次的， $\bar{f}$ 为 $f$ 在 $N_{1 i}$ （ $i = 1 ~ 4$ ）生成的函数空间上的正交投影。由式（2）推导出 $f$ 和 $\bar{f}$ 的微分方程，即

\{\begin{array}{l} {(H ψ^{'})}^{″} + {(P w^{'})}^{'} = f (x) \\ {(H \bar{ψ^{'}})}^{'} + {(P \bar{w^{'}})}^{'} = \bar{f} (x) \end{array}

（41）

考虑到荷载的等价性，有

\begin{array}{l} \int_{α}^{β} [({(H ψ^{'})}^{″} + {(P w^{'})}^{'}) - ({(H \bar{ψ^{'}})}^{'} + {(P \bar{w^{'}})}^{'})] v d x = \\ \int_{α}^{β} (f (x) - \bar{f} (x)) v d x = 0 \end{array}

（42）

\int_{α}^{β} [({(H ψ^{'})}^{″} + {(P w^{'})}^{'}) - ({(H \bar{ψ^{'}})}^{'} + {(P \bar{w^{'}})}^{'})] v d x = 0

（43）

其中： $v$ 为加权函数，属于齐次函数的解空间，即由 $N_{1 i}$ （ $i = 1 ~ 4$ ）生成的解空间。

根据式（43），得到部分积分表达式为

\begin{array}{l} [({(H ψ^{'})}^{'} + P w^{'}) - ({(H \bar{ψ^{'}})}^{'} + P \bar{w^{'}})] v |_{α}^{β} + \\ \int_{α}^{β} [(- {(H ψ^{'})}^{'} - P w^{'}) - ({(H \bar{ψ^{'}})}^{'} + P \bar{w^{'}})] v^{'} d x = 0 \end{array}

（44）

因为受到两个同等的作用，节点的位移和旋转角相等，式（44）中第一部分的和为0，因此有

\int_{α}^{β} [(- {(H ψ^{'})}^{'} - P w^{'}) - ({(H \bar{ψ^{'}})}^{'} + P \bar{w^{'}})] v^{'} d x = 0

（45）

将式（45）再次分部积分，得到

\begin{array}{l} [(H ψ^{'} + P w) - (H \bar{ψ^{'}} + P \bar{w})] v^{'} |_{α}^{β} - \\ \int_{α}^{β} [(H ψ^{'} + P w) - (H \bar{ψ^{'}} + P \bar{w})] v^{″} d x = 0 \end{array}

（46）

根据式（44），有

\int_{α}^{β} [(H ψ^{'} + P w) - (H \bar{ψ^{'}} + P \bar{w})] v^{″} d x = 0

（47）

联立式（45）~（47），得到

\{\begin{array}{l} V (x) = - H ψ^{″} - P w^{'} = K (w^{'} - ψ) - P w^{'} \\ \bar{V} (x) = - H \bar{ψ^{″}} - P \bar{w^{'}} = K (\bar{w^{'}} - \bar{ψ}) - P \bar{w^{'}} \\ M (x) = H ψ^{'} \\ \bar{M} (x) = H \bar{ψ^{'}} \end{array}

（48）

将式（48）进行简化，得到

\{\begin{array}{l} \int_{α}^{β} (V (x) - \bar{V} (x)) v^{'} d x = 0 \\ \int_{α}^{β} (V (x) - \bar{V} (x)) N_{1 i}^{'} d x = 0 \end{array}

（49）

\{\begin{array}{l} \int_{α}^{β} [(M (x) - \bar{M} (x)) + P (w - \bar{w}) +] v^{″} d x = 0 \\ \int_{α}^{β} [(M (x) - \bar{M} (x)) + P (w - \bar{w}) +] N_{1 i}^{″} d x = 0 \end{array}

（50）

梁柱每个截面 $x$ 的剪切力由表达式 $Q (x) = K (w^{'} - ψ)$ 给出。因此，式（48）中的函数 $V (x)$ 表示一种新的应力，称为修正剪切力（在截面右侧位移方向上的投影，包括 $P w^{'}$ 的贡献^［

18］）；函数

M (x)

为该截面的弯矩。

式（49）表明，因为插值相同，修正剪切差分函数 $V (x) - \bar{V} (x)$ 在区间两端为0，也与 $N_{1 i}$ ， $i = 1 ~ 4$ 生成的函数空间正交，这时产生了一个三维空间^［

2，8］。

由式（50）可得，函数 $m (x) - \bar{m} (x) = (M (x) - \bar{M} (x)) + P (w - \bar{w})$ 。因为插值相同，代表力矩差和位移差的荷载 $P$ ，在区间两端为零和零导数，也与 $N_{1 i}$ （ $i = 1 ~ 4$ ）生成的函数在二维空间正交。在每个单元中，对于 $V (x) - \bar{V} (x)$ 和 $m (x) - \bar{m} (x)$ ，正交性和插值特性是上述方法在每个有限元内对力、位移及旋转角进行近似的基础。

对于铁木辛柯柱，可以用类似的方法计算解 $\bar{U} = {[\begin{matrix} \bar{w} & \bar{ψ} \end{matrix}]}^{T}$ 和相对于 $\bar{f} (x)$ 的衍生应力 $\bar{Q} (x)$ 及 $\bar{M} (x)$ ，且不需要事先计算 $\bar{f} (x)$ ^［

8，11］。对于常数

H

和

K

的特殊情况，通过插值确定等效分布荷载

\bar{f} (x)

对应的近似解

\bar{U}

。由

H_{1} ψ^{‴} + P ψ^{'} = f

，

w^{'} = ψ - m ψ^{″}

推导出等价解

\bar{U} (z = x - α)

为

\bar{U} = [\begin{matrix} \bar{w} \\ \bar{ψ} \end{matrix}] = [\begin{matrix} 1 & z & z^{2} & z^{3} & {\bar{s}}_{1} & {\bar{c}}_{1} & z {\bar{s}}_{1} & z {\bar{c}}_{1} \\ 0 & 1 & 2 z & 3 z^{2} + 6 m & \frac{r_{1}}{ρ} \bar{c} {}_{1} & - \frac{r_{1}}{ρ} \bar{s} {}_{1} & \frac{\hat{ρ}}{ρ^{2}} {\bar{s}}_{1} + \frac{r_{1}}{ρ} z {\bar{c}}_{1} & \frac{\hat{ρ}}{ρ^{2}} {\bar{c}}_{1}_{1} - \frac{r_{1}}{ρ} z {\bar{s}}_{1}_{1} \end{matrix}] C^{- 1} \hat{u}

（51）

其中： $\hat{u} = [w (α), ψ (α), w^{'} (α^{+}), ψ^{'} (α^{+}), w (β),$ $ψ (β), w^{'} (β^{-}), ψ^{'} (β^{-})]$ 。

矩阵 $C$ 为

C = [\begin{matrix} 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 6 m & \frac{r_{1}}{ρ} & 0 & 0 & \frac{\hat{ρ}}{ρ^{2}} \\ 0 & 1 & 0 & 0 & r_{1} & 0 & 0 & 1 \\ 0 & 0 & 2 & 0 & 0 & - r_{1}^{2} / ρ^{2} & 2 r_{1} / ρ^{2} & 0 \\ 1 & h & h^{2} & h^{3} & s & c & h s_{1} & h c_{1} \\ 0 & 1 & 2 h & 3 h^{2} + 6 m & \frac{r_{1}}{ρ} c {}_{1} & - \frac{r_{1}}{ρ} s {}_{1} & \frac{\hat{ρ}}{ρ^{2}} s_{1} + \frac{r_{1}}{ρ} h c {}_{1} & \frac{\hat{ρ}}{ρ^{2}} c_{1} - \frac{r_{1}}{ρ} h s {}_{1} \\ 0 & 1 & 2 h & 3 h^{2} + 6 m & r_{1} c {}_{1} & - r_{1} s {}_{1} & s {}_{1}+ r_{1} h c {}_{1} & c {}_{1}- r_{1} h s {}_{1} \\ 0 & 0 & 2 & 6 h & - \frac{r_{1}^{2}}{ρ^{2}} s {}_{1} & - \frac{r_{1}^{2}}{ρ^{2}} c {}_{1} & \frac{2 r_{1}}{ρ^{2}} c {}_{1}- \frac{r_{1}^{2}}{ρ^{2}} h s {}_{1} & - \frac{2 r_{1}}{ρ^{2}} s {}_{1}- \frac{r_{1}^{2}}{ρ^{2}} h c {}_{1} \end{matrix}]

（52）

其中： $\bar{s} = s i n (r_{1} z)$ ； $\bar{c} = c o s (r_{1} z)$ ； $ρ = 1 + r_{1}^{2} m$ ； $m = H / K$ ； $r_{1} = \sqrt[]{P / H_{1}}$ ； $H_{1} = H (1 - P / K)$ ； $\hat{ρ} = 1 -$ $r_{1}^{2} m$ ； $h = β - α$ ； $c_{1} = c o s (r_{1} h)$ ； $s_{1} = s i n (r_{1} h)$ 。

位移和旋转角的规律由式（51）给出，而区间 $[α, β]$ 中弯矩、等效剪切力和等效分布荷载 $\bar{f} (x)$ 的规律可以通过考虑以下列因素进行推导^［

18］，即

\{\begin{array}{l} \bar{M} = H \bar{ψ^{'}} \\ \bar{Q} = - H {\bar{ψ}}^{″} \\ \bar{V} = - H {\bar{ψ}}^{″} - P \bar{w^{'}} \\ \bar{f} = H {\bar{ψ}}^{‴} + P {\bar{w}}^{″} \end{array}

（53）

上述等效分布式荷载由位移函数 $\bar{U} = {[\begin{matrix} \bar{w} & \bar{ψ} \end{matrix}]}^{T}$ 推导得到，而不是通过定义区间 $[α, β]$ 荷载 $f$ 的函数在形函数 $N_{1 i}$ 二维空间上的正交投影得到。

在每个区间或单元中，作用力 $f (x)$ 由4个函数l， $z$ ， $s i n (r_{1} z)$ 和 $c o s (r_{1} z)$ 的线性组合表达式给出。上述方法在节点处以及每个单元内部都能提供位移、旋转角和应力规律的精确解。

4 实际案例分析

4.1　算例1

试件为长6 m、矩形截面恒定的梁柱，其嵌入⁃铰接部件如图3所示。理想材料特性如下：弹性模量 $E = 2.1 \times 10^{7} k N / m^{2}$ ；泊松比 $v = 0.3$ ；剪切校正因子 $k = 5 / 6$ ；轴向加载P=15 000 kN。通过长度L=6 m的单个单元，可得到嵌入⁃铰接部件的近似解，并将这种解与精确解进行比较。为了确定精确解，根据荷载定义，将梁柱离散成4个有限元：单元1为［0，2］；单元2为［2，3］；单元3为［3，4］；单元4为［4，6］。

图3 嵌入-铰接部件

Fig.3 Embedded-hinged component

另外，还进行了线性有限元简化积分计算^［

5］。为了获得与等效分布荷载⁃有限元方法相似的计算结果，至少需要20个位移、旋转角和剪力单元以及40个弯矩单元。精确解包含4个单元的等效分布荷载，近似解包含1个单元的等效分布荷载。位移、旋转角、弯矩、剪切力以及等效分布荷载与原始荷载的比较图分别如图4~8所示。

图4 位移比较图

Fig.4 Diagram of displacement comparison

图5 旋转角比较图

Fig.5 Diagram of rotation comparison

图6 弯矩比较图

Fig.6 Diagram of bending moment comparison

图7 剪切力比较图

Fig.7 Diagram of shear force comparison

图8 等效分布荷载与原始荷载的比较图

Fig.8 Diagram of equivalent distributed load and original load comparison

由图4，5可知，等效分布荷载、简化积分之间的位移和旋转角的近似解与其精确解的数值变化趋势相一致。由图6可知，使用等效分布荷载法和1个单元，弯矩的最大误差发生在横坐标3 m处，约为3.5%，其误差比精确解的值小2.5%；使用简化积分法和40个单元，其弯矩最大误差发生在梁柱铰接部件的下端2.5 m处，约为6%。由图7可知，使用简化积分法和20个单元时，剪应力精确解与近似解的最大误差均发生在荷载的施加点位上，这是因为该近似解是通过某种方式平均了这些点位左右的剪应力值，以消除不连续的情况来形成规律，因此当近似剪应力由一个连续函数给出时，剪应力误差最小。由图8可知，等效分布荷载的最大偏差均发生在荷载的施加点位上，梁柱铰接部件底端偏差为0，顶端发生相应偏差。在点位上施加的荷载越大则两者偏差会越大，且梁柱顶端偏差也随之增大。

铁木辛柯柱和伯努利⁃欧拉柱的结果比较如图9所示。可以看出：位移的最大相对误差发生在横坐标3 m处（即梁柱铰接部件中间的位置），约为26.5%；旋转角的最大相对误差发生在梁柱铰接部件上端，约为10.5%；弯矩的最大相对误差发生在横坐标3 m处，约为5.5%；剪切力的最大相对误差发生在梁柱铰接部件上端，小于5%。

图9 铁木辛柯柱和伯努利-欧拉柱的结果比较

Fig.9 Comparison results between Timoshenko column and Bernoulli-Euler column

倾斜度（ $L / b$ ）对铁木辛柯和伯努利⁃欧拉模型在梁柱位移方面的影响较大。为此，保持试件的荷载、长度和厚度相同，改变矩形截面短边 $b$ 的长度，其位移差随倾斜度的变化如图10所示。

图10 铁木辛柯和伯努利-欧拉之间的位移差随倾斜度的变化

Fig.10 Variation of displacement difference between Timoshenko and Bernoulli-Euler with inclination

由图10可知，当倾斜度降低时，铁木辛柯和伯努利⁃欧拉模型之间的差异增大，反之则减小；当倾斜度为12时，相对差异为7%；倾斜度为3时，差异十分明显，达到45%。

4.2　算例2

嵌入底座的柱子如图11所示。该试件为长14 m、矩形截面恒定的梁柱，中间有支撑物。理想材料的特性如下：泊松比 $v = 0.3$ ，弹性模量 $E = 3 \times 10^{7} k N / m^{2}$ ，剪切校正因子 $k = 5 / 6$ 。轴向加载P=25 000 kN的值约为铁木辛柯的 $60 % P_{c r T}$ 。采用横向荷载分布，即除了点荷载外，还结合了三角形荷载和其他恒定分布荷载，以测试该方法的性能。通过长度L=14 m的单个单元，单元中分布的荷载由一个多项式部分和三角部分的函数给出，可得到嵌入底座柱子的近似解。

图11 嵌入底座的柱子

Fig.11 Column embedded in base

将算例2获得的精确解和近似解进行比较。在给定的连接条件下，使用等效分布荷载⁃有限元的近似解，将该部分离散成2个单元：单元1为［0，7］；单元2为［7，14］。此外，使用线性有限元和简化积分的近似解，将试件离散为20个位移、旋转角和剪切力单元以及40个弯矩单元。在这种情况下，所提方法的精确解需要将梁柱离散成4个有限元，根据荷载定义为：单元1为［0，3.5］；单元2为［3.5，7］；单元3为［7，10.5］；单元4为［10.5，14］。图12~15分别为等效分布荷载和简化积分法的位移、旋转角、弯矩和剪切力比较图。

图12 等效分布荷载和简化积分法的位移比较图

Fig.12 Comparison of displacements for equivalent distributed loads and simplified integral method plots

图13 等效分布荷载和简化积分法的旋转角比较图

Fig.13 Comparison of angles of rotation for equivalent distributed loads and simplified integral method plots

图14 等效分布荷载和简化积分法的弯矩比较图

Fig.14 Comparison of bending moments for equivalent distributed loads and simplified integral method Fig.

图15 等效分布荷载和简化积分法的剪切力比较图

Fig.15 Comparison of shear forces for equivalent distributed loads and simplified integral method plots

由图12，13可知，等效分布荷载、简化积分之间的位移和旋转角的精确解和近似解高度吻合。由图14可知，最大弯矩值发生在中间有支撑物处，使用等效分布荷载和2个单元（近似解）的最大弯矩误差发生在横坐标3.5 m位置（梁柱下端），约为2.3%，其误差与精确解的值相比小1.2%；使用线性插值的简化积分和40个单元的最大弯矩误差发生在横坐标7 m位置（梁柱中间支撑物处），约为2.5%。

由图15可知，精确解、等效分布荷载（近似解）及简化积分之间的剪切力的最大误差均发生在荷载的施加点与支撑点上。通过两端的等效分布作用得到的剪切力和弯矩的值与对应的精确解相吻合，即拉格朗日意义上（对于剪切力）和赫尔墨特意义上（对于弯矩）的插值属性都得到了满足。

等效分布荷载与原始荷载的比较如图16所示。由图16可知，其最大偏差均发生在荷载的施加点上（梁柱顶端与中间支撑点除外）。施加的荷载越大两者偏差也越大，梁柱顶端与底端的偏差也随之增大。等效的荷载分布在整个梁柱单元内是连续的。相比原始荷载，等效荷载近似解在梁柱单元内集中荷载加载处连续可导，具有更好的适应性。

图16 等效分布荷载与原始荷载的比较

Fig.16 Comparison of equivalent distributed load and original load

5 结束语

本研究基于铁木辛柯梁理论，提出具有精确节点位移解的等效分布荷载⁃有限元算法，并对具有剪切变形的梁柱进行了受力分析。通过与其他算法数值计算结果的对比分析，验证了该算法的精确解。在2个不同的工况算例中，所提算法分别使用1个和2个单元就可以在位移、旋转角和剪切力方面获得与待解问题精确解较为接近的近似解。倾斜度（ $L / b$ ）对铁木辛柯和伯努利⁃欧拉模型在梁柱位移方面的影响较大，柱体矩形截面短边长度越小，等效分布荷载⁃有限元算法越有优势。此外，梁柱体底部与中间在不同的连接条件下，通过该方法可以确定铁木辛柯梁柱的弯曲荷载和稳定性函数，以消除不连续的情况，有效提高了柱体结构剪切变形有限元分析解的精确性。

参考文献

HUGHES T J R. The finite element method： linear static and dynamic finite element analysis［M］. Massachusetts： Courier Corporation Press， 2003： 103-132. [百度学术]

LIU G R， NGUYEN-THOI T， LAM K Y. A novel alpha finite element method （αFEM） for exact solution to mechanics problems using triangular and tetrahedral elements［J］. Computer Methods in Applied Mechanics and Engineering， 2008， 197（45）： 3883-3897. [百度学术]

WANG B， CAI Y， LI Z C， et al. Stochastic stable node-based smoothed finite element method for uncertainty and reliability analysis of thermo-mechanical problems［J］. Engineering Analysis with Boundary Elements， 2020， 114： 23-44. [百度学术]

HIRWANI C K， PANDA S K. Nonlinear finite element solutions of thermoelastic deflection and stress responses of internally damaged curved panel structure［J］. Applied Mathematical Modelling， 2019， 65： 303-317. [百度学术]

NACHBAGAUER K， PECHSTEIN A S， IRSCHIK H， et al. A new locking-free formulation for planar， shear deformable， linear and quadratic beam finite elements based on the absolute nodal coordinate formulation［J］. Multibody System Dynamics， 2011， 26（3）： 245-263. [百度学术]

De ANGELO M， PLACIDI L， NEJADSADEHI N， et al. Non-standard timoshenko beam model for chiral metamaterial： identification of stiffness parameters［J］. Mechanics Research Communications， 2020， 103： 103462. [百度学术]

BALOBANOV V， NIIRANEN J. Locking-free variational formulations and isogeometric analysis for the timoshenko beam models of strain gradient and classical elasticity［J］. Computer Methods in Applied Mechanics and Engineering， 2018， 339： 137-159. [百度学术]

RAJASEKARAN S， KHANIKI H B. Free vibration analysis of bi-directional functionally graded single/multi-cracked beams［J］. International Journal of Mechanical Sciences， 2018， 144： 341-356. [百度学术]

王续宏，李翔宇. 轴向受载铁木辛柯梁的弯曲研究［J］. 成都大学学报（自然科学版）， 2019， 38（2）： 187-189， 205. [百度学术]

WANG Xuhong， LI Xiangyu. Study on bending of axially loaded Timoshenko beams［J］. Journal of Chengdu University（Natural Science）， 2019， 38（2）： 187-189， 205.（in Chinese） [百度学术]

MAO Q， MAZZOTTI M， DEVITIS J， et al. Structural condition assessment of a bridge pier： a case study using experimental modal analysis and finite element model updating［J］. Structural Control and Health Monitoring， 2019， 26（1）： e2273. [百度学术]

ABRAMOVICH H，程桂胜. 轴向压力作用下铁木辛柯梁的固有频率［J］. 河北建筑工程学院学报， 1994（2）： 87-93. [百度学术]

ABRAMOVICH H， CHENG Guisheng. Natural frequencies of timoshenko beams under axial pressure［J］. Journal of Hebei Institute of Architectural Engineering， 1994（2）： 87-93.（in Chinese） [百度学术]

HASHEMIAN M， FOROUTAN S， TOGHRAIE D. Comprehensive beam models for buckling and bending behavior of simple nanobeam based on nonlocal strain gradient theory and surface effects［J］. Mechanics of Materials， 2019， 139： 103209. [百度学术]

BAI R X， OU H Q， PENG K W， et al. Failure study of honeycomb sandwich structure with embedded part under axial pullout loading［J］. Thin-Walled Structures， 2020， 146： 106489. [百度学术]

刘华，叶见曙，俞博，等. 桥梁荷载横向分布系数计算方法［J］. 交通运输工程学报， 2009， 9（1）： 62-66. [百度学术]

LIU Hua， YE Jianshu， YU Bo， et al. Calculation method of transversal distribution factor for bridge load［J］. Journal of Traffic and Transportation Engineering， 2009， 9（1）： 62-66.（in Chinese） [百度学术]

ELTAHER M A， MOHAMED N， MOHAMED S A， et al. Periodic and nonperiodic modes of postbuckling and nonlinear vibration of beams attached to nonlinear foundations［J］. Applied Mathematical Modelling， 2019， 75： 414-445. [百度学术]

PALACIO-BETANCUR A， ARISTIZABAL-OCHOA J D. Statics， stability and vibration of non-prismatic linear beam-columns with semirigid connections on elastic foundation［J］. Engineering Structures， 2019， 181： 89-94. [百度学术]

EARIJ A， ALFANO G， CASHELL K， et al. Nonlinear three⁃dimensional finite⁃element modelling of reinforced⁃concrete beams： computational challenges and experimental validation［J］. Engineering Failure Analysis， 2017， 82： 92-115. [百度学术]

LI Q S. Exact solutions for buckling of non-uniform columns under axial concentrated and distributed loading［J］. European Journal of Mechanics A-Solids， 2001， 20（3）： 485-500. [百度学术]

基于荷载等效理论的柱体结构剪切变形解析 PDF

摘要

关键词

引 言