Archivi categoria: CFD

Uniform flow past a doublet…

Flow over a cylinder

A doublet alone does not give so much information about how it can be used to represent a practical flow pattern in aerodynamics. But let’s use our superposition powers: our doublet in a uniform flow turns out to be a very interesting flow pattern. Let’s first define a uniform horizontal flow.

u_inf = 1.0        # freestream speed

Now, we can calculate velocities and stream function values for all points in our grid. And as we now know, we can calculate them all together with one line of code per array.

u_freestream = u_inf * numpy.ones((N, N), dtype=float)
v_freestream = numpy.zeros((N, N), dtype=float)

psi_freestream = u_inf * Y

Below, the stream function of the flow created by superposition of a doublet in a free stream is obtained by simple addition.

The plot shows that this pattern can represent the flow around a cylinder with center at the location of the doublet. All the streamlines remaining outside the cylinder originated from the uniform flow. All the streamlines inside the cylinder can be ignored and this area assumed to be a solid object. This will turn out to be more useful than you may think.

# superposition of the doublet on the freestream flow
u = u_freestream + u_doublet
v = v_freestream + v_doublet
psi = psi_freestream + psi_doublet

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')

# calculate the stagnation points
x_stagn1, y_stagn1 = +math.sqrt(kappa / (2 * math.pi * u_inf)), 0.0
x_stagn2, y_stagn2 = -math.sqrt(kappa / (2 * math.pi * u_inf)), 0.0

# display the stagnation points
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');

Bernoulli’s equation and the pressure coefficient

A very useful measurement of a flow around a body is the coefficient of pressureCp. To evaluate the pressure coefficient, we apply Bernoulli’s equation for ideal flow and with simple mathematical steps:

\(C_{p}=1-\left(\frac{U}{U_{\infty}}\right)^{2}\)

In an incompressible flow, Cp=1 at a stagnation point. Let’s plot the pressure coefficient in the whole domain.

# compute the pressure coefficient field
cp = 1.0 - (u**2 + v**2) / u_inf**2

