准一维喷管亚声超声流动 有限差分求解

2024-03-24 13:20
文章标签 超声 差分 求解 一维 流动 有限 亚声 喷管

本文主要是介绍准一维喷管亚声超声流动 有限差分求解，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

写一个准一维喷管流动的有限差分数值解，经典算例，计算流体力学入门第七章

该运动等熵，气体在喉部达到声速，后超声速流动，本例采用非守恒方程，显式MacCormark方法求解。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from matplotlib.ticker import FixedLocator

定义格子宽度Δx = 0.1 写出无量纲长度

delta_x = 0.1 x = np.arange(0.0,3.1,0.1,dtype =  float)print("无量纲距离",x)

无量纲距离 [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.  1.1 1.2 1.3 1.4 1.5 1.6 1.71.8 1.9 2.  2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. ]

定义喷管形状A=A(x)

A = 1+2.2*(x-1.5)**2  #可以看出，在x=1.5处，喉部面积为1，故A1即为无量纲面积plt.figure(figsize=(10,5))
plt.plot(x,A,linewidth=4)
plt.xlabel('Distance',labelpad=30, size=40, weight='bold')  
plt.xticks( size=30, weight='bold')
plt.yticks( size=30, weight='bold')print("各点截面积",A)

各点截面积 [5.95  5.312 4.718 4.168 3.662 3.2   2.782 2.408 2.078 1.792 1.55  1.3521.198 1.088 1.022 1.    1.022 1.088 1.198 1.352 1.55  1.792 2.078 2.4082.782 3.2   3.662 4.168 4.718 5.312 5.95 ]

在这里插入图片描述

定义t=0时各点状态

density =1 - 0.3146*x
T = 1 - 0.2314*x
V = (0.1+1.09*x)*T**(1/2)plt.figure(figsize=(10,5))
plt.plot(x,density,linewidth=4,c='gold')
plt.plot(x,T,linewidth=4,c='tan')
plt.plot(x,V,linewidth=4,c='orange')
plt.legend(['density','T','V'],frameon=False,loc="upper left",fontsize=20)
plt.xlabel('Distance',labelpad=30, size=40, weight='bold')  
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("各点密度",density)
print("各点温度",T)
print("各点速度",V)
各点密度 [1.      0.96854 0.93708 0.90562 0.87416 0.8427  0.81124 0.77978 0.748320.71686 0.6854  0.65394 0.62248 0.59102 0.55956 0.5281  0.49664 0.465180.43372 0.40226 0.3708  0.33934 0.30788 0.27642 0.24496 0.2135  0.182040.15058 0.11912 0.08766 0.0562 ]
各点温度 [1.      0.97686 0.95372 0.93058 0.90744 0.8843  0.86116 0.83802 0.814880.79174 0.7686  0.74546 0.72232 0.69918 0.67604 0.6529  0.62976 0.606620.58348 0.56034 0.5372  0.51406 0.49092 0.46778 0.44464 0.4215  0.398360.37522 0.35208 0.32894 0.3058 ]
各点速度 [0.1        0.20656772 0.31055431 0.41191227 0.5105917  0.606540110.69970225 0.79001982 0.87743124 0.96187135 1.04327104 1.121556931.19665091 1.26846965 1.33692406 1.40191865 1.46335081 1.521109941.57507649 1.62512075 1.67110158 1.71286469 1.75024077 1.783043161.81106497 1.83407564 1.85181659 1.86399585 1.87028111 1.870290821.86358258]

在这里插入图片描述

计算时间步长

C = 0.5  #定义克朗数为0.5delta_t1=[]         
for i in range(1,31):       delta_t_i = C*(delta_x/(T[i]**(1/2)+V[i]))delta_t1.append(delta_t_i)delta_t = min(delta_t1)    #取最小时间步长
print("最小时间步长为",delta_t)

最小时间步长为 0.020134450213606162

定义速度，温度，密度，压力比，马赫数，无量纲质量流量 集合

V1=[list(V)]
density1 = [list(density)]
T1 = [list(T)]
A1=list(A)

ma=[]
for i in range(0,len(V1[0])):ma.append(V1[0][i]/T1[0][i]**(1/2))
Ma=[ma]p=[]
for i in range(0,len(V1[0])):p.append(density1[0][i]*T1[0][i])
P1=[p]mm=[]
for i in range(0,len(V1[0])):mm.append(density1[0][i]*V1[0][i]*A1[i])
M=[mm]

