Suggestions for help are requested. I have been developing a general purpose incompressible structured grid CFD code for about a year, primarily to learn CFD from the ground up. I apply an ADI-Brian solution to the momentum equations with a SOR solution of a pressure Poissons equation. The code is general purpose, easily configured and quite functional. However, i am constantly plagued with small divergences in the velocity field. This is hitting the conjugate thermal solver which generates artificial heat sources with the presence of the Div(V) terms remaining,e.g. stick a 100 Watts at a wall on the backwards facing step and one gets 500 Watts equivalent enthalpy flux at the outlet. The problem , i think, is the pressure Poissons equation solution. All solution states [u,v,w,P] reach steady conditions but the divergence of the velocity field is not zero in locals. However, the sum(Div) taken over the whole field from an inflow to outflow region is zero. I have read some papers indicating divergence problems with the PPE solutions. Further, i believe it is related to the setup of the boundary conditions. If anyone can share in this area it would be greatly appreciated. Thanks for your help.

Regards

Dean Schrage

The locale Div(v) should be withen the numerical error.

If the sum(div) cancels over the whole field, that would mean that the boundary conditions are well posed. Indeed, by the Gauss theorem sum(div) = flux(u.n) on boundaries, which should cancel for an incompressible fluid... I would rather check if your poisson equation is enough accurately solved, so that the Div(v) is small enough locally. By iterating long enough, you should be able to reach machine accuracy for the residuals of the equation (normalized by the diagonal term usually). You could also check that the div(v) which is involved in the thermal equation is correctly discretized with respect to the poisson equation. Indeed a lack in the consistency between the two equations could produce the discrepencies you observe.

Hope this help. Good luck,

franck

yes part of the problem i observe is that between two connecting nodes i can have both positive and negative divergence. The source term for the PPE dumps this information from the + div node to the - div node rather than adjusting the pressure over the field more appropriately. the analogy would be heat addition and extraction for two adjacent points. Thus, when i have divergence remaining it usually isolated to adjacent points. Flows over corners generate + div and then the CV just below the corner has a - div. Also, all interations are taken to steady state. The pressure field, given my BC's minimizes the divergence field but does not eliminate it. thanks.

regards

Dean

Hi Dean. I´ve read your question about the divergence in velocities fields. I´ve work with discretization schemes for numerical simulation. Some causes of divergence sometimes is due to the kind of scheme, such as Central, Exponencial, Power-Law or QUICK, you are using to interpolate velocities or temperature and the algorithm of coupling momentum-pressure equations, such as SIMPLE,SIMPLER or PRIME. You said about "artifical heat source", i think this is caused by some numerically diffusive scheme like Upwind or some related scheme that your are using. Try to use other kind, less diffusive like QUICK. The pressure field is dependent of them velocities field, so if the velocities are wrong, the pressure will be too. If you are using SIMPLE algorithm, try to use SIMPLER, this is easy to implement and more stable than SIMPLE. For more informations, you can read the Patankar´s text book.

Success

Carlos Vilela

Hi,

I think I have seen this one before, but I have to start with a bit of theory:

When you derive the pressure equation you start from div(U) = 0, which means (with all else consistent) that the continuity equation will be satisfied to the solver tolerance (this is crucial!). You may be experiencing one of these:

a) when calculating the face flux, the expression MUST be equivalent to a_N*(p_N - p_P), where a_N is the pressure matrix coefficient for the face and p_N and p_P are pressures around the face. Once you get this right (watch out for te part of the flux that ends up in the source!) the flux sum for the cell will be zero (for incompressible flow).

b) There's only one way of calculating the div(U) term: from the FACE FLUXES and not from the interpolated velocities:

sum_f (S . U_f) will give you the value out by the discretisation error and ot solver tolerance!!! (here, S is the face area vector and U_f = f_x U_P + (1- f_x) U_N is the interpolated face velocity. You should use sum_f flux instead.

Now, if your flow is compressible, the same applies! (our face flux is F = rho_f * S . U_f). When you need div(U) you do

sum_f F/rho_f, i.e. you interpolate the density and NOT THE VELOCITY!

I hope this helps - could you please let me know if I got it right (it's always nice to get some feedback on cases like that).

Hrv

thanks for you reply. Some more details with which to hone in on the problem:

1) non staggerred grid is applied with which it made more sense to apply a PPE rather than SIMPLER. Have not even worked up a SIMPLER scheme so cannot comment on its usage

2) yes all DIV(V) calculations are based on face interpolate velocities. This is especially crucial for CV's attached to solid impenetrable boundaries.

3) GRAD(P) terms: This has evolved considerably. Originally i coded this using central differences as you noted above. However, this lead to pressure velocity decoupling at higher Reynolds numbers. Recently, i applied a momentum interpolation scheme. GRAD(P) = W(Pe) GRAD(P, upwind) + (1-W(Pe)) GRAD(P, Downwind). The Pe is the grid Peclet number. Thus on high Peclet numbers, the code convertes to an upwind differencing on Pressure radient. This appears to be crucial to the PPE because if a node has a Div(V) = -() number then the source term for the PPE is negative. Looks like a heat equation with heat source added. Thus the pressure field curvative is negative over this node. Solving the PPE, the node pressure properly increases to attempt to eliminate the added mass indicated by the negative divergence. However, with center differencing, the pressure gradient picks up the surrounding points in differening and does not increase to blow out the added mass. Hence, pressure gradient based on velocity field.

