通过使用 Runge Kutta O(4) 方法将 2D Euler 方程的直接积分写在一个宽而短的盒子中，在该盒子中流体通过水平面进入和离开，我发现 CFL 系数，这里近似计算如下，振荡随着时间的推移，幅度增加。如果积分没有停止，速度场就会发散。

我是否面临通常被称为直接求解器对于刚性方程不稳定的情况？

$CFL_{coefficient}=\frac{max(u,v)dt}{h}$

在哪里$u,v$分别是 x 和 y 中的速度

$h$空间特征长度

$dt$整合步骤。

密度被视为常数。

一些笔记。我尝试在边界处强制执行零切向速度，作为模拟粘度的一种简单近似，而没有明确添加粘性项，它似乎会促进不稳定性。不清楚。

主要注意事项是考虑到质量守恒，如下所示，根本没有被强制执行。

$\frac{\partial u }{\partial x} + \frac{\partial v }{\partial y} = 0$

Python3 中的最小工作示例如下。请不要在已经存在名为“results”的文件夹的任何空间中运行它。要求是 numpy 和 matplotlib。

编辑：的确，答案是有一个错误。我按照接受的答案的建议简化了求解器以使用 O(1) Euler 方法。我已经更新了 MWE。现在计算时间步长以满足每个积分步长的 CLF 系数。

import numpy as np
import matplotlib.pyplot as plt
import os,sys

try:
    os.mkdir('results')
except:
    pass
try:
    os.system('rm results/*')
except:
    pass

Lx,Ly = (30,30)




RHO = 12
DT = 0.02 # Main integration timestep
H = 1
RHO_H = RHO / H
BOUNDARYVAL = 0.1


Ux = np.ones((Lx, Ly)) * BOUNDARYVAL / 2
Uy = np.zeros((Lx, Ly))

Ux[0,:] = BOUNDARYVAL
#Ux[1,:] = BOUNDARYVAL/2
#Ux[2,:] = BOUNDARYVAL/4


# Enforce the Courant-Friedrichs-Lewy condition
CFL = BOUNDARYVAL * DT / H
assert(CFL <= 0.5)


LAPS = 300
FREQ = 15
D = len(str(LAPS))
for _ in range(LAPS*FREQ):


    if _>0: 
        DT = CFL * H / max([max(Uy.flatten()), max(Ux.flatten())])
    

    # Down to up velocity vertically
    Ux[0,:] = BOUNDARYVAL
    #Ux[-1,:] = BOUNDARYVAL

    # Vertical velocity in the borders can be whatever because there is no viscosity
    #Ux[:,-1] = 0
    #Ux[:,0] = 0

    # Horizonal velocity in the horizontal borders can be whatever
    #Uy[0,:] = 0
    #Uy[-1,:] = 0

    # Horizontal velocity in the vertical borders must be zero as "fluid cant escape through the sides"
    Uy[:,0] = 0
    Uy[:,-1] = 0

    # Plot
    if _%FREQ==0:
        print(f"DT is {DT}")
        print(f"Lap {_+1} of {FREQ*LAPS}!")
        plt.imshow(Ux * Ux + Uy * Uy, cmap = 'coolwarm', vmin=min([BOUNDARYVAL,0]), vmax=max([BOUNDARYVAL,0]))
        plt.colorbar()
        plt.savefig(f'results/{str(_//FREQ).zfill(D)}.png')
        plt.close()
    
    # Compute next state
    # Mass conservation eq. is divergence(U)=0
    # Momentum conservation eqs. are dU/dt + (U \cdot \grad) U = 0
    # Where there is no gravity, viscosity, nor external pressure gradient being enforced
    
    # dU/dt = -RHO * (Ux d/dx + Uy d/dy)U
    #
    # so we compute auxiliaries 
    #
    dUx_dx = (Ux[2:,1:-1] - Ux[:-2,1:-1]) / 2 / H
    dUx_dy = (Ux[1:-1,2:] - Ux[1:-1,:-2]) / 2 / H
    dUy_dx = (Uy[2:,1:-1] - Uy[:-2,1:-1]) / 2 / H
    dUy_dy = (Uy[1:-1,2:] - Uy[1:-1,:-2]) / 2 / H
    # 
    # and evolve the field
    #
    #print(dUx_dx.shape, Ux.shape)
    #plt.imshow(dUx_dx);plt.show()
    #MX = np.max((DT * RHO * (Ux[1:-1,1:-1] * dUx_dx + Uy[1:-1,1:-1] * dUx_dy)).flatten())
    #print(f'max incoming update is {MX}')
    Ux[1:-1,1:-1] -= DT * RHO * (Ux[1:-1,1:-1] * dUx_dx + Uy[1:-1,1:-1] * dUx_dy)
    Uy[1:-1,1:-1] -= DT * RHO * (Ux[1:-1,1:-1] * dUy_dx + Uy[1:-1,1:-1] * dUy_dy)



try:
    os.remove('result.avi')
except: pass

print(f'To make a GIF showing the evolution please run ffmpeg  -framerate 5 -i "results/%0{D}d.png" -vcodec mpeg4 result.avi')