一维激波管问题——计算流体力学（Lax-Friedrichs格式）

import numpy as np import matplotlib.pyplot as plt from matplotlib.ticker import MultipleLocator, FormatStrFormatter import funclaxfriedrichs as fc # Set parameter N = 1000 # Grid number cfl = 0.5 # Courant number xmax = 1 xmin = -1 x = np.linspace(xmin, xmax, N+1) # Position coordinates dx = (xmax-xmin)/N n_dis = int(0.5*N) t_max = 0.14 # Time gamma = 1.4 param0 = [0.445, 0.311, 8.928, 0.5, 0.0, 1.4275] # Density Mass Flux Energy flag = 1 steps = 0 maxsteps = 1e5 current_time = 0 #set initial value u_vector = fc.euler_generate(N, n_dis, *param0) u_initial = u_vector u_t = np.zeros((3, N+1)) # Calculate analytical solution (Riemann sulution) initial_3 = [0.445, 0.311, 8.928] initial_2 = [0.345, 0.527, 6.570] initial_1 = [1.304, 1.994, 7.691] initial_0 = [0.5, 0., 1.4275] speed3 = -2.633 speed2 = -1.636 speed1 = 1.529 speed0 = 2.480 dis1 = int((speed3*t_max+1.)*N/2.) dis2 = int((speed2*t_max+1.)*N/2.) dis3 = int((speed1*t_max+1.)*N/2.) dis4 = int((speed0*t_max+1.)*N/2.) rho = np.zeros(N+1) mass = np.zeros(N+1) energy = np.zeros(N+1) # Density rho[:dis1+1] = initial_3[0] rho[dis2+1:dis3+1] = initial_2[0] rho[dis3+1:dis4+1] = initial_1[0] rho[dis4+1:] = initial_0[0] rho1 = np.delete(rho,np.linspace(dis1+1,dis2,dis2-dis1).astype('int64')) # Mass Flux mass[:dis1+1] = initial_3[1] mass[dis2+1:dis3+1] = initial_2[1] mass[dis3+1:dis4+1] = initial_1[1] mass[dis4+1:] = initial_0[1] mass1 = np.delete(mass,np.linspace(dis1+1,dis2,dis2-dis1).astype('int64')) # Energy energy[:dis1+1] = initial_3[2] energy[dis2+1:dis3+1] = initial_2[2] energy[dis3+1:dis4+1] = initial_1[2] energy[dis4+1:] = initial_0[2] energy1 = np.delete(energy,np.linspace(dis1+1,dis2,dis2-dis1).astype('int64')) solution = np.array([rho1, mass1, energy1]) x1 = x[:dis1+1] x2 = x[dis2+1:dis3+1] x3 = x[dis3+1:dis4+1] x4 = x[dis4+1:] x_solution = np.hstack((x1, x2, x3, x4)) # Calculate numerical solution dF = fc.fluxDiscre(dx, u_vector) while flag: steps = steps+1 u_0 = fc.consevation_to_physics(u_vector) c = np.sqrt(gamma*u_0[2, :]/u_0[0, :]) # Speed of sound dt = cfl*dx/max(abs(u_0[1, :])+c) # Step of time if steps > maxsteps: print('warning: maxsteps reach') break if current_time+dt >= t_max: dt = t_max-current_time flag = 0 current_time = current_time+dt u_t[:, 1:-1] = 0.5*(u_vector[:, :-2]+u_vector[:, 2:])-dt*dF # 推进 u_t[:, 0] = u_t[:, 1] # Boundary u_t[:, -1] = u_t[:, -2] dF = fc.fluxDiscre(dx, u_t) # Calculate dF u_vector = u_t.copy() # Plot image fig = plt.figure(figsize=(14,8)) plt.rcParams['xtick.direction'] = 'in' plt.rcParams['ytick.direction'] = 'in' ax = plt.subplot(311) plt.plot(x, u_initial[0, :], color='k', ls=':') plt.plot(x, u_t[0, :], 'k', ls='--') plt.plot(x_solution, solution[0, :], 'k') plt.scatter(x, u_t[0, :], color='', marker='o', edgecolors='k', s=80) plt.xlim((-1.0,1.0)) plt.ylim((0.2,1.4)) plt.ylabel('Density') plt.xticks(np.arange(-1.0,1.5,0.5)) xminorLocator = MultipleLocator(0.1) yminorLocator = MultipleLocator(0.05) ax.xaxis.set_minor_locator(xminorLocator) ax.yaxis.set_minor_locator(yminorLocator) plt.title('Time = '+str(t_max)+' /Lax-Friedrichs/') ax = plt.subplot(312) plt.plot(x, u_initial[1, :], color='k', ls=':') plt.plot(x, u_t[1, :], 'k', ls='--') plt.plot(x_solution, solution[1, :], 'k') plt.scatter(x, u_t[1, :], color='', marker='o', edgecolors='k', s=80) plt.xlim((-1.0,1.0)) plt.ylim((-0.5,2.5)) plt.ylabel('Mass Flux') plt.xticks(np.arange(-1.0,1.5,0.5)) xminorLocator = MultipleLocator(0.1) yminorLocator = MultipleLocator(0.1) ax.xaxis.set_minor_locator(xminorLocator) ax.yaxis.set_minor_locator(yminorLocator) plt.title('Time = '+str(t_max)+' /Lax-Friedrichs/') ax = plt.subplot(313) plt.plot(x, u_initial[2, :], color='k', ls=':') plt.plot(x, u_t[2, :], 'k', ls='--') plt.plot(x_solution, solution[2, :], 'k') plt.scatter(x, u_t[2, :], color='', marker='o', edgecolors='k', s=80) plt.xlim((-1.0,1.0)) plt.ylim((0.0,10)) plt.ylabel('Energy') plt.xticks(np.arange(-1.0,1.5,0.5)) xminorLocator = MultipleLocator(0.1) yminorLocator = MultipleLocator(0.5) ax.xaxis.set_minor_locator(xminorLocator) ax.yaxis.set_minor_locator(yminorLocator) plt.title('Time = '+str(t_max)+' /Lax-Friedrichs/') fig.tight_layout() plt.savefig('./picture/LaxFriedrichs/laxfriedrichs'+str(N)+'.eps', format='eps') plt.show()