欧拉-伯努利梁平面超大挠性变形问题变分法求解

孔祥森，李昊，宋雷鹏，沈星; KONG Xiangsen; LI Hao; SONG Leipeng; SHEN Xing

2025年4月4日 8:23 星期五

网刊加载中。。。

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

确定继续浏览么?

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

欧拉‑伯努利梁平面超大挠性变形问题变分法求解 PDF

- ORCID：
孔祥森 ^1,2
- ORCID：
李昊 ²
✉
- ORCID：
宋雷鹏 ¹
- ORCID：
沈星 ¹

1. 南京航空航天大学航空学院, 南京 210016； 2. 上海卫星工程研究所，上海 200204

中图分类号： O302

最近更新：2023-03-06

DOI：10.16356/j.1005-2615.2023.01.008

摘要

对于欧拉‑伯努利悬臂梁平面超大挠性变形问题，由于其复杂的非线性几何方程，以位移为基本变量进行求解时，通常只能采用如多重打靶、微分求积等数值方法求得梁上离散点的位移值。本文研究了欧拉‑伯努利悬臂梁平面超大挠性变形问题变分法求解理论。通过假设多项式形式的梁的曲率试函数以及常数中心线应变，基于欧拉‑伯努利悬臂梁的基本假设，推导出了相互耦合的位移函数的精确表达式，并基于变分法理论和三角函数级数展开，推导出欧拉‑伯努利梁的非线性控制方程组。利用迭代法对非线性控制方程组中的未知参数进行求解，最终得到欧拉‑伯努利悬臂梁的位移函数的解析表达式。利用有限元计算结果对提出的变分法求解理论进行验证，并分别计算了欧拉‑伯努利悬臂梁在自由端集中力及位移约束情况下的大变形。算例表明，基于本文的变分法求解理论，利用6个未知参数，即能够精确预测欧拉‑伯努利悬臂梁在自由端集中力及位移约束下的超大挠性变形，该研究成果为欧拉‑伯努利悬臂梁的超大变形问题提供了新的求解方法。

关键词

欧拉‑伯努利悬臂梁; 超大挠性变形; 变分法; 非线性; 几何方程

欧拉‑伯努利梁在忽略梁截面的变形后，在外载荷作用下其中心线位移的预测成为欧拉‑伯努利悬臂梁的基本问题。对于发生超大挠性变形的梁结构，大变形将导致几何非线性，因此需要采用几何精确梁理论进行大挠度及大转动变形分析^［

1］。基于几何精确梁理论的非线性控制方程，可以通过多重打靶法^{［参考文献 2

百度学术}2］、求积元法^{［参考文献 3

百度学术}3］、微分求积法^{［参考文献 4

百度学术}4］、几何配点法^{［参考文献 5

百度学术}5］或有限元法^{［参考文献 6

百度学术}6］进行求解。但是，以位移为基本变量，采用多重打靶或几何配点等方法求解欧拉‑伯努利梁超大挠性变形问题时，求解结果为离散节点的位移值，无法给出连续的位移函数表达式。在梁、板变形问题的研究中，变分法是求得位移函数的常用方法。变分法通过选择合适的位移试函数，代入对应科学问题的泛函中，通过对泛函求驻值，以确定位移试函数中的待定参数，从而能够得到位移函数的近似解^{［参考文献 7

百度学术}7］。以非对称复合材料薄板的大变形研究为例，众多学者采用了经典的Von‑Karman非线性几何方程。通过假设多项式形式的面外位移函数及面内应变函数，经过简单的积分即可推导得到层合板的面内位移函数，进而建立非线性控制方程^{［参考文献 8‑10}8‑10］。然而，Von‑Karman非线性几何方程在推导过程中进行了简化，并且方程中的位移是固定坐标系的单值函数。因此，Von‑Karman非线性几何方程对于超大挠性变形问题已不再适用。对于欧拉‑伯努利悬臂梁平面超大挠性变形问题，几何精确梁理论虽然给出了梁中线应变及曲率与位移函数之间的非线性几何关系^{［参考文献 1

百度学术}1， 11‑12］，但没有如小变形理论或Von‑Karman非线性几何方程一样给出位移函数之间的相容关系。因此，当采用变分法求解欧拉‑伯努利梁的超大挠性变形问题时，位移函数的假设形式成为难点。当假设的位移试函数不满足相容性时，变分法的计算结果往往误差较大。理论上虽然可以通过大量增加位移试函数中的未知参数来降低位移函数不相容的影响，但大量的未知参数很容易导致多解，并使得求解时间呈指数形式增加。因此，对于欧拉‑伯努利悬臂梁平面超大挠性变形问题，尚未见变分法求解的报道。

