Rainyboy 发表于 2010-10-28 15:28

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

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


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

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

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


任意给定{假设振型};

while( abs({修正振型} - {假设振型}) < 1e-5)
{
    {修正振型} = f({假设振型});

    {假设振型} = {修正振型} ;
}
递推式是根据欧拉梁横向振动方程得到的,以公式居多,就作为附件吧(附件[理论基础.doc]即是)。

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



还是挺像回事儿的哈。

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

源代码是:


// VibrateStyle.cpp : 定义控制台应用程序的入口点。
// 结构·第三·振型计算
#include "stdafx.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <math.h>
/********************宏,命名空间*****************/
#define MAXNUM 11
#define STOPLINE 0.000001
#define SHOW_DETAIL
#define OUTPUT(K)cout<<"****结果****"<<endl;\
cout<<K<<"阶共振 "<<w##K<<endl;for(int i=0;i<MAXNUM;++i){cout<<Shape_##K<<endl;}
using namespace std;
/***********************常量表********************/
//密度
const double Dens=2850;
//弹性模量
const double E=71.54e9;
//截面惯性矩
const double I=
{
    0.0279e-8,0.0212e-8,0.0157e-8,0.0108e-8,0.0084e-8,
    0.0061e-8,0.0045e-8,0.0037e-8,0.0032e-8,0.0030e-8,
    0.0030e-8
};
//截面面积
const double A=
{
    1.70e-4,1.46e-4,1.26e-4,1.09e-4,0.96e-4,
    0.86e-4,0.77e-4,0.73e-4,0.70e-4,0.68e-4,
    0.68e-4
};
//截面相对积分步长
const double dx=0.01;
//共振频率
double w1,w2,w3;
//各阶振型数据
double Shape_1;
double Shape_2;
double Shape_3;
/********************函数声明区*******************/
bool IfStop(double PreY,double NextY);//收敛判据
double Integral(int k,double PreY);//四重积分
void Comp(void (*f_Adjust)(double NextY),double Shape,double &w);//计算振型的函数
void Adjust_2ndorder(double NextY);//二阶振型的调整函数
void Adjust_3rdorder(double NextY);//三阶振型的调整函数
int _tmain(int argc, _TCHAR* argv[])
{
   //打开文件等等
   ofstream OutPutFile;
   OutPutFile.open("re.plt",ios::out);
   /*********计算并输出到屏幕***********/
   Comp(NULL,Shape_1,w1);OUTPUT(1)
   Comp(Adjust_2ndorder,Shape_2,w2);OUTPUT(2)
   Comp(Adjust_3rdorder,Shape_3,w3);OUTPUT(3)
   /**********结果记录到文件************/
   OutPutFile<<"第一阶"<<endl;
   for(int i=0;i<MAXNUM;++i)
       OutPutFile<<0.1*i<<""<<Shape_1<<endl;
   OutPutFile<<"第二阶"<<endl;
   for(int i=0;i<MAXNUM;++i)
      OutPutFile<<0.1*i<<""<<Shape_2<<endl;
   OutPutFile<<"第三阶"<<endl;
   for(int i=0;i<MAXNUM;++i)
      OutPutFile<<0.1*i<<""<<Shape_3<<endl;
   system("pause");
   OutPutFile.close();
   return 0;
}
void Comp(void (*f_Adjust)(double NextY),double Shape,double &w)
{
   //初始振型
   double PreY={0.1,0.1,0.1,0.1,0.1,0,1};
   //迭代振型
   double NextY={0.0};
   while(1)
   {
      #ifdef SHOW_DETAIL
      cout<<"****Test****"<<endl;
      for(int i=0;i<MAXNUM;++i)
            cout<<PreY<<endl;
      system("pause");
      system("cls");
      #endif
      //计算四重积分
      for(int i=0;i<MAXNUM;++i)
      NextY=Integral(i,PreY);
      //归一化
      for(int i=0;i<MAXNUM;++i)
            NextY=NextY/NextY;
      //根据正交条件调整
      if(f_Adjust)
            f_Adjust(NextY);
      //判断收敛
      if(!IfStop(PreY,NextY))
             break;
      //迭代
      for(int i=0;i<MAXNUM;++i)
             PreY=NextY;
   }
   //记录当前收敛解
   for(int i=0;i<MAXNUM;++i)
      PreY=NextY;
   //计算共振频率
   w=pow(E/Dens/Integral(MAXNUM-1,PreY),0.5);
   //记录到全局变量
   for(int i=0;i<MAXNUM;++i)
      Shape=NextY;
}
bool IfStop(double PreY,double NextY)
{
   double Sum=0;
   for(int i=0;i<MAXNUM;++i)
       Sum+=pow(PreY-NextY,2);
   if(pow(Sum,0.5)<STOPLINE)
       return false;//停止迭代
   return true;//继续迭代
}
double Integral(int k,double PreY)
{
    int i,j,p,q;
    double Sum_1=0,Sum_2=0,Sum_3=0,Sum_4=0;
    i=j=p=q=1;
    for(q=1;q<=k;++q)
    {
      Sum_2=0;
      for(p=1;p<=q;++p)
      {
            Sum_3=0;
            for(j=p;j<MAXNUM;++j)
            {
                     Sum_4=0;
                     for(i=j;i<MAXNUM;++i)
                     {
                            Sum_4+=A*PreY*dx;
                      }
                      Sum_3+=Sum_4*dx;
            }
            Sum_2+=Sum_3/I*dx;
      }
      Sum_1+=Sum_2*dx;
   }
   return Sum_1;
}
void Adjust_2ndorder(double NextY)
{
       double C12=0,B11=0,a1;
       for(int i=1;i<MAXNUM;++i)
       {
            C12+=A*NextY*Shape_1;
            B11+=A*Shape_1*Shape_1;
       }
       a1=C12/B11;
       for(int i=0;i<MAXNUM;++i)
       {
            NextY=NextY-a1*Shape_1;
       }
       return;
}
void Adjust_3rdorder(double NextY)
{
       double C13=0,C23=0,B11=0,B22=0,a1,a2;
       for(int i=1;i<MAXNUM;++i)
       {
       C13+=A*NextY*Shape_1;
       C23+=A*NextY*Shape_2;
       B11+=A*Shape_1*Shape_1;
       B22+=A*Shape_2*Shape_2;
       }
       a1=C13/B11;
       a2=C23/B22;
       for(int i=0;i<MAXNUM;++i)
       {
            NextY=NextY-a1*Shape_1-a2*Shape_2;
       }
       return;
}

smtmobly 发表于 2010-10-29 09:04

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

Rainyboy 发表于 2010-10-29 09:14

回复 smtmobly 的帖子

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

smtmobly 发表于 2010-10-29 09:35

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

smtmobly 发表于 2010-10-29 09:38

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

Rainyboy 发表于 2010-10-29 09:39

回复 smtmobly 的帖子

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

smtmobly 发表于 2010-10-29 09:41

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

smtmobly 发表于 2010-10-29 09:42

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

Rainyboy 发表于 2010-10-29 09:52

回复 smtmobly 的帖子

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

smtmobly 发表于 2010-10-29 12:20

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

Rainyboy 发表于 2010-10-29 12:25

回复 smtmobly 的帖子

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

smtmobly 发表于 2010-10-29 12:35

回复 Rainyboy 的帖子

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

Rainyboy 发表于 2010-10-29 12:51

回复 smtmobly 的帖子

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

smtmobly 发表于 2010-10-29 12:55

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

有梦的人 发表于 2010-11-15 17:00

不错的东西啊,帮你顶啊...
页: [1] 2
查看完整版本: 迭代法求变截面欧拉梁的固有特性(In C/C++)