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

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。

相关文章

将一维机械振动信号构造为训练集和测试集（Python）

将一维机械振动信号构造为训练集和测试集（Python）

从如下链接中下载轴承数据集。 https://www.sciencedirect.com/science/article/pii/S2352340918314124 import numpy as npimport scipy.io as sioimport matplotlib.pyplot as pltimport statistics as statsimport pandas
阅读更多...
鹅算法（GOOSE Algorithm，GOOSE)求解复杂城市地形下无人机避障三维航迹规划，可以修改障碍物及起始点（Matlab代码）

鹅算法（GOOSE Algorithm，GOOSE)求解复杂城市地形下无人机避障三维航迹规划，可以修改障碍物及起始点（Matlab代码）

一、鹅算法 鹅优化算法（GOOSE Algorithm，GOOSE)从鹅的休息和觅食行为获得灵感，当鹅听到任何奇怪的声音或动作时，它们会发出响亮的声音来唤醒群中的个体，并保证它们的安全。 参考文献 [1]Hamad R K, Rashid T A. GOOSE algorithm: a powerful optimization tool for real-world engineering
阅读更多...
SQL求解两个时间差 时间类型 时间值

SQL求解两个时间差 时间类型 时间值

sql 求解两个时间差 SELECTDATEDIFF( Second, '2009-8-25 12:15:12', '2009-9-1 7:18:20') --返回相差秒数 SELECTDATEDIFF( Minute, '2009-9-1 6:15:12', '2009-9-1 7:18:20') --返回相差分钟数 SELECTDATEDIFF( Day, '2009-8
阅读更多...
编写程序，采用辗转相除法求解两个正整数的最大公约数

编写程序，采用辗转相除法求解两个正整数的最大公约数

--编写程序，采用辗转相除法求解两个正整数的最大公约数DECLARE @a int,@b intSELECT @a=12,@b=21DECLARE @temp intprint cast(@a as varchar(5))+'和'+cast(@b as varchar(5))+'的最大公约数是'if @a<@b --或者是select @temp=@a,@a=@b,@b=@tempb
阅读更多...
软件开发-人员流动

软件开发-人员流动

软件开发 并不是一件 简单的事情。 如果程序员像白菜一样，标一个价，放在市场卖的话， 那么，软件，也不是，一堆白菜堆在一起就ok了。 人员在一个地方待久了，可能就会流动。 毕竟世界这么大。 可是，软件系统 一旦开发出来， 应该就在线上一直 运行着。 那么人员的流动，决不应该像白菜一样，一个挨一个的放着。他们彼此重叠，这样的重叠部分，应该不少于半年。这个软件系统才是可控的。 否则， 这个拼
阅读更多...
维度/标量/张量/一维/二维/三维/shape/size/索引/view理解

维度/标量/张量/一维/二维/三维/shape/size/索引/view理解

维度是对不类型的描述,有几个不同的类型就有几个不同的维度. 张量是用来描述维度更宽泛的概念 标量是一个数. 一维是一组数据. 二维是两个不同特征组合的数据. 三维是三个不同的特征组合的数据.[] 只是用来表示区间范围套更多的[]是为了表示不同的维度,要理解维度 含义而不是只看[]. 生活例子 假设我们有一个书架，书架有很多层，每层放着书： 索引 (Index)：就像书架上书的位置，比
阅读更多...
长尾式差分放大电路调零

长尾式差分放大电路调零

长尾式放大电路用了两个参数相同的三极管，但实际上并没有完全相同的三极管，所以为了提高差分放大电路的对称性(一边电流增加多少，另一边电流减小多少，即能在电阻Re上产生的压降不变(后面做虚地处理))，在下图中加入可调电阻，调节可调电阻的值，便可使输入为零时 输出也为零。 可调电阻尽量选小一些:①过大的可调电阻会影响动态的放大倍数，②在选三极管时选的是参数接近的三极管，所以用小的可调电阻微调即可。
阅读更多...
Java程序设计 第七章 一维数组

Java程序设计 第七章 一维数组

目录 7.2 数组基础 7.3 示例学习：分析数字 7.5复制数组 7.6 将数组传递给方法 7.7 方法返回数组 7.9 可变长参数列表 7.10 查找数组 7.11 数组排序（直接排序） 略 7.12 Arrays类 7.13 命令行参数 7.2 数组基础   声明数组：elementType[]  arrayRefVar;                 例
阅读更多...
【工具笔记】Microsoft数学求解器Math Solver

【工具笔记】Microsoft数学求解器Math Solver

【工具笔记】Microsoft数学求解器Math Solver 工具笔记用于记录各种有用的工具，这里记录的是一个由Microsoft提供的数学求解器Math Solver。 可以用于求解代数，三角学，微积分，矩阵等各种数学问题，并且可以获取分步解释，查看如何解决问题并获取数学概念的定义，立即画出任何公式以可视化函数并了解变量之间的关系。还会搜索出相关的视频，练习题，类似问题等。 Math S
阅读更多...
编译原理-各章典型题型+思路求解

编译原理-各章典型题型+思路求解

第2章文法和语言习题 基础知识： 思路： 基础知识： 思路： 基础知识： 编译原理之 短语&直接短语&句柄 定义与区分_编译原理短语,直接短语,句柄-CSDN博客 思路： 题目： 基础解释：   简单来说： 上下文无关文法就是这个文法中所有的产生式左边只有一个非终结： 例如：S->Abc 上下文无关文法就是第一个产生式左
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2