迭代法求变截面欧拉梁的固有特性(In C/C++)

Rainyboy

一种经典的结构动力特性数值分析方法。


总的来说，就是首先通过理论推导构造出递推式：

{实际振型} = f({实际振型});//这里，f(x)实际上是一个四次积分。

然后假设一个振型，带入递推式，具体为：

1. 
   
2. 任意给定{假设振型};
   
3. 
   
4. while( abs({修正振型} - {假设振型}) < 1e-5)
   
5. {
   
6.     {修正振型} = f({假设振型});
   
7. 
   
8.     {假设振型} = {修正振型} ;
   
9. }

复制代码

递推式是根据欧拉梁横向振动方程得到的，以公式居多，就作为附件吧（附件[理论基础.doc]即是）。

某变截面梁的前三阶振型：



还是挺像回事儿的哈。

PS:附件[计算结果.doc]是根据本文代码中硬编码的梁模型计算得到的结果。

源代码是：

1. 
   
2. // VibrateStyle.cpp : 定义控制台应用程序的入口点。
   
3. // 结构·第三·振型计算
   
4. #include "stdafx.h"
   
5. #include <iostream>
   
6. #include <iomanip>
   
7. #include <fstream>
   
8. #include <math.h>
   
9. /********************宏，命名空间*****************/
   
10. #define MAXNUM 11
    
11. #define STOPLINE 0.000001
    
12. #define SHOW_DETAIL
    
13. #define OUTPUT(K)  cout<<"****结果****"<<endl;\
    
14.   cout<<K<<"阶共振 "<<w##K<<endl;for(int i=0;i<MAXNUM;++i){cout<<Shape_##K[i]<<endl;}
    
15. using namespace std;
    
16. /***********************常量表********************/
    
17. //密度
    
18. const double Dens=2850;
    
19. //弹性模量
    
20. const double E=71.54e9;
    
21. //截面惯性矩
    
22. const double I[MAXNUM]=
    
23. {
    
24.     0.0279e-8,0.0212e-8,0.0157e-8,0.0108e-8,0.0084e-8,
    
25.     0.0061e-8,0.0045e-8,0.0037e-8,0.0032e-8,0.0030e-8,
    
26.     0.0030e-8
    
27. };
    
28. //截面面积
    
29. const double A[MAXNUM]=
    
30. {
    
31.     1.70e-4,1.46e-4,1.26e-4,1.09e-4,0.96e-4,
    
32.     0.86e-4,0.77e-4,0.73e-4,0.70e-4,0.68e-4,
    
33.     0.68e-4
    
34. };
    
35. //截面相对积分步长
    
36. const double dx=0.01;
    
37. //共振频率
    
38. double w1,w2,w3;
    
39. //各阶振型数据
    
40. double Shape_1[MAXNUM];
    
41. double Shape_2[MAXNUM];
    
42. double Shape_3[MAXNUM];
    
43. /********************函数声明区*******************/
    
44. bool IfStop(double PreY[MAXNUM],double NextY[MAXNUM]);//收敛判据
    
45. double Integral(int k,double PreY[MAXNUM]);//四重积分
    
46. void Comp(void (*f_Adjust)(double NextY[MAXNUM]),double Shape[MAXNUM],double &w);//计算振型的函数
    
47. void Adjust_2ndorder(double NextY[MAXNUM]);//二阶振型的调整函数
    
48. void Adjust_3rdorder(double NextY[MAXNUM]);//三阶振型的调整函数
    
49. int _tmain(int argc, _TCHAR* argv[])
    
50. {
    
51.    //打开文件等等
    
52.    ofstream OutPutFile;
    
53.    OutPutFile.open("re.plt",ios::out);
    
54.    /*********计算并输出到屏幕***********/
    
55.    Comp(NULL,Shape_1,w1);OUTPUT(1)
    
56.    Comp(Adjust_2ndorder,Shape_2,w2);OUTPUT(2)
    
57.    Comp(Adjust_3rdorder,Shape_3,w3);OUTPUT(3)
    
58.    /**********结果记录到文件************/
    
59.    OutPutFile<<"第一阶"<<endl;
    
60.    for(int i=0;i<MAXNUM;++i)
    
61.        OutPutFile<<0.1*i<<"  "<<Shape_1[i]<<endl;
    
62.    OutPutFile<<"第二阶"<<endl;
    
63.    for(int i=0;i<MAXNUM;++i)
    
64.       OutPutFile<<0.1*i<<"  "<<Shape_2[i]<<endl;
    
65.    OutPutFile<<"第三阶"<<endl;
    
66.    for(int i=0;i<MAXNUM;++i)
    
67.       OutPutFile<<0.1*i<<"  "<<Shape_3[i]<<endl;
    
68.    system("pause");
    
69.    OutPutFile.close();
    
70.    return 0;
    
71. }
    
