ABSTRACT

Diffusion of a solute in turbulent flows through a circular pipe or tunnel is an important aspect of environmental safety. In this study, the diffusion coefficient of turbulent flow in circular pipe has been simulated by the Discrete Tracer Point Method (DTPM). The DTPM is a Lagrangian numerical method by a number of imaginary point displacement which satisfy turbulent mixing by velocity fluctuations, Reynolds stress, average velocity profile and a turbulent stochastic model. Numerical simulation results of points’ distribution by DTPM have been compared with the analytical solution for turbulent plug-flow. For the case of turbulent circular pipe flow, the appropriate DTPM calculation time step has been investigated using a constant β, which represents the ratio between average mixing lengths over diameter of circular pipe. The evaluated values of diffusion coefficient by DTPM have been found to be in good agreement with Taylor’s analytical equation for turbulent circular pipe flow by giving β = 0.04 to 0.045. Further, history matching of experimental tracer gas measurement through turbulent smooth-straight pipe flow has been presented and the results showed good agreement.

1. Introduction

The diffusion of gas and other particulate matter in pipe or channel flows is important aspect to meet the safety requirements. It controls the longitudinal spreading and the residence time of gases or other particulate matter throughout the pipe. Diffusion occurred in turbulent flow in circular airway has been investigated for a century. Several researches were done by conducting experimental works or numerical approaches. When a pulsed substance or solute is injected into a straight pipe flow, it is advected and diffused to a relative reference moving with certain average velocity. Diffusion in the turbulent pipe flow is mainly characterized by axial velocity profile and velocity fluctuation in flow direction, because the radial gradient of solute concentration is much less than that of flow direction and also radial diffusion is limited by its pipe wall. Furthermore, turbulent mixing motions at different radial positions enhance the diffusion degree in flow direction.

Taylor [1,2] and Aris [3] have made important contributions to develop theories about longitudinal diffusion in pipe flow. Taylor [2] also analyzed longitudinal diffusion in turbulent flow. He used the Reynolds analogy assuming radial diffusivity is analogous with heat and mass transfer in turbulent flow as well as transfer of fluid momentum. He neglected the contribution of molecular diffusion in both radial and axial directions which are negligibly small compared with turbulent eddy mixing diffusion in high Reynolds number. He also evaluated velocity at a certain distance from center of pipe and radial diffusivity as a function of universal velocity profile modified from Goldstein [4]. These assumptions are valid only for higher Reynolds number (Re > 2 × 10⁴) where viscous sublayer and transition layer are negligibly thin. Taylor [2] also proposed the relationship between the longitudinal diffusion coefficient, E (m²/s), against pipe diameter and turbulent shear velocity given by following equations:

(1)

(2)

where, τ (Pa) and u* (m/s) are shear stress and friction velocity in an arbitrary sub layer of the flow, ρ (kg/m³) is fluid density, d (m) is pipe diameter, f (-) is DarcyWeisbach friction factor and U_m (m/s) is cross sectional average velocity.

There are numbers of studies especially for atmospheric pollutant dispersion using random walk as basic method. Most of researches addressed the ideal homogenous turbulence. The pioneer was Taylor [5] who proposed continuous random walk theory for ideal homogeneous turbulent. For his random walk model, Pope [6] applied the well known Langevin equation, a stochastic differential equation. He pointed out the consistency condition in the case of homogeneous turbulent satisfies linear Gaussian model while for inhomogeneous turbulent satisfies non linear Gaussian model. Milojevic [7] simulate particle dispersion in an ideal homogenous flow by incorporating both Lagrangian and Eulerian method. Further, he found high particle concentration at low fluctuation velocity and vice versa. Kroger and Drossinos [8] applied random walk method to simulate thermophoretic particles deposition in turbulent boundary layer using Lagrangian method. He considered velocity, temperature fields and thermophoretic force as Gaussian random fields of which the mean values were obtained from law-to-law wall reactions and from Knudsen number dependent expression of thermophoretic force. The root-mean-square (r.m.s) fluctuation was calculated by polynomial fit with experimental data. Luhar, Ashok and Britter [9] developed random walk for dispersion in a convective boundary layer in inhomogeneous flow. They incorporated well mixed conditions, skewness in vertical velocity and Gaussian random forcing.

