Hi, all,
I am trying to implement constant pressure boundary condition in Zou-He method for poiseuille flow. However, the velocity profile or the velocity’s norm profile never matched the velocity profile of poiseuille flow in a tube. After checking the figures on results, it was found that the velocity near the outlet of the tube is always higher (). I have already checked the code many times and had no idea on which part is wrong. Is there anyone who can help me check what’s wrong with the code and algorithm below? The code is written in Pytyon3. Thanks so much!!
[code=“python”]
“”"
File: constantpressure.py
Purpose: Test the constant pressure boundary condition algorithm in single phase
flow
“”"
import sys, os
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
nx = 32
ny = 16
tau = 1.2
nIteration = 5000
#define the wall
isWall = sp.empty([ny, nx], dtype = ‘bool’)
isWall[:, :] = False
isWall[0, :] = True
isWall[-1, :] = True
oppDirection = [0, 3, 4, 1, 2, 7, 8, 5, 6]
#define the constant for the weight coefficients
weightCoeff = sp.empty(9)
weightCoeff[0] = 4./9.; weightCoeff[1:5] = 1./9.; weightCoeff[5:] = 1./36.
#define the velocity in the different directions of lattice
microVelocity = sp.empty([9, 2], dtype = ‘int64’)
microVelocity[0] = np.array([0, 0]); microVelocity[1] = np.array([1, 0])
microVelocity[2] = np.array([0, 1]); microVelocity[3] = np.array([-1, 0])
microVelocity[4] = np.array([0, -1]); microVelocity[5] = np.array([1, 1])
microVelocity[6] = np.array([-1, 1]); microVelocity[7] = np.array([-1, -1])
microVelocity[8] = np.array([1, -1])
#define the coefficients in f_eq
eqCoeff1 = 3.0; eqCoeff2 = 9./2.; eqCoeff3 = -3./2.
#initialize distribution function on each lattice node
particleDisFunc = sp.ones([9, ny, nx])
for i in sp.arange(ny):
for j in sp.arange(nx):
particleDisFunc[:, i, j] *= weightCoeff
timeStep = 0
#define macroscopic parameters
rho = sp.zeros([ny, nx])
velocityX = sp.zeros([ny, nx])
velocityY = sp.zeros([ny, nx])
#initialize equilibrium function
equilibriumFunc = sp.empty([9, ny, nx])
specificDensityInlet = 1.1
specificDensityOutlet = 0.8
while (timeStep < nIteration):
print(“This is the %g step” % timeStep)
timeStep += 1
if (timeStep % 10 == 0 and timeStep != 0):
X, Y = np.meshgrid(np.arange(1,nx), np.arange(1,ny))
plot2 = plt.figure()
vnorm = sp.sqrt(sp.power(velocityX, 2.) + sp.power(velocityY, 2.))
plt.imshow(vnorm[1:-1, 1:-1].real, origin = 'lower')
plt.colorbar()
plt.savefig('norm%g.png' % timeStep)
#streaming
for i in sp.arange(9):
particleDisFunc[i, :, :] = np.roll(np.roll(particleDisFunc[i, :,:], \
microVelocity[i, 1], axis = 0), microVelocity[i, 0], \
axis = 1)
#Implement constant pressure boundary condition
##Inlet
velocityX[1:-1, 0] = 1. - 1./specificDensityInlet * (particleDisFunc[0, 1:-1, 0] + \
particleDisFunc[2, 1:-1, 0] + particleDisFunc[4, 1:-1, 0] \
+ 2. * (particleDisFunc[3, 1:-1, 0] + particleDisFunc[6, 1:-1, 0] + \
particleDisFunc[7, 1:-1, 0]))
particleDisFunc[1, 1:-1, 0] = particleDisFunc[3, 1:-1, 0] + 2./3. * specificDensityInlet * \
velocityX[1:-1, 0]
particleDisFunc[5, 1:-1, 0] = particleDisFunc[7, 1:-1, 0] - 1./ 2. * (particleDisFunc[2, 1:-1, 0] - \
particleDisFunc[4, 1:-1, 0]) + 1./6. * specificDensityInlet * \
velocityX[1:-1, 0]
particleDisFunc[8, 1:-1, 0] = particleDisFunc[6, 1:-1, 0] + 1./2. * (particleDisFunc[2, 1:-1, 0] - \
particleDisFunc[4, 1:-1, 0]) + 1. / 6. * specificDensityInlet * \
velocityX[1:-1, 0]
##Outlet
velocityX[1:-1, -1] = -1. + 1. / specificDensityOutlet * (particleDisFunc[0, 1:-1, -1] + \
particleDisFunc[2, 1:-1, -1] + particleDisFunc[4, 1:-1, -1] + 2. * \
(particleDisFunc[1, 1:-1, -1] + particleDisFunc[5, 1:-1, -1] + \
particleDisFunc[8, 1:-1, -1]))
particleDisFunc[3, 1:-1, -1] = particleDisFunc[1, 1:-1, -1] - 2./3. * specificDensityOutlet * \
velocityX[1:-1, -1]
particleDisFunc[6, 1:-1, -1] = particleDisFunc[8, 1:-1, -1] + 1./2. * (particleDisFunc[4, 1:-1, -1] - \
particleDisFunc[2, 1:-1, -1]) - 1. / 6. * specificDensityOutlet * \
velocityX[1:-1, -1]
particleDisFunc[7, 1:-1, -1] = particleDisFunc[5, 1:-1, -1] + 1./2. * (particleDisFunc[2, 1:-1, -1] - \
particleDisFunc[4, 1:-1, -1]) - 1./6. * specificDensityOutlet * \
velocityX[1:-1, -1]
##node on corner
##Inlet bottom
velocityX[0, 0] = 0.
velocityY[0, 0] = 0.
particleDisFunc[1, 0, 0] = particleDisFunc[3, 0, 0]
particleDisFunc[2, 0, 0] = particleDisFunc[4, 0, 0]
particleDisFunc[5, 0, 0] = particleDisFunc[7, 0, 0]
particleDisFunc[6, 0, 0] = 1./ 2. * (specificDensityInlet - (particleDisFunc[0, 0, 0] + \
particleDisFunc[1, 0, 0] + particleDisFunc[2, 0, 0] + particleDisFunc[3, 0, 0] + \
particleDisFunc[4, 0, 0] + particleDisFunc[5, 0, 0] + \
particleDisFunc[7, 0, 0]))
particleDisFunc[8, 0, 0] = 1./ 2. * (specificDensityInlet - (particleDisFunc[0, 0, 0] + \
particleDisFunc[1, 0, 0] + particleDisFunc[2, 0, 0] + particleDisFunc[3, 0, 0] + \
particleDisFunc[4, 0, 0] + particleDisFunc[5, 0, 0] + \
particleDisFunc[7, 0, 0]))
##inlet top
velocityX[-1, 0] = 0.
velocityY[-1, 0] = 0.
particleDisFunc[1, -1, 0] = particleDisFunc[3, -1, 0]
particleDisFunc[4, -1, 0] = particleDisFunc[2, -1, 0]
particleDisFunc[8, -1, 0] = particleDisFunc[6, -1, 0]
particleDisFunc[5, -1, 0] = 1./2. * (specificDensityInlet - (particleDisFunc[0, -1, 0] + \
particleDisFunc[1, -1, 0] + particleDisFunc[2, -1, 0] + \
particleDisFunc[3, -1, 0] + particleDisFunc[4, -1, 0] + \
particleDisFunc[6, -1, 0] + particleDisFunc[8, -1, 0]))
particleDisFunc[7, -1, 0] = 1./2. * (specificDensityInlet - (particleDisFunc[0, -1, 0] + \
particleDisFunc[1, -1, 0] + particleDisFunc[2, -1, 0] + \
particleDisFunc[3, -1, 0] + particleDisFunc[4, -1, 0] + \
particleDisFunc[6, -1, 0] + particleDisFunc[8, -1, 0]))
##Outlet bottom
velocityX[0, -1] = 0.
velocityY[0, -1] = 0.
particleDisFunc[3, 0, -1] = particleDisFunc[1, 0, -1]
particleDisFunc[2, 0, -1] = particleDisFunc[4, 0, -1]
particleDisFunc[6, 0, -1] = particleDisFunc[8, 0, -1]
particleDisFunc[5, 0, -1] = 1./2. * (specificDensityOutlet - (particleDisFunc[0, 0, -1] + \
particleDisFunc[1, 0, -1] + particleDisFunc[2, 0, -1] + \
particleDisFunc[3, 0, -1] + particleDisFunc[4, 0, -1] + \
particleDisFunc[6, 0, -1] + particleDisFunc[8, 0, -1]))
particleDisFunc[7, 0, -1] = 1./2. * (specificDensityOutlet - (particleDisFunc[0, 0, -1] + \
particleDisFunc[1, 0, -1] + particleDisFunc[2, 0, -1] + \
particleDisFunc[3, 0, -1] + particleDisFunc[4, 0, -1] + \
particleDisFunc[6, 0, -1] + particleDisFunc[8, 0, -1]))
##outlet top
velocityX[-1, -1] = 0.
velocityY[-1, -1] = 0.
particleDisFunc[3, -1, -1] = particleDisFunc[1, -1, -1]
particleDisFunc[4, -1, -1] = particleDisFunc[2, -1, -1]
particleDisFunc[7, -1, -1] = particleDisFunc[5, -1, -1]
particleDisFunc[6, -1, -1] = 1./ 2. * (specificDensityInlet - (particleDisFunc[0, -1, -1] + \
particleDisFunc[1, -1, -1] + particleDisFunc[2, -1, -1] + particleDisFunc[3, -1, -1] + \
particleDisFunc[4, -1, -1] + particleDisFunc[5, -1, -1] + \
particleDisFunc[7, -1, -1]))
particleDisFunc[8, -1, -1] = 1./ 2. * (specificDensityInlet - (particleDisFunc[0, -1, -1] + \
particleDisFunc[1, -1, -1] + particleDisFunc[2, -1, -1] + particleDisFunc[3, -1, -1] + \
particleDisFunc[4, -1, -1] + particleDisFunc[5, -1, -1] + \
particleDisFunc[7, -1, -1]))
#calculation of the macroscopic parameters
for i in sp.arange(ny):
for j in sp.arange(nx):
rho[i, j] = sp.sum(particleDisFunc[:, i, j])
#print(rho[i,j])
velocityX[i, j] = sp.sum(particleDisFunc[:, i, j] * microVelocity[:, 0]) /\
rho[i, j]
velocityY[i, j] = sp.sum(particleDisFunc[:, i, j] * microVelocity[:, 1]) /\
rho[i,j]
#bounceback boundary for the wall
for i in sp.arange(ny):
for j in sp.arange(nx):
if (isWall[i, j] == True):
tempDis = sp.ones(9)
for k in sp.arange(9):
tempDis[k] = particleDisFunc[k, i, j]
for k in sp.arange(9):
particleDisFunc[oppDirection[k], i, j] = tempDis[k]
#calculate the equilibrium function values
for i in sp.arange(ny):
for j in sp.arange(nx):
velocityAtLattice = np.array([velocityX[i, j], velocityY[i, j]])
for k in sp.arange(9):
equilibriumFunc[k, i, j] = weightCoeff[k] * rho[i, j] * \
(1.0 + eqCoeff1 * np.dot(microVelocity[k], velocityAtLattice) + \
eqCoeff2 * sp.power(np.dot(microVelocity[k], velocityAtLattice), 2.0) +\
eqCoeff3 * np.dot(velocityAtLattice, velocityAtLattice))
#calculate the collision
for i in sp.arange(9):
particleDisFunc[i, :, :] = particleDisFunc[i, :, :] - \
(particleDisFunc[i, :, :] - equilibriumFunc[i, :, :]) / tau