Miao Liu,Zhengming Gao,Yongjiu Yu,Zhipeng Li*,Jing Han*,Ziqi Cai,Xiongbing Huang

State Key Laboratory of Chemical Resource Engineering,School of Chemical Engineering,Beijing University of Chemical Technology,Beijing 100029,China

Keywords:PIV LES Turbulent kinetic energy Pseudo-isotropic assumption

ABSTRACT Con fined impinging jet reactor(CIJR)is a typical process intensification device used in the chemical industry.In this study,two dimensional Particle Image Velocimetry(PIV)and Large Eddy Simulation(LES)method were used to investigate the flow field in a CIJR with jets of diameter 3 mm under highly turbulent condition.The results showed LES can predictthe velocity and Turbulence Kinetic Energy(TKE)distributions in the reactorwellby comparing with the PIV results.In the CIJR,the stagnation point fluctuates with the turbulence,and its instantaneous position accords with the normaldistribution.Three methods,including s–t representation,Lumley–Newman triangle and A–G representation,were used to compare the turbulence anisotropy in the mixing chamber.It was found that the anisotropy in the impinging area and at the edge of impinging jet was strong and the maximum deviation was up to 40%.The results from 2D PIV would lead to an overestimation of the turbulent kinetic energy as much as 20%to 30%than the results from the three dimensional numerical simulation.

1.Introduction

Con fined impinging jet reactor(CIJR)is a typical process intensification device,which can be applied in absorption and desorption[1],extraction[2],liquid mixing[3–5],bioreaction[6],precipitation[7],and crystallization[8–10].The principle of CIJR is to make at least two jets colliding with each other and then the momentum,heat,and mass transfer are intensified in the impinging region.In recent years,there has been a renewed interest in this device because it can achieve excellent mixing much faster than the traditional reactors,and it has great ability to deliver mixing time less than the characteristic process time for fast reaction processes.The key performance is producing a region with highly turbulent energy dissipation and making sure that the liquid can flow through thatregion[11].Forthese reasons considerable attention has been devoted to the study of flow structure in CIJR,such as mean velocity,Root Mean Square(RMS)velocity,turbulent kinetic energy(TKE)and stagnation point offset.

Particle Image Velocimetry(PIV),as a commonly used experimental method in flow structure test,has been widely used in the research of CIJR.Different from Laser Doppler Velocimetry(LDV),PIV can conveniently be used to obtain the instantaneous velocity on an illuminated plane.Landreth and Adrian[12]measured the velocity field of a turbulent circular jet impinging on a plate for the first time by using a two dimensionalPIVsystem,thus revealing flow structures and various stages of vortex generation.Li et al.[13]employed two dimensional PIV to study the dynamic behaviors of an impingement plane in a 50 mm CIJR under inlet flow excitation by the PIV technique with Reynolds number ranging from 75 to 150,then reveal the structure and number ofvortices as well as the variation ofthe mixing effectwith the entrance lf ow excitation.Tu et al.[14]investigated the flow field at different geometric parameters in CIJR with Re ranging from 21 to 4642,observing the half de flecting oscillation behavior in CIJR for the first time and critical parameters of the half de flecting oscillation were obtained.Somashekar et al.[15]also employed micro-PIV to observe and analyze the flow field at three inlet Reynolds numbers and presented the pointwise and spatial turbulence statistics including Reynolds stresses and spatial correlations of velocity fluctuations.Gao et al.[16]investigated flow characteristics in CIJR at high Reynolds numbers ranging from 10620 to 21210 by a PIV system and found that Re had almost no in fluence on the normalized average velocity and TKE distributions in the CIJRs in fully developed turbulence,while the relationship between the stagnation point and mean velocity ratio and the effect of confined space on flow field were also discussed.

