Zai-Sha Mao ,Chao Yang2

1 CAS Key Laboratory of Green Process and Engineering,Institute of Process Engineering,Chinese Academy of Sciences,Beijing 100190,China

2 School of Chemical Engineering,University of Chinese Academy of Sciences,Beijing 100049,China

Keywords:Total drag coefficient Virtual mass force Solid particle Numerical simulation Vorticity-stream function formulation

ABSTRACT Virtual mass force is an indispensable component in the momentum balance involved with dispersed particles in a multiphase system.In this work the accelerating motion of a single solid particle is mathematically formulated and solved using the vorticity-stream function formulation in an orthogonal curvilinear coordinate system.The total drag coefficient was evaluated from the numerical simulation in a range of the Reynolds number (Re) from 10 to 200 and the dimensionless acceleration (A) between-2.0 to 2.0.The simulation demonstrates that the total drag is heavily correlated with A,and large deceleration even drops the drag force to a negative value.It is found that the value of virtual mass force coefficient(CV)of a spherical particle is a variable in a wide range and difficult to be correlated with A and Re.However,the total drag coefficient (CDV) is successfully correlated as a function of Re and A,and it increases as A is increased.The proposed correlation of total drag coefficient may be used for simulation of solid–liquid flow with better accuracy.

1.Introduction

It is recognized that numerical simulation of multiphase flow is of essential significance to the optimization of unit operations and scale-up of process equipment.For conducting the simulations,the correct and accurate constitutive equations for phase interaction are prerequisites to formulate the physical process under consideration.For solid–liquid systems the force interaction of different nature is to be included into the mathematical formulation.The forces relied on the relative motion of particles to a surrounding liquid are these acting in the direction of the relative motion,such as drag,virtual mass and Basset forces,and those perpendicular to the relative motion,such as lift,Magnus and Saffman forces [1].

The drag force is absolutely the most dominant interphase force.Many works have been devoted to the establishment of its empirical and theoretical expressions,and important progress has been achieved in using them for engineering design and numerical simulation [2],but such success does not happen in accurate description of other interphase interactions.This paper is particularly concerned with the virtual mass force,which is rather pronounced in dealing with gas–liquid flow.Gas bubbles with negligible mass are found not easily accelerated in liquid,because the acceleration of a bubble involves the accelerating and decelerating motion of surrounding liquid elements and extra force is demanded to fulfill this task.The extra force can be represented with the force to accelerate the virtual increase of mass of a bubble.This force also plays a significant role in formulation of non-steady motion of drops and solid particles.

van Winjgaarden [3] formulated the expression of extra virtual mass force as proportional to the product of the virtual particle mass ρVpand the acceleration a:

where CVis the virtual mass force coefficient,whose value should be designated before utilized in calculation.Theoretically CVis a function of particle shape and orientation [4],since they decide how the surrounding liquid is perturbed by an accelerating particle.It is proved theoretically that CVis 0.5 for a solid sphere in potential flow[5]and 1.12 for an elliptic spheroid.Zuber[6]pointed out that the coefficient was also the function of gas fraction.For deformable gas bubbles in a moving liquid,Cook and Harlow [7] found that using CV=0.25 in the simulation gave the best agreement with experimental data.For multi-bubble systems,van Winjgaarden[3] suggested

to account for the effect of gas fraction αg.Kuo and Walls [4]observed that when vortex was present in a liquid,CVwould be greater than 0.5,for the action of a particle would be propagated farther due to vortex motion.Lahey Jr.et al.[8] addressed the importance of including virtual mass force in the simulation of accelerating gas–liquid flows,since accurate and efficient prediction of transient two-phase flow is needed to understand nuclear reactor safety studies and the pressure-relief capabilities of chemical reactors,and many other processes.Sokolichin[9]and Lin et al.[10] chose CV=0.5 in their numerical simulations of gas–liquid loop reactors and bubble columns considering that the gas fractions in the reactors were relatively low.

