背景：这是本学期凝固实验课的实验之一。这节课有两个数值模拟实验，第一个是二维常物性的，只有一种介质。而第二个实验是模拟凝固过程，稍微复杂一些。这篇文章是针对第一个实验写的，实验书上是按照显示差分进行的，这里改为隐式差分以便于计算。由于本人不是学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 <tuple>
inline int indexOf(int i,int j,int cols){
return (j-1)*cols+i-1;
}
inline std::pair<int,int> 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 <QtMath>
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 <QDebug>
#include <QObject>
#include <QVector>
#include <index.hpp>
#include <parameter.hpp>
#include <eigen/Eigen/Eigen>
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<double> 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<SparseMatrix<double>> m_lu;
//初始化模拟参数
bool initLUSolver();
//将计算结果从Eigen::VectorXd中复制到QVector容器中
void copyToQVector(const VectorXd vec, QVector<double> & qvec);
};
#endif // ENGINE_H
engine.cpp
#include "engine.h"
void Engine::init(){
jumpToStart();
bool ret = initLUSolver();
m_prepared = ret;
QString output = ret?"初始化成功":"初始化失败，请检查参数";
qDebug()<<output;
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<rows;j++){
for(int i=2;i<cols;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<double> & qvec){
qvec.clear();
if (vec.rows()!=qvec.size()){
qvec.resize(vec.rows());
}
VectorXd::Map(&(qvec[0]),qvec.size()) = vec;
}
实验结果
之后的结果均是再以下参属下计算得到的
一分钟时的温度场：
中心点的温度变化曲线：
程序运行截图