Curvilinear Coordinate System for Mathematical Analysis of Inclined Buoyant Jets Using the Integral Method

The development of a local system of orthogonal curvilinear coordinates, which is appropriate to monitor the flow of an inclined buoyant jet with reference to the basic Cartesian coordinate system is presented. Such a system is necessary for the correct application of the integral method, since the well-known Gaussian profiles should be integrated on the cross-sectional area of inclined buoyant jet, where they are valid. This is the major advantage of the present work compared to all other integral methods using Cartesian coordinate systems. Consequently, the flow and mixing governing partial differential equations (PDE), i.e., continuity, momentum, buoyancy, and/or tracer conservation, are written in the local orthogonal curvilinear coordinate system and, then, the Reynolds substitution regarding mean and fluctuating components of all dependent variables is applied. After averaging with respect to time, the mean flow PDEs are taken, omitting second-order terms, as the dynamic pressure and molecular viscosity, compared to the mean flow and mixing contributions of turbulent terms. The latter are introduced through empirical coefficients.The Boussinesq’s approximation regarding small density differences is taken into consideration. The system of PDEs is closed by assuming known spreading coefficients along with Gaussian similarity profiles. Themethodology is applied in the inclined two-dimensional buoyant jet; thus, PDEs are integrated on the jet cross-sectional area resulting in ordinary differential equations (ODE), which are appropriate to be solved by applying the 4th order Runge-Kutta algorithm coded in either FORTRANor EXCEL.The numerical solution of ODEs, concerning trajectory of the inclined two-dimensional buoyant jet, as well as longitudinal variations of the mean axial velocity, mean concentration,minimum dilution, and entrainment velocity or entrainment coefficient, occurs quickly, saving computer memory and effort. The satisfactory agreement of results with experimental data available in the literature empowers the usefulness of the proposed methodology in inclined buoyant jets.


Introduction
The study of the development of the flow domain of a physical phenomenon is achieved using the equations of continuity, momentum, and conservation of buoyancy or/and tracer mass.For the mathematical description of the problem, a coordinate system must be used.The Cartesian rectangular coordinate system O(, , ) is usually used, where the unit vectors are constant at the space domain.In many physical phenomena fluid particles generally move along a curve.A Cartesian coordinate system may be preferable for use by a numerical simulation of the system of PDEs based on a fine three-dimensional mesh, structured or unstructured, covering the whole flow and mixing field.Such numerical computations of inclined buoyant discharges have been recently made using computational fluid dynamics (CFD) models [1,2].Compared to integral methods with a curvilinear coordinate system, their advantage is the generalized possibility of applications and their disadvantage is the necessary more detailed validation, which is accompanied by longer computational times (from some hours to some days) and much more computer memory and power.Alternatively, depending on the selection of the solution procedure and in case that the specific particles are in curvilinear motion, the above equations may be preferable to be expressed at a curvilinear coordinate system [3].Although, several buoyant jet phenomena can be analyzed by the governing equations written in orthogonal Cartesian or cylindrical coordinate systems [4][5][6][7][8], there exist some specific phenomena that the use of a curvilinear orthogonal coordinate system is Configuration of a curvilinear orthogonal coordinate system adapted on the flow of a buoyant jet in a calm environment.mandatory [9].Such phenomena include buoyant jet flows created by inclined discharges, in which the transverse profiles of the mean axial velocity and mean concentration are approximated with Gaussian profiles.Then, the mathematical description using a curvilinear coordinate system is necessary for the application of the integral method.This method is based on the integration of the partial derivative equations (PDE) of motion and mixing on the transverse cross-section of the main flow field of the inclined buoyant jet to get the corresponding ordinary differential equations (ODE).The latter can be solved with a 4 th order Runge-Kutta algorithm minimizing computational power, memory usage, and time.According to the authors' knowledge, there are no other studies in the literature in the field of inclined buoyant jets, which applied the integral method using a curvilinear coordinate system.
Typically, non-Cartesian coordinate systems are curvilinear [11,12].Their axes are not necessarily mutually perpendicular to each other, and the unit vectors along the axes are not constant in space.The equations that describe the problem are formulated at a local curvilinear coordinate system (O 0 , , , ) with reference to the basic Cartesian system (O, , , ), as it is defined in Figure 1.The most useful curvilinear coordinate systems are the orthogonal (O 0 , , , ), cylindrical (, , ), and spherical (, , ) systems [13][14][15].
The present work develops the mean flow and mixing PDEs of continuity, momentum, and tracer or buoyancy with respect to an orthogonal curvilinear coordinate system.This system is convenient to apply the integral method and derive ODEs, which can be easily solved to calculate the mean flow and mixing properties of buoyant jets.For the present work, an inclined two-dimensional buoyant jet is considered for the application of the mathematical integration to derive ODEs; their solution provides the mean flow and mixing properties, which are compared to experimental data available in the literature.

