2010年十月写的旧代码。
第一类边界条件是给定边界温度。
第二类是对流边界。
区域都是如下形状的:
--------------------------------
| |
| ----------------------
| |
| |
| |
| |
--------
用C++纯属蛋疼。
第一类边界条件：
1 /* 2 * 等温边界 3 */ 4 #include<iostream> 5 #include<cmath> 6 using namespace std; 7 const double out_temp=30.0;//外边界温度 8 const double in_temp=0.0;//内边界温度 9 const double accuracy=0.00000001;//精度 10 const double lambda=0.53;//导热系数 11 12 const int width=15;//上部点阵宽度 13 const int width_bottom=4;//下部点阵宽度 14 const int height=4;//上部点阵高度 15 const int height_bottom=7;//下部点阵高度 16 /* 17 //另一组参数 18 const int width=16; 19 const int width_bottom=6; 20 const int height=6; 21 const int height_bottom=6; 22 */ 23 //总共点数 24 const int num_of_points=width*height+width_bottom*height_bottom; 25 //单个点 26 class point 27 { 28 public: 29 double temp;//温度 30 int up;//数组下标 31 int down; 32 int left; 33 int right; 34 35 point(){ 36 temp=0.0;//初始化成0摄氏度 37 up=down=right=left=0; 38 } 39 }; 40 41 ostream & operator<<(ostream & src_stream,const point & src){ 42 src_stream<<"temp="<<src.temp; 43 src_stream<<" up="<<src.up<<" down="<<src.down; 44 src_stream<<" left="<<src.left<<" right="<<src.right; 45 return src_stream; 46 } 47 48 void print_grid(point * points){//输出网格各点的温度 49 cout<<endl; 50 for(int position=0;position>=0;position=points[position].down){ 51 for(int tmp=position;tmp>=0;tmp=points[tmp].right){ 52 cout.width(10); 53 cout<<points[tmp].temp; 54 } 55 cout<<endl<<endl<<endl<<endl; 56 } 57 } 58 59 double calc_direction(point * points,int direction,double &opposit_weight){ 60 //计算给定方向邻点的温度和反方向的权重。 61 switch (direction) { 62 case -1: 63 return out_temp; 64 break; 65 case -2: 66 return in_temp; 67 break; 68 case -3: 69 opposit_weight*=2; 70 return 0.0; 71 break; 72 default: 73 return points[direction].temp; 74 } 75 } 76 77 void compu_point(point * points,int now){ 78 //根据周围四个点算出指定点的温度 79 double left_temp; 80 double right_temp; 81 double up_temp; 82 double down_temp; 83 double up_wight=0.25;//上方权重,默认0.25 84 double down_wight=0.25; 85 double left_wight=0.25; 86 double right_wight=0.25; 87 left_temp=calc_direction(points,points[now].left,right_wight); 88 right_temp=calc_direction(points,points[now].right,left_wight); 89 up_temp=calc_direction(points,points[now].up,down_wight); 90 down_temp=calc_direction(points,points[now].down,up_wight); 91 points[now].temp=left_temp*left_wight+right_temp*right_wight+up_temp*up_wight+down_temp*down_wight; 92 } 93 94 void rec_walk(point * points){ 95 // "左=>右"里嵌"上=>下"遍历计算 96 for(int position=0;position>=0;position=points[position].right) 97 for(int tmp=position;tmp>=0;tmp=points[tmp].down) 98 compu_point(points,tmp); 99 } 100 101 void heat_rate(point * points){ 102 //计算导热量 103 double out_sum=0.0; 104 double in_sum=0.0; 105 //外边界 106 int pos=0; 107 double tmp=0.0; 108 for(;points[pos].right>=0;pos=points[pos].right){ 109 tmp=out_temp - points[pos].temp ; 110 if(points[pos].right<=0) 111 out_sum += 0.5*tmp; 112 else 113 out_sum+=tmp; 114 } 115 for(pos=points[0].down;points[pos].down>=0;pos=points[pos].down){ 116 tmp=out_temp - points[pos].temp ; 117 if(points[pos].down<=0) 118 out_sum += 0.5*tmp; 119 else 120 out_sum+=tmp; 121 } 122 //内边界 123 for(pos=num_of_points-1;points[pos].right<0;pos=points[pos].up){ 124 tmp = points[pos].temp - in_temp ; 125 if(points[pos].down<=0) 126 in_sum += 0.5*tmp; 127 else 128 in_sum +=tmp; 129 } 130 for(pos=points[pos].right;points[pos].right>=0;pos=points[pos].right){ 131 tmp = points[pos].temp - in_temp; 132 if(points[pos].right<=0) 133 in_sum += 0.5*tmp; 134 else 135 in_sum += tmp; 136 } 137 out_sum*=lambda; 138 in_sum*=lambda; 139 cout<<"外边界导热量="<<out_sum; 140 cout<<"\n内边界导热量="<<in_sum; 141 cout<<"\n误差="<<(in_sum-out_sum)/in_sum*100<<"%"<<endl; 142 } 143 144 void init_grid(point * points){ 145 //初始化网格 146 int position=0; 147 //初始化最上面的边 148 for(int i=0;i<width;i++){ 149 points[i].up=-1;//-1表示外边界，-2内边界,-3表示绝热 150 points[i].left=i-1; 151 points[i].right=i+1; 152 points[i].down=i+width; 153 position=i; 154 } 155 points[0].left=-1;//等温 156 points[width-1].right=-3;//绝热 157 158 //上面区域 159 for(int i=1;i<height;i++){ 160 for(int j=0;j<width;j++){ 161 points[i*width+j].up =i*width+j-width; 162 points[i*width+j].left =i*width+j-1; 163 points[i*width+j].right=i*width+j+1; 164 points[i*width+j].down =i*width+j+width; 165 position=i*width+j; 166 } 167 points[i*width+0].left=-1; 168 points[i*width+width-1].right=-3; 169 } 170 171 position-=width; 172 position++; 173 for(int j=width_bottom;j<width;j++){ 174 points[position+j].down =-2;//内边界 175 } 176 position+=width; 177 178 //下面区域 179 for(int i=0;i<height_bottom;i++){ 180 for(int j=0;j<width_bottom;j++){ 181 if(i) 182 points[position+i*width_bottom+j].up=position+i*width_bottom+j-width_bottom; 183 else 184 points[position+i*width_bottom+j].up=position+i*width_bottom+j-width; 185 186 points[position+i*width_bottom+j].left=position+i*width_bottom+j-1; 187 points[position+i*width_bottom+j].right=position+i*width_bottom+j+1; 188 points[position+i*width_bottom+j].down=position+i*width_bottom+j+width_bottom; 189 } 190 points[position+i*width_bottom+0].left=-1; 191 points[position+i*width_bottom+width_bottom-1].right=-2; 192 } 193 194 //下边界 195 position+=width_bottom*height_bottom; 196 position-=width_bottom; 197 for(int j=0;j<width_bottom;j++){ 198 points[position+j].down=-3; 199 } 200 position+=width_bottom; 201 } 202 203 int main(){ 204 point * points=new point[num_of_points]; 205 init_grid(points); 206 //初始化完。输出看看。 207 //for(int i=0;i<num_of_points;i++) 208 // cout<<i<<" => "<<points[i]<<'\n'; 209 cout<<"初始温度分布（最外层和最内层没有显示）：\n"; 210 print_grid(points); 211 //开始迭代计算 212 double diff=out_temp; 213 int times=0; 214 for(;diff>out_temp*accuracy;times++){ 215 //"上=>下" 里嵌 "左=>右"遍历计算 216 diff=0.0; 217 for(int now=0;now<num_of_points;now++){ 218 double tmp=points[now].temp; 219 compu_point(points,now); 220 tmp=abs(tmp-points[now].temp); 221 diff = diff>tmp?diff:tmp; 222 } 223 } 224 cout<<"\n迭代计算"<<times<<"次得到的温度分布（最外层和最内层没有显示）：\n"; 225 print_grid(points); 226 heat_rate(points); 227 delete [] points; 228 return 0; 229 }
第二类边界条件：
1 /* 2 * 对流边界 3 */ 4 #include<iostream> 5 #include<cmath> 6 using namespace std; 7 const double out_temp=30.0;//外边界温度 8 const double out_h=10.0;//外边界对流换热系数 9 const double in_temp=10.0;//内边界温度 10 const double in_h=4.0;//内边界对流换热系数 11 const double accuracy=0.00000001;//精度 12 const double lambda=0.53;//导热系数 13 14 const int width=15;//上部点阵宽度 15 const int width_bottom=4;//下部点阵宽度 16 const int height=4;//上部点阵高度 17 const int height_bottom=7;//下部点阵高度 18 19 const double delta_x=3.0/(2*width+1); 20 const double out_Bi=out_h*delta_x/lambda;//外表面网格Bi数 21 const double in_Bi = in_h*delta_x/lambda;//内表面网格Bi数 22 /* 23 //另一组参数 24 const int width=150; 25 const int width_bottom=40; 26 const int height=40; 27 const int height_bottom=70; 28 */ 29 //总共点数 30 const int num_of_points=width*height+width_bottom*height_bottom; 31 //单个点 32 class point 33 { 34 public: 35 double temp;//温度 36 int up;//数组下标 37 int down; 38 int left; 39 int right; 40 41 point(){ 42 temp=0.0;//初始化成0摄氏度 43 up=down=right=left=0; 44 } 45 }; 46 47 ostream & operator<<(ostream & src_stream,const point & src){ 48 src_stream<<"temp="<<src.temp; 49 src_stream<<" up="<<src.up<<" down="<<src.down; 50 src_stream<<" left="<<src.left<<" right="<<src.right; 51 return src_stream; 52 } 53 54 void print_grid(point * points){//输出网格各点的温度 55 cout<<endl; 56 for(int position=0;position>=0;position=points[position].down){ 57 for(int tmp=position;tmp>=0;tmp=points[tmp].right){ 58 cout.width(10); 59 cout<<points[tmp].temp; 60 } 61 cout<<endl<<endl<<endl<<endl; 62 } 63 } 64 65 void dump_all(point * points){//便于调试 66 for(int i=0;i<num_of_points;i++) 67 cout<<i<<" => "<<points[i]<<'\n'; 68 } 69 70 double calc_direction(point * points,int direction,double &opposit_weight,double &this_weight){ 71 //计算给定方向邻点的温度和反方向的权重。 72 switch (direction) { 73 case -1://外边界 74 this_weight*=out_Bi; 75 return out_temp; 76 break; 77 case -2://内边界 78 this_weight*=in_Bi; 79 return in_temp; 80 break; 81 case -3://绝热 82 this_weight=0.0; 83 opposit_weight*=2; 84 return 0.0; 85 break; 86 default: 87 return points[direction].temp; 88 } 89 } 90 91 void compu_point(point * points,int now){ 92 //根据周围四个点算出指定点的温度 93 double left_temp; 94 double right_temp; 95 double up_temp; 96 double down_temp; 97 double up_wight=1;//上方权重 98 double down_wight=1; 99 double left_wight=1; 100 double right_wight=1; 101 left_temp=calc_direction(points,points[now].left,right_wight,left_wight); 102 right_temp=calc_direction(points,points[now].right,left_wight,right_wight); 103 up_temp=calc_direction(points,points[now].up,down_wight,up_wight); 104 down_temp=calc_direction(points,points[now].down,up_wight,down_wight); 105 double Bi_sum=up_wight+down_wight+right_wight+left_wight; 106 points[now].temp=(left_temp*left_wight+right_temp*right_wight+up_temp*up_wight+down_temp*down_wight)/Bi_sum; 107 } 108 109 void rec_walk(point * points){ 110 // "左=>右"里嵌"上=>下"遍历计算 111 for(int position=0;position>=0;position=points[position].right) 112 for(int tmp=position;tmp>=0;tmp=points[tmp].down) 113 compu_point(points,tmp); 114 } 115 116 void heat_rate(point * points){ 117 //计算导热量 118 double out_sum=0.0; 119 double in_sum=0.0; 120 //外边界 121 int pos=0; 122 double tmp=0.0; 123 for(;points[pos].right>=0;pos=points[pos].right){ 124 tmp=out_temp - points[pos].temp ; 125 if(points[pos].right<=0) 126 out_sum += 0.5*tmp; 127 else 128 out_sum+=tmp; 129 } 130 for(pos=points[0].down;points[pos].down>=0;pos=points[pos].down){ 131 tmp=out_temp - points[pos].temp ; 132 if(points[pos].down<=0) 133 out_sum += 0.5*tmp; 134 else 135 out_sum+=tmp; 136 } 137 //内边界 138 for(pos=num_of_points-1;points[pos].right<0;pos=points[pos].up){ 139 tmp = points[pos].temp - in_temp ; 140 if(points[pos].down<=0) 141 in_sum += 0.5*tmp; 142 else 143 in_sum +=tmp; 144 } 145 for(pos=points[pos].right;points[pos].right>=0;pos=points[pos].right){ 146 tmp = points[pos].temp - in_temp; 147 if(points[pos].right<=0) 148 in_sum += 0.5*tmp; 149 else 150 in_sum += tmp; 151 } 152 out_sum*=out_h*delta_x; 153 in_sum *=in_h *delta_x; 154 cout<<"外边界导热量="<<out_sum; 155 cout<<"\n内边界导热量="<<in_sum; 156 cout<<"\n误差="<<(in_sum-out_sum)/in_sum*100<<"%"<<endl; 157 } 158 159 void init_grid(point * points){ 160 //初始化网格 161 int position=0; 162 //初始化最上面的边 163 for(int i=0;i<width;i++){ 164 points[i].up=-1;//-1表示外边界，-2内边界,-3表示绝热 165 points[i].left=i-1; 166 points[i].right=i+1; 167 points[i].down=i+width; 168 position=i; 169 } 170 points[0].left=-1;//对流 171 points[width-1].right=-3;//绝热 172 173 //上面区域 174 for(int i=1;i<height;i++){ 175 for(int j=0;j<width;j++){ 176 points[i*width+j].up =i*width+j-width; 177 points[i*width+j].left =i*width+j-1; 178 points[i*width+j].right=i*width+j+1; 179 points[i*width+j].down =i*width+j+width; 180 position=i*width+j; 181 } 182 points[i*width+0].left=-1; 183 points[i*width+width-1].right=-3; 184 } 185 186 position-=width; 187 position++; 188 for(int j=width_bottom;j<width;j++){ 189 points[position+j].down =-2;//内边界 190 } 191 position+=width; 192 193 //下面区域 194 for(int i=0;i<height_bottom;i++){ 195 for(int j=0;j<width_bottom;j++){ 196 if(i) 197 points[position+i*width_bottom+j].up=position+i*width_bottom+j-width_bottom; 198 else 199 points[position+i*width_bottom+j].up=position+i*width_bottom+j-width; 200 201 points[position+i*width_bottom+j].left=position+i*width_bottom+j-1; 202 points[position+i*width_bottom+j].right=position+i*width_bottom+j+1; 203 points[position+i*width_bottom+j].down=position+i*width_bottom+j+width_bottom; 204 } 205 points[position+i*width_bottom+0].left=-1; 206 points[position+i*width_bottom+width_bottom-1].right=-2; 207 } 208 209 //下边界 210 position+=width_bottom*height_bottom; 211 position-=width_bottom; 212 for(int j=0;j<width_bottom;j++){ 213 points[position+j].down=-3; 214 } 215 position+=width_bottom; 216 } 217 218 int main(){ 219 point * points=new point[num_of_points]; 220 init_grid(points); 221 //dump_all(points); 222 cout<<"初始温度分布（最外层和最内层没有显示）：\n"; 223 print_grid(points); 224 //开始迭代计算 225 double diff=out_temp; 226 int times=0; 227 for(;diff>out_temp*accuracy;times++){ 228 //"上=>下" 里嵌 "左=>右"遍历计算 229 diff=0.0; 230 for(int now=0;now<num_of_points;now++){ 231 double tmp=points[now].temp; 232 compu_point(points,now); 233 tmp=abs(tmp-points[now].temp); 234 diff = diff>tmp?diff:tmp; 235 } 236 } 237 cout<<"\n迭代计算"<<times<<"次得到的温度分布（最外层和最内层没有显示）：\n"; 238 print_grid(points); 239 heat_rate(points); 240 delete [] points; 241 return 0; 242 }
再见，我想我可以离开了。