精确几何模型下杆的压弯大变形

Large Deformation of a Rod with Compress and Bending Moment Load Using Accurate Geometric Model

摘要: 采用精确的几何模型，研究了一端固定一端自由约束的细长弹性杆在外部荷载作用下的压弯变形问题。在大变形过程中的精确几何关系下，建立几何方程和平衡方程，建立机械荷载作用下弹性压弯杆几何精确的大变形基本方程。采用打靶法数值对该边值问题进行求解，获得了在数值意义下问题的精确解。绘出该模型下弹性杆的压弯大变形特性曲线。研究外部载荷对弹性杆弯曲大变形的影响，分析压弯杆的力学特性。

Abstract: Based on accurate geometric model, the large deformation of a rod with free-fixed edges subjected to compress and bending moment loads was studied. Considering the precise geometric relation-ships across the deformation process, the geometric equations and equilibrium equations were es-tablished, and the basis equations of exact geometry under large deformation with linear constitu-tive relationship were established. Using shooting method to numerically solve the boundary value problem, the exact solution in the numerical sense was obtained, and the characteristic figures of the elastic large deformation with compress and bending moment loads were drawn, the influence of the large deformation elastic rod with external load was studied, and the mechanical characteris-tics of the rod were analyzed.

1. 引言

弹性杆在机械载荷作用下的过屈曲和变形问题在工程实践中具有极广泛的应用。例如作为钻杆、电子信号发射连接件、机器人臂、跳高撑杆等，近30年，在生命科学领域之中将弹性大变形杆作为生物大分子链(DNA等)的力学模型 [1] [2]。

1744年Euler等 [3] 研究了弹性细杆的平衡和稳定性问题。1927年Love等 [4] 进一步研究了细长杆的弹性曲线问题，细长弹性杆的几何非线性问题作为结构大变形的问题之一始终受到广泛重视。Kicrhhoff理论(即弹性杆静力学非线性理论)建立于1859年。Kicrhhoff理论其实质为利用欧拉-泊松(Euler-Poison)方程以及弹性杆平衡微分方程两者在数学形式中的相似性，将拉格朗日情形刚体定点转动的解析积分进行移植，以用来解决弹性杆的平衡问题 [5]。Pflǜger [6] 是第一位研究压杆在轴向压缩变形下对屈曲变形影响的人。Antanackvic [7] 等采用了分析方法研究在端部压力作用下可伸长简支梁过屈曲的分叉解，并得出了和细长比相关的临界载荷解析式。Gummadi和Palazotto [8] 对梁、拱中的大变形问题进行了研究，对现代大变形力学问题的研究产生了极大的影响。基于连续介质一般理论，Filipich和Rosales [9] [10] 系统推导出弹性细杆在轴向压力作用下的几何精确非线性控制方程，推导出了弹性压杆在9种本构关系下的过屈曲模型以及数值解，并研究在轴向载荷作用下弹性杆本构关系不同时对过屈曲行为的影响。Plaut等 [11] 基于弹性杆的轴线不可伸长几何非线性理论，采用打靶法对问题数值求解，分析了在外载荷作用下预弯曲杆在与刚性平面发生接触后的屈曲和变形问题。

薛纭等 [12] 和刘延柱等 [13] 基于弹性杆的几何非线性模型，同时考虑了横向剪切变形和杆的轴线伸长，得出了较为精确的Cosserat弹性杆模型。李世荣 [14] [15] [16] 等基于轴线可伸长杆的基本模型，采用了打靶法进行数值求解，对弹性压杆在几种不同边界条件下的过屈曲问题进行分析。郭锐 [17] 等研究了在中点附近点约束下的弹性细杆两端受压时的过屈曲变形，并给出了与轴向载荷和约束力相关的弹性杆的平衡构形以及平衡路径。李银山 [18] 等在研究弹性压杆的大变形问题时，采用了曲线积分的方法，得到了压杆过屈曲大变形的理论结果。薛纭等 [19] 于2019年基于弹性细长杆Kirchhoff动力学的比拟方法，从Kirchhoff弹性杆以及Cosserat弹性杆的拉格朗日方程出发，研究出了在弯扭度以及形心应变矢发生突变时的拉格朗日方程，并分析了弯扭度与形心应变矢突变对弹性杆挠曲线的影响。为弹性细长杆弯扭度发生突变提供了力学方法。刘荣刚等 [20] 采用能量法和挠曲线近似微分方程方法得出了临界力极限表达式，解决了细长弹性压杆临界力作用下的挠度不确定问题。

