ctovy1 发表于 2011-6-25 15:10

Wilson-θ 法解，响应无限大

在做悬壁梁简谐激振有限无分析时，分为10个单元，用Wilson-θ 法解响应，得到无限大，clear
syms ke L E I p A W
% N=;%形函数
kk=E*I/L^3*[12   6*L   -12   6*L;
             6*L    4*L^2   -6*L    2*L^2;
             -12    -6*L    12      -6*L;
             6*L    2*L^2   -6*L    4*L^2];%单元刚度阵
mm=p*A*L/420*[156 22*L 54 -13*L;
            22*L 4*L^2 13*L -3*L^2;
            54 13*L 156 -22*L;
            -13*L -3*L^2 -22*L 4*L^2];%单元质量阵

n=10;%有限单元数
KK=sym(zeros(2*(n+1),2*(n+1)));
MM=sym(zeros(2*(n+1),2*(n+1)));

%组成整体刚阵、质阵
for i=1:2:2*n
    KK(i:i+3,i:i+3)=KK(i:i+3,i:i+3)+kk;
    MM(i:i+3,i:i+3)=MM(i:i+3,i:i+3)+mm;
end

%去掉固定端约束
NK=KK(3:2*(n+1),3:2*(n+1));
NM=MM(3:2*(n+1),3:2*(n+1));

%材料常数
E=210e9;
b=0.05;
h=0.0085;
l=0.898;
I=b*h^3/12;
A=b*h;
L=l/n;
p=7800000;

%转代为数值矩阵
NK=eval(NK);
NM=eval(NM);
syms t
omg=10;
%初始条件
C=zeros(2*n,2*n);
F=sym(zeros(2*n,1));
x0=zeros(2*n,1);
F(2*n-1)=10*sin(10*t);
dx0=zeros(2*n,1);
% wilson theta 法
=WILS(NM,NK,C,F,x0,dx0);
%
m=length(x0);
for ii=m-1
    tt=(ii+1)/2;
    figure(1)
%   subplot(2,1,1)
    plot(ts,Xt(ii,:))
    title(['节点' num2str(tt) '振动曲线'])
    text(ts(1),Xt(ii,1),num2str(tt))
    xlabel('时间t');
    ylabel('位移');
    hold on
   
    figure(2)
    subplot(2,1,2)
    plot(ts,DXt(ii,:))
    title(['节点' num2str(tt) '速度曲线'])
    xlabel('时间t');
    ylabel('速度');
    hold on
end

下面是wilsontheta法函数的代码
function =WILS(M,K,C,F,x0,dx0)
dt=0.09;
n=500;
theta=1.4;

a0=6/(theta*dt)^2;
a1=3/(theta*dt);
a2=2*a1;
a3=theta*dt/2;
a4=a0/theta;
a5=-a2/theta;
a6=1-3/theta;
a7=dt/2;
a8=dt^2/6;

%初始条件
d2x0=inv(M)*(F-K*x0-C*dx0);
%等效刚度
Kv=K+a0*M+a1*C;   

t=0;
Xt(:,1)=x0;
DXt(:,1)=dx0;
D2Xt(:,1)=eval(d2x0);
Ft(:,1)=eval(F);
i=1;
ts(i)=0;
for t=dt:dt:n*dt
   i=i+1;
   ts(i)=t;
   Fdt=eval(F);
   
   Fv=Ft(:,i-1)+theta*(Fdt-Ft(:,i-1))+M*(a0*Xt(:,i-1)+a2*DXt(:,i-1)+...
       2*D2Xt(:,i-1))+C*(a1*Xt(:,i-1)+2*DXt(:,i-1)+a3*D2Xt(:,i-1));
   Xtheta=inv(Kv)*Fv;
   
   D2X=a4*(Xtheta-Xt(:,i-1))+a5*DXt(:,i-1)+a6*D2Xt(:,i-1);
   DX=DXt(:,i-1)+a7*(D2X+D2Xt(:,i-1));
   X=Xt(:,i-1)+dt*DXt(:,i-1)+a8*(D2X+2*D2Xt(:,i-1));
   
    Xt(:,i)=X;
    DXt(:,i)=DX;
    D2Xt(:,i)=inv(M)*(Fdt-K*Xt(:,i)-C*DXt(:,i));
    Ft(:,i)=Fdt;
end

ctovy1 发表于 2011-6-30 02:25

已经解决

luchun 发表于 2016-9-5 19:25

你是怎么解决的啊，我现在用wilson算法也是，出现无限大，检查不出问题啊，急急急

Catsayer 发表于 2016-9-6 12:56

怎么解决的？？？？？楼主吊胃口啊

Triste 发表于 2016-9-7 09:01

咋解决的啊？求指导
页: [1]
查看完整版本: Wilson-θ 法解，响应无限大