Computational Fluid Dynamics(CFD)has grown tremendously in recent decades.The use of CFD analysis can result in a significant time and cost saving[17].In the description of turbulence,k–ε is a classical modelbased on RANS method,such as in Liu etal.'s work[18].However,it was found that simulation using the steady-state k–ε model can only be satisfactory as a first approximation of the flow,and the turbulence kinetic energy by CFDpredictions wasunderestimated.Large Eddy Simulation(LES),as another important method in the description of turbulence,was regarded as a better way to improve the prediction by various authors.Icardi et al.[19,20]employed Large Eddy Simulation(LES)and Direct Numerical Simulation(DNS)to simulate the flow field in a 1 mm CIJR and the predictions were validated against experimental data,where the simulations generally resulted in good agreement with experimental data on both mean velocity and RMS velocity.The LES result calculated with different spatial discretization schemes and the subgrid scale(SGS)models were found not crucial for a good prediction of the turbulent behavior of the system.These results also implied that this kind of intensively turbulent mixer must be provided with unsteady boundary conditions.Li et al.[21]studied turbulent opposed jets and compared the experimental results from a hot-wire anemometer with the CFD simulation,showing that the stagnation pointoffsetwas in good agreementwith the visualization results.

In the previous study,based on the 2D PIV method,it is inevitable that the isotropic hypothesis willbe used to calculate the turbulence intensity.In fact,the researchers usually do not exactly know how much the fluid pulsation is,which is perpendicular to the target plane.In addition,the researchers do not know the deviation of this method too.In this paper,the two planes of the mixing chamber of a CIJR with 3–3 mm jets is measured by 2D PIV and compared with LES results.After obtaining the average velocity and the RMS velocity of PIV and LES respectively,the turbulent flow anisotropy is calculated and evaluated by means of three different methods,which can provide a guide for the turbulence evaluation and design of the CIJR.

2.Experimental

The experiment setup in this work consists of two parts, flowdelivery and two dimensional PIV systems.The working fluid was deionized water which was pumped from a storage tank by two multistage centrifugal pumps(CRN3-36,Grundfos,Denmark).The flow was controlled by a frequency converter(VLT5005,Danfoss,Denmark).The deionized water(995 kg·m-3,0.994 mPa·s) flowed in a PPR pipeline and the volume flow rate was measured by a liquid turbine flow meter(LW-10Z1M1SNS,Gallop,China).More information on the system can be referred in Gao et al.[16]

The geometry of the CIJR used in the presentwork is shown in Fig.1.In order to ensure the highly accurate coaxiality of the two jets,wire cutting was used in the manufacture.l/d,in which l is the length of inlet pipe and d is the diameter of jets,should be at least 10 in order to ensure that the pipe turbulence flow was fully developed.In order to make the inner of the impinging chamber visible,the top,front and back sides of the chamber were all made of perspex for the image capture and laser illumination for PIV.

A two dimensional PIV system(TSI,USA)was used to measure the tracer particle displacements in CIJR.This PIV system was conducted using a dual pulsed laser(New Wave,Nd:YAG 532 nm,200 mJ,15 Hz)and a single frame-straddling CCD camera(PowerView Plus,11 MP)with the resolution of 4008×2672 pixels.Frame and laser sequencing were controlled by a synchronizer(TSI,Laserpulse 610035).The particles seeded in the flow were hollow spherical glass whose density was 1500 kg·m-3and the average diameter was 10 μm.The particle Stokes number St that characterizes the ratio of particle response time to the flow time scale can be calculated as

where γ is the characteristic strain rate of the flow and can be approximated by 2u/d,vfand ρfare the fluid kinematic viscosity and density.For this tracer particle,the Stokes number is as low as 3.7×10-7,which is much smaller than the critical value 0.1,showing that it has excellent flow tracer fidelity with the water.

The time difference between laser pulses was chosen as 4 μs for impinging velocity u=7.07 m·s-1(Re=21230)in the experiments,which was carefully optimized to ensure that the maximum in-plane and out-of-plane displacements of seeding particles were less than one quarter of the interrogation window sizes and the thickness of light sheet.In the resolving of the flow region,a Fast Fourier Transform(FFT)cross-correlation algorithm was applied to the 32×32 pixelinterrogation windows with 50%overlap.Considering there is a resolution of 9.27 μm·pixel-1in the image,the resolution of a vector was 0.148 mm.

In the determination of sampling size in PIV experiments,the mean velocity and RMS velocity at x=0,z=1 versus sampling size are shown in Fig.2.When the sample number is over 250,the statistical result reaches basically stable.Therefore,400 samples are captured,which can ensure that the calculation of velocity is independent of the sample numbers.

