[eltenedor]

Studying Perić's CAFFA code (collocated FVM discretization and SIMPLE-type algorithm for Navier-Stokes) I came across two different implementations of the Neumann zero gradient outflow boundary condition. In the code the outflow boundary condition for scalars is implemented as

      DO IO=IOS(K)+1,IOS(K)+NOUT(K)
        IJP=IJPO(IO)
        IJB=IJO(IO)
        FI(IJB)=FI(IJP)
        DPX(IJB)=DPX(IJP)
        DPY(IJB)=DPY(IJP)
        CALL FLUXSC(IFI,IJP,IJB,IJO1(IO),IJO2(IO),FMO(IO),C,C,
     *              ONE,ZERO,FI)
      END DO

in another version of the code the same boundary condition is implemented as

      DO IO=IOS(K)+1,IOS(K)+NOUT(K)
        IJP=IJPO(IO)
        IJB=IJO(IO)
        DPX(IJB)=DPX(IJP)
        DPY(IJB)=DPY(IJP)
        CALL FLUXSC(IFI,IJP,IJB,IJO1(IO),IJO2(IO),FMO(IO),CP,CB,
     *              ONE,ZERO,FI)
        AP(IJP)=AP(IJP)-CB
        SU(IJP)=SU(IJP)-CB*FI(IJB)
      END DO

The first version implements the Neumann zero gradient boundary condition implicitly, i.e. if is a boundary control volume and corresponds to the outflow boundary face of the volume, the boundary condition implies . This would lead to a vanishing diffusive term for this boundary face. The convective contribution would depend entirely on . This is accounted for at a later stage in the code.

On the other side, the second version uses the zero gradient condition to create a "lagged" (by one outer iteration) Dirichlet boundary condition. As before but this time in an explicit manner. With respect to the implementation, this boundary condition is handled as if it was an inflow (Dirichlet) boundary condition.

I have a guess on why someone would decide to choose either one of the two previous approaches and I would appreciate getting some comments on it, maybe even a reference I could read, since I didn't find it in Ferziger/Perić nor in Versteeg:

The first approach favors convergence with respect to the number of outer iterations needed, because the boundary values are treated in an implicit way. The second approach (explicit) favors convergence with respect to the solution of the linear system at each outer iteration, because the explicit treatment of the boundary values increases the diagonal dominance of the linear system to be solved.

While I have found both described implementations used for scalar transport, I have only come across the second approach for the treatment of outflow boundary conditions for the velocities. Do you have an idea, why the explicit treatment of boundary values at Neumann zero gradient boundaries is recommendable especially for the velocities? Is it maybe, because fixing the mass flux at outflow boundaries is necessary for the velocity outflow boundary condition to be compatible with the Neumann zero gradient boundary condition for the pressure correction? (p. 206 of Ferziger/Perić Computational Methods for Fluid Dynamics (3rd Ed.) seems to imply this)