2021-04-24

拟一维喷管流动的数值解程序分享

亚声速——超声速等熵喷管流动的CFD解法（守恒型方程）

本文主要分享了安德森的《计算流体力学基础与应用》一书中：“亚声速——超声速等熵喷管流动的CFD解法 — 守恒型方程”的程序求解。水平有限，请见谅！
亚声速——超声速等熵喷管流动的CFD解法的非守恒型方程程序见之前的笔者的博客：
《计算流体力学基础与应用 》- 约翰 D. 安德森
《Computational Fluid Dynamics》- John D. Anderson
第 7 章 拟一维喷管流动的数值解
亚声速——超声速等熵喷管流动的CFD解法（非守恒型方程）

Matlab

Nozzle.m

%%
% Author: CXT
% Date: 2021/4/24 
% Theme: 拟一维亚声速-超声速等熵喷管 - 守恒型求解
%% 物理符号(无量纲量)；
% x —— 空间位置；
% A —— 喷管面积；
% rho —— 密度；
% V —— 速度；
% p —— 压力；
% N —— 计算节点数；
% m —— 质量流量；
%%
clc;
clear;
close all;
%% 喷管形状;
N = 31;% 计算节点数；
x = 0:(3-0)/(N-1):3; % 计算区间；
n = 1400; %时间步数；

A = 1 + 2.2*(x - 1.5).^2;
A = A./1;
%% 初始条件；
% 原始变量；
rho0 = 1.*(x<=0.5)...
     + (1.0 - 0.366*(x - 0.5)).*(x>0.5 & x<1.5) ...
     + (0.634 - 0.3879*(x - 1.5)).*(x>=1.5 & x<=3.0);
T0 = 1.*(x<=0.5)...
   + (1.0 - 0.167*(x - 0.5)).*(x>0.5 & x<1.5) ...
   + (0.833 - 0.3507*(x - 1.5)).*(x>=1.5 & x<=3.0);
V0 = 0.59./(rho0.*A);% 预设 U2 = 0.59; 

m0 = rho0.*A.*V0;% 恒为0.59；
%% 计算与分析；
rho_t = rho0;
T_t = T0;
V_t = V0;

tic;
for i = 1:n
    [rho_t, V_t, T_t, p_t, Ma, m] = OneStepNozzle(rho_t, V_t, T_t, N, x, A);
    
    if (i == 100)
        m1 = m;
    elseif (i == 200)
        m2 = m;
    elseif (i == 700)
        m3 = m;
    end
  
    plot(i,rho_t(16),'*');
    hold on;
end
xlabel("时间步数");
ylabel("无量纲密度");
hold off;
toc;

figure;
title("时间推进的过程中，质量流量在不同时刻的分布");
xlabel('x/L');
ylabel("无量纲质量流量");
plot(x, m0,'g');
hold on;
plot(x, m1,'m');
plot(x, m2,'b');
plot(x, m3,'c');
legend('0\Delta t','100\Delta t','200\Delta t','700\Delta t');
hold off;

OneStepNozzle.m

function [rho_t, V_t, T_t, p_t, Ma, m] = OneStepNozzle(rho, V, T, N, X, A)
%[rho_t, V_t, T_t, p_t, Ma] = OneStepNozzle(rho, V, T) 拟一维喷管流动的一个时间步长的计算
%   rho —— 初始密度；
%   V —— 初始速度；
%   T —— 初始温度；
%   N —— 计算节点数；
%   X —— 计算区间；
%   A —— 喷管形状参数；

%   rho_t —— 下一步的密度；
%   V_t —— 下一步的速度；
%   T_t —— 下一步的温度；
%   p_t —— 下一步的压力；
%   Ma —— 马赫数；
%   m —— 质量流量；

%% 默认参数；
Delta_x = (X(end)-X(1))/(N - 1); %空间步长；
gama = 1.4; %绝热指数；
C = 0.5; %柯朗数；
%% 解向量 U；
U1 = rho.*A;
U2 = rho.*A.*V; % 质量流量；
U3 = rho.*(T/(gama - 1) + gama/2*V.^2).*A;% e0 = T0;

