Projection method FVM poisson part, adding source term

The idea of the method is to decompose the Navier-Stokes equation into the solenoidal and irrotational parts.

$$\frac{\partial u}{\partial t}+u(\nabla \cdot u)=-\frac{1}{\rho}\nabla p+\nabla ^2 u$$ Discretizing the transient term and initializing the intermediate velocity $$u^*$$.
$$\frac{(u^{n+1}-u^*) + (u^*-u^n)}{\Delta t}+u(\nabla \cdot u)=-\frac{1}{\rho}\nabla p+\nabla ^2 u$$ Decomposing the equation

1. $$\frac{(u^*-u^n)}{\Delta t}+u(\nabla \cdot u)=+\nabla ^2 u$$
2. $$\frac{(u^{n+1}-u^*)}{\Delta t}=-\frac{1}{\rho}\nabla p$$ expressing the $$u^{n+1}$$ from (2) and plugging into continuity equation leads to
3. $$\nabla^2p=\frac{\rho}{\Delta t}(\nabla \cdot u^{n+1})$$

Algorithm:

1. Solve equation (1) to find $$u^*$$
2. Solve equation (3) to find the $$p$$
3. Solve equation (2) to find $$u^{n+1}$$

I was struggling with the curved no-slip boundaries in step 1 of the algorithm described here. Keep the approach of nullifying the boundary cells for now.

The new issue arose when solving the Poisson equation from step 2 of the algorithm. The discretization of Poisson equation: $$\nabla\cdot(\nabla p)=\frac{\rho}{\Delta t}(\nabla \cdot u^{n+1})$$

Integrating over the finite volume $$\int\nabla\cdot(\nabla p)dV=\frac{\rho}{\Delta t}\int(\nabla \cdot u^{n+1})dV$$

Applying Gauss divergence theorem $$\int(\nabla p)\cdot \hat ndS=\frac{\rho}{\Delta t}\int(\hat n \cdot u^{n+1})dS$$

Decomposing the left part into normal and tangential flux:

$$\sum_{f=1}^{Nf,O}\bigg( \frac{p_{N(f)}-p_O}{\delta_f} - \big( \sum_{f}^{Nf,O} \frac{p_{a(f)} - p_{b(f)}}{\delta_f A_f}\big)\hat{t_f}\cdot{l_f}\bigg)A_f=\frac{\rho}{\Delta t}\int(\hat n \cdot u^{n+1})dS$$

$$\delta_f$$ is the normal distance between the cell and its neighbour through the face $$f$$
$$p_{a(f)},{p_{b(f)}}$$ are the nodal value of $$p$$ at the face $$f$$
$$A_f$$ is the length of face $$f$$
$$l_f$$ is the vector pointing from centre of the cell into its neighbours centre
$$\hat{t_f}$$ is the tangential unit vector

RHS discretization of the source term:

$$\sum_{f=1}^{Nf,O}\bigg( \frac{p_{N(f)}-p_O}{\delta_f} - \big( \sum_{f}^{Nf,O} \frac{p_{a(f)} - p_{b(f)}}{\delta_f A_f}\big)\hat{t_f}\cdot{l_f}\bigg)A_f=\frac{\rho}{\Delta t}\int(\hat n \cdot u^{*})dS$$ $$\sum_{f=1}^{Nf,O}\bigg( \frac{p_{N(f)}-p_O}{\delta_f} - \big( \sum_{f}^{Nf,O} \frac{p_{a(f)} - p_{b(f)}}{\delta_f A_f}\big)\hat{t_f}\cdot{l_f}\bigg)A_f=\frac{\rho}{\Delta t}\sum_{f=1}^{N_{f,O}} (u^*\cdot\hat n)A_f$$

Rearranging (considering triangular mesh):

$$\bigg( \frac{A_1}{\delta_1}+\frac{A_2}{\delta_2}+\frac{A_3}{\delta_3}\bigg)p_O - \frac{A_1}{\delta_1}p_1 - \frac{A_2}{\delta_2}p_2 - \frac{A_3}{\delta_3}p_3=$$ $$=-\frac{\rho}{\Delta t}\sum_{f=1}^{N_{f,O}} (u^*\cdot\hat n)A_f + \bigg( \frac{p_a-p_b}{\delta_1|t_1|}(\hat t_1 \cdot l_1)\bigg)A_1+\bigg( \frac{p_b-p_c}{\delta_2|t_2|}(\hat t_2 \cdot l_2)\bigg)A_2 + \bigg( \frac{p_c-p_a}{\delta_3|t_3|}(\hat t_3 \cdot l_3)\bigg)A_3$$

The RHS of the equation is named SCSKEW

Used Gauss Seidel algorithm

Running the code without the source term gives perfect results (left boundary 1, right Neumann condition = 0): Adding the source term (set $$u=(1,0)$$ at the left faces):

What could be wrong with the discretization or the implementation of it into the code?
As it is expected to slightly stretch the pressure field in the direction that it is pointing to, not move the pressure and increase it by 16 times from what is currently being obtained.
(the array is named source and is added to the scskew as they both are on the RHS of the discretized equation)? (running the code on refined mesh, playing with the values of $$\Delta t$$ did not affect the result)

Can paste or email the full code and input file if needed.

Update:

After application of Pressure boundary conditions below:
Left:$$p=1$$
Right:$$p=0$$
TopBot:$$\frac{\partial p}{\partial n}=0$$
The results seem to look nicer, however, diverge, issue might be in the vector field applied improperly

• Did you consider that you are solving pure Neumann problem for pressure Poisson equation? You should set the pressure at some point in the interior. Apr 3, 2021 at 21:09
• Was thinking of L:$p=1$ R:$p=0$ TB:$\frac{\partial p}{\partial n}=0$ could you kindly correct me if im wrong Apr 5, 2021 at 19:06
• @JohntraVolta thank you, for the suggestion, I tested all Neuman and multiple combinations of Dirichlet with Neuman BC's, code started working. Thank you for the idea! Apr 7, 2021 at 16:21

Walls and inlet: $$\frac{\partial p}{\partial n}=0$$
Outlet: $$p=0$$