2021-04-26

西西西西柚

于 2021-04-26 21:23:26 发布

阅读量547

点赞数 1

分类专栏：计算流体力学 CFD 安德森文章标签：算法

本文链接：https://blog.csdn.net/ErnestChai/article/details/116173566

版权

计算流体力学同时被 3 个专栏收录

7 篇文章 3 订阅

订阅专栏

CFD

7 篇文章 0 订阅

订阅专栏

安德森

7 篇文章 0 订阅

订阅专栏

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

激波捕捉：拟一维扩张段存在正激波的拉伐尔喷管流动（守恒型方程）

本文主要分享了安德森的《计算流体力学基础与应用》一书中：“拟一维扩张段存在正激波的拉伐尔喷管流动” 的程序求解，分别给出了添加了人工粘性和未添加人工粘性的不同求解程序。水平有限，请见谅！

计算结果

喷管内的压力分布，1600个时间步
喷管内的马赫数分布
喷管内的质量流量分布

Matlab

Nozzle.m

%%
% Author: CXT
% Date: 2021/4/26
% Theme: 拟一维扩张段存在正激波的拉伐尔喷管流动 - 守恒型求解
%% 物理符号(无量纲量)；
% x —— 空间位置；
% A —— 喷管面积；
% rho —— 密度；
% V —— 速度；
% p —— 压力；
% N —— 计算节点数；
% m —— 质量流量；
%%
clc;
clear;
close all;
%% 喷管形状;
N = 61;% 计算节点数；
x = 0:(3-0)/(N-1):3; % 计算区间；
n = 1600; %时间步数；

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.702*(x - 1.5)).*(x>=1.5 & x<2.1)...
     + (0.5892 - 0.10228*(x - 2.1)).*(x>=2.1 & x<=3.0);
T0 = 1.*(x<=0.5)...
   + (1.0 - 0.167*(x - 0.5)).*(x>0.5 & x<1.5) ...
   + (0.833 - 0.4908*(x - 1.5)).*(x>=1.5 & x<2.1)...
   + (0.93968 - 0.0622*(x - 2.1)).*(x>=2.1 & x<=3.0);
V0 = 0.59./(rho0.*A);% 预设 U2 = 0.59; 

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

tic;
% 添加了人口粘性；
for i = 1:n
    [rho_t, V_t, T_t, p_t, Ma, m] = OneStepNozzleWithArtiVisc(rho_t, V_t, T_t, N, x, A);
end

figure;
plot(x, p_t,'--k');
hold on;

% 未添加人工粘性；
for i = 1:n
    [rho_t, V_t, T_t, p_t, Ma, m] = OneStepNozzleWithoutArtiVisc(rho_t, V_t, T_t, N, x, A);
end

plot(x, p_t,'-r');
title("喷管内的压力分布，1600个时间步");
legend("添加了人口粘性","未添加人工粘性");
xlabel('x/L');
ylabel("无量纲压力");
hold off;
toc;

figure;
plot(x,Ma);
xlabel('x/L'),ylabel('Ma');
title("喷管内的马赫数分布");

figure;
plot(x,m);
xlabel('x/L'),ylabel('\rho AV/\rho_0 A*a_0');
title("喷管内的质量流量分布")

OneStepNozzleWithArtiVisc.m

function [rho_t, V_t, T_t, p_t, Ma, m] = OneStepNozzleWithArtiVisc(rho, V, T, N, X, A)
%[rho_t, V_t, T_t, p_t, Ma] = OneStepNozzle(rho, V, T) 拟一维扩张段存在正激波的拉伐尔喷管流动 - 守恒型求解
% 添加了人工粘性（Artificial Viscosity）
%   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; %柯朗数；
C_x = 0.2; %人工粘性参数；
%% 解向量 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;

%% 计算初始压力 p，用于计算人工粘性;
p = rho.*T;

%% 预估步计算；
% 向前差分；
% 只计算内点；
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的预估值；
% 添加人工粘性S1,S2和S3；
S1 = zeros(1,N);
S2 = zeros(1,N);
S3 = zeros(1,N);
S = zeros(1,N); % S 为中间变量；

for i = 2:N-1
    S(i) = C_x * abs(p(i+1) - 2*p(i) + p(i-1))/(p(i+1) + 2*p(i) + p(i-1));
    
    S1(i) = S(i) * (U1(i+1) - 2*U1(i) + U1(i-1));
    S2(i) = S(i) * (U2(i+1) - 2*U2(i) + U2(i-1));
    S3(i) = S(i) * (U3(i+1) - 2*U3(i) + U3(i-1));
end

bar_U1 = U1 + partial_U1 * Delta_t + S1;
bar_U2 = U2 + partial_U2 * Delta_t + S2;
bar_U3 = U3 + partial_U3 * Delta_t + S3;

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

% 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 的校正值；
% 添加人工粘性bar_S1、bar_S2 和 bar_S3；
bar_S1 = zeros(1,N);
bar_S2 = zeros(1,N);
bar_S3 = zeros(1,N);

for i = 2:N-1
    S(i) = C_x * abs(bar_p(i+1) - 2*bar_p(i) + bar_p(i-1))/...
                (bar_p(i+1) + 2*bar_p(i) + bar_p(i-1));
    
    bar_S1(i) = S(i) * (bar_U1(i+1) - 2*bar_U1(i) + bar_U1(i-1));
    bar_S2(i) = S(i) * (bar_U2(i+1) - 2*bar_U2(i) + bar_U2(i-1));
    bar_S3(i) = S(i) * (bar_U3(i+1) - 2*bar_U3(i) + bar_U3(i-1));
end

U1_t = U1 + partial_U1_av * Delta_t + bar_S1;
U2_t = U2 + partial_U2_av * Delta_t + bar_S2;
U3_t = U3 + partial_U3_av * Delta_t + bar_S3;

% 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 和 U2 是允许变化的；
% 解向量U1、U2 和 U3；
U1_t(N) = 2*U1_t(N-1) - U1_t(N-2);
U2_t(N) = 2*U2_t(N-1) - U2_t(N-2);
V_t(N) = U2_t(N)/U1_t(N);
U3_t(N) = 0.6784*A(N)/(gama - 1) + gama/2*U2_t(N)*V_t(N);

% 原始变量；
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

OneStepNozzleWithoutArtiVisc.m

function [rho_t, V_t, T_t, p_t, Ma, m] = OneStepNozzleWithoutArtiVisc(rho, V, T, N, X, A)
%[rho_t, V_t, T_t, p_t, Ma] = OneStepNozzle(rho, V, T) 拟一维扩张段存在正激波的拉伐尔喷管流动 - 守恒型求解
%  未添加人工粘性（Artificial Viscosity）
%   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;

%% 计算初始压力 p，用于计算人工粘性;
p = rho.*T;

%% 预估步计算；
% 向前差分；
% 只计算内点；
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 和 p 的预估值；
bar_rho = bar_U1./A;
bar_T = (gama - 1) * (bar_U3./bar_U1 - gama/2 * (bar_U2./bar_U1).^2);
bar_p = bar_rho.*bar_T;

% 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 和 U2 是允许变化的；
% 解向量U1、U2 和 U3；
U1_t(N) = 2*U1_t(N-1) - U1_t(N-2);
U2_t(N) = 2*U2_t(N-1) - U2_t(N-2);
V_t(N) = U2_t(N)/U1_t(N);
U3_t(N) = 0.6784*A(N)/(gama - 1) + gama/2*U2_t(N)*V_t(N);

% 原始变量；
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