%% 计算通量向量 F; 
F1 = U2;
F2 = U2.^2./U1 + (gama - 1)/gama*(U3 - gama/2*U2.^2./U1);
F3 = gama*(U2.*U3)./U1 - gama*(gama-1)/2*U2.^3./U1.^2;

%% 计算源项；
J2 = zeros(1, N);
for i = 1:N-1
   J2(i) = 1/gama*rho(i)*T(i)*(A(i+1) - A(i))/Delta_x;
end
% 注意：J2(31)还未计算，等于 0;

%% 预估步计算；
% 向前差分；
% 只计算内点；
partial_U1 = zeros(1, N);
partial_U2 = zeros(1, N);
partial_U3 = zeros(1, N);

for i = 1:N - 1
    partial_U1(i) = -(F1(i+1) - F1(i))/Delta_x;
    partial_U2(i) = -(F2(i+1) - F2(i))/Delta_x + J2(i);
    partial_U3(i) = -(F3(i+1) - F3(i))/Delta_x;
end

%% 计算时间步长；
Delta_t = min(C * Delta_x ./ (V + sqrt(T)));

%% 计算预估值
% U1、U2和U3的预估值；
bar_U1 = U1 + partial_U1 * Delta_t;
bar_U2 = U2 + partial_U2 * Delta_t;
bar_U3 = U3 + partial_U3 * Delta_t;

% rho、T 的预估值；
bar_rho = bar_U1./A;
bar_T = (gama - 1) * (bar_U3./bar_U1 - gama/2 * (bar_U2./bar_U1).^2);

% F1、F2 和 F3 的预估值；
bar_F1 = bar_U2;
bar_F2 = bar_U2.^2./bar_U1 + (gama - 1)/gama*(bar_U3 - gama/2*bar_U2.^2./bar_U1);
bar_F3 = gama*(bar_U2.*bar_U3)./bar_U1 - gama*(gama-1)/2*bar_U2.^3./bar_U1.^2;

%% 计算校正步；
% 向后差分；
% 只计算内点；
bar_partial_U1 = zeros(1, N);
bar_partial_U2 = zeros(1, N);
bar_partial_U3 = zeros(1, N);

for i = 2:N
   bar_partial_U1(i) = -(bar_F1(i) - bar_F1(i-1))/Delta_x;
   bar_partial_U2(i) = -(bar_F2(i) - bar_F2(i-1))/Delta_x...
                      + 1/gama*bar_rho(i)*bar_T(i)*(A(i) - A(i-1))/Delta_x;
   bar_partial_U3(i) = -(bar_F3(i) - bar_F3(i-1))/Delta_x;
end

%% 计算平均时间导数；
partial_U1_av = (partial_U1 + bar_partial_U1)/2;
partial_U2_av = (partial_U2 + bar_partial_U2)/2;
partial_U3_av = (partial_U3 + bar_partial_U3)/2;

%% 校正值；
% t + Delta_t 时的参数；
% U1、U2 和 U3 的校正值；
U1_t = U1 + partial_U1_av * Delta_t;
U2_t = U2 + partial_U2_av * Delta_t;
U3_t = U3 + partial_U3_av * Delta_t;

% rho、V 、 T 和 p 的校正值；
rho_t = U1_t./A;
V_t = U2_t./U1_t;
T_t = (gama - 1)*(U3_t./U1_t - gama/2*V_t.^2);
p_t = rho_t.*T_t;
%% 边界点处的流场变量；
% 入流边界；
rho_t(1) = 1;
T_t(1) = 1;
U1_t(1) = rho_t(1) * A(1);% 定值；
U2_t(1) = 2 * U2_t(2) - U2_t(3);
V_t(1) = U2_t(1)/U1_t(1);
U3_t(1) = U1_t(1)*(T_t(1)/(gama - 1) + gama/2*V_t(1)^2);

% 出流边界；
U1_t(N) = 2*U1_t(N-1) - U1_t(N-2);
U2_t(N) = 2*U2_t(N-1) - U2_t(N-2);
U3_t(N) = 2*U3_t(N-1) - U3_t(N-2);

