Yuan-Sheng Tang

1　Introduction

In today’s manufacture,the polymer filling process possesses a crucial stage in the polymer molding,for example the form ing process of food or plastic products.It is well known that it is very important to fill these containers as rapidly as possible without spillage or deterioration in the product[Tomé et al.(1999);Hwang and Kwon(2002)].Some instability phenomena occur when the polymer mould design or parameters of fluid is inappropriate,which can lead to undesirable shape.Therefore,it is necessary to develop more effective control of the mould injection performance,especially for 3D mould filling.In addition,the capture of instability phenomena are also considered,and which are usually involving the study of complex free surface flows and entangling of interfaces,such as the splashing,spluttering,sloshing or jet buckling/coiling.Recently,the numerical simulation offers a powerful and effective means to track the 3D complex free surface and study the 3D filling process due to the difficulties of experimental way for 3D polymer melt molding.In order to better demonstrate the shear or nonlinear behavior of polymer melt in the filling process,the branched polymer melts based on the three-dimensional eXtended Pom-Pom(XPP)model[Vebeeten,Peters and Baaijens(2001);Oishi et al.(2012)]is mainly studied in this paper,which can degenerate to the Oldroyd-B fluid.

Many grid-based numerical methods have been developed and applied to treat the complex viscoelastic free surface and to simulate the filling process with the condition of low pressure or low velocity,such as finite difference method(FDM)and finite volume method(FVM)coupled with marker and cell(MAC)[Oishi et al.(2012);Tomé et al.(2002)],volume of fluid(VOF)and level set[Li et al.(2011)]methods.However,the mesh-based methods based on the Eulerian description are usually suffer from some difficulties for dealing with the extremely complex free surface behavior found in high pressure or high velocity filling process,such as the mesh generation for 3D mould filling is very expensive and so on[Hwang and Kwon(2002);Tomé et al.(2008)].For the reasons,the grid-free method is first considered and explored to study the complex 3D polymer filling process.

In recent decades,the smoothed particle hydrodynamics(SPH)method has been widely used to simulate large deformation problems in mechanics field[Li and Liu(2002);Zhou and Ge(2014)],which is a special case of the MLPG approach[Sladek et al.(2013)]and is regarded as a pure mesh-free particle method based on Lagrangian description.The main merits of SPH method over than the gridbased methods lie in:(a)It is easy to program for complex problems especially for 3D case without mesh reconstruction;(b)Complex free surfaces are modeled easily and naturally without the need of explicit surface tracking technique.Since then,The SPH method has been applied and extended to many fluid areas such as 3D viscous flows[Ferrari et al.(2009);Ferrari et al.(2009);Cleary et al.(2007)],2D multi-phase flows[Colagrossi and Landrini(2003);Hu and Adams(2006)],2D incompressible flows[Skillen et al.(2013);Zainali et al.(2013],and 2D viscoelastic fluid flows[Ellero and Tanner(2005);Fang et al.(2006);Ren et al.(2011)].In addition,the SPH method combined with the ALE formulation is well used to handle the Fluid Structure Interaction problems,which can be seen in recent work[Messahel and Souli(2013)].However,the numerical simulations of complex 3D polymer filling process based on viscoelastic constitutive model using the SPH method are rare.

On the other hand,the traditional SPH(TSPH)method is usually suffer from two major drawbacks,that is,the low accuracy and numerical/tensile instability.Hence,some corrective SPH methods based on the concept of Taylor series to restore the consistency[Chen and Beraun(2000);Liu,Xie and Liu(2005);Zhang and Batra(2004);Liu and Liu(2006);Zhang and Batra(2007);Batra and Zhang(2008)]of the kernel and gradient particle approximations of TSPH method and remedy its accuracy and stability,for example the corrected smoothed particle hydrodynamics method(CSPM)[Chen and Beraun(2000)], finite particle method(FPM)[Liu,Xie and Liu(2005)],modified SPH(MSPH)method and the symmetric SPH(SSPH)method[Zhang and Batra(2004);Liu and Liu(2006);Zhang and Batra(2007);Batra and Zhang(2008)].Unfortunately,these improved SPH methods have not extensively applied to viscous or viscoelastic free surface flows due to some disadvantages of themselves:the singularity of formed local matrix may cause serious instability;it is very complicated to extendedly apply them for complex free surface flows problems,especially for 3D problems which can be seen from their construction process.In order to find a proper improved SPH method by adopting the concept of Taylor series in the simulation of viscous fluid or free surface flows,the corrected kernel gradient scheme by adopting the partial concept of CSPM method has been introduced to the discretization schemes of TSPH and extensively applied to viscous fluid flows problems[Bonet and Lok(1999);Ren et al.(2011);Stranex and Wheaton(2011);Fetehi and Manzari(2012);Jiang,Ouyang and Ren(2012)]in recent years.However,it has not been still extendedly applied to simulate 3D viscoelastic free surface flows in our know ledge.Moreover,some other techniques are still need to remedy the tensile instability or decrease pressure oscillations even if the treatment of complex boundaries and computational efficiency of 3D particle method are also needed be enhanced.

For all the above reasons,it is necessary to seek an appropriate corrected SPH method for solving the 3D viscoelastic fluid flows problems and simulating the viscoelastic free surface flows.In this work,we mainly focus on a balanced improved SPH method that can be conveniently extended to simulate the complex filling pro-

cess based on polymer melts,and the proposed method compromises the accuracy,the stability and computational efficiency between the traditional SPH and other improved SPH methods(e.g.FPM,MSPH and SSPH).Then a 2D corrected SPH scheme combined with a density diffusive term is considered and extended to simulate the 3D viscoelastic free surface flows,which has been introduced in our previous work[Bonet and Lok(1999)]and mainly applied to 2D filling process of viscous fluid.The proposed 3D particle scheme is motivated by a coupled concept that the extended 2D kernel-gradient-corrected SPH(KGC-SPH)method without involving kernel derivative in[Ren et al.(2011);Ren et al.(2011)]is used in inner fluid region and the TSPH method is used near the boundary domain,name as the corrected 3D SPH(CSPH-3D)method.It is worth noting that the derived discretization scheme of momentum equation using the KGC-SPH method is different from the momentum discretization scheme in[Ren et al.(2011);Fetehi and Manzari(2012)].In addition,a new boundary treatment is also presented to treat complex multiple solid walls,which is effective and easier to perform than that in[Xu et al.(2012)].We can also know that a huge amount of calculation have to be handled when the proposed CSPH-3D method for solving 3D viscoelastic fluid problems due to the search of neighbor particle.For enhancing the computational efficiency,the MPI parallelization means with a dynamic cells neighbor particle search method is adopted and performed in an IBM HPC Platform with c++code.In this paper,the 3D polymer filling process based on the XPP model is mainly simulated using the proposed CSPH-3D method with density diffusive term and multiple boundary conditions and it is performed in an IBM HPC Platform with MPI parallelization means,in which the filling process with two-inlets in a rectangle container or in a ring cavity.Some interesting instability phenomena(e.g.jet coiling)are also observed.Particularly,the influences of the macroscopic rheological parameters on the complex filling process is numerically predicted and discussed.This work is organized as follows:The 3D governing equations for the XPP fluid are introduced in section 2;Section 3 describes the CSPH-3D discretization scheme of the Navier-Stokes equations based on the XPP model,including boundary treatment,density diffusive term and artificial viscosity;In section 4,the validity and ability of proposed particle scheme is tested by several benchmarks involving the 3D filling process of XPP fluid with a single inlet,and the merits of proposed CSPH-3D method combined with the density diffusive term and boundary treatment are also demonstrated;The 3D filling process with two type containers w ill be numerically predicted in section 5 and section 6,respectively;Some complex free surface phenomena have been also illustrated and discussed with different rheological parameters in sections 5 and 6.Conclusions and remarks are reported in section 7.

2　Three-dimensional governing equations for the XPP fluid

2.1　Governing equations

In a 3D Lagrangian frame,the flow of concentrated polymer solutions and melts is governed by the conservation of mass and momentum equations,together with a constitutive equation.The fluid flow is usually described by the follow ing governing equations:

whereis the velocity vector,ρis the fluid density,tis timeis the gravitational acceleration and? is the material derivative??The total Cauchy stress tensorσin Eq.(2)is usually decomposed into the ordinary isotropic pressurep,viscousand polymeric extra stress tensor

whereIrefers to identity matrix,Dis the rate of deformation tensor,which is given by

2.2　Constitutive model for the XPP

In this work,the extended pom-pom model(XPP)in multi-mode form[Vebeeten et al.(2001);Oishi et al.(2012)]is considered to study the influence of shear or viscoelastic behavior on free surface in polymer filling process.The major features of XPP model are:the dependence of melt rheology upon the polymer molecular structure;the spectrum of relaxation time to be taken into account leads to characters of orientation and stretch.The constitutive model for the XPP fluid is described as:

where the functionis given by

isthe linear relaxation modulus,Irefers to the trace of a tensor.andλ0sare the orientation and backbone stretch relaxation time respectively)shows that small values ofλrcorrespond to highly entangled backbones[Vebeeten et al.(2001)],and the anisotropy parameterThe parametervin the exponential term in Eq.(6)is incorporated into the stretch relaxation time to remove the discontinuity from the gradient of the extensional viscosity.Its value is found to be inversely proportional to the number of arms?The backbone stretchλis related to the viscoelastic stress tensor

