一、?簡介：?一維激波管內流動是一個經典的流體問題，如圖1，激波管左側壓強為p1,右側壓強為p2，中間用隔膜隔開，利用CFD技術模擬撤掉隔膜后管內的流動。

二、基礎知識：

? ? 忽略粘性影響，運用歐拉方程：

????????? ? ??

????其中U和F均為矩陣形式。

氣體內能： ? ? ? ? e=CvT=RT/(γ-1);（γ=Cp/Cv為比熱比)

氣體 ?焓 ：??????? ??h=CpT=γRT/(γ-1);

氣體總能： ?????????E=e+u2/2;

氣體總焓： ?????????H=h+u2/2; ? ??H=E+p/ρ=E+RT；

????????采用CFD方法在計算機上模擬流體流動的本質是數值求解N-S方程，或忽略粘性的歐拉方程；而CFD的本質是流動控制方程的計算機程序。

三、CFD過程

① 準備工作：

圖 2 一維激波管網格劃分

????????用差商代替偏微分形式：

上述格式為FTCS格式（時間向前，空間差分的首字母縮寫），若直接求解此方程會導致數值發散，所以需要改變形式。

所以，我們采用Lax格式求解：

展開格式為：

分析表明，Lax格式是條件穩定的，結果收斂的條件為：

CFL=|a|△t/△x<=1

②程序設計：

1.初始化各變量P,T,ρ,E,H,u,γ,?x,R，計算音速a，CFL=0.98。P的前半段為P1=10bar，后半段為P2=1bar；T=300K；ρ由P1&P2以及R和T，分別計算出ρ1和ρ2；E1=e1+u2=e1=RT/γ-1，E2同理；H=E+p/ρ=E+RT；u=0；γ=1.4；?x=0.1；R=287。

2.計算滿足條件的?t=?CFL*?x / max(a+u);

3.求解Lax格式的歐拉方程。

4.當達到所需物理時間時，停止計算，否則再次執行2、3步驟。

三、C++源碼

//#include<iostream>