rho_t(N) = U1_t(N)/A(N);
V_t(N) = U2_t(N)/U1_t(N);
T_t(N) = (gama - 1)*(U3_t(N)/U1_t(N) - gama/2*V_t(N)^2);
p_t(N) = rho_t(N)*T_t(N);
%% 马赫数与质量流量；
Ma = V_t./sqrt(T_t);
m = U2_t;
end

C++

OneStepNozzle.h

#pragma once
#include <iostream>
#include <vector>
#include <iomanip>//setprecision(3)
using namespace std;

class Nozzle
{
private:
	double X_first;
	double X_end;//计算区间[X_first, X_end]；
	double* A; //喷管面积；
	double* rho;//密度；
	double* T;//温度；
	double* V;//速度；
	int N; //计算节点数；
	int iterations;//迭代次数；

protected:
	double* p;//压力；
	double* Ma;//马赫数; 
	double Delta_x;// 空间步长；
	double Delta_t;//时间步长；
	double gama;// 绝热指数；
	double C;// 柯朗数；
	double* m;

	中间变量
	//解向量 U；
	double* U1;
	double* U2;
	double* U3;
	//通量向量 F;
	double* F1;
	double* F2;
	double* F3;
	//源项;
	double* J2;

	//中间变量；
	double* partial_U1;
	double* partial_U2;
	double* partial_U3;

	double* bar_U1;
	double* bar_U2;
	double* bar_U3;

	double* bar_rho;
	double* bar_T;

	double* bar_F1;
	double* bar_F2;
	double* bar_F3;

	double* bar_partial_U1;
	double* bar_partial_U2;
	double* bar_partial_U3;

	double* partial_U1_av;
	double* partial_U2_av;
	double* partial_U3_av;

	double* U1_t;
	double* U2_t;
	double* U3_t;

	

public:
	void _init();

	Nozzle(double x_f, double x_e, double a[], double rho0[], double T0[], 
		   double V0[], int n, int it);
	
	Nozzle(const Nozzle& N);

	void oneStepmacCormack();

	void iterativeCompute();//迭代计算；

	void print();//输出最终的计算结果；

	~Nozzle();
};

OneStepNozzle.cpp

#include "OneStepNozzle.h"

void Nozzle::_init()
{
	this->Delta_x = (X_end - X_first) / (N - 1);// 空间步长；
	this->Delta_t = -1;//时间步长；
	this->gama = 1.4;// 绝热指数；
	this->C = 0.5;// 柯朗数；

	this->U1 = new double[N];
	this->U2 = new double[N];
	this->U3 = new double[N];
	this->F1 = new double[N];
	this->F2 = new double[N];
	this->F3 = new double[N];
	this->J2 = new double[N];

	this->partial_U1 = new double[N];
	this->partial_U2 = new double[N];
	this->partial_U3 = new double[N];

	this->bar_U1 = new double[N];
	this->bar_U2 = new double[N];
	this->bar_U3 = new double[N];

	this->bar_rho = new double[N];
	this->bar_T = new double[N];

	this->bar_F1 = new double[N];
	this->bar_F2 = new double[N];
	this->bar_F3 = new double[N];

	this->bar_partial_U1 = new double[N];
	this->bar_partial_U2 = new double[N];
	this->bar_partial_U3 = new double[N];

	this->partial_U1_av = new double[N];
	this->partial_U2_av = new double[N];
	this->partial_U3_av = new double[N];

	this->U1_t = new double[N];
	this->U2_t = new double[N];
	this->U3_t = new double[N];

	this->p = new double[N];
	this->Ma = new double[N];
	this->m = new double[N];
}

Nozzle::Nozzle(double x_f,double x_e,double a[], double rho0[],double T0[],
	           double V0[], int n,int it):
	X_first(x_f), X_end(x_e),A(a),rho(rho0),V(V0),T(T0), N(n), iterations(it)
{
	//用户需要输入:计算区间[x_f,x_e];
	//           喷管形状向量 a[];
	//           初始密度 rho0[];
	//           初始温度T0[];
	//           初始速度V0[];
	//           节点数 n;
	//           迭代次数 it;
	//           出口与入口压力之比 f;

	_init();
}

