wilson.m

薄壁结构有限元计算

ct=1.37;
for i=1:500
    kx=K+3*C/ct/deltat+6*M/ct/ct/deltat/deltat;
    dpx=-M*unit*xgt(i)+M*(6*v1/ct/deltat+3*a1)+C*(3*v1+ct*a1*deltat/2);
    dxt=kx\dpx;
    da=6*dxt/ct/ct/ct/deltat/deltat-6*v1/ct/ct/deltat-3*a1/ct;
    a=a1+da;
    dv=deltat*(a+a1)/2;
    v=v1+dv;
    dx=deltat*v1+deltat*deltat*a1/3+deltat*deltat*a/6;
    x=x1+dx;
    a=-M\(M*unit*xg2(i)+C*v+K*x);
    x1=x;
    v1=v;
    a1=a;
    xt=[xt x];
    vt=[vt v];
    at=[at a];
end