72. void Comp(void (*f_Adjust)(double NextY[MAXNUM]),double Shape[MAXNUM],double &w)
    
73. {
    
74.    //初始振型
    
75.    double PreY[MAXNUM]={0.1,0.1,0.1,0.1,0.1,0,1};
    
76.    //迭代振型
    
77.    double NextY[MAXNUM]={0.0};
    
78.    while(1)
    
79.    {
    
80.         #ifdef SHOW_DETAIL
    
81.         cout<<"****Test****"<<endl;
    
82.         for(int i=0;i<MAXNUM;++i)
    
83.               cout<<PreY[i]<<endl;
    
84.         system("pause");
    
85.         system("cls");
    
86.         #endif
    
87.         //计算四重积分
    
88.         for(int i=0;i<MAXNUM;++i)
    
89.         NextY[i]=Integral(i,PreY);
    
90.         //归一化
    
91.         for(int i=0;i<MAXNUM;++i)
    
92.               NextY[i]=NextY[i]/NextY[MAXNUM-1];
    
93.         //根据正交条件调整
    
94.         if(f_Adjust)
    
95.               f_Adjust(NextY);
    
96.         //判断收敛
    
97.         if(!IfStop(PreY,NextY))
    
98.              break;
    
99.         //迭代
    
100.         for(int i=0;i<MAXNUM;++i)
     
101.              PreY[i]=NextY[i];
     
102.    }
     
103.    //记录当前收敛解
     
104.    for(int i=0;i<MAXNUM;++i)
     
105.         PreY[i]=NextY[i];
     
106.    //计算共振频率
     
107.    w=pow(E/Dens/Integral(MAXNUM-1,PreY),0.5);
     
108.    //记录到全局变量
     
109.    for(int i=0;i<MAXNUM;++i)
     
110.         Shape[i]=NextY[i];
     
111. }
     
112. bool IfStop(double PreY[MAXNUM],double NextY[MAXNUM])
     
113. {
     
114.    double Sum=0;
     
115.    for(int i=0;i<MAXNUM;++i)
     
116.        Sum+=pow(PreY[i]-NextY[i],2);
     
117.    if(pow(Sum,0.5)<STOPLINE)
     
118.        return false;//停止迭代
     
119.    return true;//继续迭代
     
120. }
     
121. double Integral(int k,double PreY[MAXNUM])
     
122. {
     
123.     int i,j,p,q;
     
124.     double Sum_1=0,Sum_2=0,Sum_3=0,Sum_4=0;
     
125.     i=j=p=q=1;
     
126.     for(q=1;q<=k;++q)
     
127.     {
     
128.         Sum_2=0;
     
129.         for(p=1;p<=q;++p)
     
130.         {
     
131.               Sum_3=0;
     
132.               for(j=p;j<MAXNUM;++j)
     
133.               {
     
134.                      Sum_4=0;
     
135.                      for(i=j;i<MAXNUM;++i)
     
136.                      {
     
137.                             Sum_4+=A[i]*PreY[i]*dx;
     
138.                       }
     
139.                       Sum_3+=Sum_4*dx;
     
140.               }
     
141.               Sum_2+=Sum_3/I[p]*dx;
     
142.         }
     
143.         Sum_1+=Sum_2*dx;
     
144.    }
     
145.    return Sum_1;
     
146. }
     
147. void Adjust_2ndorder(double NextY[MAXNUM])
     
