Stream Function
Theory and Physics
Definition of Stream Function
Professor, what is the meaning of the stream function? I hear the name often, but...
The stream function $\psi$ is a scalar function that, in two-dimensional incompressible flow, gives a velocity field that automatically satisfies the continuity equation. The definition is as follows.
From this definition, you can confirm that the continuity equation $\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0$ holds identically.
"Identically" means it holds exactly regardless of mesh coarseness or anything?
Analytically, yes. In numerical computation, discretization error enters, but if formulated with the stream function, mass conservation is structurally guaranteed. This is a major advantage of the stream function method over the primitive variable ($u, v, p$) method.
Relationship with Streamlines
What is the relationship between the stream function and streamlines?
The contour lines of the stream function coincide with streamlines. That is, fluid particles move along the curve $\psi = \text{const.}$ Furthermore, the volumetric flow rate per unit depth passing between two streamlines $\psi_1$ and $\psi_2$ is given by
This is the physical meaning of the stream function.
So the difference in stream function values directly gives the flow rate. That's convenient.
A wall is one of the streamlines, so $\psi = \text{const.}$ can be used as a boundary condition on the wall. For example, for a stationary wall, it is often set as $\psi = 0$.
Poisson Equation for Stream Function
It's used in combination with vorticity, right?
Substituting the definition of the stream function into the definition of vorticity $\omega$, $\omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}$, yields a Poisson equation.
In other words, if the vorticity distribution is known, solving this elliptic equation under boundary conditions determines the stream function (and thus the velocity field).
What happens in inviscid, irrotational flow ($\omega = 0$)?
It becomes the Laplace equation $\nabla^2 \psi = 0$. This is the stream function for potential flow and is dual to the Laplace equation for the velocity potential $\phi$. The contour lines of $\psi$ and $\phi$ are orthogonal.
Stokes Stream Function (Axisymmetric)
Can the stream function be used in 3D as well?
A stream function cannot be defined for general 3D flow. However, for axisymmetric flow, the Stokes stream function can be defined. In cylindrical coordinates $(r, z)$,
$$ u_z = \frac{1}{r} \frac{\partial \Psi}{\partial r}, \quad u_r = -\frac{1}{r} \frac{\partial \Psi}{\partial z} $$
Can the stream function be used in 3D as well?
A stream function cannot be defined for general 3D flow. However, for axisymmetric flow, the Stokes stream function can be defined. In cylindrical coordinates $(r, z)$,
Note that unlike the 2D $\psi$, the dimension of $\Psi$ becomes $[m^3/s]$. The volumetric flow rate between two stream surfaces is $Q = 2\pi(\Psi_2 - \Psi_1)$.
Seems useful for pipe flow, nozzle flow, etc.
Exactly right, it's very useful for analyzing axisymmetric jets and pipe flows. The stream function for Hagen-Poiseuille flow can be found analytically as $\Psi = \frac{U_0}{2R^2}(R^2 r^2 - r^4/2)$, so it can be used as a verification problem.
Scope of Application and Limitations
In what cases can the stream function method not be used?
The main constraints are as follows.
- General 3D flow: A scalar stream function cannot be defined. Extensions using vector potential are possible but not very practical.
- Compressible flow: Not applicable as it assumes incompressibility. A modified definition like $\rho u = \partial \psi / \partial y$ can be used for compressible flow, but it becomes complex.
- Multiphase flow: Continuity conditions for the stream function across interfaces become complex.
- 3D turbulence: Primitive variable methods (SIMPLE family, etc.) are mainstream in practical CFD.
So its main field of application is for education and research in 2D incompressible flow, right?
That understanding is correct. However, its computational robustness (structural guarantee of mass conservation) is attractive, so it is still often used for verifying 2D benchmark problems.
The Invention of the Stream Function—The Beautiful Description of 2D Flow Derived by Lagrange (1781)
The concept of the stream function ψ is said to have been introduced by the French mathematician Joseph-Louis Lagrange in a 1781 paper. It is an elegant mathematical device for two-dimensional incompressible flow, defining velocity components as u=∂ψ/∂y, v=-∂ψ/∂x, thereby automatically satisfying the continuity equation. Since the contour lines ψ=const. coincide with streamlines, it can be used directly for flow visualization. Combined with Bernoulli's theorem, it forms the core concept of potential flow theory (complex function of ψ and φ (velocity potential)) used for calculating lift on airfoils. In modern CFD, pressure-velocity coupling methods are mainstream, but plotting ψ contours remains one of the most intuitive visualization tools for deepening physical understanding of flow fields.
Physical Meaning of Each Term
- Temporal term $\partial(\rho\phi)/\partial t$: Imagine the moment you turn on a faucet. At first, water comes out in an unstable, splashing manner, but after a while, it becomes a steady flow, right? This "period of change" is described by the temporal term. The pulsation of blood flow due to heartbeats, or the flow fluctuation each time an engine valve opens and closes—all are unsteady phenomena. So what is steady-state analysis? Looking only at "after sufficient time has passed and the flow has settled down"—i.e., setting this term to zero. This significantly reduces computational cost, so solving first in steady-state is a basic CFD strategy.
- Convection term $\nabla \cdot (\rho \mathbf{u} \phi)$: What happens if you drop a leaf into a river? It gets carried downstream by the flow, right? This is "convection"—the effect where fluid motion transports things. Warm air from a heater reaching the far corner of a room is also because the "carrier," air, transports heat by convection. Here's the interesting part—this term contains "velocity × velocity," making it nonlinear. That is, as the flow becomes faster, this term rapidly strengthens, making control difficult. This is the root cause of turbulence. A common misconception: "Convection and conduction are similar" → They are completely different! Convection is carried by flow, conduction is transmitted by molecules. There's an order of magnitude difference in efficiency.
- Diffusion term $\nabla \cdot (\Gamma \nabla \phi)$: Have you ever put milk in coffee and left it? Even without stirring, after a while it naturally mixes, right? That's molecular diffusion. Now a question—honey and water, which flows more easily? Obviously water, right? Honey has high viscosity ($\mu$), so it flows poorly. When viscosity is large, the diffusion term becomes strong, and the fluid moves in a "thick" manner. In low Reynolds number flow (slow, viscous), diffusion is dominant. Conversely, in high Re number flow, convection overwhelms and diffusion plays a supporting role.
- Pressure term $-\nabla p$: When you push the plunger of a syringe, liquid shoots out forcefully from the needle tip, right? Why? Because the plunger side is high pressure, the needle tip is low pressure—this pressure difference provides the force that pushes the fluid. Dam discharge works on the same principle. On a weather map, where isobars are densely packed? That's right, strong winds blow. "Where there is a pressure difference, flow is generated"—this is the physical meaning of the pressure term in the Navier-Stokes equations. A point of confusion here: The "pressure" in CFD is often gauge pressure, not absolute pressure. When you switch to compressible analysis and suddenly get strange results, mixing up absolute/gauge pressure might be the cause.
- Source term $S_\phi$: Heated air rises—why? Because it becomes lighter (lower density) than its surroundings, so buoyancy pushes it upward. This buoyancy is added to the equation as a source term. Other examples: chemical reaction heat generated by a gas stove flame, Lorentz force acting on molten metal in a factory's electromagnetic pump... These are all actions that "inject energy or force into the fluid from the outside," expressed by the source term. What happens if you forget the source term? In natural convection analysis, if you forget to include buoyancy, the fluid doesn't move at all—a physically impossible result where warm air doesn't rise in a room with the heater on in winter.
Assumptions and Applicability Limits
- Continuum assumption: Valid for Knudsen number Kn < 0.01 (molecular mean free path ≪ characteristic length)
- Newtonian fluid assumption: Shear stress and strain rate have a linear relationship (non-Newtonian fluids require viscosity models)
- Incompressibility assumption (for Ma < 0.3): Treat density as constant. For Mach number 0.3 and above, compressibility effects must be considered.
- Boussinesq approximation (Natural convection): Consider density changes only in the buoyancy term, using constant density in other terms.
- Non-applicable cases: Rarefied gas (Kn > 0.1), supersonic/hypersonic flow (shock capturing required), free surface flow (requires VOF/Level Set, etc.)
Dimensional Analysis and Unit Systems
| Variable | SI Unit | Notes / Conversion Memo |
|---|---|---|
| Velocity $u$ | m/s | When converting from volumetric flow rate for inlet conditions, pay attention to cross-sectional area units. |
| Pressure $p$ | Pa | Distinguish between gauge pressure and absolute pressure. Use absolute pressure for compressible analysis. |
| Density $\rho$ | kg/m³ | Air: approx. 1.225 kg/m³@20°C, Water: approx. 998 kg/m³@20°C |
| Viscosity coefficient $\mu$ | Pa·s | Be careful not to confuse with kinematic viscosity coefficient $\nu = \mu/\rho$ [m²/s] |
| Reynolds number $Re$ | Dimensionless | $Re = \rho u L / \mu$. Indicator for laminar/turbulent transition. |
| CFL number | Dimensionless | $CFL = u \Delta t / \Delta x$. Directly related to time step stability. |
Numerical Methods and Implementation
Implementation of Vorticity-Stream Function Method
How do I actually implement the vorticity-stream function method in code?
The algorithm flow is as follows. Repeat the following each time step.
How do I actually implement the vorticity-stream function method in code?
The algorithm flow is as follows. Repeat the following each time step.
1. Time-integrate the vorticity equation to find $\omega^{n+1}$
2. Solve the Poisson equation $\nabla^2 \psi^{n+1} = -\omega^{n+1}$ to update the stream function
3. Compute the velocity field from $u = \partial\psi/\partial y$, $v = -\partial\psi/\partial x$
4. Update wall vorticity and return to step 1
Please teach me the specific discretization using finite difference.
On a staggered grid with uniform grid spacing $h$, the 5-point difference for the Poisson equation is
The central difference for velocity is
How is the advection term in the vorticity equation discretized?
With central difference, oscillations occur when the grid Reynolds number $Re_h = |u|h/\nu$ exceeds 2. Let's compare practical options.
| Scheme | Accuracy | Numerical Diffusion | Stability | Remarks |
|---|---|---|---|---|
| Central Difference | 2nd order | None | $Re_h < 2$ | Possibility of oscillations |
| 1st-order Upwind | 1st order | Large | Unconditional | Vortices disappear |
| QUICK | 3rd order | Small | $Re_h < 8/3$ | Good balance |
| Kawamura-Kuwahara | 3rd order | Very small | Good | Compact difference |
QUICK seems good.
Indeed, QUICK (Quadratic Upstream Interpolation for Convective Kinematics) offers a good balance for many problems. However, monotonicity is not guaranteed, so for cases with sharp discontinuities, combining it with a TVD (Total Variation Diminishing) limiter is advisable.
Iterative Solvers for Poisson Equation
Solving the Poisson equation every step becomes a cost bottleneck, right?
That's right. Let's compare the convergence speeds of typical iterative methods. For an $N \times N$ grid.
| Solver | Iteration Count ($O(\cdot)$) | Operations per Iteration | Total |
|---|---|---|---|
| Jacobi | $O(N^2)$ | $O(N^2)$ | $O(N^4)$ |
| Gauss-Seidel | $O(N^2)$ |