Numerical Study on Cavitation Generation in a Water Jet Added with

Jul 9, 2018 - For some engineering projects in submerged conditions such as ... Some theories such as mechanical effects theory, chemical effects theo...
0 downloads 0 Views 1MB Size
Subscriber access provided by UNIV OF DURHAM

Thermodynamics, Transport, and Fluid Mechanics

Numerical study on cavitation generation in water jet added with turbulent drag-reducing additives based on improved cavitation model and Cross Model Guo-Dong Li, Song-Sheng Deng, Jin-Fa Guan, and Yan Chen Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b02179 • Publication Date (Web): 09 Jul 2018 Downloaded from http://pubs.acs.org on July 10, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Numerical study on cavitation generation in water jet added with turbulent drag-reducing additives based on improved cavitation model and Cross Model Guo-Dong Li*, Song-Sheng Deng, Jin-Fa Guan,Yan Chen Department of Fuels, Army Logistics University of PLA, Chongqing 401331, China Abstract An improved cavitation model was proposed to simulate the cavitation phenomenon in water jet added with turbulent drag-reducing additives. The Cross model which expressed the shear-thinning properties of drag-reducer solution was also employed. The neglected viscous term and surface tension term in traditional cavitation models were reconsidered as viscosity and surface tension force were important factors affecting bubbles’ motion in drag-reducer solution. The numerical solutions of bubble radius in drag-reducer solutions of different concentrations were fitted and the change rates were added to the original mass transfer mechanism equations. The results showed the modified cavitation model described the cavitation characteristics in drag-reducer solutions more accurately and effectively reflected that the cavitation region’s wake became longer with the increase concentration of the drag reducing agent. With the concentration increasing, the distribution of cavitation cloud became more homogeneous and the pressure curve and velocity curve both changed more smoothly. Keywords: water jet, cavitation model, drag-reducer, numerical simulation 1 Introduction Abrasive water jet technology has been widely used in chemical, petroleum and mining industries for its better cutting performance compared with traditional machining methods in safety, heat damage and thermal stress. For some engineering projects in submerged conditions such as underwater salvage and the remove of ocean platforms, especially for some dangerous occasions like the rescue of sunken oil tankers, chemical tankers or nuclear submarines, abrasive water jet is always a preferred approach. Franz (1970) firstly pointed out that the depth penetration of a water jet would be increased due to additives, the prelude of the research on additives solution jet was opened. Since then, scholars such as Kudin, Zakin, Summerds, Howells, Annoni M 1,Wang Rui-He, Yang Yong-Yin and Zhang You-Yi 1-7 conducted valuable studies on the application performance of water jets added with polymer or surfactant. Shimizu 8 did some experimental studies on the erosion characteristics of premixed abrasive water jet under submerged conditions and found that the erosion characteristics were greatly affected by the cavitation phenomenon. Some theories such as mechanical effects theory, chemical effects theory, electrochemical effects theory and heat effects theory were gradually proposed to explain the damage mechanism of cavitation erosion. The viscoelastic properties of the solution of polymer or surfactant will significantly affect the cavitation characteristics in water jet and will further affect its cutting ability and the morphological features of the cutting section. From this perspective, it’s necessary to make a detailed study on the cavitation phenomenon in viscoelastic fluids. In 1968 9, Ellis et al. observed that the addition of a small amount of water-soluble polymer Corresponding author. E-mail address: [email protected] (Guo-Dong Li)

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

led to the weakening of cavitation to some extent, and this phenomenon attracted people's interest in studying the cavitation characteristics in non-Newtonian fluids. Later, some researchers pointed out that the addition of polymer additives weakened the cavitation intensity in water jet or rotating disc flow and the cavitation damages to equipment were then significantly affected 10-13. Brujan et al.14 (1996) conducted experimental study on bubbles’ behavior in viscoelastic fluids with shear-thinning properties (e.g. CMC, PAM solution). In their conclusions, the expansion and compression processes of the bubbles would not be obviously influenced when bubble radius was more than 0.5 mm and obvious extension of oscillation time would be observed when bubble radius was less than 0.5mm. Since some types of constitutive models were gradually proposed to describe the viscoelastic properties of non-Newtonian fluids, modern investigations of bubble dynamics in non-Newtonian fluids are always combined with an appropriate constitutive model such as Oldroyd model15, Carreau model16, Jeffreys model17 or Maxwell model18. Some investigations on bubbles’ behavior in non-Newtonian fluids are also valuable and meaningful. For example, Brujan 19 (1998) proposed the Williamson model to characterize the rheological properties of aqueous polymer solutions such as CMC and HEC in which the fluid elasticity is neglected. He concluded that when the maximum radius of the bubble was less than 0.1mm, the shear-thinning characteristics strongly affected the behavior of bubbles, and the addition of polymer reduced the velocity of the cavity walls and prolonged the collapse time. In Albernaz and Cunha’s numerical investigation20, Maxwell's constitutive model and fifth order Runge-Kutta method were used to study the nonlinear oscillation of the bubbles on the basis of considering the anisotropy of the additive solution. Conclusions that the effect of additives on bubbles showed up as the resistance of extensional viscosity in its longitudinal motion were finally got. Malvar et al. 21 (2016) conducted numerical research on the oscillations of sonic-excited spherical bubbles in non-Newtonian fluids and concluded that the key factors affecting the bubbles’ behavior were surface tension and viscosity of the liquid. When it comes to the numerical simulation of cavitation by CFD, the Singhal model 22, the Schnerr-Sauer model 23 and the Zwart–Gerber–Belamri Model 24 are available in Fluent. The Singhal model is based on the “full cavitation model” and can be used to account for the effect of non-condensable gases. The Schnerr and Sauer model regards the mixture of liquid phase and vapor phase as a mixture of a large number of spherical vapor bubbles and expresses the volume fraction of the vapor as a function of the number of bubbles in the liquid per unit volume. The Zwart-Gerber-Belamri model takes the bubble density numbers into consideration and further optimizes the mass change rate for a single bubble. The above three cavitation models have good practicality for most engineering problems and there are still some other researchers proposing some modified cavitation models by proposing different mechanism in mass transfers (e.g. Hosangadi and Ahuja, Morgut, Congedo et al and Asnaghi et al 25-27 ). The selection of cavitation models hinges on the main influencing factor that must be taken into consideration in engineering problems (e.g. viscosity, surface tension, shear stress or temperature), and those unimportant factors are often neglected. In this study, by considering the viscosity’s effects and surface tension’s effects of the non-Newtonian fluids on bubbles’ growth and collapse, an improved cavitation model was proposed to simulate the cavitation in water jet. The Cross viscosity model is also employed to reflect the shear-thinning characteristics of the liquid.