4)

I am not sure where your problem lies but a couple of pointers which may help (or not).

The pressure equation is the divergence of the momentum equation. On its own this is useless since it introduces no more information than is in the momentum equation (in fact less since it is differentiated). What makes the pressure equation usable is how continuity (that extra bit of information) is introduced into it. To simply set the divergence to zero at the current time step and start with a continuity satisfying field does not work and instabilities will develop through time (discussed in the mid 60s - see various Los Alamos publications). The solution is to include the continuity error at the previous time step in your formulation. In fact, a trivial scheme which retains only the time derivative of the continuity error on the RHS can often be made to work. This is not particularly sensible but does emphasise the role of the pressure equation. To follow this line of thinking further will, of course, lead to pressure correction schemes.

A Poisson equation without Dirichlet boundary conditions has no solution unless the rhs sums to zero. This is usually achieved by evaluating the mass fluxes conservatively and forcing the inflow to balance the outflow.

You state that you have solutions, global mass conservation but violate local mass conservation. You have introduced pressure smoothing through upwinding. The obvious question: if you include your pressure smoothing terms when evaluating your local mass imbalance does it disappear? (i.e. like a Rhie-Chow scheme).

Thanks for your reply Andy. I will annotate your comments as best i can. Confusion on your last question though.

Yes, i apply the residual divergence as a source in the PPE source term. Typically, the - rho d(div(V))/dt term is the largest of the source components. It is included in the same fashion as the MAC method, i.e. -rho Div(V(n-1)) / dt.

I observe and enforce a strict Neumann condition if the system has closed boundaries such as a cavity driven flow. The internal source is then adjusted to comply with the current-time boundary conditions. However for in-flow/out-flow problems, i apply a Dirichlet BC on the outlet (P=0) and force the normal derivative of velocity to be 0. No source adjustment is applied. The PPE converges slowly but i do not try for a time accurate solution. By lagging the PPE solution behind by one time step, i get fairly robust rapid convergence by taking 1 sweep through the PPE for each time advance of the momentum equations. There are similar schemes where the PPE is made parabolic with an artificial time derivative. In fact the relaxation factor i use makes this operation almost the same (within 10 percent). Thus, i dont think there is anything broken here.

Not sure how to comment here. Could you restate this please. Perhaps you ask, "does the Div(V) reduce itself by momentum weighted pressure gradient eval's." Yes, this did help but still does not eliminate the divergence. Always, the presence of divergence is in the vicinity of boundaries with sharp gradients. And typically it is with two nodes near each other (as if one dumps its divergence effects to the other). The divergence effects are generally not observable in the CFD simulation results but as i stated above in other posts, the divergence produces artificial heating terms in the enthalpy advection term in the energy equation. As my overall goal is a thermal solver based on the CFD solution, i need to wipe out the divergence to make this thing worthwhile.

I still focus on the BC's of the PPE. This is where the game is won or lost. Out in the bulk flow field, the PPE converges beautifully. Div(V) terms are reduces to convergence residuals. Not near boundaries.

Thanks for you thoughtful comments and assistance. Hope to continue. Thanks to everyone also.

regards

Dean

***** You state that you have solutions, global mass conservation but violate local mass conservation. You have introduced pressure smoothing through upwinding. The obvious question: if you include your pressure smoothing terms when evaluating your local mass imbalance does it disappear? (i.e. like a Rhie-Chow scheme). *****

error in post which may be confusing. The repair is:

... It is included in the same fashion as the MAC method, i.e. rho Div(V(n-1)) / dt.

A traditional Rhie and Chow scheme will count the pressure smoothing terms at cell faces as additional mass fluxes (i.e. the terms involving pressure which are not an interpolation of the velocity field). I was suggesting you did this in order to examine the faces where the positive link between "mass" flux and pressure gradient was being lost. If your scheme is conservative, the pressure smoothing "mass" fluxes will be of the opposite sign to that desired.

As a guess, I would tentatively suggest it is the treatment of the pressure smoothing at the boundary that is tripping you up.

Is your scheme conservative?

If you have both velocity components moving towards a cell face or both away from a cell face how do you evaluate the pressure gradient?

What is the numerical relationship between the pressure gradient in the momentum equations for two adjacent cells and the pressure gradient at the shared cell face in the pressure equation. That is, how do you correcting (or uncorrecting!) pressure smoothing fluxes arise.

can't comment on this. Perhaps will need to get the paper to study this. ***** A traditional Rhie and Chow scheme will count the pressure smoothing terms at cell faces as additional mass fluxes (i.e. the terms involving pressure which are not an interpolation of the velocity field). I was suggesting you did this in order to examine the faces where the positive link between "mass" flux and pressure gradient was being lost. If your scheme is conservative, the pressure smoothing "mass" fluxes will be of the opposite sign to that desired. *****

