（求问）python写的用SCMP模拟液滴验证拉普拉斯定律，有人能帮忙看看代码哪里错了吗
本帖最后由 嘿热干面 于 2018-7-25 17:29 编辑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]
e_i = [ 1 * c,0]
e_i = [ 0,1 * c]
e_i = [-1 * c,0]
e_i = [ 0, -1 * c]
e_i = [ 1 * c,1 * c]
e_i = [-1 * c,1 * c]
e_i = [-1 * c, -1 * c]
e_i = [ 1 * c, -1 * c]
w_i = 1.0/36.0 * np.ones(q) # Lattice weights.
w_i)] = 1.0 / 9.0
w_i = 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**2 + e_i**2 == 1:
omega_i = 1./9.
elif e_i**2 + e_i**2 == 2:
omega_i = 1./36.
def equilibrium(rho, u):
u_squ = u**2 + u**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 = 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 = np.sqrt(2 * (p - c_squ * rho) / (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
yp = y + e_i
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 += omega_i * psx * e_i
Ff_temp += omega_i * psx * e_i
Fint = - G * psx * Ff_temp
Fint = - G * psx * Ff_temp
F = Fint
return F
# get the quilibrium velocity in the body force term
def get_ueq(rho, u, F):
ueq = u + F * dt /rho
ueq = u + F * 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 = np.roll(np.roll(f_out, e_i, axis=0), e_i, 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)
页:
[1]