2

ACS Paragon Plus Environment

Page 2 of 18

Page 3 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

2 Simulation models 2.1 The improvement of cavitation model Generally Rayleigh-Plesset equation which is often abbreviated as R-P equation has the following form. Pb − P∞

ρL

2

= R⋅

d 2 R 3  dR  4vL dR 2 S + ⋅ ⋅ +  + R dt ρ L R dt 2 2  dt 

(1)

Here Pb denotes the bubble surface pressure, P∞ does the local far-field pressure, ρL represents the liquid density, νL stands for the liquid kinematic viscosity, S represents the liquid surface tension coefficient, R is bubble radius. As we have mentioned above, in traditional cavitation models, Rayleigh-Plesset equation is usually simplified into the following form ignoring the viscous term and surface tension force term. 2 P∞ − Pb 3 ρL

dR n = ( -1) dt

(2)

The parameter n is equal to1 during bubble’s expansion and 2 for compression. For the viscoelastic solution formed after the addition of drag-reducer, the viscosity and surface tension of the surrounding liquid medium are non-negligible factors affecting cavitation properties. Especially in the study of the cavitation phenomenon in drag-reducer solution of different concentrations, the effects of viscosity are particularly important. In these cases, the original Rayleigh-Plesset equation containing viscosity term and surface tension term needs to be considered. In Eq. (1), if the variables R and t are transformed into the non-dimensional forms as follows, β=

R n t ,τ = ( -1) R0 R0

P∞ − Pb

ρL

(3)

R0 is the initial radius of a single bubble, also the parameter n is equal to1 during bubble’s expansion and 2 for compression, in addition, if we define parameters C and D as follows, C= ( -1)

n

4µL R0 ρ L P∞ − Pb

,D= ( -1)

n

S R0 ⋅ P∞ − Pb

(4)

Eq. (1) can be transformed into the following form.

β⋅

d 2 β 3 d β 2 C d β 2D + ( ) + ⋅ + ±1 = 0 dτ 2 2 dτ β dτ β

(5)

The sign “+” is chosen during bubbles’ compression and “-” for expansion. Noting that

P∞ − Pb ρ L represents the scale of velocity. If we define U =

then the definition of parameter C can be written as follows. n 4ν L n 4 C = ( -1)

R0U

= ( -1)

Re

P∞ − Pb ρ L ,

(6)

Where ν L =µ L / ρ L represents the kinematic viscosity of surrounding liquid, Re = R0U ν L stands for Reynolds number which means the ratio of inertial forces to viscous forces. Similarly parameter D can also be expressed in the following form in which the Weber number We is included. 3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

D= ( -1)

n

S n 1 = ( -1) R0 ⋅ P∞ − Pb We

Page 4 of 18

(7)

The Weber number represents the ratio of inertia forces to surface tension forces and its definition is given below. We =

ρ LU c 2 R0 S

(8)

Then Eq. (1) can be written in the following form in which Reynolds number and Weber number are included. 2

β⋅

d 2β 3  d β  4 1 dβ 2 1 +  ⋅ ⋅ + ⋅ ±1 = 0  + dτ 2 2  dτ  Re β dτ β We

(9)