# plot the pressure coefficient field
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(1.1 * width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
contf = pyplot.contourf(X, Y, cp,
                        levels=numpy.linspace(-2.0, 1.0, 100), extend='both')
cbar = pyplot.colorbar(contf)
cbar.set_label('$C_p$', fontsize=16)
cbar.set_ticks([-2.0, -1.0, 0.0, 1.0])
pyplot.scatter(x_doublet, y_doublet,
               color='#CD2305', s=80, marker='o')
pyplot.contour(X,Y,psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');
pyplot.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')

Below, we report the complete Python code:

import math as mt
import numpy as np
from matplotlib import pyplot as pyp

N = 50                                # Number of points in each direction
x_start, x_end = -2.0, 2.0            # x-direction boundaries
y_start, y_end = -1.0, 1.0            # y-direction boundaries
x = np.linspace(x_start, x_end, N)    # creates a 1D-array for x
y = np.linspace(y_start, y_end, N)    # creates a 1D-array for y
X, Y = np.meshgrid(x, y)              # generates a mesh grid

kappa = 1.0                        # strength of the doublet
x_doublet, y_doublet = 0.0, 0.0    # location of the doublet

def get_velocity_doublet(strength, xd, yd, X, Y):
    
    u = (- strength / (2 * mt.pi) *
         ((X - xd)**2 - (Y - yd)**2) /
         ((X - xd)**2 + (Y - yd)**2)**2)
    v = (- strength / (2 * mt.pi) *
         2 * (X - xd) * (Y - yd) /
         ((X - xd)**2 + (Y - yd)**2)**2)
    
    return u, v

def get_stream_function_doublet(strength, xd, yd, X, Y):
    
    psi = - strength / (2 * mt.pi) * (Y - yd) / ((X - xd)**2 + (Y - yd)**2)
    
    return psi

# compute the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)

# compute the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
pyp.streamplot(X, Y, u_doublet, v_doublet,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyp.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o');



u_inf = 1.0        # freestream speed

u_freestream = u_inf * np.ones((N, N), dtype=float)
v_freestream = np.zeros((N, N), dtype=float)

psi_freestream = u_inf * Y

# superposition of the doublet on the freestream flow
u = u_freestream + u_doublet
v = v_freestream + v_doublet
psi = psi_freestream + psi_doublet

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
pyp.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyp.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyp.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')

# calculate the stagnation points
x_stagn1, y_stagn1 = +mt.sqrt(kappa / (2 * mt.pi * u_inf)), 0.0
x_stagn2, y_stagn2 = -mt.sqrt(kappa / (2 * mt.pi * u_inf)), 0.0

# display the stagnation points
pyp.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');

# compute the pressure coefficient field
cp = 1.0 - (u**2 + v**2) / u_inf**2

# plot the pressure coefficient field
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(1.1 * width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
contf = pyp.contourf(X, Y, cp,
                        levels=np.linspace(-2.0, 1.0, 100), extend='both')
cbar = pyp.colorbar(contf)
cbar.set_label('$C_p$', fontsize=16)
cbar.set_ticks([-2.0, -1.0, 0.0, 1.0])
pyp.scatter(x_doublet, y_doublet,
               color='#CD2305', s=80, marker='o')
pyp.contour(X,Y,psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyp.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');
pyp.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')


pyp.show()

Solving 2D Heat Equation…

Before we do the Python code, let’s talk about the heat equation and finite-difference method. Heat equation is basically a partial differential equation, it is

If we want to solve it in 2D (Cartesian), we can write the heat equation above like this:

where u is the quantity that we want to know, t is for temporal variable, x and y are for spatial variables, and α is diffusivity constant. So basically we want to find the solution u everywhere in x and y, and over time t.

We can write the heat equation above using finite-difference method like this:

If we arrange the equation above by taking Δx = Δy, we get this final equation:

where

\(\gamma = \alpha \frac{{\Delta t}}{{\Delta {x^2}}}\)

We use explicit method to get the solution for the heat equation, so it will be numerically stable whenever

\(\Delta t \leq \frac{{\Delta {x^2}}}{{4\alpha }}\)

Everything is ready. Now we can solve the original heat equation approximated by algebraic equation above, which is computer-friendly.

Let’s suppose a thin square plate with the side of 50 unit length. The temperature everywhere inside the plate is originally 0 degree (at t = 0), let’s see the diagram below:

For our model, let’s take Δx = 1 and α = 2.0.
Now we can use Python code to solve this problem numerically to see the temperature everywhere (denoted by i and j) and over time (denoted by k). Let’s first import all of the necessary libraries, and then set up the boundary and initial conditions.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

print("2D heat equation solver")

plate_length = 50
max_iter_time = 750

alpha = 2
delta_x = 1

delta_t = (delta_x ** 2)/(4 * alpha)
gamma = (alpha * delta_t) / (delta_x ** 2)

# Initialize solution: the grid of u(k, i, j)
u = np.empty((max_iter_time, plate_length, plate_length))

# Initial condition everywhere inside the grid
u_initial = 0

# Boundary conditions
u_top = 100.0
u_left = 0.0
u_bottom = 0.0
u_right = 0.0

# Set the initial condition
u.fill(u_initial)

# Set the boundary conditions
u[:, (plate_length-1):, :] = u_top
u[:, :, :1] = u_left
u[:, :1, 1:] = u_bottom
u[:, :, (plate_length-1):] = u_right

We’ve set up the initial and boundary conditions, let’s write the calculation function based on finite-difference method that we’ve derived above.

def calculate(u):
    for k in range(0, max_iter_time-1, 1):
        for i in range(1, plate_length-1, delta_x):
            for j in range(1, plate_length-1, delta_x):
                u[k + 1, i, j] = gamma * (u[k][i+1][j] + u[k][i-1][j] + u[k][i][j+1] + u[k][i][j-1] - 4*u[k][i][j]) + u[k][i][j]

    return u

Let’s prepare the plot function so we can visualize the solution (for each k) as a heat map. We use Matplotlib library, it’s easy to use.

def plotheatmap(u_k, k):
    # Clear the current plot figure
    plt.clf()

    plt.title(f"Temperature at t = {k*delta_t:.3f} unit time")
    plt.xlabel("x")
    plt.ylabel("y")

    # This is to plot u_k (u at time-step k)
    plt.pcolormesh(u_k, cmap=plt.cm.jet, vmin=0, vmax=100)
    plt.colorbar()

    return plt

One more thing that we need is to animate the result because we want to see the temperature points inside the plate change over time. So let’s create the function to animate the solution.

def animate(k):
  plotheatmap(u[k], k)

anim = animation.FuncAnimation(plt.figure(), animate, interval=1, frames=max_iter_time, repeat=False)
anim.save("heat_equation_solution.gif")

That’s it! And here’s the result:

1-D Linear Convection

The 1-D Linear Convection equation is the simplest, most basic model that can be used to learn something about CFD.

\(\frac{{\partial u}}{{\partial t}} + c\frac{{\partial u}}{{\partial x}} = 0\)

With given initial conditions (understood as a wave), the equation represents the propagation of that initial wave with speed \(c\), without change of shape.

Let the initial condition be \(u(x,0)=u_0(x)\). Then the exact solution of the equation is \(u(x,t)=u_0(x-ct)\).

We discretize this equation, in both space and time, using the Forward Difference scheme for the time derivative and the Backward Difference scheme for the space derivative.

Our discrete equation is:

\(\frac{{u_i^{n + 1} – u_i^n}}{{\Delta t}} + c\frac{{u_i^n – u_{i – 1}^n}}{{\Delta x}} = 0\)

Now let’s try implementing this in Python.

We’ll start by importing a few libraries to help us out.

  • numpy is a library that provides a bunch of useful matrix operations akin to MATLAB;
  • matplotlib is a 2D plotting library that we will use to plot our results;
  • time and sys provide basic timing functions that we’ll use to slow down animations for viewing.
# Remember: comments in python are denoted by the pound sign
import numpy as np                      #here we load numpy
from matplotlib import pyplot as plt      #here we load matplotlib
import time, sys                   #and load some utilities



nx = 41         # spatial discretization: number of grid points
dx = 2 / (nx-1) # the distance between any pair of adjacent grid points
nt = 25         # nt is the number of timesteps we want to calculate
dt = .025       # dt is the amount of time each timestep covers (delta t)
c = 1           # assume wavespeed of c = 1

We have define an evenly spaced grid of points within a spatial domain that is 2 units of length wide, i.e., \(x_i \in [0,2]\). We have define a variable nx, which represent the number of grid points we want and dx is the distance between any pair of adjacent grid points.

We also need to set up our initial conditions. The initial velocity \(u_0\) is given as \(u=2\) in the interval \(0.5 \leq x \leq 1\) and \(u=1\) everywhere else in \((0,2)\) (i.e., a hat function). Here, we use the function ones() defining a numpy array which is nx elements long with every value equal to 1.

u = np.ones(nx)                      #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s

x = np.linspace(0,2,nx)
plt.plot(x,u, 'red', label='Initial shape')

plt.show()

Now it’s time to implement the discretization of the convection equation using a finite-difference scheme.

We’ll store the result in a new (temporary) array un, which will be the solution \(u\) for the next time-step. We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has convected.

We first initialize our placeholder array un to hold the values we calculate for the \(n+1\) timestep, using once again the NumPy function ones().

Then, we may think we have two iterative operations: one in space and one in time (we’ll learn differently later), so we’ll start by nesting one loop inside the other. Note the use of the nifty range() function. When we write: for i in range(1,nx) we will iterate through the u array, but we’ll be skipping the first element (the zero-th element).

un = np.ones(nx) #initialize a temporary array

for n in range(nt):        #loop for values of n from 0 to nt, so it will run nt times
    un = u.copy()          ##copy the existing values of u into un
    for i in range(1, nx): ## you can try commenting this line and...
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
    

plt.plot(x,u, 'green', label='Final shape')

plt.xlabel('space')
plt.ylabel('value')
pplt.title('1-D Linear Convection : nx=41, nt=25, dt=.025 and c=1')


plt.legend()

plt.show()

OK! So our hat function has definitely moved to the right, but it’s no longer a hat. What’s going on?