CFD

cfd.py

import matplotlib.pyplot as plt
import numpy as np
import cv2


NACA_0018 =[[1.0000, 0.00000],
            [0.9500, 0.01210],
            [0.9000, 0.02172],
            [0.8000, 0.03935],
            [0.7000, 0.05496],
            [0.6000, 0.06845],
            [0.5000, 0.07941],
            [0.4000, 0.08705],
            [0.3000, 0.09003],
            [0.2500, 0.08912],
            [0.2000, 0.08606],
            [0.1500, 0.08018],
            [0.1000, 0.07024],
            [0.0750, 0.06300],
            [0.0500, 0.05332],
            [0.0250, 0.03922],
            [0.0125, 0.02841],
            [0.0000, 0.00000],
            [0.0125, -0.02841],
            [0.0250, -0.03922],
            [0.0500, -0.05332],
            [0.0750, -0.06300],
            [0.1000, -0.07024],
            [0.1500, -0.08018],
            [0.2000, -0.08606],
            [0.2500, -0.08912],
            [0.3000, -0.09003],
            [0.4000, -0.08705],
            [0.5000, -0.07941],
            [0.6000, -0.06845],
            [0.7000, -0.05496],
            [0.8000, -0.03935],
            [0.9000, -0.02172],
            [0.9500, -0.01210],
            [1.0000, 0.00000]]

def setup_NACA_0018(Nx, Ny, aoa=10, max_height=20):
    # find max x and y from NACA_0018
    x_max = max([x for x, y in NACA_0018])
    y_max = max([y for x, y in NACA_0018])

    scale = max_height / y_max
    x_max_scaled = int(x_max * scale)

    # scale NACA_0018 and draw it to canvas
    NACA_0018_scaled = [(int(x*scale), int(y*scale)+Ny//2) for x, y in NACA_0018]
    airfoil_image = np.zeros((Ny, x_max_scaled))
    cv2.fillPoly(airfoil_image, [np.array(NACA_0018_scaled)], 1)
    
    # rotate cw 30 degree
    M = cv2.getRotationMatrix2D((x_max_scaled//2, Ny//2), aoa, 1)
    airfoil_image = cv2.warpAffine(airfoil_image, M, (x_max_scaled, Ny))
    #cv2.imwrite(f"airfoil.png", airfoil)

    # init empty image
    canvas = np.zeros((Ny, Nx))
    print(f"canvas shape = {canvas.shape}")
    
    canvas[:, 80:x_max_scaled+80] = airfoil_image
    
    #cv2.imwrite(f"output.png", canvas)
    
    return canvas

def setup_airfoil(Nx, Ny, cross_section=None):
    if cross_section is None:
        return setup_cylinder(Nx, Ny)
    
    airfoil_image = setup_NACA_0018(Nx, Ny, -15)

    for y in range(Ny):
        for x in range(Nx):
            if airfoil_image[y, x] == 1:
                cross_section[y, x] = True
    
    return cross_section

def setup_cylinder(Nx, Ny):
    airfoil = np.full((Ny, Nx), False)

    for y in range(Ny):
        for x in range(Nx):
            if distance(Nx//4, Ny//2, x, y) < 13:
                airfoil[y, x] = True
    
    return airfoil

def distance(x1, y1, x2, y2):
    return np.sqrt((x1-x2)**2 + (y1-y2)**2)

def main():
    Nx = 600
    Ny = 300
    #rho0 = 100  # average density
    tau = 0.53  # collision timescale
    Nt = 10000   # number of timesteps

    plot_real_time = True
    plot_every = 50

    #Lattice Boltzmann parameters
    NL = 9
    cxs = np.array([0, 0, 1, 1,  1,  0, -1, -1, -1])
    cys = np.array([0, 1, 1, 0, -1, -1, -1,  0,  1])
    weights = np.array([4/9,1/9,1/36,1/9,1/36,1/9,1/36,1/9,1/36]) 

    # Initial Conditions
    F = np.ones((Ny, Nx, NL)) + 0.01*np.random.randn(Ny ,Nx, NL)

    F[:,:,3] = 2.3
    
    airfoil = setup_airfoil(Nx, Ny, np.full((Ny, Nx), False))
    # Simulation Main Loop
    for it in range(Nt):
        print(it)

        F[:, -1, [6, 7, 8]] = F[:, -2, [6, 7, 8]]   # Right wall: outflow condition
        F[:, 0, [2, 3, 4]] = F[:, 1, [2, 3, 4]]     # Left wall: outflow condition
        F[-1, :, [8, 1, 2]] = F[-2, :, [8, 1, 2]]   # Bottom wall: outflow condition
        F[0, :, [4, 5, 6]] = F[1, :, [4, 5, 6]]     # Top wall: outflow condition
        
        # Drift
        for i, cx, cy in zip(range(NL), cxs, cys):
            F[:,:,i] = np.roll(F[:,:,i], cx, axis=1)
            F[:,:,i] = np.roll(F[:,:,i], cy, axis=0)
        
        # Set reflective boundaries
        bndryF = F[airfoil,:]
        bndryF = bndryF[:,[0,5,6,7,8,1,2,3,4]]

        # Calculate fluid variables
        rho = np.sum(F, 2)
        ux  = np.sum(F*cxs, 2) / rho
        uy  = np.sum(F*cys, 2) / rho

        F[airfoil,:] = bndryF
        ux[airfoil] = 0
        uy[airfoil] = 0

        # Apply Collision
        Feq = np.zeros(F.shape)
        for i, cx, cy, w in zip(range(NL), cxs, cys, weights):
            Feq[:,:,i] = rho * w * ( 
                1 + 3 * (cx*ux + cy * uy) + 9 * (cx*ux + cy*uy)**2 / 2 - 3*(ux**2+uy**2) / 2 
            )
		
        F += -(1.0/tau) * (F - Feq)
        
        # plot in real time - color 1/2 particles blue, other half red
        if (plot_real_time and (it % plot_every) == 0):
            plt.imshow(np.sqrt(ux**2 + uy**2), cmap='viridis')
            plt.pause(0.01)
            plt.cla()
            

    # Save figure
    plt.savefig('latticeboltzmann.png',dpi=240)
    plt.show()
        
    return 0


if __name__ == "__main__":
    main()