driven.py

import numpy as np
import time
import math

# Module equivalent: Global variables
Nx, Ny = 80, 80  # Spatial grid size
current_time = 0.0  # Renamed from 'time' to avoid conflict with time module
h, dt = 0.0, 0.0  # Spatial and time step
W = np.zeros((Nx+1, Ny+1))  # Vorticity
S = np.zeros((Nx+1, Ny+1))  # Stream function
T = np.zeros((Nx+1, Ny+1))  # Temperature
U = np.zeros((Nx+1, Ny+1))  # Velocity x-component
V = np.zeros((Nx+1, Ny+1))  # Velocity y-component
X = np.zeros(Nx+1)
Y = np.zeros(Ny+1)
Ra = 1e4  # Rayleigh number
Pr = 0.7  # Prandtl number

def main():
    global current_time, h, dt
    print("**********Program of buoyancy-driven cavity flow**********")
    print(f"Related parameter: ")
    print(f"Rayleigh number: Ra={Ra}")
    print(f"Prandtl number: Pr={Pr}")
    print("Solving...")
    
    start = time.time()
    vorticity_stream_method()
    end = time.time()
    print(f"Solving time = {end - start:.4f} s")
    output_tecplot()

def vorticity_stream_method():
    global current_time, h, dt, W, S, U, V, T, X, Y
    
    ESP = 1e-4  # Convergence criterion
    iter_max = int(1e8)
    u_n = np.zeros((Nx+1, Ny+1))
    v_n = np.zeros((Nx+1, Ny+1))
    W_n = np.zeros((Nx+1, Ny+1))
    T_n = np.zeros((Nx+1, Ny+1))
    k1 = np.zeros((Nx-1, Ny-1))
    k2 = np.zeros((Nx-1, Ny-1))
    k3 = np.zeros((Nx-1, Ny-1))
    k4 = np.zeros((Nx-1, Ny-1))
    k5 = np.zeros((Nx-1, Ny-1))
    error_w = np.zeros((Nx+1, Ny+1))
    
    h = (1.0 - 0.0) / Nx  # Uniform grid
    dt = math.ceil(math.sqrt(math.sqrt(Ra/Pr))) * h**2
    
    # Grid setup
    for i in range(Nx+1):
        X[i] = 0.0 + h * i
    for j in range(Ny+1):
        Y[j] = 0.0 + h * j
    
    # Initial conditions
    for j in range(1, Ny):
        for i in range(1, Nx):
            W[i,j] = 0.0
            S[i,j] = 0.0
            U[i,j] = 0.0
            V[i,j] = 0.0
            T[i,j] = (1.0 - Y[j]) * (1.0 + 0.05 * np.sin(2.0 * np.pi * (1.0 - Y[j])))
    
    # Boundary conditions
    boundary()
    
    # Main iteration loop
    iter_count = 0
    current_time = 0.0
    while True:
        # Store previous time step values
        u_n[:,:] = U[:,:]
        v_n[:,:] = V[:,:]
        W_n[:,:] = W[:,:]
        T_n[:,:] = T[:,:]
        
        # Update equation coefficients
        for j in range(1, Ny):
            for i in range(1, Nx):
                k1[i-1,j-1] = 1.0 - 4.0 * dt / (h**2) * np.sqrt(Pr/Ra)
                k2[i-1,j-1] = dt * (1.0/(h**2) * np.sqrt(Pr/Ra) - u_n[i,j]/(2.0*h))
                k3[i-1,j-1] = dt * (1.0/(h**2) * np.sqrt(Pr/Ra) + u_n[i,j]/(2.0*h))
                k4[i-1,j-1] = dt * (1.0/(h**2) * np.sqrt(Pr/Ra) - v_n[i,j]/(2.0*h))
                k5[i-1,j-1] = dt * (1.0/(h**2) * np.sqrt(Pr/Ra) + v_n[i,j]/(2.0*h))
        
        # Solve vorticity equation (FTCS scheme)
        for j in range(1, Ny):
            for i in range(1, Nx):
                W[i,j] = (k1[i-1,j-1] * W_n[i,j] + 
                         k2[i-1,j-1] * W_n[i+1,j] + 
                         k3[i-1,j-1] * W_n[i-1,j] + 
                         k4[i-1,j-1] * W_n[i,j+1] + 
                         k5[i-1,j-1] * W_n[i,j-1] + 
                         dt/(2*h) * (T[i,j+1] - T[i,j-1]))
        
        # Solve stream function equation
        gauss_seidel()
        
        # Update velocity field
        for j in range(1, Ny):
            for i in range(1, Nx):
                U[i,j] = -(S[i,j+1] - S[i,j-1]) / (2.0 * h)
                V[i,j] = (S[i+1,j] - S[i-1,j]) / (2.0 * h)
        
        # Solve temperature field
        calculate_temperature()
        
        # Update boundary conditions
        boundary()
        
        # Check convergence
        error_w = np.abs(W - W_n)
        maxvalue_error_w = np.max(error_w)
        
        iter_count += 1
        current_time += dt
        print(f"step {iter_count} time={current_time:.6f} error={maxvalue_error_w:.6e}")
        
        if maxvalue_error_w < ESP or iter_count > iter_max:
            break