Development of a Curvilinear Orthogonal Coordinate
System.A basic Cartesian coordinate system (O, , , ) is considered with unit base vectors i, j, k corresponding to x, y, z.The buoyant jet source (either round of diameter D or slot of width D) is assumed to be at the origin O.The buoyant jet is discharged from its source, which is located at the origin O, with an initial velocity w 0 and inclination  0 to the horizontal plane xy into a calm ambient environment of density  a slightly greater than the initial jet density  0 .Without missing generality, the buoyant jet trajectory (jet centreline) is developed on the vertical plane yz, as shown in Figure 1.A local curvilinear orthogonal coordinate system (O 0 , , , ) is considered with unit vectors a, b, c, related to its coordinates , , .The local system has the  axis tangential at the jet trajectory,  axis horizontal and parallel to x axis of the basic Cartesian coordinate system (thus a ≡ i), and  axis perpendicular to the plane  and lying on the vertical plane yz.Let us consider a point M of the buoyant jet field, lying on the jet cross-section which coincides with the plane .The location of point M can be described by vector r 0 .According to Batchelor [14] and Marshall [16], the total infinitesimal motion of M to all directions can be generally described in either the basic or the local coordinate system and is expressed by the following double equality: where h 1 , h 2 , h 3 are positive scale factors, called metrical coefficients, which are defined as and |dr 0 | = √() 2 + () 2 + () 2 .Due to system orthogonality, the infinitesimal partial motion of M with respect to  axis, i.e., keeping constant the position of M from  and  axes, yields  = ,  = 0, and  = 0. Thus, The infinitesimal partial motion of M with respect to  axis, i.e., keeping constant the position of M from  and  axes, yields  = 0,  =  sin , and  = − cos .Thus, , h 2 becomes The infinitesimal partial motion of M with respect to  axis, i.e., keeping constant the position of M from  and  axes, yields motion at M  , where MM  corresponds to the infinitesimal partial change of position r 0 given as the vector variation r 0 , where r 0 = i + j + k = cℎ 3 .The vector r 0 +r 0 defines the new location M  of point M. Paying attention to the sketch in Figure 1, vectors r 0 and r 0 +r 0 can be described by the following relationships: where  0 = 0,  0 ,  0 are the coordinates of the jet centreline point O 0 .Abstracting ( 5) from ( 6) and then conducting pertinent trigonometric manipulations, the vector r 0 is calculated by the following formula: Since  → 0 ⇒ cos  ≈ 1 and sin  ≈ , (7) is further simplified as After dividing (8) by , the vector divergence r 0 / is taken as and, from (2) using ( 9), the following formula can be taken for h 3 :

Continuity Equation.
According to Batchelor [14] and Marshall [16], for the basic Cartesian coordinate system, the vector form of the continuity equation is written as where k = a+bV+c is the velocity vector with components u, v, w related to the local coordinates , , .Then, following [14] and taking into account (3), (4), and (10), ( 11) can be written with respect to the local coordinates as where ℎ ≡ ℎ 3 = 1 + (/).

