Guangchun Song,Yuxing Li*,Wuchang Wang*,Kai Jiang,Zhengzhuo Shi,Shupeng Yao

Shandong Key Laboratory of Oil-Gas Storage and Transportation Safety,College of Pipeline and Civil Engineering,China University of Petroleum,Qingdao 266580,China

Keywords:Hydrate Agglomeration Flow behavior Dynamic model Numerical simulation Population balance

ABSTRACT Dynamic modeling and numerical simulation of hydrate slurry flow behavior are of great importance to offshore hydrate management.For this purpose,a dynamic model of hydrate agglomeration was proposed in this paper.Based on population balance equation,the frame of the dynamic model was established first,which took both hydrate agglomeration and hydrate breakage into consideration.Then,the calculating methods of four key parameters involved in the dynamic model were given according to hydrate agglomeration dynamics.The four key parameters are collision frequency,agglomeration efficiency,breakage frequency and the size distribution of sub particles resulting from particle breakage.After the whole dynamic model was built,it was combined with several traditional solid–liquid flow models and then together solved by the CFD software FLUENT 14.5.Finally,using this method,the in fluences of flow rate and hydrate volume fraction on hydrate particle size distribution,hydrate volume concentration distribution and pipeline pressure drop were simulated and analyzed.

1.Introduction

Naturalgas hydrates(NGH)are crystalline solids composed ofwater and gas molecules such as methane,ethane,propane,and carbon dioxide[1].NGH are easy to form when the ambient temperature is relatively low and the ambient pressure relatively high[2].Since the first explicit hydrate plugging incident was identified in 1934[3],NGH have always been found during the low-temperature,high-pressure process of oil exploitation and transportation[4].In addition,with the development tendency of oil industry moving towards deep sea and ultra-deep sea,hazards caused by hydrate plugging are now posing a severe threat to the subsea flow assurance[5,6].The average expense for subsea hydrate prevention is approximately 1000000 USD per mile[7].Considering such situation,many investigations have been conducted on hydrate prevention strategies.Among these strategies,thermodynamic inhibition and kinetic inhibition are most commonly used both in laboratory scale and industrial scale.Thermodynamic inhibition prevents hydrate formation by injecting the thermodynamic hydrate inhibitors(THIs),which can shift the hydrate equilibrium curve to a lower temperature and higher pressure condition[8].Instead of THIs,kinetic inhibition uses kinetic hydrate inhibitors(KHIs)and antiagglomerants(AAs)to slow down the nucleation and growth process ofhydrate particles and to preventthe hydrate particles fromagglomerating,respectively[9,10].KHI and AA are collectively called the low dosage hydrate inhibitors(LDHI).Compared with THI,LDHI has the advantages of being economical and eco-friendly.More detailed information about LDHI can be found elsewhere[11–13].

The research on hydrate agglomeration is of great importance to plug prevention in offshore operations.Actually,researchers have performed a great deal of work on hydrate agglomeration phenomenon analysis,hydrate agglomeration mechanism interpretation and hydrate cohesive force measurement and made many achievements on it[14–17].However,much work still needs to be conducted to investigate the dynamic model of hydrate agglomeration.Based on the balance between shear force and cohesive force,Muhle[18]proposed a model to calculate the maximum agglomerate diameter for aggregates in the laminar flow.Using the viscosity of hydrate suspension instead of the dispersing liquid,Camargo[19]improved the model proposed by Muhle.Then,combining hydrate rheological theory with the improved model,Camargo established a phenomenological model and applied it into the calculation of hydrate slurry apparent viscosity.However,in the application of the phenomenological model,a constant cohesion force between hydrate particles was used and no explicit calculation formula for the cohesion force was given.Through a comprehensive analysis of the forces acting on hydrate particles,Wang[20]identified that capillary liquid bridge force is the main cohesive force between hydrate particles and modified the model of Camargo accordingly.Although the models above can be used to analyze hydrate agglomeration,hydrate agglomeration efficiency and the processofhydrate breakage are notfully considered.Therefore,in order to better describe actual production process,Colombel[21]proposed a preliminary dynamic model for hydrate agglomeration.This model took population balance equation as a frame and two agglomeration mechanisms were considered in it,which were known as contact-induced agglomeration mechanism[22]and shear-limited agglomeration mechanism.However,in that paper,Colombel did not give the calculating methods of the key parameters appeared in her model.The same as Colombel,Balakin[23]also built a hydrate agglomeration model based on the population balance equation and introduced the calculating methods of the key parameters.Validated with experimental data,Balakin's model showed a good accuracy.This indicates that hydrate agglomeration modeling using the population balance theory is feasible.

