nozzle.cpp

拟一维喷管流动数值解代码，运行环境vs2019，自己写的，代码段有注释。cellField类的定义在我的博客里有，直接复制粘贴到class.h就好。
https://blog.csdn.net/weixin_39124457/article/details/101038218

#include“class.h“
#include
#include

constexpr auto gamma = 1.4;//气体比热比
constexpr auto R = 287.0;//空气常数
constexpr auto cfl = 0.5;//库朗数
constexpr auto timeStep = 1000;//时间步
constexpr auto deltaX = 0.1;//x方向间距

using namespace std;


int main（） {
	cout << “开始，t=0s“ << endl;
	cellField rho V T Aa;//定义物理量（rho-密度，V-速度，T-温度，A-截面积）
	cellField drhodt dVdt dTdtrhopVpTp;//定义预估步偏导预估值
	cellField drhopdt dVpdt dTpdtdrhodtavdVdtavdTdtav;//定义校正步偏导，校正值
	cellField X deltaTime;//定义空间步、时间步
	double i jk;//定义下标
	ofstream Rfile;//定义文件类
	Rfile.open（“data.dat“）;//打开文件

	//*********************初始化**************************//
	for （i = 0; i <= cellNumber; i++）
	{
		X[i] = i * deltaX;
		A[i] = 1.0 + 2.2 * （X[i] - 1.5）* （X[i] - 1.5）;//公式7-73，p219
		rho[i] = 1.0 - 0.314 * X[i];//公式7-74a，p220
		T[i] = 1.0 - 0.2314 * X[i];//公式7-74b，p220
		V[i] = （0.1 + 1.09 * X[i]） * sqrt（T[i]）;//公式7-74c，p220
		a[i] = sqrt（T[i]）;//公式7-75，p222
	}
	//输出初始场
	Rfile << “开始，t=0s“ << endl;
	Rfile.precision（3）;
	Rfile.setf（ios::fixed）;
	Rfile << “No\tx\tA\trho\tV\tT\ta“ << endl;
	for （i = 0; i <= cellNumber; i++）
	{
		Rfile << i << “\t“ << X[i]<< “\t“ << A[i]<< “\t“ << rho[i]
			<< “\t“ << V[i]<< “\t“ << T[i]<< “\t“ << a[i]<< endl;
	}
	//******************************************************//

	for （k = 0; k <= timeStep; k++）//开始计算timeStep个时间步
	{	
		//*******************求最小时间步*****************//
		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;
		}
		//*************************************************//

		//********************预估步***********************//
		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-51p215
			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）