MacCormark 显式方法

推进1400个时间步

gamma = 1.4for j in range(0,1401):ma=[]p=[]mm=[]d_den_d_t=[]d_V_d_t=[]d_T_d_t=[]_den=[]_V=[]_T=[]_d_den_d_t=[0]_d_V_d_t=[0]_d_T_d_t=[0]av_d_den_d_t=[0]av_d_V_d_t=[0]av_d_T_d_t=[0]den_next_t=[0]V_next_t=[0]T_next_t=[0]#预测步for i in range(0,30):#用向前差分计算对时间的微分d_den_d_t.append(-density1[j][i]*(V1[j][i+1]-V1[j][i])/delta_x - density1[j][i]*V1[j][i]*(np.log(A1[i+1])-np.log(A1[i]))/delta_x - V1[j][i]*(density1[j][i+1]-density1[j][i])/delta_x)d_V_d_t.append(-V1[j][i]*(V1[j][i+1]-V1[j][i])/delta_x - 1/gamma*((T1[j][i+1]-T1[j][i])/delta_x+T1[j][i]/density1[j][i]*(density1[j][i+1]-density1[j][i])/delta_x))d_T_d_t.append(-V1[j][i]*(T1[j][i+1]-T1[j][i])/delta_x - (gamma-1)*T1[j][i]*((V1[j][i+1]-V1[j][i])/delta_x+V1[j][i]*(np.log(A1[i+1])-np.log(A1[i]))/delta_x))#得到各个预测量_den.append(density1[j][i]+d_den_d_t[i]*delta_t)_V.append(V1[j][i]+d_V_d_t[i]*delta_t)_T.append(T1[j][i]+d_T_d_t[i]*delta_t)#修正步 for i in range(1,30):#使用带横线量进行向后差分_d_den_d_t.append(-_den[i]*(_V[i]-_V[i-1])/delta_x - _den[i]*_V[i]*(np.log(A1[i])-np.log(A1[i-1]))/delta_x - _V[i]*(_den[i]-_den[i-1])/delta_x)_d_V_d_t.append(-_V[i]*(_V[i]-_V[i-1])/delta_x - 1/gamma*((_T[i]-_T[i-1])/delta_x+_T[i]/_den[i]*(_den[i]-_den[i-1])/delta_x))_d_T_d_t.append(-_V[i]*(_T[i]-_T[i-1])/delta_x - (gamma-1)*_T[i]*((_V[i]-_V[i-1])/delta_x+_V[i]*(np.log(A1[i])-np.log(A1[i-1]))/delta_x))#计算各值时间导数的平均值-averageav_d_den_d_t.append(0.5*(d_den_d_t[i]+_d_den_d_t[i]))av_d_V_d_t.append(0.5*(d_V_d_t[i]+_d_V_d_t[i]))av_d_T_d_t.append(0.5*(d_T_d_t[i]+_d_T_d_t[i]))#得到下一时刻流动参数的修正值den_next_t.append(density1[j][i]+av_d_den_d_t[i]*delta_t)V_next_t.append(V1[j][i]+av_d_V_d_t[i]*delta_t)T_next_t.append(T1[j][i]+av_d_T_d_t[i]*delta_t)#计算下一时刻第一个点的速度，密度，温度den_next_t[0]=1V_next_t[0]=2*V_next_t[1]-V_next_t[2]T_next_t[0]=1#计算下一时刻最后一个点的速度，密度，温度den_next_t.append(2*den_next_t[29]-den_next_t[28])V_next_t.append(2*V_next_t[29]-V_next_t[28])T_next_t.append(2*T_next_t[29]-T_next_t[28])for i in range(len(V1[0])):ma.append(V1[j][i]/T1[j][i]**(1/2))for i in range(0,len(V1[0])):p.append(density1[j][i]*T1[j][i])for i in range(0,len(V1[0])):mm.append(density1[j][i]*V1[j][i]*A1[i])#将下一时刻各流动参量加入流动参量大集合中V1.append(V_next_t)T1.append(T_next_t)density1.append(den_next_t)Ma.append(ma)P1.append(p)M.append(mm)
t=list(range(0,1401))
D=[]
for i in range(0,1401):D.append(density1[i][15])plt.figure(figsize=(10,5))
plt.plot(t,D,linewidth=4)
plt.xlabel('TIME(Δt)',labelpad=30, size=40, weight='bold')  
plt.ylabel('$\mathregular{ρ/ρ_0}$',labelpad=30, size=40, weight='bold') 
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("喷管喉部密度随时间步变化")