本文拟采用可伸长杆的基本理论，在精确的几何模型下推导线弹性变形杆的几何方程和平衡方程，并用打靶法数值求得压弯变形的数值解。

2. 控制方程

考虑一长为 $l$ 的矩形截面 $\left(b\times h\right)$ 压杆，一端固定一端自由，在自由端作用集中力 $F_0$ 和弯矩 $m_0$ (见图1)。

2.1. 几何方程

将杆的轴线置于直角坐标系 $oxy$，记变形前杆内任意一点为 $\left(x,y\right)$，轴线上的点 $\left(0,x\right)$， $x\in \left(0,l\right)$，杆发生弯曲变形后，轴线上点的坐标变为 $\left(w,x+u\right)$，其中 $u$ 和 $w$ 分别为 $x,y$ 方向的位移，假设变形服从Kirchhoff直线法假设，由轴线可伸长梁的大变形理论，可得如下几何关系：

$\frac{\text{d}s}{\text{d}x}=\lambda$ (1)

$\frac{\text{d}u}{\text{d}x}=\lambda \cos \theta -1$ (2)

$\frac{\text{d}w}{\text{d}x}=\lambda \text{sin}\theta$ (3)

$\frac{\text{d}\theta }{\text{d}x}=\frac{\text{d}\theta }{\text{d}s}\times \frac{\text{d}s}{\text{d}x}=\lambda \frac{1}{\rho }=\lambda k$ (4)

其中： $\theta$ 为横截面转过的角度，梁的轴线在变形后为平面曲线， $s$ 为变形后的弧长， $k$ 为轴线弯曲后的曲率， $\lambda$ 为轴线的伸长比。

Figure 1. Mechanical model of structure

图1. 结构的力学模型

轴线应变 ${\epsilon }_0$ 为

${\epsilon }_0=\frac{\text{d}s-\text{d}x}{\text{d}x}=\lambda -1$

横截面内任意一点的应变为

$\epsilon =\frac{\text{d}{s}^{\prime }-\text{d}x}{\text{d}x}=\frac{\text{d}\theta \left(\rho -y\right)}{\text{d}x}-1=\lambda -1-y\lambda k$

2.2. 物理方程

考虑应力—应变的线弹性关系，横截面上的正应力为

$\sigma =E\epsilon$

横截面上的轴力 $F_N$ 和弯矩 $m$ 为

$F_N={\displaystyle {\int }_{-h/2}^{h/2}\sigma b\text{d}y}=EA\left(\lambda -1\right)$

$m={\displaystyle {\int }_{-h/2}^{h/2}-\sigma by\text{d}y}=EI\lambda k$

其中： $A=b\times h$ 为杆的横截面面积， $I=b{h}^3/12$ 为横截面的抗弯刚度。

2.3. 平衡方程

在精确几何模型下研究变形后微段 $\text{d}x$ 上的受力平衡。用 $F_v$ 表示横截面上的内力，则 $F_v=F_0$，轴力和剪力分别为

$F_N=F_v\cos \theta ,F_S=F_v\sin \theta$

利用平衡条件可得以下平衡方程

$\frac{\text{d}{F}_N}{\text{d}x}=-F_v\sin \theta \frac{\text{d}\theta }{\text{d}x}=-F_N\lambda k\tan \theta$

