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_(), ae_(), aw_(), an_(), as_(), at_(), ab_(), b_()
{}
{
}
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 = ;
const double k_ = ;
const double A_ = ;
const double L_ = ;
const double TA_ = ;
const double TB_ = ;
double deltx = L_ / N;
using namespace std;
void InerField(vector<LinearCofs>& lcDiffusion1D)
{
for( int i = ; i < N - ; 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[].ae_ = -k_ / deltx * A_;
lcDiffusion1D[].aw_ = ;
lcDiffusion1D[].ap_ = -(lcDiffusion1D[].ae_ + lcDiffusion1D[].aw_)+ * k_ / deltx * A_;
lcDiffusion1D[N - ].ae_ = ;
lcDiffusion1D[N - ].aw_ = -k_ / deltx * A_;
lcDiffusion1D[N - ].ap_ = -(lcDiffusion1D[N - ].ae_ + lcDiffusion1D[N - ].aw_)+ * k_ / deltx * A_;
}
void CofsToMatrix(Matrix& mDiffusion1D, vector<LinearCofs>& lcDiffusion1D)
{
for(int i = ; i < N; i++)
{
for(int j = ; j < N; j++)
{
mDiffusion1D(i,j) = ;
}
}
mDiffusion1D(,) = lcDiffusion1D[].ap_;
mDiffusion1D(,) = lcDiffusion1D[].ae_;
for (int i = ; i < N- ; i++)
{
mDiffusion1D(i,i-) = lcDiffusion1D[i].aw_;
mDiffusion1D(i,i) = lcDiffusion1D[i].ap_;
mDiffusion1D(i,i+) = lcDiffusion1D[i].ae_;
}
mDiffusion1D(N-,N-) = lcDiffusion1D[N-].aw_;
mDiffusion1D(N-,N-) = lcDiffusion1D[N-].ap_;
}
void CofsToB_(vector<LinearCofs>& lcDiffusion1D)
{
lcDiffusion1D[].b_ = * k_ / deltx*A_ * TA_;
lcDiffusion1D[N-].b_ = * 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 = ; i < N; i++)
{
b[i] = lcDiffusion1D[i].b_;
}
mDiffusion1D.GaussElimination(x, b);
cout << endl;
cout << " The solution is " << endl;
for(int i = ; i < N; i++)
{
cout <<"x["<<i<<"]="<< x[i] <<" ";
}
cout << endl;
system("pause");
return ;
}
這裡第一個函數InerField()是對内部的3\4\5節點指派，BoundaryCondition()函數則是對兩端的1\5點指派，CofsToMatrix()是将得到的方程組放到矩陣中求解，矩陣求解器可以再我前面的文章中得到,大家可以直接複制，然後include 到頭檔案中。我們的矩陣求解器有兩款線性方程組求解的算法，當然後面我們對于二維問題，因為比較複雜，會有一個gauss-seidel的求解器，會在後面給出那種直接求解不需要指派到矩陣中的方法。
利用上述算法我們可以得到每個節點處的值了。