Two-equation Turbulence Model Similarity Solution of the Axisymmetric Fluid

Turbulent submerged axisymmetric jets are of interest in many engineering applications. They are formed by outflow of fluid from an axisymmetric nozzle into an unbounded space containing a similar or same fluid at rest. The distinguishing feature of submerged jets (as opposed to free jets) is that the external fluid is entrained into the jet. As a consequence of the entrainment, the jet diameter gradually increases in the axial direction – while its velocity decreases (for momentum preservation). The subject of the present investigation is steady jet flow, with parameters not dependent upon time. The usual simplifying assumptions are adopted: the external fluid is the same as the jet fluid and the nozzle exit is circular. The interest is limited to the fully developed jet – neglecting the development that takes place immediately downstream from the nozzle exit. Also neglected in the present approach (as is common in other analogous solutions – e.g., [2], [3], [4], [5], ... etc.) is the intermittence effect on the jet boundaries, where a purely turbulent flow regime alternates with laminar flow of the entrained external fluid. This, of course, influences the conditions on the extreme ends of the evaluated profiles away from the jet axis, where, however, the solution is of less practical interest, and is in any case difficult to verify experimentally. The fully developed jet possesses an important property of similarity: Profiles of various quantities – such as the velocity profile shown in the example Fig. 1 – evaluated at various downstream distances X1 are of mutually similar shape. This similarity property permits us to convert the governing partial differential equations into ordinary differential equations (which, of course, are much easier to solve). The similarity solutions based upon this idea are standard in the case of the laminar jet [1]. In most engineering applications of jet flows, however, velocities and dimensions are sufficiently large for the flow to be turbulent. Turbulence causes the governing equations to be much more complex – so much that there is a general belief that the only way to solve turbulent jet flows is numerical computation based upon the grid method. Such numerical results, of course, are only valid for a particular set of boundary conditions, and do not provide a universal solution of general validity, which is the object of interest here.


Introduction
Turbulent submerged axisymmetric jets are of interest in many engineering applications.They are formed by outflow of fluid from an axisymmetric nozzle into an unbounded space containing a similar or same fluid at rest.The distinguishing feature of submerged jets (as opposed to free jets) is that the external fluid is entrained into the jet.As a consequence of the entrainment, the jet diameter gradually increases in the axial direction -while its velocity decreases (for momentum preservation).The subject of the present investigation is steady jet flow, with parameters not dependent upon time.The usual simplifying assumptions are adopted: the external fluid is the same as the jet fluid and the nozzle exit is circular.The interest is limited to the fully developed jet -neglecting the development that takes place immediately downstream from the nozzle exit.Also neglected in the present approach (as is common in other analogous solutions -e.g., [2], [3], [4], [5], … etc.) is the intermittence effect on the jet boundaries, where a purely turbulent flow regime alternates with laminar flow of the entrained external fluid.This, of course, influences the conditions on the extreme ends of the evaluated profiles away from the jet axis, where, however, the solution is of less practical interest, and is in any case difficult to verify experimentally.
The fully developed jet possesses an important property of similarity: Profiles of various quantities -such as the velocity profile shown in the example Fig. 1 -evaluated at various downstream distances X 1 are of mutually similar shape.This similarity property permits us to convert the governing partial differential equations into ordinary differential equations (which, of course, are much easier to solve).The similarity solutions based upon this idea are standard in the case of the laminar jet [1].In most engineering applications of jet flows, however, velocities and dimensions are sufficiently large for the flow to be turbulent.Turbulence causes the governing equations to be much more complex -so much that there is a general belief that the only way to solve turbulent jet flows is numerical computation based upon the grid method.Such numerical results, of course, are only valid for a particular set of boundary conditions, and do not provide a universal solution of general validity, which is the object of interest here.
There is a universal solution of an axisymmetric turbulent jet known since 1926, derived by Tollmien [3], [1], [17].Its deficiency is the use of a simple, algebraic model of turbulence.The basic assumption of this model is that turbulence is everywhere in a state of equilibrium: local turbulence production is assumed to be equal to local turbulence dissipation rate.This is not true in real jets -mainly in the vicinity of the jet axis, where turbulent fluctuations are transported by advective and diffusion effects.As a result, Tollmien's solution is known to disagree with experimental data, in particular in the vicinity of the jet axis.In earlier publications [4] and [5], the present author succeeded in obtaining satisfactory general solution with advanced models of turbulence in the related problem of the plane turbulent jet, using one-equation [4] and two-equation [5] turbulence models.These are capable of taking into account the spatial transport of turbulence [1].The axisymmetric jet, which is of interest here, was recently also solved with the one-equation model [2].While this model can account for both advective and diffusive transport, its fails to predict the distribution of the characteristic length scale of turbulence.This was circumvented in [2] by adopting Tollmien's [3] assumption of the length scale being constant across each section of the jet.The results may be described as rather successful, when compared with experimental data -indicating that Tollmien's assumption, in spite of its simplicity, is a good approximation to real conditions.Nevertheless, a better approach, not requiring any preliminary assumptions, has been desirable also for the present axisymmetric jet problem.Such a solution is described in this paper.It uses the two-equation model of turbulence that makes it possible to evaluate, just from basic principles, the spatial distributions of all parameters of turbulence, including its length scale.