In recent years,the efforts to accurate representation of virtual mas force continue.Pudasaini and co-workers extended the reliable formulation of virtual mass effects to dense liquid–solid mixtures in studying debris flows and tsunami[11,12].Computational fluid dynamics tools are also used for more fundamental aspects of virtual mass force.Simcik et al.[13] simulated several configurations of single and multiple bubbles by VOF-based direct numerical simulation,and found that the virtual mass coefficient varied widely above and below the heuristic value of 0.5.However,their target is limited to the buoyance-driven bubbles.Earlier,the similar approach was explored by Niemann and Laurien[14]on a dense regular bubble array,and summarized the simulated virtual mass coefficients as a quadratic polynomial of gas voidage.To date,the experimental and numerical studies on establishing reliable correlations of virtual mass coefficient of single bubbles and bubble clusters are reported continuously [15,16].

However,the virtual mass force of single solid particles is less studied,perhaps it is thought the virtual mass effect is not as prominent as in the case of motion of bubbles.In this paper,numerical simulation of a single spherical solid particle is conducted as a first step for establishing a reliable correlation for particle clusters.By comparison of the drag exerted on the particle in steady motion with the total drag(including virtual mass force)in acceleration,the accurate value of virtual mass force coefficient is determined.In terms of total drag,the total drag coefficient is successfully correlated so as to facilitate thereafter accurate simulation of particulate multiphase flow.

2.Hydrodynamic Formulation of a Particle in Motion

Many works addressed the steady-state motion of single particles in an unbounded liquid medium [2],and the numerical simulation of such a situation is no longer difficult [13].However,less effort was devoted to the accelerating motion of particles,especially for ones with complex shape geometry [17].

2.1.Unsteady motion in a non-inertial reference frame

For a solid particle settles or a bubble moves rectilinearly at the terminal velocity UTin an unbounded quiescent liquid phase,the fluid flow around the particle is usually described by the equation of continuity and the transient Navier–Stokes equations:

The above momentum equation suits for liquid without volume forces other than the gravity,whose effect has been incorporated into the gradient term of dynamic pressure p.For small particles in an incompressible liquid,it is usually assumed that(i)the liquid is viscous and Newtonian;(ii) the flow is laminar and isothermal.

The liquid flow can also be viewed in a reference frame fixed on the particle,thus the reference frame moves with the same velocity,and the far boundary condition becomes-UT(Fig.1(a)).In case of steady flow,the flow relative to the particle is still governed by Eqs.(3) and (4) (noting that the time derivative irrelevant).The solution of fluid flow may be done in a cylindrical coordinate system(x,r,φ)as in Fig.1(a),or in a spherical coordinate system(r,θ,φ)in Fig.1(b).In the latter system,the computational domain may be chosen as a sufficiently large sphere with its radius robeing 100 times that of the central solid particle,R.

The axisymmetrical flow may be equivalently described by a set of partial differential equations in terms of stream function ψ and vorticity ω in a sectional plane (x,y) passing through the axis of symmetry as indicated in Fig.1(a).For the axisymmetric flow around a rectilinearly moving spherical particle,using the vorticity-stream function formulation brings some convenience by avoiding the solution of the pressure field which has no differential equation for itself.By taking the curl of Eq.(4),an equation in terms of vorticity ω and velocity u is resulted,and the pressure term is eliminated.The general governing equation of the only non-zero vorticity ω component in the continuous phase in an orthogonal coordinate system (ξ,η) in Fig.2 is [18,19]

and the distortion function f is defined as the ratio of two scale factors:

The definition of stream function ψ in a left-handed reference frame (ξ,η) leads essentially to

while Eq.(10)is simply the consequence of expressing ω in terms of the stream function ψ.Eqs.(9) and (10) are the mathematical formulation to be resolved for the external axisymmetric flow of liquid.The correspondence of the physical domain in (x,y)coordinates with the computational domain in (ξ,η) coordinates is sketched in Fig.2,where letters A through F indicate how the boundary points in two domains correspond to each other.

When transformed non-dimensional in terms of the relevant non-dimensional physical variables Ω and Ψ defined by

the governing equations set become [19,20]

with the non-dimensional time defined as

Fig.1.Sketch of a rising spherical particle.

and the Reynolds number Re=2ρRUT/μ,X=x/R,Y=Y/R,with R as the particle radius.The values of the distortion function remain unchanged after the non-dimensionalization:

For this system,the Reynolds number Re is the only nondimensional controlling parameter.

In the computational domain in Fig.2,the non-dimensional velocity components normalized with UTare related to the stream function by

When applied to the flow of liquid relative to the center particle,the above mathematical formulation remains valid in this inertial reference frame.Thus,the problem becomes a fixed particle stays stationary facing the coming liquid flow.The boundary conditions at the particle surface and at the axis of symmetry are routinely used:

At the far field boundary (ξ=0),

At the particle surface (ξ ＝1),

In case of a particle in acceleration,it will be subjected to larger resistance than that at steady motion.Not only extra force for accelerating the particle itself,liquid in a larger region would be disturbed by the accelerating particle.In the inertial reference frame fixed with the earth,the governing equations remain to be Eqs.(3) and (4).

If the liquid flow field is viewed in a reference frame fixed on the accelerating particle,namely the non-inertial reference frame moves with the same linear velocity UTand acceleration a of the center particle (Fig.1),this acceleration needs to be incorporated into the momentum balance.Based on the linear momentum equation for a non-inertial control volume[21],the differential governing Navier-Stokes equation of the relative velocity u for the present case of rectilinear acceleration becomes

with the associated acceleration a included.The far domain boundary condition remains to be

Fig.2.The grid for numerical solution,the physical reference frame and the computational coordinate system.

but the particle velocity UTis now a time dependent variable.The associated acceleration a may vary with the time,but remains a constant vector in the whole space.Since the density of liquid is constant,the continuity Eq.(3) remains valid.When nondimensionalized with the aid of Eq.(11),the governing equations remain in the same forms as Eqs.(12) and (13),and the boundary conditions Eqs.(18)–(22) work as well.When Eq.(23) is nondimensionalized,another non-dimensional number A=is introduced (see Appendix B),in addition to the already introduced Re.

2.2.Stress tensor for particle in unsteady motion

When a spherical particle moves under the gravity in an unbounded liquid at steady state,its drag coefficient CDis defined by

where FDis the drag force exerted on the particle by a surrounding liquid.It is necessary to calculate the drag on the particle to get the value of CD.

The total force FTexerted by the flow onto an accelerating particle includes the virtual mass force,which remains mathematically the same form of a surface integral of the liquid stress on the particle:

where T is the stress tensor and n=eξis the unit vector normal to the particle surface.Since the particle moves rectilinearly upwards in Fig.1,the total drag is exactly in the x-direction:

Though only the stress on the particle surface (ξ=1) is needed in Eq.(27),the correct expression of relevant stress components must be evaluated in a three-dimensional scenario.The expression of T of axisymmetric liquid flow around a closed surface in the present orthogonal coordinate system (ξ,η,φ) as that in Fig.2 is

where τξηmeans the stress acted on the plane normal to the ξ axis but in the η direction,and so on.

With the surface arc segment dS ＝（2πy）hηdη and the following relations:

In the appendices,it is shown that for a solid sphere with noslip surface(ξ=1),the non-dimensional normal stress component is

In the present case of a solid particle,not only the normalized pressure P contributes to the normal stress tensor component τξξ,but also the particle acceleration adds to the hydrodynamic pressure.The non-dimensional shear stress component is

2.3.Separation of virtual mass force from the total force of traction

For an accelerating particle,its traction from a surrounding liquid is not just equal to the drag force that a particle experiences when moving at steady state at the same linear velocity.The value of Fxis greater than FDand the excess is spent on accelerating the solid particle (or a bubble with zero mass),or in other words,to moving the liquid in the immediate neighborhood of the particle(called the virtual mass force FV).The virtual mass force coefficient CVis defined by Eq.(1).If we presume that the drag force at steady motion and the virtual mass force are linearly additive:

then the above dimensional expression for a spherical particle may be made non-dimensional to

the left hand term gives the total drag coefficient CDV,which is similarly defined as Eq.(25).The left hand term of Eq.(35),when stress non-dimensionalized withand length by R in that numerator,reduces to a simple non-dimensional form used in the latter FORTRAN program:

The virtual mass coefficient CVcan be evaluated from Eq.(35)as

