Upwind Scheme#

The Upwind Scheme is a finite difference method specifically designed to handle advection-dominated problems more effectively than symmetric schemes like FTCS. By incorporating the direction of wave propagation into the discretization of spatial derivatives, the upwind method enhances numerical stability and reduces non-physical oscillations.

In advection processes, information propagates in a specific direction determined by the flow velocity $c$. The upwind scheme leverages this directional information to bias the spatial derivative approximation, ensuring that the numerical flux aligns with the physical transport direction. This directional bias significantly improves the stability of the numerical solution, especially when dealing with sharp gradients or discontinuities.

The upwind scheme discretizes the spatial derivative based on the sign of the advection speed $c$:

• For $c>0$ (flow to the right):

  (528)#$$\begin{array}{r}\frac{\partial u}{\partial x}\approx \frac{u_i^n-u_{i-1}^n}{\mathrm{\Delta }x}\end{array}$$

• For $c<0$ (flow to the left):

  (529)#$$\begin{array}{r}\frac{\partial u}{\partial x}\approx \frac{u_{i+1}^n-u_i^n}{\mathrm{\Delta }x}\end{array}$$

Assuming $c>0$ for this implementation, the Upwind Scheme update rule becomes:

where:

• $u_i^n$ is the numerical approximation of $u$ at spatial index $i$ and time level $n$,

• $\mathrm{\Delta }t$ and $\mathrm{\Delta }x$ are the time and spatial step sizes, respectively.

The following Python code implements the upwind scheme to solve the linear advection equation for a sinusoidal initial condition.

import warnings

# Parameters
c = 1.0          # Advection speed
L = 1.0          # Domain length
T = 1.25         # Total time
nx = 100         # Number of spatial points
dx = L / nx      # Spatial step size
dt = 0.001       # Initial time step size
nt = int(T / dt) # Number of time steps

# Stability parameter
sigma = c * dt / dx
print(f"Courant number (sigma): {sigma}")

# Check CFL condition
if sigma > 1:
    warnings.warn(f"CFL condition violated: sigma = {sigma} > 1. Please reduce dt or increase dx.")

# Spatial grid
x = np.linspace(0, L, nx, endpoint=False)
u_initial = np.sin(2 * np.pi * x)  # Initial condition: sinusoidal wave

# Initialize solution array
u = u_initial.copy()

# Time-stepping loop using Upwind scheme
for n in range(nt):
    # Apply periodic boundary conditions using np.roll
    u_new = u - sigma * (u - np.roll(u, 1))
    u = u_new

# Analytical solution
u_exact = np.sin(2 * np.pi * (x - c * T))

# Plotting the results
plt.figure(figsize=(12, 6))
plt.plot(x, u_initial, label='Initial Condition', linestyle='--')
plt.plot(x, u_exact,   label='Exact Solution', linewidth=2)
plt.plot(x, u,         label='Upwind Scheme', linestyle=':', linewidth=2)
plt.xlabel('x')
plt.ylabel('u')
plt.legend()
plt.grid(True)

Courant number (sigma): 0.1

Von Neumann Stability Analysis#

Von Neumann Stability Analysis assesses the stability of numerical schemes by examining how Fourier modes are amplified or damped during the simulation.

Assume the numerical solution can be expressed as:

where:

• $G$ is the amplification factor,

• $k$ is the wave number,

• $x_i=i\mathrm{\Delta }x$.

Substitute $u_i^n=G^n{e}^{ik{x}_i}$ into the upwind update equation:

Divide both sides by $G^n{e}^{ik{x}_i}$:

Using Euler’s formula and separate the real and imaginary parts:

For stability, the magnitude of the amplification factor must satisfy $|G|\le 1$.

Calculate $|G{|}^2$:

Using ${\cos }^2(\theta )+{\sin }^2(\theta )=1$:

To find the maximum $|G|$, consider $\cos (k\mathrm{\Delta }x)=-1$ so that

