迭代法求变截面欧拉梁的固有特性(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