Nozzle::Nozzle(const Nozzle& _N)
{
	//浅拷贝；
	this->N = _N.N;
	this->X_first = _N.X_first;
	this->X_end = _N.X_end;

	this->Delta_x = _N.Delta_x;
	this->Delta_t = _N.Delta_t;

	this->gama = _N.gama;
	this->C = _N.C;

	this->iterations = _N.iterations;

	//深拷贝；
	this->rho = new double(*_N.rho);
	this->V = new double(*_N.V);
	this->T = new double(*_N.T);

	this->U1 = new double(*_N.U1);
	this->U2 = new double(*_N.U2);
	this->U3 = new double(*_N.U3);
	this->F1 = new double(*_N.F1);
	this->F2 = new double(*_N.F2);
	this->F3 = new double(*_N.F3);
	this->J2 = new double(*_N.J2);

	this->partial_U1 = new double(*_N.partial_U1);
	this->partial_U2 = new double(*_N.partial_U2);
	this->partial_U3 = new double(*_N.partial_U3);

	this->bar_U1 = new double(*_N.bar_U1);
	this->bar_U2 = new double(*_N.bar_U2);
	this->bar_U3 = new double(*_N.bar_U3);

	this->bar_rho = new double(*_N.bar_rho);
	this->bar_T = new double(*_N.bar_T);

	this->bar_F1 = new double(*_N.bar_F1);
	this->bar_F2 = new double(*_N.bar_F2);
	this->bar_F3 = new double(*_N.bar_F3);

	this->bar_partial_U1 = new double(*_N.bar_partial_U1);
	this->bar_partial_U2 = new double(*_N.bar_partial_U2);
	this->bar_partial_U3 = new double(*_N.bar_partial_U3);

	this->partial_U1_av = new double(*_N.partial_U1_av);
	this->partial_U2_av = new double(*_N.partial_U2_av);
	this->partial_U3_av = new double(*_N.partial_U3_av);

	this->U1_t = new double(*_N.U1_t);
	this->U2_t = new double(*_N.U2_t);
	this->U3_t = new double(*_N.U3_t);

	this->p = new double(*_N.p);
	this->Ma = new double(*_N.Ma);
	this->m = new double(*_N.m);

	this->A = new double(*_N.A);
}