The numerical simulation of hydrate agglomeration and hydrate slurry flow behavior in pipeline also play an important role in hydrate plug prevention,especially when experimentalsetups and experimental conditions are not available.Although there are lots of literatures on the numericalsimulation of solid–liquid flow[24,25],relative literatures especially on hydrate slurry flow are deficient.In early stage,hydrate slurry flow was simulated on the condition of constant hydrate particle size and hydrate agglomeration was not considered,just as the simulation of other traditional solid–liquid flow.Under these circumstances,flow pressure drop and the distribution of hydrate volume concentration were focused on.Then,with the development of hydrate agglomeration dynamics,changes in hydrate particle size were considered by using the model of Muhle and Camargo,respectively[26].Therefore,hydrate particle size distribution could be simulated accordingly.However,as discussed before,hydrate agglomeration efficiency and hydrate breakage could not be analyzed with the simulation methods above.Under these circumstances,a combination of computational fluid dynamics(CFD)and the population balance model(PBM)can be helpful and Balakin[27,28]have got some achievements through this way.

In the presentwork,a dynamic modelofhydrate agglomeration was proposed.Based on population balance equation,the frame of the dynamic model was established first,which took both hydrate agglomeration and hydrate breakage into consideration.Then,the calculating methods of four key parameters involved in the dynamic model were given according to hydrate agglomeration dynamics.The four key parameters are collision frequency,agglomeration efficiency,breakage frequency and the size distribution of sub particles resulting from particle breakage.After the whole dynamic modelwas built,itwas combined with some traditional solid–liquid flow models and together solved by the CFDsoftware FLUENT 14.5.Finally,using this method,the in fluences of flow rate and hydrate volume fraction on hydrate particle size distribution,hydrate volume concentration distribution and pipeline pressure drop were simulated and analyzed.

2.Hydrate Agglomeration Modeling

2.1.Population balance model

At the very beginning,population balance model is mainly used to describe the dynamic balance of population.After the introduction of Hulburt[29],itbegan to be utilized in chemical field[30,31].In chemical applications,PBM associates the microscopic behavior(nucleation,growth,agglomeration,breakage)of the discrete phase with its macroscopic properties(volume,area,particle size)by tracking the quantity density function variation of it.At present,population balance model has been widely used in the chemical processes of crystallization,flocculation,polymerization,granulation and multiphase flow.

According to the work of Ramkrishna[32],when diffusion together with nucleation and growth of the discrete phase is neglected and agglomeration and breakage are considered,PBM can be expressed as follows:

where,n(ξ，x，t)and u(ξ，x，t)are the quantity density function and velocity function of the discrete phase,respectively.ξ is a variable characterizing a macroscopic property(volume,area,particle size)of the discrete phase,x is the spatial variable and t is the time variable.aβ(ξ-ξ′，ξ′)is the agglomeration kernel of the discrete phase,where β represents the collision frequency and a represents the agglomeration efficiency after collision.b(ξ|ξ′)S(ξ′)is the breakage kernel of the discrete phase,where S(ξ′)represents the breakage frequency and b(ξ|ξ′)represents the size distribution of the discrete phase resulting from breakage.

In order to build the dynamic model for hydrate agglomeration,it is assumed that the particle size of hydrate particles is continuously distributed and the in fluences ofconvection and diffusion are negligible.Therefore,when the discrete phase is hydrate particle,Eq.(1)can be simplified to:

where,L and L′represent the size of hydrate particles.

The physical meaning of Eq.(2)can be expressed as the change in quantity density at moment t caused by hydrate agglomeration and hydrate breakage for hydrate particles whose particle size is L.According to the effect of agglomeration and breakage,the change in quantity density above can be subdivided into the four parts below:①I(mǎi)ncreasement resulting from the agglomeration between two hydrate particles with size L-L′and L′(L′is smaller than L here),respectively.The coefficient 1/2 is used to prevent the same agglomeration from being counted twice.②Decreasement resulting from the agglomeration between two hydrate particles with size L and L′(L′can be any size here),respectively. ③Increasement caused by the formation of hydrate particles with size L,where the formation results from the breakage of hydrate particles with size L′(L′is larger than L here).④ Decreasement caused by the breakage of hydrate particles with size L.

Eq.(2)is the frame of the dynamic modelfor hydrate agglomeration and it consists of two parts,the agglomeration kernel and the breakage kernel.For agglomeration kernel,it includes collision frequency and agglomeration efficiency.For breakage kernel,it includes breakage frequency and the size distribution of sub particles resulting from particle breakage.The calculating methods of the two parts are given below.

2.2.Agglomeration kernel

The agglomeration of hydrate particles in the flow field can be divided into two steps.First,two hydrate particles come close to each other and then collide under the action of agglomeration driving force.Second,two particles overcome the hydrodynamic forces in the flow field and complete the agglomeration.These two steps can be described by collision frequency and agglomeration efficiency,respectively.