Pulsed injection measurement method by using NaCl into water stream in smooth glass pipe was conducted by Sittel, Threadgill and Schnelle [10], while Taylor [2] conducted both in smooth and artificially roughened glass pipes. Furthermore, they applied least square fitting method for their measurements, and showed higher prediction values compared to Taylor’s Equation (1). Higher values of diffusion coefficient were also observed by Hull and Kent [11] by using radioactive tracer injected into a long oil pipelines. Their reason was because of pipe bends and variation in elevation due to terrain profile.

Measurements of diffusion of gas-phase in smooth pipe flows were carried out by Keyes [12], while Davidson, Farqurharson, Picken and Taylor [13] conducted in a rough and long pipeline. Widodo, Sasaki, Gautama and Risono [14] also conducted tracer gas measurement in an underground mine ventilation airways. They found higher diffusion coefficient for the airways flow than those evaluated using Equation (1) for smooth pipe. This phenomenon is well explained by several researchers [15-17]. Compared with obtained data using a solute in water flow, data provided by Keyes [12] and Taylor [2] were measured at relatively lower Reynolds number region. They concluded that the effect of molecular diffusion cannot be neglected in low Reynolds number. Furthermore, because liquid has less molecular transfer, it has a higher Schmidt number compared to gas. Tichacek et al. [17] improved Taylor’s model by considering the molecular diffusion and mean velocity profile based on averaging velocity profiles. They reported that their calculation result deviates significantly for Reynolds number about 4.2 ´ 10⁴, due to the different sets of velocity profile. Thus, the velocity profile used in calculations has a sensitive effect on evaluated values of diffusion coefficient. Figure 1 shows evaluated diffusion coefficients reported in several researches. The empirical relationship proposed by Wen and Fan [18] is also plotted in Figure 1 and shows higher

prediction compared to [10] and [17]. Taylor’s analytical solution (Equation (1)) is also plotted together with 10% error bars. All have shown the scattered result to each other. Further, it is still necessary to develop further numerical method for gas diffusion in mine airways to simulate the longer travelling time and high diffusion coefficient than those for circular tube flow as demonstrated by Sasaki, Widiatmojo, Arpa and Sugai [19].

The advantages of proposed DTPM are that the calculations of concentration gradient in space or time domain, which are commonly employed in numerical simulations, are not required. It may reduce the computational time. It is also free of grid requirement and the visualization of points’ distribution is simple.

2. Numerical Models Formulations

2.1. Point Movement

The scheme of present numerical simulations has been developed by moving points with regards to velocity profile and turbulent intensities in a turbulent circular pipe flow, depends on radial position of the point in circular cross section. Figure 2 shows schematic variable definitions on the numerical calculation scheme using Cartesian coordinate (x, y, z) to describe tracer positions in circular airway with radius R or diameter d (=2R). The x is distance from the initial position in flow (longitudinal) direction, r is radial position from tube center perpendicular to wall surface, and φ is tangential direction perpendicular to x and r. As for fully developed turbulent flow, mean

velocities in r and φ directions can be zero.

Turbulent flow is characterized by its stochastic properties. The velocity at a given specific position fluctuates around its mean value. The velocity fluctuation intensity is known as root mean square (r.m.s) values which vary as function of r. The time of averaged flow velocities and the turbulent intensities at a certain position in cylindrical coordinate system (x, r, φ) are defined as (U, 0, 0) and (, ,) respectively. The moving of points may be treated easily on Cartesian coordinate system (x, y, z), therefore, in present calculations, turbulent intensities in (x, r, φ) directions are transformed into (x, y, z) directions of Cartesian coordinate. Assumed instantaneous turbulent fluctuations are (, ,) in (x, y, z) directions, these can be obtained by transforming velocity fluctuations (, ,) as follows (see Figure 3);

(3)

Suppose the position of m^th point dosed into the flow at origin is denoted with superscript showing the elapsed time, t = 0, its moving distances (Δx, Δy, Δz) during time step Δt are given by;

(4)

Its displacement is expressed at t + Δt and t by;

(5)

Laufer [20] conducted measurements of turbulent intensities and Reynolds stress of air flowing in a straightsmooth pipe with diameter, d = 0.254 m. Measured r.m.s values of turbulent intensities expressed by (, ,) in each axis (r, φ, x) were presented as ratios against u* as function of r/R. In present numerical simulations, the polynomial approximations for Laufer’s measurements (see Figure 4) were used to calculate the turbulent intensities at a point’s position, r/R.