Thus,it seems feasible to separate the virtual mass force from the total drag exerted onto an accelerating particle.Please note that when the flow is solved in non-dimensional variables,the total drag coefficient is calculated by

3.Numerical Method

3.1.Numerical procedure

Numerical simulation of axisymmetrical liquid flow outside an accelerating solid particle in a non-inertial reference frame proceeds by solving Eqs.(12) and (13) in terms of non-dimensional vorticity Ω and stream function Ψ with pertinent boundary conditions.The computational domain is the vertical section through the axis of symmetry:the half annulus with the inner radius R and the outer radius ro=100R as sketched in Fig.2.With the help of the orthogonal boundary-fitted grid and the coordinate transform,the physical (X,Y) domain outside the particle is transformed into a square computational (ξ,η) domain [18,22].In the present case of axisymmetrical liquid flow,the vorticity vector has one nonzero component only.

The governing equations,Eqs.(12) and (13),are valid for both inertial and non-inertial coordinate systems,i.e.,for a solid sphere in steady or accelerating motion.The expressions for stress tensor component τξξand τξη,Eqs.(31)and(32),also apply to the case of A=0.

The simulation procedure is outlined as below:

(1) Specify the orthogonal grid in the physical (X,Y) domain[22,23].

(2) The liquid flow of steady state motion of a particle with designated Re is solved at first,then the pressure and tensor components at the surface are extracted,followed by the total drag.

(3) The particle is set to accelerate with designated constant acceleration A,and solve the time-dependent flow in the real time domain with small time step Δθ.

(4) Calculate Fxthe total traction on the particle by Eq.(38),which is the net force in the x-direction exerted by liquid on the particle,or CDV.The drag coefficient CDfrom literature is used to obtain the virtual mass coefficient CVfrom CDVby Eq.(37).(5) At the beginning of the(n+1)-th time step,update the flow parameters since the particle moves now at UTn+1faster than at UTn.The update involved with is

For the sake of better numerical convergence,a constant small time step Δθ is chosen.

(6) Go back to step(3)for further solution of the next time step.

The numerical program was written in the FORTRAN language as done in the previous works[22–26].The computational domain of a square 1 × 1 in the transformed ξ–η plane is discretized into(NI-1)×(NJ-1) cells of the same size of Δξ × Δη.However,in the physical domain shaped like a half annulus (in Fig.2),the discretization is done in a polar coordinate system:the NJ grid lines along the radial direction are uniformly spaced,but the NI peripheral coordinate lines are distributed with geometrically increasing spacing,with the smallest spacing at the sphere surface,which is Δr=R/(NI-1).

The governing equations are discretized by the control volume approach as described by Patankar [27]:Eq.(13) is discretized by the central difference scheme,Eq.(12) having convectional terms by the second-order QUICK scheme as detailed by Mao and Chen[22],and the temporal derivative term by the first-order explicit differencing scheme.

3.2.Grid independence and validation tests

In choosing the proper spatial mesh density and size of time step to assure the numerical simulation not being biased from the true physics,the grid independence test was conducted as shown in Table 1.As an example,the steady flow around a sphere with Re=10 is simulated with grids of different fineness and in a large but finite external flow domain.The simulated drag coefficient is compared with the simple and popularly used Schiller and Naumann correlation [28]:

When Re=10,the value of CDis 4.1510 from the present simulation.Careful examination of the results listed in Table 1,it may be concluded that a 321 × 81 grid with the outer radius roof the computational domain equal to 100 time the solid sphere radius R=1 is good enough,which gives CD=4.3078,only 3.8% above the empirical correlation.It seems that the smaller radial cell size is more important to the accurate evaluation of the shear friction included in the total drag to a particle in motion.The computational domain with ro=100R is necessary,particularly for particles at low Reynolds numbers.For particles with larger Re,the region of flow field disturbed by a moving particle shrinks gradually,especially in the lateral direction,more confidence in simulation being assured.All later simulations are conducted with a non-uniform 321 × 81 grid for the external domain with ro=100R.

Table 1 Grid independence test in term of simulated drag coefficient (ro=100R,Re=10)