The calculating method of collision frequency is investigated first.According to different agglomeration driving forces,there are mainly three different mechanisms for particle collision,Brownian motion(perikinetic),differential sedimentation and flow shear(orthokinetic).Detailed information of the three different mechanisms can be seen in Table 1[33].Combine Table 1 with hydrate properties,it can be seen that all the three collision mechanisms can lead to hydrate particle collision in the pipeline.For the collision caused by Brownian motion,the calculation formula of collision frequency proposed by Smoluchowski[34]has the highest recognition,which is also suitable for hydrate particles and is shown as follows:

where,kBis the Boltzmann constant,T is the temperature and μ is the dynamic viscosity of the fluid.

For the collision caused by differential sedimentation,the formula proposed by Camp[35]is frequently used to calculate the collision frequency:

where,V is the settling velocity,which can be obtained by Eq.(5):

As for the collision resulting from flow shear,lots of researches have been done and many different calculation formulas of the collision frequency have been proposed.In this paper,in order to improve the calculation speed of the dynamic model,the formula established by Camp[35]is finally selected to calculate the collision frequency of hydrate particles in the laminar flow,as shown below:

where,G is a local parameter referred to as the absolute velocity gradient and is defined as follows:

where,ν is the kinematic viscosity of the fluid andε is the energy dissipation rate.

The formulas established by Saffman[36]and Abrahamson[37]are selected to calculate the collision frequency of hydrate particles in the turbulent flow.First,in the turbulent dissipation region,the sizes of hydrate particles are smaller than the Kolmogorov microscale.Therefore,hydrate agglomeration in this region is dominated by the local shear in vortices.Under these circumstances,the collision frequency of hydrate particles can be calculated using the formula established by Saffman[36]:

Second,in the turbulent inertia region where the sizes of hydrate particles are larger than the Kolmogorov microscale,the motion of hydrate particles is driven by the main flow field.Under these circumstances,the collision frequency of hydrate particles can be calculated using the formula established by Abrahamson[37]:

As the collision mechanisms for Brownian motion,differential sedimentation and flow shear are different,the collision frequency of each collision mechanism should be independent of each other.Therefore,the actual collision frequency of hydrate particles in the pipeline can be calculated by:

A comparison between the calculated values ofcollision frequency is given below.Considering that the collision frequency in turbulent flow is usually higher than that in the laminar flow,so Eq.(8)is used to calculate the collision frequency of flow shear in the comparison.According to hydrate properties and actual offshore operating conditions,the size of hydrate particles used in the calculation ranges from 5 μm to 200 μm[38]and the temperature and velocity gradient used in the calculation are set as 276.15 K and 300 s-1,respectively.The calculation and comparison results are shown in Fig.1.

As seen in Fig.1a–c,under the calculation conditions of this paper,the collision frequency of flow shear is the largest while the collision frequency of Brownian motion is the smallest and the difference between them can be several orders of magnitude.From Fig.1a–c it can also be seen that,for the collision caused by Brownian motion and differential sedimentation,the collision frequency is large enough only when the difference between the two particle sizes is large enough.However,when two hydrate particles have similar particle sizes,the collision frequency of them caused by Brownian motion and differential sedimentation is very close to zero.In addition,as shown in Fig.1d,the actual collision frequency of hydrate particles is mainly in fluenced by the collision frequency of flow shear.Especially when the sizes of hydrate particles are smaller than 500 μm,the values of βTOLand βSLare almost exactly the same.

Therefore,the collision frequency caused by Brownian motion is neglected in this paper and the collision frequency used in our dynamic model is set as:

Table 1 Three different mechanisms for particle collision

Fig.1.Calculation and comparison results of the collision frequency.

The calculating method of agglomeration efficiency is studied next.At present,there are mainly two models used for the calculation of agglomeration efficiency,rectilinear model and curvilinear model[39].For the rectilinear model,it does not take flow fields and the shortrange forces between approaching particles into account.Therefore,the agglomeration efficiency in this model is 1,which is always larger than the actual value.Different from the rectilinear model, flow fields and the short-range forces between approaching particles are fully considered in the curvilinear model.Considering this,the curvilinear model is selected and used in our dynamic model to calculate the agglomeration efficiency between hydrate particles.

Eq.(12)is a typicalcurvilinear modelproposed by Ven[40],which is used in our dynamic model.

where,k is a parameter related to fluid properties,ε is a parameter which represents the ratio between van der Waals force and flow shear force and can be obtained as follows:

where,H is the Hamaker constant characterizing van der Waals force,R is the harmonic radius of the two colliding particles.