For stability, we require $|G|\le 1$. This yields two inequalities $0\le \sigma \le 1$. This implies that the Upwind Scheme is stable provided that the Courant number satisfies:

Modified Equation Analysis#

While Von Neumann Stability Analysis provides insights into the stability of numerical schemes, Modified Equation Analysis delves deeper by examining the leading-order truncation errors introduced by discretization. This analysis reveals how numerical schemes can inadvertently introduce artificial diffusion or dispersion into the solution.

To determine how the Upwind Scheme modifies the original advection equation by introducing additional terms that represent numerical errors.

1. Express the Upwind Scheme Using Taylor Series Expansions

   Start with the upwind update equation:

   (539)#$$\begin{array}{r}{u}_i^{n+1}=u_i^n-\sigma (u_i^n-u_{i-1}^n),\end{array}$$

   in addition to the grid point that we want to update, $u_i^n$, there are two additional grid points. They are $u_i^{n+1}$ and $u_{i-1}^n$.

   Expand $u_i^{n+1}$ and $u_{i-1}^n$ around the point $(x_i,t^n)$ using Taylor series:

   (540)#$$\begin{array}{rl}{u}_i^{n+1}& =u(x_i,t^n)+\mathrm{\Delta }t{\frac{\partial u}{\partial t}|}_{x=x_i,t=t_n}+\frac{\mathrm{\Delta }{t}^2}{2}{\frac{{\partial }^{2}u}{\partial t^2}|}_{x=x_i,t=t_n}+\mathcal{O}(\mathrm{\Delta }{t}^3)\\ u_{i-1}^n& =u(x_i,t^n)-\mathrm{\Delta }x{\frac{\partial u}{\partial x}|}_{x=x_i,t=t_n}+\frac{\mathrm{\Delta }{x}^2}{2}{\frac{{\partial }^{2}u}{\partial x^2}|}_{x=x_i,t=t_n}+\mathcal{O}(\mathrm{\Delta }{x}^3)\end{array}$$

   Using the shorthands $\dot{u}\equiv \partial u/\partial t$ and $u^{\prime }\equiv \partial u/partialx$, the above series can be written as

   (541)#$$\begin{array}{rl}{u}_i^{n+1}& =u_i^n+\mathrm{\Delta }t{\dot{u}}_i^n+\frac{\mathrm{\Delta }{t}^2}{2}{\ddot{u}}_i^n+\mathcal{O}(\mathrm{\Delta }{t}^3)\\ u_{i-1}^n& =u_i^n-\mathrm{\Delta }x{u^{\prime }}_i^n+\frac{\mathrm{\Delta }{x}^2}{2}{u^{″}}_i^n+\mathcal{O}(\mathrm{\Delta }{x}^3)\end{array}$$

2. Substitute the Taylor Expansions into the Upwind Scheme

   Substitute the expansions into the upwind update equation:

   (542)#$$\begin{array}{r}{u}_i^n+\mathrm{\Delta }t{\dot{u}}_i^n+\frac{\mathrm{\Delta }{t}^2}{2}{\ddot{u}}_i^n=u_i^n-\sigma [u_i^n-(u_i^n-\mathrm{\Delta }x{u^{\prime }}_i^n+\frac{\mathrm{\Delta }{x}^2}{2}{u^{″}}_i^n)]\end{array}$$

   Because the derivatives are exact (instead of numerical), we can drop the functions’ evaluations at the dicrete points $x_i$ and $t_n$. Simplify the equation, we have:

   (543)#$$\begin{array}{r}\mathrm{\Delta }t\frac{\partial u}{\partial t}+\frac{\mathrm{\Delta }{t}^2}{2}\frac{{\partial }^{2}u}{\partial t^2}=-\sigma (\mathrm{\Delta }x\frac{\partial u}{\partial x}-\frac{\mathrm{\Delta }{x}^2}{2}\frac{{\partial }^{2}u}{\partial x^2})\end{array}$$

