Most users take it for granted that CFX 5 is one of the most reliable commercial CFD programs, and the values in CFX Post are very accurate. There is no doubt that it is an excellent software package, however I have found some weird features that should be seriously examined and clarified if one cares even a little about accuracy of the results.

In many engineering tasks 10-20% inaccuracy is acceptable, but there are applications specifically in theoretical and research projects where this much inaccuracy makes a big difference. Here I do not mean the accuracy of the numerical solution compared to the measured values from experiments, because it is natural and understandable that it is very difficult to set up a mathematical model of physical processes that would solve the problems with only few percent (or even less) error. I am talking about the mathematical accuracy of the solution assuming that the mathematical model of the flow process is accurate. Let us take a primitive example to check the reliability of the derived solutions.

The geometry:

The geometry is very simple, it is a nozzle of hexahedral shape. Draw a tetragon using the following coordinates for the corners in [mm] (x,y): (0, 0); (0, 10); (40, 3); (40,0), then extrude it in z direction for 10mm width. Make the face of the drawing to be symmetry plane. The inlet is the wider face of 10mm*10mm, and the outlet is the narrowest face of 3mm*10mm. The remaining 3 faces are smooth non slip walls.

Before starting the solver: Tools/ Edit Definition file > Add expert parameter section > Add parameter > output eq flow. After starting the solver, around iteration 13 > edit run in progress > change Auto Timescale to Physical Timescale; delete Length Scale Option; and add Physical Timescale = 1 [s] > save. This change will accelerate the convergence of the Heat Transfer. This is just an example how I got the discussed weird results, and not necessarily the best; one can take even greater time step but it will not make much difference. The solver finished within 40 minutes with the following residuals:

In Objects create the following Surface Groups: Wi is the bottom surface touching the x axis Wo is the top slanted surface above Wi Wcap is the one side wall opposite to the symmetry face.

Add the following expressions:

Tain = areaAve(Temperature )@in

Pfin = -massFlowInt(Velocity ^2/2)@in

Pfout = -massFlowInt(Velocity ^2/2)@out

Phin = massFlowInt((Temperature -Tain)*2400 [J kg^-1 K^-1])@in

Phout = -massFlowInt((Temperature -Tain)*2400 [J kg^-1 K^-1])@out

Ppin = -areaInt(Pressure *Velocity u )@in

Ppout = areaInt(Pressure *Velocity u )@out

Ptot = Pfin+Pfout+Ppin+Ppout+Phout

THE PROBLEM:

Using the above expressions we can check the conservation of energy. The total energy that enters and leaves the nozzle Ptot supposed to be zero. However, we get a value of Ptot=59.56W excess energy. If we compare this value with the amount of dissipated heat energy Phout=792.7W then we can see that the error is 59.56/792.7*100%=7.5%.

If we would have used coarser mesh and if we would have achieved inferior convergence let's say about 1E-4 for the equations, then I would say that this much error is due to the inaccuracy of the mesh and insufficient convergence. But the mesh is not so bad, and the convergence is also enough good (1.5E-6 3.2E-8) to make one suspect that something is wrong with the Solver or the CFX Post that causes this much unexpected error.

I have experimented with different sizes and shapes, and different mesh quality but this characteristic error appears in all cases when using the fairly viscous Glycerol, and it can go even to 10% in some cases even if using much finer mesh than in this example, driving the residuals down to 1E-11 and using double precision. The geometry and the flow type in this case is the simplest possible. If this much error appears in the results in such simple case then what can one expect in more complicated cases and with lower convergence? I hope you understand my dilemma

So WHAT IS WRONG? Is the solution correct and there is really about 8% excess energy or is there some bug in the Solver and/or CFX Post, or are the above formulas for calculating the energy balance incorrect? Is there any way or possibility to drive the error down to the range of 1%. If yes, how to do that?

Any input on this issue is greatly appreciated. If there will be some interest, I will show some more bugs and weird errors concerning the forces and torque later.

You have obviously tried to be careful in your analysis but your expressions are inaccurate because they are not exactly what the flow solver will assemble at the boundary conditions.

Rather than using those expressions in CFX-Post to check your imbalances, why don't you check the flow solver output file instead. Is there an 8% energy imbalance reported there for the H-Energy equation? i.e look at the sum of the H-Energy flows for your inlet and outlet and look at the volume source (viscous work in your case) and see what they add up to.

You might also turn off the viscous work term in the energy equation because it's probably unecessary and adds an extra ingredient for which you are not accounting for in your expressions.

As you can see it shows an extremely good accuracy that is 4 orders of magnitude below the real accuracy of the values that we can calculate in CFX Post. I just simply do not trust these overly optimistic numbers, because I can not access a solution and results in CFX Post that would be so accurate as indicated by these figures.

Let us compare these results indicated by the solver output file with the results derived from CFX Post. Using the formulas from my previous post the power that enters the nozzle is the sum of the pressure power, kinetic (flux) power, and heat power Pin=Ppin+Pfin+Phin. Pin=-1753.45 W this is negative in my notation because it is spent, we have to invest it, or we could say: we lost it. If you compare it with the Pin=1676 W displayed by the solver (disregarding the opposite sign), then it is obvious that there is already dPin=(1753-1676)/1753*100%=4.4% difference between the two values.

