计算力学教学课件

绪论计算力学简介.........................................................3

第一章有限元法简介.......................................................3

§1-1引言..............................................................3

§1-2简例..............................................................4

§1-3有限单元法计算过程概述............................................8

第二章平面问题的有限元法.................................................9

§2-1三角形常应变单元..................................................9

§2-2形函数的性质.....................................................11

§2-3刚度矩阵.........................................................12

§2-4等效结点力.......................................................16

§2-5热应力计算.......................................................19

§2-6收敛准则.........................................................20

§2-7实施步骤.........................................................22

§2-8矩形单元.........................................................25

§2-9三角形单元面积坐标...............................................28

§2-10习题讲解........................................................30

第三章空间问题..........................................................32

§3-1四面体常应变单元.................................................32

§3-2刚度阵、等效结点力...............................................33

第四章轴对称问题........................................................33

§4-1三角形截面环单元.................................................33

§4-2单刚阵...........................................................35

§4-3等效结点力.......................................................36

§4-4精确刚度阵计算...................................................38

第五章等参数单元........................................................39

§5-1平面等参元.......................................................39

§5-2高斯求积法.......................................................46

§5-3等参元形态.......................................................48

§5-4轴对称等参元.....................................................50

§5-5空间等参元.......................................................51

第六章杆件系统有限元法................................................53

§6-1等截面梁单元....................................................53

§6-2等效结点力.....................................................55

§6-3坐标变换.......................................................58

§6-4板与杆的组合结构...............................................61

§6-5空间杆件系统...................................................63

第七章板的弯曲..........................................................65

§7-1引言.............................................................65

§7-2矩形薄板单元...................................................66

§7-3考虑横向剪切影响的平板弯曲单元.................................71

第八章壳的弯曲..........................................................77

§8-1平板壳体单元...................................................77

第九章动力学问题有限元法................................................81

§9-1动力学方程和质量阵、阻尼阵.....................................81

§9-2无阻尼自由振动.................................................83

第十章材料非线性问题的有限元法.........................................86

§10-1非线性问题求解方法............................................86

§10-2弹塑性有限元法................................................94

§10-3弹塑性问题的求解方法..........................................97

第十一章加权残数法.....................................................98

§11-1方法概述......................................................98

§11-2分类..........................................................99

§11-3最小二乘法...................................................101

各种单元的特点..........................................................102

2

绪论计算力学简介

1、计算力学概念

以数值计算方法求解力学问题的方法。

实质：将求解复杂的微分方程组问题转化为求解代数方程组问题；

优点：①可求解解析法难以解决的问题，绘出趋近于精确解的近似解。

②可利用计算机解决十分庞大复杂的问题。

2、计算力学内容

有限差分法

有限单元法

边界学元法

加权残数法

第一章有限元法简介

§1-1引言

1、起源：50年代飞机结构矩阵分析

60年代推广到弹性平面应力问题，提出“有限单元法”。

2、方法：连续体一离散体

(1)将连续体分为有限个小单元体，彼此在结点处相连联结。

如图：

(2)选函数来近似单元位移分布规律，用虚功原理建立单元结点力与结点位移间的

关系。

3

如图:

(3)单元集合，得到以结点位移为未知量的代数方程组。

-[K]⑻=[R]

(4)求解方程组，得到结点位移，再求其他量。

3、实质：把具有无限自由度的连续体化为有限个自由度的单元集合体。

4、有限元法与经典解析法的区别

经典解析法从微元体入手，求解偏微分方程。

有限单元法从有限大小的单元体的力学特性入手，求解代数方程组。

5、有限单元法的优点

(1)概念浅，易掌握

(2)应用范围广

(3)矩阵表达式，易编程计算

§1-2简例

4

【单元分析】

对①单元

{/7=&：）/+k^vt+k^u2+k^v2

匕⑴=媚〃+成匕+成%2+%如2

叫）=燔〃+成匕+噌％+成乜

匕⑴=心：％+尺加+成%2+噌丫2

可写成

5

卜1)kK⑴

