无阻尼 自由振动 暂态问题
利用有限元和直接积分求解一暂态问题(一维波动--无阻尼自由弦振动),推导了有限元的单元公式,见附图(本人推导,希望被指出错误).在集成了M矩阵和刚度矩阵K之后,利用简单的中心差分的方法,来求时间域状态.原理并不复杂.
可是MATLAB程序的图形和解析解相差太多.请好心的,热心的高手指点! :-)
================================================================================
程序如下:
clear,close all,clc;
NELES=419; %单元数
dt=0.001; %时间步长
NNODS=NELES+1; %节点数
L=1/NELES; %弦长为1,除以单元数为每个单元尺寸
x=linspace(0,1,NNODS)'; %各点x轴坐标
K=zeros(NNODS); M=zeros(NNODS); %K, M size
u=zeros(NNODS,2001); %每列为一时间状态,整体时间为0.002*2000=2秒
u(2*NNODS/7+1:3*NNODS/7,1)=sin(pi*x(2*NNODS/7+1:3*NNODS/7,1)*7); %初始位移为一小段正弦波
N11= @(x)1/L*(x-L).*(x-L)/L;%形函数N1=1/L*(x-L),此处为N1*N1
N12= @(x)-1/L*(x-L).*x/L; %形函数N2=x/L,此处为N1*N2
N22= @(x)x/L.*x/L; %此处为N1*N2
K11 = 1/L; %每个单元的2*2刚度矩阵第一个元素为1/L
K12=-K11; %此处实际为(dN1/dx)*(dN2/dx)在单元内积分
ke=; %K为对称矩阵
for ne=1:NELES %开始集成M,K
lim1=(ne-1)*L;lim2=ne*L;
M11 = quad(N11,lim1,lim2);
M12 = quad(N12,lim1,lim2);
M22 = quad(N22,lim1,lim2);
me=;%计算每个单元的m矩阵
M(ne,ne)=M(ne,ne)+me(1,1);
M(ne,ne+1)=M(ne,ne+1)+me(1,2);
M(ne+1,ne)=M(ne+1,ne)+me(2,1);
M(ne+1,ne+1)=M(ne+1,ne+1)+me(2,2);%集成m到整体M
K(ne,ne)=K(ne,ne)+ke(1,1);
K(ne,ne+1)=K(ne,ne+1)+ke(1,2);
K(ne+1,ne)=K(ne+1,ne)+ke(2,1);
K(ne+1,ne+1)=K(ne+1,ne+1)+ke(2,2); %因为每个单元的k矩阵都相等,直接集成总体K
end
M(NNODS,:)=[];M(1,:)=[];M(:,NNODS)=[];M(:,1)=[]; %因为边界条件,去掉M首,末的行与列
K(NNODS,:)=[];K(1,:)=[];K(:,NNODS)=[];K(:,1)=[]; %因为边界条件,去掉K首,末的行与列
u(2:NNODS-1,2)=2*M/dt/dt\((2*M/dt/dt-K)*u(2:NNODS-1,1)); %由推导的公式计算u2(第2时刻的位移)
for n=2:2000 %由第一时刻位移u1(初始条件)和第二时刻位移u2,开始用显式的中心差分法迭代
u(2:NNODS-1,n+1)=M/dt/dt\((2*M/dt/dt-K)*u(2:NNODS-1,n)-M/dt/dt*u(2:NNODS-1,n-1));
end
恳请各位指正! 看看！
你的matlab是哪个版本的，在6.5中出现问题？
回复 #2 无水1324 的帖子
我的matlab是7.0 问题是程序可以运行,可是把算出来的结果做了animation,明显是不对的.反复检查又找不到错在哪里.:@L
还是调试出来了
自问自答:原来问题是出在形函数积分的计算上,如果按原程序中的代码,则每个单元的对shape function的积分都是相等了,实际上是不同的.
由于一开始试算的时候,最从最简单的划分两个单元开始,这种情况下对N1和N2的积分是相等的,但不是所有情况.
希望对各位有所启示.:victory:
页:
[1]