5
$\begingroup$

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): no sauce Adding the source term (set $u=(1,0)$ at the left faces):

sauce

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

![Updated_BC

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 3 at 21:09
  • $\begingroup$ Was thinking of L:$p=1$ R:$p=0$ TB:$\frac{\partial p}{\partial n}=0$ could you kindly correct me if im wrong $\endgroup$
    – 2Napasa
    Apr 5 at 19:06
  • 1
    $\begingroup$ @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! $\endgroup$
    – 2Napasa
    Apr 7 at 16:21
2
$\begingroup$

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$

BC

$\endgroup$
1
  • 1
    $\begingroup$ +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! $\endgroup$ Apr 13 at 22:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.