2.3.4.1 Velocity Based

Pressure Poisson Equation¶

In this notebook, we provide examples of directly solving the Navier-Stokes equations in velocity space. To do this, we need to derive the pressure Poisson equation from the Navier-Stokes equations. Given 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} $$

we can take the divergence of the first equation and use the second equation to eliminate some terms:

• The divergence of the time derivative term:

$$ \nabla \cdot \frac{\partial \mathbf{u}}{\partial t} = \frac{\partial}{\partial t} (\nabla \cdot \mathbf{u}) = 0. $$

• The divergence of the pressure gradient term:

$$ \nabla \cdot \nabla p = \nabla^2 p. $$

• The divergence of the viscous term:

$$ \nabla \cdot \nu \nabla^2 \mathbf{u} = \nu \nabla^2 (\nabla \cdot \mathbf{u}) = 0. $$

The convection term and the external force term are not divergence-free, so we retain them. Thus, we have:

$$ \nabla \cdot \left[ (\mathbf{u} \cdot \nabla) \mathbf{u} - \mathbf{f} \right] = -\nabla^2 p. $$

This is the Poisson equation for the pressure. The solution of the Poisson equation gives us the pressure field, which we can use to compute the velocity field.

In TorchFSM, we directly integrate the Poisson equation into the momentum equation to obtain the velocity:

$$ \frac{\partial \mathbf{u}}{\partial t} = - \nabla \left( \nabla^{-2} \nabla \cdot \left[ (\mathbf{u} \cdot \nabla) \mathbf{u} - \mathbf{f} \right] \right) - (\mathbf{u} \cdot \nabla) \mathbf{u} + \nu \nabla^2 \mathbf{u} + \mathbf{f}. $$

The pressure is then solved by solving the Poisson equation in a postprocessing step. Since the convection term and the external force term are also involved in solving the pressure, we integrate these three terms together to avoid recomputing them. The NSPressureConvection is defined as:

$$ - \nabla \left( \nabla^{-2} \nabla \cdot \left[ (\mathbf{u} \cdot \nabla) \mathbf{u} - \mathbf{f} \right] \right) - (\mathbf{u} \cdot \nabla) \mathbf{u} + \mathbf{f}. $$

3D Taylor-Green Flow¶

In the following example, we demonstrate a 3D Taylor-Green vortex decay. The initial velocity field is given by a series of sine and cosine functions.

In [1]:

from torchfsm.operator import Laplacian,Operator,NSPressureConvection

def NavierStokes(Re:float)->Operator:
    return NSPressureConvection()+1/Re*Laplacian()

from torchfsm.operator import Laplacian,Operator,NSPressureConvection def NavierStokes(Re:float)->Operator: return NSPressureConvection()+1/Re*Laplacian()

In [2]:

from torchfsm.mesh import MeshGrid
from torchfsm.traj_recorder import CPURecorder,IntervalController
import numpy as np
import torch

mesh=MeshGrid([(0,2*np.pi,128),
               (0,2*np.pi,128),
               (0,2*np.pi,128)],
               device='cuda')
x,y,z=mesh.bc_mesh_grid()
u=torch.cat(
    [torch.sin(x)*torch.cos(y)*torch.cos(z),
     -torch.cos(x)*torch.sin(y)*torch.cos(z),
     torch.zeros_like(x)],
    dim=1)
ns=NavierStokes(1600)
traj=ns.integrate(
    u,
    dt=0.001,
    step=10000,
    mesh=mesh,
    trajectory_recorder=CPURecorder(control_func=IntervalController(200)),
    progressive=True
)

from torchfsm.mesh import MeshGrid from torchfsm.traj_recorder import CPURecorder,IntervalController import numpy as np import torch mesh=MeshGrid([(0,2*np.pi,128), (0,2*np.pi,128), (0,2*np.pi,128)], device='cuda') x,y,z=mesh.bc_mesh_grid() u=torch.cat( [torch.sin(x)*torch.cos(y)*torch.cos(z), -torch.cos(x)*torch.sin(y)*torch.cos(z), torch.zeros_like(x)], dim=1) ns=NavierStokes(1600) traj=ns.integrate( u, dt=0.001, step=10000, mesh=mesh, trajectory_recorder=CPURecorder(control_func=IntervalController(200)), progressive=True )

Integrating:   0%|          | 0/10000 [00:00<?, ?it/s]

In [3]:

from torchfsm.plot import plot_traj
plot_traj(traj)

from torchfsm.plot import plot_traj plot_traj(traj)

Out[3]:

The pressure filed can be obtained through postprocess:

In TorchFSM we use the ETDRK integrator by default. The ETDRK integrator can accurately solve the PDE with large time step while requiring high memory, which could be a server issue for 3D turbulence simulation. Thus, we also provide the classical Runge-Kutta integrator, which is more memory efficient but less accurate. For example, we can use the simple Euler integrator for the above simulation:

3D isotropic turbulence¶

To be updated...