The output power is the sum of the same 3 components at the output Pout=1813.04 W. The sign is positive because we get this back from the nozzle. After comparing this with the value indicated in the solver output file dPout=(1813-1676)/1813*100%=7.6% we see that the relative error here is even bigger than at the input. Consequently there is no reason to thrust the indicated Domain Imbalance of -8.5284E-04 ie. -0.0001 %. It seem that the values indicated by the solver output file have not much connection with the values that we can derive from the CFX Post.

You say that my equations are not accurate because they are not the same that the solver uses. Well, if they are incorrect then I would be very thankful if you would show me what equations should be used to get the same results that the solver report indicates. There is nothing mystical in these equations, one can come up with them easily by using the expressions of kinetic power, the work done by the pressure volume, and the increased energy content of the fluid of constant heat capacity if we know the temperature change. These formulas simply sum up certain products for each elementary surface to get the desired results.

If the solver achieves a very accurate solution that would produce the indicated -0.0001 % error, then the same solution supposed to be given to the CFX Post for analyzing the results too. Then if we perform the integration of the elementary power products over the boundary surfaces then we supposed to get the same excellent accuracy. The value (and usefulness) of the solution is only as good as the accuracy of the solution that I can obtain from CFX Post, and it does not matter what the solver report indicates.

There can be the following possibilities, or the combination of them:

1. The Solver arrives to a very accurate solution but the CFX Post can access only a low accuracy solution field, and consequently displays low accuracy solution.

2. The Solver does not achieve the required high accuracy and the CFX Post correctly displays low accuracy solution, consequently the solver report indicates false approximate values (using some cosmetics).

3. The Solver gives a high accuracy solution field to the CFX Post, but the CEL functions work incorrectly and do not perform accurate integration and calculations.

4. My CEL expressions used for calculating the power balance are really incorrect, thus by using the correct CEL expressions we can get the indicated very high accuracy even in CFX Post.

In the first 3 cases we can not help but be aware of the limited accuracy of the software and take the indicated values with due reservation (preferably by performing the indicated integration). If the 4th case is the only reason for the inaccuracy then let us find the correct formulas and get confidence (and practical knowledge) in the accuracy of CFX 5.

The viscous work term should be included into the energy equations because almost half of the input energy is converted into heat energy by viscous friction in the nozzle (see the above partial power values). I have detected several bugs in the CFX Post, but I am also convinced that the solver report does not reflects the complete truth but values tuned to some expectations. We can discuss the other bugs later, but let us find the solution for this problem first.

1. CFX solves the total energy equation, in which the total enthalpy is defined as

Ho = h + v^2/2 h = h_ref + Cp*(T-Tref)

You've included the 'Pv' contribution explicity, whereas CFX incorporates it into h. So it is incorrect to include the Pp contribution separately from your Ph contribution; as it is, you've double-accounted for flow work.

