Numerical Solution to Navier-Stokes

Simple numerical solution to Navier-Stokes using FreeFEM++

The Problem

This is the final post in the series of simulation using Freefem++, where I sum up how to solve Navier-Stokes equations numerically. As mentioned in previous post, I have very limited background in physics, so I cannot fully understand why and how physicist came up with the Navier-Stokes equations for incompressible fluid. Therefore this blog is presented in a mathematical perspective (which basically mean that I agree with everything the physicists said and only do the math and the coding). Now let me formally restate the famous Navier-Stokes equations: \(\begin{align} \rho \left(\frac{\partial u}{\partial t} + (u . \nabla) u \right) - \mu \Delta u + \nabla v & = \rho f \tag{1} \\ \text{div} & = 0 \tag{2} \end{align}\)

where \(u(x, t)\) is velocity of the stream at position \(x\) (\(x \in \mathcal{R}^2\) or \(x \in \mathcal{R}^3\)) at time \(t\), and \(p(x, t)\) is the pressure of the fluid. The Navier-Stokes asks to find \(u\) and $v$ which satisfies (1) and (2) with given fluid property \(\rho\) (density) , \(\mu\) (viscosity) and external force applied on each volume unit \(f\).

Note that, since \(u\) is vector value function, (1) is actually short-handed for 2 equations in \(\mathcal{R}^2\) (or 3 equations in \(\mathcal{R}^3\)).

Numerical solution

The following program describe how water flow in a pipe (as shown in the figure) given that at one end the velocity of the water is kept at constant, there is some round-shaped obstacle within the pipe and the stream at the other end is allowed to move freely. The mysterious epsilon in the weak equations is require for the stability.

int n=100;
real T=1;
real dt=0.01;
int count=0;
border s1(t=5,20){x=t;y=0;label=1;};
border s2(t=0,3){x=20;y=t;label=2;};
border s3(t=20,0){x=t;y=3;label=1;};
border s4(t=3,1){x=0;y=t;label=3;};
border s5(t=0,5){x=t;y=1;label=1;};
border s6(t=1,0){x=5;y=t;label=1;};

border c(t=0, 2*pi){x = 10 + 0.25*cos(t); y = 1.5 + 0.25*sin(t);}
mesh Th= buildmesh (s1(n)+s2(n)+s3(n)+s4(n)+s5(n)+s6(n)+c(-n));

fespace Uh(Th,P1b); 
	Uh u1,u2,v1,v2,u1old=0,u2old=0,
fespace Ph(Th,P1);
	Ph p, q;

real epsilon=1e-10;
real mu=1e-5;

real alpha =1/dt;

solve stokes([u1,u2,p],[v1,v2,q])=
	+on(1,u1=0,u2=0) +on(3,u1=0.667,u2=0) + on(c, u1=0, u2=0);

for (int i=0;i<T/dt;i++){
	plot([u1,u2],p, ps="ns_" + count +".jpg");

Ph psi,phi;
solve streamlines(psi,phi) =
	int2d(Th)(dx(psi)*dx(phi) +dy(psi)*dy(phi))
plot(psi, ps="streamline.jpg");
Navier Stokes

Another way to view the flow is through streamlines. A streamline is a curve whose point has velocity vector tangent to it. Mathematically, streamlines can be define as a function \(\psi\) satifsies: \(u = \left(\frac{\psi}{y}, -\frac{\psi}{x} \right)^T\)

Navier Stokes Streamline