U，用）K\213%14W,

后⑴

匕⑴k⑴鼠23代24V)

>—研儿22

%⑴

kK⑴k⑴

碉32长33%34%

k⑴k⑴女⑴

公？鼠42%43K44」匕，

对②单元

.一»娟3）3

6)

6

(3

3一>u；'(3i

同理可得

烧）八2）

U?'鼠36u,

匕⑵喘）女⑵

陶长46V2

无⑵

up噌婿喈鼠56“3

女⑵

M⑵一党喈长66.匕.

【单元组装】

①单元改写为

U”飞;)靖）噌0o-U】

乂⑴掰）成成碳00匕

邛00u

>—噌胤）噌噌2

VJ端理）嵋陶00V2

0000000lly

0000000*

②单元改写为

0000000M,

0000000W

00唠出噌牖M,

“(2)

'V产K

00c吠46之

uf)00噌点唠%

“⑵

M⑵00嘤嘤K66_73.

6

X」M）偿噌噌00'%、

XK⑴心）嫩噌或00W

局）成域+益）媚；媚喏媚

Xu,+u,«2

可得2

八匕⑴+匕⑵胤）爆爆+噌媚+媚媚牖丫2

X、00小）唠唠唠“3

7s.〔匕⑵J00党点牖嘤

解释:

Zx=OX2-U^-U^=0=x2=u=+up

(l)(2)

ZY=0丫2-匕⑴一匕⑵=0=>y2=v2+V2

=>{R}=[K]{6}

【刚度系数的计算】

A/=cos0

需轴向压力=r△/=^cos6

x、y方向分量为：

I7APA

结点1,=(-y-COS^)COS^=-y-COS20

内；=(岸■cos0)sin0=cos0sin0

结点2,噂=-噌=-"cos?。

[7A

&4：1；)=一招ZI)=---jcos8sing

7

§1-3有限单元法计算过程概述

1、结构离散化

2、选位移模式

位移模式：单元位移分布假定