本文提出了欧拉‑伯努利悬臂梁平面超大挠性变形问题的变分法求解方法。与经典几何精确梁理论不同的是，本方法以欧拉‑伯努利悬臂梁的曲率及中心面应变为基本变量，基于随体坐标系，通过假设欧拉‑伯努利悬臂梁的多项式形式的曲率试函数及常数中心线应变，推导出悬臂梁的位移函数，通过变分法求解，可精确预测欧拉‑伯努利悬臂梁在不同载荷作用下的超大弯曲变形。本文给出了超大弯曲变形欧拉‑伯努利悬臂梁变分法求解理论的详细推导过程，并采用有限元计算结果对变分法求解结果进行了正确性验证。

1 超大弯曲变形欧拉‑伯努利悬臂梁位移函数推导

为采用变分法求解欧拉‑伯努利悬臂梁的超大变形问题，本文基于随体坐标系，通过假设曲率方程及中线应变，逆向推导考虑相容关系的位移函数。欧拉‑伯努利梁遵循柯西霍夫假设，即变形前垂直于中面的直线在变形后依然垂直于中面，并且忽略悬臂梁厚度方向的变形。本文采用一个固定笛卡尔坐标系及一个正交随体坐标系，推导考虑相容性的超大挠性变形悬臂梁的位移函数。本文采用的坐标系如图1所示，其中x‑y‑z为笛卡尔坐标系，l‑b‑h为随体坐标系。当悬臂梁发生变形时，在任何一点随体坐标系的l轴与b轴相互垂直且与悬臂梁的中面重合，h轴垂直于悬臂梁的中面。

图1 超大弯曲变形欧拉‑伯努利悬臂梁坐标系示意图

Fig.1 Sketch map of coordinates of Euler‑Bernoulli beam with quite large deflection

基于随体坐标系，对于欧拉‑伯努利悬臂梁平面超大挠性变形问题，假设超大挠性变形悬臂梁的多项式弯曲曲率方程

\begin{array}{l} K_{1} [l] = a_{0} + \sum_{1}^{m} a_{m} l^{m} \\ K_{2} = 0 \\ K_{12} = 0 \end{array}

（1）

式中 $a_{m}$ 为未知参数。基于随体坐标系，假设悬臂梁为单轴应力状态，并假设超大弯曲悬臂梁中心线应变为常数

e [l] = c

（2）

式（1）及式（2）假设悬臂梁的横向弯曲曲率及扭曲曲率为0，悬臂梁只发生纵向弯曲变形及伸缩变形，忽略悬臂梁的横向变形。假设悬臂梁变形后在笛卡尔坐标系中的坐标为（u［l，b］，v［l，b］，w［l，b］），由于忽略横向变形，悬臂梁的坐标可简化为随体坐标l的函数（u［l］，b，w［l］）。对于任意位置l处长度为dl的悬臂梁，当dl趋于0时，在dl区间内悬臂梁的弯曲曲率可认为保持不变，即 $K_{1} [l]$ ，则悬臂梁在位置l处的曲率半径为 $R [l] = - 1 / K_{1} [l]$ ，如图2所示。图中 $θ [l]$ 为悬臂梁在随体坐标l处与固定坐标系x轴的夹角。 $θ [l]$ 可通过对曲率函数积分求得

θ [l] = - \int_{0}^{l} K_{1} [l] d l

（3）

图2 悬臂梁微单元坐标系示意图

Fig.2 Coordinate sketch map of a beam’s micro element