3. Rearrange and Substitute the Original PDE

   Taking the spatial and temporal derivatives of the original advection equation $\partial u/\partial t=-c\partial u/\partial x$, we have

   (544)#$$\begin{array}{r}\frac{{\partial }^{2}u}{\partial t^2}=-c\frac{{\partial }^{2}u}{\partial x\partial t},\frac{{\partial }^{2}u}{\partial t\partial x}=-c\frac{{\partial }^{2}u}{\partial x^2}\end{array}$$

   Combine, we obtain the two-way wave equation:

   (545)#$$\begin{array}{r}\frac{{\partial }^{2}u}{\partial t^2}=c^2\frac{{\partial }^{2}u}{\partial x^2}.\end{array}$$

   Substituting the wave equation into the modified equation, we obtain:

   (546)#$$\begin{array}{r}\mathrm{\Delta }t\frac{\partial u}{\partial t}+\frac{c^2\mathrm{\Delta }{t}^2}{2}\frac{{\partial }^{2}u}{\partial x^2}=-\sigma (\mathrm{\Delta }x\frac{\partial u}{\partial x}-\frac{\mathrm{\Delta }{x}^2}{2}\frac{{\partial }^{2}u}{\partial x^2}).\end{array}$$

4. Combine Like Terms

   Put back the definition of $\sigma$ and rearrange,

   (547)#$$\begin{array}{rl}\frac{\partial u}{\partial t}+\frac{c^2\mathrm{\Delta }t}{2}\frac{{\partial }^{2}u}{\partial x^2}& =-c(\frac{\partial u}{\partial x}-\frac{\mathrm{\Delta }x}{2}\frac{{\partial }^{2}u}{\partial x^2})\\ \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}& =\frac{c}{2}(\mathrm{\Delta }x-c\mathrm{\Delta }t)\frac{{\partial }^{2}u}{\partial x^2}\end{array}$$

   Define ${\nu }_{\text{upwind}}\equiv (\mathrm{\Delta }x-c\mathrm{\Delta }t)c/2$ be the numerical diffusion coefficient, the above equation reduces to

   (548)#$$\begin{array}{r}\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}={\nu }_{\text{upwind}}\frac{{\partial }^{2}u}{\partial x^2}.\end{array}$$

The Modified Equation Analysis reveals that the Upwind Scheme introduces an artificial diffusion term proportional to $\mathrm{\Delta }x-c\mathrm{\Delta }t$. This numerical diffusion acts to smooth out sharp gradients in the solution, which can be beneficial in suppressing non-physical oscillations that may arise from discretization errors. By damping these oscillations, the scheme enhances the overall stability of the numerical solution, ensuring that the simulation remains free from erratic and unrealistic fluctuations.

However, this introduction of numerical diffusion comes with a trade-off in accuracy. While the diffusion helps stabilize the solution, it also has the unintended effect of smearing out important features such as sharp interfaces or discontinuities within the advected quantity $u$. This attenuation can lead to a loss of fidelity in capturing precise wavefronts or shock formations, which are critical in accurately modeling phenomena like shock waves in fluid dynamics. Consequently, while the Upwind Scheme mitigates instability, it may compromise the sharpness and detail of the solution where high accuracy is required.

Furthermore, the Modified Equation Analysis reinforces the critical importance of adhering to the Courant-Friedrichs-Lewy (CFL) condition $\sigma \le 1$. Satisfying this condition is not only essential for maintaining the stability of the numerical scheme but also plays a pivotal role in controlling the extent of artificial diffusion introduced by the Upwind Scheme. When the CFL condition is violated ($\sigma >1$), not only does the numerical scheme become unstable, leading to divergent and oscillatory solutions, but the altered modified equation exacerbates the artificial diffusion. This excessive diffusion further degrades the solution quality, making it increasingly inaccurate and unreliable. Therefore, ensuring that the CFL condition is satisfied is paramount in balancing stability and accuracy, enabling the Upwind Scheme to perform optimally without introducing significant numerical artifacts.