3.Large Eddy Simulation

Although DNS is an accurate way in predicting the flow field by directly solving the Navier–Stokes equations,it costs too much machinehour.LES respectively employs DNS and the subgrid model to calculate the large scale and small scale vortices,and the basic governing equation of LES for incompressible flows after being filtered is:

Fig.1.Dimensions and solid model of a CIJR.

Fig.2.Statistical convergence of velocities at x=0,z=1 on plane XZ.

where the grid stress τijis modeled as

The subgrid scale stresses represent the interaction between the resolved larger scales and the unresolved smaller scales of motion.In the description of eddy viscosity,the Smagorinsky model[22]is always used.The unknown subgrid scale stress tensor τijand the filtered rate of strain tensorare related by

In the Smagorinsky model,the eddy viscosity is modeled as

where Csis the Smagorinsky constant which is set as 0.1 in this simulation,Δ is the filter width andis modeled as

The Smagorinsky model can be used in the numerical procedure of the Navier–Stokes equation by adding only one eddy viscosity coef ficient,which has a drawback—an excessive dissipation.In orderto overcome the defect of the eddy dissipation,Germano et al.[23]and Lilly[24]put forward a dynamic eddy viscosity model that the flow field is filtered twice and the small scale fluctuation from the coarse filtration is similar to the pulsation from the fine filtration.

The numerical simulation of the CIJR was performed by using the commercial software FLUENT 14.0.The geometric parameters in the simulation were the same as those presented in the experimental part.The total computational grid consisted of over 900000 structured hexahedral cells after the grid independence test(600000 and 1200000 cells)which is showed in Fig.3.

LES was carried out using the Smagorinsky–Lilly SGS model under a constant time step of 2×10-5s,and the maximum and average Courantnumberis about2.5 and 0.9,respectively.Inletboundary conditionswere obtained with the same procedure aswith RANS.However,it is worth reminding that in Fluent LES is carried out by superimposing stochastic components of the flow at the inlet boundaries by adding random perturbations.So,in this work turbulentintensity I is calculated according to the fully developed pipe flow:I=0.16×Re-0.125.RMS velocity profile under different inlet turbulent intensities is showed in Fig.4.As can be seen from Fig.4(a),when the turbulent intensity I ranges from 4%to 8%,the rmsU and rmsW change little,and I=4.6%is used according to the Reynolds number21230.No-slip boundary conditions were used for all solid walls.The employed spatial numerical scheme is central difference(pressure–velocity coupling solved with PISO algorithm),and the temporal discretization is second-order implicit.

In the characterization of turbulent flow,Lumley[25]proposed a method in the accurate evaluation of the local anisotropy and isotropy for turbulence.This method is based on the analysis of the normalized Reynolds stress anisotropy tensor:

Fig.3.(a)Mean velocity and(b)RMS velocity at different grid numbers.

Fig.4.RMS velocity profile of different inlet turbulent intensities.

where k is the turbulentkinetic energy.When i=j,δij=1,i≠j,andδij=0,the tensor is the basis of turbulence anisotropy analysis.

To characterize local anisotropy of turbulence by means of a single parameter,Derksen et al.[26]suggested using the distance L of the point from the origin mark O(isotropic turbulence)in Lumley's triangle:where II and III are the trace of the anisotropic matrix which is used to describe the tensor of Reynolds stress.

Fig.5.Instantaneous flow field on XZ plane of CIJRs at two moments(the red circle marks the position of the vortex;the red line marks the plane of the stagnation point;the scale of the arrow is 10 m·s-1).

Fig.6.Instantaneous stagnation offset.(a)PIV,7.07 m·s-1,1 Hz.(b)LES,7.07 m·s-1,5×104 Hz.

Many researchers showed that the turbulence in the impeller discharge region and the trailing vortex was highly anisotropic in stirred vessel[27,28].

4.Results and Discussion

4.1.Instantaneous flow field and the stagnant point

The instantaneous flow fields in the CIJR chamber are shown in Fig.5.As can be seen,the instantaneous stagnation point is unstable and oscillating with the impinging plane.In the flow field,there are lots ofsmallscale eddies,which change with the time in location and intensity.The process of the generation, flow and dissipation of these eddies can strengthen mass and heat transfer.

