有限元中用Wilson-θ法或Newmark法计算非线性动力方程结果不收敛

weizelong

Wilson-θ法在求解刚性问题时，无论你是减小时间步长还是细化有限元网格，都无法得到收敛的结果的！

hew136

能否提供有限元法的Wilson-θ法程序，谢谢

ME!

1. function [d,v,a] = Wilson( K, M, C, f, d1, v1, dt, tend )
2. % 利用Wilson-theta 法计算结构的动力响应
3. % [d,v,a] = Wilson( K, M, C, f, d1, v1, dt, tend )
4. % 输入参数
5. % K ----- 刚度矩阵
6. % M ----- 质量矩阵
7. % C ----- 阻尼矩阵
8. % d1 ----- 初始位移
9. % v1 ----- 初始速度
10. % dt ----- 时间步长
11. % tend --- 结束时间
12. % 返回值
13. % d ----- 位移
14. % v ----- 速度
15. % a ----- 加速度
16. [n,n] = size( K ) ;
17. theta = 1.4 ;
18. tao = theta*dt ;
19. alpha0 = 6/tao^2 ;
20. alpha1 = 3/tao ;
21. alpha2 = 2*alpha1 ;
22. alpha3 = tao/2 ;
23. alpha4 = alpha0/theta ;
24. alpha5 = -alpha2/theta ;
25. alpha6 = 1-3/theta ;
26. alpha7 = dt/2 ;
27. alpha8 = dt^2/6 ;
28. K1 = K + alpha0*M + alpha1*C ;
29. d = zeros( n, floor(tend/dt) + 1 ) ;
30. v = zeros( n, floor(tend/dt) + 1 ) ;
31. a = zeros( n, floor(tend/dt) + 1 ) ;
32. d(:,1) = d1 ;
33. v(:,1) = v1 ;
34. a(:,1) = M\(f(:,1)-K*d1-C*v1) ;
35. for i=2:1:floor(tend/dt) + 1
36. t = (i-1)*dt ;
37. ftheta = floor(theta) ;
38. fq = f(i-1+ftheta-1)+ (theta-ftheta)*( f(i+ftheta-1) - f(i+ftheta-2) ) ;
39. f1 = fq + M*(alpha0*d(:,i-1)+alpha2*v(:,i-1)+2*a(:,i-1)) + C*(alpha1*d(:,i-1)+2*v(:,i-1)+alpha3*a(:,i-1)) ;
40. dq = K1\f1 ;
41. a(:,i) = alpha4*(dq-d(:,i-1)) + alpha5*v(:,i-1) + alpha6*a(:,i-1) ;
42. v(:,i) = v(:,i-1) + alpha7 * ( a(:,i) + a(:,i-1) ) ;
43. d(:,i) = d(:,i-1) + dt*v(:,i-1) + alpha8 * ( a(:,i)+2*a(:,i-1) ) ;
45. end
46. return

jiangpeng

有两个原因：
1，求解大变形问题时，若采用newmark法可能需要“预测-修正”，起码我的经验是采用了“预测-修正”会好一些的。
2，即使采用预测和修正了，还有一个更根本的问题：算法本身的缺陷，参见文献：
《Energy-conserving and decaying algorithms in non-linear structural dynamics》和《On a composite implicit time integration procedure for nonlinear dynamics》
我采用了bathe的方法后，求解壳体大变形的瞬态问题后，收敛性能提高，且精度不错

建议：
1，最简单的方法，减小时间步长，或者调整参数（gamma=0.6   beta=0.3025），缺点，计算量较大，或者精度不好
2，若还不行，采用预测-修正法，这个只需在newmark-beta法上做一点改动即可
3，最根本性的解决方法，采用G-alpha法，或者bathe的方法（文献中提到的）。但是以我的经验，这些方法也只是改善，没有彻底解决。