From the discussion above it can be seen that,the agglomeration efficiency of two colliding particles is closely related to the forces acting on them.Furthermore,the forces acting on hydrate particles are dominated by the liquid phase in the pipeline.When the liquid phase in pipeline is water,the main cohesion force acting on hydrate particles is van der Waals force.However,when the liquid phase changes to oil–water mixture,the main cohesion force would then change to the capillary liquid bridge force.Therefore,Eq.(13)can only be used when the liquid phase in pipeline is water or the operating conditions accord with cold flow technology[41].Under these circumstances,a new parameterε′is proposed in this paper based on the ratio between capillary liquid bridge force and flow shear force and it is used to replace the parameter ε in Eq.(12).The calculation formula of ε′is as follows:

where,γ is the oil–water interfacial tension,μ′is the dynamic viscosity of oil–water mixture,θ is the contact angle of liquid bridge on hydrate particle surface and γ cos θ represents the capillary liquid bridge force.With this new parameter,the agglomeration efficiency of hydrate particles in oil–water mixture can be calculated then.

Based on Eqs.(12)–(15),a comparison between the agglomeration efficiency in water and the agglomeration efficiency in oil–water mixture is given below.The parameters and their corresponding values[39–41]used in the comparison are shown in Table 2.It should be noticed that,when the calculated value of a is larger than 1,just take a=1.

Table 2 Parameters and their corresponding values used in the comparison of agglomeration efficiency

Fig.2 shows the calculation and comparison results of the agglomeration efficiency.As seen in Fig.2,the variation tendency of the agglomeration efficiency in water is almost the same with that in oil–water mixture,that is,small hydrate particles usually correspond to large agglomeration efficiencies.It can also be seen from Fig.2c that,under the calculation conditions of this paper,the agglomeration efficiency ofhydrate particles in oil–water mixture(0.7–0.9)is far larger than that in water(0.02–0.05).This indicates that,compared with water system,the oil–water system has a higher tendency of hydrate agglomeration and plugging.The same conclusion has also been obtained by many other researchers[42].

In addition,through Eqs.(12)–(15),the critical harmonic radius Rcrwhich makes the agglomeration efficiency equal to 1 can be calculated.When the harmonic radius of two hydrate particles is smaller than Rcr,the agglomeration of these two hydrate particles is then destined to happen after collision.Under the calculation conditions of this paper,the critical harmonic radii for the water system and oil–water system are 0.067 μm and 9.591 μm,respectively.However,the initial size of hydrate particles formed in offshore operating conditions usually ranges from 0.2 to 50 μm[43].Therefore,once hydrate particles are formed in oil–water systems,they may agglomerate with each other quickly and dramatically,which would lead to a quick increase in particle size and pose a great threat to pipeline flow assurance.

Fig.2.Calculation and comparison results of the agglomeration efficiency.

To sumup,in ourdynamic model,Eqs.(12)–(13)are used to calculate the agglomeration efficiency in water systems and Eq.(12)together with Eq.(15)is used to calculate the agglomeration efficiency in oil–water systems.

2.3.Breakage kernel

Breakage kernel consists of breakage frequency and the size distribution of sub particles resulting from particle breakage.In the flow field,there are mainly two mechanisms for particle breakage,flow shear and particle collision.According to the work of Serra[44],the calculation formula of breakage frequency for flow shear is as follows:

And the calculation formula of breakage frequency for particle collision is given as:

where,E is a proportionality constant and m is a constant inversely proportional to the aggregate strength of hydrate particles.Breakage will happen only when the flow shear stress overcomes the aggregate strength.Both E and m can be obtained by experimental data fitting.K is a parameter represents for the possibility of breakage after particle collision.

However,according to the work of Maggi[45],the value of K is very small and it works only when the particle concentration is very high.In addition,there are always liquid bridges between hydrate particles in oil–water systems,which hinder the direct contacts between two hydrate particles during the collision process.Therefore,hydrate particles can hardly deform and break in oil–water systems due to particle collision.Under these circumstances,particle breakage caused by particle collision is neglected in this paper and Eq.(16)(m=1.90,E=800 s0.90·m-1)[23]is used to calculate the breakage frequency of hydrate particles in our dynamic model.

For the size distribution of sub particles resulting from particle breakage,it can be described mainly by binary distribution,ternary distribution and Gaussian distribution.However,according to the work of Zhang[46],the three distributions above have very little effect on the final results of particle size distribution simulation.Therefore,in order to reduce calculation,binary distribution is selected to calculate the size distribution of the sub particles in our dynamic model.The calculation formula of binary distribution is as follows:

3.Numerical Models

After the dynamic model of hydrate agglomeration is established,it is combined with some traditional solid–liquid flow models and then together solved by the CFD software FLUENT 14.5.The main numerical models used in the solving process are illustrated below.

