Hamza H — Aerospace Engineering & CFD

This solver models flow past a square cylinder using a pressure-projection method. The velocity field is advanced explicitly, the pressure Poisson equation is solved to enforce incompressibility, and the velocity is corrected so that the flow remains approximately divergence-free.

The convective terms use a second-order directional upwind stencil. The direction is chosen locally from the sign of the velocity, so the derivative uses information from the upstream side of the flow. Diffusion and pressure gradients are treated with second-order central differences. The time integration is an explicit forward-Euler update for the momentum predictor, with adaptive time stepping based on advective and viscous stability limits.

Python · NumPy · Matplotlib · Chorin projection · second-order upwind · central differences · adaptive time step

Grid spacing, numerical viscosity, and why the stencil had to improve

Before introducing the three-point upwind stencil, it is worth showing why the grid spacing matters so much. In the first-order version of the solver, the artificial diffusion from the upwind discretisation could dominate the flow field. On a coarse grid, the wake was heavily damped: instead of a clear cyclic vortex-shedding pattern, the solution looked more laminar and smeared out.

When the grid was refined, the wake structure became much more convincing. Periodic vortices became visible behind the square cylinder, showing that the physical instability was no longer being suppressed as strongly by the numerical method. This is the key point: decreasing dx reduces the artificial viscosity introduced by the scheme.

However, there is still a problem. Even if the input Reynolds number is set to 600, the calculation is not necessarily behaving like a true Re = 600 flow if the numerical viscosity is large. For a first-order upwind scheme, the leading artificial viscosity scales like

The effective Reynolds number is therefore closer to

At sufficiently large velocity, the numerical viscosity can become much larger than the physical viscosity. In that limit,

The velocity cancels. This means that increasing the inlet velocity no longer increases the effective Reynolds number in the way I intended. The scheme and grid spacing impose an upper limit. In one of these cases, the nominal Reynolds number may be 600, but the effective Reynolds number is closer to a much lower value, around 35, because the artificial viscosity is doing so much damping.

This is why I implemented a second-order upwind scheme using a three-point stencil for the advection term. The aim was to remove as much of the first-order numerical dissipation as possible. The trade-off is that lower dissipation can expose dispersive behaviour, but so far this has been controlled by using a sufficiently low CFL number.

Why divergence matters

A useful way to understand the incompressibility condition is through the volume of a moving fluid parcel. Consider a small moving control volume with fixed mass. Over a small time step, only the normal component of velocity through the surface changes the volume.

Using the divergence theorem, the surface flux can be written as a volume integral. If the small volume is shrunk down to a fluid element, this gives

So div(V) is the fractional rate of change of the volume of a fluid element. For an incompressible flow, the parcel volume should not expand or shrink, so

How I understand Chorin projection

The pressure term is what connects the momentum equation to the continuity constraint. If I advance the momentum equation without pressure, I get an intermediate velocity field, usually written as V*. This field contains the advection and diffusion effects, but it is not guaranteed to satisfy div(V*) = 0.

The pressure gradient is then used as a velocity correction. The factor dt/rho grad(p) has the units of velocity, so it makes sense to read the projection step as subtracting a pressure-generated velocity correction from the predicted field.

The corrected velocity must satisfy incompressibility:

Taking the divergence of the correction equation gives the pressure Poisson equation:

This is the core of the projection method. First compute the velocity that the momentum equation wants, then solve for the pressure field required to remove its divergence. The pressure is not just another plotted variable here. It acts like the constraint force that makes the velocity field compatible with incompressibility.

The directional upwind derivative uses a backward stencil when the local velocity is positive, and a forward stencil when the local velocity is negative:

This study models the 2D wave equation and monitors the total energy over time. The baseline energy is taken from the initial condition, and the moving energy curve is compared against that baseline. As the mesh is refined, the displacement between the initial baseline and the computed total energy reduces, which is the expected behaviour for a convergent scheme.

Small oscillations still appear when waves interact with boundaries or when separate wave packets meet and superpose. To investigate this, I used the modified-equation idea to identify leading numerical error behaviour, then tested a Lax-Wendroff-style correction designed to reduce the dominant numerical energy oscillations. FFTs of the energy error were then used to compare the frequency content of the error at different resolutions.

The continuous 2D wave equation has an energy-like quantity made from kinetic energy, through \(u_t\), and strain or gradient energy, through \(u_x\) and \(u_y\):

The energy result follows the same logic as the 1D derivation. Multiply the wave equation by \(u_t\), integrate over the domain, and integrate the spatial derivative terms by parts:

For boundary conditions where the boundary term vanishes, such as homogeneous Neumann \(\partial u/\partial n=0\) or fixed Dirichlet boundaries with no boundary motion, the energy is constant:

Energy comparison between central and modified Lax-Wendroff-style correction at N = 150.

FFT of the normalised energy error at N = 150.

Result observations

The energy comparison shows why the correction is useful. The unmodified central result has sharper energy-error excursions, especially around wave interactions and boundary reflections. The modified Lax-Wendroff-style correction reduces the largest excursions and keeps the energy curve closer to the initial baseline, although it does not remove all oscillatory behaviour.

The FFT plot shows that most of the energy-error content is concentrated at low frequencies. The correction changes the frequency content of the error and reduces the dominant energy-error response, which suggests that the remaining oscillation is not random noise but a structured numerical effect from the discretisation and boundary interactions.

2D wave equation · energy method · Lax-Wendroff correction · modified equation · FFT · mesh refinement

This study compares how different finite-difference schemes transport an initial signal. The aim is to see which schemes smear the solution through numerical dissipation and which schemes produce phase errors or oscillations through numerical dispersion.

One-dimensional solvers are a necessary building step. They make it easier to isolate the behaviour of a numerical scheme before using it in 2D or 3D. Many numerical methods are first tested on simple 1D cases because the errors are easier to interpret: diffusion shows up as smoothing, while dispersion shows up as phase error, overshoot, or oscillation.

Upwind schemes are usually robust but can be strongly diffusive. Lax-Friedrichs introduces additional smoothing. Lax-Wendroff is less diffusive, but can introduce dispersive oscillations near steep gradients. Comparing the same initial condition at several grid resolutions makes these behaviours easier to see.

Static comparison of upwind, Lax-Friedrichs, and Lax-Wendroff schemes for the 1D linear advection equation. The top row compares the transported solution against the analytical reference, while the bottom row shows the relative L2 error.

1D advection · upwind · Lax-Friedrichs · Lax-Wendroff · dispersion · dissipation · grid refinement