Python计算非恒定速度v=10+x下密度函数的二阶Lax格式方法

核心原理修正 #

原有恒定速度场景的Lax格式基于平流方程∂ρ/∂t + c ∂ρ/∂x = 0的离散，当速度改为位置相关的v(x) = 10 + x时，平流方程变为守恒形式∂ρ/∂t + ∂(v(x)ρ)/∂x = 0，需要用通量形式重构Lax离散：
ρ_i^{n+1} = 0.5*(ρ_{i+1}^n + ρ_{i-1}^n) - (dt/(2dx)) * (F_{i+1}^n - F_{i-1}^n)
其中每个位置的通量为F_i = v(x_i) * ρ_i，x_i = xmin + i*dx。

同时要满足全局CFL稳定性条件：dt ≤ CFL * dx / max(v(x))，计算域x范围是0~100时最大速度为110，按这个计算dt可避免数值发散。

完整修改后代码 #

import numpy as np
import matplotlib.pyplot as plt

# 基础参数设置
MaxIter = 800  # 迭代次数
Nx = 1000  # 空间网格点数
xmin = 0.0
xmax = 100.0
CFL = 0.2  # CFL数，与原恒定速度版本参数保持一致
v_max = 10 + xmax  # 计算域内最大速度

dx = (xmax - xmin) / Nx
dt = CFL * dx / v_max  # 按全局最大速度计算时间步长保证稳定性

# 初始密度分布
oldL = np.concatenate((np.zeros(100), [1.0]*100, np.zeros(8*100)))
newL = np.zeros(Nx)
# 边界条件
newL[0] = 0.0
newL[-1] = 0.0

# 预计算每个网格点的速度v(x) = 10 + x
x_arr = np.linspace(xmin, xmax, Nx, endpoint=True)
v_arr = 10 + x_arr

j = 0
while j < MaxIter:
    # 先计算所有点的通量
    F_arr = v_arr * oldL
    for i in range(1, Nx-1):
        # 通量形式的Lax格式更新
        newL[i] = 0.5 * (oldL[i+1] + oldL[i-1]) - (dt/(2*dx)) * (F_arr[i+1] - F_arr[i-1])
    oldL = newL.copy()
    j += 1

# 结果输出和可视化
print("网格点数 =", Nx)
print("迭代次数 =", MaxIter)
print("空间步长dx =", dx)
print("时间步长dt =", dt)

xs = [dx * i for i in range(Nx)]
plt.xlabel("x(m)")
plt.ylabel("rho(x)")
plt.title("非恒定速度二阶Lax格式计算结果")
plt.scatter(xs, newL, color='b', label='Lax格式')
plt.legend()
plt.show()

关键修改说明 #

• 新增预计算所有网格点的位置相关速度数组，避免循环内重复计算
• 时间步长按全局最大速度调整，符合CFL稳定性要求
• 将原格式中常数速度的差分项替换为通量差分项，适配位置相关速度的守恒律求解

内容的提问来源于stack exchange，提问作者HelenK