Here,the following parameters are introduced,namely the total viscosity?

In addition,the XPP model can degenerate the Oldroyd-B model whenα0=0 andin Eq.(5),and the UCM model is further obtained ifβ0=0.It is worth noting that two type constitutive models are considered for easy of comparison in this paper,i.e.the XPP and Oldroyd-B models.

2.3　Equation of state

In the application of TSPH method,the incompressible flows were usually treated as slightly compressible flows by adopting a suitable equation of state(see Monaghan[Monaghan(1994)]and Morris et al.[Morris,Fox and Zhu(1997)]).Here,the incompressible flows are also treated as weakly compressible flows using the following equation of state[Ellero and Tanner(2005)]

wherecis the speed of sound andρ0is reference density.An artificial,lower sound speed is usually used to avoid the instability and extremely small time steps.To keep the density variation of fluid less than 1%of the reference density,the Mach number?whereVis a typical reference velocity)[Zhou and Ge(2014)]must be smaller than 0.1.

3　The Corrected three-dimensional SPH scheme

In this section,a corrected 3D SPH(CSPH-3D)scheme is proposed to solve the 3D viscoelastic fluid problems based on the XPP model,that is mainly derived a coupled approach between the KGC-SPH-3D method without kernel derivative and TSPH-3D method.The CSPH-3D scheme has higher accuracy and better stability than TSPH-3D,in which the KGC-SPH-3D scheme is achieved by introducing the extended corrected scheme of first-order kernel gradient in the discretization schemes of governing equations and the construction process is based on the alliance of particles interaction.It is worth noting that the achieved KGC-SPH scheme is different from it in[Ren et al.(2011);Fetehi and Manzari(2012)].Moreover,the 2D density diffusive term(see[Ren et al.(2011);Antuono,Colagrossi,Marrone and Molteni(2010)])is extended and introduced in the CSPH-3D method for restraining the pressure oscillations.The CSPH-3D scheme for the XPP fluid flow is outlined as follows.

3.1　TSPH-3D scheme

In the constitution of TSPH-3D method,the integral interpolation theory is used for a kernel function.The fluid region ? is discretized into a finite number of particles in 3D system,where all the relevant physical quantities are approximated in terms of the integral representation over neighboring particles.Each particle carries a massm,velocityu,and other physical quantities depending on the problem.

For any functionfand its first-order derivative defined at the positionr=(x,y,z)are usually expressed by the following integral(see[Li and Liu(2002);Zhou and Ge(2014)]).

whereWdenotes the kernel function(or smoothing function)andhrepresents the smoothing length de fining the influence area ofW.TheWis usually needed to meet three properties(see[Zhou and Ge(2014)],that is the Dirac function propertythe positive propertyand the compact propertyoverwhen(whereκis a constant).

The integrating principle by parts and the divergence theorem are applied to the Eqs.(10)and(11),we can get the particle discretization scheme of TSPH for a functionf(r)and its first derivative at the positionwhich can be written in the condensed forms

The smoothing functionWis related not only with the accuracy but also with the efficiency and stability of the resulting algorithm.According to reference[Wendland(1995)],the Wendland function can produce more accurate results than the common splines kernel functions(e.g.the cubic or quintic splines function).The 3D quintic Wendland kernel function is used in this work,which is

whereuβis theβ–thcomponent of the fluid velocity,σαβis the(α,β)–thcomponent of the total stress tensor andxβis the component of 3D spatial coordinate.

The particle discretization form of total stress equation(3)can be de fined as

The discretization scheme of constitutive equation(5)-(6)for the XPP model is

Where

The above mentioned formula(16)-(17)of TSPH-3D method give expression to particle interaction of neighbor particles.

3.2　CSPH-3D scheme

The 2-D corrected symmetric kernel gradient scheme[Bonet and Lok(1999);Ren et al.(2011)]and the Taylor expansion concept[Chen and Beraun(2000);Liu,Xie and Liu(2005);Zhang and Batra(2004);Liu and Liu(2006);Zhang and Batra(2007);Batra and Zhang(2008)]are extendedly applied to remedy the first-order derivative of TSPH-3D scheme i.e.the Eqs.(15)-(17),respectively.The obtained corrected scheme can restore the particle approximations consistency and improve the numerical accuracy and stability of TSPH-3D method.The extended KGCSPH-3D schemes without involving kernel derivative for the viscoelastic fluid flow based on the XPP model,which are

where

Where

It should be noted that the corrected scheme(22)mentioned above is different from the relevant corrected scheme in[Ren et al.(2011);Fetehi and Manzari(2012)],which can accurately tracing the phenomenon of jet coiling in the simulation of polymer filling process in a rectangle container.Moreover,the corrected scheme does not include kernel derivative(see Eq.(24))which has some merits of MSPH and SSPH method(see[Xu et al.(2012);Monaghan(1994)]).

The CSPH-3D method is motivated by the coupled concept between TSPH-3D method and the KGC-SPH-3D method.Its idea is using the KGC-SPH-3D in the interior fluid field and using the TSPH-3D near the free surface(or boundary particles),and the boundary particles or free surface particles may be identified by particle densities[Zhou and Ge(2014);Monaghan(1994)]].