Parameter C depends on fluid viscosity and D depends on surface tension. In some engineering problems, if the viscosity term or surface tension term is neglected to simplify calculating, the above equation may not be appropriate for there is no corresponding Re or We exists. To obtain the numerical solution of Eq. (9), the corresponding Re and We values should be given. The values of Re and We are determined by the viscosity and surface tension coefficients of the drag reducer solution. According to Jiang Chen-Xing et al. 28 and Li Feng-Chen 29 et al.'s experimental investigations, we know that the 200ppm CTAC (cetyl trimethyl ammonium chloride) solution has a zero shear viscosity of 0.005Pa • s, a relaxation time of 0.08s and a surface tension coefficient of 0.0322N / m. CTAC solution of 500ppm, its zero shear viscosity is 0.008Pa • s, relaxation time is 0.2s, surface tension coefficient is 0.0336N / m. CTAC solution of 1000ppm, its zero shear viscosity is 0.01Pa • s, relaxation time is 1s, the surface tension coefficient is 0.036N / m. The surface tension coefficient of pure water is 0.072N / m. On the basis of the above data, the change rates of dimensionless radius with dimensionless time in CTAC solution of different concentrations (200 ppm, 500 ppm, and 1000 ppm) can be obtained by programming calculation. A Cash-Karp fourth-fifth order Runge-Kutta method with degree four interplant was used to find the numerical solution. The liquid pressure was defined by applying the boundary conditions in the bubble interface. The input data of the simulations were the dimensionless quantities: C (Re), D (We), β0 (initial radius), dβ/dτ|τ=0 (initial velocity) and the parameters which control interactions and the adaptive time steps. The numerical solutions are shown in Figure 1.

1-a

4

ACS Paragon Plus Environment

Page 5 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

1-b Figure 1. Numerical solutions of primitive Rayleigh-Plesset in drag-reducer solution of different concentrations,1-a was for the expansion process and 1-b for the compression process.

From the curves it can be seen that, if the effects of viscosity and surface tension are taken into consideration, obvious differences may be found in the variation of bubble radius in drag-reducer solution and pure water. As the concentration of drag-reducer increases, the curves become smoother, which means that the bubble radius changes more slowly. The default cavitation models in ANSYS Fluent are based on simplified Eq. (2) where the changes of the radius of a spherical bubble are mainly influenced by the pressure difference inner and outer the bubble. The change rate of bubble radius is used to describe the liquid-vapor mass transfer mechanism in transport equation. The vapor transport equation is written as follows. r ∂ (αρv ) + ∇ ⋅ ( ∂ρv uv ) = Re − Rc ∂t

(10)

Re denotes the mass transfer source term during bubbles’ growth or expansion. Rc is the mass transfer source term during bubbles’ collapse or compression. The continuity equation of the mixture model which is also employed in our simulation is shown in Eq. (11). The momentum conversation equation is given in Eq. (12). ∂ρm r + ∇ ⋅ ( ρ m um ) = 0 ∂t

(11)

∂ r r r r r ( ρ m um ) + ∇ ⋅ ( ρ m um um ) = −∇P + ∇ ⋅  µm ( ∇um + ∇umT )  ∂t

(12)

Where ρm, um and µm are the density, velocity and dynamic viscosity of the mixture respectively. Combining Eq. (11) with mass conservation equations of each respective phase in the mixture, the following equation which describes the relations between ρm and α is finally obtained. d ρm dα = − ( ρl − ρ v ) dt dt

(13)

Here α denotes the volume fraction of vapor. By substituting α with an expression with variables R (bubble radius) and n (bubble number density), Eq. (14) can be finally obtained. 1 2 d ρm dR = − ( ρl − ρv )( n4π ) 3 ( 3α ) 3 dt dt

5

ACS Paragon Plus Environment

(14)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 18

The traditional operation of building cavitation model is to substitute dR/dt with the simplification

2 P∞ − Pb 3ρ L

. However, noting that if the previous numerical solutions of R-P

equations are linearly fitted, there exists the following relationship between β and τ. d β dτ = k

(15)

k is the slope of the fitted line. According to the relationship between β and R,as well as the relationship between τ and t, the following equation is easily obtained. dR d β n = ⋅ ( -1) dt dτ

P∞ − Pb

ρL

(16)

Then by substituting dβ/dτ with k, the change rate of bubble radius can be written as follows. P∞ − Pb

dR n = ( -1) ⋅ k dt

ρL

(17)

We can also write the above equation into the following forms considering the expansion process (Eq.18) and the compression process (Eq.19) respectively. Pb − P∞  dR    = ke dt ρL  e

(18)

P∞ − Pb  dR    = kc dt ρL  c

(19)

The subscript e is the abbreviation of expansion and c is the abbreviation of compression. Similar to the traditional cavitation model, the above expressions can be written as follows. 2 Pb − P∞  dR  ⋅   = 1.2247ke 3 ρL  dt e

(20)

2 P∞ − Pb  dR  ⋅   = 1.2247kc 3 ρL  dt c

(21)

By replacing the traditional simplification of Eq.(2) with Eq.(20) and Eq.(21) in Zwart-Gerber-Belamri cavitation model, we can finally obtain the new expressions of mass transfer rates in Zwart-Gerber-Belamri model as follows. Re =Cme ⋅ Fvap

3α nuc (1 − α v ) ρv

Rc =Cmc ⋅ Fcond

ℜB

2 ( Pv − P ) ( P < Pv ) 3 ρl

2 ( P − Pv ) ( P > Pv ) 3 ρl

3α v ρ v ℜB

(22)

(23)

Pv=3540 Pa is the saturated vapor pressure at temperature of 300K, αnuc=5×10-4 denotes the nucleation site volume fraction,Fvap=50 is the evaporation coefficient and Fcond=0.01 is the condensation coefficient. Cmc=1.2247kc is given the name of modified compression coefficient and 6

