import math
import numpy as np
from matplotlib import pyplot as plt
initial_condition_negative = [0.0, 1.0, 1.0]
initial_condition_positive = [0.0, 0.125, 0.1]
gamma = 1.4
t = 0.15
delta_t = 0.0001
delta_x = 0.005
def jacobian(u, c2):
line1 = np.array([0, 1, 0])
line2 = np.array([(gamma - 3) / 2 * u * u, (3 - gamma) * u, gamma - 1])
line3 = np.array([(gamma - 2) * u * u * u - u * c2 / (gamma - 1), c2 / (gamma - 1) + (3 - gamma) / 2 * u * u, gamma * u])
jacobian = np.vstack([line1, line2, line3])
return jacobian
def eigen(matrix):
eigenvalue, eigenvector = np.linalg.eig(matrix)
return eigenvalue, eigenvector
def Roe_solution():
u1, rho1, p1 = initial_condition_negative
u2, rho2, p2 = initial_condition_positive
x = np.linspace(-0.5, 0.5, 201)
u = []
rho = []
p = []
for e in x:
if e < 0.0:
u.append(u1)
rho.append(rho1)
p.append(p1)
else:
u.append(u2)
rho.append(rho2)
p.append(p2)
u = np.array(u)
rho = np.array(rho)
mom = rho * u
p = np.array(p)
E = p / (gamma - 1) + 0.5 * np.power(u, 2) * rho
H = (E + p) / rho
time_current = 0.0
while time_current < t:
f1 = rho * u
f2 = rho * np.power(u, 2) + p
f3 = u * (E + p)
f = np.vstack([f1, f2, f3])
u_left = u[:-1]
u_right = u[1:]
rho_left = rho[:-1]
rho_right = rho[1:]
mom_left = mom[:-1]
mom_right = mom[1:]
H_left = H[:-1]
H_right = H[1:]
p_left = p[:-1]
p_right = p[1:]
E_left = E[:-1]
E_right = E[1:]
f_left = f[:, :-1]
f_right = f[:, 1:]
rho_roe = np.power((np.sqrt(rho_left) + np.sqrt(rho_right)) / 2, 2)
u_roe = (np.sqrt(rho_left) * u_left + np.sqrt(rho_right) * u_right) / (np.sqrt(rho_left) + np.sqrt(rho_right))
H_roe = (np.sqrt(rho_left) * H_left + np.sqrt(rho_right) * H_right) / (np.sqrt(rho_left) + np.sqrt(rho_right))
p_roe = (gamma - 1) / gamma * (rho_roe * H_roe - 0.5 * rho_roe * np.power(u_roe, 2))
c2_roe = (gamma - 1) * (H_roe - np.power(u_roe, 2) * 0.5)
for i in range(len(x) - 2):
jacob = jacobian(u_roe[i], c2_roe[i])
eigenvalue, eigenvector = eigen(jacob)
eigenlambda = eigenvalue * np.identity(3)
jacob_abs = np.dot(np.dot(eigenvector, abs(eigenlambda)), np.linalg.inv(eigenvector))
delta_nU = np.array([[rho_right[i] - rho_left[i]], [mom_right[i] - mom_left[i]], [E_right[i] - E_left[i]]])
delta_pU = np.array([[rho_right[i + 1] - rho_left[i + 1]], [mom_right[i + 1] - mom_left[i + 1]], [E_right[i + 1] - E_left[i + 1]]])
Kn = 0.5 * np.dot(jacob_abs, delta_nU)
kp = 0.5 * np.dot(jacob_abs, delta_pU)
fn = 0.5 * (f_left[:, i].reshape(-1, 1) + f_right[:, i].reshape(-1, 1)) - Kn
fp = 0.5 * (f_left[:, i + 1].reshape(-1, 1) + f_right[:, i + 1].reshape(-1, 1)) - kp
rho[i + 1] = rho[i + 1] - delta_t / delta_x * (fp[0][0] - fn[0][0])
mom[i + 1] = mom[i + 1] - delta_t / delta_x * (fp[1][0] - fn[1][0])
E[i + 1] = E[i + 1] - delta_t / delta_x * (fp[2][0] - fn[2][0])
u[i + 1] = mom[i + 1] / rho[i + 1]
p[i + 1] = (gamma - 1) * (E[i + 1] - 0.5 * rho[i + 1] * math.pow(u[i + 1], 2))
H[i + 1] = (E[i + 1] + p[i + 1]) / rho[i + 1]
time_current += delta_t
result = {"position":x, "velocity":u, "density":rho, "pressure":p}
return result
def main():
dict_Roe = Roe_solution()
plt.figure(6)
plt.title("Roe Solution")
plt.plot(dict_Roe["position"], dict_Roe["velocity"], label = "Velocity and Position")
plt.plot(dict_Roe["position"], dict_Roe["density"], label = "Density and Position")
plt.plot(dict_Roe["position"], dict_Roe["pressure"], label = "Pressure and Position")
plt.legend()
plt.show()
if __name__ == '__main__':
main()
- 1
- 2
前往页