The equation of the flowfield
The governing equation of the axisymmetric jet flowfield problem is the streamwise specific momentum transport balance in the form valid for cylindrical geometry.Jets are slender enough objects for Prandtl's version of diffusive transport effect (neglect of longitudinal diffusion) to be applicable.It is assumed that turbulent eddies move in a stochastic manner, without preference for any direction -so that the turbulence, being isotropic, may be characterised by a scalar quantity n t , turbulent viscosity.The longitudinal co-ordinate X 1 (Fig. 1) is assumed to coincide with the nozzle axis, while the transverse co-ordinate X 2 determines the radial direction measured away from the axis.In notation according to [1] the equation is written as: The right-hand side of Eq. ( 1) represents the effect of radial divergence.The left-hand side represents the spatial transport effect.These are expressed by means of the turbulent left-acting Prandtl transport operator [ ] P n t , which expresses the effects of two transport phenomena: advection and gradient diffusion (with neglected longitudinal diffusion) respectively.It may be decomposed into: (2) where the first term, in the present two-dimensional context equal to Ñw = Ñ + Ñ w w , is a scalar product of two vectors, Ñ and w.It represents the operator of the two mutually orthogonal components of advective transport.The second term, on the other hand, represents the operator of turbulent diffusion (with n t variable in space).
The solution of Equation ( 1) leads to evaluation of the spatial distribution ( ) (4) will take care of the other, transverse velocity component, w 2 .To evade the necessity of testing the results as to the validity of Eq. ( 4), it is expedient to transform Eq. ( 1) into an equivalent equation for the spatial distribution of the stream function Y: w X so that all solutions obey Eq. ( 4) identically.
After decomposition into terms corresponding to individual transport effects, Eq. ( 1) may be written as consisting of the following five terms:

Model of turbulence
The turbulent viscosity (or "eddy viscosity") n t , appearing in the governing Equation ( 1) is in, general, equal to In equilibrium turbulence, the two model constants c n and c z are mutually dependent, and according to [1] there is ) This relation need not hold in general turbulence.Nevertheless, it was used as a convenient simplifying assumption in the present computations.
The two quantities e f and ' must be evaluated from the simultaneously solved two equations of the turbulence model.
The first of them is the transport equation for the energy of fluctuations: Again, as in Eq. ( 1), the assumption of a slender shear region permits us here to apply the Prandtl approximation for spatial transport effects, as expressed by the Prandtl transport operator on the left-hand side of the equation.On the right-hand side, [m 2 /s 3 ] is the turbulence energy production rate ( ) while the other terms in Eq. ( 13) are the dissipation rate ' and the radial divergence effect term.
Applying the same decomposition for the Prandtl operator as in Eq. (1) and Eq.(7), it is possible to rewrite Eq. ( 13 The other equation of the turbulence model describes the spatial transport of the specific dissipation rate of turbulent fluctuation energy: where U [m ] is the destruction rate, while the last term on the right hand side, similarly as in Eq. ( 1) and ( 13), represents the effect of radial divergence.Decomposition of the Prandtl operator leads cf.Eq. ( 13) and Eq. ( 15) and [1] to: ' ' In this case, the coefficient of gradient diffusion s ' [m 2 /s] is evaluated from the eddy viscosity and from the assumed knowledge of the corresponding Prandtl number The standard recommended value [1] is Pr ' = 1 3 . .Note that the analogous Prandtl number could also be introduced in Eq. (13).There, however, its assumed value is 1.0 (this is equivalent to assuming identical transport of momentum and fluctuation energy by turbulent motions), in common with standard turbulence modelling practice.
The model constants c U and c Y introduced in Eq. (18) are mutually related in equilibrium near-wall turbulence, as shown in [1], where k is the von Kármán constant.The standard value of c U , evaluated from turbulence decay rate downstream from grids, is c U = 1.44.
Solutions of the transport Equations ( 13) and ( 16) provide the spatial distributions ( ) e f f = X and ( ) ' = f X .Inserting them into Eq.( 11) makes possible evaluation of the spatial distribution of the turbulent viscosity ( ) n t = f X .Inserting them into Eq.( 11) makes possible evaluation of the spatial distribution of the turbulent viscosity and its transverse derivative needed in all three transport equations: In the present solution, in line with previous treatments of jet flows in [2], [4], and [5]

The similarity transformation
Transformation of equations into the similarity form is based upon measurement of transverse distances by a scale which varies in the axial direction in proportion to the local jet diameter.Similarly, velocities are measured by a velocity scale which, at every jet cross section, is proportional to local maximum velocity w m .This decreases in the axial direction according to power-law with a negative exponent a.The initial task is to determine the magnitude of this exponent.The next task is determination of the the exponent in the jet diameter axial growth power-law The discussions presented e.g. in [1] and [2], supported by experimental verifications, indicate that there is and from constancy of axial momentum flow rate [1] there is The similarity transformation is achieved by introducing the following similarity variables: To remove constants that would otherwise appear in the resultant equation, it is proper -according to [3] -to introduce the definition In agreement with usage in [1] as well as [2], [4], and [5], differentiation with respect to this similarity co-ordinate is denoted by a nonindexed left-acting operator (28) d) Transformation for the relative turbulence dissipation rate: Relative stream function f is then required so that it will be possible to transform Eq. ( 5) into a corresponding expression with similarity-transformed variables u = Ñ f h .(30) Inserting the above definitions into Eq.( 5) leads to the following relation between Y and f: (31) In order to perform the transformation of Eq. ( 1), the following conversions need first to be evaluated: Analogous transformation of Eq. ( 13) requires evaluation of terms: To transform Eq. ( 16) then requires evaluation of terms: Finally, it is also necessary to transform the expressions with the eddy viscosity, Eq. ( 11) and Eq.(20): The expressions in the above boxes are inserted into the three transport Equations ( 1), (13), and ( 16) -preferably in their expanded forms -Equations ( 7), (15), and (17).
For reasons of convenience (elimination of higher derivatives in equations) and also for the interesting meaning of the introduced additional variables it is useful to introduce the following auxiliary variables: relative transverse gradient of velocity g u = Ñ, relative gradient of turbulence energy n = Ñ e , relative gradient of dissipation rate z j = Ñ.
Derivatives of the transformed stream function may be eliminated from the transformed equations by substitutions using Eq. ( 30) together with expressions derived from Eq. ( 30) by differentiation: As a result, transformation relations and introduction of auxiliary variables convert the equations into the following set of seven ordinary first-order equations:

Comparison with the plane jet case
It is instructive to compare the above resultant transformed form of the equations, Eq. (32) with the analogous similarity transformed equations obtained in [5] in the plane turbulent jet case, also solved by using the two-equation model.The equations arrived at in [5] may here be re-written as: The two systems Eq.(32) and Eq.(33) appear to be very nearly identical, the only exception being the additional 1/h terms within the first pair of brackets on the right-hand sides of the equations for the axisymmetric case, Eq. (32).Since the original partial differential equations of the problem before transformation differ only in there being the additional radial divergence term in the axisymmetric case, it would seem to be almost without doubt that the additional 1/h term in the transformed equations is the result of similarity transformation of the original radial divergence terms -the last term in Eq. ( 1), (13), and (16).It is surprising that the situation is actually far from being so simple.In spite of the similar structure before as well as after transformation, the seemingly corresponding transformed terms may be of completely different origin and meaning -the above mentioned 1/h term in fact does not stem from the radial divergence term in the equation for nÑ, but it results from the diffusion term in the equation for gÑ!This completely different origin of what seem to be corresponding and perfectly analogous terms is a very strange feature of the solved system of equations.

