大变形欧拉梁作为一种基础的力学模型，其理论研究对于分析结构几何非线性特性具有重要意义。随着工程技术及材料科学的发展，欧拉梁大变形分析需求也日益迫切，一些典型的工程结构如大展弦比柔性机翼、细长导弹等都可以简化为该模型进行分析研究。Dado与Al-Sadder^[1]总结了梁大变形问题的主要分析方法，包括椭圆积分分析方法^[2-3]、打靶法^[4-5]、有限差分方法^[6]及有限元方法^[7-8]。其中，椭圆积分分析方法可以给出半解析形式，具备结果准确度高、形式简洁及计算效率高等优势。
考虑几何非线性因素，Bisshop和Drucker^[2]最早给出了受定向载荷梁弯曲大变形的椭圆积分计算形式。Midha^[9]、Kimball^[10]、Lau^[11]等将这一形式扩展到受力矩载荷分析及受混合载荷分析。在处理变几何参数梁问题时，这一算法需要复杂的迭代流程。Shinohara和Hara^[12]将这一算法扩展至曲梁变形分析。王恩惠^[13]和李银山等^[14]利用椭圆积分方法研究压杆稳定性及变形分析问题，能够准确计算失稳后的屈曲力学特征。当考虑弹性梁受随动载荷时，复杂性提高，利用椭圆积分方法分析的研究不多^[15]。Srpcic^[6]和Argyris^[8]等分别利用有限差分方法及有限元方法给出梁端部受随动载荷的变形解。Rao等^{[4, 16]}利用改进打靶法分析了这一问题，提高了计算效率。Navaee^[17]及Shavartsman^[18]等进一步改进分析流程，通过改变积分顺序将边值问题转化为初值问题，进行Runge-Kutta数值积分求解，降低分析难度。多数现有研究局限于载荷与梁轴线垂直的情况。
本文利用椭圆积分，给出考虑弹性欧拉梁受端部定向载荷，即所谓“死载”(Dead Force)及随动载荷(Follower Force)下大变形的分析求解方法，分析定向受压杆平衡分支解问题。利用不同的变量变换关系使这两类问题控制方程具有统一形式。不限制载荷与梁轴线的角度关系，提供适用程度更高的梁大变形分析方法。本文计算结果与非线性有限元方法计算结果进行了对比。
1 梁大变形方程 考虑根部固支悬臂梁，长度L，弯曲刚度EI(s)，端部受载P。悬臂梁采用各向同性的线弹性材料组成。梁单元平衡方程满足：

(1)

式中: κ为梁上某点变形曲率; φ为梁变形角; M(s)为梁弯矩内力; EI为梁弯曲刚度; s为弧长坐标。
1.1 端部受定向载荷 考虑在悬臂梁端部给定一方向不变、大小不变的载荷。受载下梁几何关系如图 1所示。图中：α为载荷方位角，φ为梁变形角。

图 1 定向载荷作用下梁几何关系 Fig. 1 Geometric relationship of beam with dead force

图选项

将式(1)微分可得控制方程：

(2)

边界条件为

(3)

给出变换关系：

(4)

代入式(2)、式(3)可得新的控制方程及边界条件为

(5)

(6)

积分式(5)，可得

(7)

代入边界条件整理可得

(8)

式中: 。令，整理可得

(9)

引入变量θ，令：

(10)

积分式(10)可得

(11)

(12)

令相对弧长坐标，载荷因子，即有变形与载荷关系：

(13)

(14)

F(k, θ)为第一类椭圆积分^[19]。当s=1时，梁端部变形与载荷关系有

(15)

1.2 端部受随动载荷 考虑在悬臂梁端部给定一方向与梁轴向夹角不变、大小不变的载荷，即随动载荷(Follower Force)。受载下梁几何关系如图 2所示。

图 2 随动载荷作用下梁几何关系 Fig. 2 Geometric relationship of beam with follower force

图选项

将式(1)微分可得控制方程：

(16)

边界条件为

给出变换关系：

(17)

代入式(16)、式(3)可得新的控制方程为

边界条件为

(18)

积分式(5)，可得

代入边界条件整理可得

式中：z_L=α。令，整理可得

引入变量θ，令：

积分式(10)可得

(19)

即有变形与载荷关系：

(20)

同样地，梁端部变形与载荷关系有

对比可见，定向载荷与随动载荷加载后的梁控制方程形式一致，其差别表现在求解方程的边界条件差异。
2 定向受压杆问题的平衡分支解 定向受压杆的屈曲问题是梁非线性变形研究中的经典问题，超临界状态下对应一组载荷存在多种平衡分支现象。在端部受定向载荷解中，给定载荷方位角α=0°，变换关系z(s)=φ(s)，变形与载荷关系为

(21)