For the simulation of transient liquid flow around accelerating spheres,the dimensionless Δθ below 0.01 was enough to achieve the independence of solution,as suggested by our previous experience [20,24].Based on the computational procedure above,the effect of time step size on the numerical results of total traction was investigated.For a particle of Re=10 at A=5.0,Δθ ≤2.5 × 1 0-5guarantees the desired convergence and accuracy,and for a particle of Re=10 at A=0.2,Δθ ≤1 × 10-3is sufficient.This requirement on Δθ is in compliance with the general law on the Hirt stability criterion:

for a time-dependent partial differential equation such as

The present numerical simulation was validated against the Schiller–Naumann correlation,Eq.(40),for solid spheres in steady motion.As demonstrated in Fig.3,the agreement of the simulated drag coefficient with the correlation is satisfactory in Re ranging from 0.5 to 200,with an average absolute relative error of 2%.

4.Results and Discussion

In this work,the steady and accelerating motion of solid spheres in the range of Re from 10 to 200 with the nondimensional acceleration number A between-2.0 and 2.0 are simulated.The flow field structures,and the predicted total drag coefficient and virtual force coefficient are presented and analyzed.Moreover,an empirical correlation of total drag coefficient is proposed for its application in numerical simulation of particle motion in complicated situations.

4.1.Non-steady flow structure

Intuitively it is anticipated that at the same value of Re the fluid flow around a spherical particle in acceleration would be different from that around a steadily moving sphere due to the inertia of fluid.It is found that the difference is quite obvious for a sphere between steady motion and transient motions at the same Re=50.Fig.4 illustrates the flow structures of an accelerating or decelerating sphere temporarily at Re=50 are quite different from that of a steadily moving sphere.At the nose(front stagnant point)the compression of streamline spacing by acceleration (A=0.5) or the relaxation by deceleration (A=-0.5) is marginally observed,but the length of the circulating wake behind a steady sphere is obviously larger than that of a sphere in acceleration and muchshorter than that in deceleration.More difference in the velocity components and pressure across the flow field can also be observed.

Fig.3.Comparison of simulated drag coefficient of a steady moving sphere with the empirical correlation (321 × 81 grid).

4.2.Behavior of total drag coefficient

The total drag of a sphere at steady motion at Re=50 is extracted from the simulated flow field,and expressed as the total drag coefficient CDVin Eq.(38).In steady motion it is essentially the conventional drag coefficient CD,having been shown in good accordance with the well-known Schiller–Naumann correlation in Fig.3.When the sphere moves at A=0.2 from Re=10 to 100,the total drag is enhanced significantly to give a larger CDVin Fig.5,while when the sphere is decelerated with A=-0.2 from Re=100 to 10,the total drag is reduced by a large amount.

The virtual mass force coefficient CVevaluated from Eq.(37) is presented in Fig.6,showing that CVdepends on both Re and the algebraic sign of A.Also its value is largely biased from 0.5 which is the theoretical result for a sphere in an inviscid fluid.The trend of change of CVseems irregular.Part of the reason might be due to the partly empirical nature of CDfrom Eq.(40),but it is believed that the main reason is relied on other complex interaction of acceleration with the flow field.

The effects of acceleration on CDVand CVare presented in Fig.7.It is observed that CDVis increased by acceleration and decreased by deceleration,both in a monotonic way.Acceleration and deceleration seem to generate their unique patterns in changing CV.The more regular and roughly parallelism of lines in Fig.7(a) suggest that the job of establishing an empirical correlation for CDVis relatively easy,while the correlation of CDVdata in Fig.7(b) would be much difficult in view that the particle behaviors are quite different in the cases of acceleration and deceleration.

Fig.4.Comparison of the streamline maps of particle with Re=50(on 321×81 grid):(a)at steady state;(b)in accelerating motion with constant A=0.5 from Re=10;(c)decelerated with constant A=-0.5 starting from Re=100.Stream function values plotted:0,0.0025,0.005,0.01,0.02,0.04,...&-1 × 10-3,-2 × 10-3,-4 × 10-3,...

Fig.5.Total drag coefficient of a sphere in steady and non-steady motion(with A at 0.2 and -0.2).

Fig.6.Virtual mass force coefficient of a sphere in steady and non-steady motion(with A at 0.2 and -0.2).

