這篇博文將剩下的問題解決完。第一篇在：

https://mp.csdn.net/postedit/101038218

在初始化步驟之后，就到了計算下一時間步的步驟了。計算之前先講一講這里用到的計算方法：顯式麥考馬克方法。簡單地來講，顯式意味著可以根據(jù)當前時間步的流場量去算下一個時間步的流場量。相對的，隱式則要用到下一個時間步的流場量，換言之要求解方程組，計算及其復(fù)雜。麥考馬克方法是比較早期的十分簡單的方法，由預(yù)估-校正兩步組成，每一步分別只有一階精度，但合在一起就有了二階精度，具體可以看看書中介紹。

一、求最小時間步

求最小時間步的目的是讓各個格點的庫朗數(shù)不高于預(yù)設(shè)值（這里是0.5），顯式格式中，庫朗數(shù)若大于1則會出現(xiàn)不穩(wěn)定的現(xiàn)象，因此這里最好選擇最小時間步來統(tǒng)一時間的推進。

cellField Va = V + a; deltaTime.inverse(Va, cfl * deltaX);//公式7-67，p217 j = 0; for (i = 0; i <= cellNumber; i++) {if (deltaTime[j] > deltaTime[i])j = i;//求最小時間間隔，將下標儲存為j，即deltaTime[j]為最小時間步//cout << deltaTime[j] << endl; }

deltaTime.inverse是在類定義中設(shè)置的一個函數(shù)，用來算公式：

將算到的最小時間步對應(yīng)的下標存在j里面，下面的計算會用到。

二、預(yù)估步

for (i = 0; i < cellNumber; i++) {drhodt[i] = -V[i] * (rho[i + 1] - rho[i]) / deltaX- rho[i] * (V[i + 1] - V[i]) / deltaX- rho[i] * V[i] * (log(A[i + 1]) - log(A[i])) / deltaX;//7-51,p215dVdt[i] = -V[i] * (V[i + 1] - V[i]) / deltaX- 1 / gamma * ((T[i + 1] - T[i]) / deltaX + T[i] / rho[i] * (rho[i + 1] - rho[i]) / deltaX);//7-52dTdt[i] = -V[i] * (T[i + 1] - T[i]) / deltaX- (gamma - 1) * T[i] * ((V[i + 1] - V[i]) / deltaX + V[i] * (log(A[i + 1]) - log(A[i])) / deltaX);//7-53rhop[i] = rho[i] + drhodt[i] * deltaTime[j];//7-54,p215Vp[i] = V[i] + dVdt[i] * deltaTime[j];//7-55Tp[i] = T[i] + dTdt[i] * deltaTime[j];//7-56 }

對應(yīng)公式：

7-52等號右端第一項應(yīng)該帶負號，書中有誤，代碼里已更正。

三、校正步

for (i = 1; i < cellNumber; i++) {drhopdt[i] = -Vp[i] * (rhop[i] - rhop[i-1]) / deltaX- rhop[i] * (Vp[i] - Vp[i - 1]) / deltaX- rhop[i] * Vp[i] * (log(A[i]) - log(A[i - 1])) / deltaX;//7-57,p216dVpdt[i] = -Vp[i] * (Vp[i] - Vp[i - 1]) / deltaX- 1 / gamma * ((Tp[i] - Tp[i - 1]) / deltaX + Tp[i] / rhop[i] * (rhop[i] - rhop[i - 1]) / deltaX);//7-58dTpdt[i] = -Vp[i] * (Tp[i] - Tp[i - 1]) / deltaX- (gamma - 1) * Tp[i] * ((Vp[i] - Vp[i - 1]) / deltaX + Vp[i] * (log(A[i]) - log(A[i - 1])) / deltaX);//7-59drhodtav[i] = (drhodt[i] +drhopdt[i])/2;//7-60,p215dVdtav[i] = (dVdt[i] +dVpdt[i])/2;//7-61dTdtav[i] = (dTdt[i] +dTpdt[i])/2;//7-62rho[i] += drhodtav[i] * deltaTime[j];//7-63V[i] += dVdtav[i] * deltaTime[j];//7-64T[i] += dTdtav[i] * deltaTime[j];//7-65 }

對應(yīng)公式：

到此，計算的大部頭已完成，接下來算頭和尾。

四、計算頭尾

除了密度、溫度在x=0處設(shè)定為1之外，其他頭尾的值都用插值法估算。

rho[0] = 1; V[0] = 2 * V[1] - V[2]; T[0] = 1;//7-70 ,p219 x=0處賦值 rho[cellNumber] = 2 * rho[cellNumber - 1] - rho[cellNumber - 2]; V[cellNumber] = 2 * V[cellNumber - 1] - V[cellNumber - 2]; T[cellNumber] = 2 * T[cellNumber - 1] - T[cellNumber - 2];//7-72, 結(jié)尾處賦值

對應(yīng)公式：

總結(jié)

以上是生活随笔為你收集整理的持续不定期更新：CFDC++之拟一维喷管流动的数值解（2）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。