The location of stagnation point,including the average and instantaneous position,is an important factor in the design of CIJR,which is closely related to the long-term operation ofCIJR.And the position ofinstantaneous stagnation point also re flects the vibration of flow filed.

The instantaneous stagnation points in experimentaland simulation results operated at u1/u2=1.0 were statistically computed,as shown in Fig.6.It can be seen thatthe instantaneous stagnation points oscillate in the range of X=±2.5 mm.

Then,the frequency distribution histogram of instantaneous stagnation position is shown in Fig.7.Both ofthe experimentaland simulation data can be expressed in the normal distribution,as the following equation.

Fig.7.Distribution histogram of stagnation offset.

The position of stagnation point in PIVresultis more dispersed,indicating the randomness of instantaneous stagnation point is strong.However,the correlation coefficient in the regression of LES result is high.The fluctuation ofexperimentdata is largerbecause ofthe instability ofpumping even in such a multi-stage centrifugalpump.In CIJR fluid high-speed impact process is strongly turbulent,with instability that can lead to pipeline characteristic curve changes over time,and make the centrifugal pump work under unstable conditions,making the pump's flow unstable.In addition,the centrifugal pump used in the experiment itself has a 4%error,which is one of the reasons for the instability of the pump outlet flow.

4.2.Velocity profiles

In the literature ofGao etal.[16],the mean velocity and TKE distributions under different Reynolds numbers were investigated and it was found that the normalized flow fields in fully turbulent flow are similar.Therefore,Re=21210 is chosen as the operation condition in the comparison of experimental and simulation results.

As can be observed in Fig.8,experimentalresults show a typical flow characteristic of CIJR on planes XZ and XY.It represented that two impinging jets collide at the center of the cubic mixing chamber.Mean velocity decreases sharply along the x-axis,and the average velocity at the stagnation pointis nearly zero.Aftertwo impinging jetscollide,the fluid nose-dives and comes out of the outlet port of the cubic vessel.

The comparison of mean velocity and RMS velocity in x,y,and z directions between PIV and LES is shown in Fig.9.In mean velocities U and W,LES result is in a good agreement with PIV data,but there is a slight difference in the peak of V velocity.The RMS velocities rmsV and rmsW in LES and PIV are similar,but the fluctuation of rmsU velocity in the experiments is larger than the LES data.In X＞2 mm and X＜-2 mm along x-axis,velocity U remains stable.But in the range of-2＜X＜2 mm,the velocity U decreases sharply and reaches zero at the center of the mixing chamber.In Fig.9,the RMS velocity is large in the area where the mean velocity is low.In the main area of impinging jet,the mean velocity distribution in two jets is consistent with turbulence characteristics of the fluid,which indicates two impinging jets are fully turbulent.

Fig.8.Mean velocity and TKE distributions(k,m2·s-2)on XY and XZ plane.(The scale of arrow is 10 m·s-1.)

Fig.9.Mean velocity and RMS velocity at different axes using PIV and LES.

In the range of Y＜-2 and Y＞2 mm,V increases as Y decreases,and reaches a maximumof4 m·s-1atthe point Y=±2 mm,which is about 57%of the maximum velocity of impinging jet.Near the center of the mixing chamber,the velocity decreases and approached to zero at the stagnation point.After the fluids collide in the mixing chamber,mean velocity decreases to zero soon.Then,the fluid is extruded and pushed by the high pressure in the chamber,leading to the sharp increase of U.Along the z axis,velocity W and rmsW distributions are similar to that of V and rmsV velocity.

4.3.TKE distribution and the calculation

The full expression of turbulent kinetic energy(k)is known as follows:

where rmsU,rmsV,and rmsW are fluctuating velocity in the x,y,and z directions,respectively.Since the PIV system used in the present work is two dimensional,the fluctuating velocity in the third direction cannot be obtained directly,thus,pseudo-isotropic assumption is often used and then,