[门-单元内任一点位移，{5}一单元结点位移，[N]-形函数阵

如：

3、分析单元特性

（1）几何方程：{£}=[B]{b},

（2）物理方程：{<T}=[D][B]{J}e

（3）平衡方程：[R「=[k]"L

[k]-单刚阵[R--结点载荷阵

4、等效结点力

将外力移置到结点上。

5、集合单刚体，形成结构平衡方程

[K]"]=[R]

6、引入边界条件后，求解结点位移，计算单元应力

8

第二章平面问题的有限元法

§2-1三角形常应变单元

一、离散化

结点位移表示为：

｛盯=[%V,.UjVjun,vj=[邛或]7

假设代入i,/,机的结点坐标,

[v=a4+%x+。6y

Uj=ax+。2玉+戊3%匕-。4+二5巧

=><%=%+a2Xj+a3yjvy.=a4+a5x.+a(iyj位移模式

u

,n=%+a/m+a3ymvnt=%++a6ym

=>,,%，%,%,%q

将药,%,%代入位移模式,

"(也+3+a〃,+(%+㈠+3飒+a+〃x+，/)〃"』

其中A,q，e，q为:

9

%

2A=

a,j,机)

1

4=XjJ,“一xy;b；=y-y„,;q=-(x-x)

m五m

N

.-(q+6jX+qy)

/_

N

.-L

令形函数4J2(%+/?/X+c.y)贝lj：u=NM+N秒j+N“M

NA±

w-

2

同理：v=NR+AN匕+Nmvni

MON,。N0

其中［N］为形函数阵，［N］=m为二阶单位阵

oM0Nj0

三、应变

一、

£

X,〃0耳0bm0

l卜

-<£=

v

»

I：

£

z

)

+

a-v

ax

n{e}=[8]{3『

也O'

其中n[6]=困BjB,„],闻0%("M

cb

L(t.

四、应力

{(T}=<ay”=[D]{e}=回[8]{b『，[D]见书本。

巴.

i〃o

平面应力问题固=占〃10

1一Ni

O0〜

2

io

10

1一〃

平面应变问题：固=£(「〃)A10

L」(1-2〃)(1+〃)1-〃

001*

§2-2形函数的性质

形函数N,=」-(%+4x+qy)

2A

1、形函数在结点上的值

N，a，%)=亢3+%+纱，)=1

M(弓,刀)=*(《++q刀)=0

对乂，bixj

N,G“，》,“)=[(4

+3“+。,几)=0

[2A

Nj(x,,y,)=O'N〃G，X)=O

同理：＜Nj(xj,yj)=l,■NQj,y)=O

巴区,,%)=0乂(%%)=1

2、单元任一点上三个形函数之和

Ni(x,y)+Nj(x,y)+Nm(x,y)=l即三个形函数中只有两个独立。

3、在三角单元的一个边上

0

如手边，"边的方程式为：y-~~~—(x-x,.)+y;=—~—(x-x(.)+y;

Xi-xtXj—Xj

对于①单元，bm=y,.-y.cm=-(x.-x.)

11

：.y^--(x-xi}+yi

c,“

N,“(x,)')=am+b,„x+c,“一F(X一苫，)+'J[=2(《，，+既,£+%y,)=0

x—xx—X

N,(x,y)=l---------,N,(x,y)=——

Xixrx.,

对于②单元，Nj(x,y),N/x,y)与①相同，而N,“(x,y)=O。

在相邻三角形公共边上，如〃边，

N„,(x,y)=N“(x,y)=。

u=NM+NjUj

*

v=NM+NM

因此在公共边ij边上，位移U,V完全由公共边的(/的位移确定，相邻单元位移是

连续的。

§2-3刚度矩阵

一、单元刚度阵

设单元结点力：｛R『=［q匕U,v｝umvj

单元结点位移：3｝''=［%匕勺,

单元结点虚位移：｛b*『=［阴阳8Uj8vj8um8vJ

则单元内各点的虚位移：｛门=［；,=皿］｛6*｝'

12

*

J

单元内各点的虚应变：k}=£：

.4.

外力虚功为：

风"+见匕+8ujUj+6vyj+8unfJm+8vnym

=应叫8uj8vj3umSvm][UiViUjVjUmVm]'

=({bT)[R}'

单元内应力在虚应变上的虚功为：

/、

巴

[[，(£；%+£；%+匕％)”=[小，(£；£；匕)o-.vdV=^[e]{cy}tdxdy

虚功原理：

({bT)T[/?r=JJ{£*}{b}fdxdy

n{/?『=(

令[幻=\\[B][D][B]tdxdy——单元刚度矩阵

则位}'=内⑹'

对于三角形单元，伏]=[8『[0[8]也

原%kim

网=k'k/k》“

"刈%」6x6

JJ1—M7

„bb.+------cc.pbc.+23

Etrs2r

(r=i,j,ms=i,j,m)

4(1-〃2)△卜N卜

4c也+2"Jcc+—j,b,

rs2_

发M的向

.；做营贺号苧被腐d核‘4支水》'向夕"力

.•.一

-$.......~'

7金@------------

cmMNL

13

二、整体刚度阵

举例说明：

[幻=\\[B][D][B]tdxdy

乂⑴

单元①<

匕⑴

匕⑴

匕⑵

《2)

单元②

匕⑵

乂⑵

14

■U?

匕⑴

匕

K[噌o]

«2

〃⑴+“(2)>(1),n(2)卜⑵

匕⑴+匕⑵V,

单元组装：•>=

“⑴:去⑵"⑴+"⑵“(2)

u；）+u?矶）“32十凡32“33十”33*34«3

匕⑴+匕⑵“(2)L(2)7,(2)V3

0“42”43儿44J

《2）(

.匕⑵%

所以{R}=[K]{6}

Y

其中{R}={R【R；R；{/?,}=(i=l,2,3,4)

M}=⑻8[si3J=削G=1,2,3,4)

三、整刚阵性质

（1）刚度阵每一列元素的物理意义：要是弹性体某一结点在坐标方向发生单位位移,

而其他结点位移都保持为零的变形状态，在所有结点上需要施加的结点力。

（2）刚度阵主元素总是正的。

（3）刚度阵是对称阵。

（4）刚度阵是稀疏阵。

（5）刚度阵在排除刚体位移后是正定阵。

15

§2-4等效结点力

1、集中力等效载荷阵

等效结点力｛a'=出几F—耳」

虚功原理：集中载荷在虚位移上做的功等于等效结点力在结点虚位移上的功。

集中载荷做虚功：Gx6u+Gy6v=｛f｝'｛G｝

等效结点力做虚功：

FM+耳+与@〃)+以阳+工,科,“+Fmy8vm=(M*｝］｛月"

由虚功原理：｛F｝e=｛f｝T｛G｝

而｛门=网肛，即｛/*『=(｛町),［町

=俯｝『｛尸『=(㈤)阿｛G｝

={丹=['/{G}

（M）」G}

其中：但｝'=［此外原5七小/<（M）c{G}>

.（NJ{G}.

（.\Fix\G

{用=/=（M）cG=（M）c{G}（i,j,m）

(N)为形函数N,在集中力作用点处的值。

16

2、表面力等效载荷阵｛。｝"

设表面力（单位面积上的力）｛4｝=/

4,

等效结点力｛。『=[。八2,Qjx0,。““Q,”J

表面力做虚功：j（q@"+q0）tds=｛q/ds

等效结点力做虚功：（桓*｝，｛。｝“

由虚功原理：（•*｝'『｛。『=[｛f｝T｛q｝tds

=｛。｝'=口呵｛。｝"

{板

Q.JM

,[N.q)tds'

其中：QiyQixQiyQmx2,J=0，

Qm.

Q.

{2}={0'}=jM⑷心(bj,m)

匕5J

特例：

3、体积力等效载荷阵｛P｝‘

17

等效结点力{尸}'=出p,pjxpjyp心

体积力做虚功：+py3v）dV=J，{/"}'{p}tdxdy

等效结点力做虚功：（W}丁{P『

由虚功原理：（B}『{P},=\\{f}r{q}tdxdy

={0}'=[p}tdxdy

{p}tdxdy

P.

其中：{P『=出P,PjxPjyPmx%]'=1><JjNj{p*dxdy>

P",.JjN„,{p}fdxdy

7

{《}=,M>=JJM{p}tdxdy（i，j，m）

%

特例：自重p=<>

1-/J

若三种载荷同时作用在单元上，则等效结点力为：{R}'={F『+{Q『+{P『

18

§2-5热应力计算

1、应力应变关系

设a为热膨胀系数，T为温度改变值，卜。}为温度变化造成的应力。

对于平面应力问题：{/}=aT[l10]r

对于平面应变问题：{%}=(l+〃)ani1Of

则应力应变关系为：Q}=。](㈤-{%})

=>{"=⑵(闻⑶{岛})

如：平面应力问题，温度变化和载荷共同作用

2、热载荷阵

考虑热应力，弹性体内应力在虚应变上做的虚功为：

e

价*}'{cr}tdxdy=JJ{Zf[D]([B]{S}-{s0})tdxdy

=^*}C)r(\\[B]T[D][B]{^}etdxdy⑹-JJ网,[。电。}汕办)

外力在结点虚位移上做虚功：

e

由虚功原理得：{R}'+\\[BY[D]{£0}tdxdy^[k]{^}

令{//}'=JJ网,[£>]{£()}冠xdy---热载荷阵，

则:{??}'+{"}'=伙]{阡

对于平面应力问题，如三角形单元：

19

{"}'=0W，［0aTU1Oftdxdy二费/,qh

j%bm%/［，。哂

【说明】

§2-6收敛准则

一、收敛准则

收敛性要求满足：

1、位移模式需包含单元刚体位移

,一工V…一伍=%+%x+a3y

如二角形单兀，1123)

v=a4+%冗+。6y

【从应变角度分析】

单元发生刚体运动且无外力时，应变邑与右，应为零，即:

du-dvdudv八

—=0;—=0;—+—=0

dxdydydx

得：a2=a6=0;%+%=0

[u=a,、

/J13)(1)

[v=a4+a5x

【从刚体位移角度分析】

东若三角形单元绕z轴做刚体转动跳角，则：

s，、〜口«=sin0=-co^y

任后、一点A(x,y)的位移是：，

v=rcoQcos0=co[}x

20

*若三角形还有沿x轴和y轴的刚体平移位移“0,%

则总位移为：…。一阳⑵

+为x

-&,

对比（1）、（2）可得：a}-M（）；%=%；«3-0；。5=。0

最终可得：/=〃0；。2=0；a3=-COQ；%=%；%=利）；。6=0

可见：«,%4%提供刚体位移。

2、位移模式须包含常应变

应变包括常应变和变应变，

dudvdudv

£*==%；f=—=<z-，r,=—+—=a+a

oxoyv6•onyoxi5

可见：a?aba3%提供常应变。

结论：位移模式中应包含一次项和常数项（对于三角单元，4%不仅是刚体位移的需

要，也是常应变的需要）。

3、位移模式须在单元内位移连续、相邻单元间位移协调

单元内连续（单元位移不间断）——选多项式构造位移模式即可满足。

单元间协调（单元间不开裂，也不挤入）——单元交界面上位移由该交界面上的结点位

移确定。

总结：满足条件1、2——完备单元

满足条件3——协调单元

满足条件1、2、3——完备协调单元（收敛性要求）。

二、多项式位移模式的选择

21

要求：1、几何各向同性——不应有一个偏惠的坐标方向

1

x:y

x2xyy2

x3x2y:xy2y3

x4x2y2xy3y4

2、多项式项数等于或稍大于单元结点自由度数

§2-7实施步骤

一、求解步骤

二、注意事项

1、对称性的利用

22

结构对称、受力反对称——变形反对称

2、结点选择、单元划分

(1)以集中载荷作用点，分布载荷起终点，约束支承点作为结点。

(2)物体厚度有变，或有不同材料组成，则不要把厚度不同或材料不同的区域划在

一个单元内。

(3)在应力变化剧烈的区域，单元应划分的密一些。

(4)单元各条边长度不要悬殊太大。

3、结点编号

同一单元相邻结点号之差应尽量小。

23

4、边界条件处理

采取指定位移处理方法。如:

方法一：修改矩阵法

方程改写为:

1000[〃[B、

。七2夫2-后2圈-氏23夕3

0010%A

。”44_了2七一”414一火43A

0%42

求出匕后代入原式，可求/?1,尺2。

方法二：置大数法

方程改写为：

%「十%12匕3%14MrW5

“21及22^23%24

匕>—<

儆31。」

2。5

&32攵34〃2

人攵42上43攵44_

近似求出匕,岭后代入原式，可求飞,凡。

24

5、应力计算结果处理

（1）绕结点平均法

=/4加+(%)⑵+(%)⑶+(%)(4)+(q)⑸+(q)⑹]

(2)单元平均法

=，9儿）+C

3）=；[（%）⑵+（/）（”]

A、B、C为公共边界中点。

W）c=g3、）⑶+。）（4）】

§2-8矩形单元

25

设任一点坐标：］x=%+殆一点位置映射

其中（了0，汽）为矩形单元中点坐标。

设该点位移为：

।2：3/4：/―一点位移映射（位移模式）

由结点位移，如（-1,-1）—＞（%,匕）

4

Z

〃=

E

nw

r-//

l=

4lAI

NV

V=z匕

/=i

0NQN.0N0

其中[N]=<24

N、0N20华0%

N,=(l+4)(1+%)/4,4。=基,%=mi=l,2,3,4

总结:

u=u{x,y)

v=v(xj)

v=v(f,〃)

，13

〃

a3g

£史

X1aN

由

于㈤

-£>--<

y/axb18-

a-v

£I

Z史1a1

ay〃

/++

a-w-1一451

私b8"S

ayk

=>{£}=网.『

=>{<7}=[0网同

n{k}=^[BY[D][B]tdxdy=ff[BY[D][B]tabd^d?j

26

等效结点力计算同三角形单元。

【说明】

(1)位移模式中有常数项、一次项——完备元。

(2)边界位移协调——协调元。

例如图中单元：在kj边上,

对于①单元，J=1

u-ay+a2+a3rj+a4rl=@+a2)+(%+%)〃=4+

_uj+uk

,1-2

由「

uk=At+A2._Uk-Uj

.2-2

对于②单元，^=-1

u=a；-a2+a3rl-%〃=(a；-&)+Q-%j〃=A；+A2rj

由<J_=><

限=4+44=七乜=4

I2-

所以〃在以边上协调，同理v在灯边上也协调。

(3)矩形单元不是常应力常应变单元。

q沿y方向线性分布，o■，沿x方向线性分布，%沿小y两方向线性分布。

(4)矩形单元的缺点：

1)不能适应斜边界，曲线边界；

2)不便于对不同部位采用不同大小的单元。

27

§2-9三角形单元面积坐标

ATAijin的面积;

A.->\pjm的面积；Apini的面积；即炉的面积

令L=幺，L=当,L=4

AAA

4,4,L,“为p点面积坐标。

A+匕+Lm=14，乙，4中只有两个独立。

可以证明：L：=Ni,Lm=Nm

x=xiLi+xjLj+xmLmuL+UjLj+u“,L“,

y=+刀4+y,„L„,v=v.L.+v,L.+vmLm

【说明】

我为弟无

点.我船

注：Lm=\-L-Lj

28

X=xL+XjLj+xJ,“

单元内任意一点坐标:

y=M,+yjLj+y“L“

u=uL+uL+uL

而位移模式：iiJJmm

VVLVLVL

^ii+jj+,nm

{/}=：=卬]的

0L,0L0

其中W]=m

L,0儿0Lm

可见:

点坐标（见”,疝,协一|点位移（U,V）

变换飞

位置映於

\/

局部坐标&也）

注：点位置映射的变换与位移模式的变化在形式上可以相同，也可以不相同。

史

&

V加

--一»=网{亚

办

包dv

办十一

dx

Q}=⑵网⑶

{k}=\\[Bir[D][B]tdxdy

三角形面积上积分：口"阻9=31省产-，⑼

a邛、

三角形边上积分：'以4ds■I(ijm)

(a+〃+l)!

其中a;/?7为整常数。

1!1!…1

如：jj^Ljdxdy----------2A=-----------2A=A(a=1;£=l;y=0)

(1+1+2)!4x3x2xl

29

§2-10习题讲解

一

力

彳

一

N

已知：平板尺寸如图所示，厚度为f,受载荷集度为4的均布载荷作用，划分2个单元

求：单元应力。

解：

1、划分单元，结点编号

2、求单刚阵

对单元①：对单元②:

噌K7

伙］⑴=噌成

图）成

3、组装总刚阵

7-4-4-2-3-200

413-2-12-2-100

靖）k⑴

鼠12〜30-4-27004-3-2

码）成+磴掰）+成）唠_3Et-2-1201340-2-1

［勺=

6⑴4〃2)~12

研假32十人32“33十代33媾-3-20470-4-2

0娉)-2-140013-2-12

00-3-2-4-274

00-2-1-2-12413

4、计算等效结点力：

[01

q=\\

｛。｝'=]阿｛加5=[&Qy&Q3y*Q,J=[0000_乎了

原结构变为：

30

7-4-4-2-3-20