Implementation of a simulation of an incompressible Newtonian fluid with uniform density

Let $d\in\left\{2,3\right\}$ and $\Omega\subseteq\mathbb R^d$ be a bounded domain. I want to simulate an incompressible Newtonian fluid with uniform density $\rho$ and viscosity $\nu$.

The evolution up to time $T>0$ of such a fluid is given by the instationary Navier-Stokes equations $$\left\{\begin{matrix}\displaystyle\left(\frac\partial{\partial t}+\boldsymbol u\cdot\nabla\right)\boldsymbol u&=&\displaystyle\nu\Delta\boldsymbol u-\frac 1\rho\nabla p+\boldsymbol f&&\text{in }\Omega\times (0,T)\\\nabla\cdot \boldsymbol u&=&0&&\text{in }\Omega\times (0,T)\end{matrix}\right.\;,\tag 1$$ where $\boldsymbol u:\Omega\times [0,T]\to\mathbb R^d$ and $p:\Omega\times [0,T]\to\mathbb R$ are the velocity field and pressure we're looking for and $\boldsymbol f:\Omega\times (0,T)\to\mathbb R^d$ is the sum of all external forces.

Suppose we're given $\boldsymbol u_0=\boldsymbol u(\;\cdot\;,0)$ and $p_0=p(\;\cdot\;,0)$. How do we obtain $\boldsymbol u_\tau$ and $p_\tau$ for some delta time $\tau\in (0,T)$?

Let's restrict the question to $d=2$. Using the implicit Euler method, we need to solve the stationary problem $$\left\{\begin{matrix}\displaystyle\frac{\boldsymbol u_\tau-\boldsymbol u_0}\tau+\left(\boldsymbol u_\tau\cdot\nabla\right)\boldsymbol u_\tau&=&\displaystyle\nu\Delta\boldsymbol u_\tau-\frac 1\rho\nabla p_\tau+\boldsymbol f_\tau&&\text{in }\Omega\\\nabla\cdot \boldsymbol u_\tau&=&0&&\text{in }\Omega\end{matrix}\right.\;.\tag 2$$ This problem can be numerically solved using a Oseen iteration: $$\left\{\begin{matrix}\displaystyle\frac{\boldsymbol u_\tau^{(n)}-\boldsymbol u_0}\tau+\left(\boldsymbol u_\tau^{(n-1)}\cdot\nabla\right)\boldsymbol u_\tau^{(n)}&=&\displaystyle\nu\Delta\boldsymbol u_\tau^{(n)}-\frac 1\rho\nabla p_\tau^{(n)}+\boldsymbol f_\tau&&\text{in }\Omega\\\nabla\cdot \boldsymbol u_\tau^{(n)}&=&0&&\text{in }\Omega\end{matrix}\right.\;.\tag 3$$

Let's draw our attention to a simple example: Let

• $\boldsymbol u=(u,v)$ and $\boldsymbol f=(f,g)$
• $\Omega=(a,b)\times (c,d)$, $$x_i:=a+i\Delta x\;\;\;\text{for }i=0,\ldots,m:=\frac{b-a}{\Delta x}$$ and $$y_j:=c+j\Delta y\;\;\;\text{for }j=0,\ldots,n:=\frac{d-c}{\Delta y}$$
• $\varphi^{ij}:=\varphi(x_i,y_j)$ for any function $\varphi$ defined on $\Omega$

Using $$\frac{\partial\varphi}{\partial x}(x,y)\approx\frac{\varphi(x+h,y)-\varphi(x-h,y)}{2\Delta x}$$ and $$\frac{\partial^2\varphi}{\partial x^2}(x,y)\approx\frac{\varphi(x+h,y)-2\varphi(x,y)+\varphi(x-h,y)}{{\Delta x}^2}$$ (you can imagine how we approximate $\partial/\partial y$ and $\partial^2/\partial y^2$) we obtain a system

\begin{equation} \begin{split} \left(\frac{{\Delta x}^2{\Delta y}^2}\tau+4\nu\left({\Delta x}^2+{\Delta y}^2\right)\right)\color{red}{\tilde u_\tau^{ij}}&+{\Delta y}^2\left(\Delta xu_\tau^{ij}-2\nu\right)\color{red}{\tilde u_\tau^{(i+1)j}}\\&-{\Delta y}^2\left(\Delta x u_\tau^{ij}+2\nu\right)\color{red}{\tilde u_\tau^{(i-1)j}}\\&+{\Delta x}^2\left(\Delta y v_\tau^{ij}-2\nu\right)\color{red}{\tilde u_\tau^{i(j+1)}}\\&-{\Delta x}^2\left(\Delta y v_+^{ij}+2\nu\right)\color{red}{\tilde u_\tau^{i(j-1)}}\\&+\frac{\Delta x{\Delta y}^2}\rho\color{blue}{\tilde p_\tau^{(i+1)j}}\\&-\frac{\Delta x{\Delta y}^2}\rho\color{blue}{\tilde p_\tau^{(i-1)j}}=\frac{2{\Delta x}^2{\Delta y}^2}\tau u_0^{ij}+2{\Delta x}^2{\Delta y}^2f_\tau^{ij}\\ \left(\frac{{\Delta x}^2{\Delta y}^2}\tau+4\nu\left({\Delta x}^2+{\Delta y}^2\right)\right)\color{green}{\tilde v_\tau^{ij}}&+{\Delta y}^2\left(\Delta xu_\tau^{ij}-2\nu\right)\color{green}{\tilde v_\tau^{(i+1)j}}\\&-{\Delta y}^2\left(\Delta x u_\tau^{ij}+2\nu\right)\color{green}{\tilde v_\tau^{(i-1)j}}\\&+{\Delta x}^2\left(\Delta y v_\tau^{ij}-2\nu\right)\color{green}{\tilde v_\tau^{i(j+1)}}\\&-{\Delta x}^2\left(\Delta y v_+^{ij}+2\nu\right)\color{green}{\tilde v_\tau^{i(j-1)}}\\&+\frac{{\Delta x}^2\Delta y}\rho\color{blue}{\tilde p_\tau^{(i+1)j}}\\&-\frac{{\Delta x}^2\Delta y}\rho\color{blue}{\tilde p_\tau^{(i-1)j}}=\frac{2{\Delta x}^2{\Delta y}^2}\tau v_0^{ij}+2{\Delta x}^2{\Delta y}^2g_\tau^{ij}\\ \frac{1}{2\Delta x}\color{red}{\tilde u_\tau^{(i+1)j}}-\frac{1}{2\Delta x}\color{red}{\tilde u_\tau^{(i-1)j}}&=-\frac{1}{2\Delta y}\color{green}{\tilde v_\tau^{i(j+1)}}+\frac{1}{2\Delta y}\color{green}{\tilde v_\tau^{i(j-1)}} \end{split}\tag 4 \end{equation}

of $(d+1)(m-1)(n-1)$ linear equations (and just as many unknows), where I've written $\tilde\varphi$ instead of $\varphi^{(n)}$ and $\varphi$ instead of $\varphi^{(n-1)}$, for $\varphi\in\left\{u,v,p\right\}$, in order to beautify the notation.

If we introduce the more suitable indexing $$k\equiv i(n+1)+j\;\;\;\text{for }i=0,\ldots,m\text{ and }j=0,\ldots,n\;,$$ we're able to write $(4)$ in the form $$Ax=b\;.\tag 5$$

However, I'm unsure whether or not my attempt so far is sensible or not. Moreover, $A$ in $(5)$ is sparse and somehow regular, but it's still not easy to implement a solver.

To be explicit: Given $d=2$ and $\Omega=(a,b)\times (c,d)$, I want to implement a simple solver for $(5)$. How can I do so? How should I generate $A$ and $b$? Should I use a Gauss–Seidel iteration to solve $(5)$? And how would such an iteration look like algorithmically?

• If it's not clear enough: The "unknows" are the values of $u_\tau$, $v_\tau$ and $p_\tau$ at the "inner" grid points. Moreover, you can assume that the values of these quantities at the boundary grid points are given. – 0xbadf00d Dec 10 '15 at 16:44
• What you've done is just a standard first-order implicit central-differencing scheme. What language are you going to write your code in? In Matlab, you can just use A\b to solve the matrix problem. The matrix M isn't too hard to generate, diag will be useful, just look at some examples. You need to include your boundary conditions in the matrix problem as well. – David Dec 13 '15 at 22:59
• @David I need to write the code in C++ and I don't want to use an external library. – 0xbadf00d Dec 14 '15 at 0:19
• I don't really know anything about C++, but I'm suer you'll need some kind of linear algebra/math libraries. If you can write the matrix problem you need to solve down on paper, it should be easy enough to code up. Maybe you should start with a simple PDE rather than the Navier-Stokes equations for your first time doing numerics? – David Dec 14 '15 at 0:38
• @David This is not my first time doing numerics ;) And I don't want to solve the Navier-Stokes equations cause they are interesting PDEs, but because I need to simulate a fluid. I wasn't sure, if my approach is sensible or not. (and I'm aware of the literature in CFD. I just want to look what I can achieve with simple considerations). – 0xbadf00d Dec 14 '15 at 11:53