3.1.Geometric model

In this paper,a set of freon(R11)hydrate slurry flow experiments conducted by Balakin[47]are selected as the object of our numerical simulation and the experimental results are used to verify our simulation results.In order to simulate hydrate agglomeration and hydrate slurry flow precisely,a three-dimensional geometric model of the horizontal pipeline used in Balakin's experiments is established,as shown in Fig.3.

Fig.3.The three-dimensional geometric model and its mesh structure used in the simulation.

The same as Balakin's experimental flow loop,the diameter of the horizontal pipeline model is set to 45.2 mm and a length of 3 m is selected in order to get a fully developed hydrate slurry flow.When generating the mesh,hexahedrons and a step length of 1 mm are used.In the near wall region,8 layers of mesh are generated to deal with the boundary layer effect.In general,the mesh of the geometric model consists of 267696 hexahedrons and the mesh quality is 0.913.

Fig.4 shows the mesh dependency ofthe geometric model.As can be seen in Fig.4,the velocity gradients in the near wall region for two different hexahedron numbers are almost the same,with a relative deviation of 6.3%.This indicates that the mesh dependency of the geometric model is good enough for numerical simulation.Therefore,in order to reduce calculation,the mesh with 267696 hexahedrons is finally selected.

Fig.4.Flow rate profiles of the continuous phase on a midline cross-section ofthe pipeline for two different hexahedron numbers.Mean flow rate is 2 m·s-1.

3.2.Multiphase model and viscous model

In this paper,multiphase model and viscous model are two main numerical models used in the simulation and the numerical simulation processes are based on the following assumptions:①The flowing process is isothermal and no interphase mass transfer is considered.②The formation and dissociation of hydrates are neglected.③The hydrate slurry used in the simulation consists of water phase and hydrate phase and no oil phase and gas phase are considered.The numerical simulation of hydrate slurry flow using oil–water mixture will be investigated in another paper of our research group.④All the hydrate particles are continuous medium.⑤The process of hydratewall adhesion is neglected.

For multiphase model,the Eulerian–Eulerian two fluid model is selected and it consists of two governing equations(continuity equations and momentum equations)and several constitutive equations used for the closure of the equation set.According to the assumptions above,the continuity equations(Eq.(19))and momentum equations(Eq.(20))in the multiphase model are written as follows:

where,i represents the water phase or hydrate phase,ρ is the density,α is the volume fraction,? is the Laplasse operator,u is the velocity vector,p is the pressure,τ is the stress tensor and M is the interphase momentum exchange.

The liquid–solid coupling ofhydrate slurry mustbe focused on when establishing the multiphase model.In numerical simulations,liquid–solid coupling is achieved by the interphase momentum exchange.In this paper,the interphase momentum exchange M is calculated mainly considering the interphase drag force Mdand the turbulent diffusion force Mt.Then,an interphase force model is established accordingly and it is used for the closure of the equation set,as can be seen below:

where,l and s represent water and hydrate respectively,uris the interphase relative velocity,μt,mis the turbulent viscosity,σdis the Planck diffusion coefficient and klsis the momentum transfer coefficient.

In this paper,klsis calculated by the model of Ding[48].For this model,when αs≤ 20%,it can be written as:

When αs＞ 20%,it can be written as:

where,CDis the drag coefficient.

Besides the interphase force model,the dynamic viscosity ofhydrate phase is also necessary for the closure of the multiphase model,which can be obtained by the following equation:

where,μmis the dynamic viscosity of the hydrate slurry and it can be calculated by Roscoe–Brinkman equation[49],given as Eq.(25).

In the experimental work of Balakin[47],which is the object of our simulation,the apparent viscosity of R11 hydrate slurry can be well described by Eq.(25)under the conditions in Table 4.Therefore,although the in fluence of hydrate particle diameter on the dynamic viscosity of hydrate slurry[50]is not considered in Eq.(25),Eq.(25)is still selected considering its high accuracy in the simulation scenario of this paper and for model simplification.

As for the viscous model,the standard k–ε model is selected and standard wall functions are used to deal with the near wall treatment.More detailed information of this part can be found elsewhere[27]and no repetition will be listed here.

3.3.Model solution

In this paper,the dynamic model of hydrate agglomeration together with the multiphase model and the viscosity model are all solved using the CFD software FLUENT 14.5.As the frame of PBM and solid–liquid flow models is originally inserted in FLUENT 14.5,only parameter settings and callings according to Eqs.(4)–(18)and(24)–(25)are needed before solving the models.For boundary conditions,the inlet boundary condition of the geometric model is set as the velocity inlet and the outlet boundary condition of the geometric model is set as the pressure outlet.Besides,a second order upwind difference scheme is used in momentum equation discretization and a phase coupled SIMPLE scheme is used in pressure–velocity coupling.In addition,PBM is solved by the discrete method in this paper.When the maximum residual value is no more than 1×10-5,the model solution process ends.