喷管喉部密度随时间步变化

在这里插入图片描述

t=list(range(0,1401))
TTT=[]
for i in range(0,1401):TTT.append(T1[i][15])plt.figure(figsize=(10,5))
plt.plot(t,TTT,linewidth=4)
plt.xlabel('TIME(Δt)',labelpad=30, size=40, weight='bold')  
plt.ylabel('$\mathregular{T/T_0}$',labelpad=30, size=40, weight='bold') 
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("喷管喉部温度随时间步变化")

喷管喉部温度随时间步变化

在这里插入图片描述

t=list(range(0,1401))
MMA=[]
for i in range(0,1401):MMA.append(Ma[i][15])plt.figure(figsize=(10,5))
plt.plot(t,MMA,linewidth=4)
plt.xlabel('TIME(Δt)',labelpad=30, size=40, weight='bold')  
plt.ylabel('Ma',labelpad=30, size=40, weight='bold') 
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("喷管喉部马赫数随时间步变化")

喷管喉部马赫数随时间步变化

在这里插入图片描述

t=list(range(0,1401))
PPP=[]
for i in range(0,1401):PPP.append(P1[i][15])plt.figure(figsize=(10,5))
plt.plot(t,PPP,linewidth=4)
plt.xlabel('TIME(Δt)',labelpad=30, size=40, weight='bold')  
plt.ylabel('$\mathregular{P/P_0}$',labelpad=30, size=40, weight='bold') 
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("喷管喉部压力比马赫数随时间步变化")

喷管喉部压力比马赫数随时间步变化

在这里插入图片描述

plt.figure(figsize=(10,10))
plt.plot(x,Ma[1400],linewidth=4)
plt.scatter(x,Ma[1400],marker = "^",linewidths=6)
plt.xlabel('DISTANCE (x)',labelpad=30, size=40, weight='bold')  
plt.ylabel('Ma',labelpad=30, size=40, weight='bold')  
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("数值解：喷管内马赫数随随无量纲距离x的变化")

数值解：喷管内马赫数随随无量纲距离x的变化

在这里插入图片描述

plt.figure(figsize=(10,10))
plt.plot(x,density1[1400],linewidth=4)
plt.scatter(x,density1[1400],marker = "^",linewidths=6)
plt.xlabel('DISTANCE (x)',labelpad=30, size=40, weight='bold')  
plt.ylabel('$\mathregular{ρ/ρ_0}$',labelpad=30, size=40, weight='bold')  
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("数值解：喷管内无量纲密度随无量纲距离x的变化")

数值解：喷管内无量纲密度随无量纲距离x的变化

在这里插入图片描述

plt.figure(figsize=(8,10))
plt.plot(x,M[0],linewidth=4,c='#D9CBDF')
plt.plot(x,M[50],linewidth=4,c='#C19CC4')
plt.plot(x,M[100],linewidth=4,c='#FAE9BD')
plt.plot(x,M[150],linewidth=4,c='#F4D6CE')
plt.plot(x,M[200],linewidth=4,c='#E3E2E1')
plt.plot(x,M[700],linewidth=4,c='#C1BDBC')
#plt.scatter(x,density1[1400],marker = "^",linewidths=6)
plt.legend(['0Δt','50Δt','100Δt','150Δt','200Δt','700Δt'],frameon=False,loc="upper right",fontsize=20)
plt.xlabel('DISTANCE (x)',labelpad=30, size=40, weight='bold')  
plt.ylabel('$\mathregular{ρVA/{ρ_0}{a_0}{A^*}}$',labelpad=30, size=40, weight='bold')  
plt.xticks( size=25, weight='bold')
plt.yticks( size=25, weight='bold')print("数值解：喷管内无量纲密度随无量纲距离x的变化")

数值解：喷管内无量纲密度随无量纲距离x的变化