The CSPH-3D schemes for the XPP fluid flow can be obtained by using Eqs.(15)-(18)for the boundary fluid particles and using Eqs.(18),(21)–(23)for the interior fluid particles,which are

where

At the end of each time-step,the position of each particle is updated using

3.3　Numerical diffusive and dissipative terms

It is well known that the pressure oscillations easily occur and happen to numerical noise when the density values are update by using the continuity discretization equation(25).In order to avoid this problem,a 2D density diffusive term is extended and introduced in the continuity equation which can effectively eliminate a large part of the numerical noise(see[Ren et al.(2011);Antuono,Colagrossi,Marrone and Molteni(2010);Marrone et al.(2011)]).On the other hand,a dissipative term(name as the “artificial viscosity”,see[Zhou and Ge(2014);Monaghan(1994)])is also added to the discrete momentum equation(26)for the purpose of increasing the stability of numerical simulations.As a result,the discretization schemes(25)and(26)can be written as

where theφijand Πijare chosen to be

A ll the coefficients mentioned above are positive and have to be properly tuned.The coefficients?,χcontrol the order of the magnitude of the diffusive term and viscous term,respectively,which are usually chosen as 0＜?＜0.1,0＜χ＜0.2.Theη=0.1hterm is included to prevent numerical divergence when two particles get too close to each other.TheαIandβIare usually chosen approximately equal to 1.

3.4　Time integration scheme

In order to solve the system of ordinary differential equations(30),(31),(18)and(29),the predictor-corrector scheme is used in this paper which possesses secondorder accuracy and better stability.The predictor step consists of an Eulerian explicit evaluation of all quantities for each particle

The time step and space step must satisfy the well-known Courant-Friedrichs-Lewy(CFL)condition for ensuring the numerical stability.According to[Fang et al.(2006)],we may choose the following stability condition

whereis the kinematic viscosity.

3.5　Boundary treatment

In SPH simulations,the treatment of physical boundary including the free surface and rigid wall is very important.The boundary can be stationary or in motion.Fortunately,the boundary condition of free surface is satisfied naturally by the SPH or improved SPH methods according to references[Zhou and Ge(2014);Fang et al.(2006);Monaghan(1994)].Therefore,two types of boundary conditions needed to be properly treated in the simulations of polymer filling process which are the rigid wall boundary and the in flow boundary except for the free surface.

3.5.1Rigid wall boundary

In this paper,the multi-wall boundaries are involved due to the complex filling process in a closed three-dimensional container(see Fig.2).The multi-wall boundaries mainly consist of the straight and curve boundaries which can be seen in Fig.1.Several methods for treating rigid wall boundary conditions have been presented in previous work[Zhou and Ge(2014);Monaghan(1994)].There are mainly two methods i.e.the artificial repulsive force method and the image particles method.Considering the advantages and disadvantages(see[Zhou and Ge(2014);Fang et al.(2006)])of mentioned above methods and the complexity of closed three dimensional filling container,a virtual particles refinement method is presented to implement the multi-wall boundary conditions to prevent fluid particles from penetrating rigid walls in this work,and its convenience and validity are demonstrated in the following simulations.

Figure 1:The pro file graphics of different rigid particles distribution.

Figure 2:Initial sketches of filling process:(a) filling with two-inlets in rectangle container;(b) filling of die casting in ring container.

The proposed solid wall boundary treatment is feasible to prevent fluid particles from penetrating rigid wall(see section 4)and possesses some advantages overt than that of proposed wall boundary treatment in[Fang et al.(2006);Xu et al.(2012)].

3.5.2In flow boundary

The enforcement of in flow boundary condition is also taken into account in the complex filling process,which lies in the initial distribution of fluid particles before entering the container and the treatment of inlet velocity near the in flow boundary(see Fig.2).The numerical accuracy,stability and the capture of some phenomena in simulations of filling process are all affected by the in flow boundary treatment.For the purpose of obtaining a desired numerical simulation,the in flow boundary treatment can be enforce as with different containers(see Fig.2):(a)At the first time-step of numerical simulation,many reserved fluid particles are uniformed distribution and away from the nozzle entry,and only one layer particles just on the inlet nozzle;(b)The pro file sketch of reserving fluid particles is a disc for the injection process in a rectangle container with two-circle inlets(see Fig.2(a)),and the pro file of reserving fluid particles distribution is a rectangle for the filling process of in a ring container(see Fig.2(b));(c)These reserved fluid particles are set the uniform injection velocity i.e.the inlet velocity.

It is worth noting that the mode of initial fluid particles distribution for the injection process in a rectangle container(see Fig.2(a))is different from the initial particles distribution mode in[Xu et al.(2012)].The capture of jet coiling using the initial fluid particles distribution mode mentioned above is more obvious than that using the initial particles distribution mode in[Xu et al.(2012)],which can be observed from the present numerical results in subsection 4.2.

3.6　Computational performance

In the simulations of using present particles method,the ordinary differential equations(30),(31),(18)and(29)need to be solved and the physical quantities are evolved using the predictor-corrector scheme(32)-(33)and pressure updated by formulation(9).Meanwhile,the multi-wall boundaries are considered and the physical variations of rigid particles are calculated by using the proposed boundary treatment in subsection 3.5.

As well as known,a key point have to be considered and handled which is the huge calculate caused by the search of neighbor particles in the 3D computational performance.In order to enhance the computational efficiency,a neighbor particle searching technique based on dynamic background mesh following the motion of fluid particles is adopted that likes the dynamic cell method in the molecule simulation,which is detailed in[Zhou and Ge(2014);Ferrari et al.(2009);Shen and Vassalos(2011);Cleary et al.(2007)].In the whole numerical simulations,our code is written in visual C++program and the data are organized in a flexible way of searching technique of linking lists following the fluid motion.We can also know that the simulation procedure is usually suffer from the storage problems due to required million particles in 3D polymer filling process,and the dynamic arrays method of C++program is used for overcoming the storage problems.Moreover,the present particle method has some characters the same as the molecule method in which the particle method is easy to perform the Message Passing Interface(MPI)parallelization.Therefore,the MPI parallelization combined with the dynamic cell particle searching method is performed in an IBM HPC Platform with number of CPUs for better reducing the computer consumes in this paper.The detail and merits of using MPI system is described in[Ferrari et al.(2009)].

4　Tests for the validity of present scheme

In this section,the validity and capacity of the proposed CSPH-3D scheme combined with the present boundary treatment is demonstrated by solving several benchmarks and comparing with other numerical results,and the effect of density diffusive term is also be tested.Meanwhile,the merits of presented particle method over than the other methods for simulating the free surface flows are illustrated.Noting that the Oldroyd-B fluid model is considered in this section for comparisons whenα0=0,f(λ,τ)=1 in the constitutive model of XPP model(5).

