无阻尼 自由振动 暂态问题

godfane

利用有限元和直接积分求解一暂态问题(一维波动--无阻尼自由弦振动),推导了有限元的单元公式,见附图(本人推导,希望被指出错误).
在集成了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=[K11 K12;K12 K11];          %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=[M11 M12;M12 M22];  %计算每个单元的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

恳请各位指正!

无水1324

看看！
你的matlab是哪个版本的，在6.5中出现问题？

godfane

我的matlab是7.0 问题是程序可以运行,可是把算出来的结果做了animation,明显是不对的.
反复检查又找不到错在哪里.:@L

godfane

自问自答:
原来问题是出在形函数积分的计算上,如果按原程序中的代码,则每个单元的对shape function的积分都是相等了,实际上是不同的.
由于一开始试算的时候,最从最简单的划分两个单元开始,这种情况下对N1和N2的积分是相等的,但不是所有情况.
希望对各位有所启示.:victory: