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

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

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

一、求最小時間步

求最小時間步的目的是讓各個格點的庫朗數不高于預設值（這里是0.5），顯式格式中，庫朗數若大于1則會出現不穩定的現象，因此這里最好選擇最小時間步來統一時間的推進。

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是在類定義中設置的一個函數，用來算公式：

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

二、預估步

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 }

對應公式：

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

三、校正步

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 }

對應公式：

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

四、計算頭尾

除了密度、溫度在x=0處設定為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, 結尾處賦值

對應公式：

總結

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

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