ACS Paragon Plus Environment

Page 7 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Cme=1.2247ke is given the name of modified expansion coefficient. The values of calculated Cmc and Cme are given in Table 1.

Media properties kc ke Cmc Cme

water

200ppm solution 500ppm solution

1000ppm solution

µ=0.001Pa.s, σ=0.072N/m

µ=0.005Pa.s, σ=0.0322N/m

µ=0.008Pa.s, σ=0.0336N/m

µ=0.01Pa.s, σ=0.036N/m

0.1213 0.7283 0.1486 0.8920

0.0237 0.3343 0.0290 0.4094

0.0167 0.1481 0.0205 0.1814

0.0141 0.0899 0.0173 0.1101

Table 1. The values of Cmc and Cme in drag-reducer solution of different concentrations.

For the three default cavitation models in ANSYS Fluent, Zwart-Gerber-Belamri model and Schnerr and Sauer model have better stability and robustness than Singhal et al. model. The Zwart-Gerber-Belamri model and Schnerr and Sauer model are also compatible with all the turbulence models available in ANSYS Fluent. The difference between Zwart-Gerber-Belamri model and Schnerr and Sauer model is only in the density terms in the mass transfer rate. In Zwart-Gerber-Belamri model, the unit volume mass transfer rate is only related to the vapor phase density. Unlike Schnerr and Sauer model, mass transfer rate has no relation with the liquid phase and mixture densities in this model. Besides, the Zwart-Gerber-Belamri model already contains coefficients Fvap and Fcond , so we can realize the modification either by UDFs(user-defned functions ) or just by changing the coefficients via constant pre-multiplier in the Zwart-Gerber-Belrami model. That’s why we chose the Zwart-Gerber-Belrami model as an example. 2.2 Cross model Both the polymer solution and the surfactant solution are viscoelastic fluids and the numerical simulation of viscoelastic fluids needs constitutive equations that can accurately describe the flow characteristics. The commonly used constitutive equations include the generalized Newtonian fluid models, the linear viscoelastic fluid models and the nonlinear viscoelastic fluid models. Four generalized Newtonian fluid models are available in Fluent, they are the power law model, the Carreau model, the Cross model and the Herschel-Bulkley model, and all of them can well characterize the shear-thinning characteristics of viscoelastic fluid. The linear viscoelastic fluid models can describe the unsteady flow behavior of some viscoelastic fluids, but they can not describe the characteristics of viscosity changing with shear rate. The typical linear viscoelastic constitutive models are the Maxwell model and the Jeffreys model. The commonly used nonlinear viscoelastic fluid models include differential viscoelastic fluid models (such as Oldroyd-B model, Giesekus model, FENE-P model and PTT model) and integral fluid models (such as Lodge model, K-BKZ model and Rivlin Sawyers model), both the two kinds of models have good generality and can well describe the flow characteristics of viscoelastic fluids. If the nonlinear viscoelastic constitutive models are employed for numerical simulation in Fluent, the UDS function should firstly be opened to calculate the molecular deformation rate tensor of the drag-reducer solution, and the calculated results will be added into N-S equation as a 7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

source term for further calculation. However, both the differential constitutive models and the integral constitutive models can not express the stress as an explicit function of velocity or velocity gradient when the Weissenberg number reaches a certain value. This may lead to a high degree of nonlinearity in numerical calculation 30. So the solving difficulty and computational load of the nonlinear viscoelastic constitutive models are significantly higher than the generalized Newtonian fluid model. According to the analysis above, the Cross model which is one of the generalized Newtonian fluid models was chosen for the cavitation simulation in drag reducer solution. The expression is as follows. m µ = µ0 1+ ( λγ& ) 





(24) & γ µ0 denotes the zero shear viscosity, λ is the characteristic time, is the shear rate, m denotes the empirical index. λ=0.4µsM[µ]2c/R0T, M denotes the average molar mass of additive molecules, µs denotes the dynamic viscosity of the solvent, c is the solute concentration, R=8.31J/mol•K, T is the absolute temperature, [µ]=1.03×10-2M0.78 denotes the intrinsic viscosity. Jiang Chen-Xing 31 measured the zero-shear viscosity and characteristic time of 200ppm, 500ppm and 1000ppm PAM (polyacrylamide) solution and CTAC / NaSal solution and concluded that it was feasible to use Cross model to characterize the viscosity features of PAM and CTAC solution. The calculated results generally matched well with the measured values. However, there is a deviation between the experimental value and the calculated value at the 1000ppm concentration when the shear rate is large. According to Jiang’s explanation in his doctoral thesis, this may be caused by the error in the process of CTAC solution preparation or measurement error in the experiment. In Jiang’s published papers, the Cross model has been proved to be feasible to describe the viscosity characteristics of the drag-reducer solutions mentioned above. 3 Simulation and analysis 3.1 Geometry, grid, turbulence model and boundary conditions Based on the improved cavitation model and Cross viscosity model, the cavitation phenomenon in CTAC solution of the above three conditions issuing from a nozzle was numerically simulated. The geometry of the nozzle is shown in Figure 2 and the parameters are given as D/d= 2.88, L/d=4 where D, d, L denote inlet diameter, the nozzle diameter and nozzle length, respectively.