$\frac{\text{d}m}{\text{d}x}=F_v\frac{\text{d}w}{\text{d}x}=F_N\lambda \tan \theta$

将物理方程带入以上平衡方程中得

$\frac{\text{d}\lambda }{\text{d}x}=-\left(\lambda -1\right)\lambda k\tan \theta$ (5)

$\frac{\text{d}k}{\text{d}x}=\left(\lambda -1\right)k^2\tan \theta +\frac{\left(\lambda -1\right)A\tan \theta }{I}$ (6)

式(1)~(6)构成了弹性杆在精确几何模型下非线性大变形问题下的基本方程。考虑具体边界条件，就构成杆件压弯大变形的边值问题，对于一端固定另一端自由的弹性杆，相应的边界条件为

$固定端\left(x=0\right):u=w=\theta =s=0$ (7)

$自由端\left(x=l\right):F_v=-F_0,m=m_0$ (8)

文献 [9] [10] 研究了弹性杆在端部压力和剪力作用下的过屈曲问题。基于Kirchhoff假设推导出的非线性基本方程不适用于弯曲变形过大的情况。

3. 对控制方程无量纲化

引入下列无量纲量

$\left(X,S,U,W\right)=\left(x,s,u,w\right)/l$

$K=k\cdot l,{\rm T}=\theta ,\Lambda =\lambda ,\delta =h/l$

$N=F_N/EA{\delta }^2,M=ml/EI\delta$

$V_0=F_0/EA{\delta }^2,M_0=m_{0}l/EI\delta$

其中： $\delta$ 为杆的细长比。将无量纲量代入物理方程可得无量纲轴力和弯矩的表达式

$N=\left(\Lambda -1\right)/{\delta }^2,M=\Lambda K/\delta$ (9)

将无量纲量代入式(1)~(8)，可得以下无量纲基本方程

$\frac{\text{d}S}{\text{d}X}=\Lambda$ (10)

$\frac{\text{d}{\rm T}}{\text{d}X}=\Lambda K$ (11)

$\frac{\text{d}U}{\text{d}X}=\Lambda \cos {\rm T}-1$ (12)

$\frac{\text{d}W}{\text{d}X}=\Lambda \text{sin}{\rm T}$ (13)

$\frac{\text{d}\Lambda }{\text{d}X}=-\left(\Lambda -1\right)\Lambda K\tan {\rm T}$ (14)

$\frac{\text{d}K}{\text{d}X}=\left(\Lambda -1\right)K^2\tan {\rm T}+\frac{\left(\Lambda -1\right)A\tan {\rm T}}{I}$ (15)

无量纲边界条件为

$固定端(x=0):S={\rm T}=U=W=0$ (16)

$自由端\left(x=1\right):\left(\Lambda -1\right)+V_0{\delta }^2\cos {\rm T}=0,\Lambda K-M_0\delta =0$ (17)

4. 弯曲问题的数值分析

由于微分方程的强非线性以及多未知量间的耦合效应，难以获得问题的解析解，因而采用打靶法对此边值问题进行数值求解。取 $\delta =0.1$，研究自由端的位移随端部无量纲载荷的变化规律。

采用打靶法求解，令：

$Y={\left\{y_1,y_2,y_3,y_4,y_5,y_6\right\}}^{\text{T}}={\left\{S,U,W,{\rm T},\Lambda ,K\right\}}^{\text{T}}$ (18)

问题的控制方程可写为一阶常微分方程组：

$\frac{\text{d}Y}{\text{d}X}=F\left(X;Y;\delta ,V_0,M_0\right)=\left\{y_5,U,W,{\rm T},\Lambda ,K\right\}$ (19)

$\frac{\text{d}{y}_1}{\text{d}x}=y_5,\frac{\text{d}{y}_2}{\text{d}x}=y_5\cdot \cos \left(y_4\right)-1,\frac{\text{d}{y}_3}{\text{d}x}=y_5\cdot \sin \left(y_4\right),\frac{\text{d}{y}_4}{\text{d}x}=y_5\cdot y_6$ (20)

