背景：这是本学期凝固实验课的实验之一。这节课有两个数值模拟实验，第一个是二维常物性的，只有一种介质。而第二个实验是模拟凝固过程，稍微复杂一些。这篇文章是针对第一个实验写的，实验书上是按照显示差分进行的，这里改为隐式差分以便于计算。由于本人不是学CS的，因此代码的质量可能不是很高。

简要说明：

二维非稳态传热、常物性、第一类边界条件、无内热源、

网格的划分

计算原理概述

直角坐标系内二维导热过程温度场控制微分方程：

若时间项采用向后差分，空间项采用中心差分，则方程可离散为：

令

​，则上式可简化为：

整理可得：

通过隐式差分，该问题转化成线性方程组的求解问题。只需要给出系数矩阵A和此时刻的温度​即可求出下一时刻的温度分布。 由于该线性方程组有一定有确切的解，因此其系数矩阵A一定为满秩的方阵，基于此可以选择LU分解法来快速求解该方程组。

为了方便计算机内的存储与计算，将某时刻各点的温度存储在一个列向量中T中，具体的存储方式为

​

不妨假设在划分网格的时候X轴方向一共有cols个格点，Y轴方向一共有rows个格点。 由于在C++语言中，数组的下标都是从0开始的，因此有以下的关系：

为了方便表示，定义一个的函数,且始终满足​,则系数矩阵A满足以下关系:

为了方便表示，定义一个的函数indexof(i,j),且始终满足indexof(i,j) = (j-1)*cols+i-1,则系数矩阵A满足以下关系:

当然，以上的表达式仅仅只考虑了差分方程，而没有考虑边界条件。

对于边界点(即i=1或i=cols或j=1或j=rows的点)，应满足:

此时系数矩阵A已经得到了，就可以根据初始温度场向后迭代计算了。由于LU分解以及其解方程的原理较复杂，因此这里不重新写了，而是直接使用一个名为Eigen的C++矩阵计算库来完成此过程。

程序

整个程序主要由以下部分组成index.hpp 实现上文中的indexof函数

parameter.hpp 保存、计算以及传递模拟时需要的参数

engine.h engine.cpp 调用Eigen库进行数值模拟

record.hpp 存储计算结果

plot.h plot.cpp 计算结果的绘制

widget.h widget.cpp widget.ui 主程序的绘制

main.cpp 程序的入口

Eigen线性代数运算库、Qwt绘图库以及Qt界面库

为了节省篇幅，这里只给出与计算有关的代码(前3项)

index.hpp

#ifndef INDEX_HPP#define INDEX_HPP#include

inline int indexOf(int i,int j,int cols){

return (j-1)*cols+i-1;

}

inline std::pair positionOf(int index,int cols){

//ret.first = i, ret.second = j //i和j从1开始 return std::make_pair(index%cols+1,index/cols+1);

}

#endif// INDEX_HPP

parameter.hpp

#ifndef PARAMETER_HPP#define PARAMETER_HPP#include

struct Parameter{

int xNums,yNums,stop;

//分别为：x方向上格点数、y方向上格点数、计算多少个时间步长的时候停下来 double t0,tw,a,dt,dx,dy,l1,l2,rx,ry;

//分别为：初始温度、边界温度(第一类边界条件)、热扩散系数、时间步长 //x方向网格步长、y方向网格步长、x方向总长度、y方向总长度、x方向网格傅里叶数、y方向网格傅里叶数 Parameter() = default;

Parameter(double t0,double tw,double a,

double dt,double l1,double dx,int stop){

init(t0,tw,a,dt,l1,dx,stop);

}

Parameter(double t0,double tw,double a,double dt,

double l1,double l2,double dx,double dy, int stop){

init(t0,tw,a,dt,l1,l2,dx,dy,stop);

}

inline void init(double t0,double tw,double a,

double dt,double l1,double dx, int stop){

init(t0,tw,a,dt,l1,l1,dx,dx,stop);

}

inline void init(double t0,double tw,double a,double dt,

double l1,double l2,double dx,double dy,int stop){

this->t0 = t0;

this->tw = tw;

this->a = a;

this->dt = dt;

this->l1 = l1;

this->l2 = l2;

this->dx = dx;

this->dy = dy;

this->stop = stop;

this->rx = (a*dt)/(dx*dx);

this->ry = (a*dt)/(dy*dy);

this->xNums = qRound(l1/dx)+1;

this->yNums = qRound(l2/dx)+1;

}

};

