[selig5576]

When operating on a collocated grid when solving the incompressible Navier-Stokes, do you need to do any special treatment for the pressure on the inlet? When solving the PPE I just apply unless I have an outflow in which I use . The method in which I am doing my collocation is the Rhie-Chow interpolation.

[FMDenaro]

To be honest I don't understand what are you setting...In inflow you set the known velocity profile and write the boundary condition dp/dn=(1/dt)(u*-u_inflow)

[selig5576]

I agree with the boundary condition you wrote. What I meant was the way I define the outlet BC for a collocated grid is . With this you plug it into . (for a projection method)

[FMDenaro]

Originally Posted by selig5576

I agree with the boundary condition you wrote. What I meant was the way I define the outlet BC for a collocated grid is . With this you plug it into . (for a projection method)

You first wrote about " a special treatment of the pressure at the inlet", no one being need other than prescribing the velocity profile.

At the outlet, prescribing du/dx=0 you get the relation for the second derivative of the pressure. We already discussed about that, in 2D such condition is equivalent to set v=0 at the outlet. The condition allow to write a 1D poisson equation according to d2p/dy^2 = (1/dt) v* you can solve with the BC.s at the two wall points. After solving this 1D equation you have a pressure profile at the outlet that can be prescribed as Dirichlet condition for the 2D Poisson equation.

[selig5576]

Yes you are right. I just wanted to make sure that the way I was treating the inlet in the PPE was indeed correct (not forgetting something), which it is. This is how I am doing it.

[FMDenaro]

Originally Posted by selig5576

Yes you are right. I just wanted to make sure that the way I was treating the inlet in the PPE was indeed correct (not forgetting something), which it is. This is how I am doing it.

In your chosen approximation d/dx=0, you could also simplify Eq.(1) as d2p/dx^2=0. That will produce the 1D Poisson equation at the outlet

[selig5576]

But that would mean that the R_{nx-1,j,k} would then retain the du^{*}/dx?.

[FMDenaro]

Originally Posted by selig5576

But that would mean that the R_{nx-1,j,k} would then retain the du^{*}/dx?.

I mean your physical approximation is d/dx=0 at the outlet for all the variables. Therefore du/dx=du*/dx=0 -> d2p/dx^2=0

[FMDenaro]

There is no a unique and exact way to prescibe the outlet condition, you always introduce an approximation. Its validity depends on the type of flow and the distance of the numerical outlet section.

[selig5576]

I entirely agree on what you say. The reason I posted this thread was I was leading up that my solver (pressure-free PM, AB2-CN integration) for the lid driven cavity in R^3 works perfectly. Even with a jet it works quite well. Though when I test my solver on a channel flow, after 70 iterations of SOR my PPE solver diverges. For sanity I have chosen a domain 80 x 40 x 40 with a kinematic viscosity of 0.01. If I make my inlet smaller then my solver under the same parameters does not diverge. Why would an inlet that spans u(1,2:ny-1,2:nz-1) cause a problem? It seems strange to me.

[FMDenaro]

Originally Posted by selig5576

I entirely agree on what you say. The reason I posted this thread was I was leading up that my solver (pressure-free PM, AB2-CN integration) for the lid driven cavity in R^3 works perfectly. Even with a jet it works quite well. Though when I test my solver on a channel flow, after 70 iterations of SOR my PPE solver diverges. For sanity I have chosen a domain 80 x 40 x 40 with a kinematic viscosity of 0.01. If I make my inlet smaller then my solver under the same parameters does not diverge. Why would an inlet that spans u(1,2:ny-1,2:nz-1) cause a problem? It seems strange to me.

divergence in the SOR applied to the pressure equation means you do not fulfill the discrete compatibility condition. That should correspond to an error in setting the BC.s

[selig5576]

