# Projection algorithm

###### Tags: `CFDLab` `modmesh`

[](https://www.youtube.com/embed/FO3IDwShGMA)

<iframe width="100%" height="350" src="https://www.youtube.com/embed/oaSOAeCJnuY" title="Lid Driven cavity-Projection algorithm" frameborder="0" allow="accelerometer; style="max-width: 100%; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

## Govern equations for 2D incompressible flow
For simplicity, suppose $\rho=1$.
+ Momentum equation in x-dir$$
\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}=-\frac{\partial p}{\partial x}+ \nu (\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2})
\tag{1}
$$
+ Momentum equation in y-dir$$
\frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}=-\frac{\partial p}{\partial y}+ \nu (\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2})
\tag{2}
$$
+ Continuity $$
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0
\tag{3}
$$

### Vector form
+ Momentum equation$$
\mathbb{u}_t+\mathbb{u}(\nabla \cdot \mathbb{u})=-\nabla p+\nu \nabla^2 \mathbb{u}
\tag{4}
$$
+ Continuity$$
\nabla \cdot \mathbb{u}=0
\tag{5}
$$

## Projection algorithm
Velocity at $n^{th}$ time-step is denoted as $\mathbb{u}^n$. Convection term is denoted as $\mathbb{C}$, and diffusion term is denoted $\mathbb{D}$.$$
\frac{\mathbb{u}^{n+1}-\mathbb{u}^{n}}{\Delta t}+\mathbb{C}(\mathbb{u}^{n})=-\nabla p+\mathbb{D}(\mathbb{u}^{{n}})
\tag{6}
$$
Solve for $\mathbb{u}^{n+1}$, we need to decide $\mathbb{P}$ beforehand. Projection algorithm is implemented with intermediate term $\mathbb{u}^{*}$, the momentum equation is split into 2 parts: $$
\frac{\mathbb{u}^{*}-\mathbb{u}^{n}}{\Delta t}+\mathbb{C}(\mathbb{u}^{n})=\mathbb{D}(\mathbb{u}^{{n}})
\tag{6a}
$$$$
\frac{\mathbb{u}^{n+1}-\mathbb{u}^{*}}{\Delta t}=-\nabla p
\tag{6b}
$$

take divergence, the equation becomes$$
\frac{\nabla \cdot(\mathbb{u}^{n+1}-\mathbb{u}^{*})}{\Delta t}=-\nabla ^2 p
\tag{7}
$$

A distinguishing feature of the projection method is that the velocity field is forced to satisfy the continuity as a part of the algorithm, where$$
\nabla \cdot(\mathbb{u}^{n+1})=0
\tag{8}
$$

and therefore the equation becomes, $$
\frac{\nabla \cdot(\mathbb{u}^{*})}{\Delta t}=\nabla ^2 p
\tag{9}
$$
##

Step 1- Solve $\mathbb{u}^{*}$$$
\frac{\mathbb{u}^{*}-\mathbb{u}^{n}}{\Delta t}+\mathbb{C}(\mathbb{u}^{n})=\mathbb{D}(\mathbb{u}^{{n}})
\tag{6a}
$$

Step 2- Solve $p$$$
\frac{\nabla \cdot(\mathbb{u}^{*})}{\Delta t}=\nabla ^2 p
\tag{9}
$$

Step 3- Solve $\mathbb{u}^{n+1}$$$
\frac{\mathbb{u}^{n+1}-\mathbb{u}^{*}}{\Delta t}=-\nabla p
\tag{6b}
$$

