import math
import decimal
import numpy as np
from decimal import *
wnx=31#最后一个网格记数31
dx=3/(wnx-1)#网格的大小
bzb=1.4#比热比
dtb=200
C = 0.5
at=0.0201
x=np.zeros((wnx,1))
A=np.zeros((wnx,1))
A1=np.zeros((wnx,1))
midu=np.zeros((wnx,dtb+1))
T=np.zeros((wnx,dtb+1))
V=np.zeros((wnx,dtb+1))
t=np.zeros((wnx,dtb+1))
dt=np.zeros((1,wnx+1))
midut=np.zeros((wnx,dtb+1))
Tt=np.zeros((wnx,dtb+1))
Vt=np.zeros((wnx,dtb+1))
midut1=np.zeros((wnx,dtb+1))
Tt1=np.zeros((wnx,dtb+1))
Vt1=np.zeros((wnx,dtb+1))
dtbshu=np.zeros((1,dtb))
y1=np.zeros((wnx,dtb+1))
ma=np.zeros((wnx,dtb+1))
for i in range(0,wnx):
x[i]=round(3*(i)/(wnx-1),3)
A[i]= round(float(1+2.2*(x[i]-1.5)**2),5)
midu[i,0]=round(float(1-0.3146x[i]),3)#出始密度
T[i,0]=round(float(1-0.2314x[i]),3)
V[i,0]=round(float(0.1+1.09*x[i])*math.sqrt(T[i,0]),3)

print(V[15,0]),print(V[16,0]),print(A[15]),print(A[16])

print(T[15,0]),print(T[16,0])

print(midu[15,0]),print(midu[16,0])

print(“x方向”,[i[0]for i in x])
print(“面积”,[i[0]for i in A])
print(“密度”,[i[0]for i in midu])
print(“速度”,[i[0]for i in V])
print(“能量”,[i[0]for i in T])
print(V[16,0])

print(“压力”,[i[0]for i in y1])

print(“马赫数”,[i[0]for i in ma])

midu[0,0]=1
T[0,0]=1
V[0,0]=0.1

midu[wnx-1,0]=2*midu[wnx-2,0]-midu[wnx-3, 0]

T[wnx-1,0] = 2*T[wnx-2,0]-T[wnx-3,0]

V[wnx-1,0] = 2*V[wnx-2,0]-V[wnx-3,0 ]

for i in range(0,30):
A1[i]=math.log(A[i])
for i in range(1,30):
midut[i,0] =round( -V[i,0] * (midu[i+1,0]-midu[i,0])/dx-midu[i,0](V[i+1,0]-V[i,0])/dx-midu[i,0]V[i,0](math.log(A[i+1])-math.log(A[i]))/dx,4)
Vt[i,0]=round(-V[i,0](V[i+1,0]-V[i,0])/dx-((T[i+1,0]-T[i,0])/dx+T[i,0]*(midu[i+1,0]-midu[i,0])/(midu[i,0]*dx))/bzb,4)
# Tt[i, 0] = round(-V[i, 0] * (T[i + 1, 0] - T[i, 0]) / dx - (bzb - 1) * T[i, 0] * (
# (V[i + 1, 0] - V[i, 0]) / dx + V[i, 0] * (math.log(A[i + 1]) - math.log(A[i])) / dx), 4)
ss=A1[i + 1] - A1[i]
Tt[i, 0] =np.round(-V[i, 0] * (T[i + 1, 0] - T[i, 0]) / dx - (bzb - 1) * T[i, 0] * ((V[i + 1, 0] - V[i, 0]) / dx + V[i, 0] *ss / dx),3)

print(midut[15,0]),print(Vt[15,0]),print(Tt[15,0])
print([i[0]for i in midut])
print([i[0]for i in Vt])
print([i[0]for i in Tt])

s=len(T)
print(s)
# 各网格点∆t的最小值
# for i in range(0, 30):
# t [i] = C * (dx / (V[i,j] + math.sqrt(T[i,j])))
# at = min(t[i]) # 找出了各网格点∆t的最小值并赋值给at
# print(at)
for i in range(30):
midu[i,0] = round(midu[i,0] + midut[i,0]at,3)
V[i,0]= round(V[i,0] + Vt[i,0] * at,3)
T[i,0]=round(T[i,0] + Tt[i,0] * at,3)
print(midu[15,0]),print(V[15,0]),print(T[15,0])
for i in range(0,30):
midut1[i,0] =round(-V[i,0](midu[i,0]-midu[i-1,0])/dx-midu[i,0](V[i,0]-V[i-1,0])/dx-V[i,0]midu[i,0](math.log(A[i])-math.log(A[i-1]))/ dx,3)
Vt1[i,0]= round(-V[i,0] * (V[i,0]-V[i-1,0])/dx-((T[i,0]-T[i-1,0])/dx+T[i,0](midu[i,0]-midu[i-1,0])/dx/midu[i,0])/bzb,1)
Tt1[i,0]= round(-V[i,0](T[i,0]-T[i-1,0])/dx-(bzb-1)T[i,0]((V[i,0]-V[i-1,0])/dx+(V[i,0]((math.log(A[i])-math.log(A[i-1]))/dx))),5)
print(midut[15,0]),print(Vt[15,0]),print(Tt[15,0])
print(midut[15,0]),print(Vt[15,0]),print(Tt[15,0])
print([i[0]for i in midut1])
print([i[0]for i in Vt1])
print([i[0]for i in Tt1])
s=len(Tt1)
print(s)
s=len(Tt1)
print(s)
for i in range(1, 30):
midut1[i,0]=round(0.5*(midut1[i,0]+midut[i,0]),3)
Vt1[i,0]=round(0.5*(Vt1[i,0]+Vt[i,0]),3)
Tt1[i,0]=round(0.5*(Tt1[i,0]+Tt[i,0]),4)
print(midut1[15,0]),print(Vt1[15,0]),print(Tt1[15,0])
print(midut[15,0]),print(Vt[15,0]),print(Tt[15,0])
s=len(Tt1)
print(s)
for i in range(1, 30):
midu[i,0]=round(midu[i,0]+midut1[i,0]*at,3)
V[i,0]=round(V[i,0]+Vt1[i,0]*at,3)
T[i,0]=round(T[i,0]+Tt1[i,0]*at,3)

