Euler-Bernoulli beam theroy(欧拉伯努利梁理论与有限元计算)

Euler-Bernoulli beam

一、理论部分

Euler-Bernoulli beam 假设

()()

'000u u yv v v x u u x ⎧=-⎪

=⎨⎪

=⎩ (1.1)

由式可得

'

''00

xx yy xy u yv εεγ⎧=-⎪

=⎨⎪

=⎩ (1.2)

虚位移原理 0int ext P P δδ+=

(1.3) 其中

d d ext P f u V T u S δδδ=+⎰⎰

(1.4)

()()()/2

/20/2

'''

'''

/20

''''

''

e

e

e int h l xx

xx

h h l h l P dV

b dxdy

b E u yv u y v dxdy EAu u EIv v dx

δσδεσδεδδδδ--=-=-=---=-+⎰⎰

⎰⎰⎰⎰

令1

2

e x l ξ+=，e l 为单元长度，则上式成为

()1

''

''''0012e int l

P EAu u EIv v d δδδξ-=-+⎰

(1.5)

单元节点位移取为 {}()

''01110222,,,,,T

q u v v u v v =

(1.6)

令

{}(){}{}(){}

01212340000u u u N N q v N N N N q = 0 0 ⎧⎪⎨

= ⎪⎩ (1.7)

其中形函数

1212

12

u u N N ξξ-⎧=⎪⎪⎨

+⎪=⎪⎩ (1.8)

()()()()()()()()

2122232

41124118

1124118e e N l N N l N ξξξξξξξξ⎧=-+⎪⎪

⎪=-+⎪⎨

⎪=+-⎪⎪⎪=+-+⎩

(1.9)

分别对式、求一阶和两阶导数得

'

1'212

12u u N N ⎧=-⎪⎪⎨

⎪=⎪⎩ (1.10)

''1''2''3''432314

4323144e e

N N l N N l ξξξξ⎧=⎪⎪

⎛⎫⎪=- ⎪⎪⎪⎝⎭⎨

⎪=-⎪⎪⎛⎫⎪=+ ⎪⎪⎝⎭⎩

(1.11)

由

2

e

x l ξ∂=∂，可得 '2ui ui e N N x l ∂=∂ (1.12)

2

2

2'''''2

2i i

i i i e N N N N N x x x

x x x l ξξ⎛⎫∂∂∂∂∂∂⎛⎫⎛⎫⎛⎫==== ⎪ ⎪ ⎪ ⎪

∂∂∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭

(1.13)

将、式代入式，可得

{}

()'11''

12'

12''1''2''''''''1234''3''4002000000000u T

int u u e u N P q EA N N d l N N N EI N N N N N N δδξ-⎧⎛⎫⎪ ⎪

⎪ ⎪

⎪ ⎪ ⎪ ⎪=- ⎨ ⎪⎪

⎪⎪ ⎪⎪ ⎪ ⎪⎝⎭⎩

⎛⎫

⎪

⎪

⎪ ⎪ +0 ⎪ ⎪ ⎪ ⎪ ⎝⎭⎰(){}3112e d q l ξ-⎫⎪⎪⎪⎛⎫⎪ ⎬ ⎪⎝⎭⎪

⎪⎪⎪⎭

⎰ (1.14)

刚度矩阵

e e e

u v K K K =+ (1.15)

()0'11''12'12002000000u e

u u u e u N K EA N N d l N ξ-⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪

⎪ ⎪ ⎪ ⎝⎭⎰

(1.16)

()''13''120''''''''12341''3''40200v e N N K EI N N N N d l N N ξ- ⎛⎫ ⎪ ⎪

⎪ ⎛⎫ ⎪=0 ⎪ ⎪⎝⎭

⎪ ⎪ ⎪ ⎝⎭

⎰

(1.17)

二、算例

本文所举算例为悬臂梁端部受竖向载荷P ，挠度理论曲线为

23

26Pl P y x x EI EI

=-+ 具体数据如下：

0.02b m =，0.05h m =，1l m =，82.110EA N =⨯，424.37510EI N m =⨯⋅，10P kN =

划分10个单元，11个结点。

00.10.20.30.4

0.50.60.70.80.91

x(m)

t h e d e f l e c t i o n o f t h e b e a m (m )

the deformation of the beam

由上图可以看出，数值模拟结构跟理论解吻合非常好。

源代码(Matlab)

% program_main initial_statistic;

le=pretreatment(ND,NE,lamda,x,y);

[Dis,K]=Euler_Bernoulli_beam(ND,NE,lamda,EA,EI,le,P,Node_c,Dis_c);

post(Dis,x,y,ND,NE);

% initial statistic; ND=11;NE=10;

EA=*ones(NE,1);EI=*ones(NE,1); x=(0::1)';y=zeros(11,1);

lamda=[1,2;2,3;3,4;4,5;5,6;6,7;7,8;8,9;9,10;10,11;]; P=zeros(3*ND,1);P(32)=;

Node_c=[1,2,3]';Dis_c=[0,0,0]';

function le=pretreatment(ND,NE,lamda,x,y) % the program of pretreatment; for i=1:NE

P