#endif// PARAMETER_HPP

engine.h

#ifndef ENGINE_H#define ENGINE_H

#include #include #include #include #include #include using namespace Eigen;

class Engine : public QObject

{

Q_OBJECT

public:

explicit Engine(QObject *parent = nullptr);

signals:

//准备好开始计算的信号 void prepared(bool ok);

//计算完成的信号 void finished(bool ok);

//完成一步计算的信号，通知主程序当前计算状态、步数、计算结果 void stepFinished(bool ok, int curStep, QVector data);

public slots:

//初始化 void init();

//计算下一时刻温度场 void next();

//反初始化 void unInit();

//将计算结果复位 void jumpToStart();

//进行计算前的准备 void prepare(bool ok);

//设置模拟参数 void setParameter(const Parameter parm);

private:

//记录是否初始化 bool m_prepared;

//现在计算到什么时候了 int m_currentStep;

//存储模拟参数 Parameter m_parm;

//存储温度场 VectorXd m_tFirst,m_tSecond;

//LU求解器 SparseLU> m_lu;

//初始化模拟参数 bool initLUSolver();

//将计算结果从Eigen::VectorXd中复制到QVector容器中 void copyToQVector(const VectorXd vec, QVector & qvec);

};

#endif// ENGINE_H

engine.cpp

#include "engine.h"

void Engine::init(){

jumpToStart();

bool ret = initLUSolver();

m_prepared = ret;

QString output = ret?"初始化成功":"初始化失败，请检查参数";

qDebug()<

emit prepared(ret);

}

void Engine::jumpToStart(){

int rows = m_parm.yNums;int cols = m_parm.xNums;

if(rows*cols<=0){

return ;

}

m_currentStep = 1;

m_tFirst.resize(rows*cols);

m_tSecond.resize(rows*cols);

m_tFirst.fill(m_parm.t0);

m_tSecond.fill(m_parm.t0);

auto mat = MatrixXd::Map(&m_tFirst[0],rows,cols);

mat.row(0).fill(m_parm.tw);

mat.row(rows-1).fill(m_parm.tw);

mat.col(0).fill(m_parm.tw);

mat.col(cols-1).fill(m_parm.tw);

}

bool Engine::initLUSolver(){

int rows = m_parm.yNums;int cols = m_parm.xNums;

double onePlusTwoR=1+2*(m_parm.rx+m_parm.ry);

MatrixXd A(MatrixXd::Identity(rows*cols,rows*cols));

for(int j=2;j

for(int i=2;i

int index = indexOf(i,j,cols);

A(index,index) = onePlusTwoR;

A(index,indexOf(i-1,j,cols)) = -m_parm.rx;

A(index,indexOf(i+1,j,cols)) = -m_parm.rx;

A(index,indexOf(i,j-1,cols)) = -m_parm.ry;

A(index,indexOf(i,j+1,cols)) = -m_parm.ry;

}

}

m_lu.compute(A.sparseView());

return m_lu.info()==Success;

}

void Engine::next(){

if(!m_prepared){

qDebug()<

emit finished(false);

return;

}

if(m_currentStep>m_parm.stop){

qDebug()<

emit finished(true);

return;

}

m_tSecond = m_lu.solve(m_tFirst);

m_tFirst = m_tSecond;

copyToQVector(m_tSecond,m_cache);

emit stepFinished(m_lu.info()==Success,m_currentStep,m_cache);

m_currentStep ++;

}

Engine::Engine(QObject *parent) : QObject(parent)

{

}

void Engine::prepare(bool ok){

if(ok){

init();

}else{

unInit();

}

}

void Engine::unInit(){

m_prepared = false;

emit prepared(false);

}

void Engine::setParameter(Parameter parm){

m_parm = parm;

}

void Engine::copyToQVector(const VectorXd vec, QVector & qvec){

qvec.clear();

if (vec.rows()!=qvec.size()){

qvec.resize(vec.rows());

}

VectorXd::Map(&(qvec[0]),qvec.size()) = vec;

}

实验结果

之后的结果均是再以下参属下计算得到的

一分钟时的温度场：

中心点的温度变化曲线：

程序运行截图