The pressure gradient is comprised as forward and backward components. If the Peclet number is -10 it is weighted more heavily in the forwards direction than backwards. If Pe = + 10 then more heavily backwards than forwards. The Peclet number is based on the cell centered velocity. I did this for simplicity. I see some schemes apply the momentum weighting on the face interpolated Peclet numbers. However, in doing this i felt that there would be a disconnect with the method. Example if Pe(i-1/2) = +10 and Pe(i+1/2) = -10 then the scheme would give a much larger value for the pressure gradient (basically doubling the gradient). Perhaps this is more correct, but at the time of coding it didn't seem so (just my gut feel).

***** see above If you have both velocity components moving towards a cell face or both away from a cell face how do you evaluate the pressure gradient? *****

the PPE solve the boundary nodes by solving the residual equation:

GRAD(P,en) = mu DELTA(V.en)

where DELTA(V.en) is evaluatate at the face of the CV's just inside of the fluid domain in contact with a pressure boundary.

Thus at a boundary, i take one step into the solution domain and press, up or down, the wall pressure based on the above equation. I think the answer to your question is forward differencing is applied to the BC's of the PPE and the revised Peclet weighting scheme i described above is used in the momentum equations. With this scheme, the momentum equation is always in contact with the boundary to some degree.

***** see above What is the numerical relationship between the pressure gradient in the momentum equations for two adjacent cells and the pressure gradient at the shared cell face in the pressure equation. That is, how do you correcting (or uncorrecting!) pressure smoothing fluxes arise. ****

thanks

regards

Dean

I am not sure we are quite connecting (the medium is not conducive to this sort of discussion) and I do not fully understand the differencing operators you are expressing. However, it may not be important to the problem (assuming I understand it correctly, your scheme is conservative and the pressure equation is derived from the discretized momentum equations).

If your want your scheme to numerically locally conserve mass then this must be reflected in the differencing of the pressure equation because this is where continuity is imposed. If your pressure field is not to decouple you must introduce some form of pressure smoothing (by averaging, one-sided differencing, time stepping or whatever is not important). Where you have a mass error across a cell face you will find that your smoothing terms (however your differenced pressure equation differs from a simple centred scheme) are generating mass fluxes which increase rather than decrease the mass error. To see this, take your solution, identify the worst offending face, perturb the pressure gradient across the face and work through the resulting change in the velocity field. This process will identify what is different about the regions which are behaving and those that are not.

PS What stops the pressure field decoupling for Peclet numbers less than 10?

(1). I don't know what you are doing. And I guess I don't want to know it. (2). But I do have a suggestion here for you. Take a fully developed 2-D channel flow problem as a test case and check out your code thoroughly first. (3). First prepare a long 2-D channel, at the inlet prescribe the fully-developed parabolic velocity profile. And then solve this problem using your code. (4). you can also check out the energy equation by specifying the wall temperatures. (5). This problem has the analytical solution which you can use for the validation of your code. (Do not specify the exact solution as the initial flow field.) (5). Then you can move on to more complex flow problems, such as the developing 2-D channel flow, flow over a back-step,...etc..(6). If you still have some problems, try to follow existing methods, to see whether you can get it to work. If you are not satisfied with the mass conservation, you can still try the stream function-vorticity method.

You need to buy, or at least read, our book on the subject,

"INCOMPRESSIBLE FLOW AND THE FINITE ELEMENT METHOD:, By P.M.Gresho and R.L.Sani, John Wiley & Sons; 1999 for the extant hardcopy version, and 2000 for the SOON-TO-BE-RELEASED 2-part paperback version....

While emphasing FEM, the book covers ALL saspect of CFD for incompressible flows, with VERY MUCH EMPHASIS on the pressure Poissom eqn!!!

Andy, I think I have solved the problem. The code was applying the corner pressure boundary condition on a wall to both N and W facing nodes. The north node was derived using the proper boundary condition equation (Neumann in V.ey) while the W facing node applied (erroneously) the same condition. I am restructing the preprocessor for this.

Second, the maximum divergence is a function of the time step. I found a paper which analyzed non staggered grids showing that the divergence cannot be zero - only small and proportional to the group (dx^2 dt). This was exactly observed. Intersting, the pressure field becomes very spiky and oscillatory when the time step is reduced. The flow appears to become incompressible in simulation (note that i lagged the pressure and velocity solutions without update at each time step).

Thanks for your help. You seemed to be the most consistent in your attempts at help. Some of the other postings were downright goofy.

best regards

Dean

Thessaloniki 11/04/2000 Hi Mr College I am postgraduate student at the national technical university of Athens, Greece and i have worked on Pessure Correction Techniques(Simple Algorithm) on unstructured Meshes that consist of 2-D triangular elements. Although the code converges i have great problems as regards the arising pressure oscillations.I'd like to know your opinion