Figure 2. The diagram of the geometric structure of

Figure 3. The mesh of computational domain

the nozzle

2D structure grid in axisymmetric coordinates was chosen for the simulation. An overview of the generated grid is shown in Figure 3. The Coupled algorithm was used for pressure-velocity 8

ACS Paragon Plus Environment

Page 8 of 18

Page 9 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

coupling and a PRESTO! scheme was used for pressure interpolation. Discretization schemes for momentum and volume fraction were selected as QUICK. Convergence criteria were equal to 10-6 for all cases. The choice of turbulence model will play a significant role in predictions of cavitation phenomenon. Coutier-Delgosha O et al. (2003) 32 studied the influence of different models (standard k-ε RNG; modified k-ε RNG; k-w with and without compressibility effects) on the numerical simulation of cavitating flows. They concluded that the standard models k-ε RNG and k-w without compressibility effects poorly described the vapor cloud shedding behavior. Their modified k-ε RNG model and the k-w model including compressibility effects reliably described Venturi unsteady cavitation behavior. In our investigation, we aimed to use our improved cavitation model to study the effects of liquid viscosity and surface tension on the stabilized cavitation morphology. The phenomenon of cloud shedding during the development is not within the scope of this study. Therefore, for ease of calculation and convergence, we chose realizable k-ε model as the turbulence model. The pressure inlet was set 600kPa and the pressure outlet was set 200kPa.The temperature was set 300K. UDF function in Fluent was employed to embed the modified coefficients into the traditional Zwart-Gerber-Belamri model to get the improved model described in Eq.(22) and Eq.(23). 3.2 Vapor volume fraction Numerical simulations of 200ppm, 500ppm and 1000ppm CTAC solutions were carried out based on traditional cavitation model (the Zwart-Gerber-Belamri model), the combination of Cross constitutive model and Zwart-Gerber-Belamri model, the combination of Cross constitutive model and improved Zwart-Gerber-Belamri model, respectively. The contours of vapor volume fraction simulated by different models are shown in Figure 4.

Zwart-Gerber-Belamri model

Cross model and Zwart-Gerber-Belamri model

Cross model and improved Zwart-Gerber-Belamri model

a

b

c

200 ppm

500 ppm

1000 ppm

Figure 4. Contours of vapor volume fraction for three kinds of additive solutions simulated by different models

From the simulation results in Figure 4-a, it can be seen that the scale of the cavitation region decreases as the concentration increases from 200ppm to 1000ppm. Results in Figure 4-a can not 9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

match the changing trend of the cavitation area in Jiang Chen-Xing’s experimental results shown in Figure 5, for the scale of the cavitation region increases as the concentration increases in the experiments. So, conclusions can be obtained that it is not appropriate to simulate the cavitation phenomenon in drag-reducer solution just based on cavitation model. The simulation results based on the Cross constitutive model and the original cavitation model can make up for this problem (shown in Figure 4-b). In Figure 4-b, the cavity lengths at each concentration are 18.24mm, 19.20mm and 20.35mm respectively. A significantly increasing trend of the cavitation region scale was showed as the concentration increases, which can partially match the experiments by Jiang. But the phenomenon that the cavity’s wake becomes longer with the increase of concentration in Figure 5 was not reflected in Figure 4-b. However, the simulation results based on the combination of the Cross constitutive model and the improved cavitation model in Figure 4-c predict the cavity wake’s characteristics very well. As is shown in Figure 4-c, the lengths of the whole cavitation region at each concentration are 25.72mm, 26.21mm and 26.55mm respectively. What’s more, the lengths of the area where the vapor volume fraction is less than 0.2 at each concentration are 1.31mm, 1.93mm and 2.21mm respectively. It is obvious that the wake of cavitation cloud gradually becomes longer with the increase of concentration. The wake of the supercavity shown in Figure 5 has a length of 2.3 cm, 4.2cm and 5.2cm at 3 different concentrations. Comparing the two sets of data, we find that the simulation results are consistent with experimental phenomena to some extent in predicting the morphology of the cavitation wake. Here we define the cavitation area where the vapor volume fraction is less than 0.2 as the weak cavitation region (shown in Figure 6), and the variations of the length of weak cavitation region with concentration are shown in Figure 7.

Figure 6. Schematic of the so-called weak cavitation region

Figure 5. Experimental results of supercavity in drag-reducer solution of different concentrations

10

ACS Paragon Plus Environment

Page 10 of 18

Page 11 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

30

0.1 0.08

20

0.06 0.04

10

0.02 0

0 200ppm

500ppm

1000ppm

The length of the whole cavitation region The length of the weak cavitation region The ratio of the two lengths Figure 7. Variations of the scales of weak cavitation region and the whole cavitation area as the concentration of additive changes