## Discretization
### Taylor expansion
$$
f(x+h)=f(x)+f'(x)\frac{h}{1!}+f''(x)\frac{h^2}{2!}+f'''(x)\frac{h^3}{3!}+... \tag{10}
$$ $$
f(x-h)=f(x)-f'(x)\frac{h}{1!}+f''(x)\frac{h^2}{2!}-f'''(x)\frac{h^3}{3!}+... \tag{11}
$$
+ First derivatives via (10)-(11)
$$
f(x+h)-f(x-h)=2\times f'(x)\frac{h}{1!}+2 \times f'''(x)\frac{h^3}{3!}+...
$$$$
f(x+h)-f(x-h)=2\times f'(x)\frac{h}{1!}+O(h^3)
$$$$
f'(x)=\frac{f(x+h)-f(x-h)}{2h}
\tag{12}
$$
+ Second derivative via (10)+(11)
$$
f(x+h)+f(x-h)=2\times f(x)+2 \times f''(x)\frac{h^2}{2!}+4 \times f''''(x)\frac{h^4}{4!}
$$$$
f(x+h)+f(x-h)=2\times f(x)+2 \times f''(x)\frac{h^2}{2!}+O(h^4)
$$$$
f''(x)=\frac{f(x+h)+f(x-h)-2\times f(x)}{h^2}
\tag{13}
$$
##
Step 1- Solve $\mathbb{u}^{*}$$$
\frac{\mathbb{u}^{*}-\mathbb{u}^{n}}{\Delta t}+\mathbb{C}(\mathbb{u}^{n})=\mathbb{D}(\mathbb{u}^{{n}})
\tag{6a}
$$
where$$
\mathbb{C}(u^n_{i,j})=u^n_{i,j} \cdot \frac{u^n_{e}-u^n_{w}}{\Delta x}+v^n_{p} \cdot \frac{u^n_{n}-u^n_{s}}{\Delta y}
$$$$
\mathbb{D}(u^n_{i,j})=\nu (\frac{u^n_{e}+ u^n_{w}-2u^n_{i,j}}{(\Delta x/2)^2}+\frac{u^n_{n}+ u^n_{s}-2u^n_{i,j}}{(\Delta y/2)^2})
$$

Step 2- Solve $p$$$
\frac{\nabla \cdot(\mathbb{u}^{*})}{\Delta t}=\nabla ^2 p
\tag{9}
$$
where$$
\nabla \cdot(\mathbb{u}^{*})=\frac{u^{*}_{e}-u^{*}_{w}}{\Delta x}+\frac{v^{*}_{n}-v^{*}_{s}}{\Delta y}
$$$$
\nabla ^2 p=\frac{p_{i+1,j}+p_{i-1,j}-2p_{i,j}}{\Delta x^2}+\frac{p_{i,j+1}+p_{i,j-1}-2p_{i,j}}{\Delta y^2}
$$

Step 3- Solve $\mathbb{u}^{n+1}$$$
\frac{\mathbb{u}^{n+1}_{i,j}-\mathbb{u}^{*}_{i,j}}{\Delta t}=-\nabla p
\tag{6b}
$$
where$$
\nabla p=\frac{p_{i+1,j}-p_{i,j}}{\Delta x}
$$

## Programming
Code can be found [here](https://github.com/EN-Chou/Projection_algorithm-Lid_driven_cavity).

Projection algorithm

int main(){
    while(velocity_not_converge()){
        //step 1
        cal_u_star();
        cal_v_star();
        //step 2
        cal_p();
        //step 3
        cal_u();
        cal_v();

        timestep++;
        moniter();
    }
    return 0;
}

Step 1- Solve $\mathbb{u}^{*}$$$
\frac{\mathbb{u}^{*}-\mathbb{u}^{n}}{\Delta t}+\mathbb{C}(\mathbb{u}^{n})=\mathbb{D}(\mathbb{u}^{{n}})
\tag{6a}
$$

C=u[i][j]*(u_e-u_w)/h+v_p*(u_n-u_s)/h;
D=1.0/Re*(u_e+u_w+u_n+u_s-4*u[i][j])/(h*h/4.0);
u_star[i][j]=(D-C)*delta_t+u[i][j];

C=u_p*(v_e-v_w)/h+v[i][j]*(v_n-v_s)/h;
D=1.0/Re*(v_e+v_w+v_n+v_s-4*v[i][j])/(h*h/4.0);
v_star[i][j]=(D-C)*delta_t+v[i][j];

Step 2- Solve $p$$$
\frac{\nabla \cdot(\mathbb{u}^{*})}{\Delta t}=\nabla ^2 p
\tag{9}
$$

// Iteratively compute until converges
poisson_left=(u_star_e-u_star_w+v_star_n-v_star_s)/h/delta_t;
p[i][j]=0.25*(p[i][j-1]+p[i][j+1]+p[i-1][j]+p[i+1][j]-poisson_left*h*h);

Step 3- Solve $\mathbb{u}^{n+1}$$$
\frac{\mathbb{u}^{n+1}-\mathbb{u}^{*}}{\Delta t}=-\nabla p
\tag{6b}
$$

u[i][j]=u_star[i][j]-delta_t*(p[i+1][j]-p[i][j])/h;

v[i][j]=v_star[i][j]-delta_t*(p[i][j+1]-p[i][j])/h;