The key parameters involved in the solution process and their corresponding values[23,26,39]are listed in Table 3.

Table 3 Parameters involved in the solution process and their corresponding values

4.Numerical Simulation

4.1.Simulation conditions

In this paper,the numerical simulation mainly focuses on the in fluences of flow rate and hydrate volume fraction on hydrate particle size distribution,hydrate volume concentration distribution and pressure drop.In order to compare with the experimental results,the flow rates and hydrate volume fractions used in the simulations are all from the work of Balakin[47],as seen in Table 4.

Table 4 Simulation condition list

Forcases 1–4,they are used to investigate the in fluences of flow rate.For cases 4–7,they are used to analyze the in fluences ofhydrate volume fraction.

In addition,the initial size distribution of hydrate particles should be given when using PBM.In order to be comparable,the size distribution shown in Table 5 is used for all the seven cases.That is to say,in each case,the sizes of all the hydrate particles in pipeline are 5 μm before simulation.

According to the work of Cheng[51],particle size distribution is sensitive to the number of bins when the discrete method is used to solve the PBM.Therefore,the analysis of the bin number independence is necessary here.Taking case 4 as an example,Fig.5 illustrates the particle size distribution ofcase 4 simulated with differentbin numbers.As seen in Fig.5,hydrate particle size distributions are almost the same when bin number is larger than 4.Therefore,bin number 7 is selected in this work in order to reduce calculation.

Table 5 Initial size distribution used in the simulation

Fig.5.Particle size distribution ofcase 4 simulated with different bin numbers.(Minimum size is 5 μm and maximum size is 640 μm).

4.2.Model validation

According to the experiments of Balakin[47],the experimental data of pressure drop and hydrate particle size distribution are selected to verify the simulation results.

Table 6 shows the comparison of pressure gradient between experimental results and simulation results.As seen in Table 6,the variation tendency of pressure gradient for simulation cases 5–7 is almost the same with that for corresponding experiments.In addition,it can be seen that the relative errors between experimental results and simulation results are all smaller than 18%.Therefore,the simulation method and numerical models used in this paper can well describe the pressure drop characteristics of hydrate slurry flow.

Table 6 Comparison of pressure gradient between experimental results and simulation results

Fig.6.Comparison of hydrate particle size distribution between simulation results(cases 5–7)and experimental results.

Fig.5 gives the comparison of hydrate particle size distribution between simulation results(cases 5–7)and experimental results.It can be seen from Fig.6 that the law of hydrate particle size distribution obtained from numerical simulation approximates log-normal distribution and is almost the same with that obtained from corresponding experiments.For both simulation results and experimental results,hydrate particles with a particle size smaller than 50 μm account for the largest proportion and hydrate particles whose particle sizes are larger than 150 μm account for a very small proportion.Thus,the simulation method and numerical models used in this paper can also well describe the particle size distribution of hydrate slurry.

As for the error of hydrate particle size distribution between simulation results and experimentalresults,itis considered to be derived from the following aspects.Firstly,in experiments,the sampling processes were stochastic and discontinuous.Therefore,the samples could not well re flect the real situation of hydrate particle size distribution in the pipeline.Secondly,for both experiments and simulations,the hydrate particle size distribution was calculated using the data of several cases(cases 5–7)rather than one single case.Thus,the simulation results and experimental results are thought to be better matching for one single case.

4.3.In fluences of flow rate

The in fluences of flow rate on hydrate particle size distribution in the outlet radial cross section of the pipeline are analyzed first.

Fig.7 shows the hydrate particle size distribution in the outlet radial cross section atdifferent flow rates when hydrate volume fraction in the pipeline is 10%.As can be seen from Fig.7,in the radial cross section of the pipeline,hydrate particles in the near wall region possess larger diameters while hydrate particles in the center region possess smaller diameters.According to hydrate agglomeration dynamics, flow shear is the main reason for the collision and agglomeration of hydrate particles.In a pipeline,the flow shear in the near wall region is the strongest.Therefore,hydrate particles in the near wall region are more likely to collide and agglomerate with each other,which further leads to the increase of particle diameter.However,in the center region of a pipeline, flow shear is relatively weak and hydrate particles are homogeneously distributed(as seen in Fig.8).Thus,the probability ofhydrate particle collision and agglomeration is relative small.The conclusion above is contrary to the results of the simulations using the model of Muhle and Camargo.The explanations are as follows.Firstly,in the simulations using the model of Muhle and Camargo,hydrate cohesion force at different places of the pipeline was assumed to be a same constant.Secondly,hydrate agglomeration efficiency was neglected when using the model of Muhle and Camargo.Therefore,hydrate particles would agglomerate right after formation and the agglomeration process would notstop untilthe maximum particle diameter is reached.However,these two circumstances would never happen in actual production processes.Therefore,the simulation results in this paper are more believable.