As is shown in Figure 7, the length of the weak cavitation region accounts for an increasing proportion of that of the whole cavitation region with the increase of concentration. At the same time, as shown in Figure 4-c, the red area which represents a region of high vapor volume fraction gradually becomes smaller, but the area where the vapor volume fraction is low becomes relatively larger. Combined with the legend, conclusions may be drawn that cavitation clouds distribute more uniformly with increase in the concentration. The improved cavitation model in which the effects of the fluid’s viscosity and surface tension on the bubbles’ growth and collapse are taken into consideration is proved to be more accurate to predict the cavitation characteristics in drag-reducer solution. The simulation results based on our modified cavitation model and Cross model successfully predicted the increase in the length of cavitation wake as the concentration of CTAC solution increases, which is in good agreement with the experimental results shown in Figure 5. Actually, when the viscosity of surrounding fluid increases, the maximum bubble radius decreases while the minimum bubble radius increases. This is because the rise of liquid viscosity increases the resistance of bubbles’ oscillation and increases the loss of bubble energy 21. Moreover, the more viscous the liquid becomes, the faster the bubble energy is consumed and the momentum diffused. According to the definition of Reynolds number, when Re decreases or C increases, such effects of viscosity will become more distinct in bubbles’ oscillation. Some scholars pointed out in their experimental research that surface tension force at the bubble wall slows down the bubble growth process but speeds up the collapse process33. What’s more, the influence of surface tension on bubble’s expansion is more than that on collapse. Increasing surface tension will decrease the maximum and minimum bubble radius and shorten the oscillation times. With the increase of the concentration of the CTAC solution, the inhibitory effects of the liquid viscosity and surface tension on bubbles’ motion are more significant, which microscopically leads to the reduction of cavitation intensity. As shown in Figure 4, non-zero values of vapor volume fraction were predicted close to the maximum nozzle length, and as such gradients in the solution variables in the normal direction will persist at the outlet boundary if its position coincides with the nozzle length. So it’s necessary to conduct a detail investigation on the cavitation in the area near the outlet of the nozzle. As a 11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 18

result, water jets issuing from the nozzle into submerged conditions were simulated and a pressure inlet of 5000kPa was selected. The transient development of cavitation is shown in the Figure 8 and the stable states are shown in Figure 9.

Figure 8. The transient development of cavitating flows

Pure water 200 ppm

500 ppm

1000 ppm

Figure 9. The stabilized state of cavitation in submerged conditions at inlet of 5000kPa based on our improved cavitation model

Obvious cavitation area can be seen at the nozzle exit due to the shear effect during the development of cavitating water jet. At the given pressure, the stabilized cavitation regions extend beyond the nozzle into submerged conditions and reach a length of 46.08 mm, 50.688 mm and 57.60m respectively at each concentration. At higher pressure drop the velocity across the periphery of nozzle would be higher and thus more cavities are generated which ultimately results into the higher vapor volume fraction. Some studies have also shown the inlet pressure is a significant factor that affects the intensity of cavitation and the length of the cavitation zone34-35. Here inlet pressures of 500kPa, 600kpa and 700kPa were chosen to see the effects of inlet pressure on the length of the cavitation zone within the nozzle. The concentration of the drag reducer solution selected here was 500ppm. Also the improved cavitation model and Cross model were employed. The results are shown in Figure 10.

500kPa

600kPa 12

ACS Paragon Plus Environment

700kPa

Page 13 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 10. Contours of vapor volume fraction for three inlet pressures by Cross model and improved cavitation model

It can be seen from the simulation results that as the inlet pressure increases, the length of the entire cavitation region increases significantly. When the inlet pressure reaches 700kPa, the length of the cavitation zone has exceeded the nozzle. Figure 5 appears to show a separated wake region downstream of the vapor cavity, which is an unsteady phenomenon. Actually cavitation phenomenon along with the submerged water jet bears unique unsteady characteristics. Subsequent images taken by Hai-Xia Liu et al. provide a good proof of this phenomenon36. As for the cavitation inception mechanism, available knowledge is far from sufficient. According to the investigation by Manolis Gavaises 37 on the cavitation cloud formation and collapse in an axisymmetric geometry based on Visualisation and Large Eddy Simulations (LES) and traditional cavitation model, the cavities formed at the primary inception point at the turn are unstable and the cavity detachment was well reflected by LES simulation. This research and the study by Coutier-Delgosha, O. mentioned above remind us to combine our improved cavitation model and LES or other suitable turbulence model in the simulation of unsteady cavitation phenomenon in our future research. 3.3 Pressure and velocity distribution The changes of the pressure and velocity of the flow fields are plotted in Figure 11 and Figure 12.

Figure 11. The curves of pressure changes at the

Figure 12. The curves of velocity changes at the

location of Y=0.035m

location of Y=0.035m

T stands for the traditional cavitation model, C represents the cross constitutive model, and I represents the improved cavitation model. As can be seen from the figures, the pressure and the velocity both show a gradual slowing trend along the nozzle with the increase of the CTAC concentration. This phenomenon is particularly obvious for the change of velocity. As is shown in Figure 12, the calculation results based on Cross model and the improved cavitation model show that the curvature of the velocity curve of 1000ppm is obviously less than that of 200ppm. It means that the velocity becomes less fluctuating along the axis as the concentration of CTAC increases. From another point of view, it is confirmed that the increase of concentration leads to 13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the increase of viscosity, and then leads to a greater effect on the momentum diffusion.

200 ppm

500 ppm

1000 ppm