Boundary conditions
The system (Eq.( 32)) of simultaneous equations is subject to the following obvious boundary conditions on the jet axis, at h = 0: … for zero velocity profile slope n ax = 0 … for zero e f profile slope z ax = 0 … for zero e profile slope zero values of the nondimensional gradients due to profile symmetry.
These are, unfortunately, only five conditions instead of the seven required for solution of the seven Equations (32).This is a typical example of a problem of split boundary conditions: the sixth and seventh values, e ax and j ax at h = 0 need be adjusted so as to fulfill the required values ' = 0 and j = 0 at the other end of the integration region, at h ® ¥.
The problem of finding the proper values e ax and j ax at h = 0 is not easy.Since the unknown valeus are two and the integrated equations are extremely sensitive even to minute changes, the use of the usual "shooting method" is hardly practicable.Fortunately, its is possible to base the starting values upon the results of the earlier one-equation model solution [2].It may be useful to recall that the solution found in [2] is not unique: the equations are fulfilled for any value of the dimensionless parameter c Z /ks dependent on jet diameter growth factor s and mixing length factor k. The proper value of the parameter was evaluated by comparison with experiments.In [2], the value c ks These boundary values Eq. ( 35) and Eq. ( 36) were now used as starting values in the present solution with the two--equation model.

Solution
There is still another problem associated with integration of the Equation system (32).It is caused by the fact that what seems to be the natural starting point of the integration procedure, h = 0 cannot actually be used.The transformed radial co-ordinate h appears in Equations (32) in denominator position, leading to infinite values of some terms there.This is no substantial problem and was easily solved by shifting the starting position a small distance away from the axis and evaluating the shifted off-axis boundary conditions using Taylor series expansion (in fact quite simple, because the most important first derivatives are zero on the axis due to the symmetry).
Standard values of turbulence model constants c Z and c n from Eq. ( 12), as well as c Y corresponding to Pr .
n 144 (37) with von Kármán constant value k = 0.41 -were used.The only deviation from the standard practice was to insert slightly higher Pr .' = 13837 in place of the usual Pr .' = 13(for the latter there being, after all, no physical justification anyway).The higher value leads to better fulfillment of the boundary conditions on the outer boundary (an effect which can also be achieved by adjusting the conditions e ax and j ax ).
Integration of Eq. ( 32) was performed by a simple computer program TES-AXI2EQ , using the standard Runge-Kutta integration method.A complete listing of the program is attached in the Appendix.Fig. 2 shows the integration result: the half-profiles of the transformed axial velocity u as well as its transformed transverse gradient g and Fig. 2: Solution of the governing Eq. ( 1) in terms of the three similarity variables that determine velocities.This solution was obtained using the boundary conditions evaluated in [2] for one-equation solution.
transformed stream function, f (an integral of velocity multiplied by radius).These quantities are plotted in Fig. 2 against co-ordinate h, the similarity-transformed radius.Fig. 3 then presents the analogous half-profiles of the results obtained from simultaneous integration of the transformed relative energy of turbulent fluctuations: there is the relative energy e as well as its (transformed) transverse gradient n, again plotted against the similarity transformed radius h.
In a similar manner, Fig. 4 presents the evaluated half--profiles of the turbulence dissipation rate ', the spatial distribution of which is described by Eq. ( 16).
In the present solution this quantity is computed mainly as a means for evaluating the turbulence length scale using Eq. ( 10).

