薄板弯曲的变分原理及有限元素法

第三章 薄板弯曲的变分原理及有限元素法

3.1 基本问题

基本认识：板作为承力的结构元件，主要通过弯曲起作用。如果垂直于板面的挠度与板的厚度相

比很小的话1<

w

，则由弯曲而引起的板中面的拉伸作用就可以忽略不及，这是所谓的小挠度问题，一般认为4.0

以下。 反之，H

w

越大，弯曲引起的中面拉伸的影响越来越大，就不能忽略不计，导致所谓大挠度

问题。

除板的弯曲变形之外，还伴随有剪切变形，剪切作用的影响一方面取决于材料的剪切模量，另一方面取决于厚度/跨度（l H ）之比，即横向剪切随l H 的增大而增大。通常把不考虑剪切作用（横向剪应变无穷大）的板理论叫做薄板理论，把考虑剪切作用的板理论叫做厚板理论。本章仅考虑小挠度薄板问题。

基本假设：取板的中面为xy 平面，取z 轴与y x ,轴垂直，设板的厚度为h ，可以是()y x ,的函数。 ① 变形假设：变形前垂直于中面的直线段在变形后没有伸缩，并且继续垂直变形后的中面。

由此得： ② 内力假设：板内应力的6个分量的大小不是同一量级，一般xy y x τσσ,,最大，yz xz ττ,约小

一个量级，而z σ又小一个量级；在静力学分析中，0=z σ。 控制方程（内力平衡方程及物理方程）

① 由弹性力学方法，对于均质材料构成的薄板，应力分量yz xz xy y x τττσσ,,,,可用5个内力

()()()()()y x Q y x Q y x M y x M y x M y x xy y x ,,,,,,,,,表示，即：

x x zM h 312=

σ y y zM h 312=σ xy xy

zM h 3

12=τ （矩定义为单位宽度上的矩） Note ：上述的弯距及剪力代表单位宽度上的，而不是整个板侧面的。 ② 用内力表示的平衡方程：

0=+??+??p y

Q x Q y

x ()y x p p ,= 分布的横向载荷 在薄板理论中，内力y x Q Q ,不产生应变，因而也不做功，可在以后的分析中不计算它们，在上式中消去y x Q Q , 即得：

③ 几何关系： ④ 物理关系：（各向同性体）

点应力应变关系：??

???????????????????

?--=??

?

??

?????xy y x xy y x v v v v E εεετσσ2100010112

内力与应变关系：

注意： ()v E G +=12 ()

2

3

112v Eh D -=

⑤ 单位面积上的应变能及余应变能（密度）

应变能密度（曲率作为自变量）变分：

[]xy xy y y x x k M k M k M U δδδδ2++= ??

?

??????=

x w x k x (单位长度上转角的变化) ∴ x x k U

M ??=

y y k U M ??= xy

xy k U M ??=21（这也是一种物理关系） 代入关于内力矩的物理关系，有：

注意：上式中都是关于曲率的二次项，而且从物理上对于任意的曲率U >0，故U 称为正定的二次齐次函数。

余应变能密度：（可看作是内力矩的函数） 当板的变形由一种状态变到相邻的另一种状态时， V 的变分为：

在以前的研究中，我们把内力表示成变形的函数，从而构造和研究关于变形的泛函；在这里，我们把体系中的两类量都看成独立可变的，故有上式V δ，（上式中由于V 的表达式，变分的结果又可只认为只有内力矩变化）。这是一种认识观点，对于后面理解广义变分原理有利的。当然，还应当注意这两类变量之间存在着物理关系的约束。 由V δ的表达式可知： x x M V

k ??=

y y M V k ??= xy

xy M V k ??=2（研究其中一项时，让其它两项的变分为0） 代入内力矩关系，可知V 是xy y x M M M ,,的正定二次齐次函数。 ⑥ 用挠度表示的平衡关系：

把内力矩的物理关系及曲率的几何关系代入原平衡方程，即得：

p w D p y w y x w x w D =??

=???

? ????+???+??444224442 （双调和方程）

坐标旋转引起的变换

在研究板问题时，经常用到不同坐标系表示的包括法向导数等，因此在不同坐标系下，板弯

曲的基本量之间有什么联系是我们经常要遇到的计算。 取两个不同的坐标系，如右图 ① 坐标变换关系： ② 函数的方向导数：

对于一个函数()()()[]ηξηξ,,,,y x F y x F = 由求导的链式规则：