4.1 fluid flows in 2D coordinate system

4.1.1Example 1:viscoelastic flow past a periodic array of cylinders

In the previous works[Ellero and Tanner(2005);Fang et al.(2006);Ren et al.(2011)],the validity and accuracy of the conventional SPH method and earlier improved SPH method for solving the New tonian or viscoelastic fluid flow problem is tested by simulating the benchmark i.e.three dimensional or two dimensional Poiseuille flow.The pro file of three dimensional problems is usually as the key studied point in the simulation of 3D bench problem,which is actually translated into the study of 2D problem.Here,the problem of Oldroyd-B fluid flow around a periodic array of cylinders con fined in a channel is simulated in 2D coordinate system,for the purpose of further illustrating the capacity and merit of the proposed CSPH method with a new boundary treatment for solving the viscoelastic fluid flow.Meanwhile,the influences of elastic stress parameters on the velocity pro files or phenomenon of vortex flow are discussed and compared with other numerical or experimental results in[Quesada and Ellero(2012);Liu et al.(1998)].

Figure 3:The 2D geometric con figuration of periodic array of cylinders:(a)periodic array of cylinders;(b)analog con figuration of single cylinder.

Figure 4:Comparisons of different numerical results for the flow at stable state along different paths(see Fig.3(b)):(a)u velocity;(b)pressure p.

(II)Subsequently,the complex rheological behavior of the Odroyd-B fluid flow past a periodic array of cylinders con fined in a channel is studied using the proposed CSPH method,and some phenomena have been numerically predicted which have not been discussed in[Quesada and Ellero(2012);Liu et al.(1998)].All the physical quantities are set as:LH=0.04m,rc=0.02m,F=5×10?5ms?2,ρ=103kgm?3,η=0.1Pa·s,theU=1.2×10?4ms?1,c=0.02ms?1,and corresponding to theRe=0.024;the fluid particles are uniform ly distributed on the whole domain,and the particle numberNy=71 along the y-axis;the ratioβ0=0.2,Weissenberg numberWe=0.06,time step Δt=5×10?3and smoothing lengthh=1.3d0.

The phenomena of velocity overshooting and vortex flow are mainly shown using the CSPH method with the ratioLr=2.5,which can be observed from the Fig.5 and Fig.6.According to the References[Quesada and Ellero(2012);Liu et al.(1998)],the vortex flow phenomenon occurs between two neighbor cylinders for the Newtonian or polymeric fluid through a periodic array of cylinders,whenLr=2.5.In addition,the overshooting phenomenon ofu-component velocity also appear for the Oldroyd-B fluid before achieving stable stage(see Fig.5),which has not been illustrated in[Quesada and Ellero(2012);Liu et al.(1998)].

Figure 5:Velocity pro files for different fluid flows along different paths at different times(Lr=2.5):New tonian fluid( first column);Oldroyd-B fluid(second column).

Figure 6:The phenomenon of vortex flow for the Oldroyd-B fluid with different Weissenberg numbers

Fig.5 shows the velocity pro files for the Newtonian and Oldroyd-B fluid flow along different paths(with differentRowhich is the ratio ofxcand cylinder radius)at different times.From Fig.5,we can find that:(1)the consumed CPU time for the Oldroyd-B fluid(t≈8s)is more than that for the Newtonian fluid flow(t≈4s)before flow achieving steady;(2)the overshooting phenomenon for the Oldoryd-B fluid is more evident than the Newtonian fluid,which is similar the case of Poiseuille flow in[Ellero and Tanner(2005)].The mentioned above remarks have just given expression to the complex non-linear character of viscoelastic fluid.Fig.6 shows the vortex flow phenomenon for the viscoelastic fluid with different Weissenberg number,which is also discussed in[Quesada and Ellero(2012);Liu et al.(1998)].The size of shaping vortex flow between two neighbor cylinders is enlarging with the increase of Weissenberg number,due to that the influence of the elastic stress on the fluid flow.Moreover,the contour distributions of shear elastic stress and first normal stress difference have been illustrated in Fig.7 withLr=2.5.The peak values ofτ′xyorcan be observed,and they appear near the cylinder at different positions,which is in agreement with the experimental results in[Liu et al.(1998)].

Figure 7:Numerical convergence for different quantities along different paths with increasing the particle numbers(Lr=6,We=0.06):(a)u velocity;(b)shear elastic stress;(c) first normal stress difference.

In order to further test the reliability of the proposed particle method for simulating the Oldroyd-B fluid flow around a periodic array of cylinders,the simple numerical convergence is investigated by increasing the particle numbersNyor decreasing the initial distance of particlesd0,which can be seen in Fig.7.Fig.7 shows the Numerical convergence for the velocity,shear elastic stress and first normal stress difference along different paths with increasing the particle numbersNyatLr=6,We=0.06,in which theLxy(see Fig.7(c))represents the path over the cylinder centerline that of a 45oto the horizontal line.The numerical results ofu,τ′xyandN′1is change little with the increase of particle numbers from observing the Fig.7,and we can believe that the CSPH results is valid and credible.Moreover,from the Fig.7,we can deduce that the proposed approach to simulate the problem of flow past a cylinder at low Reynolds number has desirable numerical convergence.

Example 2:Spin-down problem

Through the simulation of spin-down problem obtained using the present CSPH method with boundary treatment,the reliability of present boundary treatment is further verified.Meanwhile,the influences of particles distribution mode on the numerical accuracy or stability are also studied by this simulation with circle or square particles distribution.Spin-down is that the fluid con fined in a cylindrical container initially rotating with the same angular velocity then suddenly stopped,which involves a critical boundary problem and has been studied in[Ren et al.(2011)].

A ll the parameters are the same as in[Ren et al.(2011)],the number of rings is 23(3 rings boundary particles are included)and the centre of the rings isx=0,y=1.2.Fig.8 shows that the particles positions and pressure field for the spin-down problem obtained using the proposed method with different initial particles distribution modes att=0 andt=0.8 after the boundary particles were set to rest.It can be seen the result of adopting the circle participles distribution mode(Fig.8(a))att=0.8 agrees well with that in[Ren et al.(2011);],but the result of employing the square particles distribution mode(Fig.8(b))is unacceptable and undesired.From observing the Fig.8,we can get:(1)the proposed boundary treatment is credible;(2)the numerical result using the circle mode of particles distribution has higher accuracy and better stability than that using the case of square particles distribution when fluid flow con fined in a circle container.It is worth noting that the circle mode of initial particles distribution is only considered for the jet bulking problem of filling process in 3D coordinate system in this paper.

The validity and capacity of the proposed particle method have been tested by two challenge examples,and the feasibility of present boundary treatment is also verified in this subsection.The ability and merits of the proposed CSPH-3D method for simulating the surface flows are demonstrated in the follow ing subsection 4.2.

Figure 8:Particles positions and pressure field for the spin-down problem obtained using the proposed method with different particles distribution modes at t=0 and t=0.8 after the boundary particles were set to rest.

