持续不定期更新：CFD&C++之拟一维喷管流动的数值解（2）

这篇博文将剩下的问题解决完。第一篇在：

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,p215
	dVdt[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-52
	dTdt[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-53
	rhop[i] = rho[i] + drhodt[i] * deltaTime[j];//7-54,p215
	Vp[i] = V[i] + dVdt[i] * deltaTime[j];//7-55
	Tp[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,p216
	dVpdt[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-58
	dTpdt[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-59
	drhodtav[i] = (drhodt[i] +drhopdt[i])/2;//7-60,p215
	dVdtav[i] = (dVdt[i] +dVpdt[i])/2;//7-61
	dTdtav[i] = (dTdt[i] +dTpdt[i])/2;//7-62
	rho[i] += drhodtav[i] * deltaTime[j];//7-63
	V[i] += dVdtav[i] * deltaTime[j];//7-64
	T[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, 结尾处赋值

对应公式：