# 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

The Gauss Seidel algorithm:

res = 1;
epsilon = 1e-5;
ctr = 0;

phi = zeros(ncells,1);

% LOOP HERRE
while res > epsilon
for ifc = 1:nfaces
c1 = link_face_to_cell(ifc,1);
c2 = link_face_to_cell(ifc,2);
if c2 ~= 0
uface(ifc,:) = wf(ifc) * ucell(c1,:) + (1-wf(ifc)) * ucell(c2,:);
end
end
% end calculating inner face values

calculate boundary face values in U DIRERCTION
for ifc = 1:nbfaces
if ismember(ifc,botf)
uface(ifc,1) = 0;
elseif ismember(ifc,topf)
uface(ifc,1) = 0;
elseif ismember(ifc,leftf)
uface(ifc,1) = 1;
elseif ismember(ifc,rightf)
ic = link_face_to_cell(ifc,1);
uface(ifc,1) = ucell(ic,1);
end
end
ustar = ucell;
res = 1;
epsilonjac = epsilon;
ctrpsn = 0;
phiv = zeros(nnodes,1);

phib = zeros(nbfaces,1); % ??? ADDED FOR STH? DUNNO WHY !?!?!?!!!!!!!!!!!!!!??!!??!?!?!?!?!?!?!?!?
for i = 1:length(leftf)
phib(leftf(i)) = 1;
end

ap = zeros(ncells,1);
sc = zeros(ncells,1);
anb = zeros(ncells,3);
for ic = 1:ncells
ap(ic) = 0;
sc(ic) = 0;
for j = 1:3
ifc = link_cell_to_face(ic,j);
if link_face_to_bface(ifc) == 0 % if interior
ap(ic) = ap(ic) + areaf(ifc)/deltaf(ifc); % Ao,O
anb(ic,j) = - areaf(ifc)/deltaf(ifc); % Aj,O
elseif  ismember(ifc,rightf) || ismember(ifc,leftf)% boundary
ifb = link_face_to_bface(ifc);
ap(ic) = ap(ic) + areaf(ifc)/deltaf(ifc); % dirichlet. Ao,O THIS IS Anb at boundary
anb(ic,j) = 0;
sc(ic) = sc(ic) + phib(ifb) * areaf(ifc) / deltaf(ifc); % THIS IS BOUNNDARY VALUE RHS OF MATMUL=X
elseif ismember(ifc,topf) || ismember(ifc,botf)
ifb = link_face_to_bface(ifc);
ap(ic) = ap(ic) + 0; % no dirichlet, nullifying. Ao,O THIS IS Anb at boundary
anb(ic,j) = 0;
sc(ic) = sc(ic) + 0; % dp/dn = 0  RHS FLUX
end
end

end

% LOOP HERRE
while res > epsilonjac
% boundaries NEUMANN
for ifc = 1:nbfaces
ic = link_face_to_cell(ifc,1);
if ismember(ifc,botf)
phib(ifc) = phi(ic);
elseif ismember(ifc,topf)
phib(ifc) = phi(ic);
elseif ismember(ifc,leftf)
phib(ifc) = 1;
elseif ismember(ifc,rightf)
phib(ifc) = 0;
end
end

%compute vertex values
phinode = zeros(nnodes,1);

for iv = 1:nnodes
phinode(iv) = 0;
weight(iv) = 0;
for ic = 1:sum(wv(iv,:)~=0) %over nonzero wv elements
if bnode(iv) == 0 % FOR INTERIOR
phinode(iv) = phinode(iv) + phi(link_node_to_cell(iv,ic)) * wv(iv,ic);
else % FOR BOUNDARY
phinode(iv) = phinode(iv) + phi(link_node_to_cell(iv,ic)) * wv(iv,ic);
end

end
end

%%%%%%%%%%%%%%%%%%%%%%% HARDWIRE BOUNDARIES AT THE LEFT
for ifc = 1:nbfaces
if ismember(ifc,leftf)
phinode(link_face_to_node(ifc,1)) = 1;
phinode(link_face_to_node(ifc,2)) = 1;
elseif ismember(ifc,rightf)
phinode(link_face_to_node(ifc,1)) = 0;
phinode(link_face_to_node(ifc,2)) = 0;
end
end % check with find(phinode) of initial