2.2. Average Velocity Profile

If dimensionless value of longitudinal average velocity is defined as:

(6)

Dimensionless distance from the wall is given by:

(7)

where υ (m²/s) is kinematic viscosity. According to Kenyon [21], the relationships of U⁺ and y⁺ have been presented by Nikuradse with equations for three regions, that are viscous sub layer (y⁺ ≤ 5), buffer layer (5 < y⁺ ≤ 30) and turbulent zone (y⁺ > 30).

(8)

(9)

(10)

Equations (8) to (10) have been confirmed to agree well with equations proposed by Reichardt [22]. In present simulations, these equations were used to calculate axial velocity, U(r) of each point.

2.3. Turbulent Stochastic Model

A stochastic approach was applied to determine points’ diffusion process in the turbulent flow. In present numerical model, it is supposed that the point movements were based on turbulent eddy motion, which satisfies Gaussian probability density function (hereinafter GPDF) with a standard deviation equal to the turbulent intensities or the r.m.s value of velocity fluctuations in each direction (see Rouse [23]). In the simulations, three pseudo-random numbers follow GPDF were generated using Box Muller algorithm to calculate turbulent velocity vector (, ,).

As described previously, Laufer’s measurement results of turbulent intensities were presented in normalized values over shear velocity, u^*, which is function of cross sectional average velocity, U_m, and friction factor, f. The relationship between f and Re was presented with empirical equation by Colebrook [24];

(11)

2.4. Turbulent Reynolds Stress Model

In turbulent shear flow, fluid particles are translated from slower region to faster one and Reynolds stress is generated. It shows a time-averaged correlation between longitudinal and radial velocity fluctuations in the shear flows. In present study, effects of Reynolds stress correlations have been investigated by giving relationship between x axis velocity fluctuation, , with velocity fluctuation in r direction, , expressed as;

(12)

In the simulations, the value of including its sign was firstly given as a random number following GPDF described previously in preceding chapter, then absolute values of || were generated using similar method, but the signs of u’ were decided to satisfy Equation (12).

2.5. Boundary Condition and Initial Condition

In fact, point’s displacement near wall and its wall interaction has not been well understood to simulate breaking sub-layer. The boundary condition for Lagrangian random walk still needs a calculation model. Several boundary models have been proposed in previous studies. Those methods depend on the physical and numerical factors considered in the simulation [25].

In present DTPM simulations, the reflection boundary condition at the airway wall is satisfied by a repositioning numerical treatment if points are out of the flow region. Thus, by this boundary condition the zero-flux condition at wall can be satisfied. A reflection boundary scheme is modeled as shown in Figure 5 and is used in present simulation. Similar boundary model was also proposed by Szymczak and Ladd [26,27].

Another boundary condition namely rearrangement model has been investigated. In this model, random number is continuously generated until the point is repositioned within the flow regime. Szymczak and Ladd [26] also reviewed this kind of boundary condition as multiple rejection method. Figure 6 shows the comparison of two models. It can be seen that the reflection model resulted in higher dispersion than the rearrangement model. The reason is because the reflection model allows higher possibilities for points to be repositioned at near wall region. These movements make higher “retaining effect” in low velocity layer near wall. Further simulations indicated that the reflection model shows closer results to Taylor’s analytical solution. Thus, reflection model is adapted for present DTPM numerical simulations.

Other zero-flux boundary conditions also have been proposed. Drazer and Koplik [28]; Kurowski, Ippolito,

Hulin, Koplik and Hinch [29] proposed rejection boundary where the points do not change its position in the given time step and other random diffusivity is recalculated until the movement satisfy r < R. Salles et al. [30] and Maier et al. [31] suggested interruption boundary condition at which points stop at wall and its time is incremented by modified time step, λDt where λ = 0 ~ 1 is adjustment factor. However, only the reflection method imposes the zero-flux boundary condition while other boundary conditions lead to incorrect concentration profiles near wall boundary [26].

3. Calculation Results and Evaluation of Diffusion Coefficient

3.1. One-Dimensional Diffusion