在这里插入图片描述

给出1400步以后的流场参数

data_1400 = {'x':x, 'A':A, 'ρ':density1[1400], 'V':V1[1400], 'T':T1[1400], 'P':P1[1400], 'Ma':Ma[1400],'m':M[1400]}df_1400 = pd.DataFrame(data_1400)
df_1400
xAρVTPMam
00.05.9501.0000000.0986141.0000001.0000000.0986140.586751
10.15.3120.9975350.1119350.9990100.9965470.1119910.593135
20.24.7180.9966730.1252570.9986680.9953450.1253410.588998
30.34.1680.9939750.1425780.9975830.9915730.1427500.590683
40.43.6620.9913140.1622660.9965170.9878600.1625490.589057
50.53.2000.9870270.1865330.9947900.9818850.1870210.589163
60.62.7820.9814900.2154330.9925570.9741850.2162390.588241
70.72.4080.9734700.2507490.9893070.9630600.2521010.587785
80.82.0780.9623040.2935560.9847590.9476380.2958190.587014
90.91.7920.9463550.3457600.9782090.9257330.3495900.586364
101.01.5500.9238240.4089980.9688500.8950470.4155210.585655
111.11.3520.8921970.4850040.9554880.8524840.4961720.585036
121.21.1980.8488890.5747380.9367400.7951880.5938270.584490
131.31.0880.7917980.6779720.9111470.7214440.7102600.584056
141.41.0220.7207190.7924840.8776880.6325670.8459030.583724
151.51.0000.6384590.9139580.8363430.5339710.9993880.583525
161.61.0220.5508221.0365920.7884940.4343201.1673710.583540
171.71.0880.4648061.1545370.7367000.3424231.3451250.583860
181.81.1980.3862021.2633250.6839260.2641341.5276030.584503
191.91.3520.3182541.3604950.6327440.2013731.7103420.585393
202.01.5500.2617331.4454830.5848670.1530791.8900990.586412
212.11.7920.2158351.5188600.5411910.1168082.0646310.587458
222.22.0780.1790041.5819690.5019220.0898462.2329560.588445
232.32.4080.1495981.6360700.4669570.0698562.3942160.589363
242.42.7820.1260511.6828400.4358620.0549412.5489920.590130
252.53.2000.1071721.7229890.4083450.0437632.6963030.590901
262.63.6620.0918291.7584550.3837480.0352392.8386280.591327
272.74.1680.0794501.7885590.3620600.0287662.9724380.592280
282.84.7180.0690201.8166030.3422330.0236213.1052660.591549
292.95.3120.0609351.8391380.3252390.0198183.2248770.595301
303.05.9500.0528501.8616740.3082450.0162913.3531730.585412

参考资料：
[1] 计算流体力学基础及其应用. [美]John D.Anderson. 著；吴颂平, 刘赵淼译. 北京：机械工业出版社，2007. 6

这篇关于准一维喷管亚声超声流动 有限差分求解的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/841702。

相关文章

2024年流动式起重机司机证模拟考试题库及流动式起重机司机理论考试试题

2024年流动式起重机司机证模拟考试题库及流动式起重机司机理论考试试题

题库来源：安全生产模拟考试一点通公众号小程序 2024年流动式起重机司机证模拟考试题库及流动式起重机司机理论考试试题是由安全生产模拟考试一点通提供，流动式起重机司机证模拟考试题库是根据流动式起重机司机最新版教材，流动式起重机司机大纲整理而成（含2024年流动式起重机司机证模拟考试题库及流动式起重机司机理论考试试题参考答案和部分工种参考解析），掌握本资料和学校方法，考试容易。流动式起重机司机考试技
阅读更多...
poj 3159 （spfa差分约束最短路） poj 1201

poj 3159 （spfa差分约束最短路） poj 1201

poj 3159: 题意： 每次给出b比a多不多于c个糖果，求n最多比1多多少个糖果。 解析： 差分约束。 这个博客讲差分约束讲的比较好： http://www.cnblogs.com/void/archive/2011/08/26/2153928.html 套个spfa。 代码： #include <iostream>#include <cstdio>#i
阅读更多...
poj 3169 spfa 差分约束

poj 3169 spfa 差分约束