You can estimate the enthalpy advection through the boundaries simply by doing a massFlowInt of Total Enthalpy. If you want to see it separated into thermal and kinetic effects, do a massFlowInt of Static Enthalpy and of kinetic energy (which you'll have to calculate yourself).

2. Viscous work also gives an energy flow which you've left out of your analysis. It represents the work done on the system by viscous stresses on the boundary. These are a little harder to evaluate in Post because Post doesn't have direct access to the stress tensor. A first sanity check would be to turn off viscous work and see that the balances make sense then.

where massFlowInt is the sum of the integration point massflows times the scalar value (in this case, total enthalpy). This is not what you have. eg, for the inlet you have used

Phin = massFlowInt((Temperature -Tain)*2400 [J kg^-1 K^-1])@in

Pfin = -massFlowInt(Velocity ^2/2)@in

Ppin = -areaInt(Pressure *Velocity u )@in

The flow solver includes Ppin in the enthalpy already. and there will be small differences in how CFX-Post calculates Phin+Pfin when compared to how the flow solver substitutes boundary conditions.

Also, as Phil and Neale say, you also do not account for the viscous work term in any of your expressions. This has an effect.

You also mention that you expect the sign of the flow through the inlet to be negative. This is not how CFX-5 works. Flows coming into the domain are considered positive, flows leaving the domain are negative.

Another thing you have misunderstood is the concept of accuracy vs. conservation. These are two compeletly different things. The degree to which the flow solver is conserving flows has nothing to do with accuracy (eg: under a mesh refinement or timestep refinement). The imbalance and residuals only indicate how well the flow solver has solved the equations on the mesh and boundary conditions you have thrown at it.

Thank you for the suggestions. I am convinced that I am not the only one on this forum who will benefit from our discussions, so your efforts to help finding the solution is very valuable.

First of all let me clarify that my intention with these calculations was not to find fault or play the nitpicking game. The very simple reason is that my work demands that all types of energy components should be known separately with high accuracy. From this demand naturally follows that when all these components are summed up, they supposed to give a residual as close to zero as possible. By summing up the energy components one can also check or at least estimate the accuracy of the energy components. One can feel a certain confidence in the solution if besides knowing the values, one also knows the accuracy of these values.

I wanted to avoid using enthalpy for my calculations because it is a "slippery" stuff that represents the sum of certain energy components. Now from this "conglomerate" one can separate each component only by calculating some components separately and subtracting from the total enthalpy. Why complicate the case when not necessary? Why can I not calculate each energy component separately as I did with the formulas in the first message? Why would these separate components be incorrect? And if each one is correct then why supposed to be their sum different from the result calculated through the total enthalpy, or the value indicated by the solver report?

OK, let me try to comply with your suggestions separately and see what we get, then we can come back to the clarification of each energy component separately.

Phil, You say that I should not calculate the energy component of flow work ie. kinetic energy or in my notation: Pfin = -massFlowInt(Velocity ^2/2)@in (and the same at the outlet) because this is already included in the Htot total enthalpy, and thus I have double accounted for this component. This objection might be correct only if I would have used Htot for calculating the energy balance. But this is not the case; I have absolutely not used enthalpy for calculating the energy balance.

In this case of incompressible liquid with constant heat capacity and laminar flow I can conceive only 3 different types of energy components: 1. kinetic energy 2. the energy that is represented by pressure volume flowing through the surface 3. heat energy. All these components were calculated separately at both inlet and outlet. If we add them together they supposed to give zero in theoretical case. As you see since I have not used enthalpy for my calculations that includes the kinetic energy component, I had to calculate it separately as Pf (energy of the flux) at the inlet and outlet.

You suggest that I should calculate the total energy flowing through the inlet and outlet by doing a MassFlowInt of Total Enthalpy. I did it, here is the result: Pein= 1682.31 [W] Peout= -1716.8 [W] if you compare these results with the values from solver report Boundary: in 1.6760E+03 Boundary: out -1.6760E+03, they do not match (and make no sense at all for me).

In your second point you say that I have left out the viscous work from my analysis. As far as I know in incompressible laminar liquid flow the viscous friction can generate two kinds of energy. One is the heat energy developed by the viscous friction between the fluid layers moving relative to each other, and the other is mechanical work that the wall shear forces can do if they are allowed to move the wall in the direction of the flow. Since the walls of the nozzle are stationary the second term is not present, the wall shear forces do not perform any work. Thus, only the heat component is present that has been taken into account in my analysis with the formula: Phout = -massFlowInt((Temperature -Tain)*2400 [J kg^-1 K^-1])@out This power is generated by the viscous friction.

You say further that "A first sanity check would be to turn off viscous work and see that the balances make sense then." I have run the solver again with the viscous work term not included and got the following results: I have changed the Auto timescale to physical 1 s at time step 11, and the solved finished within 32 minutes by reaching below 1E-6 RMS residual in time step 81:

In CFX Post the 3 energy components (the kinetic, heat, and pressure) are calculated for in and out: Pfin= -63.528 [W], Pfout= 1014.78 [W], Phin= -0.0328285 [W], Phout= 757.569 [W], Ppin=-1689.95 [W], Ppout=5.55364 [W], Ptot= 24.4261 [W] After adding up the components we get that Pin= -1753.51 [W], power enters the nozzle and Pout=1777.9 [W] leaves the nozzle. These values are still much different from 1676.6 indicated in the solver report as seen above.

For the sake of easy comparition let me reiterate here again the original values got from the first run with viscous work term included: Pfin= -63.528 [W]; Pfout=1014.78 [W]; Phin=0.0275742 [W]; Phout=792.705 [W]; Ppin= -1689.94 [W]; Ppout=5.55357 [W]. Ptot=59.56 [W]

You can see that only the power converted into heat has changed, and honestly I do not understand why did it change as it did. To understand better what is meant by viscous work term I looked into the manual, but found only some meager explanations:

Reference Guide CCL Content P189 Include Viscous Work Term : Indicates whether the shear stress term dot velocity term is included in the energy equation.

Solver Modelling Basic Capabilities Modelling P7 A heat transfer model is used to predict the temperature throughout the flow. Heat transfer by conduction, convection, and (where appropriate) turbulent mixing and viscous work are included.

As just said, according to my present understanding, from viscosity we can get heat and mechanical work if the walls are moving. Since the walls are not moving we can talk here only about converting energy into heat by viscous friction. Now, if we do not include this component (viscous work term) then according to my understanding either no heat increase should be indicated (if its meaning is heat generation) or the generated heat should be the same in both cases (if the mechanical work of moving the walls is understood under viscous work term). I would appreciate if someone would explain this thing in details.

Michael,

You also suggested to calculate the massFlowInt(H_total) and you can see the results above. Let me try to understand your formulas here a bit better. massFlowInt(H_total) = massFlowInt(h_static+0.5*vel*vel) = assFlowInt(cp*T+0.5*vel*vel)

The first part is clear, it says that the total enthalpy is the sum of static enthalpy and the kinetic energy component. I have calculated the kinetic energy component at the out as

Pfout=-massFlowInt(Velocity ^2/2)@out.

Further you write that the specific static enthalpy is cp*T. If I understand correctly then cp in your notation is the specific heat capacity at constant pressure and T represents the increase of temperature. If this is correct then the component massFlowInt(cp*T) represents only the energy component that is converted into heat (I calculated at out as:

Phout=-massFlowInt((Temperature -Tain)*2400 [J kg^-1 K^-1])@out

and does not include the work of pressure volume that passes through the boundaries. In Phyl's formula: h = h_ref + Cp*(T-Tref) h_ref seems to represent the energy pressure volume flow component.

Looking into the manual found the following definition:

Solver Theory Basic Solver Capability Theory P7 Specific Static Enthalpy is a measure of the energy contained in a fluid per unit mass. enthalpy is defined in terms of the internal energy of a fluid and the fluid state: hstat =ustat+pstat/ï²stat

If I understand correctly then here ustat supposed to be the heat energy content of the flow and pstat/ï²stat is the energy represented by the pressure volume flow, that I calculated at in with the formula:

Ppin = areaInt(Pressure *Velocity u )@in

You can see that I have taken into account all the 3 components that make up the total enthalpy with the 3 kinds of CEL expressions. Now if my formulas for the partial energy components are correct and there are no other components left out, then their sum supposed to correctly represent the energy balance of the system. If there is some additional energy component not included in my CEL expressions then what is it (what CEL expression could calculate it?). If the CEL expressions are incorrect then what is wrong with them, and what is the correct way of getting them from basic variables of temperature, pressure, and velocity. The enthalpy supposed to be calculated from the same 3 variables, so let us do it ourselves (not using enthalpy) in order to understand better what we are dealing with.

In my previous message I did not say that I expect the sign of the flow through the inlet to be negative. I am aware that a flow into a closed surface is considered positive and flow out of it negative. I just wanted to explain, that In my CEL functions I have used opposite signs, not as this convention dictates, and also gave the reason and logic for this modification. If one keeps this in mind then it has no significance. If you like I can turn around the signs to comply with the convention.

You say: "Another thing you have misunderstood is the concept of accuracy vs. conservation. These are two compeletly different things". I understand that the real accuracy of the solution depends on the mesh quality, convergence, and the accuracy of the physical model built into the software. As said earlier, now I do not care about the accuracy of the physical model, only about mathematical accuracy that the residuals supposed to indicate. Better to say that for now I care only about the accuracy of the derived variables (like pressure, velocity, temperatire in the res. file) that determine the accuracy of the partial energy equations, that in turn determines the residual that remains if these components are summed up. This is why I did not rely on the residuals indicated in the solver report, but wanted to get it by using the variables from the results file through CFX Post.

Under point 1, you are correct to include Pfin (kinetic energy advection). It is Ppin (flow work) which is incorrect because this effect is automatically included in the enthalpy. I suppose you can split enthalpy into internal energy and pv if you really want to, but you'll need to be careful to use reference conditions for internal energy which is consistent with CFX's reference conditions for enthalpy.

Viscous work refers to the work done on the fluid by viscous forces. Refer to Eq. 5 in the 5.7 solver modelling documentation. For an independent source, see the 1960 version of Bird, Stewart, and Lightfoot (1960), Eq. 10.1-9 (last term). Viscous work is nonzero at inlets and outlets as well as moving walls because both stress and velocity are nonzero.

Yes, the Ppin flow work term is indeed included into the static enthalpy. However, it is irrelevant whether it is included or not, because I have absolutely not used the enthalpy in my original analysis to check the energy balance. I have also given the reason why do I want to avoid using enthalpy at all, and why do I want to calculate each energy component separately from basic variables (velocity, pressure, and temperature); not through enthalpy.

You have seen also that I have calculated the MassFlowInt of Total Enthalpy and it still did not show values that would make much sense (specially at the outlet). So before starting to bother with enthalpy I would suggest, let us concentrate on each single energy component separately and calculate them using the basic variables by CEL expressions. I have still not found any error in my formulas, so I have no reason to throw them away. If you find any error in them, please let me know what is it, and show the correct formula for the same energy component (without using enthalpy).

You say: "Viscous work refers to the work done on the fluid by viscous forces." If I am not mistaken here under "work" we understand mechanical work that represents movement under force ie. W=F*s where F is the force that is pushing the object and s is the length of displacement of the object under the effect of force.

The case of a moving no-slip wall is clear: the object to be moved is the wall and the force that is moving it is the viscous friction of fluid against the wall. So here we talk about work done by viscous forces on the wall, or work that the fluid performs when moving the wall. The power of this work term can be calculated if we integrate the Wall Shear stress * elementary Surface Area on the moving wall* Velocity in the direction of the wall shear stress. This has a solid physical meaning.

In the case of our nozzle there is no moving wall. The boundary surface of inlet and outlet is imaginary, and it does not move, there is no object that could be moved by the fluid on these surfaces. So I can not conceive any physical meaning of " work done on the fluid by viscous forces". The viscous forces are not separate and independent from the fluid. So how can the fluid perform work on itself by moving itself?

I have checked the formula and the relevant explanation in your reference: Solver Theory Basic Solver Capability Theory P22 If viscous work is significant, then an additional term is added to the right hand side of the above energy equation to account for the effect of viscous shear, and the energy equation becomes,"

Unfortunately this explanation and the differential equation still do not clarify the physical meaning of this term any better. According to your suggestions I have run the solver for the same problem first with viscous work term included and later without including it. The only significant difference between the partial energy components of these two run was that when the viscous work term was included we have got Phout=792.705 [W] power converted into heat; while when the viscous work term was not included we have got only Phout= 757.569 [W]. That means the viscous work term seems to be heat energy component dissipating energy into heat, and not performing mechanical work. This is the only sensible explanation for the indicated heat (temperature) increase when viscous work term is included. So now we have fallen into a nice confusion, because this additional heat energy increase at the output does not seem to diminish the other energy components at the outlet, as one would expect based on the law of energy conservation. So I suppose this viscous work term must be some imaginary energy component that has no real physical meaning, and consequently it is nonsense to include and use it when we want to model the real physical processes.

Unfortunately I have no access to the book issued in 1960 that you referred to, so I would be glad if you could quote the relevant definitions found there or explain with your own words what is written about this issue.

After reading your questions and several of your answers to others, I would like to express my views. By the way, I am going to use Enthalpy despite your dislike to it.

From first law for an open system in steady state with one inlet and one outlet (Introduction to Thermodynamics: Classical and Statistical, by R. Sonntag and G. Van Wylen, any edition any year, usually by John Wiley):

So, at this point everybody's equations are correct. Everybody happy I hope.. The problem is how we compute the integrals and what reference conditions to use. You seem to compute your own reference conditios which may not be the ones CFX-5.7 is using.

The CFX-5 flow solver uses enthalpy for multiple physical reasons beyond this discussion. The effective flow reported in the outfile is equivalent to a mass flow integral of the term in the equation above, that is

1 - Using area integral for the P*V term, when it should be massFlowInt (P/density)

2 - You must use the hybrid values (read manual for detailed explanation) since the solver uses the boundary conditions for this calculations as anyone would expect.

3 - To avoid CFX-Post from using the conservative value for the inlet calculations, I would hardcode 10 [m s^-1] into your CEL expressions. Perhaps calling CFX helpdesk will be better here (have you not done so already, anyways?)

Now that the quantities are clear with/without enthalpy. You can repeat Neale/Phil's suggestion running without the viscous work term. The energy flows should be within round off between CFX-Post and solver output file.

Now the mysterous (to you at least) viscous work term energy flow you cannot account for at the outlet. This term as indicated in previous post and the documentation is scalar product of Shear_stress * Velocity vector. At your outlet, the velocity vector is not completely normal to the outlet (do a vector plot to verify for yourself); therefore, the term is not zero and the flow solver energy flow reflects that. Can it be calculated in CFX-Post? Not easily I believe, or not at all.

How come then I believe this difference is not something else? Well, add an extension to your domain with uniform cross section to allow the flow to develop and be fully normal at the outlet. Repeat your energy flow check, and you will see that tau*v is zero at the outlet and the energy flows should be round off compared to the CFX-5.x flow solver regardless of using massFlowInt(Total Enthalpy) or your split version of energy components.

I hope this clear your interpretation of the CFX-5.7 flow solver results which should have been explained to you in the training course. If it does not, my best advice is to contact CFX helpdesk directly and someone there must be able to help you.

Juan Carlos

PS. I am amazed you do not know about Phil's reference. It is perhaps the most published book in transport phenomena, Edition 62 just came out (kind of the VW Beetle or the Energizer bunny) . See http://www.engr.wisc.edu/che/bsl.html

Thank you for the excellent demonstration and very clear formulas. Your analysis of the problem and suggestions how to find the solution are very eye-opener. I am also glad that all our formulas are correct, and include all the 3 energy flow components (shown for the inlet):

Heat Energy:

In my original formulation:

Phin = massFlowInt((Temperature -Tain)*2400 [J kg^-1 K^-1])@in

Where Cp=2400 [J kg^-1 K^-1], and Tain = areaAve(Temperature )@in supposed to represent a conservative reference temperature measured at the inlet. Naturally this can be as well substituted with the hybrid reference value of 298.15 K at the inlet, and I have already tried both cases earlier, but it makes not much difference.

Your corresponding formula of this heat energy component is:

m_dot * Cp * (T - Tref) which is basically the same.

Flow Work:

Ppin = -areaInt(Pressure *Velocity u )@in

In your notation:

V_in * area * (P - Pref) where Pref=0 in our case, is in good agreement with my version if one takes into account that (V_in * area) is a scalar product of two vectors, that means we should multiply the area only with the velocity component that is normal to the surface (in my notation this is Velocity u).

Kinetic Energy:

Pfin = -massFlowInt(Velocity ^2/2)@in

Your version:

0.5 * m_dot * V_in^2 which is basically again the same.

The only component in your formula that I have not calculated is the reference energy content term: m_dot * h_ref, because I have taken it to be zero. We are interested only in the change of energy levels and not in the absolute energy content of the flow. Since this reference term is the same at the inlet and outlet we can simply take it to be zero and still get correct energy balance (ie. Valid energy conservation) when summing up all the components for the whole domain.

After agreeing about the basic components let's proceed to your suggestions. You write: "The effective flow reported in the outfile is equivalent to a mass flow integral of the term in the equation above, that is

1 - Using area integral for the P*V term, when it should be massFlowInt (P/density)"

Although theoretically massFlowInt(P/density)@in supposed to give the same result as Ppin=areaInt(Pressure*Velocity u )@in and as your formula of V_in*area*(P-Pref) when integrated over the inlet surface, the discrepancy between the real and originally declared boundary condition values causes significant difference in this case.

In CFX Post we can calculate the massflow through the inlet with the following CEL expressions (or their equivalents in the calculator) when the Density=1262 [kg m^-3] and the surface Area of the inlet=1E-4 [m^2] and the declared boundary condition is vin=10 [m s^-1]:

1. massFlow()@in= 1.262 [kg s^-1]

2. massFlowInt(1)@in= 1.262 [kg s^-1]

3. areaInt(Velocity*Density)@in= 1.25662 [kg s^-1]

4. areaInt(Velocity u*Density)@in= 1.24663 [kg s^-1]

The results are the same independent whether they are calculated with hybrid or conservative values in the calculator. If the velocity through the inlet boundary would be really normal to the surface everywhere and exactly 10 m/s everywhere, then all the 4 expressions would give the same 1.262 [kg s^-1] result. Obviously the massFlowInt calculates with these idealistic (and unrealistic) declared conditions.

But if you take a look at the velocity vector plot on the inlet boundary, then it is visible that there is a significant area near the no slip walls (due to strong viscosity) where the velocity is NOT normal to the inlet surface and its intensity is also much different from the declared 10m/s. It does not matter what do we declare as the boundary condition at the inlet or outlet when defining the setup, if it is unrealistic (non-physical) ie. it can not be realized in reality, then it is incorrect to calculate the massFlowInt as it is calculated: using unrealistic variables.

The expressions 1 and 2 give the most unrealistic value for the massflow because they ignore both the change in intensity and change in direction of the velocity. The expression 3 takes into account the more realistic changing intensity of velocity near the walls, but still ignores the deviation of it direction from the normal. The expression 4 takes into account both deviations of intensity and direction from the declared and physically impossible boundary condition. So the most accurate value is calculated using areaInt and only the velocity component that is normal to the inlet boundary ie. Velocity u. I think it is obvious that if the velocity is not normal to the inlet surface, then the massflow can not be calculated accurately as massFlowInt does, taking it to be constant and normal to the surface.

If the solver indeed calculates this energy flow component with massFlowInt (P/density) and takes the unrealistic boundary conditions, then I would say, that the fault is in this way of calculation and not in my CEL expression using areaInt and Velocity u. I have calculated the discussed flow work with both approaches. With :

-areaInt(Pressure *Velocity u )@in= -1689.94 [W]

-massFlowInt(Pressure/Density)@in= -1714.86 [W]

The difference is significant 24.92 [W] when compared to the original residual of Ptot=59.56W and this causes about half of the error. You have hit the nail on the head, this helped me a great deal to understand the source of a major error, and also the way to avoid it.

My conclusion is that my way of calculating this energy term through aeaInt is the correct way and the solver builds an error into the calculated variables if it does not takes into account that the realistic variables at the inlet and outlet differ from the unrealistic originally declared boundary conditions. The only way I see at the moment to eliminate this inaccuracy is to run the solver again, but this time using the realistic variables at the inlet and outlet taken from the first run. Depending on the accuracy requirement one might need to repeat this rerun process several times.

You say in next two points:

"2 - You must use the hybrid values (read manual for detailed explanation) since the solver uses the boundary conditions for this calculations as anyone would expect.

3 - To avoid CFX-Post from using the conservative value for the inlet calculations, I would hardcode 10 [m s^-1] into your CEL expressions. Perhaps calling CFX helpdesk will be better here (have you not done so already, anyways?)"

I would not say that anyone would expect that the solver should use physically impossible boundary condition variables for these calculations. I would certainly not expect this, rather that it should use the real velocity values at the boundaries as it is refined in each timestep. I would expect that it should model the physical reality as close as possible, and this is not the closest possible if it sticks to an unrealistic boundary condition.

I have tried to implement your last 2 suggestions, that have basically the same aim to use the originally declared boundary values instead of conservative values as follows:

The so calculated Ptot1= 34.7921 [W] (34.79/797.97*100=4.36%) is closer to zero than the original Ptot=59.56W but I think this supposed to be still more accurate, if everything would be OK.

The next step was to make a new run with an extension at the outlet and by implementing your other suggestions to see whether we would really get the same residual with our formulas as indicated in the solver report within round off errors.

The geometry:

As suggested by Juan I have added a straight 30mm long (10 times the height of the outlet) extension in x direction having uniform cross section (the same as the original nozzle outlet) at the outlet of the original hexahedral nozzle (dimensions of the outlet: y=3mm* z=10mm). The one symmetry plane (XY Plane) has included one side face of the extension as symmetry face.

The mesh:

The meshing parameters in CFX-Mesh: Default body spacing 0.7mm, Default face spacing 0.7mm. A user specified face spacing includes 3 faces: the new outlet now at the end of the extension and the two narrow side faces touching the outlet face (3mm high, and extending only to the end of the original nozzle end). This face spacing has a constant edge length of 0.7mm. One might think that this will not improve the mesh since the default body and face spacing have the same 0.7mm edge length. But, there is indeed difference between the mesh quality when using this user face spacing, or without using it.

The 3 composite walls (the symmetry plane not included) have inflated boundaries, Number of inflated layers 5, Expansion factor 1.2, total thickness, maximum thickness 1mm. No edge and no surface proximity. The not mentioned parameters are left as default. This produces a mesh of :

Total Number of Elements = 382153

Total Number of Tetrahedrons = 328813

Total Number of Prisms = 53340

CFX Pre:

The setup is analogous to the previous two runs (described in my first and previous message), except that the domain now includes the newly added extension; and the 3 walls parts on this extension are Free Slip Walls (unlike the 3 walls on the nozzle that are No Slip Walls). My intention with using free slip walls on the extension is to minimize the unnecessary head loss caused by the extension, because the only purpose of this extension is to make the velocity parallel with the x axis and normal to the outlet surface.

The RESULTS:

I have run the solver first not including Viscous Work Term, as suggested. Starting with Aggressive Auto Time scale and changing to Physical Time scale of 2 s at time step 7, the solver finished after 37 minutes with residuals:

When calculating the energy flow components using the last shown modified CEL expressions with hybrid values and "hardcored" inlet velocity and adding them up we get the energy flow into the inlet Pin1= -1555.21 [W] whereas the solver report indicates 1453.8W. The same for the outlet Pout1= 1556.21 [W] instead of 1453.8W. I suppose that the difference 1555.21-1453.8=101.41 W is some kind of constant reference energy flow term, that I have not included, and I have no idea how does the solver chooses exactly this value.

The residual when the Pin1 and Pout1 is added together is Ptot1=1.00006 [W] (1 [W]/842.149 [W]*100%=0.12%) this is still much more than the indicated 3.6621E-04, but it is already much better than in previous examples, what more it is already an excellent accuracy. If Ptot is calculated with my original formulas we get Ptot=20.1108 [W]. Almost all this difference comes from the difference in calculation of Ppin with areaInt and conservative values: Ppin=-1470.48 [W] and with mssFlowInt hybrid values: Ppin1=-1492.11 [W]; Ppin-Ppin1=21.63. I someone wants to perform the same CEL calculations using once hybrid and once conservative values might run into difficulties, due to some bugs in the CFX Post, but let us discuss these things under another heading later.

Next I have run the solver for the same problem but with included Viscous Work Term. Starting with Aggressive Auto Time scale and changing to Physical Time scale of 2 s at time step 11, and later to 20 s at step 52 the solver finished after 1 hours 45 minutes with residuals:

Pin1=-1550.02 [W], Pout1=1588.58 [W], Ptot1=33.3706 [W] this is already a significant error. When calculated with the original CEL expressions Ptot=49.3151 [W] that is even greater. The increased error comes from the increase of the fluid's temperature and heat energy. The modified (hybrid) heat power in the previous run excluding viscous work term was Phout1=842.149 [W], and the same in this run including viscous work term Phout1=874.496 [W]. The difference is 32.35W that almost fully covers the increased Ptot1=33.3706 [W] compared to the previous run.

It is obvious that when including viscous work term then the temperature of the fluid increases, so this term supposed to be a dissipative term that increases heat generation, but the so generated heat energy is not subtracted from the other energy components at the output to conserve energy. This leads to seeming excess energy and bad energy conservation.

You say about viscous work term: " This term as indicated in previous post and the documentation is scalar product of Shear_stress * Velocity vector. At your outlet, the velocity vector is not completely normal to the outlet (do a vector plot to verify for yourself); therefore, the term is not zero and the flow solver energy flow reflects that. Can it be calculated in CFX-Post? Not easily I believe, or not at all. How come then I believe this difference is not something else? Well, add an extension to your domain with uniform cross section to allow the flow to develop and be fully normal at the outlet. Repeat your energy flow check, and you will see that tau*v is zero at the outlet and the energy flows should be round off compared to the CFX-5.x flow solver regardless of using massFlowInt(Total Enthalpy) or your split version of energy components." I understand that the velocity vector at the outlet was not normal in the original 2 runs and this might have caused a nonzero Shear_stress * Velocity vector work term. But in the last 2 runs (done as you suggested) using an extension, the velocity vector at the outlet is normal to the outlet boundary surface with very good accuracy. That means this tau*v component supposed to have disappeared to allow good energy conservation even if the viscous work term is included, as you have predicted. But this did not happen, and the significant difference in accuracy between the cases run with- and without viscous work term does not comes from a mechanical work component that intends to move something parallel to the outlet surface, but it comes from temperature increase.

I have already indicated in my previous post that the viscous work term contributes an additional heat energy portion to the normal heat generation, and increases the temperature of the fluid, while not diminishing the other energy components to keep the conservation valid. So your explanation of the nature and effect of this viscous work term at the outlest seems to be incorrect, and does not explains the heat increase, and the consequent energy imbalance.

Why is the temperature increased by this term, how is this increase calculated by the solver, and what really does this included viscous work term in the solver is still not understood by me, and we can say it is still a mystery. But the most important question to be answered is: whether the effect of this term on the solution (additional temperature increase) really exists in reality, or it is only an imaginary term that physically does not exist? If we want to model the physical reality should we include it or not. When should we include and when not, and why?

So we have got much closer to the truth thanks to all of you who contributed to this discussion, but this viscous work term is still not demystified

By the way for the last two runs: MassFlowInt(Total enthalpy)@in= 1455.85 [W] viscous work term NOT included MassFlowInt(Total enthalpy)@out= -1454.86 [W] viscous work term NOT included and MassFlowInt(Total enthalpy)@in= 1459.02 [W] viscous work term included MassFlowInt(Total enthalpy)@out= -1487.21 [W] viscous work term included Please compare these values with the valuse indicated in the solver report (see above) and with the values derived with CEL expressions and hopefully you will realize why am I so much reluctant to use Enthalpy for any practical purpose in CFX 5 (because it make no sense at all).

Finally you say that these things supposed to be explained to me in the training course, and I should contact CFX helpdesk. Well, I have never attended any CFD training course but learning myself (and with your help). I can not use CFX helpdesk, nor their private website because I only use this program (as a favor), when it is not occupied (not the owner of the program) thus have no access to any support.

And I do not know the quoted famous bible of transport phenomena, because I have studied electrical engineering (not mechanical engineering and CFD). I still "dare" to learn and use the software because we have got sufficient education in mechanics and physics to not scare me away, and I need this knowledge and the program for my research work.

Under this subject we have discussed two sources of inaccuracy in the derived solution.

1. The massflow and velocity are treated independently in the solver, ie. the massflow is not calculated by the solver using the velocity components stored at the mesh nodes (that are saved to the results file). The massflow is calculated from the initial boundary condition at the inlet no matter if this condition can be realized in physical reality or not, and then kept constant through the whole domain to get the same massflow at the outlet too.

This takes place for example when a liquid of high viscosity (ex. Glycerol) is used and the walls near the inlet are no slip walls and using a constant, relatively high inlet velocity normal to the inlet surface. It is evident that the velocity near the no slip wall can not be the same as near the center of the inlet (but much less) and consequently its direction will also have to point slightly towards the center. In such case the initial condition of constant high velocity normal to the inlet is physically impossible, and the massflow that is calculated from such impossible condition is also inaccurate. The massflow that is closer to reality should be calculated from conservative velocity components.

While the solver very intelligently changes the conservative values (and direction) of the inlet velocity to obtain a physically possible solution, the massflow is not updated accordingly but kept constant as initially. Since the velocities have changed but the massflow not modified accordingly, but both variables are used to calculate the equations, an unnecessary error is built into the solution that will not disappear by refining the mesh and by obtaining very small indicated residuals.

The simplest and best way to avoid this in-built source of error is to insure by using appropriate geometry and walls, that the initial boundary condition at the inlet (and possibly also at the outlet) should be physically possible, so that the massflow calculated from this condition remains correct even after convergence. This is accomplished so that an extension with constant cross section and FREE SLIP walls is added in front of the real inlet with no slip walls (and preferably the same at the outlet too). This way high accuracy solution is possible (considering the energy flow balance) when the walls are stationary and the option "Include Viscous Work Term" is switched OFF, but I suspect that the temperature distribution is not correct.

The reason for this suspicion is that when the Viscous Work Term is included, we get a completely different temperature distribution in the solution compared to the case when the term is not included. This difference is obvious from temperature contour plots for both cases. The heat generation and temperature distribution seems to be closer to reality when the Viscous Work Term is included.

2. The second source of inaccuracy comes from including the Viscous Work Term. If this term is not included, then probably the solver does not calculates separately the generated heat using the velocity gradients, but simply takes it to be the necessary amount of energy flow that is required to conserve energy for each finite volume. Consequently the it is natural that we obtain a very accurate energy flow balance, and valid energy conservation.

But when at least one of the no slip walls (or faces) is moving and the Viscous Work Component is not included, then the work component performed by the liquid on this wall is completely ignored by the solver when the energy balance is calculated (taken as if they were stationary). From this fact and from the above described way of calculating the dissipated heat energy and the temperature distribution, we get completely wrong energy balance in the results file. The velocities and pressures are correct but the temperature of the liquid is wrong.

Therefore if there is a moving wall, and we want to know the real temperature and dissipated heat, then we must Include Viscous Work Term. When this is included, then the solver most probably calculates the dissipated heat and the resulting temperature change from the velocity gradients and viscosity for each finite volume, and consequently obtains a more realistic temperature distribution. The energy flow balance is correct this way with moving walls, but there is about 1% error that is almost impossible to eliminate. This inaccuracy naturally occurs due to the fact that the temperature change is really calculated separately from velocity gradients as it should be, and the small errors in the interpolated velocity values (and velocity gradient values) produce an even higher error when used in the energy equations.

The value of this error strongly depends on the quality of the inflated mesh near the walls where the velocity gradient is the highest. The mentioned 1% has been achieved using 10 inflated layers with the thinnest layer of 0.05mm and 1.1 expansion factor. For coarser mesh we get greater error (about 3.5% when using 5 layers, max. thickness 1mm and 1.2 expansion factor). Although in the energy flow balance calculations it seems that the solution has greater error when the Viscous Work Term is included with all stationary walls, I am convinced that this solution and heat distribution is closer to physical reality and consequently more accurate.

In order to check whether the solver really calculates the temperature change and heat dissipation form the velocity gradients and viscosity I have tried to find a formula and a way to calculate the same using CEL expressions from the results file in CFX Post. This issue will be discussed under a new heading "VISCOUS heat dissipation". I will appreciate your contribution to the clarification of this issue.

> Why is the temperature increased by this term, how is this increase calculated by the solver, and what really does this included viscous work term in the solver is still not understood by me, and we can say it is still a

The viscous work term in the Total Energy equation is obtained simply by dotting the viscous stress term in mometum by velocity = viscous work term. It adds energy so the flow so will act to increase the total temperature.

The rule and basic properties of "Include Viscous Work Term" option (and when and why should we use it) has been clarified under "Summary and Solution" subtitle in this discussion, and continued under the new heading "VISCOUS heat dissipation". So now we can say that this term has been demystified, and the only thing that remains is to find a correct formula for the calculation of viscous heat dissipation that uses velocity gradients and dynamic viscosity from the results file. This is the subject of the new heading.

For your request here is the CCL for glycerol (you might need to reformat it if you want to load it directly into CFX 5 because the Post Editor of this forum alters the text formatting: