I’ve not updated the site in a month or so. Unfortunately, despite my prior article about carving out a study day, life sometimes gets in the way! I’ve been spending a great deal of time finishing off the writing of a book for QuantStart, an entrepreneurial venture of mine in the quantitative finance realm. Despite a busy writing schedule I have spent some time on physics, if only reading in a more sedate manner. I’ve been trying to follow Quantum Field Theory For The Gifted Amateur, which is an absolutely fantastic higher level, albeit still quite mathematically technical, overview of Quantum Field Theory. More on that in later posts.

## My Compressible Computational Fluid Dynamics Research

This post is concerned with Computational Fluid Dynamics (CFD), which I believe to be relevant here because it is fundamentally computational physics, and thus deserves space on the blog! In particular it is a continuation of my postgraduate aeronautics research at Imperial College. I’ve not looked at my PhD work for some time but I thought I would take some of the content I researched and wrap it up into a paper, which I will attempt to submit to a journal. This particular work concerned idealised two-dimensional axisymmetric supersonic air-breathing engine (“ramjet”) intake simulations under conditions that are “off-design”, namely where the compression process is not optimised. This leads to some interesting physical phenomena in the inviscid Euler regime, including strong-shock Mach stems and shock-expansion fan interactions.

Of particular interest to me is a hysteresis effect that occurs in a certain domain of intake cone-angle and free-stream Mach number. My supervisor and I were able to show that such a hysteresis effect can indeed be simulated numerically. Further, we were able to map out the regime in the cone-angle/Mach number parameter space where such hysteresis occurred. The domain was bounded by the so-called Detachment Criterion and Von-Neumann Criterion.

The goal of this article is to describe how I plan to replicate this result with more modern tools, e.g. Python and C++, and then extend the result to multiple engine geometries, potentially including curved surfaces within the intakes. This would allow an isentropic intake cone as well as a curved cowl to replicate the shock-on-lip and shock-on-shoulder conditions.

There are two components to the simulation. The first will be the design and implementation of a finite-volume method (FVM) compressible inviscid solver, using a Godunov-type method and a Riemann solver. My supervisor previously had a sophisticated Fortran RANS-based turbulence code that utilised the Generalised Riemann Problem of Ben Artzi and Falcovitz^{[1]}. I plan to build my own in Python, due to reasons of development speed and testing. Once I am happy with the results, I will replicate the code in C++ and optimise, making sure to compare results with the Python code at every step in order to ensure validity. This code will solve the two-dimensional axisymmetric Euler equations in an inhomogeneous “stretched” block-structured fashion, which makes the assumption that the tangent free-stream flow vector is parallel to the cone axis of symmetry, i.e. zero angle of attack. The output will consist of a velocity, density and pressure fields. A suitable visualisation tool will also be needed in the absence of the commercial TecPlot code!

The second code will utilise the Method of Characteristics (MoC) to produce a highly performant analytic solution to flow-field values at suitable domains within an engine. I had previously written such a code in C++, known as the Axisymmetric Method of Characteristics Solver, or ARMOCS. I aim to re-implement this in Python, since execution speed will be more than sufficient for this task. The goal is then to produce the domain of duality that sits between the Detachment Criterion and the Von Neumann Criterion in free-stream Mach number/cone-angle space. Once these boundaries have been established, the FVM code can perform the actual simulation to check whether the hysteresis effect manifests itself. I have already established that this is the case for the particular geometry I studied for my PhD, but I plan to consider more sophisticated regimes and make the code more general.

Much has changed since I carried out the research in 2008 (six years ago!). In the interim more Riemann Solvers have been developed. Hence I have more scope to “swap out” the underlying solver being used in the FVM simulation. In addition, computational speed and random-access memory have moved on substantially. Thus I will able to try more highly resolved meshes, albeit with problems identified by Quirk^{[3]}. I have not yet carried out an up-to-date literature review on the latest shockwave/expansion-fan literature, so the first step will be to bring myself up to speed.

## Python FVM Godunov-Type Roe-Based Euler Solver

I’ve already begun the process by constructing a one-dimensional Roe-based Godunov solver for the inviscid Euler equations in one-dimension. An interval $[0,1] \in \mathbb{R}$ is divided into $N$ cells, each of which contains a discretised pressure, velocity and density field. Such a grouping of the three fields is known as a *primitive* formulation. The Godunov procedure of updating the relations is recast into the *conservative* formulation, which essentially provides a (discretised) density field, a momentum field and an internal energy field. The conservative formulation of the Euler equations is given by:

