2
$\begingroup$

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}$$


1) The first step is to update time (current time = $\Delta t$ + current time) and enforce boundary conditions on $U$ and $V$.


2) 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^{*}= u^{n}-\Delta t\left(\vec{U}^{n}\cdot\nabla\vec{U}^{n}-\nu\nabla^{2}\vec{U}^{n}\right)$$


3) 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^{*}$$


4) 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)$$


5) 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()

Editted:

Boundary conditions are all zero except for u=1 at the top (v left=v right=v top=v bottom=0, u left=u right=u bottom=0). I use neuman boundary conditions for pressure so pressure derivative normal to grid walls = 0.

I've screenshot a plot. It does not show any influence from V, thus there's no swirling of the fluid in the cavity.

enter image description here

$\endgroup$
  • 2
    $\begingroup$ A couple of clarifications: a) what is $V$ boundary condition entering sentence in 1)? $V$ does not seem to participate in the Euler momentum equation. b) Can you add a plot of the non-converging results? $\endgroup$ – Anton Menshov Aug 5 '19 at 17:05
  • $\begingroup$ Could you also confirm which velocity component corresponds to which physical coordinate, and on which grid? at the moment your $d/dx$s seem to be happening over the variable you've labelled as $y$, and vice versa, which is somewhat confusing. $\endgroup$ – origimbo Aug 5 '19 at 17:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.