In numerical simulations with three dimensional mesh,all the velocity components can be derived.But we deliberately used Eq.(12)to calculate the TKE in LES and compared it with the experiment results.Therefore,LES data is used to analyze the accuracy of TKE calculation between two dimensionalpseudo-isotropic assumption(showed in Fig.10,marked as LES)and 3D results(showed in Fig.10,marked as LES-3D).

The TKE along the x-and z-axes is shown in Fig.10.It can be noted that there is an obvious difference in the calculation of TKE between pseudo-isotropic assumption and the three-dimension results.The maximum value of TKE in the three-dimension calculation is 11 m2·s-2,while that in two-dimension is 16 m2·s-2,with the relative difference of about 45%.

In the impinging zone,strong turbulence anisotropy makes the pseudo-isotropy assumption not applicable,because the value of TKE is directly related to the RMS velocities in three directions.It's showed in Fig.11 that the value of rmsU is larger than rmsV and rmsW in-2 mm＜X＜2 mm.Especially,at X=0 mm and Z=0 mm the difference between rmsU and rmsV or rmsW reaches maximal,leading to the TKE calculation value of Eq.(12)being larger than Eq.(11).

To show the difference between Eqs.(11)and(12)more convenient,a RATIO representing the proportion of Eq.(12)to Eq.(11)is defined,

The contour of the RATIO on plane XZ is shown in Fig.12.

The RATIO deviates from 1.0,indicating that the two dimensional approximation of TKE is larger than the actual value.In the impinging zone,the maximum deviation ofTKE reaches more than 40%.Following that,in the discharging area and at the edge of the impinging body,the maximum deviation of TKE is in the range of 20%–30%.This indicates that the turbulence anisotropy in the mixing chamber is not negligible.In addition,when using the pseudo-isotropic assumption in the calculation of the TKE in a 2D-PIV measurement in the present coordinate system,rmsV2=0.5(rmsU2+rmsW2)is assumed in Eqs.(11)and(12).According to the result shown in Fig.11,the difference of TKE by 3D and 2D calculation is expressed by the RATIO(Eq.(13)).If the value of RATIO needs approaching 1,which means the calculation of TKE by 3D and 2D is the same,a parameter γ is defined as follows,

Substitute Eq.(14)into Eq.(13),

where γ is the correction factor in the calculation of TKE in 2D-PIV measurement.

Fig.10.TKE distribution of different calculation methods.

Fig.11.Distribution of RMS velocity.

Fig.12.The contour of the TKE RATIO.

For example,as can be seen from Fig.12 and Eq.(15),the values of RATIO for the impinging area and the edge of impinging area region are 1.4 and 1.2,soγin these two areas should be 0.071 and 0.25 respectively.It indicates that in different areas the parameter γ should be different.

4.4.Turbulence anisotropy analysis

The above contents indicate that the turbulent anisotropy in the mixing chamber is very strong,and three different methods are used to evaluate the characteristics in this section by utilizing the results from LES.

Lumley and Newman[29]proposed a method forevaluating the turbulence anisotropy,and used the normalized Reynolds stress invariants represented the inherent characteristics of the Reynolds stress.Besides Lumley and Newman's method,Escudié and Liné [30]developed another two evaluation methods:(1)s–t representation based on the eigenvalues of the anisotropy and(2)A–G representation based on the axisymmetric parameter A and two-dimensional parameter G.

Reynolds stress is the second order symmetric tensor,with the normalized Reynolds stress invariants partial in the representation of inherent characteristics of the Reynolds stress.In an anisotropic turbulence,the normalized Reynolds stress bijcan be used to describe the difference between three dimensional isotropic and anisotropic.

When i=j,δijshould be 1.And when i≠ j,δijshould be zero.After Gauss transformation,the matrix eigenvalues γij∣ψcan be derived and substituted into the matrix:

(1)s–t representation

s and t are the diagonal elements of the matrix.Since bii=0,the trace ofthe anisotropy tensormatrix is zero.The firstand second principal stress is set as s and t,and the third principal stress is–(s+t),then

when s≥t≥ -(s+t),

(2)Lumley–Newman Curved edge triangle method

The parameters I,II and III,which are trace of the matrix,square trace and cubic trace,are defined as:

For the second order symmetric tensor,they are invariant.The trace of bijis zero,so I=0

