# -*- coding: utf-8 -*-
"""
Created on Tue Mar 2 20:21:39 2021
@author: Rolle
"""
import numpy as np
import matplotlib.pyplot as plt
import math
gama=1.4;c=0.9;t=0;flag=1;
xmax=1;nx=100;
U = np.zeros((3,103))#定义二维零数组
F =U.copy();W =U.copy();#同样定义F W为二维零数组
E = np.zeros(103) ;H = E.copy();U1 = E.copy();U2 = E.copy();U3 = E.copy()
roh = np.zeros(103);u = roh.copy(); p = roh.copy()
roh_L=np.zeros(102);u_L=roh_L.copy();E_L=roh_L.copy();H_L=roh_L.copy();p_L=roh_L.copy()
roh_R=np.zeros(102);u_R=roh_R.copy();E_R=roh_R.copy();H_R=roh_R.copy();p_R=roh_R.copy()
WL=np.zeros((102,3));WR=WL.copy();
FL=np.zeros((3,102));#FR=FL.copy();F=FL.copy()
ua=np.zeros(102);Ha=ua.copy();aa=ua.copy();alf1=ua.copy()
alf2=ua.copy();alf3=ua.copy();
K1=np.zeros((102,3));K2=K1.copy();K3=K1.copy();A=K1.copy()
lamd1=ua.copy();lamd2=ua.copy();lamd3=ua.copy()
dx=xmax/100;
rohL = 1;uL = 0.75;pL = 1;rohR = 0.125;uR = 0;pR =0.1;tmax = 0.2;xp = 30;#test1
#rohL=1;uL=-2;pL=0.4;rohR=1;uR=2;pR=0.4;tmax=0.012;xp=50;#test2
#=roh1;uL=0;pL=1000;rohR=1;uR=0;pR=0.01;tmax=0.012;xp=40;#test3
#rohL=5.99924;uL=19.5975;pL=460.894;rohR=5.99242;uR=-6.19633;pR=46.0950;tmax=0.035;xp=40;#test4
#rohL=1;uL=-19.59745;pL=1000;rohR=1;uR=-19.59745;pR=0.01;tmax=0.012;xp=80;#test5
#赋初值
for j in range(0,103):
if j<=xp:
roh[j] = rohL;u[j] = uL;p[j] = pL;#用方括号，不用圆括号,且二维数组要用两个方括号
else:
roh[j] = rohR;u[j] = uR;p[j] = pR;
#利用for循环得到E H的计算
#E[0][j] = 0.5*roh[0][j]*u[0][j]**2+p[0][j]/(gama-1);
#H[0][j] = (E[0][j]+p[0][j])/roh[0][j];
#U1[0][j]=roh[0][j];U2[0][j]=roh[0][j]*u[0][j];U3[0][j]=E[0][j];
#可直接利用数组进行计算
E = 0.5*roh*u**2+p/(gama-1);
H = (E+p)/roh;
U1 =roh.copy();U2=roh*u;U3=E.copy();
U=np.c_[U1,U2,U3] #将三个行向量纵向合并
while(flag):
roh_L=roh[0:102];u_L=u[0:102];E_L=E[0:102];H_L=H[0:102];p_L=p[0:102]
roh_R=roh[1:103];u_R=u[1:103];E_R=E[1:103];H_R=H[1:103];p_R=p[1:103]
WL=np.c_[roh_L,u_L,p_L];WR=np.c_[roh_R,u_R,p_R]
dt=c*dx/np.nanmax((np.abs(u)+np.sqrt(gama*p/roh)))#np.nanmax()用于返回非空值的最大值
t=t+dt
if t>=tmax:
flag=0
#roe平均
for j in range(0,102):
ua[j]=(math.sqrt(roh_L[j])*u_L[j]+math.sqrt(roh_R[j])*u_R[j])/(math.sqrt(roh_L[j])+math.sqrt(roh_R[j]));
Ha[j]=(math.sqrt(roh_L[j])*H_L[j]+math.sqrt(roh_R[j])*H_R[j])/(math.sqrt(roh_L[j])+math.sqrt(roh_R[j]));
aa[j]=math.sqrt((gama-1)*(Ha[j]-0.5*ua[j]**2));
alf2[j]=(gama-1)/aa[j]**2*((U1[j+1]-U1[j])*(Ha[j]-ua[j]**2)+ua[j]*(U2[j+1]-U2[j])-(U3[j+1]-U3[j]));
alf1[j]=0.5/aa[j]*((U1[j+1]-U1[j])*(ua[j]+aa[j])-(U2[j+1]-U2[j])-aa[j]*alf2[j]);
alf3[j]=(U1[j+1]-U1[j])-alf1[j]-alf2[j];
K1=np.c_[np.ones(102),ua-aa,Ha-ua*aa];
K2=np.c_[np.ones(102),ua,0.5*ua**2];
K3=np.c_[np.ones(102),ua+aa,Ha+ua*aa];
lamd1=ua-aa;
lamd2=ua;
lamd3=ua+aa;
for j in range(0,102):
#熵修正
eps=0.4
if abs(lamd1[j])<eps:
lamd1[j]=(lamd1[j]**2+(2*eps)**2)/(2*eps)
A=0.5*(alf1*np.abs(lamd1)*K1.T+alf2*np.abs(lamd2)*K2.T+alf3*np.abs(lamd3)*K3.T)
WL=WL.T;WR=WR.T
FL=np.c_[WL[0]*WL[1],WL[0]*WL[1]**2+WL[2],(WL[2]+E_L)*WL[1]]
FR=np.c_[WR[0]*WR[1],WR[0]*WR[1]**2+WR[2],(WR[2]+E_R)*WR[1]]
F=0.5*( FL+ FR )-A.T
for j in range(1,102):
U[j]=U[j]-dt/dx*(F[j]-F[j-1])
U[0]=U[1];U[102]=U[101]#边界条件
roh=U[:,0].copy()
u=U[:,1]/U[:,0]
p=(gama-1)*(U[:,2]-0.5/U[:,0]*U[:,1]**2)
E=U[:,2].copy()
H=(E+p)/U[:,0]
e=p/(gama-1)/roh
U1=U[:,0].copy();U2=U[:,1].copy();U3=U[:,2].copy()
x=np.arange(0,1.03,0.01)
plt.figure(dpi=30)
plt.figure(figsize=(8,15))
plt.subplot(421)
plt.plot(x,roh,'bo',ms=1)
plt.xlabel('x')
plt.ylabel('Denisty')
plt.subplot(422)
plt.plot(x,u,'b')
plt.xlabel('x')
plt.ylabel('Velocity')
plt.subplot(423)
plt.plot(x,p,'b')
plt.xlabel('x')
plt.ylabel('Pressure')
plt.subplot(424)
plt.plot(x,e,'b')
plt.xlabel('x')
plt.ylabel('Internal energy')
plt.savefig('激波管图片.jpg')
评论10