void Nozzle::oneStepmacCormack()//使用Mac Cormack法计算；
{
	for (int i = 0;i < N;i++)
	{
		// 解向量 U；
		U1[i] = rho[i] * A[i];
		U2[i] = rho[i] * A[i] * V[i];// 质量流量；
		U3[i] = rho[i] * (T[i] / (gama - 1) + gama / 2 * pow(V[i] , 2)) * A[i];// e0 = T0;

		// 计算通量向量 F;
		F1[i] = U2[i];
		F2[i] = pow(U2[i] , 2) / U1[i] + (gama - 1) / gama * (U3[i] - gama / 2 * pow(U2[i] , 2) / U1[i]);
		F3[i] = gama * (U2[i] * U3[i]) / U1[i] - gama * (gama - 1) / 2 *pow( U2[i] , 3) / pow(U1[i] , 2);
	}

	// 计算源项；
	for (int i = 0;i <= N - 2;i++)
	{
		J2[i] = 1 / gama * rho[i] * T[i] * (A[i + 1] - A[i]) / Delta_x;
	}
	J2[N - 1] = 0;// 注意：J2[N-1]设为 0：计算内点时不会用到此值;

	 预估步计算；
	// 向前差分；
	// 只计算内点;
	for (int i = 0; i <= N - 2; i++)
	{
		partial_U1[i] = -(F1[i + 1] - F1[i]) / Delta_x;
		partial_U2[i] = -(F2[i + 1] - F2[i]) / Delta_x + J2[i];
		partial_U3[i] = -(F3[i + 1] - F3[i]) / Delta_x;
	}
	partial_U1[N - 1] = 0;
	partial_U2[N - 1] = 0;
	partial_U3[N - 1] = 0;

	计算最小时间步长；

	Delta_t = C * Delta_x / (V[0] + sqrt(T[0]));//将首节点计算的时间步长赋值给 Delta_t;

	for (int i = 1;i < N;i++)//从第个二节点计算到最后一个节点；
	{
		if (C * Delta_x / (V[i] + sqrt(T[i])) <= this->Delta_t)
		{
			this->Delta_t = C * Delta_x / (V[i] + sqrt(T[i]));
		}
	}

	// 计算预估值
	for (int i = 0;i < N;i++)
	{
		// U1、U2和U3的预估值；
		bar_U1[i] = U1[i] + partial_U1[i] * Delta_t;
		bar_U2[i] = U2[i] + partial_U2[i] * Delta_t;
		bar_U3[i] = U3[i] + partial_U3[i] * Delta_t;

		// rho、T 的预估值；
		bar_rho[i] = bar_U1[i] / A[i];
		bar_T[i] = (gama - 1) * (bar_U3[i] / bar_U1[i] - gama / 2 * pow(bar_U2[i] / bar_U1[i], 2));

		// F1、F2 和 F3 的预估值；
		bar_F1[i] = bar_U2[i];
		bar_F2[i] = pow(bar_U2[i],2)/bar_U1[i] + (gama - 1)/gama * (bar_U3[i] - gama/2*pow(bar_U2[i],2)/bar_U1[i]);
		bar_F3[i] = gama * (bar_U2[i]*bar_U3[i])/bar_U1[i] - gama * (gama-1)/2*pow(bar_U2[i],3)/pow(bar_U1[i],2);
	}

	 计算校正步；
	// 向后差分；
	// 只计算内点；
	for (int i = 1;i <= N - 1;i++)
	{
		bar_partial_U1[i] = -(bar_F1[i] - bar_F1[i - 1]) / Delta_x;
		bar_partial_U2[i] = -(bar_F2[i] - bar_F2[i - 1]) / Delta_x
			              + 1 / gama * bar_rho[i] * bar_T[i] * (A[i] - A[i - 1]) / Delta_x;
		bar_partial_U3[i] = -(bar_F3[i] - bar_F3[i - 1]) / Delta_x;
	}
	bar_partial_U1[0] = 0;
	bar_partial_U2[0] = 0;
	bar_partial_U3[0] = 0;

	 计算平均时间导数；
	for (int i = 0;i < N;i++)
	{
		partial_U1_av[i] = (partial_U1[i] + bar_partial_U1[i]) / 2;
		partial_U2_av[i] = (partial_U2[i] + bar_partial_U2[i]) / 2;
		partial_U3_av[i] = (partial_U3[i] + bar_partial_U3[i]) / 2;
	}
	
	 校正值；
	// t + Delta_t 时的参数；

	for (int i = 0;i < N;i++)
	{
		// U1、U2 和 U3 的校正值；
		U1_t[i] = U1[i] + partial_U1_av[i] * Delta_t;
		U2_t[i] = U2[i] + partial_U2_av[i] * Delta_t;
		U3_t[i] = U3[i] + partial_U3_av[i] * Delta_t;

		// rho、V 、 T 和 p 的校正值；
		rho[i] = U1_t[i] / A[i];
		V[i] = U2_t[i] / U1_t[i];
		T[i] = (gama - 1) * (U3_t[i] / U1_t[i] - gama / 2 * pow(V[i],2));
		p[i] = rho[i] * T[i];
	}

	// 边界点处的流场变量；
	// 入流边界；
	rho[0] = 1;
	T[0] = 1;
	U1_t[0] = rho[0] * A[0];// 定值；
	U2_t[0] = 2 * U2_t[1] - U2_t[2];
	V[0] = U2_t[0] / U1_t[0];
	U3_t[0] = U1_t[0] * (T[0] / (gama - 1) + gama / 2 * pow(V[0],2));

	// 出流边界；
	U1_t[N - 1] = 2 * U1_t[N - 2] - U1_t[N - 3];
	U2_t[N - 1] = 2 * U2_t[N - 2] - U2_t[N - 3];
	U3_t[N - 1] = 2 * U3_t[N - 3] - U3_t[N - 3];

	rho[N - 1] = U1_t[N - 1] / A[N - 1];
	V[N - 1] = U2_t[N - 1] / U1_t[N - 1];
	T[N - 1] = (gama - 1) * (U3_t[N - 1] / U1_t[N - 1] - gama / 2 * pow(V[N - 1],2));
	p[N - 1] = rho[N - 1] * T[N - 1];
	
	// 马赫数与质量流量；
	for (int i = 0;i < N;i++)
	{
		Ma[i] = V[i] / sqrt(T[i]);
		m[i] = U2_t[i];
	}
}