在笛卡尔坐标系中，初始长度为dl的悬臂梁单元变形后如图 2所示，其长度变为 $d l (1 + e [l])$ ，悬臂梁单元的起点与笛卡尔坐标系x轴的夹角为 $θ [l]$ 。将薄板单元绕着y轴顺时针旋转 $θ [l]$ 后，悬臂梁单元在其起始端与x轴相切，旋转后其末端坐标可表示为

\begin{array}{l} d u [l] = R [l] s i n [\frac{d l (1 + e [l])}{R [l]}] \\ d w [l] = R [l] - R [l] c o s [\frac{d l (1 + e [l])}{R [l]}] \end{array}

（4）

则初始长度为dl悬臂梁单元在旋转前的末端坐标，可通过旋转变换得到

\begin{array}{l} [\begin{matrix} d u [l] \\ b \\ d w [l] \end{matrix}] = [\begin{matrix} c o s [- θ [l]] & 0 & s i n [- θ [l]] \\ 0 & 1 & 0 \\ - s i n [- θ [l]] & 0 & c o s [- θ [l]] \end{matrix}] • \\ [\begin{matrix} R [l] s i n [\frac{d l (1 + e [l])}{R [l]}] \\ b \\ R [l] - R [l] c o s [\frac{d l (1 + e [l])}{R [l]}] \end{matrix}] \end{array}

（5）

根据导数的定义，悬臂梁上任意一点在笛卡尔坐标系坐标的导数可表示为

d [\begin{array}{l} d u [l] \\ b \\ d w [l] \end{array}] / d l = [\begin{matrix} \underset{d l \to 0}{l i m} \frac{d u [l]}{d l} \\ 0 \\ \underset{d l \to 0}{l i m} \frac{d w [l]}{d l} \end{matrix}]

（6）

将式（1~5）代入式（6）并求极限，最终将得到悬臂梁上任意一点在笛卡尔坐标系中坐标的导数 $u^{'} [l]$ 和 $w^{'} [l]$ 。

\begin{array}{l} u^{'} [l] = (1 + e [l]) c o s [\int K_{1} [l] d l] \\ w^{'} [l] = - (1 + e [l]) s i n [\int K_{1} [l] d l] \end{array}

（7）

将式（7）对坐标l进行积分，将得到悬臂梁在笛卡尔坐标系中的坐标，并能够得到悬臂梁在笛卡尔坐标系中的位移方程

\begin{array}{l} U [l] = \int (1 + e [l]) c o s [\int K_{1} [l] d l] d l - l \\ W [l] = - \int (1 + e [l]) s i n [\int K_{1} [l] d l] d l \end{array}

（8）

基于式（8），可推导得到悬臂梁基于随体坐标系的非线性几何方程