4.2　free surface flows in 3D coordinate system

Example 1:Stretching of a sphere droplet

A ll physical quantities of stretching of a sphere droplet are:the reference densityρ0=103kgm?3,the viscosityμ=10?3kgm?1s?1,and the speed of soundc=1.4×103ms?1.The initial geometry of the water drop is a circle of radiusR=1 m.There is no external forces but a initial velocity fieldand the initial pressure field??The number of fluid particles is 113081 and corresponding to the initial distanced0=0.0167 m,the smoothed lengthh=1.3d0,and the time stepdt=10?5s.

Fig.9 illustrates the particles distributions obtained using different method for the stretching of a sphere droplet att=0.01 s.The distribution of particles obtained by the CSPH-3D method are more uniformly than the TSPH-3D method,especially the external surface particles simulated by the proposed method are far smoother than the TSPH-3D.Moreover,the comparisons of consumed CPU time using different methods for a sphere droplet stretching problem with different number of CPUs are shown in Tab.1(corresponds to the case of Fig.9).From seeing the Tab.1,we can know that:(1)the computational efficiency has been enhanced for three different methods with increasing the number of CPUs;(2)the more number of CPUs the longer communication time between computer nodes in the simulation,so the reduced proportion of consumed CPU time is less than the increased proportion of number of CPUs;(3)the present CSPH-3D method has larger calculated amount than the TSPH-3D method.

Table 1:Comparisons of consumed CPU time using different methods for the stretching problem with different CPU number(corresponds to the case of Fig.9).

Example 2:Jet coiling of filling process in a rectangle container with a single inlet

The jet coiling problem of 3D filling process confined in a rectangle container with a single inlet has been numerically or experimentally studied in[Toméetal.(1999);Hwang and Kwon(2002);Tomé et al.(2008);Oishi et al.(2008)]for the Newtonian or Oldroyd-B fluid,and it is a typical and challenge example for testing the feasibility of the presented particle method applied to viscoelastic free surface flow.Moreover,the numerical simulations of 3D filling process for the polymer fluid usually suffer from two difficulties(see[Hwang and Kwon(2002);Oishi et al.(2008)]):(a)the treatment of complex multi-boundaries;(b)the enhancement of computational efficiency,especially using the particle method.Here,the merits and ability of the proposed particle method combined with presented boundary treatment and MPI parallelization technique are further illustrated by simulating the jet coiling of 3D polymer filling process based on the XPP fluid,and the results for the case of Newtonian fluid are also shown for comparison.

Figure9:The particles distributions obtained using different method for the stretching of a sphere droplet at t=0.01 s.

Figure 10:The shape and w?velocity contour of the filling process with single inlet obtained using present CSPH-3D method based on the Newtonian fluid( first column)and the XPP fluid(second column).

Figure 11:The sketch of filing process with three different containers.

Fig.10 shows the shape andw?velocity contour of the filling process with single inlet obtained using present method based on the Newtonian fluid and the XPP fluid(see[Oishi et al.(2012)]).In the simulation,the jet widthD=4 mm and the inlet is located in the center of container,the size of rectangle containerLx=Ly=5cm andLz=7cm,inlet velocityU=?1 ms?1,the stable time step is 5×10?6.The reference densityρ0=1.1×103kgm?3,kinematic viscosity?speed of soundc=8 ms?1,λ0b=0.01 s,λr=0.8,β0=0.3 andα0=0.1.The values of parameter mentioned above correspond toand the ratioIt can be seen that the deformation process of jet coiling for the Newtonian fluid(see Fig.10( first column))is basically agree with the results obtained using FDM method in[Tomé et al.(2008);Oishi et al.(2008)],and the phenomenon of jet coiling for the XPP fluid also appears which can be observed from the second column of Fig.10.There are some small differences between the results of Fig.10 and those in[Tomé et al.(2008);Oishi et al.(2008)]due to itself instability of jet coiling and the differences of numerical techniques.Through the illustration of Fig.10,we can get that:(1)it is credible for the proposed method to simulate the complex 3D viscoelastic free surface;(2)the present boundary treatment is feasible to treat multi-boundaries,and it is easier to enforce than the boundary method in[Fang et al.(2006)];(3)the tendency of coiling is more obvious than that in[Fang et al.(2006)]at the short time of fluid impact on the container bottom.Moreover,the variation ofw-velocity can also be observed by the demonstration of velocity contours in Fig.10.In this simulation,the CPUs(30 CPUs are used)consumed time is about 30 hours and 37 hours for simulating Newtonian fluid and XPP fluid jet,respectively(with 307749 fluid particles reaching 90000 steps).

In order to further demonstrate the ability of extended application of proposed CSPH-3D method to simulate the 3D polymer filling process,the polymer filling processes based on XPP fluid with different containers(see Fig.11)are numerically investigated and some instability phenomena are illustrated in the following section 5 and section 6.The sketch of filing process with three different containers is shown in Fig.11.In the following numerical investigations,two typical containers are mainly considered i.e.the rectangle container with two-inlets(see Fig.11(a))and the ring container(see Fig.11(b)-(c)).Moreover,some size parameters are introduced and labeled in Fig.11.

5　Numerical investigation of filling process in a rectangle container

5.1　Jet coiling of filling process with two-in lets

Fig.12 shows the shape andw?velocity contour of the filling process with two inlets obtained using present CSPH-3D method based on the Newtonian fluid.Meanwhile,the deformation process of the filling process based on the Newtonian fluid and XPP fluid with different times are also illustrated in Fig.13(corresponding to the case of Fig.12).In this simulations of Fig.12 and Fig.13,the size of rectangle containerLx=5 cm,Ly=4cm andLz=7cm,theRe=0.4,We=2.5,and the other parameters are the same as in Fig.10.From observing Fig.12 and Fig.13,we can get that:(1)Under the effect of gravity,the bigger value region ofw?velocity appears on a certain part between the inlet and bottom of container(see Fig.12);(2)For the reason of shear thinning of the XPP fluid,the time of fluid jets reaching the bottom is shorter for the XPP fluid than that for the Newtonian fluid(see Fig.13 att=0.075s);(3)The jet coiling phenomenon occurs for two jets;(4)Although two jets are an identical fluid and alike inject condition,the deformation process of two jets is asymmetric att=0.36s andt=0.41s due to that itself stability of jet coiling phenomenon.The asymmetric phenomenon also occurs for two XPP fluid jets after impacting on container bottom which can be seen in Fig.13 and Fig.14( first column),and it is not observed for two Newtonian fluid jets in 2D coordinate system(see[Ren et al.(2011)]).

Figure 12:The shape and w?velocity contour of the filling process with two-inlets obtained using present CSPH-3D method based on the Newtonian fluid.

5.2　Influences of macroscopic rheological parameters on the filling process