V[0, 0] = round(2 * V[1, 0] - V[2, 0],3)
midu[29, 0] =round( 2 * midu[28, 0] - midu[27,0],3)
T[29, 0] = round(2 * T[28, 0] - T[27, 0],3)
midu[30, 0] =round( 2 * midu[29, 0] - midu[28,0],3)
T[30, 0] = round(2 * T[29, 0] - T[28, 0],3)
V[30, 0] = round(2 * V[29, 0] - V[28, 0],3)

print(midut[15,0]),print(Vt[15,0]),print(Tt[15,0])

for i in range(31):
y1[i,0]=round(midu[i,0]*T[i,0],4)
ma[i,0]=round(V[i,0]/(math.sqrt(T[i,0])),3)
print(“x方向”,[i[0]for i in x])
print(“面积”,[i[0]for i in A])
print(“密度”,[i[0]for i in midu])
print(“速度”,[i[0]for i in V])
print(“能量”,[i[0]for i in T])
print(“压力”,[i[0]for i in y1])
print(“马赫数”,[i[0]for i in ma])

用Python实现拟一维喷管流动的数值解相关推荐

持续不定期更新：CFDC++之拟一维喷管流动的数值解（1）
前言: 学习openfoam,或者准确来说,CFD,断断续续也有5个多月.3月底在本科学校写下的第一篇openfoam学习的博文,到现在在研一的学校里,尝试自己写代码计算cfd问题,慢慢一点一点地构建 ...
拟一维喷管流动的数值解——全亚声速等熵喷管流动的非守恒型CFD解法（MacCormack方法）
一.Matlab代码片 %全亚声速等熵喷管流动非守恒型麦考马克方法数值求解 clear; %清理内存变量 clc; %清理工作窗中的所有显示内容 r=1.4; %比热比 L=3; %喷管长度 i=3 ...
拟一维喷管流动的数值解——亚声速-超声速等熵喷管流动的非守恒型CFD解法（MacCormack方法）
一.Matlab代码片 %亚声速-超声速等熵喷管流动非守恒型麦考马克方法数值求解 clear; %清理内存变量 clc; %清理工作窗中的所有显示内容 r=1.4; %比热比 L=3; %喷管长度 ...
拟一维喷管流动的数值解——亚声速-超声速等熵喷管流动的守恒型CFD解法（MacCormack方法）
一.Matlab代码片 %亚声速-超声速等熵喷管流动守恒形CFD解法 MacCormack方法 clear; %清理内存变量 clc; %清理工作窗中的所有显示内容 r=1.4; %比热比 L=3; ...
持续不定期更新：CFDC++之拟一维喷管流动的数值解（2）
这篇博文将剩下的问题解决完.第一篇在: https://mp.csdn.net/postedit/101038218 在初始化步骤之后,就到了计算下一时间步的步骤了.计算之前先讲一讲这里用到的计算方法 ...
亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）
入门CFD,主要参考书目<计算流体力学基础及其应用>(John D.Anderson 著,吴颂平等译) 实现了第 7.3 节亚声速-超声速等熵喷管流动的CFD解法的代码,采用的是 ...
亚声速 – 超声速等熵喷管流动数值模拟(文字）
亚声速 – 超声速等熵喷管流动问题描述流过拉伐尔喷管的定常等熵流动在喷管,流动经过等熵地加速已达到超声速喷管出口处流体的压力.温度.速度和马赫数分别记为在喷管收缩段,流动是亚声速的:在喷管的 ...
喷管流动的守恒型CFD解法及激波捕捉（附完整代码）
入门CFD,主要参考书目<计算流体力学基础及其应用>(John D.Anderson 著,吴颂平等译) 实现了第 7.6 节激波捕捉的代码,采用的是 MacCormack 方法, ...
python中的一维卷积conv1d和二维卷积conv2d
先来看二维卷积conv2d conv2d(input, filter, strides, padding, use_cudnn_on_gpu=True, data_format="NHWC& ...

用Python实现拟一维喷管流动的数值解