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

 --------------------------------
|                               |
|                               |
|                               |
|          ----------------------
|         |     
|         |     
|         |     
|         |
|         |     
|         |
|         |     
|         |
 ---------

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