I will have to do some testing on this matter. With that said, the way I implement the outlet BC is exactly how I did it in the PDF doc, but also at the 4 corners of my control volume (I have 8 corners because I'm doing 3D.)

!Left BC and right BC in the x direction  
!ufs, vfs, wfs refer to the face velocties of u^*, v^*, w^*.

do k=3,nz-2
  do j=3,ny-2
    !Left BC no outflow
     R(2,j,k) = (rho/dt)*(idx*(ufs(2,j,k) - un(1,j,k)) + idy*(vfs(2,j,k) - vfs(2,j-1,k)) + idz*(wfs(2,j,k) - wfs(2,j,k-1)))
     P(2,j,k) = (1.0-omega)*P(2,j,k) + (omega)/(idxx + 2.0*idyy + 2.0*idzz)*(idxx*P(3,j,k) + idyy*(P(2,j+1,k) + P(2,j-1,k)) + idzz*(P(2,j,k+1) + P(2,j,k-1)) - R(2,j,k))
!The outflow BC is only imposted at i=nx-1
     R(nx-1,j,k) = (rho/dt)*(idy*(vfs(nx-1,j,k) - vfs(nx-1,j-1,k)) + idz*(wfs(nx-1,j,k) - wfs(nx-1,j,k-1)))
     P(nx-1,j,k) = (1.0-omega)*P(nx-1,j,k) + (omega)/(2.0*idyy + 2.0*idzz)*(idyy*(P(nx-1,j+1,k) + P(nx-1,j-1,k)) + idzz*(P(nx-1,j,k+1) + P(nx-1,j,k-1)) - R(nx-1,j,k))
  end do
end do

As far as I can see and worked out on paper I don't see any issue in this implementation, but maybe another set of eyes will say otherwise

[FMDenaro]

Check the conditions at the corners where you have both walls and outlet. See also that the magnitude of the (div v) function is the same in each cell or you have some cell having greater magnitude. I assume you do not satisfy exactly the divergence-free constraint but only satisfy up to the local truncation error magnitude.

[selig5576]

I've spent considerable double checking there were no errors in the implementation of my pressure BCs. You are correct in that I am using an approximate projection method (APM), so continuity is only being satisfied up to the LTE.

The only observation I can make is if I reduce my kinematic viscosity to 0.1 from 0.01 my pressure solver does not diverge. At first I thought my cell Reynolds number was being violated but I ensured that I had a cell Reynolds number proportional to O(1).

For clarity I am working with a simple forward Euler scheme to reduce the 'uncertainity' in my code. I've attached two images: one image after 1 time step and the other after ~400 time steps.

EDIT: I also looked at the pressure profile and it seems consistent with what we should expect, pressure decreases as it gets farther from the inlet. I've also included a third image of the velocity profile at t = 1.5

Many thanks!

[FMDenaro]

In the first figure, I see two blue spots at the inlet that seem unphysical...

[selig5576]

I was about to comment on those. In my lid driven cavity test those do not appear. Only when I prescribe an outlet. When doing my SOR routine, I first iterate the interior then the boundaries, finally I prescribe the corners. This is all done in 1 while loop that iterates my PPE. I'm gonna assume I don't need a different while loop for the boundaries than the interior.

[FMDenaro]

Originally Posted by selig5576

I was about to comment on those. In my lid driven cavity test those do not appear. Only when I prescribe an outlet. When doing my SOR routine, I first iterate the interior then the boundaries, finally I prescribe the corners. This is all done in 1 while loop that iterates my PPE. I'm gonna assume I don't need a different while loop for the boundaries than the interior.

Differences in processing the nodes have effects in terms of convergence rate when using SOR. But in your case the problem appears in the final solution...I can suggest to do a simple test in a plane channel test prescribing the Poiseuille exact solution as initial condition and checking if it is conserved.

[selig5576]

I actually created a 'toy' 2D code for the same problem to maybe help illuminate a fundamental problem in my pressure solver. For 1 time step my pressure profile is correct however these 2 points gradually form near the inlet. To my understanding these are unphysical, but I am not quite sure.

[FMDenaro]

The inlet profile is prescribed as a plug flow? do you see the vertical velocity at the two corners?

[selig5576]

Yes, my inlet is that of a plug flow, in this case u = 1. My inlet profile seems consistent with that of my pressure profile. EDIT: This velocity profile is at step = 30 (not 1)