一维梁的MATLAB有限元法分析

一维梁的M A T L A B有限元法分析
问题如下：
梁A B在A和B两端固定，中间点表示为C，在中间区域承受均匀分布的载荷q，如图所示。

梁A B的抗弯刚度为E I。

1，使用R i t z法确定点C处的位移和弯矩，并讨论随着包含更多基本函数的准确性。

提示：根据偏转曲线的形状，可以选择基函数的形式，其中系数，并且应该由点A或B处的边界条件确定。

2，采用一维有限元法解决问题，并讨论网格越细时的准确性。

提示：使用1-D梁单元。

1R i t z法
一维欧拉-伯努利梁的势能如下：
设选择基函数，容易看出基函数满足边界条件
设，
代入势能表达式得到
由于三角函数是正交函数系，所以得到
令q=10N/m m，E=200000M P a，I=10000，L=200m m
在M A T L A B中计算A k前十项
得到
A=
-0.008400754770396-0.000320811945459-0.000049922673823 -0.000020050746591-0.000009258470170-0.000003960641302 -0.000001943426801-0.000001253171662-0.000000837688763 -0.000000513299113
计算C点位移,使用1-10个试函数结果如下：
0.0084007547703960.008400754770396
0.008500600118041
0.0085006001180410.008519117058380
0.008519117058380
0.0085230039119830.008523003911983
0.0085246792895100.008524679289510
计算C点弯矩,,使用1-10个试函数结果如下：
0.8291212625436840.7024697829907620.746814316705690 0.7151514468174600.7379958063005830.723923419683592 0.7333220380018130.7254063205297550.732103122461494 0.727037063279377
可以看到，位移收敛是很快的，弯矩收敛速度慢于位移。