Comparisons with experiments
The following several figures present comparisons of the profiles resulting from the present solution with experimental data.For plotting the present results in common co-ordinates with experimental and other data, it is useful to use another relative transverse co-ordinate, e.g., the co-ordinate x related  The evaluation of boundary conditions in [2] used just a single number from Trüpel's experiment, its radius growth factor s 0.5 .The agreement in the shape of the velocity profiles is, therefore, a result of good modeling of the flow and its turbulence.This may perhaps be proved by inspecting another comparison with another experiment in Fig. 6, where there is no such relation.
The apparent disagreement visible in this case near the jet outer boundaries is of no importance, for several reasons.First, the measurements the results of which are shown in Fig. 6 were made using a hot wire anemometer.It is a typical property of a hot wire sensor that it lacks the capability to discriminate positive and negative velocity directions.Since near the jet boundary the time-mean velocity is small and the fluctuations are relatively very large, negative values occur but are interpreted by the anemometer as positive values.As a result, the apparent mean velocity data in this region tend to be higher than the actual mean velocity -a trend which is in line with what is seen in Fig. 6.Secondly, the present solution is not intended to be exactly valid on the jet boundary, because -as was already mentioned in the introduction -it does not take into account the effects of intermittence.
Anemometric measurements performed by the present author and his co-workers made possible not only evaluations of time-mean velocity, but also velocity fluctuation.From the data on the latter, it is possible to calculate local energy of turbulent fluctuations e f -under the somewhat simplifying assumption of isotropy, which is necessitated by the only single-channel available anemometer.The experimental profiles of this fluctuation energy, convertedaccording to Eq. ( 28) -to the nondimensional similarity variable for which the transport equation is here solved, are compared in Fig. 7 with the present theoretical solution (corresponding to profiles from Fig. 3).
It is evident that the anisotropy assumption is not restrictive, and the experimental and theoretical profiles are in quite good agreement.The surprising fact, however, is that this agreement is arrived at, after gradual development, only at a very large distance downstream from the nozzle.Note that in the simultaneously made velocity measurements shown in Fig. 6, the fully developed shape of the profile is attained already at a distance several times shorter.The extreme length of turbulence parameter development, equal to around 60 nozzle exit diameters, does not seem to be known or discussed in the available literature, where it is commonly expressed (no doubt on the grounds of comparing only veloc-ity data) that a fully development state is already attained at 10 or perhaps 15 nozzle diameters.
In order to investigate this interesting finding in detail, experiments were performed during which the velocity fluctuation was measured by a hot wire anemometer along the jet axis.The results of these experiments are shown in Fig. 8. Increasing time scale of fluctuation with the downstream distance caused a rather large scatter of data at larger values of X 1 .Nevertheless, the dependence demonstrates quite clearly the asymptotic character of the approach to the value Eq. ( 35), found by the similarity solutions both from [2] and from the present analysis.

Turbulence length scale
Apart from the above described comparisons with experimental data, it may be of interest to compare the present results with what was obtained earlier in [2] with a less sophisticated, one-equation model of turbulence.Fig. 9 compares the half-profiles of the axial time-mean velocity computed from the present solution (cf.Fig. 2) and from [2].
Both velocity profiles are nearly coincident -to the degree that there seems to be no possibility to decide between them on the basis of any conceivable experiment.
The main difference between the present two-equation model solution and the earlier solution [2], which used the one-equation model, is the closed character of the present approach, not dependent upon outside information or assumptions.The one-equation model required an assumption on the spatial distribution of turbulence length (taken as constant in [2]) and, to achieve unicity, some information about the jet spreading rate.It may be of interest to see what was accomplished by the present two-equation approach in these two respects.
As for the length scale, its distribution across the jet cross section may now be calculated from the present solution, using the results on the distributions of fluctuation energy (Fig. 3) and dissipation rate (Fig. 4), with the use of Eq. ( 10) and of the values of constants Eq. ( 12).In Fig. 10, values obtained this way are shown in the form of their ratio to the constant turbulence length of the earlier solution [2].There was, in [2], while the present 2EQ is given by Eq. ( 10).Using the definitions Eq. ( 28) and Eq. ( 29), the ratio of the two lengths is In Fig. 10, it is immediately apparent that the length scale distribution resulting from the present solution changes only insignificantly across the jet profile.Indeed, an indication of the small changes of length scale across the profiles was visible from the small difference between the two solutions in Fig. 9 and, in fact, already from the similar shapes of the curves in Figs. 3 and 4 (with constant c Z cf.Eq. ( 10)).This insignificant variation provides an explanation of the apparent success of earlier jet solutions, beginning with [3] and including [2], which assumed the length scale to be simply constant across the jet.

