背景：这是本学期凝固实验课的实验之一。这节课有两个数值模拟实验，第一个是二维常物性的，只有一种介质。而第二个实验是模拟凝固过程，稍微复杂一些。这篇文章是针对第一个实验写的，实验书上是按照显示差分进行的，这里改为隐式差分以便于计算。由于本人不是学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;
}
}
pute(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;
}
实验结果
之后的结果均是再以下参属下计算得到的
一分钟时的温度场：
中心点的温度变化曲线：
程序运行截图