M A T L A B代码如下
%量纲为mm，N，k g
%定义弹性模量泊松比,定义截面惯性矩
E=200000;v v v v v=0.3;I=10000;
%定义梁的长度
L e n g t h=200;L=L e n g t h;
%定义载荷集度
P=10;
%定义基函数数目
n=10;
s y ms x;
%计算A j
A=z e r o s(1,n);
f o r j=1:n
A(j)=(L^3/(8*j^4*p i^4*E*I))*i n t(P*r e c t a n g u l a r P u l s e(0.25*L e n g t h,0.75*L e n g t h,x)*(c o s(2*j*p i*x/L) -1),0,L);
e n d
A
%计算C点位移与弯矩
U C=z e r o s(n,1);x C=L/2;
U C(1)=A(1)*(c o s(2*p i*1*x C/L));
f o r j=2:n
U C(j)=U C(j-1)+A(j)*(c o s(2*p i*j*x C/L)-1);
e n d
U C
M C=z e r o s(1,n);
M C(1)=-A(1)*(-(2^2*p i^2*1^2/L^2)*c o s(2*p i*1*x C/L));
f o r j=2:n
M C(j)=M C(j-1)-A(j)*(-(2^2*p i^2*j^2/L^2)*c o s(2*p i*j*x C/L));
e n d
M C
2有限元法
有限元法理论不再叙述，M A T L A B代码如下
%量纲为mm，N，k g
%定义弹性模量泊松比,定义截面惯性矩
E=200000;v v v v v=0.3;I=10000;
%定义梁的长度,定义单元数
L e n g t h=200;E l e m e n t N u m=40;
%定义载荷集度
P=10;
%定义符号变量
s y ms x;
L=L e n g t h/E l e m e n t N u m;
%代码段建立单元刚度矩阵
K e=E*I/(L*L*L)*[126*L-126*L;6*L4*L*L-6*L2*L*L;-12-6*L12-6*L;6*L2*L*L-6*L4*L*L]; %代码段建立总体刚度矩阵K
i=1;d o f i=1;d o f j=1;
K=z e r o s(2*(E l e m e n t N u m+1));
K g=z e r o s(2*(E l e me n t N u m+1));%K g是工具矩阵
f o r i=1:E l e m e n t N u m
f o r d o f i=1:4
f o r d o f j=1:4
K g(d o f i+2*(i-1),d o f j+2*(i-1))=K e(d o f i,d o f j);
e n d
e n d
K=K+K g;
K g=z e r o s(2*(E l e me n t N u m+1));
e n d
K K=K;
%建立局部坐标系下的单元形函数矩阵
e=x/L;
N1=1-3*e*e+2*e*e*e;
N2=L*(e-2*e*e+e*e*e);
N3=3*e*e-2*e*e*e;
N4=L*(e*e*e-e*e);
N=[N1N2N3N4];
%r e c t a n g u l a r P u l s e(x,0.25*L e n g t h,0.75*L e n g t h);
%代码段建立总体外载荷向量P G,G~g l o b a l
P G=z e r o s(2*(E l e m e n t N u m+1),1);
P G g=z e r o s(2*(E l e m e n t N u m+1),1);
i=1;d o f i=1;
f o r i=1:E l e m e n t N u m
e=(x-(i-1)*L)/L;
N1=1-3*e*e+2*e*e*e;
N2=L*(e-2*e*e+e*e*e);
N3=3*e*e-2*e*e*e;
N4=L*(e*e*e-e*e);
N=[N1N2N3N4];
P g=i n t(P*r e c t a n g u l a r P u l s e(0.25*L e n g t h,0.75*L e n g t h,x)*N,(i-1)*L,i*L);
f o r d o f i=1:4
P G g(2*(i-1)+d o f i)=P g(d o f i);
e n d
P G=P G+P G g;
P G g=z e r o s(2*(E l e me n t N u m+1),1);
e n d
%代码段求解线性方程组，采用主对角线置大数法加载约束,将节点位移存入向量q
K(1,1)=99999999999999;
K(2,2)=99999999999999;
K(2*E l e m e n t N u m+1,2*E l e m e n t N u m+1)=9999999999999;
K(2*E l e m e n t N u m+2,2*E l e m e n t N u m+2)=9999999999999;
q=K\P G;
%输出C点挠度U C到控制台
U C=q(E l e m e n t N u m+1);
i f m o d(E l e me n t N u m,2)
q e=[q(E l e me n t N u m);q(E l e me n t N u m+1);q(E l e me n t N u m+2);q(E l e me n t N u m+3)];
e=x/L;N1=1-3*e*e+2*e*e*e;N2=L*(e-2*e*e+e*e*e);N3=3*e*e-2*e*e*e;N4=L*(e*e*e-e*e);
N=[N1N2N3N4];
U=N*q e;
U C=50000*i n t(U,0.5*L-0.00001,0.5*L+0.00001);
e n d
U C
网格与中点位移关系如下表
单元数目C点位移
20.016927089638169
50.016893756304203
100.016927089638188
150.016926678114924
200.016927089638213
图线如下结果呈振荡收敛。

本章附录：有限元法程序正确性检验
将以下两行注释掉，边界条件即转化为悬臂梁
K(2*E l e m e n t N u m+1,2*E l e m e n t N u m+1)=9999999999999;
K(2*E l e m e n t N u m+2,2*E l e m e n t N u m+2)=9999999999999;
计算得到端点位移为
a n s=
0.437500018452276
A N S Y S软件计算得到0.44326，十分接近。

（A N S Y S中建立正方形截面，其惯性矩略小于10000）
由于两端固定支撑梁问题的特殊性，不以此检验。