题意： 给n只牛，这些牛有些关系。 ml个关系：fr 与 to 牛间的距离要小于等于 cost。 md个关系：fr 与 to 牛间的距离要大于等于 cost。 隐含关系： d[ i ] <= d[ i + 1 ] 解析： 用以上关系建图，求1-n间最短路即可。 新学了一种建图的方法。。。。。。 代码： #include <iostream>#include
阅读更多...
POJ 1364差分约束

POJ 1364差分约束

给出n个变量，m个约束公式 Sa + Sa+1 + .... + Sa+b < ki or > ki ,叫你判断是否存在着解满足这m组约束公式。 Sa + Sa+1   +   .+ Sa+b =  Sum[a+b] - Sum[a-1]  . 注意加入源点n+1 。 public class Main {public static void main(Strin
阅读更多...
简单的Q-learning|小明的一维世界(3)

简单的Q-learning|小明的一维世界(3)

简单的Q-learning|小明的一维世界(1) 简单的Q-learning|小明的一维世界(2) 一维的加速度世界 这个世界，小明只能控制自己的加速度，并且只能对加速度进行如下三种操作：增加1、减少1、或者不变。所以行动空间为： { u 1 = − 1 , u 2 = 0 , u 3 = 1 } \{u_1=-1, u_2=0, u_3=1\} {u1​=−1,u2​=0,u3​=1}
阅读更多...
简单的Q-learning|小明的一维世界(2)

简单的Q-learning|小明的一维世界(2)

上篇介绍了小明的一维世界模型 、Q-learning的状态空间、行动空间、奖励函数、Q-table、Q table更新公式、以及从Q值导出策略的公式等。最后给出最简单的一维位置世界的Q-learning例子，从给出其状态空间、行动空间、以及稠密与稀疏两种奖励函数的设置方式。下面将继续深入，GO！ 一维的速度世界 这个世界，小明只能控制自己的速度，并且只能对速度进行如下三种操作：增加1、减
阅读更多...
2024 年高教社杯全国大学生数学建模竞赛题目——2024 年高教社杯全国大学生数学建模竞赛题目的求解

2024 年高教社杯全国大学生数学建模竞赛题目——2024 年高教社杯全国大学生数学建模竞赛题目的求解

2024 年高教社杯全国大学生数学建模竞赛题目 （请先阅读“ 全国大学生数学建模竞赛论文格式规范 ”） 2024 年高教社杯全国大学生数学建模竞赛题目 随着城市化进程的加快、机动车的快速普及， 以及人们活动范围的不断扩大，城市道 路交通拥堵问题日渐严重，即使在一些非中心城市，道路交通拥堵问题也成为影响地方经 济发展和百姓幸福感的一个“痛点”，是相关部门的棘手难题之一。 考虑一个拥有知名景区
阅读更多...
Python中差分进化differential_evolution的调用及参数说明

Python中差分进化differential_evolution的调用及参数说明

在场景应用中，要求我们的函数计算结果尽可能的逼近实际测量结果，可转化计算结果与测量结果的残差，通过最小化残差，便可求出最优的结果。但使用最小二乘等方法来计算时，常常会使迭代的结果显然局部最优点而导致结算错误。 差分进化原理 差分进化（Differential Evolution，DE）是一种基于群体差异的进化算法，其计算思想主要包括以下几个方面： 一、初始化种群 首先，随机生成一个初始种群
阅读更多...
正规式与有限自动机例题

正规式与有限自动机例题

答案：D 知识点： 正规式 正规集 举例 ab 字符串ab构成的集合 {ab} a|b 字符串a,b构成的集合 {a,b} a^* 由0或者多个a构成的字符串集合 {空,a,aa,aaa,aaaa····} (a|b)^* 所有字符a和b构成的串的集合 {空,a,b,ab,aab,aba,aaab····} a(a|b)^* 以a为首字符的a,b字符串的集
阅读更多...
c语言——用一维数组输出杨辉三角形

c语言——用一维数组输出杨辉三角形

一.代码 #include <stdio.h>int Num[100];int Hang;int Lie;int a;int Flag;int main() {Lie = 1;Hang = 1;a = 0;while (1) {//列1为1if (Lie == 1) {Num[1] = 1;Lie++;}//数据存到数组里面while (Hang >= Lie && Hang !=
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2