Fig.7.Contours of hydrate particle size(in m)distribution in the outlet radial cross section of the pipeline(cases 1–4).

Fig.8.Log-normal fittings of the data on hydrate particle size distribution for cases 1–4.

Table 7 lists the characteristic hydrate particle sizes in each simulation case.Combine Table 7 with Fig.7 it can be seen that,in the radial cross section of the pipeline,the proportion of hydrate particles with large diameters increases with the increase of the flow rate and so is the mean particle diameterin pipeline.This can be illustrated asfollows.As the flow rate increases, flow shear in the pipeline also increases.Thus,hydrate particles are more likely to agglomerate and large particle diameters are more easily to form.In addition,it can also be concluded from Table 7 and Fig.7 that,the maximum hydrate particle diameter gradually decreases with the increase of flow rate.This is because the strong flow shear resulting from high flow rates reduces the maximum diameter which hydrate particles could maintain.Besides,in simulation processes,the initial hydrate particle diameter(5 μm)is set to be the minimum diameter which could not be further broken.Therefore,the minimum hydrate particle diameter is about 5 μm in all the seven cases.

Table 7 Characteristic hydrate particle sizes in each case

Probability distribution function is frequently used when describing the distribution of hydrate particle size,such as Gaussian distribution,log-normal distribution,Gamma distribution and Weibull distribution.The fittings of the data on hydrate particle size distribution for the seven cases indicate that only log-normal distribution can well describe the hydrate particle size distribution in this paper,with the determination coefficients all above 0.98.This conclusion accords with the research results of Balakin[47]and Clarke[52].Fig.7 shows the lognormal fittings of the data on hydrate particle size distribution for cases 1–4.As seen in Fig.8,hydrate particle size distribution in the flowing system is well fitted by the log-normal distribution.In addition,it can also be seen from Fig.8 that the range of hydrate particle size distribution gradually decreases from case 4 to case 1.This indicates that at the same hydrate volume fraction,the maximum hydrate particle diameter gradually decreases with the increase of flow rate.This conclusion is the same as the one obtained by analyzing Fig.6 and Table 7.

The in fluences of flow rate on hydrate volume concentration distribution in the outlet radial cross section of the pipeline is analyzed next.

Fig.9 shows the hydrate volume concentration distribution in the outlet radial cross section at different flow rates when hydrate volume fraction in the pipeline is 10%.As seen in Fig.9,two kinds of solid–liquid flow patterns,homogeneous suspension flow(case 4)and heterogeneous suspension flow(cases 1–3),can be observed in the simulation processes.According to the analysis before,hydrate particle diameters in case 4 are relatively small and the mean diameter is about 39 μm.In this case,the solid–liquid flow pattern in the pipeline is homogeneous suspension flow,where hydrate volume concentration in the near wall region is relatively high and hydrate volume concentration in the center region is relatively low.This can be illustrated by the difference in flow rate.In the center region of the pipeline,the flow rate is higher than that in the near wall region and so is the dispersion coefficient of hydrate particles.Therefore,hydrate particles can be dispersed homogeneously in the center region.However,with the increase of flow rate,the mean diameter of hydrate particles in the pipeline grows sharply.The proportion of hydrate particles whose diameters are larger than 150μm also increases greatly.Underthese circumstances(cases 1–3),the flow conditions in the pipeline could not maintain the homogeneous suspension of hydrate particles any more.Thus,the homogeneous suspension flow in pipeline gradually turns into heterogeneous suspension flow.Considering that the density of R11 hydrate is higher than that of water,so hydrate particles tend to concentrate at the bottom of the pipeline in heterogeneous suspension flow.On the contrary,hydrate volume concentration at the top of the pipeline is relatively low.In addition,it can also be concluded from Fig.9 that the hydrate volume concentration at the bottom of the pipeline gradually decreases with the increase of flow rate.This phenomenon can also be explained by the strong ability of carrying and dispersing resulting from high flow rates.

Fig.9.Contours of hydrate volume concentration distribution in the outlet radial cross section of the pipeline(cases 1–4).

Finally,the in fluences of flow rate on pressure drop ofhydrate slurry are given below.

Table 8 shows the pressure gradient in each simulation case.Compare the pressure gradient of cases 1–4,it can be found that the pressure gradient increases with the increase of flow rate.According to the discussion before,the increase ofpressure drop is notonly related to flow rate increase itself,but also related to the change of the flow pattern.After the flow pattern converts into heterogeneous suspension flow,the heterogeneous distribution of hydrate particles would lead to the increase of flow friction,which could cause a further increase of the pressure drop.