Figure 13. Contours of pressure for three kinds of additive solutions simulated by Cross model and improved cavitation model

The distribution of pressure for all the three concentrations is shown in Figure 13. It can be seen that the length of the transition zone of pressure near the exit of the nozzle gradually increases with the increase of the concentration. This coincides with the change trends of the weak cavitation zone at different concentrations. The velocity is an important factor for the application in water jet. The changes of jet velocity at different concentrations were reflected in Figure 12. Figure 14 shows the effects of different inlet pressures on the velocity distribution in the nozzle.

500 kPa

600 kPa

700 kPa

Figure 14. Contours of velocity for three inlet pressures simulated by Cross model and improved cavitation model

14

ACS Paragon Plus Environment

Page 14 of 18

Page 15 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

It can be seen from the above figures that the velocity within the nozzle significantly increased with the increase in inlet pressure. The velocity at the exit of the nozzle under the three pressures is 18.04m/s, 19.71m/s and 22.20m/s, respectively. 4 conclusions 1 Taking into consideration the effects of liquid viscosity and surface tension on bubbles’ growth and collapse processes, the original cavitation model was improved. The simulation results based on our improved cavitation model and Cross constitutive model are proved to be more accurate to reflect the cavitation characteristics in the drag-reducer solution. 2 The entire cavitation area and the weak cavitation area both expand with increase in the concentration. The weak cavitation region accounts for an increasing proportion of the whole cavitation region. 3 With the increase of the concentration, the pressure and the velocity both showed a gradual slowing trend along the nozzle. Acknowledgements This work is supported by the National Natural Science Foundation of China (Grant Nos. 51475469) Abbreviations

Pb : pressure within the bubble P∞: external pressure far from bubble ρL: the density of surrounding liquid νL: Kinematic viscosity of surrounding liquid,νL=µL/ρL µL: dynamic viscosity of surrounding liquid µ0 : the zero shear viscosity S: surface tension of bubble R: radius of bubble R0: the initial radius of bubble. t: time ρm: the density of mixture um: the velocity of mixture µm: the viscosity of mixture α : the volume fraction of vapor phase ρv: the density of vapor phase uv: the velocity of vapor phase Re: mass transfer source term connected to the growth of the vapor bubbles Rc: mass transfer source term connected to the collapse of the vapor bubbles n: bubble number density Fvap: evaporation coefficient Fcond: condensation coefficient ℜ B : bubble radius in mass transfer equation in Fluent Compliance with ethical standards 15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Conflict of interest: The authors declare that they have no conflict of interest. Reference (1) Kudin, A. M.; Barenblatt, G. I.; Kalashnikov, V. N.; etc. Destruction of Metallic Obstacles by a Jet of Dilute Polymer Solution[J]. Nature (London), Phys. Sci., 1973, 245,95-96. (2)Zakin,J.L.; Summers, D.A. The Effect of Viscoelastic Additives on Jet Structures[J]. Proceedings of the Third International Symposium on Jet Cutting Technology, BHRA Group, Jan 1976. (3) Howells, W.G. Super-water jetting applications from 1974 to 1999,Proc.10th America Water jet Conference,14th-17th August, Houston,USA. (4) Annoni, M.1. The role of polymeric additives in water jet cutting[J]. International Journal of Machining & Machinability of Materials, 2009, 6. (5) Rui-he, Wang; Shu-jiang, Wang; Yu-huan, Bu. Structural features of jet with high molecular polymer additive [J]. J. China Univ. Pet. 1998, 22,38-40. (6) Yong-yin, Yang; Rui-he, Wang; Wei-dong, Zhou. Study on the Effect of PAM on Abrasive Jet Flow[J]. Pet. Drill. Tech.2001, 29,4-6. (7) You-yi, Zhang; Xiao-chuan,Wang; Dong, Hu. Influence of Adding a Small Amount High Molecular Polymer on Cutting Performance of Abrasive Water Jet[J]. J. Sichuan Univ. 2014, 46,181-187. (8) Shimizu,S. Erosion due to Premixed Abrasive Water Jet under Submerged Condition[J]. Nihon Kikai Gakkai Ronbunshu B Hen/transactions of the Japan Society of Mechanical Engineers Part B, 1993, 59,2964-2968. (9) Ellis, A.T.; Hoyt, W.J.; Robertson, J.M. Some Effects of Macromolecules on Cavitation Inception, ASME Cavitation Forum, 1968, pp.2–3. (10) Hoyt, J.W. Effect of Polymer Additives on Jet Cavitation[J].J. Fluids Eng.1976, 98,106-112. (11) Ting, R.Y. Characteristics of flow cavitation in dilute solutions of polyethylene oxide and polyacrylamide[J]. Phys. Fluids. 1978, 21,898-901. (12) Shima, A.; Tsujino,T.; Nanjo, H.; etc. Cavitation Damage in Polymer Aqueous Solutions[J]. J. Fluids Eng. 2014, 107,134-138. (13) Tsujino,T.; Shima,A.; Nanjo,H. Effects of Various Polymer Additives on Cavitation Damage[J]. Archive Proceedings of the Institution of Mechanical Engineers Part C Mechanical Engineering Science 1983-1988 (vols 197-202), 1986, 200,231-235. (14) Brujan,E.A.; Ohl,C.D.; Lauterborn,W.; etc. Dynamics of Laser-Induced Cavitation Bubbles in Polymer Solutions[J]. Acta Acust. Acust. 1996, 82,423-430.

