added 440 characters in body
Source Link

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

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

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

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

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

edited formulas to be formulas, not images
Source Link
Anton Menshov
  • 8.1k
  • 5
  • 34
  • 90

$$\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. FirstThe first step is to update time (current time = dt$\Delta t$ + current time) and enforce boudaryboundary conditions on U$U$ and V$V$.

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

  1. The calculated is$U^*$ is used to find the pressure values $P^{n+1}$ using the poisson'sPoisson's equation below:

$$\Delta P^{n+1}=-\frac{\rho}{\Delta t}\nabla \cdot U^{*}$$

  1. 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)$$

  1. Enforce boundary conditions, plot velocities and pressure and repeat 1-5 again until the end time.
  1. First step is to update time (current time = dt + current time) and enforce boudary conditions on U and V.
  1. The calculated is used to find the pressure values using the poisson's equation below:
  1. Use values to find :
  1. Enforce boundary conditions, plot velocities and pressure and repeat 1-5 again until end time.

$$\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$.

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

  1. 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^{*}$$

  1. 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)$$

  1. Enforce boundary conditions, plot velocities and pressure and repeat 1-5 again until the end time.
Source Link

Verification on pressure predictor method for CFD code

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:


  1. First step is to update time (current time = dt + current time) and enforce boudary conditions on U and V.

  1. Next step is to discretize the above formula to obtain predicted velocities:

  1. The calculated is used to find the pressure values using the poisson's equation below:

  1. Use values to find :

  1. Enforce boundary conditions, plot velocities and pressure and repeat 1-5 again until 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()