In order to illustrate the influences of macroscopic parameters on the phenomena of jet coiling and the asymmetric deformation case of two fluid jets.Figs.14-16 show that the deformation process of two XPP fluid jets with different macroscopic parametersWeandRe.Meanwhile,the variations of the jet length obtained using present method based on the XPP fluid versus dimensionless time with different macroscopic parameters are demonstrated in Fig.17(the other parameters corresponds to the case in Fig.12).For the Figs.14-16,the same parameters are:the jet widthD=0.4 cm;the height of rectangle containerLz=7cm;the inlet veloc-ity;the reference densityρ0=1.1×103kgm?3;the speed of soundc=8 ms?1;theλr=0.8,β0=0.3 andα0=0.1;the ratioThe other parameters are:(1)in Fig.14,theLx=Ly=5cm,the time-stepdt=5×10?6,the kinematic viscosity? andRe=0.4,theλ0b=0.01 s(We=2.5)orλ0b=0.4 s(We=100);(2)in Fig.15,theLx=Ly=4cm,the time stepdt=5×10?6(corresponds toυ0=5×10?3m2s?1)ordt=1×10?5(corresponds toυ0=1×10?3m2s?1),the kinematic viscosityυ0=5×10?3m2s?1(Re=0.8)orυ0=1×10?3m2s?1(Re=4),theλ0b=0.01 s(We=2.5)orλ0b=0.1 s(We=25);(3)in Fig.16,theLx=Ly=4cm,the time-stepdt=1×10?5,the kinematic viscosityυ0=5×10?4m2s?1(Re=8),theλ0b=0.01 s(We=2.5).

Figure 13:The deformation process of the filling process with two-inlets based on the Newtonian fluid( first column)and XPP fluid(second column).(Corresponds to the case of Fig.12).

Figure 14:The shape and w?velocity contour of the filling process with two-inlets obtained using present CSPH-3D method based on the XPP fluid.18CPUs We=2.5( first column),We=100(second column).

Figure 15:The shape of filling process obtained using present CSPH-3D method based on the XPP fluid with two-inlets:Re=0.8,We=25( first column),Re=4,We=2.5(second column).

Figure 16:The shape and w?velocity contour of the filling process with two-inlets obtained using present CSPH-3D method based on the XPP fluid(Re=8,We=2.5).

Figure 17:The variations of the jet length obtained using present method based on the XPP fluid versus dimensionless time with different model parameters(the other parameters corresponds to the case in Fig.12).

It can be seen that the phenomenon of jet coiling doesn’t occur when the Weissenberg number is increased to 100 from Fig.14,due to that the higherWelead to the bigger effect of shear thinning of XPP fluid(see[Oishi et al.(2012)]),and the deformation of two jets become symmetric atWe=100.With increasing Weissenberg number is,the longer jet length is achieved before fluid jets reaching the container bottom,which also can seen in Fig.17(a).Meanwhile,the variations ofw?velocity are illustrated at different times in Fig.14,and some fluid particles near the bottom after jets impacting on the container bottom are not moving.

Through observing the Fig.15 and Fig.16,we can get that:(1)the phenomenon of jet coiling is disappear when the Reynolds number is increased;(2)the thinner disk near the container bottom is gained with biggerRe,which can be seen in Fig.15 att=0.12s(Re=4)and Fig.16 att=0.09s(Re=8);(3)the filling process is stable whenRe=0.8 andWe=25(see Fig.15( first column));(4)the phenomenon of concave near the fluid jet appear after two fluid jets collided and the filling process is unstable withRe=4 andRe=8;(5)the bigger value region ofw?velocity is decreased versus time after two jets impacting on the bottom withRe=8(see Fig.16);(6)with increasing the Reynolds number(it increased from 0.8 to 8),the shape of two fluid jets has well symmetric at different times even if the filling process is unstable;(7)the longer jet length is gained before fluid jets reaching the container bottom when the Reynolds number is increased,which also can seen in Fig.17(b).Moreover,we can know that the shorter jet length is achieved with increasing the ratioβ0due to that the shear thinning become weak.

In a word,the influences ofWeandReon the deformation of filling process are important.The phenomenon of jet coiling may be occur when the physical parameters are adopted appropriately for the XPP fluid,and the deformation process of two fluid jets is more complex than that of a single jet.The asymmetric phenomenon of two fluid jets is observed in the simulations of 3D filling process at lower Reynolds number(for exampleRe=0.4).Moreover,Comparing the Oishi’s works[Oishi et al.(2012);Figueiredo et al.(2013)]with the present work,the main differences or novelties lie in two points:1)the deformation process of two XPP jets coiling is asymmetric at some time which can be not observed in the case of a single jet process;2)the required time of container filled for two jets is less than that for single jet in the stable filling process.

6　Further numerical investigation of filling process in a ring container

In this section,the 3D filling processes of die casting with two different form ring containers(see Fig.11(b)-(c))are numerically predicted using the proposed CSPH-3D method.The closed ring container has more complex multi-boundaries than that of the closed rectangle container in section 5,so the proposed particle method combined with present boundary treatment can be further tested by simulating the 3D filling process of die casting with ring container.Indeed,the only multi-straight-boundaries(see Fig.1(a))are treated in the simulations of filling process in a rectangle container(see Fig.11(a)),and the multi-straight-boundaries and multi-curve-boundaries(see Fig.1)are needed to treat in the simulations of filling process in a ring container(see Fig.11(b)-(c)).In the previous work[Ren et al.(2011)],the 2D filling process in a ring cavity based on the Cross fluid is numerically investigated using an improved particle method.Here,the 3D filling process in a ring container based on the XPP fluid is numerically studied,and the case of Newtonian fluid is also considered for comparisons.It is noting that the main differences of proposed CSPH-3D method over than the improved SPH method in[Ren et al.(2011);Fetehi and Manzari(2012)]lie in:(1)the discretized form of momentum equation is more adapted to simulate unstable filling process than the case in[Ren et al.(2011);Fetehi and Manzari(2012)];(2)the proposed boundary treatment is more easier performed than that in[Fang et al.(2006);Ren et al.(2011);Xu et al.(2012)].

Figure 18:The deformation process of die casting in ring container obtained using the present CSPH-3D method based on the Newtonian fluid( first column)and XPP fluid(second column)with Re=90.

Figure 19:The u?velocity and pressure distribution of die casting in ring container based on the Newtonian fluid( first column)and XPP fluid(second column)with different times,which corresponding to the case of Fig.18.

