无网格法大作业

二、利用伽辽金无网格法实现悬臂梁弯曲问题

算例描述：悬臂梁自由度受集中力作用如图1所示，编写2D FEGM 程序计算位移、应变、应力。

梁的尺寸：L = 42m , D = 12m , 厚度默认为1m ; 弹性常数：E = 30M Pa, ν = 0.3;

不考虑梁的自重，在自由端施加的集中载荷为P = -1000 N ;

图1自由端受集中力作用的悬臂梁

弹性力学解析解

x 方向位移：22

(,)[(63)(2)(+)(]642Py D D u x y L x x y D y EI ν=--++-+ y 方向位移：22

2(,)[3()(45)(3)]64

P D x v x y y L x L x x EI νν=-+++- 梁截面的法向应力：()(,)xx P L x y

x y I

σ-=- y 方向正应力：0yy σ= 梁截面剪应力：2

2(,)()24

xy P D x y y I τ=-

应变能：1 4.452T

E D d εεΩ

=

Ω≈⎰ 1、 节点间距收敛性分析

为研究节点间距对计算结果的影响，考虑如下5种等间距节点分布情况： 3X3(轴向X 竖向)、6X3、11X3、11X5、11X7。依次增加轴向和竖向的节点数，考察应变能残量和悬臂梁中心轴y=0的计算位移的收敛性，以及某一梁截面(x=23.667)的计算应力收敛性。此部分程序与相应结果见附录1。

不同节点网格数目对能量误差范数的影响

不同节点数目计算得到中心轴y=0的梁的位移与解析解比较

0102030

40

50607080

0.10.20.30.40.50.60.70.8

0.9场节点数目

能量误差范数

5

10

15

20

25

30

35

40

45

50

-3

x

位移

不同节点数目计算得到x=23.667的截面法向应力与解析解比较

不同节点数目计算得到x=23.667的截面剪应力与解析解比较

从上述结果可以发现，随着节点数的增加，FEGM 计算收敛于解析解。考虑应变能和位移以及应力的收敛性可以发现，节点数从3X3增加至6X3，应变能残量迅

-6

-4-20246

-1500

-1000

-500

500

1000

1500

y

s t r e s s (x x )

-6

-4

-2

0246

-140

-120-100-80-60-40-200

204060y

s t r e s s (x y

)

速减小，中心轴位移迅速向解析解逼近，截面正应力始终非常接近于解析解，而截面剪应力则不然，虽然有向解析解逼近的趋势，但是与解析解相距仍甚远。节点数进一步增加至11X3，发现不管是应变能残量还是截面应力都没有明显进一步收敛至解析解，究其原因是仅增加轴向节点数，导致不能进一步精确的描述截面剪应力。当增加竖向节点数后，可以发现截面剪应力迅速逼近解析解，从而应变能残量随之减小，再进一步增加竖向节点数，不管是位移还是应力已基本没有太大变化，且与解析解吻合很好。

2、 高斯积分点数的影响

为了探讨不同高斯积分点对计算结果的影响，分别进行2、3、4、5、6积分点的高斯积分与解析解比较，此时固定场节点为11X5。计算的实现通过调整子程序pguass 的参数。此部分程序与相应结果见附录1。

不同高斯积分点数目对能量误差范数的影响

高斯积分点数

能量误差范数

不同高斯积分点数目计算的梁中心轴的位移与解析解间的关系

不同高斯积分点数目计算的梁x=24截面法向应力与解析解间的关系

5

10

15

20

25

30

35

40

45

50

x 10

-3

x

位移

y

n o d e s -s t r e s s (x x )

不同高斯积分点数目计算的梁x=24截面法向应力与解析解间的关系

通过改变高斯积分点数目的计算结果可以发现，高斯积分点数为2、3、4、5、6的计算结果差异很小，且与解析解吻合很好。说明当高斯积分点数不小与2时，增加高斯积分点数对自由端受集中载荷悬臂梁计算问题影响甚小，2点高斯积分已足够精确。

3、 使用罚函数方式施加边界条件，并于原拉氏乘子法比较 罚函数法介绍

将位移边界条件引入伽辽金弱形式

()()0t

T T T u u f d u td δδεσδδΩ

Γ∏=-+Ω+Γ=⎰⎰

()()()0u

T p u u u u u d δδαδΓ∏=∏+-Γ=⎰

式中α为罚系数，一般取为103E~107E ，E 为弹性模量。 则得系统总体的离散方程：

()u

u

p p T p T p K K d P P K N Nd P N ud ααΓΓ+=+=Γ=Γ

⎰⎰

根据上述理论修改原主程序(见附录 1 gbmm2Dpf)，此时仍保持场节点数位11X4，四点高斯积分，罚系数α取103E 。此部分程序与结果见附录1。

-6

-4-20246

y

n o d e -s t r e s s (x y )