def boundary():
    global U, V, S, W, T
    
    # Bottom boundary
    for i in range(Nx+1):
        U[i,0] = 0.0
        V[i,0] = 0.0
        S[i,0] = 0.0
        W[i,0] = -3.0/h**2 * (S[i,1] - S[i,0]) - 0.5 * W[i,1]
        T[i,0] = 1.0
    
    # Top boundary
    for i in range(Nx+1):
        U[i,Ny] = 0.0
        V[i,Ny] = 0.0
        S[i,Ny] = 0.0
        W[i,Ny] = 3.0/h**2 * (S[i,Ny] - S[i,Ny-1]) - 0.5 * W[i,Ny-1]
        T[i,Ny] = 0.0
    
    # Left boundary
    for j in range(Ny+1):
        U[0,j] = 0.0
        V[0,j] = 0.0
        S[0,j] = 0.0
        W[0,j] = -3.0/h**2 * (S[1,j] - S[0,j]) - 0.5 * W[1,j]
        T[0,j] = (4.0 * T[1,j] - T[2,j]) / 3.0
    
    # Right boundary
    for j in range(Ny+1):
        U[Nx,j] = 0.0
        V[Nx,j] = 0.0
        S[Nx,j] = 0.0
        W[Nx,j] = 3.0/h**2 * (S[Nx,j] - S[Nx-1,j]) - 0.5 * W[Nx-1,j]
        T[Nx,j] = (4.0 * T[Nx-1,j] - T[Nx-2,j]) / 3.0
    
    # Corner points (average method)
    W[0,0] = 0.5 * (W[0,1] + W[1,0])
    W[Nx,0] = 0.5 * (W[Nx-1,0] + W[Nx,1])
    W[Nx,Ny] = 0.5 * (W[Nx,Ny-1] + W[Nx-1,Ny])
    W[0,Ny] = 0.5 * (W[0,Ny-1] + W[1,Ny])

def gauss_seidel():
    global S
    Alpha = 1.50  # Relaxation factor
    ESP = 1e-2    # Convergence criterion
    S_n = np.zeros((Nx+1, Ny+1))
    error_S = np.zeros((Nx+1, Ny+1))
    
    while True:
        # Store previous iteration values
        S_n[:,:] = S[:,:]
        
        # Update stream function
        for j in range(1, Ny):
            for i in range(1, Nx):
                S[i,j] = 0.25 * (S_n[i+1,j] + S[i-1,j] + S_n[i,j+1] + S[i,j-1]) + h**2 * W[i,j] / 4.0
                S[i,j] = Alpha * S[i,j] + (1.0 - Alpha) * S_n[i,j]
        
        # Compute error
        error_S = np.abs(S - S_n)
        maxvalue_error = np.max(error_S)
        
        # Check convergence
        if maxvalue_error < ESP:
            break

def calculate_temperature():
    global T, U, V
    for j in range(1, Ny):
        for i in range(1, Nx):
            dTxx = (T[i+1,j] - 2*T[i,j] + T[i-1,j]) / h**2
            dTyy = (T[i,j+1] - 2*T[i,j] + T[i,j-1]) / h**2
            dTx = (U[i,j] * (T[i+1,j] - T[i-1,j])) / (2*h)
            dTy = (V[i,j] * (T[i,j+1] - T[i,j-1])) / (2*h)
            T[i,j] = T[i,j] + dt * (np.sqrt(1.0/Pr/Ra) * (dTxx + dTyy) - dTx - dTy)

def output_tecplot():
    filename = "result.plt"
    with open(filename, 'w') as f:
        f.write(f'TITLE = " Ra= {Ra:.2f},Pr= {Pr:.2f}"\n')
        f.write('VARIABLES = "X","Y","U","V","T","w","S"\n')
        f.write(f'ZONE I={Nx+1},J={Ny+1} ,DATAPACKING=POINT\n')
        for j in range(Ny+1):
            for i in range(Nx+1):
                f.write(f"{X[i]:16.8E} {Y[j]:16.8E} {U[i,j]:16.8E} {V[i,j]:16.8E} {T[i,j]:16.8E} {W[i,j]:16.8E} {S[i,j]:16.8E}\n")

if __name__ == "__main__":
    main()