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

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;
}

第二类边界条件:

  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
231
232
233
234
235
236
237
238
239
240
241
242
243
/*
 * 对流边界
 */
#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;
}