$\frac{\text{d}{y}_5}{\text{d}x}=y_6\cdot y_5\cdot \left(1-y_5\right)\cdot \tan \left(y_4\right),\frac{\text{d}{y}_6}{\text{d}x}=\left({\delta }^2\cdot y_6{}^2+12\right)\cdot \left(1-y_5\right)\cdot \tan \left(y_4\right)/{\delta }^2$ (21)

边界条件为：

$固定端:y_1=0,y_2=0,y_3=0,y_4=0$ (22)

$自由端:\left(y_5-1\right)-V_0\cdot \text{cos}\left(y_4\right)\cdot {\delta }^2=0,y_5\cdot y_6-M_0\cdot \delta =0$ (23)

当端部载荷增加时，可以求解出弯曲程度较大的情形，说明方程(1)~(6)采用了精确的几何模型，对问题的求解适应大变形的弹性变形问题。文献 [18] 采用了微分求积法，建立大变形时弹性压杆的数学模型，并采用Maple编程对此模型进行了求解，得到自由端压力，转角，挠度和轴向位移的对应数值关系，其数值结果与本文结果基本一致。

图2为端部弯矩 $M_0$ 为0.1时，不同轴向压力 $V_0$ 的平衡构形。随着杆端轴向压力 $V_0$ 增加，弹性杆弯曲程度逐渐增大，并发生回弯(横向位移减小)。杆件的变形过程与文献 [18] 的结果是一致的。

图3给出了不同 $V_0$ 时，轴向位移 $U_1$ 随 $M_0$ 的变化规律。图中可以看出，随着 $M_0$ 的增加，轴向位移不断(负向)增加，在 $M_0$ 较大时，不同 $V_0$ 的变化规律是类似的，而在 $M_0$ 较小时，差异较大。 $V_0\le 0.2$ 三条曲线，轴向位移 $U_1$ 几乎从0开始，当 $V_0=0.25$ 时， $U_1$ 变为从-0.36开始。压杆稳定的无量纲临界载荷的Euler公式为

$V_{cr}=\frac{{\pi }^{2}EI}{{\left(\mu l\right)}^2}\cdot \frac{1}{EA{\delta }^2}=0.\text{2}0\text{56}$

在 $M_0=0$ 时， $V_0$ 达到临界载荷前，压杆的轴向位移几乎为0，而达到临界载荷后，轴向位移明显增加，变形从轴向压缩变形过渡到过屈曲变形。

Figure 2. Equilibrium configurations at $M_0=0.1$

图2. $M_0=0.1$ 的平衡构形

Figure 3. Variation curves of $U_1$ with $M_0$

图3. $U_1$ 随 $M_0$ 的变化曲线

图4给出了不同 $V_0$ 时，横向位移 $W_1$ 随 $M_0$ 的变化规律，当 $V_0$ 小于临界载荷时，曲线从0附近开始增加，当 $V_0$ 大于临界载荷时处于过屈曲平衡，此时 $M_0=0$ 时已经出现了横向位移。 $M_0$ 较小时，曲线在 $V_0=0.4$ 附近出现变化，此时变形开始回弯(见图2)。当 $V_0$ 较大时，过屈曲变形较大，挠度随 $M_0$ 增加而减小，回弯现象更为明显。

Figure 4. Variation curves of $W_1$ with $M_0$

图4. $W_1$ 随 $M_0$ 的变化曲线

图5给出了不同 $V_0$ 时，轴线长度随 $M_0$ 的变化规律。图中可以看出，在轴向压力作用下，随着 $M_0$ 的增加，轴线的长度是增加的，图像说明随着弯曲变形的增加， $V_0$ 引起轴线变形量下降。

Figure 5. Variation curves of $S_1$ with $M_0$

