用位移变分法求解问题

用位移变分法求解问题
一 求解问题：
铅直平面内的正方形薄板，边长为2a ，四边固定，图1，只受重力作用。

设μ=0，试取位移分量的表达式为
2222
1232222u (1)1(),
a a a a a a
x y x y x y A A A ⎛⎫=--+++
⎪⎝⎭
2222
1232222v (1)(1)(),
a a a a
x y x x B B B =--+++
用里兹法或伽辽金法求解（在u 的表达式中， 布置了因子x 和y ，因为按照问题的对称条件， u 应为x 和y 的奇函数）。

二 求解步骤：
取坐标如图1所示，位移分量为
2222
1232222u (1)1(),
a a a a a a
x y x y x y A A A ⎛⎫=--+++
⎪⎝⎭
2
2
2
2
1232222v (1)(1)(),
a a a a
x y x x B B B =-
-+++
上述位移分量满足位移边界条件
()a u 0x =±= a (u )0y =±=
a (v )0
y =±= a (v )0x =±=
现在，式（1）取各式前两个系数为待定系数，也就是取
221
22u (1)1a a a a
x y x y
A ⎛⎫=-- ⎪⎝⎭
22
1
22v (1)(1)a a x y B =--
将式（2）代入式（3），即
()222ε1-μE u v u v v u V 2μd d 2(1-μ)2x y x y x y x y ⎡⎤⎛⎫⎛⎫∂∂∂∂∂∂⎛⎫
=++++⎢⎥ ⎪ ⎪ ⎪∂∂∂∂∂∂⎝⎭⎢⎥⎝⎭⎝⎭⎣⎦⎰⎰
其中，μ=0，得
221111ε646464E V -+217522515A B A B ⎛⎫
= ⎪
⎝⎭
图1
（1）
（2）
（3）
（4）
因为在这里f f 0,f 0,f ρg,m=1x y x y --
====可得
1V 0A ε∂=∂ 22
221V ρg(1)(1)d d B a a x x x y
ε∂=--∂⎰⎰
由式（5）得
1132E
(18-7)01575A B =
2
11-32E 16ρga (-30)=2259A B
解式（6）得，
2
1175
=
ρga 1066E A 2
1225=ρga 533E B
将式（7）代入式（2）得位移函数
22
22175u=(1)(1)ρg 1066E a a x y xy --
222
22255v=(1)(1)ρga 533E a a x y xy --
将位移函数式（8）代入位移-应力函数中，可得
22
x 221753σ(1)(1)ρg 1066a a y x y
=-- 2
y 2450σ(1)ρg 533a x y
=--
2222
24225τ(2915217)ρg
2132a a a xy y x y x x =---+
三 MATLAB 求解
应力空间图
这里取，
22
ρ=1KN/m g=1m/s a=1m 由应力函数可到以下图2，图3，图4。

（5）
（6）
（7）
（8）
（9）
分别代表x 方向正应力，y 方向正应力，剪应力
图3
y
σ方向正应力 图2 x σ方向正应力
x
σ
y
σ
由应力图很容易观察到图像关于x 轴y 轴的对称性，对比应力的函数的奇偶对称性正好相对应。

图2知x 方向的正应力峰值出现在薄板上（±a,±0.5a ）；由图3知y 方向正应力最大出现在（0，±a ）;由图4切应力最大值出现在（±a,0）。

其中，max σ 0.1231MP x a = max σ 0.8443MP y a = max τ 0.4221MP xy a =
三 总结
创新实践将在弹性力学的不容易求解的问题借助于软件得以实现，一方面，从思维上培养了自己做事情认真严谨的作风；另一方面，从动手能力方面得到提升。

创新实践不单单是做一个题目的问题，更加重要的是从做题中在，做到了对理论的复习和提升；做到了对软件的全面掌握，并熟练度得到提高；做到了一种分析问题借，助工具解决问题的能力。

图4
xy
τ剪应力 xy
τ
附录：
clear
clc
syms x y u v a A B dux dvy duy dvx V p g E;
u=(1-x^2/a^2)*(1-y^2/a^2)*(x/a)*(y/a)*A;
v=(1-x^2/a^2)*(1-y^2/a^2)*B;
dux=diff(u,x);
dvy=diff(v,y);
dvx=diff(v,x);
duy=diff(u,y);
duy=SIMPLIFY(duy);
dux=SIMPLIFY(dux);
dvy=SIMPLIFY(dvy);
dvx=SIMPLIFY(dvx);
dux2=dux^2;
dvy2=dvy^2;
dxy=(dvx+duy)^2;
dux2=SIMPLIFY(dux2);
dvy2=SIMPLIFY(dvy2);
dxy=SIMPLIFY(dxy);
ji=dux2+dvy2+0.5*dxy;
ji=SIMPLIFY(ji);
V=int(int((ji),x,-a,a),y,-a,a);
V=SIMPLIFY(V);
V=0.5*E*V;
m=p*g*int(int((1-x^2/a^2)*(1-y^2/a^2),x,-a,a),y,-a,a);
m=SIMPLIFY(m);
dva=diff(V,A);
dvb=diff(V,B);
dva=SIMPLIFY(dva);
dvb=SIMPLIFY(dvb);
[A,B]=solve('32/1575*E*(18*A-7*B)=0','-32/225*E*(A-30*B)=16/9 *p*g*a^2','A','B');
A=A
B=B
u=(1-x^2/a^2)*(1-y^2/a^2)*(x/a)*(y/a)*A;
v=(1-x^2/a^2)*(1-y^2/a^2)*B;
fx=E*diff(u,x);
fy=E*diff(v,y);
fxy=0.5*E*(diff(v,x)+diff(u,y));
fx=simplify(fx)
fy=simplify(fy)
fxy=simplify(fxy)
p=1;
g=1;
a=1;
fx=175/1066*(a^2-y^2)*y*p*g*(-3*x^2+a^2)/a^4
fy=-450/533*(a^2-x^2)*y*p*g/a^2
fxy=-25/2132*x*p*g*(29*a^4-15*y^2*a^2-21*y^2*x^2+7*x^2*a^2)/a ^4
ezmeshc(fx,[-1,1,-1,1])
figure
ezmeshc(fy,[-1,1,-1,1])
figure
ezmeshc(fxy,[-1,1,-1,1])
（3程序运行主要结果
A=175/1066/E*p*g*a^2
B=225/533/E*p*g*a^2
fx=175/1066*(a^2-y^2)*y*p*g*(-3*x^2+a^2)/a^4
fy=-450/533*(a^2-x^2)*y*p*g/a^2
fxy=-25/2132*x*p*g*(29*a^4-15*y^2*a^2-21*y^2*x^2+7*x^2*a^2)/a ^4。