（求问）python写的用SCMP模拟液滴验证拉普拉斯定律，有人能帮忙看看代码哪里错了吗

嘿热干面

python写的代码，用SCMP模型模拟液滴处在蒸汽中，只考虑了液体-液体之间的作用力，四周都采用周期性边界，但是算了几步之后，密度就变得超级大，有人能帮我看下代码哪里可能有问题吗


import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

# Parameters definition
nx = 300
ny = 200
maxIter = 10000
dx = 1.0
dt = 1.0
RR = 30.0
c = dx/dt  # lattice velocity
c_squ = c**2/3.
G = -1.   # fluid interaction strength
tau = 1.0
nu = c_squ * dt * (tau - 0.5)

# parameters in the P-R equation of state
a = 3./49.
b = 2./21.
R = 1.0
omega = 0.344    # acentric factor in P-R EOS, 0.344 for water
Tc = 0.0778 * a / (0.45724 * b * R)    # critical temperature
Ts = 0.86 * Tc    # saturated temperature

rho_l = 6.50    # density of liquid
rho_g = 0.38    # density of vapor

# D2Q9 Lattice constants
# e_i gives all 9 velocity vectors.
q = 9      #D2Q9
e_i = np.zeros((q, 2), dtype=np.int64)
e_i[0] = [ 0,  0]
e_i[1] = [ 1 * c,  0]
e_i[2] = [ 0,  1 * c]
e_i[3] = [-1 * c,  0]
e_i[4] = [ 0, -1 * c]
e_i[5] = [ 1 * c,  1 * c]
e_i[6] = [-1 * c,  1 * c]
e_i[7] = [-1 * c, -1 * c]
e_i[8] = [ 1 * c, -1 * c]
w_i = 1.0/36.0 * np.ones(q)       # Lattice weights.
w_i[np.asarray([np.linalg.norm(e)<1.1 for e in e_i])] = 1.0 / 9.0
w_i[0] = 4.0 / 9.0

f_eq = np.zeros((q, nx, ny))
ueq = np.zeros((2, nx, ny))
p = np.zeros((nx, ny))
psx = np.zeros((nx, ny))
Fint = np.zeros((2, nx, ny))    # interaction force between fluid nodes
omega_i = np.zeros((q,1))
Fg = np.zeros((2, nx, ny))
for i in range(q):
    if e_i[i, 0]**2 + e_i[i, 1]**2 == 1:
        omega_i = 1./9.
    elif e_i[i, 0]**2 + e_i[i, 1]**2 == 2:
        omega_i = 1./36.

def equilibrium(rho, u):
    u_squ = u[0]**2 + u[1]**2
    cu = np.dot(e_i, u.transpose(1,0,2))
    for i in range(q):
        f_eq = (rho * w_i * (1. + cu/c_squ + 0.5 * cu**2/c_squ**2  - 0.5 * u_squ/c_squ))
    return f_eq

def initialize():
    ''' the domains are filled with vapor and a spherical liquid droplet is
    placed at the center'''
    rho = rho_g * np.ones((nx, ny))
    for y in range(ny):
        for x in range(nx):
            if ((x - nx/2)**2 + (y - ny/2)**2 <= RR**2):
                rho[x, y] = rho_l
    u = np.zeros((2, nx, ny))
    f_eq = equilibrium(rho, u)
    f = f_eq.copy()
    return f

# obtain the macro variables such as velocity and density
def get_var(f):
    rho = np.sum(f, axis = 0)
    u = np.dot(e_i.transpose(), f.transpose((1,0,2)))/rho
    return rho, u

# get pressure and psx using the P-R equation of state
def get_p(rho):
    epl = (1 + (0.37464 + 1.54226 * omega - 0.26992 * omega**2) * (1 - (Ts/Tc)**0.5))**2
    # The temperature in the domain is set as the saturated temperature
    p = rho * R * Ts/(1 - b * rho) - a * epl * rho**2/(1 + 2 * b * rho - (b * rho)**2)
    C0 = 6.0
    for y in range(ny):
        for x in range(nx):
            psx[x, y] = np.sqrt(2 * (p[x, y] - c_squ * rho[x, y]) / (C0 * G))
    return p, psx

def calcu_Force(rho, psx):
    Ff_temp = np.zeros((2, nx, ny))
    for y in range(ny):
        for x in range(nx):
            for i in range(q):
                xp = x + e_i[i, 0]
                yp = y + e_i[i, 1]
                if (xp < 0):   xp = nx-1
                if (xp > nx-1):   xp = 0
                if (yp < 0):   yp = ny-1   # !!!!
                if (yp > ny-1):   yp = 0   # !!!!
               
                # terms behind the summation symbol
                Ff_temp[0, x, y] += omega_i * psx[xp, yp] * e_i[i, 0]
                Ff_temp[1, x, y] += omega_i * psx[xp, yp] * e_i[i, 1]
    Fint[0, :, :] = - G * psx * Ff_temp[0, :, :]
    Fint[1, :, :] = - G * psx * Ff_temp[1, :, :]
    F = Fint
    return F

# get the quilibrium velocity in the body force term
def get_ueq(rho, u, F):
    ueq[0, :, :] = u[0, :, :] + F[0, :, :] * dt /rho
    ueq[1, :, :] = u[1, :, :] + F[1, :, :] * dt /rho
    return ueq

def collision(rho, u, ueq, f):  
    df = equilibrium(rho, ueq) - equilibrium(rho, u)    #force term in the LBE, calculated using the exact difference method
    f_eq = equilibrium(rho, u)
    for i in range(q):
        f_out = f - (f - f_eq) / tau + df
    return f_out

def stream(f_out):
    for i in range(q):
        f[i, :, :] = np.roll(np.roll(f_out[i, :, :], e_i[i, 0], axis=0), e_i[i, 1], axis=1)
    return f
     
# initialization
f = initialize()

print('Begin iteration')
# main loop
for t in range(maxIter):
    print('Iteration loop:' , t)
    rho, u = get_var(f)
    p, psx = get_p(rho)
    F = calcu_Force(rho, psx)
    ueq = get_ueq(rho, u, F)
    f_out = collision(rho, u, ueq, f)
    f = stream(f_out)