The unstable filling processes are demonstrated in Figs.18-19,in which the deformation process of die casting in a ring container(see Fig.11(b))based on the Newtonian fluid and XPP fluid withRe=90(see Fig.18).The center of internal circle is located in the center of big ring in the ring containers of Fig.11(b)-(c).The physical parameters of Fig.18 are:the size of container?? andrb=1.6L(L=4.5cm);the inject velocityU=4 ms?1,the stable time-step is 1×10?5;the reference densitykgm?3,kinematic viscosity?,speed of soundc=8U,λ0b=0.02 s,? andthe values of parameter mentioned above correspond to??From observing the Fig.18 and 19,it can be seen that:(1)the filling process is unstable,and the numerical results are basically similar to the results in[Ren et al.(2011)]but there are some differences because that the present inject fluid is also affected by the upper-down rigid walls in a 3D cavity;(2)the fracture phenomena near the front surface are observed in the larger inject velocity for the Newtonian and XPP fluid at aboutt=0.05s;(3)the re flux phenomenon for the XPP fluid is more obvious than that for the Newtonian fluid due to the effect of shear thinning of XPP fluid,and the required time is shorter for the XPP fluid than that for the Newtonian fluid when the container is full;(4)the thickness(which alongz?axis)near front fluid region is thinner than the initial thickness of fluid due to the effect of upper-down rigid walls,which is agree with the results in[Hwang and Kwon(2002)],and it also can be observed in Fig.20;(5)the largeru?velocity value region happen to near the middle part of container along thex?axis before the fluid impacting on the internal cylinderrs(see Fig.19(a)at);(6)theu?velocity distribution is nonuniform after the fluid impacting on the external ringrbdue to the instability of filling process;(7)the larger pressure value is increased at short time after the fluid impacting on the internal cylinder(see Fig.19(b)),and it is decreased versus time.of Fig.20 are:the size of containerand rb=1.15L(L=2cm);the inject velocity U=6 ms?1,the stable time-step is 2×10?6;the reference density ρ0=1.0×103kgm?3,kinematic viscosity10?3m2s?1,speed of soundand α0=0.1;the values of parameter mentioned above correspond toThe 3D numerical results of Fig.20 is similar with the 2D results in[Ren et al.(2011)],but the feature of point(4)mentioned above in Fig.18 is observed only for3D case.Though the deformation process for the Newtonian fluid is basically similar with that for the XPP fluid,the moving front free surface for the XPP fluid flow is faster than that for the Newtonian fluid after the fluid impacting on the internal cylinder due to the effect of shear thinning for the XPP fluid.Moreover,the front free surface likes a convexity,and the track effect is obviously observed about t=0.025s in Fig.20.

Fig.20 shows that the deformation process of ring filling container(see Fig.11(c))based on the Newtonian fluid and XPP fluid with Re=12.The physical parameters

Through the numerical investigations mentioned above,we can know that:(1)the present boundary treatment is valid and convenient to treat the complex multiboundaries;(2)some complex phenomena are observed in the numerical simulations of 3D polymer filling process,but not seen in the 2D case,for example the jet coiling and the thickness of front surface region is thinner(see Fig.20 att=0.006s);(3)the polymer filling process is more complex than that for the Newtonian fluid filling due to the characteristics of the fluid flow itself.

7　Some conclusions

In this paper,a corrected 3D parallel SPH(CSPH-3D)method is first proposed to simulate the polymer filling process in a rectangle or ring container for the XPP fluid,and some complex instability phenomena are also numerically investigated.The present CSPH-3D method is based on a coupled concept that an extended kernel-gradient-corrected SPH method is used in the interior of fluid region and the traditional SPH(TSPH)method is used near the boundary domain,which has higher accuracy and better stability than the TSPH-3D.In order to restrain the pressure oscillations and treat the complex multi-boundaries,a density diffusive term is introduced and a new boundary treatment is proposed,respectively.Moreover,the MPI parallelization means with a dynamic cells neighbor particle search method is adopted for enhancing the computational efficiency.Subsequently,several bench tests are solved using the proposed particle method for demonstrating the validity and merits of proposed CSPH method,and some remarks can be obtained which are:(1)it is reliable and accurate for the CSPH method to solve the viscoelastic fluid flow;(2)the present CSPH method has higher accuracy and better stability than the TSPH method;(3)it is it is validate and feasible for the proposed CSPH-3D method to capture the viscoelastic free surface;(4)the proposed boundary treatment is feasible and convenient to treat straight or curve boundary.

Subsequently,the deformation processes in a rectangle container filling with two inlets and in a ring cavity molding for the XPP fluid are numerically predicted using the CSPH-3D method,and the influences of rheological parameters on the filling process are studied.All the numerical results show that:(1)the present boundary treatment is valid and convenient to treat complex multi-boundaries;(2)some complex phenomena are observed in the numerical simulations of 3D polymer filling process,but not seen in the 2D case,for example the jet coiling and the asymmetric character of two jets;(3)the polymer filling process based on the XPP fluid is more complex than that for the Newtonian fluid filling due to the characteristics of the fluid flow itself;(4)the influences ofWeandReon the deformation of filling process are important;(5)the phenomenon of track effect is observed in the 3D polymer filling process in a ring container.The proposed CSPH-3D method is further expected to be extended and applied to more complex fluid dynamic problems.

Acknowledgement:All the authors acknowledge the support from he Postdoctoral Science Foundation of China(Grant No.2014M 550310),the Natural Science Foundation of Jiangsu Province,China(Grant No.BK20130436),and the Innovation Cultivation Funds of Yangzhou University,China(Grant No.2013CXJ003).

Antuono,M.;Caolagrossi,A.;Marrone,S.;Molteni,D.(2010):Free-surface flows solved by means of SPH schemes with numerical diffusive terms.Computers Physics Communications,vol.181 pp.532–549

Bonet,J.;Lok,T.(1999):Variational and momentum preservation aspects of smooth particle hydrodynamic formulations.Computer Methods in Applied Mechanics and Engineering,vol.180,no.(1–2),pp.97–115.

Batra,R.C.;Zhang,G.M.(2008):SSPH basis functions for meshless methods,and comparison of solution of solutions with strong and weak formulations.Computational Mechanics,vol.41,pp.527–545

Cleary,P.W.;Ha,J.;Prakash,M.;Nguyen,T.(2007):3D SPH flow predictions and validation for high pressure die casting of automotive componentsApplied Mathematical Modelling,vol.30,pp.1406–1427.

Colagrossi,A.;Landrini,M.(2003):Numerical simulation of interfacial flows by smoothed particle hydrodynamics.Journal of Computational Physics,vol.191,pp.448–475.

Chen,J.K.;Beraun,J.E.(2000):A generalized smoothed particle hydrodynamics method for nonlinear dynam ic problems.Computer Methods in Applied Mechanics and Engineering,vol.190,pp.225–239.

Ellero,M.;Tanner,R.I.(2005):SPH simulations of transient viscoelastic flows at low Reynolds number.J.Non-Newtonian Fluid Mech.,vol.132,pp.61–72.

Ferrari,A.;Dumbser,M.;Toro,E.F.;Armanini,A.(2009):A new 3D parallel SPH scheme for free surface flows.Computers&Fluids,vol.38 pp.1203–1217

Fetehi,R.;Manzari,M.T.(2012):A consistent and fast weakly compressible smoothed particle hydrodynamics with a new wall boundary condition.Int.J.Num.Meth.Flu.,vol.68,pp.905–921.

Figueiredo,R.A.;Oishi,C.M.;Cuminato,J.A.;Alves,M.A.(2013):Three dimensional transient complex free surface flows:Numerical simulation of XPP fluid.Journal of Non-Newtonian Fluid Mechanics,vol.195,pp.88–98.

Fang,J.N.;Owens,R.G.;Tacher,L.;Parriaux,A.(2006):A numerical study of the SPH method for simulating transient viscoelastic free surface flows.J.Non-Newtonian Fluid Mech.,vol.13,pp.68–84.

