Navier-Stokes Equation Simulation with SIMPLE algorithm (transcient state)

gistfile1.txt

# Import necessary libraries
import numpy as np

# Define fluid properties and boundary conditions
density = 1.0
viscosity = 1.0
pressure_inlet = 1.0
velocity_inlet = 1.0
pressure_outlet = 0.0

# Define grid size and time step
nx = 100
ny = 100
nz = 100
dt = 0.01

# Initialize velocity and pressure fields
u = np.zeros((nx, ny, nz))
v = np.zeros((nx, ny, nz))
w = np.zeros((nx, ny, nz))
p = np.zeros((nx, ny, nz))

# Iterate over time steps
for t in range(num_timesteps):
  # Calculate intermediate velocity field
  u_star = u + dt * ((1.0 / density) * (p[1:,:,:] - p[:-1,:,:]) - viscosity * ((u[1:,:,:] - u[:-1,:,:]) / dx + (v[:,1:,:] - v[:,:-1,:]) / dy + (w[:,:,1:] - w[:,:,:-1]) / dz))
  v_star = v + dt * ((1.0 / density) * (p[:,1:,:] - p[:,:-1,:]) - viscosity * ((u[:,1:,:] - u[:,:-1,:]) / dx + (v[1:,:,:] - v[:-1,:,:]) / dy + (w[:,:,1:] - w[:,:,:-1]) / dz))
  w_star = w + dt * ((1.0 / density) * (p[:,:,1:] - p[:,:,:-1]) - viscosity * ((u[:,1:,:] - u[:,:-1,:]) / dx + (v[:,:,1:] - v[:,:,:-1]) / dy + (w[1:,:,:] - w[:-1,:,:]) / dz))


  # Apply boundary conditions
  u_star[0,:,:] = velocity_inlet  # Inlet
  u_star[-1,:,:] = 0  # Outlet
  v_star[:,0,:] = 0  # Wall
  v_star[:,-1,:] = 0  # Wall
  w_star[:,:,0] = 0  # Wall
  w_star[:,:,-1] = 0  # Wall


  # Calculate pressure
  p_new = p.copy()
  for i in range(nx):
    for j in range(ny):
      for k in range(nz):
        p_new[i,j,k] = (1.0 - omega) * p[i,j,k] + omega * (p[i+1,j,k] + p[i-1,j,k] + p[i,j+1,k] + p[i,j-1,k] + p[i,j,k+1] + p[i,j,k-1] - (dx**2) * (u_star[i,j,k] - u_star[i-1,j,k]) / dt - (dy**2) * (v_star[i,j,k] - v_star[i,j-1,k]) / dt - (dz**2) * (w_star[i,j,k] - w_star[i,j,k-1]) / dt) / (2.0 / (dx**2) + 2.0 / (dy**2) + 2.0 / (dz**2))
  p = p_new
  
  # Calculate final velocity field
  u = u_star - (dt / density) * (p[1:,:,:] - p[:-1,:,:]) / dx
  v = v_star - (dt / density) * (p[:,1:,:] - p[:,:-1,:]) / dy
  w = w_star - (dt / density) * (p[:,:,1:] - p[:,:,:-1]) / dz
  
  # Update time step
  t += dt

# Output results
# Save velocity and pressure fields to file
np.savetxt('velocity_x.txt', u)
np.savetxt('velocity_y.txt', v)
np.savetxt('velocity_z.txt', w)
np.savetxt('pressure.txt', p)

# End simulation
print("Simulation complete.")