同理：θθθθη22222

222

2cos sin cos 2sin y

F y x F x F F ??+???-??=?? ③ 只要将F 换成w 即得曲率在不同坐标系下的转轴公式。 ④ 弯、扭矩的转轴规律

应注意到y xy x M M M ,,是单位宽度上的矩，推导时要依据平衡关系。 由此得：

剪力的转轴规律符合矢量投影规律，即： 典型的边界条件 （如右图）

Ω：板在xy 平面上所占区域；

c ：板的边界； n ：边界外法线；

s ：边界切线方向；s n ,符合右手定则

典型边界条件有三种：

① 固支边：如在1c 上，w w =，

n n

w

ψ=?? 这里，n w ψ,

为1c 上已知的关于弧长s 的函数。

② 简支边：如在2c 的部分边界上简支，则有：

n M w ,为2c 上已知的关于弧长的函数。

③ 自由边：在自由边上(指无位移限制)已知作用在边界上的力（即所谓的自然边界条件）。 从内力和内力矩角度看，边界上能反映出来的有n ns n Q M M ,,三个，但不能取：n n M M =，

ns ns M M =，n n Q Q =

因从做功角度讲，ns M 和n Q 并不完全独立，分析如下：

给自由边界3c 上的挠度有一变分w δ，则n ns Q M ，在w δ上所作的功为： 可见，ns M 相当于

s

M ns

??线分布载荷以及作用在自由边两端的集中载荷ns M ±。由于在自由边

两端总是有支持的，所以在该两端上的ns M 对板变形不产生影响，所以分析时略去不计。 ∴wds Q s

M U c n ns δδ???

?

??+??=

3 由上分析，

n ns

Q s

M +??是与w δ相应的广义力，故自由边的边界条件应为： 在3c 上，n n M M =，

n n ns

q Q s

M =+??（与点应力的力边界条件相似） ()s q n 是已知作用在3c 上的线布载荷。

若在自由边界的某一点0s s =上有一集中载荷p ，那么有：

()()p s M s M ns ns =--+0000 p 表示单位宽度上的力。

若没有集中载荷，ns M 应是s 的连续函数，ns M 的跳跃量相应于集中载荷。 ﹡自由边界上有尖角的情况： -n 为A 点前面一段终点的法向量；

+n 为A 点后面一段起点的法向量。

命两个法向量与x 轴的夹角为+-θθ,。 在A 点前后的两个扭矩：

由于要求ns M 连续，即+

-=ns

ns M M 由此得：

当A 点尖角的两条边平行于y x ,轴，

0,90==-+θθ 故得： 0=xy M 3.2

最小势能原理

考虑如图的板受载系统：

1c 上：n n

w

w w ψ=??=,

2c 上：n n M M w w ==,

3c 上：q Q s

M M M n ns

n n =+??=, 整个系统的势能包括两部分：

C 2 3

① 板的应变能??Ω

=

∏Udxdy 2 U 为应变能密度，是曲率的二次函数；

② 外载荷（包括边界力）的势能： ds n

w

M wds q pwdxdy c c c n

??+--

=∏????

Ω

+3

3

2'

