Wilson-θ 法解，响应无限大

ctovy1

在做悬壁梁简谐激振有限无分析时，分为10个单元，用Wilson-θ 法解响应，  得到无限大，

1. clear
   
2. syms ke L E I p A W
   
3. % N=[1-3*ke^2+2*ke^3,L*ke*(1-2*ke+ke^2),3*ke^2-2*ke^3,L*ke^2*(ke-1)];%形函数
   
4. kk=E*I/L^3*[12     6*L     -12     6*L;
   
5.              6*L    4*L^2   -6*L    2*L^2;
   
6.              -12    -6*L    12      -6*L;
   
7.              6*L    2*L^2   -6*L    4*L^2];%单元刚度阵
   
8. mm=p*A*L/420*[156 22*L 54 -13*L;
   
9.               22*L 4*L^2 13*L -3*L^2;
   
10.               54 13*L 156 -22*L;
    
11.               -13*L -3*L^2 -22*L 4*L^2];%单元质量阵
    
12. 
    
13. n=10;%有限单元数
    
14. KK=sym(zeros(2*(n+1),2*(n+1)));
    
15. MM=sym(zeros(2*(n+1),2*(n+1)));
    
16. 
    
17. %组成整体刚阵、质阵
    
18. for i=1:2:2*n
    
19.     KK(i:i+3,i:i+3)=KK(i:i+3,i:i+3)+kk;
    
20.     MM(i:i+3,i:i+3)=MM(i:i+3,i:i+3)+mm;
    
21. end
    
22. 
    
23. %去掉固定端约束
    
24. NK=KK(3:2*(n+1),3:2*(n+1));
    
25. NM=MM(3:2*(n+1),3:2*(n+1));
    
26. 
    
27. %材料常数
    
28. E=210e9;
    
29. b=0.05;
    
30. h=0.0085;
    
31. l=0.898;
    
32. I=b*h^3/12;
    
33. A=b*h;
    
34. L=l/n;
    
35. p=7800000;
    
36. 
    
37. %转代为数值矩阵
    
38. NK=eval(NK);
    
39. NM=eval(NM);
    
40. syms t
    
41. omg=10;
    
42. %初始条件
    
43. C=zeros(2*n,2*n);
    
44. F=sym(zeros(2*n,1));
    
45. x0=zeros(2*n,1);
    
46. F(2*n-1)=10*sin(10*t);
    
47. dx0=zeros(2*n,1);
    
48. % wilson theta 法
    
49. [ts,Xt,DXt]=WILS(NM,NK,C,F,x0,dx0);
    
50. %
    
51. m=length(x0);
    
52. for ii=m-1
    
53.     tt=(ii+1)/2;
    
54.     figure(1)
    
55. %     subplot(2,1,1)
    
56.     plot(ts,Xt(ii,:))
    
57.     title(['节点' num2str(tt) '振动曲线'])
    
58.     text(ts(1),Xt(ii,1),num2str(tt))
    
59.     xlabel('时间t');
    
60.     ylabel('位移');
    
61.     hold on
    
62.    
    
63.     figure(2)
    
64.     subplot(2,1,2)
    
65.     plot(ts,DXt(ii,:))
    
66.     title(['节点' num2str(tt) '速度曲线'])
    
67.     xlabel('时间t');
    
68.     ylabel('速度');
    
69.     hold on
    
70. end
    
71. 
    
72. 下面是wilsontheta法函数的代码
    
73. function [ts,Xt,DXt]=WILS(M,K,C,F,x0,dx0)
    
74. dt=0.09;
    
75. n=500;
    
76. theta=1.4;
    
77. 
    
78. a0=6/(theta*dt)^2;
    
79. a1=3/(theta*dt);
    
80. a2=2*a1;
    
81. a3=theta*dt/2;
    
82. a4=a0/theta;
    
83. a5=-a2/theta;
    
84. a6=1-3/theta;
    
85. a7=dt/2;
    
86. a8=dt^2/6;
    
87. 
    
88. %初始条件
    
89. d2x0=inv(M)*(F-K*x0-C*dx0);
    
90. %等效刚度
    
91. Kv=K+a0*M+a1*C;   
    
92. 
    
93. t=0;
    
94. Xt(:,1)=x0;
    
95. DXt(:,1)=dx0;
    
96. D2Xt(:,1)=eval(d2x0);
    
97. Ft(:,1)=eval(F);
    
98. i=1;
    
99. ts(i)=0;
    
100. for t=dt:dt:n*dt
     
101.    i=i+1;
     
102.    ts(i)=t;
     
103.    Fdt=eval(F);
     
104.    
     
105.    Fv=Ft(:,i-1)+theta*(Fdt-Ft(:,i-1))+M*(a0*Xt(:,i-1)+a2*DXt(:,i-1)+...
     
106.        2*D2Xt(:,i-1))+C*(a1*Xt(:,i-1)+2*DXt(:,i-1)+a3*D2Xt(:,i-1));
     
107.    Xtheta=inv(Kv)*Fv;
     
108.    
     
109.    D2X=a4*(Xtheta-Xt(:,i-1))+a5*DXt(:,i-1)+a6*D2Xt(:,i-1);
     
110.    DX=DXt(:,i-1)+a7*(D2X+D2Xt(:,i-1));
     
111.    X=Xt(:,i-1)+dt*DXt(:,i-1)+a8*(D2X+2*D2Xt(:,i-1));
     
112.    
     
113.     Xt(:,i)=X;
     
114.     DXt(:,i)=DX;
     
115.     D2Xt(:,i)=inv(M)*(Fdt-K*Xt(:,i)-C*DXt(:,i));
     
116.     Ft(:,i)=Fdt;
     
117. end

复制代码

ctovy1

已经解决

luchun

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

Catsayer

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

Triste

咋解决的啊？求指导