【旧代码】传热过程数值模拟

传热过程数值模拟（《传热学》实验指导书第四部分第一题，第一，第二类边界条件）

2010年十月写的旧代码。

第一类边界条件是给定边界温度。

第二类是对流边界。

区域都是如下形状的:

1 2 3 4 5 6 7 8 9 10 11 12 13 14	-------------------------------- \| \| \| \| \| \| \| ---------------------- \| \| \| \| \| \| \| \| \| \| \| \| \| \| \| \| ---------

用C 纯属蛋疼。

第一类边界条件：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230

/* * 等温边界 */ #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; }

double int position sum width

0 人点赞