148. {
     
149.        double C12=0,B11=0,a1;
     
150.        for(int i=1;i<MAXNUM;++i)
     
151.        {
     
152.               C12+=A[i]*NextY[i]*Shape_1[i];
     
153.               B11+=A[i]*Shape_1[i]*Shape_1[i];
     
154.        }
     
155.        a1=C12/B11;
     
156.        for(int i=0;i<MAXNUM;++i)
     
157.        {
     
158.               NextY[i]=NextY[i]-a1*Shape_1[i];
     
159.        }
     
160.        return;
     
161. }
     
162. void Adjust_3rdorder(double NextY[MAXNUM])
     
163. {
     
164.        double C13=0,C23=0,B11=0,B22=0,a1,a2;
     
165.        for(int i=1;i<MAXNUM;++i)
     
166.        {
     
167.        C13+=A[i]*NextY[i]*Shape_1[i];
     
168.        C23+=A[i]*NextY[i]*Shape_2[i];
     
169.        B11+=A[i]*Shape_1[i]*Shape_1[i];
     
170.        B22+=A[i]*Shape_2[i]*Shape_2[i];
     
171.        }
     
172.        a1=C13/B11;
     
173.        a2=C23/B22;
     
174.        for(int i=0;i<MAXNUM;++i)
     
175.        {
     
176.               NextY[i]=NextY[i]-a1*Shape_1[i]-a2*Shape_2[i];
     
177.        }
     
178.        return;
     
179. }

复制代码

smtmobly

做数值计算也好几年了，怎么除了我毕业论文写了一个通用的计算椭圆、抛物、双曲方程的那个程序超过百行，其他的就纯数值处理上从没有超过百行的。

Rainyboy

回复 smtmobly 的帖子

呃……我写的这些个数值程序可都是200-300行的……有限元就更长了……

smtmobly

有限元也没那么长吧！
当然你不能把矩阵运算的代码算进去！

smtmobly

一般有限元一次元是比较好的！我都用一次元，在一次元基础上用外推和插值，结果会比高次元好很多的，当然比如是那些奇异性问题效果不见得好，但是特殊问题，高次也没啥用啊！
个人观点

Rainyboy

回复 smtmobly 的帖子

当然要算矩阵运算，一般是变半带宽压缩，还有方程组的求解也要算……
老师说了，不许用MATLAB，要自己写！

smtmobly

是自己写，方程组的求解，有标准算法啊！你用什么

smtmobly

下降法和超松弛就足够了吧

Rainyboy

回复 smtmobly 的帖子

对，我用的是列主元高斯消去法和超松弛迭代

smtmobly

一次元最简单，而且光滑性要求也较低，只比常元高。
关于有限元、差分等算法的外推可以看《分裂外推与组合技巧》吕涛的书。
外推和组合是基于模型的一种后处理方法，这本书里说的比较详细。我的所有的外推上的
知识都是这本书里的！效果很好。

Rainyboy

回复 smtmobly 的帖子

我一般也倾向于线性单元，但是会构造简单的等效模型换用高次的单元试算一下，确认差别不大再放心使用线性单元。谢谢你的指点，这是我没有听说和涉及过的方面，谢谢你！

smtmobly

回复 Rainyboy 的帖子

客气了！高次的效果好，说明光滑性好，那么低次的外推与组合效果一定很好的。
现在很多算法，比如什么多尺度在一定程度上是一种组合或者外推，组合和外推的差别有时候很难区分。（形式上），理论上就完全不一样。
关键还是看你计算的模型是什么类型的。
考虑基于模型的这种组合方法，会让你在模型层面上去考虑，使用有限元还是差分，或者在某些区域使用有限元，另一个区域使用差分。甚至边界元、谱方法等。比您点评的那个四阶方程的例子实际上就是谱方法的一个特例。
基于模型的考虑，思路会更宽阔，并且将目光放在主流算法的基本点上，不至于陷入到一些特别艰深的问题中去

Rainyboy

回复 smtmobly 的帖子

说得是，也许我常见的工作是迅速实现一种计算方法（这种方法只要在计算空间上可行、时间上可接受就差不多了），然后就开始调整参数，分析对比结果，算个寿命算个应力什么的，反而对方法本身没有那么敏感了……呵呵

smtmobly

有效就好哈！其他就当做玩了

有梦的人

不错的东西啊，帮你顶啊...