On Boundary Conditions
 

I wrote a two-phase, 1-D CFD code. It has a segregated solver and a full coupled solver. I just can't seem to get the residuals to behave. They keep oscillating all over the place. Furthermore, my pressure profile looks completely unrealistic for most of the cases I test. I'm trying to model some experimental tests of two-phase flow in a pipe that discharges into a vessel full of water. The boundary conditions I give to my code are:

1) mass flow rate at the pipe inlet
2) pressure at the pipe outlet
3) volume fraction of the vapor at the inlet (the experiment gives the flow quality at the inlet)
4) volume fraction of the vapor at the outlet (since the pipe discharges to a large vessel filled with water, I set this value to zero)

I'm wondering if I'm allowed to specify the volume fraction on both ends of the pipe like I have. It looks like the void profile tries to take a shape but then suddenly jumps down at the end to satisfy the boundary condition. Could this be the cause of the unrealistic pressure profile and the ill-behaved residuals?

I can't speak to this specific problem but typically it is not a good idea to use Dirichlet conditions on both the inflow and outflow boundaries. You may be over constaining the problem or at the very least making it difficult to converge towards the correct solution. Try instead a Neumann condition at the outlet. Provided the physics would have the fluid leaving as pure liquid, you should still get the result that you want on this boundary.

I agree with cdegroot's solution. Setting a fixed value both in inlet and outlet leads always to inestabilities. Try fixing a gradient (first or second derivative) of the vapour concentration and see if it improves.

I imagine it can even have more physical sense.

I never really thought of specifying a Neumann boundary condition. Out of my set of governing equations (vapor momentum, liquid momentum, vapor continuity, and liquid continuity), what derivative would I be specifying and to what value?

The experiment I'm modeling gives the mass flow rate and the flow quality at the inlet. It also gives pressure at the outlet. The facility is setup such that the pipe discharges into a large vessel filled with water, so volume fraction is zero at the outlet.

It's risky to say it without knowing in detail what the experiment is about, but I would say derivative of the vapour concentration in space. If the vapour concentration is really zero, then it will also be true that it does not change in space. You will close the equation system without the need of forcing the mass balance.

I'm not sure if I'm right, but I have a hinch. ; )

Yes I can see this is a bit tricky. Certainly it makes sense to specify your velocity (mass flux) at the inlet and simply extrapolate your pressure to the inlet boundary. This is a very robust boundary condition and will surely be stable. Now, this type of boundary conditions is most robust when coupled with some sort of "outflow" boundary condition where the derivative of velocity and the value of pressure is set.

I would set the derivative of velocity to zero at the outlet. Presumably the volume fraction of zero is attained within the pipe, before it enters the vessel so this should accomplish the same thing, i.e. your solution with the Neumann outlet condition should return a volume fraction of zero at the outlet. If it doesn't, that probably means the liquid enters the vessel as a mixture and becomes liquid within the vessel. Hopefully this is not this case because your 1D model would not handle this. You would need to also model the vessel.

In summary:

At the inlet: specify velocity to satisfy the desired mass flow rate (Dirichlet condition) and extrapolate the pressure (Neumann condition).
At the outlet: specify pressure as the measured pressure (Dirichlet condition) and set the derivatives of the velocities to zero (Neumann condition).

Thank you both for your recommendations. I will have to play around with these ideas today and report back when I get some results.

I imagine it's possible that the vapor concentration goes to zero before actually discharging to the vessel. The experiment facility contains a gamma densomiter for measuring the void fraction about 1.4 m before the vessel junction. The volume concentration is very far from zero at that point, but the void gradient in space would be very steep after that.

Considering the actual experiment, the velocity gradient would not be zero at the junction because the flow area expands from the pipe area to the much larger vessel area; however, I do not model the vessel region, so for my code's sake, the velocity shouldn't change in that last momentum mesh cell so long as the volume fraction also does not change.

I wonder, though, that, while a zero velocity gradient at the exit will ensure that the volume fraction does not change across that cell, what's to force the volume fraction to go to zero at that point?

Presumably it would be the physics of whatever is going on in the pipe. For example, if there is a wall heat flux cooling the fluid. Outlet boundary conditions are a tricky subject. It is typically better when they do not try to influence the flow, rather just let it exit gracefully.

But you are right in a sense that there is no real reason that it has to go to zero, other than the fact that the tank is at zero and I am presuming it doesn't reach that state immediately at the tank. Hopefully that is true, otherwise you are going to have to devise something a bit innovative to specify a non-zero derivative at that boundary.

The problem I'm seeing, the more I think about this, is that there is not really any physics in the pipe that could account for the volume fraction going to zero. Actually, there is not heat flux since this is adiabatic. As far as I see, the volume fraction goes to zero because it discharges into a tank filled with water. If I modeled the tank, then I believe the volume fraction would naturally be driven to zero, but since I don't actually model that tank, I can't see another option but to specify that the volume fraction must go to zero at that point. I wonder what the result would be to not specify the inlet void fraction.

Ah very interesting. Maybe you can specify the volume fraction at the outlet and let it float at the inlet. Possibly a zero derivative condition would work at the inlet if it is far enough from the tank and is fully developed.

Yes, that's what I was thinking of trying. I'd imagine that, in the experiment, inlet volume fraction was not something they directly controlled, but was a secondary effect of modifying the pumping rate and the amount of enthalpy put into the fluid before entering the test section, so maybe if I approach this the same way, that will work.

I think you are right it should work. As long as you do not specify on both boundaries it should be okay.

I tried setting the volume fraction at the outlet and letting the inlet volume fraction drift with the solution. Interestingly, some of the cases saw the inlet void drift to the measured inlet value for volume fraction. Others saw it go far from it. The residual behavior didn't improve at all, nor did the pressure profile, so it seems there's something else causing that. There are some cases that have a much harder time converging and if those are the same cases that saw the voids drift to the wrong value, then I wonder if there is something wrong with the way I'm modeling them numerically in the code. So while my problem is not yet solved, this may point me towards some new insights to pursue.

On another note, I have been playing with the void fraction I use for the momentum cell faces (I use a staggered mesh) and I saw some improvements by doing that, which prompts me to believe that there may be something wrong with the interpolation approach I use (first order upwind scheme) to calculate mesh cell face values. By increasing the mesh refinement a whole bunch, I see quite different solutions emerge and, notably, I see much steeper void gradients at the end of the pipe. I believe that, with the coarser mesh, I was seeing a more gradual gradient because of the diffusion effect of using first order UDS. I'm currently implementing a second order upwind interpolation scheme (QUICK) to see what effects that has.

Thanks for all posts