I have developed a python code for a lid-drive cavity model. However, my results are not converging. The algorithm of my code looks like this:
Euler Momentum Equation looks like this:
$$\frac{u^{n+1}-u^{n}}{\Delta t}=-\frac{1}{\rho}\vec{P}^{n+1}-\vec{U}^{n}\cdot\nabla\vec{U}^{n}+\nu\nabla^{2}\vec{U}^{n}$$
- The first step is to update time (current time = $\Delta t$ + current time) and enforce boundary conditions on $U$ and $V$.
- Next step is to discretize the above formula to obtain predicted velocities:
$$\frac{u^{*}-u^{n}}{\Delta t}=-\vec{U}^{n}\cdot\nabla\vec{U}^{n}+\nu\nabla^{2}\vec{U}^{n}\quad\Rightarrow\quad u^{*}=\Delta u^{n}+t\left(\vec{U}^{n}\cdot\nabla\vec{U}^{n}+\nu\nabla^{2}\vec{U}^{n}\right)$$
- The calculated $U^*$ is used to find the pressure values $P^{n+1}$ using the Poisson's equation below:
$$\Delta P^{n+1}=-\frac{\rho}{\Delta t}\nabla \cdot U^{*}$$
- Use $P^{n+1}$ values to find $U^{n+1}$ :
$$u^{n+1}=u^{n}+\Delta t\left(-\frac{1}{\rho}\vec{P}^{n+1}-\vec{U}^{n}\cdot\nabla\vec{U}^{n}+\nu\nabla^{2}\vec{U}^{n}\right)$$
- Enforce boundary conditions, plot velocities and pressure and repeat 1-5 again until the end time.
I think I missed something that it does not show correct results. The source code is shown below:
import numpy as np
import math
import matplotlib.pyplot as plot
#--------------------------------------------------------------------------------------------------------------------------
#1.0 INPUTING MESH VALUES AND CREATING MESH GRID
Lx=1 # length at x direction
dx=0.05 # x step size
Ly=1 # length at y direction
dy=0.05 # y step size
i=math.ceil(Lx/dx) #1.3 Indexing ith cell faces
j=math.ceil(Ly/dy) #1.3 Indexing jth cell faces
#1.4 initializing velocity matrices
v=np.zeros((i+1,j+2))
u=np.zeros((i+2,j+1))
#--------------------------------------------------------------------------------------------------------------------------
#2.0 INPUT INITIAL VALUES
current_time=0
end_time=200
dt=0.025
#--------------------------------------------------------------------------------------------------------------------------
#3.0 INPUT CONSTANTS
mu=1/400 #Input kinematic viscosity
rho=1 # Input density
#--------------------------------------------------------------------------------------------------------------------------
# 4.0 FUNCTIONS
def b_conditions(uu,vv): # TO ESTABLISH BOUNDARY CONDITIONS
vv[0]=0 #top
vv[(len(vv)-1)]=0 #bot
uu[:,0]=0 #left
uu[:,(len(uu[0])-1)]=0 #right
#3.2 Force boundary conditions at face
u_top=1
uu[0]=-uu[1]+2*u_top
u_bot=0
uu[(len(uu)-1)]=-uu[i]+2*u_bot
v_left=0
vv[:,0]=-vv[:,1]+2*v_left
v_right=0
vv[:,(len(vv[0])-1)]=-vv[:,j]+2*v_right
return uu,vv
def velo_tem(u,v,j,i,dx,dy,dt,mu): # TEMPORAL DISCRETIZER FOR VELOCITIES
vs=np.zeros((i+1,j+2))
us=np.zeros((i+2,j+1))
for y in range(1,j):
for x in range(1,i+1):
dudx2=(u[x][y-1]-2*u[x][y]+u[x][y+1])/(dx*dx)
dudy2=(u[x+1][y]-2*u[x][y]+u[x-1][y])/(dy*dy)
ududx=u[x][y]*(u[x][y+1]-u[x][y-1])/(2*dx)
vdudy=0.25*(v[x][y]+v[x][y+1]+v[x-1][y]+v[x-1][y+1])*(u[x-1][y]-u[x+1][y])/(2*dy)
us[x][y]=u[x][y]+dt*(mu*(dudx2+dudy2)-ududx-vdudy)
for y in range(1,j+1):
for x in range(1,i):
dvdx2=(v[x][y-1]-2*v[x][y]+v[x][y+1])/(dx*dx)
dvdy2=(v[x+1][y]-2*v[x][y]+v[x-1][y])/(dy*dy)
udvdx=0.25*(u[x+1][y-1]+u[x][y-1]+u[x+1][y]+u[x][y])*(v[x][y+1]-v[x][y-1])/(2*dx)
vdvdy=v[x][y]*(v[x-1][y]-v[x+1][y])/(2*dy)
vs[x][y]=v[x][y]+dt*(mu*(dvdx2+dvdy2)-udvdx-vdvdy)
b_conditions(us,vs)
return us,vs
#--------------------------------------------------------------------------------------------------------------------------
#5.0 INITIATE LHS FOR LINEAR EQUATION (PRESSURE POISSON'S EQUATION)
#5.1 x direction
lhs_x=-2*np.eye(j,j)+np.eye(j,j,k=1)+np.eye(j,j,k=-1)
lhs_x[0,0]=lhs_x[j-1,j-1]=-1
lhs_x=np.kron(np.eye(i,i),lhs_x)/(dx*dx)
#5.2 y direction
lhs_y=-2*np.eye(i,i)+np.eye(i,i,k=1)+np.eye(i,i,k=-1)
lhs_y[0,0]=lhs_y[i-1,i-1]=-1
lhs_y=np.kron(np.eye(j,j),lhs_y)/(dy*dy)
#5.3 Initiate LHS
ln=j*i
lhs=np.zeros((ln,ln))
xx=np.arange(ln).reshape((i,j))
xx=np.transpose(xx).reshape(ln)
for x in range(ln):
for y in range(ln):
lhs[x,y]=lhs_x[x,y]+lhs_y[xx[x],xx[y]]
#5.4 i dont know the significance of this but was told to put value 1 at P(0,0)
lhs[ln-1]=0
lhs[ln-1,ln-1]=1
#--------------------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------------------
# CALCULATION STARTS HERE
#--------------------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------------------
# Initializing boundary conditionS
b_conditions(u,v)
# Conditions to stop when time reaches end_time
while current_time<=end_time:
#--------------------------------------------------------------------------------------------------------------------------
#CFL control of dt
print("current time: ",current_time)
cfl=dt*(np.max(abs(u))/dx+np.max(abs(v))/dy)
print("cfl: ",cfl)
if cfl>=1:
dt=cfl/(3*(np.max(abs(u))/dx+np.max(abs(v))/dy))
#--------------------------------------------------------------------------------------------------------------------------
#Time update
current_time=current_time+dt
#--------------------------------------------------------------------------------------------------------------------------
#Initializing predictor velocities cells - step 2
us,vs=velo_tem(u,v,j,i,dx,dy,dt,mu)
#--------------------------------------------------------------------------------------------------------------------------
#5.5 Creating Tridiagonal Matrix (RHS) for Poisson Equation
n=0
rhs=np.zeros(ln)
for x in range(i):
for y in range(j):
rhs[n]=((us[x+1,y+1]-us[x+1,y])/dx)+((vs[x,y+1]-vs[x+1,y+1])/dy)
rhs=-rho*rhs/dt
#--------------------------------------------------------------------------------------------------------------------------
#5.6 Solving pressure field
p1=np.linalg.solve(lhs,rhs)
p=np.reshape(p1,(i,j))
#--------------------------------------------------------------------------------------------------------------------------
#6.0 UPDATE VELOCITY FIELD
vn=np.copy(v)
un=np.copy(u)
for y in range(1,j):
for x in range(1,i+1):
dudx2=(u[x][y-1]-2*u[x][y]+u[x][y+1])/(dx*dx)
dudy2=(u[x+1][y]-2*u[x][y]+u[x-1][y])/(dy*dy)
ududx=u[x][y]*(u[x][y+1]-u[x][y-1])/(2*dx)
vdudy=0.25*(v[x][y]+v[x][y+1]+v[x-1][y]+v[x-1][y+1])*(u[x-1][y]-u[x+1][y])/(2*dy)
un[x,y]=u[x,y]-dt*((p[x-1,y]-p[x-1,y-1])/(dx*rho)+ududx+vdudy-mu*(dudx2+dudy2))
for y in range(1,j+1):
for x in range(1,i):
dvdx2=(v[x][y-1]-2*v[x][y]+v[x][y+1])/(dx*dx)
dvdy2=(v[x+1][y]-2*v[x][y]+v[x-1][y])/(dy*dy)
udvdx=0.25*(u[x+1][y-1]+u[x][y-1]+u[x+1][y]+u[x][y])*(v[x][y+1]-v[x][y-1])/(2*dx)
vdvdy=v[x][y]*(v[x-1][y]-v[x+1][y])/(2*dy)
vn[x,y]=v[x,y]-dt*((p[x-1,y-1]-p[x,y-1])/(dy*rho)+udvdx+vdvdy-mu*(dvdx2+dvdy2))
u=un
v=vn
b_conditions(u,v)
#--------------------------------------------------------------------------------------------------------------------------
# VISUALIZATION & PLOTTING
# Interpolation Function
def inter(q11,q21):
return (q11+q21)/2
#// Generation of new loop
#Interpolation and Coordinates of U
#U Visualization
u_i=np.zeros((i,j))
for x in range(i):
for y in range(j):
u_i[x,y]=inter(u[x+1,y],u[x+1,y+1])
#V Visualization
v_i=np.zeros((i,j))
for x in range(i):
for y in range(j):
v_i[x,y]=inter(v[x,y+1],v[x+1,y+1])
#MeshGrid
xx,yy=np.meshgrid(np.linspace(0.5,j-0.5,num=j),np.linspace(0.5,i-0.5,num=i))
yy=np.flip(yy,0)
fig, ax = plot.subplots()
ax.contourf(xx,yy,p)
ax.quiver(xx,yy,u_i,v_i)
plot.show()