Momentum Equation.
The vector form of the momentum equation is written as where  is the local density of the jet fluid, t is the time variable, F is the fluid acceleration due to gravity,  is the molecular viscosity of jet fluid, and p is the total pressure, which may be considered as the sum of the hydrostatic and dynamic pressure; i.e.,  =  ℎ +   .It is well known that the flow of a buoyant jet, with initial Reynolds number Re>1000, after a short distance from its exit ( ≈ 10) becomes fully turbulent.In order to take into account the turbulent contribution to the mean flow and buoyancy fluxes, the Reynolds substitution is considered according to the form =  +   , where f represents any dependent variable of the buoyant jet with a mean value  and turbulent fluctuation  .In such a condition, after averaging the terms of the governing equations, the contribution of the last term of the right part of ( 13) becomes negligible compared to the mean turbulent contribution and may be omitted [5,9].A common approximation for buoyant jet flows is the Boussinesq's approximation made for small initial density differences, i.e., when   / 0 ≈ 1 and Δ 0 =   −  0 <<  0 ≤   .Thus, the local fluid density  is considered nearly constant and equal to a reference density  r , where it is taken as either  r = 0 or  r = a .The usual assumptions and approximations are described in more detail elsewhere [17][18][19], while [5,9,20] have included the second-order terms of the turbulence contribution to momentum and buoyancy fluxes.Following Batchelor [14] and considering (3), (4), and (10), the term k ⋅ ∇k of ( 13) can be written with respect to the local coordinates as where 3 / = a/ + b/ + cℎ −1 / by substitution of V with u, v, and w, the terms k ⋅ ∇, k ⋅ ∇V, k ⋅ ∇ are given by the following relationships: The hydrostatic pressure at the point M is where p 0 is the hydrostatic pressure at O 0 (0,y 0 ,z 0 ).Using ( 16), the partial changes of the hydrostatic pressure with respect to , , and  are given by the following equations: since  0 / = 0,  0 / = 0, and  0 / = −   sin .The components of the vector of gravity acceleration related to , ,  are   = 0,   =  cos , and   = − sin , correspondingly.Consequently, the quantity (F − ∇) using ( 17) becomes The components of momentum vector on , , and  directions can be found by (13) and the vectors calculated by ( 14), (15), and (18).As justified above, the last vector of the right part of ( 13) is omitted in the present analysis, while  of inertial terms is substituted by the reference density  r .

Tracer Equation.
The tracer equation can be extracted from the mass conservation equation given by Batchelor [14] and replacing the fluid density  by the relative concentration definition given above, as Using the continuity equation, performing the Reynolds' substitutions, and averaging, the following final form for the tracer equation is taken: where overbars declare average values and S declares the source term.It is noted that the transient term is included in the momentum and tracer equations for the sake of generality.Nevertheless, in the analysis followed this term is omitted, since it concerns a quasi-steady-state phenomenon and it is not clear whether the assumptions of Gaussian similarity are valid in the initial stages of the transient evolution.The case of a transient buoyant jet phenomenon occurring in the staring time steps and/or when the flow field is confined, is beyond the scope of the present study and could be investigated in a future work.

Inclined Two-Dimensional Buoyant
Jet.An inclined twodimensional turbulent buoyant jet is considered, as shown in Figure 2, which is described by a local curvilinear orthogonal coordinate system (O 0 , , , ), where  axis is horizontal and parallel to  axis;  and  axes lay on the vertical plane  with  axis to follow the buoyant jet trajectory.Under these conditions, the PDEs describing the mean flow and mixing properties are ( 12), ( 22), ( 23), (24), and (26).These equations can be integrated on the transverse cross-sectional area of the jet, provided that the transverse profiles of the mean axial velocity and mean concentration are known functions.Experimental evidence has shown that the cross-sectional area occupied by a turbulent buoyant jet is circular and the aforementioned profiles may be well approximated by Gaussian type functions of the form: where  ≡  or c and   =    is the nominal width of the buoyant jet, which is equal to the distance  where  =  −1   , and   is the spreading coefficient of the property .
Index  indicates that the related value is at the  axis of the jet centreline.Equations ( 12), ( 22), ( 23), (24), and (26) must satisfy the following boundary conditions: where   =   V  is the shear stress due to turbulence,   is the physical half width of the buoyant jet, where   = /2 +     ,   = 2 1/2 , and V  is the entrainment velocity.The velocity  = 0 everywhere in the buoyant jet field due to the two-dimensionality of the flow.Therefore, (22) remains redundant and all partial changes of the variables with respect to  are equal to zero.The tracer is considered conservative and, therefore, the source term S equals zero.12), ( 23), (24), and (26) are integrated on the cross-sectional area of the twodimensional buoyant jet under the boundary conditions and approximations defined above.The integration of the terms of PDEs is managed using the Leibnitz rule for the terms regarding gradient to the main flow direction and the Green's theorem for the terms regarding the transverse gradients [15].The contribution of the fluxes of the turbulent terms of momentum and tracer with respect to the main flow direction to the mean flow fluxes is empirically introduced by coefficients based on available experimental findings [5,9].The cross-sectional area A occupied by the two-dimensional buoyant jet, where integration takes place, is spanned from  = -B w up to  = B w .In the initial region after the jet exit, which is called core region or zone of flow establishment (ZFE) because the turbulent character of the flow is not established in the whole jet width, the profiles are not of Gaussian type.These profiles are synthesized according to the methodology described elsewhere [20,21].To enable integration, the fluxes regarding volume , momentum m, weight deficit , and buoyancy  are defined and computed based on the profiles given by ( 27), as shown in Table 1.