(15) Oldroyd, J.G. On the Formulation of Rheological Equations of State[J]. Proc. R. Soc. London.1950, 200,523-541. (16) Carreau,P.J.; Macdonald,I.F.; Bird, R.B. A nonlinear viscoelastic model for polymer solutions and melts—II[J]. Chem. Eng. Sci.1968, 23,901-911. (17) Shima, A.; Tsujino,T.; Oikawa, Y. The collapse of bubbles in viscoelastic fluids (the case of Jeffreys model fluid), Rep. Inst. High Speed Mech. Tohoku Univ. 1988,55, 1–15. (18) Kim, C. Collapse of spherical bubbles in Maxwell fluids[J]. J. Non-Newtonian Fluid Mech.1994, 55,37-58. (19) Brujan, E.A. Bubble dynamics in a compressible shear-thinning liquid[J]. Fluid Dynamics Research, 1998, 23,291-318. 16

ACS Paragon Plus Environment

Page 16 of 18

Page 17 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(20) Cunha, F. R.; Albernaz, D.L. Oscillatory motion of a spherical bubble in a non-Newtonian fluid[J]. J. Non-Newtonian Fluid Mech.2013, 191,35-44. (21) Malvar, S.; Gontijo, R.G.; Cunha, F.R.Pattern identifcation on the nonlinear radial motion of an oscillating spherical bubble using neural networks[J]. Lat Am J Solids Struct.2016,13,2464-2489. (22) Singhal, A. K.; Athavale,M.M.; H, Li.; etc. Mathematical Basis and Validation of the Full Cavitation Model[J]. J. Fluids Eng.2002, 124,617-624. (23) Schnerr, G. H.; Sauer, J. Physical and numerical modeling of unsteady cavitation dynamics[C]. In Fourth International Conference on Multiphase Flow, New Orleans, USA. 2001. (24) Zwart, P. J.; Gerber, A. G.; Belamri, T. A two-phase flow model for predicting cavitation dynamics[C]. In Fifth International Conference on Multiphase Flow, Yokohama, Japan. 2004. (25) Hosangadi, A.; Ahuja, V. A. New unsteady model for dense cloud cavitation in cryogenic fluids[C].ASME 2005 Fluids Engineering Division Summer Meeting. 2005,623-629. (26) Congedo, P. M.; Goncalves, E.; Rodio, M. G. Uncertainty quantification of turbulence and cavitation models in cavitating flows simulations[J]. European Journal of Mechanics B/Fluids, 2015, 53,190-204. (27) Asnaghi, A.; Feymark, A.; Bensow, R E. Improvement of cavitation mass transfer modeling based on local flow properties[J]. Int. J. Multiphase Flow, 2017,93,142–157. (28) Jiang C X; Shuai Z J; Zhang X Y; etc. Numerical study on evolution of axisymmetric natural supercavitation influenced by turbulent drag-reducing additives[J]. Appl. Therm. Eng. 2016, 107,797-803. (29) Li F C; Dong Y; Kawaguchi Y; etc. Experimental study on swirling flow of dilute surfactant solution with deformed free-surface[J]. Exp. Therm. Fluid Sci. 2009, 33,161-168. (30) Zhang M; Zhang L; Jiang B. Progress in the numerical simulation for agitated flow field of non-Newtonian fluid[J]. Chem. Ind. Eng. Prog. (Beijing, China), 2009, 28,1296-1301. (31) Jiang C X. Study on the mechanism of the supercavitating flow influenced by turbulent drag-reducing additives [D].Harbin: Harbin Institute of Technology,2015. (32) Coutier-Delgosha, O.; Reboud, J. L.; Delannoy, Y. Numerical simulation of the unsteady behaviour of cavitating flows[J]. Int. J. Numer. Methods Fluids.2003, 42,527-548. (33) He Jie; Liu Xiumei; Lu Jian; etc. Surface Tension Effects on the Growth and Collapse of Cavitation Bubbles near a Rigid Boundary[J]. Chin. J. Lasers (Chin. Ed.). 2009, 36,342-346. (34) Bashir, T. A.; Soni, A. G.; Mahulkar, A. V.; etc. The CFD driven optimisation of a modified venturi for cavitational activity[J]. Can.J.Chem.Eng. 2011, 89,1366-1375. (35) Kuldeep, Saharan,V. K.; Computational study of different venturi and orifice type hydrodynamic cavitating devices[J]. Journal of Hydrodynamics, Ser. B. 2016,28, 293-305. (36) Liu H; Kang C; Zhang W; etc. Flow structures and cavitation in submerged waterjet at high jet pressure[J]. Exp. Therm. Fluid Sci.2017, 88. (37) Gavaises, M.; Villa, F.; Koukouvinis, P.; etc. Visualisation and les simulation of cavitation cloud formation and collapse in an axisymmetric geometry[J]. Int. J. Multiphase Flow, 2015, 68,14-26.

17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

For Table of Contents Only

18

ACS Paragon Plus Environment

Page 18 of 18