一维理想气体Euler方程 有限体积求解 附python代码

一维理想气体Euler方程 有限体积求解 附python代码

Euler方程记为
$$U_t+F(U{)}_x=0$$
经过有限体积离散之后，守恒量更新格式为，
$$U_i^{n+1}=U_i-\frac{\Delta t}{\Delta x}(F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}})$$
其中，$F_{i-\frac{1}{2}}$是数值通量，这里使用local Lax Friedriches数值通量。

数值算例

在$[0,1]$的区间内，网格数量为1000，CFL数为0.5.

下图是sod激波管问题在 $t=0.2$s的密度，速度和压强的分布结果。

代码

import matplotlib.pyplot as plt
import numpy as np


class Fluid_data:
    DIM = 1
    U = 1e-8 * np.ones(DIM + 2, dtype=float)
    gamma = 1.4
    P = 0.0
    rho = 1e-8
    u = 0

    def __init__(self, DIM=1, U=1e-8 * np.ones(DIM + 2, dtype=float), gamma=1.4):
        self.DIM = DIM
        self.U = U
        self.gamma = gamma
        self.update()

    def set_U(self, U):
        self.U = U

    def set_gamma(self, gamma):
        self.gamma = gamma

    def update(self):
        self.rho = self.U[0]
        self.u = self.U[1] / self.U[0]
        self.P = (self.gamma - 1) * (self.U[-1] - 0.5 * np.power(np.linalg.norm(self.U[1]), 2) / self.U[0])

    def is_fisable(self):
        if (self.P >= 0 and self.U[0] >= 0):
            return True
        else:
            return False

    def sound_speed(self):
        return np.sqrt(self.gamma * self.P / self.U[0])

    def max_speed(self):
        self.update()
        return abs(self.U[1] / self.U[0]) + self.sound_speed()

    # python 中没有 Switch
    def eig_value(self, i):
        self.update()
        if i == 1:
            return self.u - self.sound_speed()
        elif i == 2:
            return self.u
        elif i == 3:
            return self.u + self.sound_speed()

    def flux(self):
        self.update()
        Fu = np.zeros(DIM + 2)
        Fu[0] = self.U[1]
        Fu[1] = self.U[1] * self.U[1] / self.U[0] + self.P
        Fu[2] = (self.U[2] + self.P) * self.U[1] / self.U[0]
        return Fu


# 初值函数
def init_fluid_data(U_data, Ele_c):
    for i in range(len(U_data)):
        if Ele_c[i] < 0.5:
            rho = 1
            u = 0
            P = 1.0
            U_temp = np.zeros(U_data[i].DIM + 2)
            U_temp[0] = rho
            U_temp[1] = rho * u
            U_temp[2] = P / (U_data[i].gamma - 1) + 0.5 * rho * u * u
            U_data[i].set_U(U_temp)
            U_data[i].update()
        else:
            rho = 0.125
            u = 0
            P = 0.1
            U_temp = np.zeros(U_data[i].DIM + 2)
            U_temp[0] = rho
            U_temp[1] = rho * u
            U_temp[2] = P / (U_data[i].gamma - 1) + 0.5 * rho * u * u
            U_data[i].set_U(U_temp)
            U_data[i].update()


def numerical_flux(Ul, Ur, label="LLF"):
    if label == "LLF":
        max_speed = max(Ul.max_speed(), Ur.max_speed())
        return 0.5 * (Ul.flux() + Ur.flux()) - 0.5 * max_speed * (Ur.U - Ul.U)

DIM = 1
CFL = 0.5
U_init = 1e-8 * np.ones(DIM + 2, dtype=float)
gamma = 1.4
x_min = 0
x_max = 1
Nx = 1000
max_t = 0.2
point = np.linspace(0, 1, Nx + 1)
n_edge = np.shape(point)[0]
n_ele = n_edge - 1
ele_vol = (x_max - x_min) / n_ele
edge_mark = np.zeros(n_edge, dtype=int)
edge_mark[0] = 1
edge_mark[-1] = 1

ele2edge = np.array([np.arange(0, n_edge - 1), np.arange(1, n_edge)])
ele2edge = ele2edge.transpose()
ele2point = ele2edge
# edge2cell = np.array([n_edge,2])
edge2cell = np.array([np.arange(-1, n_ele), np.arange(0, n_ele + 1)])
edge2cell[-1, -1] = -1
edge2cell = edge2cell.transpose()
# 将向量组合在一起 可能不再是array了
ele_centroid = 0.5 * (point[0:-1] + point[1:])

# 对数据进行初始化
U_data = []
for i in range(n_ele):
    U_data.append(Fluid_data())

# 守恒量 初始值
init_fluid_data(U_data, ele_centroid)

t = 0
while (t < max_t):
    # 计算时间步长
    max_speed = 0
    for i in range(n_ele):
        max_speed = max(max_speed, U_data[i].max_speed())

    dt = CFL * ele_vol / max_speed

    # 计算数值通量
    Flux_vec = np.zeros([DIM + 2, n_edge])
    # lax-Friedriches
    for i in range(n_edge):
        if edge2cell[i, 0] == -1:
            Ur = U_data[edge2cell[i, 1]]
            Ul = Ur
        elif edge2cell[i, 1] == -1:
            Ul = U_data[edge2cell[i, 0]]
            Ur = Ul
        else:
            Ul = U_data[edge2cell[i, 0]]
            Ur = U_data[edge2cell[i, 1]]

        Flux_vec[:, i] = numerical_flux(Ul, Ur)
    # 更新守恒量
    for i in range(n_ele):
        U_data[i].U = U_data[i].U - (dt / ele_vol) * (Flux_vec[:, ele2edge[i, 1]] - Flux_vec[:, ele2edge[i, 0]])
        U_data[i].update()

    t = t + dt


Rho = np.zeros(n_ele)
V = np.zeros(n_ele)
P = np.zeros(n_ele)

for i in range(n_ele):
    Rho[i] = U_data[i].U[0]
    V[i] = U_data[i].U[1] / U_data[i].U[0]
    P[i] = U_data[i].P
plt.figure()
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.subplot(2, 2, 1)
plt.plot(ele_centroid, Rho)
plt.xlabel('x')
plt.ylabel(r'$\rho$')
plt.subplot(2, 2, 2)
plt.plot(ele_centroid, V)
plt.xlabel('x')
plt.ylabel(r'$u$')
plt.subplot(2, 2, 3)
plt.plot(ele_centroid, P)
plt.xlabel('x')
plt.ylabel(r'$p$')
plt.show()