Hu,X.;Adam s,N.(2006):A multi-phase SPH method for macroscopic and mesoscopic flows.Journal of Computational Physics,vol.213,pp.844–861

Hwang,C.J.;Kwon,T.H.(2002):A Full 3D Finite Element Analysis of the Powder Injection Molding Filling Process Including Slip Phenomena.Polymer Engineering and Science,vol.42,pp.33–50.

Monaghan,J.J.(1994):Simulating free surface?ows with SPH.J.Comp.Phys.,vol.110,pp.399–406.

Jiang,T.;Ouyang,J.;Ren,J.L.(2012):A mixed corrected symmetric SPH(MCSSPH)method for computational dynamic problems.Comp.Phys.Comm.,vol.183,pp.50–62.

Liu,A.W.;Bornside,D.E.;Arm strong,R.C.;Brown,R.A.(1998):Viscoelastic flow of polymer solutions around a periodic,linear array of cylinders:comparisons of predictions for microstructure and flow fields.Journal of Non-Newtonian Fluid Mechanics,vol.77,pp.153–190.

Liu,M.B.;Liu,G.R.(2006):Restoring particle consistency in smoothed particle hydrodynamics.Applied Numerical Mathematics,vol.56,pp.19–36.

Li,S.F.;Liu,W.K.(2002):Mesh-free particle methods and their applications.Applied Mechanics Review,vol.54,pp.1-34.

Li,Q.;Ouyang,J.;Li,X.J.;Wu,G.R.;Yang,B.X.(2011):Numerical simulation of gas-assisted injection molding process for a door handle.CMES:Computer Modeling in Engineering&Sciences,vol.74,pp.247–267.

Liu,M.B.;Xie,W.P.;Liu,G.R.(2005):Modeling incompressible flows using a finite particle method.Applied Mathematical Modelling,vol.29,pp.1252–1270.

Marrone,S.;Antuono,M.;Colagrossi,A.;olicchio,G.;Touze,D.L.;Graziani,G.(2011):δ-SPH model for simulating violent impact flows.Comput.Meth.Appl.Mech.Engrg.,vol.200,pp.1526–1542.

Morris,J.P.;Fox,P.J.;Zhu,Y.(1997):Modeling low Reynolds number incompressible flows using SPH.J.Comput.Phys.,vol.136,pp.214–226.

Messahel,R.;Souli,M.(2013):SPH and ALE Formulations for Fluid Structure Coupling.CMES:Computer Modeling in Engineering&Sciences,vol.96,pp.435–455.

Oishi,C.M.;Martins,F.P.;Tomé,M.F.;Alves M.A.(2012):Numerical simulation of drop impact and jet buckling problems using the eXtended Pom-Pom model.J.Non-Newtonian Fluid Mech.,vol.169-170,pp.91–103.

Oishi,C.M.;Tome′M.F.;Cuminato,J.A.;McKee,S.(2008):An implicit technique for solving 3D low Reynolds number moving free surface flows.Journal of Non-Newtoian Fluid Mechanics,vol.227 pp.7446–7468

Quesada,A.V.;Ellero,M.(2012):SPH simulations of a viscoelastic flow around a periodic array of cylinders con fined in a channel.Journal of Non-Newtonian Fluid Mechanics,vol.168,pp.1–8.

Ren,J.L.;Ouyang,J.;Jiang,T.;Li,Q.(2011):Simulation of complex filling process based on the generalized New tonian model using a corrected SPH scheme.Computational Mechanics,vol.49,pp.643–665.

Ren,J.L.;Ouyang,J.;Yang,B.X.;Jiang,T.;M ai,H.Y.(2011):Simulation of container filling process with two inlets by improved smoothed particle hydrodynamics(SPH)method.International Journal of Computational Fluid Dynamics,vol.25,pp.365–386.

Skillen,A.;Lind,S.;Stansby,P.K.;Rogers,B.D.(2013):Incompressible smoothed particle hydrodynamics(SPH)with reduced temporal noise and generalised Fickian smoothing applied to body–water slam and efficient wave–body interaction.Computer Methods in Applied Mechanics and Engineering,vol.265,pp.163–173.

Sladek,J.;Stanak,P.;Han,Z-D.;Sladek,V.;Atluri,S.N.(2013):Applications of the MLPG Method in Engineering&Sciences:A Review.CMES:Computer Modeling in Engineering&Sciences,vol.92,pp.423–475.

Shen,L.;Vasaalos,D.(2011):Applications of 3D Parallel SPH for Sloshing and Flooding.Contemporary Ideas on Ship Stability and Capsizing in Waves Fluid Mechanics and Its Applications,vol.97,pp.709–721.

Stranex,T.;Wheaton,S.(2011):A new corrective scheme for SPH.Computer Methods in Applied Mechanics and Engineering,vol.200,pp.392–402.

Tomé,M.F.;Castelo,A.;Ferreira,V.G.;M cKee,S.(2008):A finite difference technique for solving the Oldroyd-B model for 3D-unsteady free surface flows.J.Non-Newtonian Fluid Mech.,vol.154,pp.179–206.

Tomé,M.F.;M ckee,S.;Barratt,L.;Jarvis,D.A.;Patrick,A.J.(1999):An experimental and numerical investigation of container filling with viscous liquids.Int.J.Numer.Meth.Fluids,vol.31,pp.1333–1353.

Tomé,M.F.;Mangiavacchi,N.;Castelo,A.;Cuminato,J.A.;McKee,S.(2002):A finite difference technique for simulating unsteady viscoelastic free surface flows.J.Non-Newtonian Fluid Mech.,vol.106,pp.61–106.

Vebeeten,W.M.H.;Peters,G.W.M.;Baaijens,F.P.T.(2001):Differential constitutive equations for polymer melts:The extended Pom-Pom model.J.Rheol.,vol.45,pp.823-843.

Wend land,H.(1995):Piecewise polynomial,positive definite and compactly supported radial functions of minimal degree.Advances in Computational Mathemat-ics,vol.4,pp.389–396.

Xu,X.Y.;Ouyang,J.;Jiang,T.;Li,Q.(2012):Numerical simulation of 3D-unsteady viscoelastic free surface flows by improved smoothed particle hydrodynamics method.J.Non-Newtoian Fluid Mechanics,vol.177-178 pp.109–120

Zhang,G.M.;Batra,R.C.(2004):Modified smoothed particle hydrodynamics method and its application to transient problems.Comput.Mech.,vol.34,pp.137–146.

Zhang,G.M.;Batra,R.C.(2007):Wave propagation in functionally graded materials by modified smoothed particle hydrodynamics(MSPH)method.J.Comput.Phys.,vol.222,pp.374–390.

Zhou,G.Z.;Ge,W.(2014):Progress of smoothed particle hydrodynamics in complex flows.CIESC Journal,vol.65,pp.1145–1161.

Zainali,A.;Tofighi,N.;Shadloo,M.S.;Yildiz,M.(2013):Numerical investigation of Newtonian and non-Newtonian multiphase flows using ISPH method.Computer Methods in Applied Mechanics and Engineering,vol.254,pp.99–113