The illusory "plane jet/round jet anomaly"
While in the simpler previous solution [2] it was necessary to evaluate the jet diameter growth using information obtained from experimental data, the closed character of the present two-equation solution makes possible direct evaluation of the jet spreading rate.In fact, the information about jet spreading resulting from the present solution made possible one of important results of this paper.It is a clarification of what for a considerable time has been one of the basic problems of turbulence modelling: the unresolved problem of the "round-jet/plane-jet spreading anomaly".According to assertions in the most authoritative literature on the subject -e.g.[11] -computations with the present standard two-equation turbulence models, while predicting plane jets properly, are believed to fail when applied to axi-symmetric jets.In particular, computed spreading rates are said to be between 25 % and 40 % larger than experimental spreading rates -as shown in Fig. 11.
Such a discrepancy is believed to be a manifestation of a fact that there are physical mechanisms participating in the dissipation process, that have so far not been taken into account in usual turbulence models.Attempts to remove the anomaly have led to various suggested alterations to the model equation for the axisymmetric jet without, however, any convincing success.Usually, as is the typical case of Pope's [12] adaptation of the two-equation model, the suggested modifications remove the "round-jet/plane-jet anomaly" only to exchange it for a spreading rate anomaly in other situations.The reported absence of the anomaly in the Wilcox k -w [11] model was at least partly the reason behind its popularity.The comparisons in Fig. 11, however, are not really convincing.Furthermore, recent investigations (Speziale and Abid, 1995) show that the k -w model is hardly perfect, as it has certain difficulties associated with assumed basic physics of turbulence.
The present results, on the other hand, show that what is known as the "round-jet/plane-jet spreading anomaly" is only a fictitious effect, caused by improper numerical values of boundary conditions (and perhaps model constants) used by earlier authors in their earlier computations.In Fig. 12, the velocity profiles from Fig. 11 are compared with the present solution, which is re-plotted using the same transverse co--ordinate (which directly corresponds to the value of the jet diameter growth constant) as in Fig. 11.Obviously, no spreading anomaly is seen, although the present solution shown in Fig. 12 uses standard model constants (with the exception of the small adjustment in the value of Pr ' , which is irrelevant in the present context) as they are used also for the plane jet solution.Fig. 13 shows the error that was probably responsible for the wrong behaviour of two-equation solution in Fig. 11: without changes in model constants but with a different set of boundary conditions e ax and j ax , the curve that should represent the "standard two-equation solution" in [11] is quite well computed by the present solution, using the computer program TES-AXI2EQ as listed in the Appendix.In fact, by slightly decreasing e ax below the value e ax = 0.1213 indicated in Fig. 13, it is possible to achieve almost perfect coincidence of the two curves -the slight mismatch was preferred for presentation in Fig. 13 so that both curves may be individually recognisable.Obviously, the bad reputation of the two-equation model was due to a similar erroneous set of boundary conditions used in some earlier computations.It should be stressed that such an error is extremely difficult to Fig. 11: This illustration is taken from [11].It shows the common belief in the literature on axisymmetric jets, that computation using the standard two-equation turbulence model leads to a wrong spreading rate.The model due to Wilcox [11] is proposed as a remedy.
Fig. 12: The present results re-plotted in the same co-ordinates as used in Fig. 11 show that the spreading anomaly is a fictitious effect [10]: no anomalous spreading is found for the present solution with proper constants and proper boundary conditions identify in usual computations based upon a grid solution of partial differential equations -in particular in view of another of the present findings, the extremely long way towards equilibrium, as shown in Fig. 8 (computations are seldom performed on regions as long as is now known to be necessary).It is an advantage of the similarity transformation approach that all these details may be easily clarified.

Sensitivity analysis
A serious problem with the standard two-equation model of turbulence, which actually forms the background of enigmatic features such as the discussed fictitious anomaly, is its extreme sensitivity to exact values of model constants and boundary conditions [14].
The sensitivity problem becomes really serious in common applications of the model to computations using grid methods, where the problem of boundary conditions often appears in a disguised form and may lead to conclusions that the system of solved equations is unstable or unsolvable.On the other hand, the present similarity approach is particularly suitable for clarifying this question.After the transformation, the boundary condition appears in the form of just a single numerical value to be determined.Using the present solution, it is tractable and in fact quite instructive to study the dependence of the solution errors on changes of the boundary condition values.In the sensitivity analysis, the magnitude of the solution error was judged by the magnitudes of the two end boundary conditions which are to be fulfilled at the end of integration: the values of relative velocity u and relative fluctuation energy e outside the jet boundaries -at h = 0.8 (Fig. 15) where they should be zero.The results of repeated integrations with varied values of boundary conditions at the starting point (on the jet axis) are shown in Fig. 16.Fig. 13: Demonstration of the fact that the improper spreading rate, believed to be an inherent drawback of the two-equation model, is just a result of using improper boundary conditions on the jet axis: with a change in the value of the conditions, the present solution may be changed to match the curve from [11] Fig. 14: Classical sensitivity analysis is concerned with finding the slope of linear (or locally linearised) dependence between the deviations of parameters from the correct value and the resultant error in the solution They indicate that, in the present problem, there occurs the worst conceivable situation known by classical sensitivity analysis: the linearised dependence is here a vertical line with an infinite slope.This means that the solution of the jet with the two-equation model exhibits infinite sensitivity even to the smallest variations of the starting point values.It is only due to the nonlinearity of the solved system -the fact that the dependence is not linear and there is some curvature -that any solution at all is arrived at.Note that some small variation of the starting value (leading, however, to an error of the order of tens percent) is admissible only in one sense: only an increase in e ax or a decrease in j ax is tolerable.Even a small change in the wrong direction causes the whole solution to collapse.
Of course, simultaneous changes of both values in their acceptable directions are possible without collapse of the integration -these are the compensatory variations of [15].Nevertheless, they lead to wrong predictions: the wrong "classical" results from Fig. 11 are a consequence of such compensated deviations from the proper values without collapse of the solution, as demonstrated in Fig. 13. 12 Conclusions 1) A general solution was derived for the axisymmetric submerged jet flow with a two-equation model of turbulence, valid for any conditions (but excluding, as usual, the effects of flow development and intermittence).With the exception of Pr ' -which it is advisable to insert slightly larger than the standard value -the solution uses standard values of model constants.
2) The solution is based upon the similarity transformation approach, whereby the governing partial differential equations are transformed into a set of seven ordinary differential Equations (32).Their integration is straightforward and presents no stability or stiffness problems.
3) The problem has split boundary conditions: of the seven required conditions at the starting point, only five are known.The success then depends critically upon proper insertion of the two unknown boundary conditions at the starting point of the integration -the knowledge obtained in an earlier solution [2] with a simpler, one-equation model was of substantial help here.4) The profiles of velocity and other variables computed from the present solution are in excellent agreement with available experimental data -a new result, however, is that especially some parameters of turbulence require a substantially longer downstream distance from the nozzle to attain a fully developed condition than was previously considered necessary.5) The main advantage of the two-equation model is the internal closedness of the solution which -in contrast to earlier solutions such as [3] or [2] -is not dependent upon a priori assumptions or information from experimental data.Otherwise the resultant profiles of many variables (such as velocity) are almost coincident with earlier one-equation solution results with , indicating that this assumption is a rather good approximation.6) A stain upon the reputation of the two-equation model has been its alleged inability to predict the proper spreading rate of axisymmetric jets.The present analysis shows convincingly that this "spreading rate anomaly" is a fiction: there is nothing wrong with the model, the erroneous results have been merely a consequence of wrongly inserted numerical values of boundary conditions.7) The system of equations describing the turbulent axisymmetric jet with turbulence modelled by the two-equation model is extremely sensitive (in fact, the performed sensitivity analysis shows theoretically infinite sensitivity) to the values of inserted constants and boundary conditions, which must be used with at least four (and preferably five-) digit accuracy to obtain proper results.
The program has only 100 lines -of which a substantial proportion have only an auxiliary role: for example, lines 10-110 serve only to print the identification heading.Lines 120 to 180 contain the definitions of the seven solved equations Eq. ( 32).Lines 190 to 240 are there to draw the grid used for plotting the resultant profiles.Lines 190 to 260 insert numerical values of boundary conditions on the jet axis, i.e., at h = 0, as well as the value of the parameter CKS = c Z / k×s.In the present form, the program runs with values of the solution parameter and fifth boundary condition as evaluated for Trüpel's jet.Adjustment to different experimental data is possible by changes in lines 19 and 200.In line 250, the small stepsize DETA = Dh and in line 260 the starting radial position are inserted.The following lines 270 to 420 serve to insert constants and boundary conditions of the solution.The actual solution begins in line 430.It is performed in a loop (loop parameter I), which closes at line 990.The integration method used is the standard Runge-Kutta fourth-order algorithm, written in the repeated seven lines corresponding to the five similarity variables: F = f, U = u, EPS = e, G = g, N = n, J = j and S = z.Line 790 then advances the solutio along the transverse co-ordinate ETA = h.The solution is terminated if the value of this similarity transverse co-ordinate becomes larger than h = 0.6 in line 800.The remaining lines 810 to 980 plot the results in a diagrammatic form.Among the plotted variables there is also the turbulence production rate PROD = and the ratio of turbulence length scales l.