4.3.Existing correlations on accelerating particles

Fig.7(a) demonstrates that the virtual mass force is rather important when a particle is in non-steady motion.However,it is still ignored in some current numerical simulations of particulate-fluid flow,because there is a lack of reliable and applicable correlations of virtual mass coefficient.

Zhang et al.[29] proposed a correlation for the total drag coefficient of bubbles in accelerating and steady motion in highly viscous liquids in the form of

where A is the acceleration number (A ＝) and Ar is the Archimedes number(Ar ＝).The bubble acceleration effect due to the added mass force was separated from Eq.(40) as

This strategy also worked for non-Newtonian liquids [30].Eq.(43) for spherical bubbles was used together with the literature drag force correlations to simulate rising bubbles released from a nozzle in a quiescent liquid,and the simulated rising process was found in good agreement with experimental measurements for a wide range of liquid properties and bubble size [29].

So far,no similar reports on solid particles and liquid drops are available.It is also desired to get similar correlations for solid spheres and liquid drops for reliable and efficient simulation of two-phase flow with dispersed particles.

Fig.7.The effects of acceleration on total drag coefficient and virtual mass force coefficient.

4.4.New correlation of total drag in non-steady motion

It is interesting to see if a correlation for the total drag coefficient and virtual mass force coefficient can be established,by incorporating the acceleration number A into the Schiller-Naumann correlation:

Unfortunately the trial fails when correlating in this form.

After many trial-and-error efforts,the following functional form is found a good candidate:

where the symbol ± stands for +1 for acceleration and -1 for deceleration.

The data in Fig.7(a) were fitted to Eq.(46) to decide the undetermined parameters.A complex algorithm in FORTRAN for optimization was used for data-fitting.The final result is

where the average absolute error is 0.060.Fig.8 shows the satisfactory coincidence between the simulated and correlated total drags.

It is desirable to know the mechanisms by which the acceleration influences the total drag onto a particle.Eq.(47) suggests superficially that the value and algebraic sign of a contribute largely to the change of CDV,because firstly,non-dimensional A appears in the expression of surface pressure,Eq.(31),and the term AX will contribute to the integral of Eq.(38);secondly,the pressure distribution having been altered by A will further change the velocity vectors due to the interaction embodied in the Navier–Stokes equation,but this mechanism is rather implicit.It is interesting to note that the term of contribution of A,±（δAε）,looks as if irrelevant of the particle Reynolds number,and Re and A seem to act separately.This observation awaits a thorough interpretation.

Fig.8.The contrast of the numerical simulation with the proposed correlation on the total drag coefficient.(Re:10 to 200,A:-1.0 to 1.0).

5.Conclusions

Though the virtual mass force is indispensable in the momentum balance involved with non-steady motion of particles in multiphase systems,the numerical simulation on this topic is still very scarce.In this work the accelerating motion of a solid sphere in an orthogonal curvilinear coordinate system is mathematically formulated with the vorticity-stream function formulation and solved numerically with the control volume scheme.The simulation provides numerical data on total drag coefficient and virtual mass force coefficient.It is found that the acceleration increases the total drag and the deceleration decreases it significantly at almost the same extent.The total drag coefficient behaves much regularly,while the virtual mass force coefficient shows a complex trend of variation.Though the effort of correlating the virtual mass coefficient fails,the total drag summing up the steady drag force and virtual mass force is successfully correlated as a function of the particle Reynolds number Re and the acceleration number A.The proposed correlation represents well the results of total drag,and is hopeful in making contribution to better accuracy in numerical simulation of particulate multiphase flow.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was financially supported by the National Key Research and Development Program (2020YFA0906804),the National Natural Science Foundation of China (22035007,91934301),External Cooperation Program of BIC,Chinese Academy of Sciences (122111KYSB20190032),and Chemistry and Chemical Engineering Guangdong Laboratory,Shantou (No.1922006).Dedicated to deep-in-heart memory of prof.CY Chen,our caring teacher and guide.

Nomenclature

A dimensionless acceleration

a acceleration,m·s-2

CDVtotal drag coefficient

CVvirtual mass force coefficient

d particle diameter,m