\begin{eqnarray}

\frac{\partial U}{\partial t} = \frac{\partial F(U)}{\partial x}

\end{eqnarray}

Where $F(U) = (\rho, \rho u, )$ are the fluxes. The numerical evolution of the system is governed by the Godunov updating procedure:

\begin{eqnarray}

U_i^{n+1} = \frac{\Delta t}{\Delta x} (F(U_{i+\frac{1}{2}}^n) – F(U_{i-\frac{1}{2}}^n))

\end{eqnarray}

For this simple situation we are considering a modified version of the classic Sod Test, which is found in Chapter 6 of Toro^{[4]}. We choose a CFL number of 0.9, a diaphragm point of the “Riemann shock tube” of $x = 0.3$, for the interval $[0,1]$, a primitive left state $W_L$ of $(\rho_L, u_L, p_L) = (1.0, 0.75, 1.0)$ and a primitive right state $W_R$ of $(\rho_R, u_R, p_R) = (0.1, 0.125, 0)$.

The first step in the code is to import the necessary Python libraries. We simply require NumPy and os.path:

import os.path import numpy as np

It is useful to have functions that can convert between primitive states $W$ and conservative variable states $U$:

def w2u(w, gamma): """ Convert the primitive to conservative variables. """ u = np.zeros(3) u[0] = w[0] u[1] = w[0]*w[1] u[2] = w[2]/(gamma-1.0)+0.5*w[0]*w[1]*w[1] return u def u2w(u, gamma): """ Convert the conservative to primitive variables. """ w = np.zeros(3) w[0] = u[0] w[1] = u[1]/u[0] w[2] = (gamma-1.0)*( u[2] - 0.5*w[0]*w[1]*w[1] ) return w

The flux of the conservative variables is calculated using the Euler equations:

def euler_flux(w, gamma): """ Calculate the conservative Euler fluxes. """ rho = w[0] u = w[1] p = w[2] a2 = gamma*p/rho f_1 = rho*u f_2 = rho*u*u + p f_3 = rho*u*( a2/(gamma-1.0) + 0.5*u*u ) return np.array([f_1, f_2, f_3])

The majority of the code is concerned with an implementation of a Roe-based approximate Riemann solver:

def roe_flux(wL, wR, gamma): """ Use the Roe approximate Riemann solver to calculate fluxes. """ uL = w2u(wL, gamma) uR = w2u(wR, gamma) # Primitive and other variables. # Left state rhoL = wL[0] vL = wL[1] pL = wL[2] aL = np.sqrt(gamma*pL/rhoL) HL = ( uL[2] + pL ) / rhoL # Right state rhoR = wR[0] vR = wR[1] pR = wR[2] aR = np.sqrt(gamma*pR/rhoR) HR = ( uR[2] + pR ) / rhoR # First compute the Roe Averages RT = np.sqrt(rhoR/rhoL); rho = RT*rhoL v = (vL+RT*vR)/(1.0+RT) H = (HL+RT*HR)/(1.0+RT) a = np.sqrt( (gamma-1.0)*(H-0.5*v*v) ) # Differences in primitive variables. drho = rhoR - rhoL du = vR - vL dP = pR - pL # Wave strength (Characteristic Variables). dV = np.array([0.0,0.0,0.0]) dV[0] = 0.5*(dP-rho*a*du)/(a*a) dV[1] = -( dP/(a*a) - drho ) dV[2] = 0.5*(dP+rho*a*du)/(a*a) # Absolute values of the wave speeds (Eigenvalues) ws = np.array([0.0,0.0,0.0]) ws[0] = abs(v-a) ws[1] = abs(v) ws[2] = abs(v+a) # Modified wave speeds for nonlinear fields (the so-called entropy fix, which # is often implemented to remove non-physical expansion shocks). # There are various ways to implement the entropy fix. This is just one # example. Try turn this off. The solution may be more accurate. Da = max(0.0, 4.0*((vR-aR)-(vL-aL)) ) if (ws[0] < 0.5*Da): ws[0] = ws[0]*ws[0]/Da + 0.25*Da Da = max(0.0, 4.0*((vR+aR)-(vL+aL)) ) if (ws[2] < 0.5*Da): ws[2] = ws[2]*ws[2]/Da + 0.25*Da # Right eigenvectors R = np.zeros((3,3)) R[0][0] = 1.0 R[1][0] = v - a R[2][0] = H - v*a R[0][1] = 1.0 R[1][1] = v R[2][1] = 0.5*v*v R[0][2] = 1.0 R[1][2] = v + a R[2][2] = H + v*a # Compute the average flux. flux = 0.5*( euler_flux(wL, gamma) + euler_flux(wR, gamma) ) # Add the matrix dissipation term to complete the Roe flux. for i in range(0,3): for j in range(0,3): flux[i] = flux[i] - 0.5*ws[j]*dV[j]*R[i][j] return flux