Integration of Equations. Equations (
In Table 1, the coefficient   ≈ 1.1 introduces the turbulence contribution to the mean flux of momentum; b c is the spreading rate variable for the concentration field; the coefficient   introduces the turbulence contribution to the mean flux of buoyancy;  =   /  =   /  ;  = √(1 + )/(2);   ,   ,   ,   are shape factors introducing the contribution of the core region of the buoyant jet, where the profiles are not of Gaussian type; and  0 = 1.The coefficients   and  are computed by the empirical formula [9]: ; and  0 =  0 /√  0  is the initial Froude number of the buoyant jet [5].
Taking into consideration all data given above, ( 12), ( 23), (24), and (26) are integrated with respect to  on the actual cross-section of the two-dimensional buoyant jet, to give the following ODEs: Momentum Buoyancy conservation  =  0 (33)

Solution of ODEs and Comparisons.
Initially, the system of (31) throughout (33), along with aforementioned boundary conditions and flux definitions, is solved applying the 4 th order Runge-Kutta algorithm [15] and using either FORTRAN or EXCEL.The input parameters used in the computational code for the solution of the inclined twodimensional buoyant jet are the following: the effective acceleration of gravity  0 ' = 0.0301×9.81ms −1 , the width of the buoyant jet exit  = 0.002 m, the initial inclination angle  0 = 90 ∘ , 45 ∘ , 0 ∘ , or -45 ∘ , the initial exit velocity  0 = 0.365 or 0.486 ms −1 or alternatively the initial value of the Froude number  0 = 15 or 20, the initial value of relative concentration  0 = 1 at the exit, the exit coordinates with respect to the basic Cartesian coordinate system  0 =  0 =  0 = 0, the constant coefficients   = 0.132,   = 1.5,   = 1.21,   = 1.1,   = 1.04,   = 1.18, the limiting Richardson value   = 0.2907, the parameter  = 0.144, and the total spreading coefficient for the velocity field   /  = 1.9 and for the concentration field   /  = 2.5.The marching step assigned was  = 0.1D, for a slot width  = 0.002 m, and the computations stopped at the axial distance  = 1500.The computational time needed using FORTRAN was 0.19 s, irrespective of inclination angle.There was not any meaningful difference using EXCEL.Using the 4 th order Runge-Kutta algorithm, the convergence is obtained at each marching step and there is no need for any iteration.To analyze the sensitivity of the computational procedure to the results obtained, a tentime finer marching step was tested without taking practically different results.Consequently, the buoyant jet trajectory, the centreline mean axial velocities, and mean concentrations are computed.Then, the entrainment velocity or, alternatively, the entrainment coefficient defined as  = V  /  can be computed from (30).For a concise presentation of the results, all independent and dependent variables have been normalised by pertinent quantities coming from dimensional analysis.The distances , , , , , ,   ,   ,   are normalised by dividing them by ( 4/3 0 ); the velocities , w m , v B by ( 0  −2/3 0 ); the concentrations ,   by ( 0  −2/3 0 ); and the minimum dilutions (c 0 /c m ) by  2/3 0 .The same normalised nondimensional results are obtained for different values of g 0 ', D, and w 0 or F 0 .The buoyant jet centreline for initial inclination angles  0 = -45 ∘ , 0 ∘ , and 45 ∘ is shown in Figure 3(a), while in Figures 3(b) and 3(c) the trajectories of horizontal two-dimensional buoyant jets are compared successfully with measurements reported by Voustrou et al. [10].The data available by Cederwall [22] and Cuthbertson et al. [23], regarding this trajectory of a horizontal buoyant jet, have not be plotted in Figure 3(b), because these researchers have used a slot source in their experimental apparatuses without a protruding nozzle; thus, the Coanda wall effect occurring at the jet exit altered significantly the trajectory location.Nevertheless, the concentration or dilution remained unaffected.Paying attention to Figure 3(b), regarding the trajectory of the horizontal two-dimensional buoyant jet, it is obvious that the numerical results are in excellent agreement with the experimental data [10].However, regarding the trajectory of the negatively inclined two-dimensional buoyant jet shown in Figure 3(c), it is observed that the numerical results underestimate the trajectory.This fact may be attributed to the intense deformation of the buoyant jet field in the region of high curvature variation, where the flow and mixing symmetry with respect to the centreline plane is not maintained.This justification can be verified by the observation of approximately constant differences between model and experimental data after the high curvature region.
The variations of the normalised centreline mean axial velocities and the mean concentrations for initial inclination angles  0 = -45 ∘ , 0 ∘ , 45 ∘ , and 90 ∘ are shown in Figures 4(a) and 4(b), correspondingly.It is observed that near the jet exit and at large axial distances Ξ the plots in each figure are approximately in coincidence, while they differ in the region 0.5 ≤ Ξ ≤ 4.This behaviour is expected as all buoyant jets are initiated with the same velocity and relative concentration near their exit, while tending to obtain the same behaviour at large Ξ due to buoyancy, as they are gradually approaching the vertical direction.The variations of the normalised minimum dilutions for initial inclination angles  0 = -45 ∘ , 0 ∘ , 45 ∘ , and 90 ∘ are shown in Figure 5(a), while the normalised minimum dilution for the horizontal buoyant jet is shown also in Figure 5(b), where it is compared to available experimental data.For the case of  0 = -45 ∘ , Voustrou et al. [10] have reported the values of minimum dilutions measured at the maximum height of rise and at the return point, which are 0.53 and 1.16, respectively, while the model predicts the values 0.59 and 0.83.Despite the appreciable scatter of data, it is evident that the model predicts satisfactorily the minimum dilutions of a horizontal two-dimensional buoyant jet.The entrainment all cases examined and it is shown in Figure 6.
It is observed that the positively two-dimensional buoyant jets with  0 ≥ 45 ∘ approximately the same entrainment.The negatively buoyant jet shows the largest variations of a in the region 0 ≤ Ξ ≤ 2, where initially from Ξ = 0 and up to Ξ = 0.7 the flow is decelerating, while for Ξ > 0.7 the flow is accelerating, as it is implied from Figures 3(a) and 4(a).The minimum value of a = 0.0399 occurs at Ξ ≈ 0.7 and the maximum a = 0.1375 occurs at Ξ ≈ 2.

Conclusions
The main target of the present work is the mathematically correct application of the integral method, according to which integration of the system of PDEs, i.e., continuity, momentum, and tracer and/or buoyancy conservation, which governs the flow and mixing, should take place on the transverse cross-sectional area of the inclined buoyant jet, where the profiles of mean axial velocity and mean concentration are Gaussian.This necessity can be satisfied by developing a convenient system orthogonal curvilinear will monitor the flow motion with reference to the basic Cartesian system.Although the work of development of the curvilinear coordinate system was complicated, the integration of the governing equations resulted in quite simple ODEs, which are appropriate for numerical solution using the 4 th order Runge-Kutta algorithm programmed in FORTRAN or EXCEL.The numerical solution takes place very quickly, saving memory and computational effort.The mean flow and mixing properties taken from the solution concern the normalised trajectories of the two-dimensional buoyant jet for initial inclination angles  0 = -45 ∘ , 0 ∘ , and 45 ∘ , as well as the variation of the normalised mean axial velocities, mean concentrations, minimum dilutions, and entrainment coefficient, with respect to the normalised axial or vertical distance from the jet exit.These properties are presented in diagrams in analogy with the corresponding already confirmed results of the vertical two-dimensional buoyant jet.Additional comparisons were made for the cases with available experimental data in the literature.Finally, it is concluded that the whole model predictions are in satisfactory agreement with the available experimental data, thus justifying the complicated mathematical methodology described herein.

Figure 2 :
Figure 2: Configuration of a positively inclined two-dimensional buoyant jet in a calm environment.

Figure 5 : 3 Figure
Figure 5: Normalised centreline minimum dilutions (a) for several initial inclination angles and (b) for a horizontal two-dimensional buoyant jet in a calm environment compared to available experimental data.

Table 1 :
Mean flow kinematic fluxes and their initial values. d = √