Fxprojection of total drag in the x-direction

f(ξ,η) distortion function

g gravity acceleration (=9.812 m·s-2)

Hξ,Hηnon-dimensional scaling factor

hξ,hηscaling factor,m

p pressure,Pa

R particle radius,m

Re Reynolds number (=UTdρ/μ)

roouter domain radius,m

U dimensionless velocity (=u/UT)

UTterminal velocity,m·s-1

u velocity component,m·s-1

X,Y dimensionless coordinate in physical plane by R

x,y coordinate in physical plane,m

θ dimensionless time,(=tUT/R)

μ viscosity,Pa·s

ν kinematic viscosity,m2·s-1

ξ,η coordinate in computational plane,0 ≤ξ,η ≤1

ρ liquid density,kg·m-3

τ stress tensor component

φ azimuthal angle

Ψ dimensionless stream function (=ψ/(R2UT))

ψ stream function,m3·s-1

Ω dimensionless vorticity (=ωR/UT)

ω vorticity,s-1

Subscripts

DV total drag

V virtual mass

Appendix A

Stress tensor components for a solid particle in steady motion

There are only two independent components inT,τξη=τηξand τξξbeing relevant to the total drag force in axisymmetric flow[31]:

since τηηand τφφdo not contribute to the force on the particle surface.

Eq.(A1) may be simplified further.Since liquid is incompressible,

Because the flow is axisymmetric,the third term in the bracket is zero,and with hφ=y,the above equation reads

One the solid particle surface uξ=0,the first term disappears,leading to

Substituting this expression into Eq.(A1) results in

For a solid particle,both velocity components uξand uηat the surface are zero,thus leading to

With uξ=uη=0 on the surface,Eq.(A2) is reduced to

Thus the total traction force at the surface in the x-direction is the projection of FTin the x-direction,i.e.,

or in terms of tensor components as

When transformed to non-dimensional in terms of relevant non-dimensional vorticity Ω and stream function Ψ as defined in Eq.(11),and the stress components,Eqs.(A7) and (A8),together with pressure p normalized withthe relevant dimensional stress components may be converted to Ref.[20]:

It is necessary to get the expression for the static pressure p in Eq.(A9).This may start from the time-dependent momentum equation in the inertial reference frame:

As being non-dimensionalized by divided withit reads

The non-dimensionalizing factors are UTfor u,for p,and R/UTfor t(θ=tUT/R).In axisymmetric flow the only non-zero components of U and Ω are

It is noted that the vorticity vector Ω has only one component in the φ–direction.In the above,eξ，eηand eφare unit vectors in respective axis.

Then,the expanded expressions for some terms in Eq.(A12)are

Substitute the above into Eq.(A12) and then dot multiplied with eηresults in

Because only the pressure at the surface(ξ=1)is needed in Eq.(30),for a solid sphere in steady motion with=0,the integration of Eq.(A19)along the surface from η=0 gives the pressure on the particle surface over there with Uξ=Uη=0,and we have

There is an undetermined constant C in Eq.(A20),but on integration over the whole surface of a particle to get the total traction in the x-direction by Eq.(30),it makes no contribution to Fxat all,and thus can be omitted from Eq.(A20).

Eqs.(A21) and (A10) are inserted into Eq.(38) to get the total drag on a solid particle.

Appendix B

Stress tensor components for accelerating solid particle

The governing equations remain to be Eqs.(3) and (4) when a particle moves with acceleration a.Now the flow field becomes time-dependent.If the liquid flow field is viewed in a reference frame fixed on the accelerating particle,namely the non-inertial reference frame moves with the same linear velocity UTand acceleration a as the center particle(Fig.1),the Navier-Stokes equation of the relative velocity u becomes

However,the acceleration a is involved and the flow field becomes time-dependent.The counterpart of Eq.(A12) becomes

where the non-dimensional acceleration number is A=.Eq.(A18) is now

Eq.(A19) now reads

With the no-slip condition at the surface (ξ=1),we have Eq.(31) in Section 2:

Similar,the integration constant becomes irrelevant when integrated over a closed surface.

The tangential stress component τξηremains to be Eq.(A10),i.e.,Eq.(32).