Velocity profile disturbance due to loss coefficient
 

I wrote a code that simulates a two-phase flow through a horizontal pipe. The boundary conditions are pressure at the outlet and mass flow at the inlet. My code uses the SIMPLE algorithm to couple pressure and velocity. It only solves the momentum and pressure correction equation at this point.

When I run a simulation with no loss at the pipe exit, the velocity profiles and pressure profile looks realistic. The liquid phase velocity decreases linearly from inlet to outlet and the vapor phase velocity increases linearly from inlet to outlet. But when I apply an exit loss coefficient to the last momentum cell (cell N), a ripple forms in the velocity profiles. At cell N-1, the liquid velocity suddenly increases instead of continuing its decent. In cell N, the liquid velocity drops back below the value in cell N-2. The vapor phase velocity does the opposite, decreasing in cell N-1. Since one phase zigs while the other zags, mass is conserved and the pressure correction equation allows the behavior to stay.

I ran this using only single-phase liquid by actually getting rid of the vapor phase equations and there is no such ripple, with or without the loss coefficient. But when I run it as two-phase, but with a very tiny vapor void fraction, the ripple still persists, so it seems to be due to the presence of the vapor equations and not just due to the value of void (which is fixed throughout the pipe since I don't solve phase continuity equations right now). I also shut off the mass and momentum transfer terms at the vapor/liquid interface to rule that out and the ripple still persists. So it looks like it's occurring through the phase coupling in the pressure correction equation.

I've seen this problem before, but I can't figure out what causes it. As far as I know, such a ripple in the velocity profile is not physical. So does this mean that I have improperly applied the loss coefficient or is there a problem with my formulation of the conservation equations?

Since a picture is worth a thousand words, I plotted the liquid velocity profile (normalized) along the horizontal pipe length. The plot shows the momentum cell number on the x axis and the normalized liquid velocity on the y axis. Since the inlet velocity is the given, the velocities are normalized to inlet velocity.

The velocity profile is shown for the first 4 iterations and then the 500th iteration, which is where the solution converges. There is an exit loss at the end of the pipe, which is why the velocity drops so much in the last cell during the first iteration. But then the velocity profile slumps in the middle of the pipe next iteration. By the third iteration, the profile goes to a fairly steady decrease until cell N-1, which pops up. This pop holds up until convergence and you can see the converged profile has a ripple in it.

I've looked over my code and I think I applied the algorithm correctly... I'm wondering if my boundary conditions are not applied right.

I also plotted up the velocity profile behavior through the iterations.

http://i898.photobucket.com/albums/a...d_velocity.jpg

http://i898.photobucket.com/albums/a...r_velocity.jpg

I did some void sensitivity studies, adjusting vapor volume fraction from 0.1 to 0.9 in 0.1 increments. It seems my wall shear model may be causing some problems because the vapor phase is being driven to zero flow when vapor volume fraction is high. I'm using the Lockhart-Martinelli two-phase pressure drop model. I distribute the resulting wall drag force to the phases by multiplying the force by the phase velocities. See the results of the void sensitivity study:

http://i898.photobucket.com/albums/a...71/liq_vel.png

http://i898.photobucket.com/albums/a...71/vap_vel.png

The pressure drop in the pipe looks bizarre, as I expect it to be linear, making me again think there is something wrong with my wall shear formulation.

http://s14.postimage.org/t9vuov1q9/pressure.png
gif hosting

Now, if I get rid of the form loss at the end of the pipe:

http://s19.postimage.org/mm8ztop5f/liq_vel.png
free image hosting

http://s19.postimage.org/4gw1fmpn7/vap_vel.png
picture hosting

The pressure distribution still looks messed up:

http://s19.postimage.org/5x7jxrsk3/pressure.png
upload pictures

I checked and for all these cases, the two-phase mixture mass flow rate is conserved throughout the pipe and is equal to the inlet boundary condition. My void fraction equations are on, but void is predicted not to change throughout the pipe, it stays equal to the inlet value. Since the last post, I also normalized the continuity equations to reference densities. This is supposed to lead to better conservation of the lighter (vapor) fluid. These equations are used both in determining the pressure correction and for determining the phase void fractions.

It seems to me that the problem is the way the mass fluxes are defined at momentum cell faces. I use the upwind difference scheme. Thus the convection term:

d(rho*u*u)/dx

is expressed as,

(rho*u*u)_e - (rho*u*u)_w

or,

F_e*u_e - F_w*u_w,

where F_e is the east face mass flux (rho*u)_e and F_w is the west face mass flux (rho*u)_w. We solve our momentum equation for u_e and u_w. The F terms are explicit from the last iteration. This is the linearized form of the momentum equation. This is what this control volume looks like:

http://s19.postimage.org/vzsx71wsz/cv_patankar.jpg
picture hosting

Note, however, that velocity is not defined on the faces of the momentum cell, but in the center. The upwind scheme says that:

u_e = u_P if F_e>0 and
u_e = u_E if F_e<0

Similarly,

u_w = u_W if F_w>0 and
u_w = u_P if F_w<0

This basically means we take the upwind velocity as the velocity to convect. The problem is, how do we define the u in the F terms? I thought it should be defined as u_e = 0.5*(u_P + u_E). In other words, a linear interpolation between the adjacent momentum cell values. From the books I read on the upwind scheme, this is how it should be done. But it confuses me because then we are still reaching into the downwind side of the solution to get velocity information. So I changed my code a bit to define the u_e in F_e as:

u_e = u_P if u_P > 0 and
u_e = u_E if u_P < 0

The same way as the upwind difference scheme when formulating the momentum equation. I didn't think this was the right way, so somebody please correct me on the right useage of the upwind difference scheme if this is not correct. However, my results look normal with that implementation. Note that in these runs, I moved the form loss to momentum cell number 10 to more clearly show how velocity reacts before and after the loss.

Old implementation of UDS:
http://s19.postimage.org/aovd2seoj/liq_vel.png
free image hosting

New implementation of UDS:
http://s19.postimage.org/6egp779lf/liq_vel.png
upload pictures

As for the vapor phase:

before:
http://s19.postimage.org/z7xek3j2r/vap_vel.png
upload pictures

after:
http://s19.postimage.org/kq07c3rrn/vap_vel.png
adult upload image

The liquid velocity still looks weird for the higher void cases, which I think is attributable to my wall drag still not being formulated correctly. It shows also in my pressure drop in the pipe. When void goes over 0.5, the relative pressure in the pipe goes negative after the form loss and actually rises as moving to the exit.