(1 + U^{'} {[l]}^{2}) + W^{'} {[l]}^{2} = (1 + e {[l])}^{2}

（9）

基于式（8），可推导得到悬臂梁的位移相容方程

\frac{1 + U^{'} [l]}{c o s [\int K_{1} [l] d l]} + \frac{W^{'} [l]}{s i n [\int K_{1} [l] d l]} = 0

（10）

针对欧拉‑伯努利悬臂梁，式（9）与几何精确梁理论给出的结果完全相同，但通过直接假设位移函数U［l］和W［l］进行变分法求解时，所假设的位移函数往往难以满足式（10）。针对悬臂梁，本文通过假设曲率函数K₁［l］及应变函数ε₁［l］，可以直接通过式（8）推导位移函数U［l］和W［l］的精确解析形式。

2 超大变形欧拉‑伯努利悬臂梁变分法求解

基于K₁［l］、e［l］、U［l］和W［l］，并基于几何精确梁的一维本构定律，可以推导出欧拉‑伯努利悬臂梁的位移变分方程

\begin{array}{l} δ P = ∭ ε [l] E δ (ε [l]) d V - \iint (F_{x} [l] δ (U [l]) + \\ F_{z} [l] δ W [l] + M_{y} δ θ [l]) d b d l = 0 \end{array}

（11）

式中： $P$ 表示悬臂梁的总势能； $E$ 为弹性模量；ε［l］为应变

ε [l] = e [l] + h K_{1} [l]

（12）

对于悬臂梁，由于其应变及曲率假设为多项式形式，则应变能密度变分ε［l］Eδ（ε［l］）也为多项式，可通过积分得到应变势能变分的精确形式。而对于其外力虚功部分，由于U［l］及W［l］的解析表达式中包括三角函数与多项式的乘积，难以直接进行积分运算。因此，为实现变分法求解，本文将对式（7）和（8）进行级数展开简化。

三角函数在 $x = 0$ 处的幂级数展开式为

\begin{array}{l} c o s [x] \approx \sum_{0}^{n} c o s^{(n)} [0] \frac{x^{n}}{n!} \\ s i n [x] \approx \sum_{0}^{n} s i n^{(n)} [0] \frac{x^{n}}{n!} \end{array}

（13）

将式（13）代入式（7）后，将坐标的导数u'［l］和w'［l］变为多项式形式

\begin{array}{l} u^{'} [l] = (1 + e [l]) \sum_{0}^{n} c o s^{(n)} [0] \frac{{(\int K_{1} [l] d l)}^{n}}{n!} \\ w^{'} [l] = - (1 + e [l]) \sum_{0}^{n} s i n^{(n)} [0] \frac{{(\int K_{1} [l] d l)}^{n}}{n!} d l \end{array}

（14）

将式（14）代入式（8），可得到多项式形式的位移函数

\begin{array}{l} U [l] = \\ \int (1 + e [l]) \sum_{0}^{n} c o s^{(n)} [0] \frac{{(\int K_{1} [l] d l)}^{n}}{n!} d l - l \\ W [l] = \\ - \int (1 + e [l]) \sum_{0}^{n} s i n^{(n)} [0] \frac{{(\int K_{1} [l] d l)}^{n}}{n!} d l \end{array}

（15）

将式（15）代入式（11），得到一组非线性方程组

f_{m} = \frac{\partial P}{\partial a_{m}}, f_{m + 1} = \frac{\partial P}{\partial c^{}}

（16）

对于悬臂梁，在悬臂梁上施加已知载荷 $F_{x}$ 、 $F_{z}$ 时，非线性方程组中方程数量等于式（1）及式（2）中未知参数的总和。当在悬臂梁上施加位移边界条件时，式（11）中的 $F_{x}$ 、 $F_{z}$ 、 $M_{y}$ 变为未知的支座反力。因此，为求得悬臂梁的变形，需将支座反力作为未知变量，并将式（16）与边界条件方程进行联合求解，能够同时得到式（1）及式（2）中的未知参数及支座反力 $F_{x}$ 、 $F_{z}$ 、 $M_{y}$ 的值。

本文采用了数学软件Mathematica进行上述公式的推导及数值计算，采用FindRoot函数进行非线性方程组的数值求解。

3 计算结果分析讨论

本文分别采用有限元方法和上述变分求解法，对欧拉‑伯努利悬臂梁的平面超大挠性变形进行预测。悬臂梁长度100 mm，截面宽度2 mm，厚度0.5 mm，材料为钢。计算中采用的材料参数为 $E = 200 G P a$ ， $μ = 0.3$ 。有限元仿真采用ABAQUS，悬臂梁采用2结点梁单元B31进行建模，同时在计算过程中开启几何非线性选项Nlgeom。

在悬臂梁末端施加集中力或弯矩，其变形预测结果对比如图3、4所示。在式（1，15）中，分别取参数 $m = 6, n = 10$ 。在末端集中载荷作用下，悬臂梁产生了超大挠性变形。图3、4中的计算结果显示，本文所提出的变分法求解结果与有限元方法预测结果几乎完全吻合。

图3 悬臂梁末端加集中力载荷变形预测结果(m=6,n=10)

Fig.3 Deformation prediction of cantilever beam with concentrate force loads at its end (m=6, n=10)

图4 悬臂梁末端加弯矩及集中力载荷变形预测结果(m=6,n=10)

Fig.4 Deformation prediction of cantilever beam with a concentrate force and a moment load at its end (m=6,n=10)

从式（13~15）可以看出，位移函数 $U [l]$ 及 $W [l]$ 为多项式形式，但其阶数随着三角函数的幂级数阶数n呈指数增加。为验证式（13）中三角函数展开阶数n对计算结果的影响，图5中给出了不同展开阶数n对应的结果。分析结果显示，三角函数幂级数展开阶数 $n \geq 8$ 时，n继续增加对计算结果精度的影响有限。为降低计算成本，可取 $n = 8$ 。

图5 三角函数幂级数展开阶数n影响分析(m=6)

Fig.5 Parametric study on influence of series number n of power series expansion (m=6)

为验证多项式曲率函数 $K_{1} [l]$ 的阶数对计算结果的影响，改变m的取值并求解，计算结果如图6所示。在图3及图4中的末端载荷作用下，当多项式的阶数 $m < 3$ 时，变分法的计算结果将出现明显误差。因此，取m=3并假设中心线应变e［l］为常数，则模型中只包含5个未知参数，即可精确预测悬臂梁在末端载荷作用下的超大变形。

图6 多项式曲率函数阶数m影响分析(n=8)

Fig.6 Parametric study on influence of order m of polynomial curvature (n=8)

在悬臂梁末端施加位移约束 $U [0.1] = - 0.05, W [0.1] = 0$ ，变分法计算结果与有限元计算结果对比见图7及表1所示。由图7显示，当曲率函数 $K_{1} [l]$ 的阶数为m=6时，变分法计算的悬臂梁形状与有限元预测结果几乎完全重合。当 $K_{1} [l]$ 的阶数m低于4时，变分法预测结果与有限元预测结果存在较为明显的误差。表1中支座反力的预测结果对比显示，当 $K_{1} [l]$ 的阶数m高于2时，变分法预测的支座反力与有限元预测结果相比误差小于1%。因此，在位移约束条件下，采用四阶多项式曲率函数 $K_{1} [l]$ 即能获得较高的预测精度，此时模型中仅包括6个未知参数，即 $a_{0} ~ a_{4}$ 与c。

图7 位移约束下悬臂梁变形的变分法与有限元预测结果对比(n=8)

Fig.7 Comparison of results predicted by variational method and FEA of cantilever beam deformation under displacement constraints (n=8)

表1 悬臂梁支座反力变分法和有限元预测结果对比

Table 1 Comparison of results predicted by variational method and FEA on the end reaction forces of the restrained cantilever beam

支座反力/N	FEA	变分法
支座反力/N	FEA	$m = 6$	$m = 5$	$m = 4$	$m = 3$	$m = 2$
$F_{x}$	-9.564	-9.543	-9.548	-9.603	-9.632	-9.925
$F_{z}$	5.584	5.594	5.598	5.616	5.626	6.084

上述变分法与有限元预测结果对比说明，本文所提出的变分法，使用极少的未知参数即能够精确预测悬臂梁在不同末端边界条件下的超大挠曲变形。对于更加复杂的边界条件，可通过适当增加多项式曲率函数 $K_{1} [l]$ 的阶数来保证计算精度。基于多项式形式的曲率函数（1）、常应变函数（2）以及位移函数（15），可以利用Hamilton原理进一步推导得到大挠性欧拉‑伯努利悬臂梁的非线性动力学方程组^［

10］，进而对其大挠性低频动力学响应进行求解。在后续研究中，将基于本文的研究成果进一步预测欧拉‑伯努利悬臂梁低频挠性振动。

4 结论

采用变分法求解欧拉‑伯努利悬臂梁的平面超大挠性变形问题，常规思路是假设位移函数 $U [l]$ 及 $W [l]$ ，并通过几何方程求得应变函数，代入位移变分方程（11），可得到一组非线性方程组。对非线性方程组进行求解，将得到悬臂梁在外力作用下的变形。但是，直接假设位移函数 $U [l]$ 及 $W [l]$ 时，基于有限的未知参数将无法正确表征位移函数 $U [l]$ 及 $W [l]$ 之间复杂的耦合关系，进而导致计算结果不准确。本文所提出的变分求解法，其核心是基于随体坐标系，通过假设多项式形式的曲率函数及常数中心线应变，推导得到相互耦合的位移函数 $U [l]$ 及 $W [l]$ 。本文将变分法的计算结果与有限元仿真结果进行对比，计算结果显示，基于所假设的多项式曲率函数、常数中心线应变以及所推导的耦合位移函数，本文提出的变分法采用6个未知参数即能精确预测悬臂梁在末端载荷作用下的超大变形及悬臂梁末端约束反力，既有一定理论前沿性也有较好的工程应用前景。

参考文献

PAI P F. Highly flexible structures： Modeling， computation， and experimentation［M］. ［S.l.］：American Institute of Aeronautics and Astronautics（AIAA）， 2007. [百度学术]

PAI P F. Three kinematic representations for modeling of highly flexible beams and their applications［J］. International Journal of Solids and Structures， 2011， 48： 2764-2777. [百度学术]

ZHONG H， GAO M. Quadrature element analysis of planar frameworks［J］. Archive of Applied Mechanics， 2010， 80（12）： 1391-1405. [百度学术]

张琼，杜永峰，朱前坤. 基于微分求积法的 Euler-Bernoulli 梁的大变形力学行为研究［J］. 工程力学， 2014（S1）： 1-5. [百度学术]

ZHANG Qiong， DU Yongfeng， ZHU Qiankun. Study on large deformation mechanical behavior of Euler-Bernoulli using DQM［J］. Engineering Mechanics， 2014（S1）： 1-5. [百度学术]

黄正. 基于等几何配点法的几何精确Euler-Bernoulli梁几何非线性分析［D］. 武汉：华中科技大学， 2017. [百度学术]

HUANG Zheng. Geometrically nonlinear analysis of geometrically exact Euler-Bernoulli beams based on isogeometric collocation method［D］. Wuhan： Huazhong University of Science &Technology， 2017. [百度学术]

古雅琪，王海龙，杨怀宇.一种大变形几何非线性Euler-Bernoulli梁单元［J］.工程力学， 2013， 30（6）： 11-15. [百度学术]

GU Yaqi， WANG Hailong， YANG Huaiyu. A large deformation geometric nonlinear Euler-Bernoulli beam element［J］. Engineering Mechanics， 2013， 30（6）： 11-15. [百度学术]

徐芝纶. 弹性力学简明教程［M］. 北京：高等教育出版社， 2012. [百度学术]

XU Zhilun. Elastic mechanics brief tutorial［M］. Beijing： Higher Education Press， 2012. [百度学术]

HYER M W. Calculations of the room-temperature shapes of unsymmetric laminates［J］. Composite Materials， 1981， 15： 296-310. [百度学术]

WU Z， LI H， CHEN Y. An improved model for unsymmetric plates［J］. Composite Structures， 2020， 252： 1-14. [百度学术]

WU Z， LI H， FRISWELL M I. Advanced nonlinear dynamic modelling of bi-stable composite plates［J］. Composite Structures， 2018， 201： 582-596. [百度学术]

李世荣，孙云，刘平. 关于 Euler-Bernoulli 梁几何非线性方程的讨论［J］.力学与实践， 2013， 35： 77-80. [百度学术]

LI Shirong， SUN Yun， LIU Ping. Discussion on geometrically nonliear equation of Euler-Bernoulli beam［J］. Mechanics in Engineering， 2013， 35： 77-80. [百度学术]

吴坛辉，洪嘉振，刘涛永.非线性几何精确梁理论研究综述［J］.中国科技论文， 2013， 8（11）： 1126-1130. [百度学术]

WU Tanhui， HONG Jiazhen， LIU Taoyong. Advances of geometrically exact 3D beam theory［J］. China Sciencepaper， 2013， 8（11）： 1126-1130. [百度学术]

请使用 Firefox、Chrome、IE10、IE11、360极速模式、搜狗极速模式、QQ极速模式等浏览器，其他浏览器不建议使用!