Hence,Reynoldsdeviatoric stresstensorcan be expressed by these two parameters(II,III)in evaluating the turbulence anisotropy.

(3)A–G representation

On plane II–III,the axial symmetry boundary condition is i:

Then the parameter A is defined as

On plane II-III,the two dimensional boundary condition is:

Then the parameter G is defined as

By summarizing above algorithm,the property ofanisotropy parameters is show in Table 1.

The contour of turbulent anisotropy in CIJR by those three methods is shown in Fig.13.It can be noted the region with high anisotropic turbulence is mainly in the impinging zone and the edge of impinging jet,shown in Table 2.

Table 1 Property of anisotropy parameters

Fig.13.Contour of turbulent anisotropy.

Table 2 Value of anisotropy parameters in CIJR

The k–ε model is based on the turbulence viscosity concept and the assumption of isotropy.As can be seen in Table 2,in these two regions,the three methods show that the turbulence approaches 1D isotropy and deviation from3Disotropic.Therefore,the deviation ofTKE calculation in these areas using 2D method is very large.That's why the k–ε model is not recommended when describing the TKE and turbulence dissipation rate of CIJR.In Liu and Fox's work[18],they obtained the similar results in CIJR at Re in the range of 211–1003 and predicted that the value of k is higher than the experimental results at low Re.

There is no shear stress in the isotropic flow field.Actually,the shear stress cannot be ignored in the impinging zone and the edge of impinging jet,where velocity gradient was very high.

5.Conclusions

In this work,the flow characteristics in CIJR were investigated by two dimensional PIV experiment and LES.Through extensive data collection and analysis,it is found that the stagnation point in CIJR is not fixed,but will fluctuate with the time and flow,which is also observed in the LES.The stagnation pointin this mixing chamber is approximately within-2.5 mm＜X＜2.5 mm.The position of the stationary point in the CIJR mixing chamber can be described by a normal distribution,and the distribution of the location of the stationary points obtained by the experiment and simulation is generally matched.

The results of comparison experiment and numerical simulation show that the velocity field is consistent.In addition,using pseudoisotropic assumption to calculate the three dimensional LES results,the velocity and the mean square velocity are almost equal to the two dimensional PIV,which proves that LES is credible.

The turbulence anisotropy of CIJR is calculated using three methods.By comparing differentanisotropy parameters,the anisotropy ofthe impinging area and edge of impinging jet in CIJR was found to be stronger and the maximum deviation was 40%.These results show that the isotropic hypothesis used in the k–ε model and two dimensional PIV is notapplicable to the description ofturbulentkinetic energy distribution in CIJR,which results in an error of 20%to 30%for the estimation of turbulence intensity.When using 2D-PIVresults for TKE calculation,different correction factors(γ)need to be adopted for different areas of concern,which can help improve the accuracy of TKE calculation.

The nature of the stagnation point fluctuation could assist in choosing the optimal mixing chamber width for designing the CIJR,reducing the reactor volume and ensuring that the impact will always occur in the mixing chamber.The mean velocity and TKE distribution could help the design of the internal shape of CIJR,such as manufacture the flow channel which anastomoses with TKE distribution to avoid dead zones.The results about turbulence anisotropy give rise to provide advice of TKE calculation method for the design of CIJR.

Nomenclature

A Axisymmetric parameter

d Inner diameter of jet,mm

dpDiameter of the particle,μm

G Two-dimensional parameter

I Turbulence intensity

II-III Trace of the anisotropic matrix

k Turbulent kinetic energy,m2·s-2

l Inner length of jet,mm

Re Reynolds number in the jets

rmsU Root Mean Square velocity in the x-direction,m·s-1

rmsV Root Mean Square velocity in the y-direction,m·s-1

rmsW Root Mean Square velocity in the z-direction,m·s-1

s–t Eigenvalues of the anisotropic matrix

U Mean velocity in the x-direction,m·s-1

u Mean velocity magnitude of jet,m·s-1

V Mean velocity in the y-direction,m·s-1

W Mean velocity in the z-direction,m·s-1

γ Correction factor for TKE calculation

νfKinematic viscosity,m2·s-1

ρpDensity of the particle,kg·m-3

ρfDensity of flow,kg·m-3