5.2.4 二维无内热源稳态导热程序

周星星9527

发布于 2019-08-19 11:23:41

6040

发布于 2019-08-19 11:23:41

文章被收录于专栏：javascript趣味编程

理论可以查看前文视频，或者编程基础。2D显式温度场内部节点迭代格式如下：

算例如下图矩形，左侧和下侧边界上温度分别是某一函数分布；右侧和上侧温度为0.

假设网格取5×3，即共6×4=24个节点。温度场求解程序如下：

Array.create = function(dimension, initialValue){
    var arr = [];
    for(var i=0;i<dimension;i++){
        arr[i]=initialValue;
    }
    return arr;
};

var Solution=function(tmpr){
  this.tmpr=tmpr;

  this.SetUpGeometryAndMesh=SetUpGeometryAndMesh;
  this.SetUpBoundaryCondition=SetUpBoundaryCondition;
  this.Solve=Solve;
};

function SetUpGeometryAndMesh(dx,dy){
  for(var i=0;i<5+1;i++){
    tmpr[i]=Array.create(3+1,0);
  }
}

function SetUpBoundaryCondition(dx,dy){
  for(var j=0;j<3+1;j++){
    tmpr[0][j]=4*j*(3.0-j);//左侧
    tmpr[5][j]=0;//右侧
  }

  for(var i=0;i<5+1;i++){//
    tmpr[i][0]=10*Math.sin(Math.PI*i/5.0)
    tmpr[i][3]=0;//
  }
}

function Solve(){
  var err,iter=0,maxIter,tmp,N=5*3/*节点总数*/;
  do{
    err=0/*误差*/,maxIter=100/*最大迭代次数*/;
    for(var i=1;i<5;i++){
      for(var j=1;j<3;j++){
        tmp=0.25*(tmpr[i+1][j]+tmpr[i-1][j]+tmpr[i][j+1]+tmpr[i][j-1]);
        err+=(tmp-tmpr[i][j])*(tmp-tmpr[i][j]);
        tmpr[i][j]=tmp;
      }
    }
    iter++;//迭代次数
    err=Math.sqrt(err/N);
    errors.push(err);iters.push(iter)
  }while(err>1E-4)

  console.info('最终误差：',err);
}

var tmpr=[],iters=[],errors=[];

function getHeat2DProfile(){
  var solution=new Solution(tmpr);
  var dx=0.1; var dy=0.1;//本特例网格dx与dy必须一致
  solution.SetUpGeometryAndMesh(dx,dy);//网格
  solution.SetUpBoundaryCondition(dx,dy);//设置边界条件
  solution.Solve();//求解

  var results=[];
  for(var j=0;j<3+1;j++){
    for(var i=0;i<5+1;i++){
      var arr=[j,i,tmpr[i][j].toFixed(1)];
      results.push(arr);
    }
  }
  return results;
}

module.exports = {
  getHeat2DProfile: getHeat2DProfile
}

上述程序是小程序新工科课程设计HTML5小程序里自带程序24个节点的温度场，使用百度echarts.heatmap可视化，虽然丑了点，但信息很全面，运行效果如下：

加密网格到51×31个，有兴趣的读者可以自行将温度场可视化！此处给计算的关键代码，仅供参考：

Array.create = function(dimension, initialValue){
    var arr = [];
    for(var i=0;i<dimension;i++){
        arr[i]=initialValue;
    }
    return arr;
};

var Solution=function(tmpr){
  this.tmpr=tmpr;

  this.SetUpGeometryAndMesh=SetUpGeometryAndMesh;
  this.SetUpBoundaryCondition=SetUpBoundaryCondition;
  this.Solve=Solve;
  this.HeatfluxEval=HeatfluxEval;
  this.GetContourData=GetContourData;
  this.ShowResults=ShowResults;
};

function SetUpGeometryAndMesh(dx,dy){
  for(var i=0;i<50+1;i++){
    tmpr[i]=Array.create(30+1,0);
    Qx[i]=Array.create(30+1,0);
    Qy[i]=Array.create(30+1,0);
  }
}

function SetUpBoundaryCondition(dx,dy){
  for(var j=0;j<30+1;j++){
    tmpr[0][j]=0.05*j*(30.0-j);//左侧
    tmpr[50][j]=0;//右侧
  }

  for(var i=0;i<50+1;i++){//
    tmpr[i][0]=10*Math.sin(Math.PI*i/50.0)
    tmpr[i][30]=0;//
  }

  tabulateResults(tmpr);
}

function Solve(){
  var err,iter=0,tmp,N=51*31/*节点总数*/;
  do{
    err=0/*误差*/,maxIter=1000/*最大迭代次数*/;
    for(var i=1;i<50;i++){
      for(var j=1;j<30;j++){
        tmp=0.25*(tmpr[i+1][j]+tmpr[i-1][j]+tmpr[i][j+1]+tmpr[i][j-1]);
        err+=(tmp-tmpr[i][j])*(tmp-tmpr[i][j]);
        tmpr[i][j]=tmp;
      }
    }
    iter++;//迭代次数
    err=Math.sqrt(err/N);
    errors.push(err);iters.push(iter)
  }while(err>1E-4)
  tabulateResults(tmpr);
  console.info('最终误差：',err);
}

function ShowResults(){}

var tmpr=[],Qx=[],Qy=[],iters=[],errors=[];

function onSolve(){
  var solution=new Solution(tmpr);
  var dx=0.1;  var dy=0.1;//本特例网格dx与dy必须一致
  solution.SetUpGeometryAndMesh(dx,dy);//网格
  solution.SetUpBoundaryCondition(dx,dy);//设置边界条件
  solution.Solve();//求解
  solution.HeatfluxEval(dx,dy);//由流函数计算流体速度
  solution.ShowResults();//绘制流函数及速度云图
}

最终结果应当如下：