亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）

入门CFD，主要参考书目《计算流体力学基础及其应用》（John D.Anderson 著，吴颂平等译）

实现了 第 7.3 节亚声速-超声速等熵喷管流动的CFD解法 的代码，采用的是 MacCormack 方法。代码中增加了动态显示无量纲温度和无量纲密度的功能，参考的是 python中plot实现即时数据动态显示方法，最终结果如下：图中红点代表喉口所在位置。

不足之处，欢迎指正！

完整代码如下：

# -*- coding: utf-8 -*-
"""
Created on Wed Aug  5 19:11:15 2020@author: PANSONG
"""#### Laval 喷管 拟一维计算 ###############import numpy as np
import matplotlib.pyplot as plt#%% 参数，采用无量纲量
#喷管截面形状
def section(x):return 1+2.2*(x-1.5)**2# 初始条件
def initial_condition(x):rho = 1-0.314*xT = 1-0.2314*xV = (0.1+1.09*x)*np.sqrt(T)return rho,T,V# 理想气体比热比
gamma = 1.4## 离散化
max_x = 3
dx = 0.1
num_point = int(max_x/dx + 1)
x = np.linspace(0,3,num_point)max_Iteration = 10000#%% 动态绘图
plt.ion() #开启interactive mode 成功的关键函数
plt.figure()
plt.title('Dimensionless Parameters Evolution')
plt.xlabel('Iteration')iterations = []#%%  CFD 计算
#### 初始化，无量纲量
A = section(x)## 这里给两个初值，一个 0，一个线性，避免迭代过程中 残差判断时 重复判断 i
array0 = np.ones(shape=x.shape)*0.01 rho,T,V = initial_condition(x)
rho_history = np.vstack([array0,rho])
T_history = np.vstack([array0,T])
V_history = np.vstack([array0,V])####### MacCormack 方法 ###################################
# 向前，向后差分计算，输入 N 个数据，输出 N-2 个偏导数
def forward_partial(y,dx):return (y[2:]-y[1:-1])/dxdef backward_partial(y,dx):return (y[1:-1]-y[:-2])/dxdef new_rho_T(original,partial,dt):new = original + np.hstack([0,partial,0])*dt # P_rho_t 首尾补 0 new[-1] = 2*new[-2]-new[-3]  # 出口处又外插值确定，入口处保持不变return newdef new_V(V,P_V_t,dt):new_V = V + np.hstack([0,P_V_t,0])*dtnew_V[0] = 2*new_V[1] - new_V[2] #入口的速度可变，外插值确定new_V[-1] = 2*new_V[-2]-new_V[-3]return new_VC = 1 # 柯朗数,小于等于 1;从精度考虑，尽可能接近 1
# 关于柯朗数，通过对线性双曲型偏微分方程的稳定性分析，要求 C<=1，对应数值解的依赖区域包含精确解依赖区域
# 从解的精度考虑，应使数值解的依赖区域尽可能等于精确解依赖区域，即 C 尽可能接近 1
# 对于非线性偏微分方程，起指导性作用，参见参考书目 P221for i in range(max_Iteration):# 绘制 喉口处 无量纲温度 和 无量纲密度iterations.append(i)plt.plot(iterations,rho_history[1:,15],'r-',linewidth=1.0)plt.plot(iterations,T_history[1:,15],'b-',linewidth=1.0)plt.legend(['Density','Temperature'])plt.pause(1e-3)##########################################   rho_error = max(np.abs((rho_history[i]-rho_history[i-1])/rho_history[i-1]))V_error = max(np.abs((V_history[i]-V_history[i-1])/V_history[i-1]))T_error = max(np.abs((T_history[i]-T_history[i-1])/T_history[i-1]))error = max(rho_error,V_error,T_error) ## 取最大相对误差if error < 1e-4:breaka = np.sqrt(T) #无量纲当地声速        dt = min(C*dx/(V+a)) # 计算最大允许时间推进步长   ###### 预估步 ########P_rho_t = -V[1:-1]*forward_partial(rho, dx) - rho[1:-1]*forward_partial(V, dx) - rho[1:-1]*V[1:-1]*forward_partial(np.log(A), dx)rho_pred = new_rho_T(rho,P_rho_t,dt)P_V_t = -V[1:-1]*forward_partial(V, dx) - (1.0/gamma)*(forward_partial(T, dx)+(T[1:-1]/rho[1:-1])*forward_partial(rho, dx))V_pred = new_V(V,P_V_t,dt)P_T_t = -V[1:-1]*forward_partial(T, dx) - (gamma-1.0)*T[1:-1]*(forward_partial(V, dx)+V[1:-1]*forward_partial(np.log(A), dx))T_pred = new_rho_T(T,P_T_t,dt)###### 校正步 ########P_rho_t_pred = ( -V_pred[1:-1]*backward_partial(rho_pred, dx) -rho_pred[1:-1]*backward_partial(V_pred, dx) -rho_pred[1:-1]*V_pred[1:-1]*backward_partial(np.log(A), dx) )P_rho_t_av = 0.5*( P_rho_t_pred + P_rho_t )P_V_t_pred = ( -V_pred[1:-1]*backward_partial(V_pred, dx)-1.0/gamma*backward_partial(T_pred, dx)-1.0/gamma*T_pred[1:-1]/rho_pred[1:-1]*backward_partial(rho_pred, dx) )P_V_t_av = 0.5*( P_V_t_pred + P_V_t )P_T_t_pred = ( -V_pred[1:-1]*backward_partial(T_pred, dx)-(gamma-1.0)*T_pred[1:-1]*backward_partial(V_pred, dx)-(gamma-1.0)*T_pred[1:-1]*V_pred[1:-1]*backward_partial(np.log(A), dx) )P_T_t_av = 0.5*( P_T_t_pred + P_T_t )rho = new_rho_T(rho,P_rho_t_av,dt)rho_history = np.vstack([rho_history,rho])V = new_V(V,P_V_t_av,dt)V_history = np.vstack([V_history,V])T = new_rho_T(T,P_T_t_av,dt)T_history = np.vstack([T_history,T])plt.ioff() #关闭动态绘图#%% 计算结果可视化
def plot_results(x,y,x0,y0,title='value'):plt.figure()plt.plot(x,y)plt.scatter(x0,y0,color='red')plt.xlabel('Dimensionless distance')plt.title('Steady '+title+' in different position')plot_results(x,rho_history[-1,:],1.5,rho_history[-1,15],title='density')
plot_results(x,T_history[-1,:],1.5,T_history[-1,15],title='temperature')p = rho_history[-1,:]*T_history[-1,:]
plot_results(x,p,1.5,p[15],title='pression')Ma = V_history[-1,:]/np.sqrt(T_history[-1,:])
plot_results(x,Ma,1.5,Ma[15],title='Mach number')plt.show()

亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）相关推荐

拟一维喷管流动的数值解——亚声速-超声速等熵喷管流动的非守恒型CFD解法（MacCormack方法）
一.Matlab代码片 %亚声速-超声速等熵喷管流动非守恒型麦考马克方法数值求解 clear; %清理内存变量 clc; %清理工作窗中的所有显示内容 r=1.4; %比热比 L=3; %喷管长度 ...
拟一维喷管流动的数值解——亚声速-超声速等熵喷管流动的守恒型CFD解法（MacCormack方法）
一.Matlab代码片 %亚声速-超声速等熵喷管流动守恒形CFD解法 MacCormack方法 clear; %清理内存变量 clc; %清理工作窗中的所有显示内容 r=1.4; %比热比 L=3; ...
亚声速 – 超声速等熵喷管流动数值模拟(文字）
亚声速 – 超声速等熵喷管流动问题描述流过拉伐尔喷管的定常等熵流动在喷管,流动经过等熵地加速已达到超声速喷管出口处流体的压力.温度.速度和马赫数分别记为在喷管收缩段,流动是亚声速的:在喷管的 ...
『迷你教程』识别人类活动的一维卷积神经网络模型，附完整代码
文章目录使用智能手机数据集进行活动识别开发一维卷积神经网络拟合和评估模型总结结果调整一维卷积神经网络过滤器数量内核大小多头卷积神经网络一些拓展的想法人类活动识别是将专用线束或智能手 ...
喷管流动的守恒型CFD解法及激波捕捉（附完整代码）
入门CFD,主要参考书目<计算流体力学基础及其应用>(John D.Anderson 著,吴颂平等译) 实现了第 7.6 节激波捕捉的代码,采用的是 MacCormack 方法, ...
流过平板的超声速流动的CFD计算（附完整代码）
入门CFD,主要参考书目<计算流体力学基础及其应用>(John D.Anderson 著,吴颂平等译) 实现了第 10 章流过平板的超声速流动的代码.利用有限差分法求解二维 Nav ...
【图像分割】基于布谷鸟算法实现二维Tsallis熵、kapur、oust多阈值图像分割附matlab代码
1 内容介绍本文介绍了一种基于布谷鸟算法的多级阈值(MT)算法.布谷鸟优化算法[CuckooSearch (CS)],也叫杜鹃搜索,是智能算法的其中一种,于2009年由剑桥大学Xin-SheYang ...
Matlab 编程《计算流体力学基础及应用（约翰D安德森）》全亚声速等熵喷管流动CFD解法拉瓦尔喷管非守恒形式方程解法
Matlab 编程 <计算流体力学基础及应用(约翰D安德森)> 全亚声速等熵喷管流动CFD解法拉瓦尔喷管非守恒形式方程解法问题之全亚声速等熵喷管流动的CFD解法初始化参数迭代过 ...
拟一维喷管流动的数值解——全亚声速等熵喷管流动的非守恒型CFD解法（MacCormack方法）
一.Matlab代码片 %全亚声速等熵喷管流动非守恒型麦考马克方法数值求解 clear; %清理内存变量 clc; %清理工作窗中的所有显示内容 r=1.4; %比热比 L=3; %喷管长度 i=3 ...

亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）

亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）相关推荐

最新文章

热门文章