式中: θ₁=mπ, m=0, -1, -2, …，方程均成立，即存在多平衡解分支。需要说明的是，由边界条件及变换关系可知，，因此m≤0。
当m=0时，称为平衡Ⅰ态，其临界载荷因子。载荷超过该临界载荷时，平衡解为平衡Ⅰ态解曲线，由式(11)及式(15)可知其解曲线由积分区间0≤θ≤单独确定。
当m=-1时，称为平衡Ⅱ态，其临界载荷因子。载荷超过该临界载荷时，平衡解出现新的分支解即平衡Ⅱ态解曲线，由式(11)、式(15)及引入变量变换关系可知，解曲线由积分区间确定。
当m=-2时，称为平衡Ⅲ态，其临界载荷因子。载荷超过该临界载荷时，平衡解在平衡Ⅰ态/Ⅱ态解曲线基础上出现新的分支解即平衡Ⅲ态解曲线，由引入变换关系可知，解曲线由积分区间-2π≤θ≤确定。以此类推，当载荷因子β介于平衡n态及平衡n+1态的临界载荷因子之间时，梁变形挠曲线共有n种形式，即n个平衡分支解。
值得说明的是，当梁端部随动载荷受压时，由式(20)可知k=0，解曲线为0，即随动载荷压杆不存在静不稳定解。
3 计算流程 梁端部受定向载荷时，求解式(13)~式(15)。由于椭圆积分参数k与梁末端变形角φ_L相关，无法直接给出。首先在式(15)中可以得到不同k值对应的载荷因子β。找到与工况指定的载荷因子对应的k即可确定椭圆积分参数及梁末端变形角。由式(13)即可得到弧长坐标s与变量z的关系。根据变换关系式(4)可恢复梁实际变形角。挠曲线结果可通过沿弧长坐标积分计算求解：

(22)

(23)

当s=1即s=L时，梁端部y方向变形量δ=y(1)，x方向变形量Δ=L-x(1)，梁挠曲线几何关系如图 3所示。计算流程如图 4所示。

图 3 梁变形挠曲线 Fig. 3 Deflection curves of beam

图选项

图 4 梁端部受定向载荷计算流程 Fig. 4 Flowchart of calculation with dead force at the end of beam

图选项

梁端部受随动载荷时，求解式(13)、式(15)、式(20)。椭圆积分参数仅与载荷方位角α相关，可以直接给出。在式(15)中可以得到不同φ_L值对应的载荷因子。找到与工况指定的载荷因子对应的φ_L即可确定梁末端变形角及中间变量θ₁。由式(13)即可得到弧长坐标与变量的关系。根据变换关系式(17)可恢复梁实际变形角。挠曲线结果同样可通过沿弧长坐标积分式(13)、式(21)计算求解。计算流程如图 5所示。

图 5 梁端部受随动载荷计算流程 Fig. 5 Flowchart of calculation with follower force at the end of beam

图选项

4 数值算例及结果 4.1 大变形梁模型 给定一柔性梁模型，梁截面圆形，不考虑沿轴向伸长/缩短，模型参数如表 1所示。
表 1 梁模型参数 Table 1 Parameters of beam model

参数	数值
长度L/m	1
截面半径r/m	0.001 7
弹性模量E/(N·m-2)	7×1010
泊松比	0.29

表选项






4.2 梁变形计算结果 给定不同的载荷方位角α后，指定不同的载荷因子β=PL²/(2EI)后，计算梁挠曲线变形。为验证计算结果的准确性，变形结果与非线性有限元计算结果进行对比。图 6给出了梁端部受定向载荷时不同α/β组合下的梁挠曲线计算结果对比。随着载荷因子增大，梁变形逐渐增大，本文方法计算结果与非线性有限元方法计算结果保持高度一致，且计算精度没有随载荷因子的增大及载荷方位角的改变而降低。

图 6 梁端部受定向载荷挠曲线 Fig. 6 Deflection curves with dead force at the end of beam

图选项

图 7给出了梁端部受随动载荷时不同α/β组合下的梁挠曲线计算结果对比。同样的，随着载荷因子增大，梁变形逐渐增大。计算结果与非线性有限元方法计算结果保持高度一致。非线性有限元计算通过MSC.Nastran软件实现，计算精度没有随载荷因子的增大及载荷方位角的改变而降低，计算方式适用性高。

图 7 梁端部受随动载荷挠曲线 Fig. 7 Deflection curves with follower force at the end of beam

图选项

图 8给出了梁端部变形在载荷方位角α=90°时受载荷因子影响的变化曲线，同时给出非线性有限元方法计算结果及文献[2]计算结果，三者保持高度一致。与相关文献结果^[2]的对比进一步验证了本文方法对于梁变形计算的准确性。表 2和表 3给出了定向载荷及随动载荷下，不同α/β组合下梁端部变形δ的具体数值计算结果对比。可以发现，其与非线性有限元方法计算结果偏差值很小，已有状态下，最大偏差值不超过0.75%。

