5
【旧代码】传热过程数值模拟
source link: https://byronhe.com/post/2011/11/15/heat-grid/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
【旧代码】传热过程数值模拟
2011-11-15
传热过程数值模拟（《传热学》实验指导书第四部分第一题，第一，第二类边界条件）
2010年十月写的旧代码。
第一类边界条件是给定边界温度。
第二类是对流边界。
区域都是如下形状的:
--------------------------------
| |
| |
| |
| ----------------------
| |
| |
| |
| |
| |
| |
| |
| |
---------
用C++纯属蛋疼。
第一类边界条件：
/*
* 等温边界
*/
#include<iostream>
#include<cmath>
using namespace std;
const double out_temp=30.0;//外边界温度
const double in_temp=0.0;//内边界温度
const double accuracy=0.00000001;//精度
const double lambda=0.53;//导热系数
const int width=15;//上部点阵宽度
const int width_bottom=4;//下部点阵宽度
const int height=4;//上部点阵高度
const int height_bottom=7;//下部点阵高度
/*
//另一组参数
const int width=16;
const int width_bottom=6;
const int height=6;
const int height_bottom=6;
*/
//总共点数
const int num_of_points=width*height+width_bottom*height_bottom;
//单个点
class point
{
public:
double temp;//温度
int up;//数组下标
int down;
int left;
int right;
point(){
temp=0.0;//初始化成0摄氏度
up=down=right=left=0;
}
};
ostream & operator<<(ostream & src_stream,const point & src){
src_stream<<"temp="<<src.temp;
src_stream<<" up="<<src.up<<" down="<<src.down;
src_stream<<" left="<<src.left<<" right="<<src.right;
return src_stream;
}
void print_grid(point * points){//输出网格各点的温度
cout<<endl;
for(int position=0;position>=0;position=points[position].down){
for(int tmp=position;tmp>=0;tmp=points[tmp].right){
cout.width(10);
cout<<points[tmp].temp;
}
cout<<endl<<endl<<endl<<endl;
}
}
double calc_direction(point * points,int direction,double &opposit_weight){
//计算给定方向邻点的温度和反方向的权重。
switch (direction) {
case -1:
return out_temp;
break;
case -2:
return in_temp;
break;
case -3:
opposit_weight*=2;
return 0.0;
break;
default:
return points[direction].temp;
}
}
void compu_point(point * points,int now){
//根据周围四个点算出指定点的温度
double left_temp;
double right_temp;
double up_temp;
double down_temp;
double up_wight=0.25;//上方权重,默认0.25
double down_wight=0.25;
double left_wight=0.25;
double right_wight=0.25;
left_temp=calc_direction(points,points[now].left,right_wight);
right_temp=calc_direction(points,points[now].right,left_wight);
up_temp=calc_direction(points,points[now].up,down_wight);
down_temp=calc_direction(points,points[now].down,up_wight);
points[now].temp=left_temp*left_wight+right_temp*right_wight+up_temp*up_wight+down_temp*down_wight;
}
void rec_walk(point * points){
// "左=>右"里嵌"上=>下"遍历计算
for(int position=0;position>=0;position=points[position].right)
for(int tmp=position;tmp>=0;tmp=points[tmp].down)
compu_point(points,tmp);
}
void heat_rate(point * points){
//计算导热量
double out_sum=0.0;
double in_sum=0.0;
//外边界
int pos=0;
double tmp=0.0;
for(;points[pos].right>=0;pos=points[pos].right){
tmp=out_temp - points[pos].temp ;
if(points[pos].right<=0)
out_sum += 0.5*tmp;
else
out_sum+=tmp;
}
for(pos=points[0].down;points[pos].down>=0;pos=points[pos].down){
tmp=out_temp - points[pos].temp ;
if(points[pos].down<=0)
out_sum += 0.5*tmp;
else
out_sum+=tmp;
}
//内边界
for(pos=num_of_points-1;points[pos].right<0;pos=points[pos].up){
tmp = points[pos].temp - in_temp ;
if(points[pos].down<=0)
in_sum += 0.5*tmp;
else
in_sum +=tmp;
}
for(pos=points[pos].right;points[pos].right>=0;pos=points[pos].right){
tmp = points[pos].temp - in_temp;
if(points[pos].right<=0)
in_sum += 0.5*tmp;
else
in_sum += tmp;
}
out_sum*=lambda;
in_sum*=lambda;
cout<<"外边界导热量="<<out_sum;
cout<<"\n内边界导热量="<<in_sum;
cout<<"\n误差="<<(in_sum-out_sum)/in_sum*100<<"%"<<endl;
}
void init_grid(point * points){
//初始化网格
int position=0;
//初始化最上面的边
for(int i=0;i<width;i++){
points[i].up=-1;//-1表示外边界，-2内边界,-3表示绝热
points[i].left=i-1;
points[i].right=i+1;
points[i].down=i+width;
position=i;
}
points[0].left=-1;//等温
points[width-1].right=-3;//绝热
//上面区域
for(int i=1;i<height;i++){
for(int j=0;j<width;j++){
points[i*width+j].up =i*width+j-width;
points[i*width+j].left =i*width+j-1;
points[i*width+j].right=i*width+j+1;
points[i*width+j].down =i*width+j+width;
position=i*width+j;
}
points[i*width+0].left=-1;
points[i*width+width-1].right=-3;
}
position-=width;
position++;
for(int j=width_bottom;j<width;j++){
points[position+j].down =-2;//内边界
}
position+=width;
//下面区域
for(int i=0;i<height_bottom;i++){
for(int j=0;j<width_bottom;j++){
if(i)
points[position+i*width_bottom+j].up=position+i*width_bottom+j-width_bottom;
else
points[position+i*width_bottom+j].up=position+i*width_bottom+j-width;
points[position+i*width_bottom+j].left=position+i*width_bottom+j-1;
points[position+i*width_bottom+j].right=position+i*width_bottom+j+1;
points[position+i*width_bottom+j].down=position+i*width_bottom+j+width_bottom;
}
points[position+i*width_bottom+0].left=-1;
points[position+i*width_bottom+width_bottom-1].right=-2;
}
//下边界
position+=width_bottom*height_bottom;
position-=width_bottom;
for(int j=0;j<width_bottom;j++){
points[position+j].down=-3;
}
position+=width_bottom;
}
int main(){
point * points=new point[num_of_points];
init_grid(points);
//初始化完。输出看看。
//for(int i=0;i<num_of_points;i++)
// cout<<i<<" => "<<points[i]<<'\n';
cout<<"初始温度分布（最外层和最内层没有显示）：\n";
print_grid(points);
//开始迭代计算
double diff=out_temp;
int times=0;
for(;diff>out_temp*accuracy;times++){
//"上=>下" 里嵌 "左=>右"遍历计算
diff=0.0;
for(int now=0;now<num_of_points;now++){
double tmp=points[now].temp;
compu_point(points,now);
tmp=abs(tmp-points[now].temp);
diff = diff>tmp?diff:tmp;
}
}
cout<<"\n迭代计算"<<times<<"次得到的温度分布（最外层和最内层没有显示）：\n";
print_grid(points);
heat_rate(points);
delete [] points;
return 0;
}
第二类边界条件：
/*
* 对流边界
*/
#include<iostream>
#include<cmath>
using namespace std;
const double out_temp=30.0;//外边界温度
const double out_h=10.0;//外边界对流换热系数
const double in_temp=10.0;//内边界温度
const double in_h=4
.0;//内边界对流换热系数
const double accuracy=0.00000001;//精度
const double lambda=0.53;//导热系数
const int width=15;//上部点阵宽度
const int width_bottom=4;//下部点阵宽度
const int height=4;//上部点阵高度
const int height_bottom=7;//下部点阵高度
const double delta_x=3.0/(2*width+1);
const double out_Bi=out_h*delta_x/lambda;//外表面网格Bi数
const double in_Bi = in_h*delta_x/lambda;//内表面网格Bi数
/*
//另一组参数
const int width=150;
const int width_bottom=40;
const int height=40;
const int height_bottom=70;
*/
//总共点数
const int num_of_points=width*height+width_bottom*height_bottom;
//单个点
class point
{
public:
double temp;//温度
int up;//数组下标
int down;
int left;
int right;
point(){
temp=0.0;//初始化成0摄氏度
up=down=right=left=0;
}
};
ostream & operator<<(ostream & src_stream,const point & src){
src_stream<<"temp="<<src.temp;
src_stream<<" up="<<src.up<<" down="<<src.down;
src_stream<<" left="<<src.left<<" right="<<src.right;
return src_stream;
}
void print_grid(point * points){//输出网格各点的温度
cout<<endl;
for(int position=0;position>=0;position=points[position].down){
for(int tmp=position;tmp>=0;tmp=points[tmp].right){
cout.width(10);
cout<<points[tmp].temp;
}
cout<<endl<<endl<<endl<<endl;
}
}
void dump_all(point * points){//便于调试
for(int i=0;i<num_of_points;i++)
cout<<i<<" => "<<points[i]<<'\n';
}
double calc_direction(point * points,int direction,double &opposit_weight,double &this_weight){
//计算给定方向邻点的温度和反方向的权重。
switch (direction) {
case -1://外边界
this_weight*=out_Bi;
return out_temp;
break;
case -2://内边界
this_weight*=in_Bi;
return in_temp;
break;
case -3://绝热
this_weight=0.0;
opposit_weight*=2;
return 0.0;
break;
default:
return points[direction].temp;
}
}
void compu_point(point * points,int now){
//根据周围四个点算出指定点的温度
double left_temp;
double right_temp;
double up_temp;
double down_temp;
double up_wight=1;//上方权重
double down_wight=1;
double left_wight=1;
double right_wight=1;
left_temp=calc_direction(points,points[now].left,right_wight,left_wight);
right_temp=calc_direction(points,points[now].right,left_wight,right_wight);
up_temp=calc_direction(points,points[now].up,down_wight,up_wight);
down_temp=calc_direction(points,points[now].down,up_wight,down_wight);
double Bi_sum=up_wight+down_wight+right_wight+left_wight;
points[now].temp=(left_temp*left_wight+right_temp*right_wight+up_temp*up_wight+down_temp*down_wight)/Bi_sum;
}
void rec_walk(point * points){
// "左=>右"里嵌"上=>下"遍历计算
for(int position=0;position>=0;position=points[position].right)
for(int tmp=position;tmp>=0;tmp=points[tmp].down)
compu_point(points,tmp);
}
void heat_rate(point * points){
//计算导热量
double out_sum=0.0;
double in_sum=0.0;
//外边界
int pos=0;
double tmp=0.0;
for(;points[pos].right>=0;pos=points[pos].right){
tmp=out_temp - points[pos].temp ;
if(points[pos].right<=0)
out_sum += 0.5*tmp;
else
out_sum+=tmp;
}
for(pos=points[0].down;points[pos].down>=0;pos=points[pos].down){
tmp=out_temp - points[pos].temp ;
if(points[pos].down<=0)
out_sum += 0.5*tmp;
else
out_sum+=tmp;
}
//内边界
for(pos=num_of_points-1;points[pos].right<0;pos=points[pos].up){
tmp = points[pos].temp - in_temp ;
if(points[pos].down<=0)
in_sum += 0.5*tmp;
else
in_sum +=tmp;
}
for(pos=points[pos].right;points[pos].right>=0;pos=points[pos].right){
tmp = points[pos].temp - in_temp;
if(points[pos].right<=0)
in_sum += 0.5*tmp;
else
in_sum += tmp;
}
out_sum*=out_h*delta_x;
in_sum *=in_h *delta_x;
cout<<"外边界导热量="<<out_sum;
cout<<"\n内边界导热量="<<in_sum;
cout<<"\n误差="<<(in_sum-out_sum)/in_sum*100<<"%"<<endl;
}
void init_grid(point * points){
//初始化网格
int position=0;
//初始化最上面的边
for(int i=0;i<width;i++){
points[i].up=-1;//-1表示外边界，-2内边界,-3表示绝热
points[i].left=i-1;
points[i].right=i+1;
points[i].down=i+width;
position=i;
}
points[0].left=-1;//对流
points[width-1].right=-3;//绝热
//上面区域
for(int i=1;i<height;i++){
for(int j=0;j<width;j++){
points[i*width+j].up =i*width+j-width;
points[i*width+j].left =i*width+j-1;
points[i*width+j].right=i*width+j+1;
points[i*width+j].down =i*width+j+width;
position=i*width+j;
}
points[i*width+0].left=-1;
points[i*width+width-1].right=-3;
}
position-=width;
position++;
for(int j=width_bottom;j<width;j++){
points[position+j].down =-2;//内边界
}
position+=width;
//下面区域
for(int i=0;i<height_bottom;i++){
for(int j=0;j<width_bottom;j++){
if(i)
points[position+i*width_bottom+j].up=position+i*width_bottom+j-width_bottom;
else
points[position+i*width_bottom+j].up=position+i*width_bottom+j-width;
points[position+i*width_bottom+j].left=position+i*width_bottom+j-1;
points[position+i*width_bottom+j].right=position+i*width_bottom+j+1;
points[position+i*width_bottom+j].down=position+i*width_bottom+j+width_bottom;
}
points[position+i*width_bottom+0].left=-1;
points[position+i*width_bottom+width_bottom-1].right=-2;
}
//下边界
position+=width_bottom*height_bottom;
position-=width_bottom;
for(int j=0;j<width_bottom;j++){
points[position+j].down=-3;
}
position+=width_bottom;
}
int main(){
point * points=new point[num_of_points];
init_grid(points);
//dump_all(points);
cout<<"初始温度分布（最外层和最内层没有显示）：\n";
print_grid(points);
//开始迭代计算
double diff=out_temp;
int times=0;
for(;diff>out_temp*accuracy;times++){
//"上=>下" 里嵌 "左=>右"遍历计算
diff=0.0;
for(int now=0;now<num_of_points;now++){
double tmp=points[now].temp;
compu_point(points,now);
tmp=abs(tmp-points[now].temp);
diff = diff>tmp?diff:tmp;
}
}
cout<<"\n迭代计算"<<times<<"次得到的温度分布（最外层和最内层没有显示）：\n";
print_grid(points);
heat_rate(points);
delete [] points;
return 0;
}
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK