目录
杆件是最常用的承力构件,特点是:
- 连接杆件的两端一般都是铰接接头
- 主要承受沿轴线的轴向力(两个连接的构件在铰接结构处可以转动,不传递和承受弯矩)
1 问题描述
下图是一个左端固定的拉杆,其右端承受一外力 P P P。拉杆的长度为 l l l,横截面积为 A A A,弹性模量为 E E E。
1.1 1D问题的基本变量
由于该问题是沿X轴方向的一维问题,因此只有沿X轴方向的基本变量,即
- 定义沿X轴方向的移动为位移: u ( x ) u(x) u(x)
- 定义沿X轴方向的相对伸长量(或缩短)量为应变: ε x ( x ) \varepsilon_x (x) εx(x)
- 定义沿X轴方向的单元横截面上的受力为应力: σ x ( x ) \sigma_x (x) σx(x)
1.2 1D问题的基本方程
(1)取杆件的任意一个截面,可得到平衡方程(无体力)为
【无体力】指不考虑体积力(也称为体力或体载荷)的影响,例如重力、惯性力等
σ x = c 1 或 d σ x d x = 0 \begin{equation} \begin{aligned} \sigma_x = c_1 \ 或 \ \cfrac{\mathrm{d} \sigma_x}{\mathrm{d} x} = 0 \end{aligned} \end{equation} σx=c1 或 dxdσx=0
式中: c 1 c_1 c1为待定的常数。
(2)取出杆件 x x x位置处的一段长度 d x \mathrm{d} x dx,设其伸长量为 d u \mathrm{d} u du,则其相对伸长量为
ε x = d u d x \begin{equation} \begin{aligned} \varepsilon_x = \cfrac{\mathrm{d} u}{\mathrm{d} x} \end{aligned} \end{equation} εx=dxdu
称为几何方程。
(3)由该材料的拉伸试验,可由虎克定律得到
ε x = σ x E \begin{equation} \begin{aligned} \varepsilon_x = \cfrac{\sigma_x}{E} \end{aligned} \end{equation} εx=Eσx
(4)边界条件(BC, boundary condition)
位移边界条件 B C ( u ) u ( x ) ∣ x = 0 = 0 力边界条件 B C ( p ) σ x ( x ) ∣ x = l = P A = P ˉ x \begin{equation} \begin{aligned} &位移边界条件 \mathrm{BC}(u) \quad u(x)|_{x=0} = 0 \\ &力边界条件 \mathrm{BC}(p) \quad \sigma_x(x)|_{x=l} = \cfrac{P}{A} = \bar{P}_x \end{aligned} \end{equation} 位移边界条件BC(u)u(x)∣x=0=0力边界条件BC(p)σx(x)∣x=l=AP=Pˉx
对于上面的力边界条件,只能作为一种近似,因在 x = l x=l x=l端面处, σ x ( x ) \sigma_x(x) σx(x)不应是均匀分布的。
由圣维南原理(Saint-Venant’s principle),在远离 x = l x=l x=l的截面,力的边界条件才能较好地满足。
1.3 求解思路
(1) 直接求解方法
对于该简单问题,可以由3个方程来直接求解3个变量
(2) 基于试函数的间接方法
可先选取一个变量(如位移)作为最基本的待求变量,将其他变量都有它来表达,并采用间接的近似求解方法,具体的做法为:
- 先对待求的位移变量假设一种事先满足位移边界条件的可能解(其中有一些待定的系数),称为试函数(trail function)
- 让该受力系统的【势能取最小值】 / 【内部变形虚功=外部施加力的虚功】 来最后确定处可能解(试函数)中的那些待定系数(unknown constant)
2 1D问题的直接求解
σ x ( x ) = c 1 或 d σ x ( x ) d x = 0 ε x ( x ) = d u d x ε x ( x ) = σ x E } \begin{equation} \begin{aligned} \left.\begin{matrix} \sigma_x(x) = c_1 \ 或 \ \cfrac{\mathrm{d} \sigma_x(x)}{\mathrm{d} x} = 0 \\ \varepsilon_x(x) = \cfrac{\mathrm{d} u}{\mathrm{d} x} \\ \varepsilon_x(x) = \cfrac{\sigma_x}{E} \end{matrix}\right\} \end{aligned} \end{equation} σx(x)=c1 或 dxdσx(x)=0εx(x)=dxduεx(x)=Eσx⎭ ⎬ ⎫
求解上述方程组,得到
σ x ( x ) = c 1 ε x ( x ) = c 1 E u ( x ) = c 1 E x + c 2 } \begin{equation} \begin{aligned} \left.\begin{matrix} \sigma_x(x) = c_1 \\ \varepsilon_x(x) = \cfrac{c_1}{E} \\ u(x) = \cfrac{c_1}{E}x + c_2 \end{matrix}\right\} \end{aligned} \end{equation} σx(x)=c1εx(x)=Ec1u(x)=Ec1x+c2⎭ ⎬ ⎫
式中: c 1 , c 2 c_1, c_2 c1,c2为待定的常数。由边界条件式(4)可求得 c 2 = 0 c_2 = 0 c2=0, c 1 = P A c_1 = \cfrac{P}{A} c1=AP,故有
σ x ( x ) = P A ε x ( x ) = P E A u ( x ) = P E A x } \begin{equation} \begin{aligned} \left.\begin{matrix} \sigma_x(x) = \cfrac{P}{A} \\ \varepsilon_x(x) = \cfrac{P}{EA} \\ u(x) = \cfrac{P}{EA}x \end{matrix}\right\} \end{aligned} \end{equation} σx(x)=APεx(x)=EAPu(x)=EAPx⎭ ⎬ ⎫
3 1D问题的虚功原理求解
3.1 虚功原理
如图为一个平衡力系,由于该系统处于平衡状态,故有
P A P B = l B l A \begin{equation} \begin{aligned} \cfrac{P_A}{P_B} = \cfrac{l_B}{l_A} \end{aligned} \end{equation} PBPA=lAlB
假想在该平衡力系上作用有微小的扰动(不影响原平衡条件),且外力所作用的位置产生了微小的位移变化,即 Δ A \Delta_A ΔA, Δ B \Delta_B ΔB
该假想的位移如果不影响原平衡条件,应满足
Δ B Δ A = l B l A \begin{equation} \begin{aligned} \cfrac{\Delta_B}{\Delta_A} = \cfrac{l_B}{l_A} \end{aligned} \end{equation} ΔAΔB=lAlB
这就是任意扰动的位移应满足的条件,称为许可位移条件。将满足许可位移条件的、任意微小的假想位移称为虚位移(virtual displacement)
联立式(8)和(9),有
P A Δ A − P B Δ B = 0 \begin{equation} \begin{aligned} P_A \Delta_A - P_B \Delta_B = 0 \end{aligned} \end{equation} PAΔA−PBΔB=0
P A Δ A P_A \Delta_A PAΔA为 A A A点的虚功(virtual work),即外力 P A P_A PA在 A A A点虚位移 Δ A \Delta_A ΔA上所做的功
− P B Δ B - P_B \Delta_B −PBΔB为 B B B点的虚功,即外力 P B P_B PB在 B B B点虚位移 Δ B \Delta_B ΔB上所做的功,负号表示力与位移的方向相反
式(10)即为基于虚位移的虚功原理(principle of virtual work)
对于一个处于平衡状态的系统,作用于系统上的所有外力在满足许可位移条件的虚位移上所做的虚功总和恒为0
对于弹性力学中含有变形体的虚功原理,此时的虚功应包括外力虚功 δ W \delta W δW和内力虚功 − δ U -\delta U −δU, δ U \delta U δU叫做虚应变能(virtual strain energy)。由于弹性体在变形过程中,内力是抵抗变形所产生的,其方向总是与变形的方向相反,所以内力虚功取负。
由于虚功总和为零,则有
δ W − δ U = 0 → δ W = δ U \begin{equation} \begin{aligned} \delta W - \delta U = 0 \to \delta W = \delta U \end{aligned} \end{equation} δW−δU=0→δW=δU
弹性力学中的虚功原理可表述为
在外力作用下处于平衡状态的变形体,当给物体以微小虚位移时,外力所做的总虚功等于物体的总虚应变能(即应力在由虚位移所产生虚应变上所做的功)
注意:此处的虚位移指仅满足位移边界条件 B C ( u ) BC(u) BC(u)的许可位移
3.2 虚功原理求解
对于上文中的一端固定的拉杆问题,设满足位移边界条件式(4)的位移场为
u ( x ) = c x \begin{equation} \begin{aligned} u(x) = cx \end{aligned} \end{equation} u(x)=cx
式中: c c c为待定系数,代表杆件在单位长度上的位移变化。 u ( x ) = c x u(x) = cx u(x)=cx为待定函数,也称为试函数
基于试函数,该系统的应变、虚位移以及虚应变为
ε ( x ) = c δ u ( x ) = δ c ⋅ x δ ε ( x ) = δ c } \begin{equation} \begin{aligned} \left.\begin{matrix} \varepsilon(x) = c \\ \delta u(x) = \delta c \cdot x \\ \delta \varepsilon(x) = \delta c \end{matrix}\right\} \end{aligned} \end{equation} ε(x)=cδu(x)=δc⋅xδε(x)=δc⎭ ⎬ ⎫
式中: δ c \delta c δc为待定系数的增量,即待定系数 c c c发生微小变化。
虚位移 δ u ( x ) = δ c ⋅ x \delta u(x) = \delta c \cdot x δu(x)=δc⋅x表示在试函数中引入了一个微小的变化 δ c \delta c δc后,位移场的变化
虚应变 δ ε ( x ) = δ c \delta \varepsilon(x) = \delta c δε(x)=δc表示由于 δ c \delta c δc引起的应变变化
计算得到拉杆结构的虚应变能为
δ U = ∫ Ω σ x δ ε x d Ω = ∫ 0 l ∫ A E ⋅ ε x δ ε x d A d x = E ⋅ c ⋅ δ c ⋅ A ⋅ l \begin{equation} \begin{aligned} \delta U &= \int_{\Omega} \sigma_x \delta \varepsilon_x \mathrm{d \Omega} \\ &= \int_{0}^{l} \int_{A} E \cdot \varepsilon_x \delta \varepsilon_x \mathrm{d} A \mathrm{d} x \\ &= E \cdot c \cdot \delta c \cdot A \cdot l \end{aligned} \end{equation} δU=∫ΩσxδεxdΩ=∫0l∫AE⋅εxδεxdAdx=E⋅c⋅δc⋅A⋅l
外力虚功为
δ W = P ⋅ δ u ( x = l ) = P ⋅ δ c ⋅ l \begin{equation} \begin{aligned} \delta W &= P \cdot \delta u(x=l) \\ &= P \cdot \delta c \cdot l \end{aligned} \end{equation} δW=P⋅δu(x=l)=P⋅δc⋅l
由虚功原理,有
E ⋅ c ⋅ δ c ⋅ A ⋅ l = P ⋅ δ c ⋅ l c = P E A \begin{equation} \begin{aligned} E \cdot c \cdot \delta c \cdot A \cdot l &= P \cdot \delta c \cdot l \\ c &= \cfrac{P}{EA} \end{aligned} \end{equation} E⋅c⋅δc⋅A⋅lc=P⋅δc⋅l=EAP
4 1D问题的最小势能原理求解
考虑在变形过程中,材料从无应变状态逐渐增加应变到某一值 ε \varepsilon ε,在这个过程中,外力做功,并且这部分功被存储在材料中作为应变能
在应变从0增加到 ε \varepsilon ε的过程中,瞬时应力 σ \sigma σ随应变线性增加。对于一个微小的应变增加 d ε \mathrm{d} \varepsilon dε,瞬时应力 σ \sigma σ也会变化。因此,微小增量下的做功(即能量增量)为
d U = σ d ε \begin{equation} \begin{aligned} \mathrm{d} U = \sigma \mathrm{d} \varepsilon \end{aligned} \end{equation} dU=σdε
为了得到从无应变状态到最终应变 ε \varepsilon ε过程中存储的总能量,需要对上述能量增量进行积分
U = ∫ 0 ε σ d ε = ∫ 0 ε E ε d ε = E ∫ 0 ε ε d ε = E ⋅ 1 2 ε 2 = 1 2 σ ⋅ ε \begin{equation} \begin{aligned} U &= \int_0^{\varepsilon} \sigma \mathrm{d} \varepsilon \\ &= \int_0^{\varepsilon} E \varepsilon \mathrm{d} \varepsilon \\ &= E \int_0^{\varepsilon} \varepsilon \mathrm{d} \varepsilon \\ &= E \cdot \cfrac{1}{2} \varepsilon^2 \\ &= \cfrac{1}{2} \sigma \cdot \varepsilon \end{aligned} \end{equation} U=∫0εσdε=∫0εEεdε=E∫0εεdε=E⋅21ε2=21σ⋅ε
4.1 最小势能原理求解
设有满足位移边界条件 B C ( u ) \mathrm{BC} (u) BC(u)的许可位移场 u ( x ) u(x) u(x),计算该系统的势能(potential energy)为
∏ ( u ) = U − W \begin{equation} \begin{aligned} \prod(u) = U - W \end{aligned} \end{equation} ∏(u)=U−W
式中: U U U为应变能; W W W为外力功。对于上文的一维拉杆结构,有
U = 1 2 ∫ Ω σ x ( u ( x ) ) ⋅ ε x ( u ( x ) ) d Ω W = P ⋅ u ( x = l ) \begin{equation} \begin{aligned} U &= \cfrac{1}{2} \int_{\Omega} \sigma_x(u(x)) \cdot \varepsilon_x(u(x)) \mathrm{d \Omega} \\ W &= P \cdot u(x=l) \end{aligned} \end{equation} UW=21∫Ωσx(u(x))⋅εx(u(x))dΩ=P⋅u(x=l)
对于包含有待定系数的试函数 u ( x ) u(x) u(x)而言,真实的位移函数 u ( x ) u(x) u(x)应使得该系统的势能取极小值(能量处于最低状态),即
m i n u ( x ) ∈ B C ( u ) [ ∏ ( u ) = U − W ] \begin{equation} \begin{aligned} \mathop{min}\limits_{u(x) \in \mathrm{BC}(u)}[\prod(u) = U - W] \end{aligned} \end{equation} u(x)∈BC(u)min[∏(u)=U−W]
同样取满足位移边界条件的位移场为
u ( x ) = c x \begin{equation} \begin{aligned} u(x) = cx \end{aligned} \end{equation} u(x)=cx
式中: c c c为待定系数。计算得到应力和应变为
ε x ( x ) = d u d x = c σ x ( x ) = E ⋅ ε x ( x ) = E ⋅ c \begin{equation} \begin{aligned} \varepsilon_x(x) &= \cfrac{\mathrm{d} u}{\mathrm{d} x} = c \\ \sigma_x(x) &= E \cdot \varepsilon_x(x) = E \cdot c \end{aligned} \end{equation} εx(x)σx(x)=dxdu=c=E⋅εx(x)=E⋅c
计算该系统的势能为
∏ ( u ) = U − W = 1 2 E ⋅ c 2 ⋅ A ⋅ l − P ⋅ c ⋅ l \begin{equation} \begin{aligned} \prod(u) = U - W = \cfrac{1}{2} E \cdot c^2 \cdot A \cdot l - P \cdot c \cdot l \end{aligned} \end{equation} ∏(u)=U−W=21E⋅c2⋅A⋅l−P⋅c⋅l
求极值有
∂ ∏ ( u ) ∂ c = 0 E ⋅ c ⋅ A ⋅ l − P ⋅ l = 0 c = P E A \begin{equation} \begin{aligned} \cfrac{\partial \prod(u)}{\partial c} &= 0 \\ E \cdot c \cdot A \cdot l - P \cdot l &= 0 \\ c &= \cfrac{P}{EA} \end{aligned} \end{equation} ∂c∂∏(u)E⋅c⋅A⋅l−P⋅lc=0=0=EAP
5 总结
从上面的计算可以看出,基于试函数的方法,包括虚功原理和最小势能原理(priciple of minimum potential energy),仅计算系统的能量,实际上是计算积分(式(14)和式(20)),然后转化为求解线性方程,不需求解微分方程,这样就极大地降低了求解难度。
试函数方法的关键是如何构造出适用于所求问题的位移试函数,并且该构造方法还应具有规范性以及标准化,基于“单元”的构造方法就可以完全满足这些要求。
参考
- 曾攀. 有限元分析基础教程[M]. 北京: 清华大学出版社, 2008.