Hi experts! I need your opinion on the following; a colleague says that the k-e turbulent model was calibrated long time ago, when due to computing power the mesh density it was tested for was not as fine as it can be at present. This means that for different mesh densities different solutions will be obtained due to incompatible constants in the k-e model. Is it so? because I think that after certain mesh density the solution should become mesh independent. But he says that it is only true for laminar flows and not true when turbulent model is used.

When using the one layer wall function and k-e model, the solution change a lot with the mesh density. I did that test. It seems that the more fine grid near the wall will give the more friction force which will cause the near wall velocity decrease. Then, I remove the wall function, the mesh sensitivity will better. Good luck!

For many cases you will obtain a mesh-independent solution using the k-epsilon model, and as you refine your mesh, your solution will "converge" to a "correct" result, "correct" meaning a correct solution according to the model, which of course might be wrong. Hence, your colleauge is wrong, there is no inherent incompatibility in the k-epsilon model which makes it impossible to obtain mesh-independent results.

However, there are many cases where you will have problems to obtain this mesh independence. First, like sheng said, you must of course always make sure that if you use wall-functions you should have the first mesh point at around y+ 50 or something like that - if you start resolving the viscous sublayer your wall-laws will not work and you have to use something else (use a low-Re model). There are also many other cases, where the k-epsilon model simply isn't good enough to model the physics of the flow - as you refine your grid, and thereby reduce the numerical viscosity, you might, for instance, suddenly get an unsteady flow which will not converge to a steady solution. The result might be that you will get a solution of poor numerical quality, which then of course will not be mesh-independent. I´m sure that there are many similar examples where the k-epsilon model is not good enough to model the flow and thereby will not give mesh-independent results.

Thanks Jonas for your explanation. If I understood correct; the reason of not getting mesh independent result is case dependent i.e. if k-e model is not good enough to model the physics of the flow it is difficult to get mesh independent solution. My friend says it is same to say that the constants in the k-e model are not tested for very fine meshes of today. Is it correct argument? By the way we are having trouble getting mesh independent solution for an axixsymmetric jet.

(1). I don't have problem with k-epsilon model, high Re model or low Re model. (2). For high Re model, normally 30 points across a boundary layer is needed. For low Re model, you need 60 points across the boundary layer in order to obtain reliable results. (3). To check if a turbulence model is a function of the mesh density (this includes both low Re and high Re models), you can check the terms in the model and the coefficient functions in the model and the wall function implementation in the model, to see if the mesh spacing (delta(y), or y(j+1)-y(j) ) is used as a parameter. (4). The quality of the mesh must also be checked to make sure that it is acceptable. And sometimes convergence is a problem with some models. Unable to obtain mesh independent solution? Never heard of it.

(1). Is your jet initial condition (inlet) properly defined and is independent of the mesh system? (2). To avoid the initial condition problem, move it further upstream into the nozzle or the pipe and compute the flow development inside the nozzle, the flow around the lip, the wake and so forth. For example, to compute the jet-in-a-cross-flow problem, it is essential to move the condition upstream the hole exit into the region on the other side of the wall. (3). The jet problem is normally parabolic and the solution depends mainly on the initial conditions. So you must be very careful about the definition of the initial conditions. (4). If you still have problem and the flow is turbulent, try to model it by setting the viscosity to a constant value. In this way, it is easier to isolate the problem.

Your colleague probably has a typical flow solver where the k-e wall boundary conditions are set such that they only correspond to the log-layer part of the boundary layer. These values typically hold for y+ between 50 and 300. If the first point of the mesh is say at y+ = 60, and you increase the resolution by a factor of four, while insisting that you would continue to use the log-law expressions or wall-functions, you will only make the situation worst. Ofcourse, I am talking about worst compared to experimentally observed results, mathematically you would expect to get precisely this sort of behavior.

The typically imposed conditions for k and epsilon in the log-layer are:

k = u_star**2 / sqrt(C_mu)

epsilon = u_star**4 / (kappa * y_plus * molecular_visc)

By using these values at y+ = 15, you will have your k value a factor of 2 higher, then it ought to be, and your epsilon will be about 20% higher, than it ought to be. The net effect will be that your eddy viscosity will be about a factor of 3 higher, which represents a terrible situation. All of this assumes that you have a resolution independent value of u_star calculated, which may not be the case for coarser resolution. Your colleague may have a slightly more ellaborate model for setting the boundary conditions, but if the inherent assumption is that y+ for the first grid point ought to be in a certain range, you can rest assured that you will get bad results as you successively increase the resolution. In the core region, you will have better accuracy, but these regions will be masked by propagation of effect from the boundary conditions.