List of symbols:
c Z [1] turbulence dissipation rate coefficient c n [1] eddy viscosity coefficient c U [1] turbulence dissipation rate generation coefficient c Y [1] turbulence dissipation rate destruction coefficient ] specific energy of turbulent fluctuations f [1] relative stream function g [1] relative transverse gradient of velocity j [1] relative dissipation rate of turbulence j ax [1] relative dissipation rate on jet axis k [1] coefficient of the integral length growth ] destruction rate of turbulence dissipation z [1] relative transverse gradient of turbulence, dissipation rate a [1] exponent of velocity decrease b [1] exponent of the jet diameter growth d [1] transverse dimension of shear region d 0.5 [1] conventional jet diameter e [1] relative energy of fluctuation e ax [1] relative fluctuation energy on jet axis e 0.8 [1] energy boundary condition error ' [m 2 /s 3 ] turbulence dissipation rate k [1] von Kármán constant l [1] ratio of turbulence length scales h [1] similarity transverse co-ordinate h 0.5 [1] value of h at a position where u = 0.5 s ' [m The nomenclature as well as the form of the equations closely follows the usage in textbook [1].

Fig. 1 :
Fig. 1: Schematic representation of the computed jet flow and its typical profile of time-mean axial velocity.The illustration also serves to define the orientation of the used co-ordinate axes X 1 and X 2 .
the axial component w 1 of the time-mean velocity vector ) where w t is the velocity scale of turbulent fluctuations, local maximum velocity w m (Eq.(21)) on the jet axis.Introducing proportionality constant m (needed for dimensional reasons), we can write Relative kinetic energy:Using maximum velocity on the jet axis as the reference, a suitable definition for the transformed similarity variable is adjusting the jet growth factor s 0.5 to the magnitude evaluated for the classical Trüpel (1915) experiment.The value Eq. (34) led to the relative specific energy of turbulent fluctuations on jet axis e ax = 007614 .(35) and to relative specific energy dissipation rate on the axis