Derivation#

We begin by expanding $u(x,t+\mathrm{\Delta }t)$ in time around $t$ using a Taylor series up to second order:

Substitute the advection equation and the wave equation into the Taylor series expansion:

and approximate the spatial derivatives using centered finite differences:

we obtain:

Simplify the equation to isolate $u_i^{n+1}$, we obtain the Lax-Wendroff Scheme for the linear advection equation

# Parameters
c = 1.0          # Advection speed
L = 1.0          # Domain length
T = 1.25         # Total time
nx = 100         # Number of spatial points
dx = L / nx      # Spatial step size
dt = 0.005       # Time step size
nt = int(T / dt) # Number of time steps

# Courant number
sigma = c * dt / dx
print(f"Courant number (sigma): {sigma}")

# Stability condition (for Lax-Wendroff, no strict CFL condition, but accuracy improves with sigma <=1)
if sigma > 1:
    warnings.warn(f"Courant number sigma = {sigma} > 1. Accuracy may decrease.")

# Spatial grid
x = np.linspace(0, L, nx, endpoint=False)
u_initial = np.sin(2 * np.pi * x)  # Initial condition: sinusoidal wave

# Initialize solution array
u = u_initial.copy()

# Time-stepping loop using Lax-Wendroff scheme
for n in range(nt):
    u_new = np.zeros_like(u)
    for i in range(nx):
        u_new[i] = (u[i]
                    - 0.5 * sigma    * (u[(i+1)%nx] -          u[(i-1)%nx])
                    + 0.5 * sigma**2 * (u[(i+1)%nx] - 2*u[i] + u[(i-1)%nx]))
    u = u_new

# Analytical solution
u_exact = np.sin(2 * np.pi * (x - c * T))

# Plotting the results
plt.figure(figsize=(12, 6))
plt.plot(x, u_initial, label='Initial Condition', linestyle='--')
plt.plot(x, u_exact,   label='Exact Solution', linewidth=2)
plt.plot(x, u,         label='Lax-Wendroff Scheme', linestyle=':', linewidth=2)
plt.xlabel('x')
plt.ylabel('u')
plt.legend()
plt.grid(True)

Courant number (sigma): 0.5

Purpose of Non-Dimensionalization#

Non-dimensionalization serves several critical purposes in the analysis and solution of PDEs:

1. Simplifying Equations for Analysis

   By scaling variables and parameters, non-dimensionalization reduces the number of independent variables and parameters in the governing equations. This simplification often transforms complex, dimensional equations into a more manageable, dimensionless form. For instance, consider the non-dimensionalization of the Navier-Stokes equations:

   (562)#$$\begin{array}{r}\frac{\partial \overline{\mathbf{u}}}{\partial \overline{t}}+(\overline{\mathbf{u}}\cdot \overline{\mathrm{\nabla }})\overline{\mathbf{u}}=-\overline{\mathrm{\nabla }}\overline{p}+\frac{1}{\text{Re}}{\overline{\mathrm{\nabla }}}^2\overline{\mathbf{u}}+\overline{\mathbf{f}}\end{array}$$

   Here, the introduction of dimensionless variables $\overline{\mathbf{u}}$, $\overline{t}$, and $\overline{p}$ has reduced the complexity of the original equations by consolidating physical parameters into dimensionless groups like the Reynolds Number (Re).

2. Identifying Dominant Physical Effects in Specific Regimes

   Non-dimensionalization reveals the relative importance of various physical phenomena within a system. By examining the dimensionless numbers that emerge from the process, one can determine which terms in the equations are significant and which can be neglected under certain conditions. This identification is crucial for developing approximate solutions and understanding the behavior of the system in different regimes.

   For example, in high Reynolds number flows (Re $\gg 1$), inertial forces dominate over viscous forces, simplifying the Navier-Stokes equations by reducing the influence of the viscous term. Conversely, in low Reynolds number flows (Re $\ll 1$), viscous forces are predominant, and inertial terms can often be neglected.