(一次泛函，固支边界位移变分为0） 令w 是问题的精确解，k w 是可能的挠度（在1c 上满足n n

w

w w ψ=??=,；在2c 上满足w w =） 最小势能原理指出：

与精确解w 相应的总势能()w ∏达到最小值，即小于任何其它可能位移k w 相应的总势能。 证明的过程与梁的步骤全同。 令：w w w w w w k k ?+=-=?,

w ?满足齐次的边界条件：?????

=?=???=?0

,

0,0,

21

w c n

w

w c 上上 与k w 相应的总势能

可以证明：()0,''=?∏w w ，且若w ?不是刚体运动，则()02>?∏w 。 因此有， ()()w w k ∏>∏

将最小势能原理写成变分形式：0=∏δ (留给同学自己完成)

()022322222=??--??

?

??-+??+???? ??-??-???-??-?????Ωds n w M M wds q Q s M wdxdy p y M y x M x M c n n c n ns y

xy x δδδ（这里是Galerken 虚功的形式）

?平衡方程： 022

2222=+??+???+??p y

M y x M x M y xy x

3c 上

q Q s

M n ns

=+?? 2c 上 n n M M = 3.3

最小余能原理

① 概说

余能原理是一个数学原理，并不对应物理上的保守势场中的能量守恒原理。它的着眼点是固定位移，看可能的力系统变化（这些力系统一定满足平衡关系和力边界条件，这是一种狭义理解

方式）。那么从此观点上看，显然只有满足位移协调关系的那一个力系才是精确解。

余能原理：在满足平衡关系的所有力系中，只有满足位移协调关系的那组力系，使系统余能取极值。如果这个系统是稳定的，则取最小值。 ② 数学证明：

考虑与上小节相同的一块薄板弯曲问题，命y x xy y x Q Q M M M w ,,,,,是问题的精确解；再取

s

y

s x s xy s y s x Q Q M M M ,,,,为一组可能的内力。 按要求可能的力系满足：

a ． 平衡条件：??

???????

??=+??+??=-??+??=-??+??000p y Q x

Q Q y M x M Q y M x M s y

s x s

y s

y s xy s x s xy

s x b ． 力边界条件：在2c 上：n s

n M M =

在3c 上：q Q s

M M M

s

n s

ns n s ns

=+??=,

首先，要做的工作是分别给出精确解和可能的内力的余能大小。系统余能包括两部分，一部分是域内的余应变能2

Γ；另一部分是已知位移的部分边界上的余能1

Γ。 xy xy y y x x k M k M k M V 2++= （在线性应力应变关系下） 整个系统的总余能： ds M ds w Q s M Vdxdy n c n c c n ns ψ????

+??

?

??+??-=

Γ+Γ=ΓΩ

+1

212

1

(边界上已知力变分为0) 固定了某种位移形态，则自变函数即为{}

,,,,,M =y x xy y x Q Q M M M 与精确解相应的余能记为：()()()M M M 321,,ΓΓΓ；与可能内力相应的余能为：()()()

s s s M M M 321,,ΓΓΓ。

取：

x s

x x Q Q Q -=? y s

y y Q Q Q -=? （内力，而不是外载荷）

由于两套力系都满足平衡关系及力的边界条件，这些关系又都是线性的，所以，可能内力增量满足下列齐次方程和齐次边界条件： 在2c 上：0=?n M 在3c 上：0,

0=?+???=?n ns

n Q s

M M

上两套式子表明：与内力增量相应的外载荷为零。将系统余能细致给出表达式，有：

由功的互等定理，知上式等于零。 ()()()[]

dxdy M M d M d M d M y x y x ??Ω

??+?+?=

?Γ122

222112221 （大于零） 若0≠?M 则有：()

()M M s Γ≥Γ 此即最小余能原理。

余能原理的变分表达即：0=Γδ

由余能原理要求，可能内力必须满足平衡条件，所以y x M M δδ,及xy M δ不能独立变分（即任意变化），即是说，这种变分是条件变分（原理）。其二，可以想象，因为可能力系的变化要求满足平衡关系，所以变分的结果肯定是用力表达的变形协调关系。 3.4.2 9个位移参数的三角形部分协调元 1．三角形划分

每个节点的位移参数： y

w

x w w y x ??=??=ψψ,, 节点的位移列阵

T T T T

q q q ][32

1

=

2．面积坐标系 ① 面积坐标的定义：

i i A

A

ξ=（i = 1,2,3）

2

2

33

111223

3

111111x y x y x y x y x y x y ξ=

，其他两个类推。 ② 插值函数阶次与系数的基本概念

? 完备的单（多）项式； ? Pascal 图；

? 冗余坐标标架下的多项式 3．位移插值函数：

一般选择完备的单项式基，不完备的单项式基一般不能保证收敛性（特别是少常数项和1次项）。这里，选择全三次多项式（三次项是完备的，Note : 用节点参数计算三次多项式系数时少一个）

各基函数的形状：

i

j

ξ2

ξ1 ξ3

k

32110ξξξα

，称为一种内部自由度。 决定10α的办法：（Zenkiewicz ）

)(9876542110ααααααα+++++=, (可满足收敛性要求) 非物理参数的挠度插值表达式：α?T w =,

物理参数变换：

在三角形的三个顶点，挠度及其导数应等于设定的值（待求参数），代入各节点的挠度及其导数，由此可得9个方程：

T yk xk k yj xj j yi xi i e w w w ][ψψψψψψ=q （节点位移参数列阵）

协调性分析：

○1 单元边上的挠度函数是一条三次曲线，需四个参数，现有6个参数（每个节点上的两个转角参数并不独立，独立的参数只有4个），完全可以由该边上选择相邻单元的共同一组位移参数决定，而无需第3个节点上的位移参数。故若两个三角形元素公用一条边，则在这条边上的挠度是连续的。

○

2 边上的法向斜率：二次曲线，需三个参数，但仅可知道两节点的n

w

??，故不能确定，要借助第三个节点的n

w

??。所以，单元的连接边上有硬折。 4．刚阵计算

○

1 计算曲率与扭率 公式：

T

w w w w

w y x w y w x w ]2[2

123123222322222

1

2222

22ξξξξξξξξξ???????????????=????

?

????

????????????????Ω（具体计算Ω！）

再对α?

T

w =求偏导，有：

② 单元弯曲应变能及刚阵计算：

????

??????-=γγγ

1000101d D , )