Table 8 Pressure gradient in each simulation case

4.4.In fluences of hydrate volume fraction

Fig.10 shows the hydrate particle size distribution in the outlet radial cross section at different hydrate volume fractions when the flow rate in pipeline is 1.5 m·s-1.The same as cases 1–4,in the radial cross section of the pipeline,hydrate particles in the near wall region possess larger diameters while hydrate particles in the center region possess smaller diameters.In addition,combine Fig.10 and Table 7 it can be seen that in the radial cross section of the pipeline,the proportion of hydrate particles with large diameters increases with the increase of hydrate volume fraction and so are the mean particle diameterand maximum particle diameter in pipeline.This is because hydrate particles are more likely to collide and agglomerate with each other at high hydrate volume fractions.

Fig.11 shows the log-normal fittings of the data on hydrate particle size distribution for cases 4–7.The same as cases 1–4,hydrate particle size distribution in cases 4–7 can also be well described by log-normal distribution.Besides,the range of hydrate particle size distribution gradually increases from case 4 to case 7.This indicates that at the same flow rate,the maximum hydrate particle diameter gradually increases with the increase of hydrate volume fraction.This conclusion is the same as the one obtained by analyzing Fig.10 and Table 7.

Fig.12 shows the hydrate volume concentration distribution in the outlet radial cross section at different hydrate volume fractions when the flow rate in pipeline is 1.5 m·s-1.As seen in Fig.12,only homogeneous suspension flow is observed in cases 4–7.The same as case 4,hydrate volume concentration in the near wall region is relatively high and hydrate volume concentration in the center region is relatively low.In addition,with the increase ofhydrate volume fraction,the range of the low hydrate volume concentration region in the pipeline center gradually decreases.However,the corresponding concentration value of the low hydrate volume concentration region gradually increases.This indicates that at the same flow rate,the carrying and dispersing ability resulting from liquid flow gradually decreases with the increase of hydrate volume fraction.

Fig.10.Contours of hydrate particle size(in m)distribution in the outlet radial cross section of the pipeline(cases 4–7).

Fig.11.Log-normal fittings of the data on hydrate particle size distribution for cases 4–7.

As for pressure gradient,from Table 8 it can be concluded that the pressure gradient increases with the increase of hydrate volume fraction.This is due to the increase of flow friction caused by hydrate volume fraction increasing.In addition,it can also be seen from Table 8 that the in fluence of hydrate volume fraction on pressure gradient is relatively small than that of flow rate.This indicates that flow rate is the key parameter which effects pressure gradient most in this paper.

5.Conclusions

A dynamic model of hydrate agglomeration was proposed based on the population balance theory and was then used to simulate the in fluences of flow rate and hydrate volume fraction on hydrate particle size distribution,hydrate volume concentration distribution and pipeline pressure drop.

Fig.12.Contours of hydrate volume concentration distribution in the outlet radial cross section of the pipeline(cases 4–7).

When building the dynamic model,the conclusions below were revealed.Population balance model could be used to describe the agglomeration process of hydrate particles.Brownian motion could be neglected when calculating the collision frequency of hydrate particles and only flow shear and differential sedimentation should be mainly considered.A typical curvilinear model was used to calculate the agglomeration efficiency of hydrate particles in water systems.Based on the ratio between capillary liquid bridge force and flow shear force,the curvilinear model above was modified and then used to calculate the agglomeration efficiency of hydrate particles in oil–water systems.Flow shear was focused on when calculating the breakage frequency of hydrate particles and binary distribution was used to describe the size distribution of sub hydrate particles resulting from particle breakage.

The conclusions obtained by analyzing the numerical simulation results are as follows.Hydrate particle size distribution in the pipeline can be well described by log-normal distribution.In the radial cross section of the pipeline,hydrate particles in the near wall region possess larger diameters while hydrate particles in the center region possess smaller diameters.At the same hydrate volume fraction and in the radial cross section of the pipeline,the proportion of hydrate particles with large diameters increases with the increase of the flow rate and so is the mean particle diameter in pipeline.However,the maximum hydrate particle diameter gradually decreases with the increase of flow rate.At the same flow rate and in the radial cross section of the pipeline,the proportion of hydrate particles with large diameters increases with the increase of hydrate volume fraction and so are the mean particle diameter and maximum particle diameter in pipeline.However,the carrying and dispersing ability resulting from liquid flow gradually decreases with the increase of hydrate volume fraction.Compared with hydrate volume fraction, flow rate is the key parameter which effects pressure gradient most in our simulations.