We need to enforce a CFL condition to maintain numerical stability:

def calc_time_step(cfl, dx, gamma, bcells, U): """ Calculates the maximum wavespeeds and thus the timestep via an enforced CFL condition. """ max_speed = -1.0 for i in range(1,bcells-1): w = u2w(U[i], gamma) u = w[1] c = np.sqrt(gamma*w[2]/w[0]) max_speed = max(max_speed, abs(u)+c) dt = cfl*dx/max_speed # CFL condition return dt

We then carry out the Godunov temporal updating procedure:

def update_solution(U, fluxes, dt, dx, gamma, bcells): """ Updates the solution of the equation via the Godunov procedure. """ # Create fluxes for i in range(0, bcells-1): wL = u2w(U[i], gamma) wR = u2w(U[i+1], gamma) fluxes[i] = roe_flux(wL, wR, gamma) # Update solution for i in range(1, bcells-1): U[i] = U[i] + (dt/dx) * (fluxes[i-1]-fluxes[i]) # Boundary Conditions U[0] = U[1] U[bcells-1] = U[bcells-2]

Finally we wrap the previous functions in a `__main__`

block and output the primitive variables to disk:

if __name__ == "__main__": gamma = 1.4 cells = 1000 bcells = cells + 2 dx = 1.0/cells cfl = 0.8 t = 0.0 tf = 0.2 nsteps = 0 U = np.zeros((bcells,3)) fluxes = np.zeros((bcells,3)) test_case_sod(U, bcells, gamma) for n in range(1, 50000): if (t==tf): break dt = calc_time_step(cfl, dx, gamma, bcells, U) if (t+dt > tf): dt = tf - t update_solution(U, fluxes, dt, dx, gamma, bcells) t += dt nsteps += 1 out_csv = open("output/1d_roe_godunov_euler_sod.csv", "w") for elem in U: new_elem = u2w(elem, gamma) out_csv.write("%s,%s,%s\n" % tuple(new_elem))

We now plot a three-pane of density, velocity and pressure. Comparison with the Sod Test in Toro^{[4]} gives full agreement for $N=1000$. The matplotlib code for the output is as follows:

import sys import matplotlib.pyplot as plt import pandas as pd if __name__ == "__main__": filename = sys.argv[1] data = pd.io.parsers.read_csv( filename, delimiter=",", header=None, names=["rho", "u", "p"], ) # Plot three charts fig = plt.figure() fig.patch.set_facecolor('white') # Set the outer colour to white # Plot the density ax1 = fig.add_subplot(311, ylabel='Density', ylim=[0.0, 1.2]) data['rho'].plot(ax=ax1, color="black", lw=2.) # Plot the velocity ax2 = fig.add_subplot(312, ylabel='Velocity', ylim=[0.0, 1.5]) data['u'].plot(ax=ax2, color="black", lw=2.) # Plot the pressure ax3 = fig.add_subplot(313, ylabel='Pressure', ylim=[0.0, 1.2]) data['p'].plot(ax=ax3, color="black", lw=2.) # Plot the figure plt.show()

The Sod Test visualisation is given below:

While this is a "good start" there are many improvements needed before being able to simulate a full isentropic supersonic intake:

- I need to consider the two-dimensional axisymmetric Euler equations on a homogeneous grid. For this I will employ Strang's operator splitting
^{}technique, that allows a one-dimensional "sweep" in each spatial direction. - I need to use a more robust second-order Riemann solver, such as the GRP. This will take some time to implement as it is a non-trivial algorithm. Thankfully the Ben-Artzi paper
^{[]}and monograph^{}provides a detailed overview of how to carry it out. I would also like to try the HLLE/HLLC solvers of Harten, Lax, van Leer, Einfeldt and Toro. - While I will maintain a structured grid, I would like to add grid inhomogeneity, such that certain areas can receive a higher resolution than others. This is appropriate at the cowl lip and cowl shoulder. Eventually I would like to use a fully unstructured adaptive mesh refinement code, but this is some time away!
- The grid will of course need to have a non-Cartesian geometry to simulate the isometric intake and curved cowl.

In the next post I will discuss the development of the 2D axisymmetric solver using operator splitting.