1一维热传导问题与C++描述

1.1理论推导

扩散问题的通用方程可以写为：

∇∇Γ+S=0(1)

一维热传导问题可以讲上式子简化为来：

ddx(Γdϕdx)+S=0(2) 将其在网格上离散可以得到下图的控制体单元表示：

将1.2式在网格上积分并利用Gauss定理就可以得到：

∫ΔVddx(Γdϕdx)dV+∫ΔVSdV=(ΓAdϕdx)e−(ΓAdϕdx)w+S¯ΔV=0(3) Γw 与 Γe 可以用线性插值的方法得到。1.3再在w,e面上的值即可以写为 (ΓAdϕdx)e=ΓeAeϕE−ϕPδEP(4) (ΓAdϕdx)w=ΓwAwϕP−ϕWδPW(5) 源项可以用 S¯ΔV=Su+Spϕp 代替，则式3可以写为： ΓeAeϕE−ϕPδEP−ΓwAwϕP−ϕWδPW+Su+Spϕp=0(6) 整理就得到： aPϕP=aWϕW+aEϕE+Su(7)

其中:

aW aE aP ΓWδxWPAW ΓEδxPEAE aW+aE+Su

1.2一维传热的小case

其中 δx 为 0.1m ， k=1000w/m⋅s ，截面积 A=0.01m2 。 根据公式7我们可以直接得到2、3、4节点处的方程组的系数。而对于边界1点处，可以类比公式6:

ΓAT2−T1δ21−ΓATP−TAδ1A=0(8)

同理对于5点我们也有:

ΓATB−T5δ5B−ΓAT5−T4δ54=0(8)

对于这个case我们所需的全部都已经得到，可以用C++编程语言进行具体的求解。

1.3 C++描述

1.3.1一个有用的数据结构

观察公式7、8、9我们发现实际上对于结构化网格，实际上网格任意一点所需求解的值都仅与相邻网格或者边界上的值相关，我们也就可以用

aPϕP=aWϕW+aEϕE+aNϕN+aSϕS+aTϕT+aBϕB+Su(7)

来表示任意一个网格节点处的离散化的控制方程。在这里就可以定义一个非常哟用的类或者结构体来储存这些系数。

#ifndef _LinearCofs_ #define _LinearCofs_ using namespace std; class LinearCofs { public: double ap_, ae_, aw_, an_, as_, at_, ab_; double b_; LinearCofs() :ap_(0), ae_(0), aw_(0), an_(0), as_(0), at_(0), ab_(0), b_(0) {} { } LinearCofs(double ap, double ae, double aw, double an, double as, double at, double ab, double b) :ap_(ap), ae_(ae), aw_(aw), an_(an), as_(as), at_(at), ab_(ab), b_(b) {} }; #endif

即使在后面的三维问题的求解中我们也可以利用上述类。

1.3.2 求解过程

#include <iostream> #include <vector> #include "LinearCofs.h" #include "Matrix.h" const double N = 5; const double k_ = 0.5; const double A_ = 1; const double L_ = 0.02; const double TA_ = 100; const double TB_ = 500; double deltx = L_ / N; using namespace std; void InerField(vector<LinearCofs>& lcDiffusion1D) { for( int i = 1; i < N - 1; i++) { lcDiffusion1D[i].ae_ = -k_ / deltx * A_; lcDiffusion1D[i].aw_ = -k_ / deltx * A_; lcDiffusion1D[i].ap_ = -(lcDiffusion1D[i].ae_ + lcDiffusion1D[i].aw_); } } void BoundaryCondition(vector<LinearCofs>& lcDiffusion1D) { double deltx = L_ / N; lcDiffusion1D[0].ae_ = -k_ / deltx * A_; lcDiffusion1D[0].aw_ = 0.0; lcDiffusion1D[0].ap_ = -(lcDiffusion1D[0].ae_ + lcDiffusion1D[0].aw_)+2 * k_ / deltx * A_; lcDiffusion1D[N - 1].ae_ = 0.0; lcDiffusion1D[N - 1].aw_ = -k_ / deltx * A_; lcDiffusion1D[N - 1].ap_ = -(lcDiffusion1D[N - 1].ae_ + lcDiffusion1D[N - 1].aw_)+2 * k_ / deltx * A_; } void CofsToMatrix(Matrix& mDiffusion1D, vector<LinearCofs>& lcDiffusion1D) { for(int i = 0; i < N; i++) { for(int j = 0; j < N; j++) { mDiffusion1D(i,j) = 0.0; } } mDiffusion1D(0,0) = lcDiffusion1D[0].ap_; mDiffusion1D(0,1) = lcDiffusion1D[0].ae_; for (int i = 1; i < N-1 ; i++) { mDiffusion1D(i,i-1) = lcDiffusion1D[i].aw_; mDiffusion1D(i,i) = lcDiffusion1D[i].ap_; mDiffusion1D(i,i+1) = lcDiffusion1D[i].ae_; } mDiffusion1D(N-1,N-2) = lcDiffusion1D[N-1].aw_; mDiffusion1D(N-1,N-1) = lcDiffusion1D[N-1].ap_; } void CofsToB_(vector<LinearCofs>& lcDiffusion1D) { lcDiffusion1D[0].b_ = 2 * k_ / deltx*A_ * TA_; lcDiffusion1D[N-1].b_ = 2 * k_/ deltx*A_ * TB_; } int main() { Matrix mDiffusion1D(N,N); vector<LinearCofs> lcDiffusion1D; AllocateDim(lcDiffusion1D, N); InerField(lcDiffusion1D); BoundaryCondition( lcDiffusion1D ); CofsToMatrix(mDiffusion1D, lcDiffusion1D); CofsToB_( lcDiffusion1D ); cout << "Matrix A is: " << endl; cout << mDiffusion1D; vector<double> x(N),b(N); for(int i = 0; i < N; i++) { b[i] = lcDiffusion1D[i].b_; } mDiffusion1D.GaussElimination(x, b); cout << endl; cout << " The solution is " << endl; for(int i = 0; i < N; i++) { cout <<"x["<<i<<"]="<< x[i] <<" "; } cout << endl; system("pause"); return 0; }

这里第一个函数InerField()是对内部的3\4\5节点赋值，BoundaryCondition()函数则是对两端的1\5点赋值，CofsToMatrix()是将得到的方程组放到矩阵中求解，矩阵求解器可以再我前面的文章中得到,大家可以直接复制，然后include 到头文件中。我们的矩阵求解器有两款线性方程组求解的算法，当然后面我们对于二维问题，因为比较复杂，会有一个gauss-seidel的求解器，会在后面给出那种直接求解不需要赋值到矩阵中的方法。 利用上述算法我们可以得到每个节点处的值了。