代码如下：
%量纲为mm，N，k g
%定义弹性模量泊松比,定义截面惯性矩
E=200000;v v v v v=0.3;I=10000;
%定义梁的长度,定义等分数
L e n g t h=200;E l e m e n t N u m=100;
%定义载荷集度
P=10;
%定义符号变量
s y ms x;
L=L e n g t h/E l e m e n t N u m;
%代码段建立单元刚度矩阵
K e=E*I/(L*L*L)*[126*L-126*L;6*L4*L*L-6*L2*L*L;-12-6*L12-6*L;6*L2*L*L-6*L4*L*L]; %代码段建立总体刚度矩阵K
i=1;d o f i=1;d o f j=1;
K=z e r o s(2*(E l e m e n t N u m+1));
K g=z e r o s(2*(E l e me n t N u m+1));%K g是工具矩阵
f o r i=1:E l e m e n t N u m
f o r d o f i=1:4
f o r d o f j=1:4
K g(d o f i+2*(i-1),d o f j+2*(i-1))=K e(d o f i,d o f j);
e n d
e n d
K=K+K g;
K g=z e r o s(2*(E l e me n t N u m+1));
e n d
%建立局部坐标系下的单元形函数矩阵
e=x/L;
N1=1-3*e*e+2*e*e*e;
N2=L*(e-2*e*e+e*e*e);
N3=3*e*e-2*e*e*e;
N4=L*(e*e*e-e*e);
N=[N1N2N3N4];
%r e c t a n g u l a r P u l s e(x,0.25*L e n g t h,0.75*L e n g t h);
%代码段建立总体外载荷向量P G,G~g l o b a l
P G=z e r o s(2*(E l e m e n t N u m+1),1);
P G g=z e r o s(2*(E l e m e n t N u m+1),1);
i=1;d o f i=1;
f o r i=1:E l e m e n t N u m
e=(x-(i-1)*L)/L;
N1=1-3*e*e+2*e*e*e;
N2=L*(e-2*e*e+e*e*e);
N3=3*e*e-2*e*e*e;
N4=L*(e*e*e-e*e);
N=[N1N2N3N4];
P g=i n t(P*r e c t a n g u l a r P u l s e(0.25*L e n g t h,0.75*L e n g t h,x)*N,(i-1)*L,i*L);
f o r d o f i=1:4
P G g(2*(i-1)+d o f i)=P g(d o f i);
e n d
P G=P G+P G g;
P G g=z e r o s(2*(E l e me n t N u m+1),1);
e n d
%代码段求解线性方程组，采用主对角线置大数法加载约束,将节点位移存入向量q
K(1,1)=999999999999999;
K(2,2)=999999999999999;
%K(2*E l e me n t N u m+1,2*E l e me n t N u m+1)=99999999999999;
%K(2*E l e me n t N u m+2,2*E l e me n t N u m+2)=99999999999999;
q=K\P G;
%输出C点挠度U C到控制台
U C=q(E l e m e n t N u m+1);
i f m o d(E l e me n t N u m,2)
q e=[q(E l e me n t N u m);q(E l e me n t N u m+1);q(E l e me n t N u m+2);q(E l e me n t N u m+3)];
e=x/L;N1=1-3*e*e+2*e*e*e;N2=L*(e-2*e*e+e*e*e);N3=3*e*e-2*e*e*e;N4=L*(e*e*e-e*e);
N=[N1N2N3N4];
U=N*q e;
U C=50000*i n t(U,0.5*L-0.00001,0.5*L+0.00001);
e n d
U C
%输出端点位移
q(2*(E l e me n t N u m+1)-1)
边界元法求解弹性波散射问题
2.考虑无限扩展各向异性弹性固体介质中的波散射问题。

在x3方向上具有无限长度的圆柱形空腔切割介质内，半径为r。

现在，入射平面波以正x1方向传播，置换场：入射波被腔散射。

事件波在腔中散射。

问题：请用B E M求解散射问题，并表达散射和总波场沿腔边界的位移。

提示：
（1）散射波是典型的S H波场，满足控制方程：
（2）总波场边界上的合力为零。

即.
（3）使用常量元素。

1问题的转化
根据弹性波的相关理论，我们假设散射波在x1与x2方向无位移，在x3方向做如下形式振动：
由于入射波的平移不变性，u3在x3方向恒定，代入波动方程即可得到
由于问题中空腔内部不存在任何物质，因此
该问题归结为求解亥姆霍兹方程第二边值问题。

，
直接离散得到。