图 8 定向载荷下梁变形对比(α=90°) Fig. 8 Contrast of beam deflection with dead force (α=90°)

图选项

表 2 定向载荷下梁端部变形δ对比 Table 2 Contrast of beam end deflection δ with dead force

载荷方位角/(°)	载荷因子	椭圆积分解/m	有限元解/m	偏差/%
90	0.5	0.301 63	0.301 81	0.06
90	1.0	0.493 5	0.494 96	0.3
90	2.0	0.670 65	0.671 49	0.12
90	4.0	0.787 85	0.788 11	0.03
45	0.5	0.295 45	0.295 58	0.04
45	1.0	0.582 65	0.584 92	0.39
45	2.0	0.797 67	0.798 19	0.65
45	4.0	0.847 02	0.843 12	0.46
135	0.5	0.178 43	0.178 44	0.01
135	1.0	0.280 25	0.281 05	0.28
135	2.0	0.388 11	0.388 32	0.05
135	4.0	0.480 46	0.480 58	0.02

表选项






表 3 随动载荷下梁端部变形δ对比 Table 3 Contrast of beam end deflection δ with follower force

载荷方位角/(°)	载荷因子	椭圆积分解/m	有限元解/m	偏差/%
90	0.5	0.324 8	0.324 64	0.05
90	1.0	0.576 22	0.579 23	0.52
90	2.0	0.785 65	0.784 63	0.13
90	4.0	0.623 69	0.623 01	0.11
45	0.5	0.221 02	0.221 0	0.01
45	1.0	0.386 93	0.385 97	0.25
45	2.0	0.584 41	0.583 68	0.12
45	4.0	0.652 39	0.651 96	0.15
135	0.5	0.252 09	0.252 06	0.01
135	1.0	0.492 37	0.496 02	0.74
135	2.0	0.793 19	0.792 17	0.13
135	4.0	0.495 42	0.495 06	0.07

表选项






图 9展示了受定向载荷时，梁端部变形在不同载荷方位角下，受载荷因子影响的变化曲线。梁端挠度随载荷增大呈现逐渐增大趋势，最终趋于平缓。图 10展示了受随动载荷时，梁端部变形在不同载荷方位角下，受载荷因子影响的变化曲线。与受定向载荷情况不同，梁端部挠度随载荷的增大呈现先增大后减小的趋势。这与两类载荷加载中梁局部剪力变化趋势不同密切相关。

图 9 定向载荷下梁变形 Fig. 9 Deflection of beam with dead force

图选项

图 10 随动载荷下梁变形 Fig. 10 Deflection of beam with follower force

图选项

4.3 定向受压杆平衡分支解计算结果 给定载荷方位角α=0°计算该算例模型前3项临界载荷因子及平衡态解曲线。由椭圆积分性质可得，有

(24)

结合载荷因子定义，临界载荷为

(25)

这与材料力学中心压杆失稳欧拉解一致。类似的，通过计算给定参数椭圆积分极值，可以得到第2及第3项临界载荷因子：

(26)

(27)

图 11~图 13给出了在不同载荷因子下，前3个平衡态解曲线计算结果，可以发现越靠后出现的平衡态其梁曲线越复杂。结合前述理论分析，平衡Ⅰ态解曲线由0≤θ≤π/2积分区间确定，其解曲线包括1个驻点、0个反弯点；平衡Ⅱ态解曲线由-π≤θ≤π/2积分区间确定，其解曲线包括2个驻点、1个反弯点；平衡Ⅲ态解曲线由-2π≤θ≤π/2积分区间确定，其解曲线包括3个驻点、2个反弯点；对于平衡n态，其解曲线将由-(n-1)π≤θ≤π/2积分区间确定，其解曲线包括n个驻点、n-1个反弯点。本文方法适用于计算定向压杆各级平衡态临界载荷及给定载荷下的解曲线。

图 11 平衡Ⅰ态解曲线 Fig. 11 Solution curves of balance branchⅠ

图选项

图 12 平衡Ⅱ态解曲线 Fig. 12 Solution curves of balance branch Ⅱ

图选项

图 13 平衡Ⅲ态解曲线 Fig. 13 Solution curves of balance branch Ⅲ

图选项

5 结论 1) 利用椭圆积分推导了悬臂梁端部受集中定向载荷及随动载荷的大变形分析方程，通过变量代换将两者形式进行了统一，使两者拥有一致的控制方程，仅在中间变量计算中存在不同。
2) 计算了给定载荷因子及载荷方位角组合下的梁大变形，展示了大变形的变化规律，并将结果与非线性有限元分析结果进行对比，两者一致性很高。当载荷逐渐增大时，非线性有限元方法收敛性逐渐变差，而椭圆积分解不会出现该问题。
3) 对定向受压杆平衡分支解问题进行了分析计算，给出了各级平衡态临界载荷及解曲线的计算方法，并简要分析了各级平衡态解曲线性质。后续可对不同载荷方位角下的平衡分支解问题进行深入研究。