#include<fstream>
#include<math.h>
#include<string>
#include<algorithm>
using namespace std;
double max(double *B, int n)//求數組B[n]中的最大值
{
double Max = *B;
for (int i = 0; i<n; ++i)
{
if (Max<*(B + i))
Max = B[i];
}
return Max;
}
int main()
{
const int nnod = 101;//節點數
ofstream out("out.txt");//將結果輸出到該文件
double timeOut = 0.005;//要計算到該時間
double lendom = 10.0;//激波管長度
double dx = lendom / (nnod - 1);
double X[nnod];
X[0] = 0.0;
double Tini = 300.0;//溫度
double PiniL = 10.0;//壓力P1
double PiniR = 1.0;//壓力P2
double uini = 0.0;//初始速度
double gamma = 1.4;//比熱比
double R = 287.0;//理想氣體常數J/kg/K
double CFL = 0.98;
double DiniL = 100000.0*PiniL / (R*Tini);//左邊密度
double DiniR = 100000.0*PiniR / (R*Tini);//右邊密度
double Eini = R*Tini / (gamma - 1.0);//初始內能
double Hini = Eini + R*Tini;//初始總焓
double aini = sqrt(gamma*R*Tini);//初始聲速

double Den[nnod] = { 0 };
double Tem[nnod] = { 0 };
double Pre[nnod] = { 0 };
double Vel[nnod] = { 0 };
double M[nnod] = { 0 };
double Ene[nnod] = { 0 };
double Ent[nnod] = { 0 };
double a[nnod] = { 0 };
double a_v[nnod] = { 0 };
fill_n(Den, 101, 1.0);
fill_n(Tem, 101, 1.0);
fill_n(Pre, 101, 1.0);
fill_n(Vel, 101, 1.0);
fill_n(M, 101, 1.0);
fill_n(Ene, 101, 1.0);
fill_n(a, 101, 1.0);
double dt = CFL*dx / aini;
double U1[101] = { 0 };
double U2[101] = { 0 };
double U3[101] = { 0 };
fill_n(U1, 50, DiniL);
fill_n(U1+50, 51, DiniR);
fill_n(U2, 101, 0);
fill_n(U3, 50, DiniL*Eini);
fill_n(U3+50, 51, DiniR*Eini);


double Un[3][101] = { 0 };
double Un1[3][101] = { 0 };
memcpy(Un[0], U1, nnod*sizeof(double));
memcpy(Un[1], U2, nnod*sizeof(double));
memcpy(Un[2], U3, nnod*sizeof(double));
memcpy(Un1[0], Un[0], 101 * sizeof(double));
memcpy(Un1[1], Un[1], 101 * sizeof(double));
memcpy(Un1[2], Un[2], 101 * sizeof(double));
double F1[nnod] = { 0 };
double F2[nnod] = { 0 };
double F3[nnod] = { 0 };
double Fn[3][nnod];
fill_n(F2, 50, PiniL);
fill_n(F2+50, 51, PiniR);
memcpy(Fn[0], F1, nnod*sizeof(double));
memcpy(Fn[1], F2, nnod*sizeof(double));
memcpy(Fn[2], F3, nnod*sizeof(double));
double time = 0.0;
double DX(0);
for (int i(0); i<101; ++i)
{
X[i] = DX;
DX += dx;


}
while (time<timeOut)
{
time += dt;
for (int j(0); j<nnod; ++j)
{
if (j !=0 && j != (nnod-1))
{
Un1[0][j] = 0.5*(Un[0][j - 1] + Un[0][j + 1]);
Un1[1][j] = 0.5*(Un[1][j - 1] + Un[1][j + 1]);
Un1[2][j] = 0.5*(Un[2][j - 1] + Un[2][j + 1]);


Un1[0][j] = Un1[0][j]-dt/dx*(Fn[0][j + 1] - Fn[0][j - 1])*0.5;
Un1[1][j] = Un1[1][j]-dt/dx*(Fn[1][j + 1] - Fn[1][j - 1])*0.5;
Un1[2][j] = Un1[2][j]-dt/dx*(Fn[2][j + 1] - Fn[2][j - 1])*0.5;
}
}
for (int i(0); i<101; ++i)
{
Den[i] = Un1[0][i];
Vel[i] = Un1[1][i] / Den[i];
Ene[i] = Un1[2][i] / Den[i];
Tem[i] = (Ene[i] - 0.5*Vel[i] * Vel[i])*(gamma - 1.0) / R;
Pre[i] = R*Den[i] * Tem[i];
Ent[i] = Ene[i] + R*Tem[i];
a[i] = sqrt(gamma*R*Tem[i]);
M[i] = Vel[i] / a[i];
a_v[i] = a[i] + Vel[i];
}
dt = CFL*dx / (max(a_v, nnod));
//memcpy(Un1[0], Un[0], 3 * 101);
memcpy(Un[0], Un1[0], 101 * sizeof(double));
memcpy(Un[1], Un1[1], 101 * sizeof(double));
memcpy(Un[2], Un1[2], 101 * sizeof(double));
for (int j(0); j<101; ++j)
{
Fn[0][j] = Den[j] * Vel[j];
}
for (int j(0); j<101; ++j)
{
Fn[1][j] = Den[j] * Vel[j] * Vel[j]+Pre[j];
}
for (int j(0); j<101; ++j)
{
Fn[2][j] = Den[j] * Vel[j] * Ent[j];
}
}
out << "nnod"<< " "<<"dencity"<<" " << "pressure" << " " << "velocity" << " " << "Ma" << endl;
for (int i(0); i<101; ++i)
{
out << i << " " <<Den[i]<<" "<< Pre[i]<<" "<<Vel[i]<<" "<<M[i] << endl;

}

return 0;

}

總結

以上是生活随笔為你收集整理的一维激波管内流动CFD实现（附C++源码）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。