图5. $S_1$ 随 $M_0$ 的变化曲线

图3~5反映出： $M_0$ 较小时对杆变形影响在前屈曲阶段 $\left(V_0<V_{cr}\right)$ 和过屈曲阶段 $\left(V_0\ge V_{cr}\right)$ 是不同的，此时 $V_0$ 对变形过程影响明显。而当 $M_0$ 较大时，前屈曲阶段和过屈曲阶段的变形规律基本一致， $V_0$ 对变形过程的影响减弱。

Figure 6. Variation curves of $U_1$ with $V_0$

图6. $U_1$ 随 $V_0$ 的变化曲线

图6给出不同端部外力偶 $M_0$ 作用时，轴向位移随 $V_0$ 的变化曲线，图中可以看出，随着 $V_0$ 的增加，端部轴向位移是(负向)增加的。在 $V_0$ 较小时，不同的 $M_0$ 对应的变化曲线差异较大， $M_0<3$ 时，能明显的看出由弯曲变形到过屈曲阶段的过程，临界载荷在 $V_0=0.2$ 附近；而 $M_0>3$ 时， $M_0$ 引起的变形较大，前屈曲到过屈曲的变化并不明显。当端部压力 $V_0>1.5$ 时，不同 $M_0$ 引起的轴向位移差异变得很小。

Figure 7. Variation curves of $W_1$ with $V_0$

图7. $W_1$ 随 $V_0$ 的变化曲线

图7给出了不同 $M_0$ 时，横向位移 $W_1$ 随 $V_0$ 的变化规律，图中可以看出，随 $V_0$ 的增大，横向位移先增后减， $W_1$ 存在最大值。随 $M_0$ 增大最大挠度点横向偏移，对应的 $V_0$ 在数值上减小。当 $M_0=0.1$ 时，能明显看出随 $V_0$ 的增加，变形由前屈曲到过屈曲阶段的变化过程，当 $M_0$ 较大时，此过程不明显。

Figure 8. Variation curves of $S_1$ with $V_0$

图8. $S_1$ 随 $V_0$ 的变化曲线

图8给出了不同 $M_0$ 时，轴线长度 $S_1$ 随 $V_0$ 的变化规律。图中可以看出，随着 $V_0$ 增大， $S_1$ 呈现先减后增的趋势。 $M_0$ 较小时， $S_1$ 减小到增大的转折点较明显(临界载荷附近)，当 $M_0$ 较大时， $S_1$ 从减到增，变化比较平缓。

图6~8反映出：当 $V_0$ 较大时，不同 $M_0$ 时的变形规律基本相同，当 $V_0$ 较小时，尤其是在临界载荷附近，不同 $M_0$ 时的变形过程差别较大。

5. 结论

本文导出了一端固定一端自由约束的弹性杆在外荷载下压弯变形问题的精确几何模型的基本方程，根据求得的数值解，绘出了弹性杆在端部外荷载作用下的压弯构形和非线性平衡路径，并精确分析了外载荷对弹性杆变形的影响。得出了以下结论：

1) 本文基本方程采用了大变形下精确的几何模型，既可求解压杆的过屈曲问题，也可精确求解杆件的大变形弯曲问题。

2) 由精确几何模型下的基本方程求得的小变形时的数值结果与参考文献中的结果基本一致，而本文的基本方程更为精确。

3) 弹性杆压弯变形时，杆件变形随 $M_0$ 的变化规律受到 $V_0$ 的影响，在 $M_0$ 较小时， $V_0$ 的影响比较明显。同时杆件变形随 $V_0$ 的变化规律也受到 $M_0$ 的影响，当在 $V_0$ 临界载荷附近以及小于临界载荷时， $M_0$ 对变形的影响较为明显；随 $V_0$ 的增大， $M_0$ 对变形的影响逐渐减弱。

基金项目

本文得到国家自然科学基金11862012的资助。