scskew = zeros(ncells,1); % TANGENTIAL FLUX SOURCE
for ic = 1:ncells
scskew(ic) = 0;
sumf = 0;
sumsource = 0;
for j = 1:3 % over faces of each cell
ifc = link_cell_to_face(ic,j);
if link_face_to_bface(ifc) == 0 % skip boundary faces
c1 = link_face_to_cell(ifc,1);
c2 = link_face_to_cell(ifc,2);
v1 = link_face_to_node(ifc,1);
v2 = link_face_to_node(ifc,2);
dxl = xc(c2) - xc(c1);
dyl = yc(c2) - yc(c1);
tdotl = st(ifc,1) * dxl + st(ifc,2) * dyl;
sumf = sumf + tdotl * (phinode(v2) - phinode(v1)) * snsign(ic,j)/deltaf(ifc); % check order of subtraction
sumf = sumf; %  - uface(ifc,:) * sn(ifc,:)' * snsign(ic,j) * (rho/dt) * areaf(ifc)
sumf = sumf -  uface(ifc,:) * sn(ifc,:)' * snsign(ic,j) * (rho/dt) * areaf(ifc);
end
%%% - uface(ifc,:) * sn(ifc,:)' * snsign(ic,j) * (rho/dt) * areaf(ifc)

end
scskew(ic) = sumf;

end

for ic = 1:ncells % DIVERGENCE ON THE RIGHT AS PREDICTED?
sumS = 0;
for j = 1:3
ifc = link_cell_to_face(ic,j);
sumS = sumS + uface(ifc,:) * sn(ifc,:)' * snsign(ic,j) * dt/rho * areaf(ifc);
end
scskew(ic) = scskew(ic)  - sumS;
end

for ic = 1:ncells
sumf = 0;
for j = 1:3 % over neighbours
ifc = link_cell_to_face(ic,j);
icn1 = link_face_to_cell(ifc,1);
icn2 = link_face_to_cell(ifc,2);
if icn1 ~= ic && icn2 == ic
icn = icn1; % need to sum over neighbours,
elseif icn2 ~=ic && icn1 == ic % so if not me, then its my neighbour
icn = icn2;
end
if icn ~= 0 % skip boundary faces
sumf = sumf + anb(ic,j) * phi(icn);
end
end

% GAUSS JORDAN method updates some
% vertices based on the currently solved, that is why cells
% that are not boundary appearing to be nonzero at iter 1
phi(ic) = (sc(ic)+scskew(ic) - sumf)/ap(ic);
if phi(ic) > 0 && printflag == 1
fprintf("ic %d: %6.3f %6.3f %6.3f %6.3f\n",ic,sc(ic),scskew(ic) ,-sumf,ap(ic));
end
end

% RESIDUAL CALCULATION CHECK MAX
sumr = 0; %
for ic = 1:ncells
sumf = 0;
for j = 1:3 % over neighbours
ifc = link_cell_to_face(ic,j);
icn1 = link_face_to_cell(ifc,1);
icn2 = link_face_to_cell(ifc,2);
if icn1 ~= ic && icn2 == ic
icn = icn1; % need to sum over neighbours,
elseif icn2 ~=ic && icn1 == ic % so if not me, then its my neighbour
icn = icn2;
end
if icn ~= 0
sumf = sumf + anb(ic,j) * phi(icn);
end
end
sumr = sumr + (sc(ic) + scskew(ic) - ap(ic) * phi(ic) - sumf) ^ 2;
end
res = sqrt(max(0,sumr));

ctrpsn = ctrpsn + 1;
end

end


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. – Johntra Volta Apr 3 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 – 2Napasa Apr 5 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! – 2Napasa Apr 7 at 16:21

## 1 Answer

The smooth solution turned out to have BC's applied in the following way:

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

• +1 Good work figuring out the answer, and thanks for coming back and writing an answer in case it's useful for someone in the future! – user1271772 Apr 13 at 22:39