Taylor [1,2] proposed that concentration distribution of solute after certain time, t, is symmetrically distributed follows the partial differential equation. He proposed concentration gradient at x and r direction which move to x direction with constant average cross sectional velocity, U_m;

(13)

Several studies (see Wen and Fan [18]) also proposed similar numerical expression of solute diffusion plug flow model.

The solution of Equation (13) can be obtained by assuming that molecular diffusion is neglected and concentration gradient in radial direction is negligible. The variable, E shows effective diffusion coefficient in axial direction. Since the center of dispersed solute is assumed to be at x = U_mt after elapsed time, t, the solution of Equation (13) can be given with an equation similar with Gaussian distribution:

(14)

where G (m³) is total volume of released solute at x = 0, and A (m²) is cross sectional area of pipe.

In the DTPM simulation, concentration of diffused points, C, can be calculated by;

(15)

where Δm is number of point counted in the numerical control volume which located at a certain downstream position, given by ΔV = ΔXπd²/4 and ΔX = U_mΔt. Suppose the total number of points released from the origin is M_D, The normalized concentration of points is expressed as C_TA/M_D.

3.2. Simulation of Turbulent Flow through Circular Airway

The DTPM simulations on turbulent flow of a straightsmooth airway have been carried out. Figure 7 shows the evolution of point’s distribution after several elapsed time since released into the flow (U_m = 1.5 m/s, d = 0.5 m). It can be observed that the asymmetry of points distribution is gradually diminished as the travelled distance increases.

Figures 8 and 9 show the results of simulation for flow with d = 2 m, U_m = 4 m/s (f = 0.01315) and d = 4 m, U_m = 5 m/s (f = 0.0112) respectively after t = 100 s. It can be inferred that different calculation time step, Δt, gives different evaluated value of effective diffusion coefficient. To consider this effect, the evaluated value of effective diffusion coefficient is plotted against dimensionless value, β (-) defined as:

(16)

where is r.m.s value of streamwise velocity fluctuation in the center of pipe. It may be possible to

regard β as ratio of the average mixing length to pipe diameter.

Figures 10(a)-(d) shows evaluated values of effective diffusion coefficient of different time step for Re = 5 × 10⁵, 7.5 × 10⁵, 9.4 × 10⁵ and 1.25 × 10⁶. The values of E calculated by Equation (1) are also presented as comparisons. From the results, it can be seen that the evaluated values of E from present DTPM for β = 0.04 ~ 0.045 intersect with those Taylor analytical solution given by Equation (1).

The results of evaluated E of DTPM showing relatively non-linear and inversely proportional correlation to the value of β, before gradually attain constant value as higher value of Δt is applied. Further simulations were done by setting β = 0.040, 0.043 and 0.045 at different flow conditions. As shown in Figure 11, the results indicate linear relationship between Taylor’s analytical equation and ones evaluated by DTPM.  The value of β = 0.043 used for DTPM simulations is appropriate in order to get longitudinal diffusion coefficient, E, which has good agreement with Taylor’s analytical solution.

4. Simulation of Tracer Gas Experiment through Smooth-Circular Pip

Tracer gas diffusion experiment was conducted at a laboratory scale by Widodo [32]. The experimental apparatus mainly consists of smooth circular pipe as scaled pipe, gas injection apparatus, and gas measurement apparatus.

For the straight airway, he used a pipe with smooth lining 0.025 m diameter, 30 m length and placed horizontally. Tracer gas was released by breaking balloons filled with methane (CH₄) and measured arriving concentration at the end of pipe. Methane concentration was measured by an original infrared adsorption gas detector

using 3.4 mm He-Ne laser with an infrared (IR) sensor, air pump, air mass flow meter and amplifier as shown in Figure 12. Sampling rate was set as the maximum flow rate at which the gas sensor could still detect the gas concentration correctly. The data were recorded with a data logger connected to a PC.

Beforehand, the conversion of voltage reading by data logger to gas concentration was calibrated based on calibration data with the standard methane-air mixture. The validity of concentration reading from IR sensor was crosschecked using gas chromatograph and showed good linear fit.

The measurement was conducted for Reynolds number, Re = 6085 (U_m = 3.87 m/s, d = 2R = 0.025 m). Measurements for higher Re were also done; however, the data were inadequately acquired due to insufficient sampling rate to compensate higher flow velocity.