Fig. 3 :
Fig. 3: Solution of Equation (13) for turbulent fluctuation energy in similarity transformed coordinates: relative magnitude of fluctuation energy e and its transverse gradient n.Note the extended vertical scale compared with Fig. 2.

Fig. 4 :Fig. 5 :Fig. 5
Fig. 4: Solution of Equation (16) for turbulence dissipation rate in similarity transformed coordinates: again, the transverse gradient of the evaluated quantity is plotted (on a 10-times decreased scale) in the bottom part of the diagram

Fig. 6 :Fig. 7 :
Fig. 6: Another comparison of the velocity profile obtained by the present solution -as shown in Fig. 2 -with experimental velocity profile data measured by a hot wire anemometer in the author's laboratory

ActaFig. 8 :
Fig. 8: Experimentally determined relative values of fluctuation energy on the jet axis as a function of the downstream distance from the nozzle exit

Fig. 9 :Fig. 10 :
Fig. 9: Comparison of the velocity profile computed from the present two-equation model solution (round symbols) with an earlier solution [2] based on the one-equation model of turbulence (continuous line)

Fig. 15 :Fig. 16 :
Fig. 15: Meaning of the solution error in the sensitivity analysis of the present solution: deviations of the solution for some of the basic nondimensional quantities (velocity and kinetic energy of fluctuations) from their correct zero value outside the jet

4 ] 1
] one-equation model turbulence length scale 2EQ [m] two-equation model turbulence length scale m [m 2 /s] coefficient of the maximum velocity n [1] relative transverse gradient of fluctuation energy generation of turbulence dissipation rate w m [m/s] maximum velocity in the profile w t [m/s] characteristic velocity of turbulent motions w [m/s] velocity vector w [m/s] time-mean velocity vector w 1 [m/s] time-mean axial velocity component w 2 [m/s] time-mean transverse (radial) velocity component X