void Nozzle::iterativeCompute()//迭代计算；
{
	for (int it = 0;it < this->iterations;it++)
	{
		oneStepmacCormack();
	}
}

void Nozzle::print()
{	
	iterativeCompute();

	cout << "节点" << "\t" << "密度" << "\t" << "速度" << "\t" 
		 << "温度" << "\t" << "压力" << "\t" << "马赫数" << "\t"
		 << "质量流量" << endl;

	for (int i = 0;i < N;i++)
	{
		cout << i + 1 << "\t" 
			<< setprecision(4) << rho[i] << "\t" 
			<< setprecision(4) << V[i] << "\t"
			<< setprecision(4) << T[i] << "\t" 
			<< setprecision(4) << p[i] << "\t" 
			<< setprecision(4) << Ma[i] << "\t"
			<< setprecision(4) << m[i] << endl;
	}
}

Nozzle::~Nozzle()
{
	delete[] U1;
	delete[] U2;
	delete[] U3;
	delete[] F1;
	delete[] F2;
	delete[] F3;
	delete[] J2;

	delete[] partial_U1;
	delete[] partial_U2;
	delete[] partial_U3;

	delete[] bar_U1;
	delete[] bar_U2;
	delete[] bar_U3;

	delete[] bar_rho;
	delete[] bar_T;

	delete[] bar_F1;
	delete[] bar_F2;
	delete[] bar_F3;

	delete[] bar_partial_U1;
	delete[] bar_partial_U2;
	delete[] bar_partial_U3;

	delete[] partial_U1_av;
	delete[] partial_U2_av;
	delete[] partial_U3_av;

	delete[] U1_t;
	delete[] U2_t;
	delete[] U3_t;

	delete[] p;
	delete[] Ma;
	delete[] m;
}

拟一维喷管流动的数值解.cpp

/*
 Author: CXT
 Date : 2021 / 4 / 24
 Theme : 拟一维亚声速-超声速等熵喷管 - 守恒型求解
 %物理符号(无量纲量)；
 x —— 空间位置；
 A —— 喷管面积；
 rho —— 密度；
 V —— 速度；
 p —— 压力；
 N —— 计算节点数；
 */

#include <iostream>
#include <cmath>//sqrt();
#include <ctime>
#include "OneStepNozzle.h"
using namespace std;

void test()
{
	double x_f = 0;//计算区间[x_f,x_e];
	double x_e = 3;
	int n = 31;//节点数；
	int it = 1400;//迭代次数；

	double* x = new double[n];
	double* a = new double[n];
	double* rho0 = new double[n];
	double* T0 = new double[n];
	double* V0 = new double[n];

	//准备节点；
	for (int i = 0; i<n; i++)
	{
		x[i] = x_f + (x_e - x_f) / (n - 1) * i;

		a[i] = 1 + 2.2 * pow(x[i] - 1.5, 2);//喷管形状；

		// 初始条件；
		rho0[i] = 1 * (x[i] <= 0.5)
			+ (1.0 - 0.366 * (x[i] - 0.5))*(x[i] > 0.5 && x[i] < 1.5) 
			+ (0.634 - 0.3879 * (x[i] - 1.5))*(x[i] >= 1.5 && x[i] <= 3.0);
		T0[i] = 1 * (x[i] <= 0.5)
			+ (1.0 - 0.167 * (x[i] - 0.5))*(x[i] > 0.5 && x[i] < 1.5)
			+ (0.833 - 0.3507 * (x[i] - 1.5))*(x[i] >= 1.5 && x[i] <= 3.0);
		V0[i] = 0.59 / (rho0[i] * a[i]);// 预设 U2 = 0.59;

	}

	Nozzle N1(x_f,x_e,a,rho0,T0,V0,n,it);
	N1.print();

	delete[] x;
	delete[] a;
	delete[] rho0;
	delete[] T0;
	delete[] V0;
}

int main()
{
	clock_t start = clock();
	test();
	clock_t end = clock();

	cout << "Total time: " << end - start << "ms" << endl;

	system("pause");
	return 0;
}