For the DTPM simulation, calculation parameters were decided based on the experimental properties such as pipe’s length, average cross sectional velocity and pipe’s diameter. Value of friction factor, f, was calculated using Equation (11) and calculation time step was defined using Equation (16) as Dt = 0.0092 s. To verify the effect of initial points’ distribution, three different conditions were considered; 1) point source (x₀ = 0, r₀ = 0); 2) uniformly distributed at x₀ = 0 (0 ≤ r₀ < R); 3) uniformly distributed at x₀ ± 2d.

Figure 13 shows the evolution of radial and longitudinal point’s distribution in five consecutive times for point source initial condition. The results shown in Figure 14 implied that the difference of initial points’ distribution has no significant effect on result of DTPM. The length of flow domain is long enough for the radial diffusivity to attain radial homogeneity of flow. It may also imply that the sudden break of balloons during release has less effect on the arriving concentration as the measurement distance is long compared to diameter d pipe length. In general, the results confirm that present DTPM is able to simulate the turbulent in straight and smooth circular pipe.

In the real case, the utilization of tracer gas measurement is not merely straight pipe, but also tunnels network [19,32-35]. This kind of network is assembled in the form of interconnecting tunnels which allow airflow to be separated or rejoined at junction. The mechanism of points tracking as proposed in this study can be developed to consider network flow by combining with a scheme to treat flow separation. The developed scheme is supposed to be able to simulate point’s distribution at arbitrary position in the network and allow easy dispersion evaluation of gas or other particulates spreading.

5. Conclusions

In this study, effective diffusion coefficients of turbulent flow in circular pipe or channel have been evaluated by the Discrete Tracer Point Method (DTPM), a Lagrangian numerical simulation method. The present study is summarized as follows:

1) DTPM simulation procedures have been presented to

simulate the displacement of points released into straightsmooth circular airway flow by giving turbulent average velocity, intensity of velocity fluctuation and Reynolds stress;

2) A simple procedure to represent diffusion of points in the airway flow has been employed by generating random number which satisfies Gaussian probability functions with value of turbulent intensities in each direction as standard deviations;

3) Appropriate calculation time step was proposed by matching with Taylor’s analytical equation by considering the ratio between the average mixing lengths to airway’s diameter, β @ 0.043;

4) The result of matching curve of DTPM with experimental result by considering different initial points distribution has shown that the present DTPM can be used to simulate turbulent advection-diffusion in smoothstraight pipe representing the ideal model of straight pipe or channel;

5) Present DTPM simulation has promising possibility to simulate the network flow by combining with a scheme to treat flow separation at junction.

6. Acknowledgements

Authors would like to thank Japan Ministry of Education, Culture, Sports, Science and Technology and Kyushu University Global COE program for the financial support.

REFERENCES

Nomenclature

A = Cross-sectional area of pipe (m²)

C = Gas concentration (−)

C_T = Number of points located within numerical control volume (−)

d = Diameter of pipe (m)

E = Longitudinal diffusion coefficient (m²/s)

f = Darcy-Weisbach friction factor (−)

M_D = Total number of released points (−)

r = Radial position of a point (m)

R = Radius of pipe (m)

Re = Reynolds number (−)

t = Elapsed time (s)

U = Time-averaged flow velocity in x direction u^* = Friction velocity (m/s)

U⁺ = Dimensionless value of flow velocity (−)

y⁺ = Dimensionless distance from the wall (−)

Um = Cross-sectional average velocity (m/s)

() = Turbulent intensities in cartesian coordinate (x, y, z) (m/s)

= Stream-wise velocity fluctuation in the pipe center (m/s)

(, ,) = Turbulent intensities in cylindrical coordinate (r, φ, x) (m/s)

(, ,) = r.m.s value of turbulent intensities in cylindrical coordinate (r, φ, x) (m/s)

(Δx, Δy, Δz) = Displacement vector of a point (m)

β = ratio between average mixing in stream-wise direction to the diameter of pipe

Γ = Released volume of tracer gas (m³)

Δt = Numerical time step (s)

ΔV = Numerical control volume (m³)

λ = Time adjustment factor for interruption boundary condition (−)

ρ = Density of fluid (kg/m³)

τ = Shear stress (Pa)

υ = Kinematic viscosity of fluid (m²/s)

NOTES

^*Corresponding author.