Key Dimensionless Numbers#

Several dimensionless numbers play pivotal roles in fluid dynamics and the study of PDEs. Among the most important are the Reynolds Number (Re), Mach Number (Ma), and Prandtl Number (Pr). Each of these numbers encapsulates the ratio of different physical effects, providing insight into the system’s behavior.

Applications of Dimensionless Numbers#

Dimensionless numbers like Re, Ma, and Pr are instrumental in analyzing and predicting the behavior of physical systems across various domains:

1. Flow Control and Design

   In engineering, these numbers guide the design of systems involving fluid flow. For example, ensuring that aircraft operate efficiently at desired Reynolds and Mach numbers is crucial for performance and safety.

2. Scale Modeling and Similitude

In experimental studies, ensuring that the dimensionless numbers of a model match those of the real system (dynamic similarity) allows for accurate predictions and validations of theoretical models through physical experiments.

3. Astrophysical Fluid Systems

   In astrophysics, dimensionless numbers help in scaling and understanding fluid behaviors in vastly different environments:

   • Stellar Convection: High Reynolds numbers indicate turbulent convection currents within stars, affecting energy transport and mixing of stellar material.

   • Accretion Disks: Mach numbers determine the compressibility of gas in accretion disks around black holes, influencing the formation of shock waves and jet structures.

   • Interstellar Medium: Prandtl numbers aid in modeling the thermal and momentum diffusion in the diffuse interstellar gas, impacting star formation rates and the dynamics of molecular clouds.

Historical Anecdote: Fermi’s Estimation of Atomic Bomb Energy Release#

One of the most illustrative examples of the power of non-dimensionalization and dimensional analysis in physics is Enrico Fermi’s remarkable estimation of the energy released during the first atomic bomb test, known as the Trinity test, conducted on July 16, 1945. This anecdote not only highlights Fermi’s ingenuity but also demostrate the practical utility of dimensionless numbers and scaling laws in tackling complex, real-world problems with limited empirical data.

During World War II, the Manhattan Project was a top-secret initiative aimed at developing nuclear weapons. As the project progressed, scientists faced the daunting challenge of predicting the yield of the atomic bombs they were designing. Precise calculations were hindered by the lack of comprehensive empirical data, necessitating innovative approaches to estimate explosive power accurately.

Enrico Fermi, a renowned physicist known for his ability to make quick, accurate estimates with minimal data, was tasked with predicting the energy release of the impending Trinity test. His method exemplified dimensional analysis—a technique that uses the fundamental dimensions (such as mass, length, and time) to derive relationships between physical quantities.

Fermi’s estimation process involved the following steps:

1. Observing the Shock Wave: After the detonation, Fermi and his colleagues observed the speed at which the shock wave propagated through the air. Let’s denote this speed as $U$.

2. Estimating the Shock Wave Radius: By timing how long it took for the shock wave to reach a certain distance from the explosion site, Fermi estimated the radius $L$ of the blast wave after a specific duration $t$.

3. Applying the Sedov-Taylor Scaling: The Sedov-Taylor solution describes the propagation of a strong shock wave from a point explosion in a homogeneous medium. According to this theory, the energy $E$ of the explosion is related to the shock wave’s speed $U$, radius $L$, and the density $\rho$ of the surrounding air by the following relationship:

   (566)#$$\begin{array}{r}E\sim \rho U^2{L}^3\end{array}$$

   This equation is derived by balancing the kinetic energy imparted to the air by the shock wave with the energy of the explosion.

4. Calculating the Energy: By substituting the estimated values of $\rho$, $U$, and $L$ into the above equation, Fermi arrived at an approximate value for the energy $E$ released during the explosion.

There is an arXiv paper on this topic.