1(1223γ-=Eh d ③ 等效载荷列阵计算

ξi

ξi 2ξj

ξi ξj ξk

x y x 2 xy y 2 x 3 x 2y x y 2

y 3 ??-

=e

dxdy pw e Ω

∏1

（p 是分布载荷）

④ 单位刚阵/载荷列阵集合

⑤ 计算刚度方程（先置入边界条件）

K q = F

3.4.3 12个位移参数的矩形部分协调元 1．矩形域中的无量纲坐标

ηξb y a ==, 坐标变换方法2：)1(2

),1(2

,12,12ηξηξ+=

+=

-=

-=b

y a

x b

y a

x 现用1f 映射： 定义导数算子关系：

η

ξ??=????=??b y a x 1,1 2．12个位移参数的矩形部分协调元

○

1 单元划分：若干个有限的矩形单元（如上图）每个角点3个位移参数 y x θθ,定义与常规转角的定义不同，可看作广义转角位移，主要为计算方便。

○

2 位移选择 3

1231131029283726524321ξηαηξαηαξηαηξαξαηαξηαξαηαξαα+++++++++++=w 注意：这个插值式包含了完全的三次多项式，而四次项只适当地选用了两项。选择这样的多项式

是关于坐标对称的，因此，能够保证关于对称轴LL ，MM 对称或反对称的四种变形都能用3个自由度来描述(关于LL , MM 都是对称的；关于LL 对称，关于MM 反对称；关于LL 反对称，关于MM 对称；关于都是反对称的四种变形)。 协调性分析：

* 在每条边上，w 是三次函数（x 或y ）， 恰好4个参数（因为是直角），w 连续。 * 在每条边上，

η??w 是x 和y 的二次多项式（如1=η的边）。在两端点上只有两个参数与η

?w 有关，不能确定，需要其它边上的参数，故而

η

??w

不协调。 ○

3 挠度函数与节点自由度的关系: P 278 (17.4式) L

1

2

Homework : 用矩阵形式推导出矩形板单元的刚度矩阵计算式及等效的载荷列阵。

○

5 2

∏ P276, 16.4式，转角代数式P279,17.6式 ○

6 刚阵 P280 3、16个位移参数的矩形过分协调元素分析

仍以12个参数的单元矩形ABCD 为例，取局部坐标x , y 和局部无量纲坐标ξ ，η。对于每一个节点赋予四个位移参数：

这样，每一个矩形共有16个位移参数，挠度w 用一个包含16个未知数的多项式进行插值：

∑∑===303

m n n m mn w ηξα (对两个坐标变量都是完全的三次多项式)

上述插值方法，能够保证在每条公用边上，w 及

n

w

??的连续性，因为在每条边上（如ξ＝0），w 是三次曲线，完全可以由两点的w 和θ 所决定（如A 、D 两点的w 和θy ）；同理，关于转角也是三次函数（如ξ＝0边上的

ξ??w ），同样可以由A 、D 两点的ξ??w 和???

?

??????ξηw 所决定。故这种单元保证了挠度和转角的连续性。但过分协调了，因为此位移函数还导致了公共节点上ξ

η???w

2的连续

性，从单元的能量积分角度而言，这种连续性是不必要的。 3.4.4 建立协调元素的方法

1. 二次分片插入法

2. 条件极值法

仍以9 在这个由9法向斜率并不连续，且有一

个任意的自由度。。关于每条边上的斜联立这个条件求解系统能量的极值，就成为一个有条件的函数极值问题。可以用Lagrange 乘子

法进行计算，但一般说，这种计算有一个大的问题：求解的方程组不再保持对称性，使求解复杂起来。

○

4曲率与节点自由度的关系 P 278-279, 17.5式