2.3.4.1 Vorticity Based
Introduction to the Navier-Stokes Equations¶
The Navier-Stokes equations describe the motion of fluid. In this context, we focus on the incompressible form of the Navier-Stokes equations:
$$ \begin{align} \frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} &= -\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}, \\ \nabla \cdot \mathbf{u} &= 0, \end{align} $$where:
- $\mathbf{u}$ is the velocity field,
- $p$ is the pressure field,
- $\rho$ is the fluid density,
- $\nu$ is the kinematic viscosity, and
- $\mathbf{f}$ is the external force field.
The first equation is the momentum equation, and the second is the continuity equation.
Solving the Navier-Stokes equations is a challenging task due to the non-linearity of the equations. The non-linearity arises from the convective term $(\mathbf{u} \cdot \nabla)\mathbf{u}$, which makes the equations difficult to solve analytically. However, the equations can be solved numerically using various methods. Here, we introduce two common methods used in Fourier Spectral Methods:
Vorticity-Based Formulation¶
The vorticity-streamfunction formulation is a common method for solving the Navier-Stokes equations in two dimensions. The vorticity, $\boldsymbol{\omega}$, is defined as the curl of the velocity field:
$$ \boldsymbol{\omega} = \nabla \times \mathbf{u}. $$For two-dimensional, incompressible flow, the vorticity is a scalar quantity, and its transport equation is given by:
$$ \frac{\partial \omega}{\partial t} + (\mathbf{u} \cdot \nabla) \omega = \nu \nabla^2 \omega + \nabla \times \mathbf{f}. $$The derivation of the above vorticity transport equation can be found here.
Here, we provide two examples of solving the vorticity transport equation:
Decaying Flow¶
In this example, we consider a decaying flow in a periodic domain. There is no external force, and the initial flow is randomly generated. The flow decays over time due to viscosity.
from torchfsm.operator import Laplacian,Operator,VorticityConvection
from typing import Optional
def NavierStokesVorticity(Re:float=100,force:Optional[Operator]=None)->Operator:
ns_vorticity=-VorticityConvection() + 1/Re*Laplacian()
if force is not None:
ns_vorticity+=force
return ns_vorticity
import torch
from torchfsm.mesh import MeshGrid
from torchfsm.plot import plot_traj,plot_field
from torchfsm.field import kolm_force,diffused_noise
from torchfsm.traj_recorder import AutoRecorder,IntervalController
device='cuda' if torch.cuda.is_available() else 'cpu'
L=torch.pi*2; N=128;
mesh=MeshGrid([(0,L,N),(0,L,N)],device=device)
x,y=mesh.bc_mesh_grid()
decayed=NavierStokesVorticity(Re=1000)
u0=diffused_noise(mesh,device=device,n_batch=2,
zero_centered=False,
unit_magnitude=True,
unit_variance=False,
diffusion_coef=0.05)
traj=decayed.integrate(
u_0=u0,
mesh=mesh,
dt=0.1,
step=100*50,
trajectory_recorder=AutoRecorder(IntervalController(interval=50)),
progressive=True
)
Integrating: 0%| | 0/5000 [00:00<?, ?it/s]
plot_traj(traj)
Kolmogorov flow¶
Another good example is the Kolmogorov flow. Compared to previous decayed flow, kolmogorov flow is driven by a sinusoidal force and exhibit a more complex flow structure.
kolm=NavierStokesVorticity(Re=100,force=kolm_force(y),)
u0=diffused_noise(mesh,device=device,n_batch=2)
traj=kolm.integrate(
u_0=u0,
mesh=mesh,
dt=0.5/50,
step=100*50,
trajectory_recorder=AutoRecorder(IntervalController(interval=50))
)
Note that the force term in NavierStokesVorticity is actually the curl of the original force term in the Navier-Stokes equations, i.e., $\nabla \times \mathbf{f}$
plot_traj(traj)
Obtain velocity and pressure from vorticity¶
In the vorticity based formulation, we only solve the vorticity transport equation. To obtain the velocity, we can do some simple post-processing.
For a 2d velocity field $\mathbf{u}=(u,v)$, there exists a streamfunction $\psi$ such that: $$ \begin{align} u &= -\frac{\partial \psi}{\partial y},\\ v &= \frac{\partial \psi}{\partial x}. \end{align} $$
and the vorticity can be then expressed as
$$ \omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} = \nabla^2 \psi. $$Thus, we could first obtain the streamfunction $\psi$ from the vorticity $\omega$ by solving the Poisson equation and then obtain the velocity field $\mathbf{u}$ from the streamfunction $\psi$. Here, we directly provide an operator for all these post-processing steps.
from torchfsm.operator import Vorticity2Velocity
v2v=Vorticity2Velocity()
ori_shape=traj.shape
traj=traj.view(-1,*ori_shape[2:]) # move the time dimension to the batch dimension for parallelism
velocity=v2v(traj,mesh=mesh).view(*ori_shape[0:2],2,*ori_shape[3:])
plot_traj(velocity,channel_names=['$u_x$','$u_y$'])
Meanwhile, the pressure can also be solved through the Pressure Poisson equation, which we will discuss in the next section. But here we can show how we do it using our operator.
from torchfsm.operator import Vorticity2Pressure
v2p=Vorticity2Pressure(kolm_force(y))
pressure=v2p(traj,mesh=mesh).view(*ori_shape[0:2],1,*ori_shape[3:])
plot_traj(pressure,channel_names=['$p$'])