How to Use NumPy for Computational Fluid Dynamics

Introduction

Computational Fluid Dynamics (CFD) is a branch of fluid mechanics that utilizes numerical analysis and algorithms to solve and analyze problems involving fluid flows. NumPy, a fundamental package for scientific computing with Python, is widely used in various fields, including CFD, due to its powerful n-dimensional array object and a collection of routines for fast operations on arrays. This tutorial explores the essentials of using NumPy in CFD applications.

Getting Started with NumPy

Before delving into CFD, you need to install NumPy. Use the following command to install NumPy if you haven’t already:

pip install numpy

Once installed, import NumPy in your Python script as follows:

import numpy as np

NumPy Basics for CFD

CFD calculations often involve dealing with multidimensional arrays representing physical quantities over a spatial domain. Let’s start by creating a 2D grid using NumPy:

# Define the domain dimensions
x_len = 10
y_len = 10

# Generate a 2D grid
x = np.linspace(0, x_len, num=x_len+1)
y = np.linspace(0, y_len, num=y_len+1)
XX, YY = np.meshgrid(x, y)

This code snippet creates a 10×10 grid using linspace and meshgrid functions which is the foundation for our CFD domain.

Discretization and the Finite Difference Method

The finite difference method (FDM) is a numerical technique for approximating derivatives, which is fundamental in simulating fluid dynamics. Let’s discretize the 2D Laplace equation as an example:

# Initialize the scalar field
phi = np.zeros((y_len+1, x_len+1))

# Boundary conditions
phi[:, 0] = 0   # Left wall
phi[:, -1] = 1  # Right wall

# Discretization using the finite difference
for k in range(1000): # Iterations
    for i in range(1, x_len):
        for j in range(1, y_len):
            phi[j, i] = 0.25 * (phi[j, i-1] + phi[j, i+1] +
                               phi[j-1, i] + phi[j+1, i])

# Visualizing the results with matplotlib
import matplotlib.pyplot as plt
plt.contourf(XX, YY, phi)
plt.colorbar()
plt.title('2D Laplace Equation Solution')
plt.show()

This example demonstrates a simple implementation of the FDM to solve the Laplace equation.

Advanced CFD Examples with NumPy

As the complexity of fluid dynamics problems increases, so does the complexity of the numerical methods required to solve them. Let’s look at an example of the Navier-Stokes equations, which govern the motion of fluid substances.

The Navier-Stokes equations describe the motion of fluid substances and are a set of nonlinear partial differential equations. A common approach to solve them, especially in a learning or experimental context, is to use a method like the finite difference method for discretization.

Assuming this is a 2D, incompressible flow, and we’re looking at a simple scenario (like flow in a lid-driven cavity), a basic implementation could look like the following. Keep in mind that this is a very simplified version and real-world fluid dynamics problems are often much more complex and require more sophisticated methods and boundary conditions:

import numpy as np
import matplotlib.pyplot as plt

# Initialization of flow field variables
nx, ny = 50, 50
nt, nu, dt = 100, 0.1, 0.001
dx, dy = 2 / (nx - 1), 2 / (ny - 1)

u = np.zeros((ny, nx))
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))
b = np.zeros((ny, nx))

def build_up_b(b, rho, dt, u, v, dx, dy):
    b[1:-1, 1:-1] = (rho * (1 / dt * 
                    ((u[1:-1, 2:] - u[1:-1, 0:-2]) / 
                     (2 * dx) + (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
                    ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
                      2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
                           (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
                          ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2))

    return b

def pressure_poisson(p, dx, dy, b):
    pn = np.empty_like(p)
    pn = p.copy()
    
    for q in range(nit):
        pn = p.copy()
        p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 + 
                          (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
                          (2 * (dx**2 + dy**2)) -
                          dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * 
                          b[1:-1, 1:-1])

        # Boundary conditions for pressure
        p[:, -1] = p[:, -2] # dp/dx = 0 at x = 2
        p[0, :] = p[1, :]   # dp/dy = 0 at y = 0
        p[:, 0] = p[:, 1]   # dp/dx = 0 at x = 0
        p[-1, :] = 0        # p = 0 at y = 2
        
    return p

rho = 1
nit = 50 

for n in range(nt):
    un = u.copy()
    vn = v.copy()
    
    b = build_up_b(b, rho, dt, u, v, dx, dy)
    p = pressure_poisson(p, dx, dy, b)
    
    u[1:-1, 1:-1] = (un[1:-1, 1:-1]-
                     un[1:-1, 1:-1] * dt / dx *
                    (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
                     vn[1:-1, 1:-1] * dt / dy *
                    (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
                     dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +
                     nu * (dt / dx**2 *
                    (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                     dt / dy**2 *
                    (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))

    v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                     un[1:-1, 1:-1] * dt / dx *
                    (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
                     vn[1:-1, 1:-1] * dt / dy *
                    (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
                     dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +
                     nu * (dt / dx**2 *
                    (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
                     dt / dy**2 *
                    (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))

    # Boundary conditions for velocity
    u[0, :], u[-1, :], u[:, 0], u[:, -1] = 0, 0, 0, 0
    v[0, :], v[-1, :], v[:, 0], v[:, -1] = 0, 0, 0, 0

# Plotting the velocity field (u, v)
plt.figure(figsize=(11, 7))
plt.quiver(X, Y, u, v)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

This code provides a basic framework for solving the Navier-Stokes equations in a 2D space using finite difference methods. It includes:

• The build-up of the source term b.
• The iterative solution of the pressure Poisson equation.
• The update of the velocity field u and v.
• Simple boundary conditions.

Performance Optimization in NumPy

In CFD simulations, performance is key. NumPy’s ability to handle vectorized operations can significantly speed up computation. Always use built-in functions and avoid loops where possible for optimal performance:

# Bad performance
for i in xrange(y_len):
    for j in xrange(x_len):
        phi[j, i] = some_operation(phi[j, i])

# Good performance
phi = some_operation_vectorized(phi)

Using vectorized functions takes advantage of NumPy’s optimized C backends.

Conclusion

NumPy offers powerful tools to handle the numerical computations required for CFD. By leveraging NumPy arrays and optimized operations, you can perform complex fluid flow simulations with ease. Continue exploring various numerical methods and performance optimizations to further enhance your CFD projects.