迭代法求变截面欧拉梁的固有特性(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 的帖子
呃……我写的这些个数值程序可都是200-300行的……有限元就更长了…… 有限元也没那么长吧!
当然你不能把矩阵运算的代码算进去! 一般有限元一次元是比较好的!我都用一次元,在一次元基础上用外推和插值,结果会比高次元好很多的,当然比如是那些奇异性问题效果不见得好,但是特殊问题,高次也没啥用啊!
个人观点 回复 smtmobly 的帖子
当然要算矩阵运算,一般是变半带宽压缩,还有方程组的求解也要算……
老师说了,不许用MATLAB,要自己写!
{:{17}:} 是自己写,方程组的求解,有标准算法啊!你用什么 下降法和超松弛就足够了吧 回复 smtmobly 的帖子
对,我用的是列主元高斯消去法和超松弛迭代 一次元最简单,而且光滑性要求也较低,只比常元高。
关于有限元、差分等算法的外推可以看《分裂外推与组合技巧》吕涛的书。
外推和组合是基于模型的一种后处理方法,这本书里说的比较详细。我的所有的外推上的
知识都是这本书里的!效果很好。 回复 smtmobly 的帖子
我一般也倾向于线性单元,但是会构造简单的等效模型换用高次的单元试算一下,确认差别不大再放心使用线性单元。谢谢你的指点,这是我没有听说和涉及过的方面,谢谢你! 回复 Rainyboy 的帖子
客气了!高次的效果好,说明光滑性好,那么低次的外推与组合效果一定很好的。
现在很多算法,比如什么多尺度在一定程度上是一种组合或者外推,组合和外推的差别有时候很难区分。(形式上),理论上就完全不一样。
关键还是看你计算的模型是什么类型的。
考虑基于模型的这种组合方法,会让你在模型层面上去考虑,使用有限元还是差分,或者在某些区域使用有限元,另一个区域使用差分。甚至边界元、谱方法等。比您点评的那个四阶方程的例子实际上就是谱方法的一个特例。
基于模型的考虑,思路会更宽阔,并且将目光放在主流算法的基本点上,不至于陷入到一些特别艰深的问题中去 回复 smtmobly 的帖子
说得是,也许我常见的工作是迅速实现一种计算方法(这种方法只要在计算空间上可行、时间上可接受就差不多了),然后就开始调整参数,分析对比结果,算个寿命算个应力什么的,反而对方法本身没有那么敏感了……呵呵 有效就好哈!其他就当做玩了 不错的东西啊,帮你顶啊...
页:
[1]
2