参考文献

[1]	DADO M, AL-SADDER S. A new technique for large deflection analysis of non-prismatic cantilever beams[J]. Mechanics Research Communications, 2005, 32(6): 692-703. DOI:10.1016/j.mechrescom.2005.01.004
[2]	BISSHOP K E, DRUCKER D C. Large deflection cantilever beams[J]. Quarterly of Applied Mathematics, 1945, 3(3): 272-275. DOI:10.1090/qam/13360
[3]	TIMOSHENKO S P, GERE J M. Theory of elastic stability[M]. New York: McGraw-Hill, 1963: 125-127.
[4]	RAO B N, RAO G V. Applicability of static and dynamic criterion for the stability of a cantilever column under a tip-concentrated subtangential follower force[J]. Journal of Sound and Vibration, 1987, 120(1): 197-200.
[5]	WANG C M, KITIPORNCHAI S. Shooting optimization technique for large deflection analysis of structural members[J]. Engineering Structure, 1992, 14(4): 231-240. DOI:10.1016/0141-0296(92)90011-E
[6]	SRPCIC S, SAJE M. Large deflections of thin curved plane beam of constant initial curvature[J]. International Journal of Mechanical Science, 1986, 28(5): 275-287. DOI:10.1016/0020-7403(86)90041-X
[7]	WANG G, SHAHINPOOR M. Design prototyping and computer simulations of a novel large bending actuator made with a shape memory alloy contractile wire[J]. Smart Materials Structure, 1996, 6(2): 214-221.
[8]	ARGYRIS J H, SYMEONIDS S P. Non-linear finite element analysis of elastic systems under nonconservative loading-natural formulation.Part 1.Quasi-static problems[J]. Compute Methods of Applied Mechanical Engineering, 1981, 26(1): 75-123. DOI:10.1016/0045-7825(81)90131-6
[9]	HOWELL L L, MIDHA A. Parametric deflection approximations for end-loaded large deflection beams in compliant mechanic[J]. Journal of Mechanical Design, 1995, 117(1): 156-165. DOI:10.1115/1.2826101
[10]	KIMBALL C, TSAI L W. Modeling of flexural beams subjected to arbitrary end loads[J]. Journal of Mechanical Design, 2002, 124(2): 223-234. DOI:10.1115/1.1455031
[11]	LAU J H. Closed-form solutions for the large deflections of curved optical glass fibers under combined loads[J]. Journal of Electron Package, 1993, 115(3): 337-339. DOI:10.1115/1.2909337
[12]	SHINOHARA A, HARA M. Large deflection of a circular c-shaped spring[J]. International Journal of Mechanical Science, 1979, 21(3): 179-186. DOI:10.1016/0020-7403(79)90022-5
[13]	王恩惠. 弹性悬臂杆的平衡分支[J]. 固体力学学报, 1983, 4: 143-151. WANG E H. Balance branch of elastic cantilever rod[J]. Acta Mechanica Solida Sinica, 1983, 4: 143-151. (in Chinese)
[14]	李银山, 刘波, 潘文波, 等. 弹性压杆的大变形分析[J]. 河北工业大学学报, 2011, 40(5): 31-35. LI Y S, LIU B, PAN W B, et al. Analysis of large deflection of flexible compression bars[J]. Journal of Hebei University of Technology, 2011, 40(5): 31-35. DOI:10.3969/j.issn.1007-2373.2011.05.007 (in Chinese)
[15]	FRISCH-FAY R. Flexible bars[M]. London: Butterworths, 1962: 119-151.
[16]	RAO B N, RAO G V. Applicability of static or dynamic criterion for the stability of a non-uniform cantilever column subjected to a tip-concentrated subtangential follower force[J]. Journal of Sound and Vibration, 1988, 122(1): 188-191. DOI:10.1016/S0022-460X(88)80017-8
[17]	NAVAEE S, ELLING R. Equilibrium configurations of cantilever beams subjected to inclined end loads[J]. Journal of Applied Mechanics, 1992, 59(3): 572-579. DOI:10.1115/1.2893762
[18]	SHAVARTSMAN B S. Large deflections of cantilever beam subjected to a follower force[J]. Journal of Sound and Vibration, 2007, 304(3-5): 969-973. DOI:10.1016/j.jsv.2007.03.010
[19]	王竹溪, 郭敦仁. 特殊函数概论[M]. 北京: 北京大学出版社, 2000: 520-556. WANG Z X, GUO D R. Special functions[M]. Beijing: Peking University Press, 2000: 520-556. (in Chinese)