There is nothing fancy about the boundary layer equations that stops the code developers from adding in more accurate boundary conditions implementations that hold in all regions, including the viscous sub-layer. One can in fact assume that one has the compete law-of-the-wall specified (say to y+=500) in inner-layer coordinates (u+ vs y+, say using the Spalding expression), and for the boundary conditions assume that the corresponding k+, epsilon+ equations satisfy the above expression on one side, and are zero at the other end of y+=0. These two coupled equations can be numerically solved to give a consistent set of k,e profiles. This is the sort of procedure that I used to come up with the numbers quoted for y+=15. The procedure results show that k+ is nearly a constant [1./sqrt(C_mu)] from about y+=50 and beyond and has a monotonic drop-off to zero in the 0-50 range. The value of epsilon+ varies like 1./(kappa*y_plus) from about y+=30 and beyond, but it does not monotonically drop to zero in the 0-30 range. In fact, there is a peak in the vicinity of y+=10, and then a drop-off in the lower range. There is a tremendous departure on the epsilon+ behavior, from the log-layer in the y+=10-30 range.

The above profiles have been coded in our PowerFlow code, in the model which supports the k-e equations. The results based on such approaches will be presented next month in the ASME conference in Boston (August 1, 1999).

(1). I am getting the feeling that wall function is becoming the black art of turbulence modeling. (2). The problem is related to the exact concept of flow field matching using wall function (functions). It is important to know that there is a region ( a gap) where the solution is replaced by the analytical wall functions. This gap is under the user control, and there is no mesh point inside. The mesh and mesh density does not apply in this region at all. So, the mesh size does not affect the wall function region. (3). The problem occurs because in some formulations (especially the finite-volume approaches), this gap region was also included in the finite volume formulation. I have mentioned the proper method to handle this region a couple of times here, and I am not going to repeat it. (4). So, it is all right to say that wall function is a black art.

The requirement that the first cell be at a specified y+ location is an historic one that has been 'rectified' by evolution of wall functions employed by the more 'advanced' commercial CFD codes.

Robin.

The fact about the grid dependence of k-e model is related to the transport equation of epsilon. In the model by Jones and Launder (1972) (low Reynolds number k-e model) the epsilon is a dissipation function rather than dissipation itself as is considered in High Reynolds number version of k-e model. The dissipation rate itself varies a lot near the wall with respect to cross-stream distance as can be seen from the DNS data for many turbulent flows. Therefore the model using dissipation rate itself would be grid dependent. On the other hand the variation of the dissipation function used by Jones and Launder does not vary much with respect to cross-stream distance, so their model is better than all those models using dissiaption rate itself. Sorry it has become very academic. But I think this is the real cause.

the values of k-e model coefficients are determined according to our knowledge of simple turbulent flows (nothing to do with mesh density testing). For example the decay of isotropic turbulence suggests Ceplison2 has a value of around 1.8 and near wall turbulence dictates that (Ceplison2-Cepsilon1) is fixed. Also Cmu/Sigmaeplison=0.069.

Numerically, it's a different matter - with wall function, you can not refine the near wall cell too much because of y+ requirement. You may remove the problem by adoption two-layer or low reynolds number model.

(1). I like this answer. I think it is getting closer to the reality. (2).The basic problem is the implementation of the Wall Function numerically. (3). The first grid point is on the wall. The second grid point is the matching point, which should be in the law of the wall REGION for the matching to be carried out properly. The region between the wall point and the matching point is not solved separately. (4). The computational domain is defined from the matching point outward. In this computational domain, fine mesh can be used.(which is independent of the spacing between the wall and the first matching point). (5). The basic assumption of the law of the wall is the total shear stress near the wall is constant, which is in the function itself. (6). It is another story when using control volume and TKE to get the shear stress.(this part is the numerical implementations of the wall function, which can be confusing depending on the methods used) Following the leader does not mean that you will get the right answer all the time. (7). To avoid this problem, try to solve the region near the wall also using two-layer, or low Reynolds number models. This is the right suggestion!

(8). Just trying to make it a little bit clear. When one uses the wall function approach, the flow field is divided into two regions ( or two problems) by the user specified matching point location. On one side, you obtain a solution from the analytical law of the wall formulation. On the other side, you carry out the normal CFD analysis, using the matching point as the boundary condition. In this region, you can refine the mesh in any way you want. There is no limitation there.

Well, then I assume that you are using a hybrid method since wall-functions are only valid in the logarithmic region and not in the linear viscous sublayer. My limited experience with hybrid methods is that it is very difficult to get these to work well in general. A simple example - If you go from a poorly resolved wall to a more well resolved wall the upstream region is invetiably not capable of providing the downstream region with enough information to give justice to the good resolution. If you have separations etc. these things can also cause problems.

I'd be very interested in hearing a bit from the commercial side about your experience with these "hybrid methods" - I assume that that is what you are referring to. Could you perhaps elaborate a bit more about what kind of hybrid methods you are using and how they work, if there are any problems with them, etc?