† Corresponding author. E-mail:

Project supported by the Natural Science Foundation of Jiangsu Province, China (Grant Nos. BK20130436 and BK20150436) and the Natural Science Foundation of the Higher Education Institutions of Jiangsu Province, China (Grant No. 15KJB110025).

In this work, the behavior of the three-dimensional (3D) jet coiling based on the viscoelastic Oldroyd-B model is investigated by a corrected particle scheme, which is named the smoothed particle hydrodynamics with corrected symmetric kernel gradient and shifting particle technique (SPH_CS_SP) method. The accuracy and stability of SPH_CS_SP method is first tested by solving Poiseuille flow and Taylor–Green flow. Then the capacity for the SPH_CS_SP method to solve the viscoelastic fluid is verified by the polymer flow through a periodic array of cylinders. Moreover, the convergence of the SPH_CS_SP method is also investigated. Finally, the proposed method is further applied to the 3D viscoelastic jet coiling problem, and the influences of macroscopic parameters on the jet coiling are discussed. The numerical results show that the SPH_CS_SP method has higher accuracy and better stability than the traditional SPH method and other corrected SPH method, and can improve the tensile instability.

The complex jet coiling phenomena,^{[1]} although they play an important role in many industrial processes, such as extrusion, and container filling with polymers, are still less well understood in terms of fundamental fluid dynamics. In order to predict the complex rheological phenomena in a real three-dimensional (3D) viscoelastic jet buckling process, a number of grid-based numerical methods have been proposed, such as the finite difference method (FDM), finite element method (FEM), finite volume method (FVM), etc. The mesh-based methods mentioned above have been applied to many fluid mechanic fields.

However, these grid-based methods also suffer numerical difficulties, such as severe mesh distortion, and the expensive cost in the mesh generation process. As one may expect, the various mesh-free methods^{[2–4]} in a Lagrangian framework have been proposed to solve the fluid dynamic problems, such as the smoothed particle hydrodynamics (SPH) method,^{[4,5]} element-free Galerkin (EFG) method,^{[3]} reproducing kernel particle method (RKPM),^{[6–8]} and so on.

The SPH method^{[5]} is a pure meshless method which employs Lagrangian description of motion. Recently, the SPH method has been extensively applied to the fluid dynamic areas due to its merits, such as viscous flows,^{[9]} free surface flows,^{[5,10,11]} incompressible fluids,^{[12,13]} multi-phase flows,^{[14]} and non-Newtonian fluid flows.^{[13,15]} However, the application of the SPH method to the viscoelastic fluid mainly focuses on the two-dimensional (2D) viscoelastic jet buckling, and the application of the SPH method to the 3D viscoelastic jet coiling has received a little attention recently. The reason is that the SPH method usually suffers two major drawbacks, that is, the particle inconsistency^{[16,17]} for uninformed particle distributions/near the boundary domain and the numerical/tensile instability.^{[18]} In order to restore the consistency of the kernel and gradient particle approximations of SPH, some improved SPH methods based on the Taylor series expansion have been proposed, such as the corrected smoothed particle hydrodynamics method (CSPM),^{[19]} finite particle method (FPM),^{[11,20]} the modified SPH (MSPH) method,^{[21]} and the symmetric SPH (SSPH) method.^{[22]} Although these improved methods have higher accuracy than the traditional SPH (TSPH) method, some disadvantages still exist and make it more difficult to be extended to the complex fluid dynamic problems,^{[22]} the choice of the kernel function is limited for solving the partial differential equations (PDEs) with higher derivatives; the local matrix is asymmetric, which can cause numerical instability; its computational efficiency is lower than that of TSPH, such as the MSPH method and SSPH method. Therefore, a corrected kernel gradient scheme^{[23]} has been introduced into the TSPH method and extensively applied to Newtonian fluid^{[24,25]} in recent years. However, the corrected kernel gradient scheme has been rarely extended to simulate the viscoelastic jet buckling, especially the 3D simulations of the viscoelastic jet coiling.

In the present paper, the 3D viscoelastic jet coiling process based on the Oldroyd-B model is simulated and analyzed by a parallel coupled corrected particle scheme named the SPH_CS_SP (SPH with a corrected symmetric kernel gradient scheme and shifting particle technique) method. Some special phenomena of the 3D jet coiling are observed and the physical explanations of these complex phenomena are also analyzed. The proposed particle method combines three improved concepts, that is, introducing the incompressible condition in the momentum equation to achieve better robustness,^{[26]} extending the corrected symmetric kernel gradient scheme^{[27]} to further improve the numerical accuracy and employing the shifting particle technique to overcome the tensile instability.^{[28]} It should be noted that the main difference between the proposed SPH_CS_SP method and other improved SPH methods lies in the following respects: (i) the obtained discretized form of momentum equation is different from those in Refs. [24] and [25]; (ii) the shifting particle technique is first applied to the viscoelastic fluid flows; (iii) the presented corrected SPH method is different from the particle methods in Refs. [24] and [29], which have not been extended to simulate the 3D polymer free surface flow.

The rest of this paper is organized as follows. The governing equations and the viscoelastic constitutive model are introduced in Section 2. In Section 3, the details of the proposed SPH_CS_SP scheme for the viscoelastic fluid are described. The accuracy and robustness of SPH_CS are tested by simulating the Poiseuille flow in Section 4. Moreover, the necessity of enforcing the shifting particles technique is also investigated by simulating the Taylor–Green flow in Section 4. In Section 5, the capacity for SPH_CS_SP to simulate the viscoelastic flow is verified by the polymer flow around a periodic array of cylinders. In Section 6, the proposed particle method with MPI parallelization technique is extended to simulate the 3D viscoelastic jet coiling problem in the filling process. Some conclusions are drawn from the present study in Section 7.

In the 3D Lagrangian frame, the governing equations of a weakly compressible fluid can be described as follows:

*υ*= (

*u*,

*v*,

*w*) is the velocity vector,

*ρ*is the fluid density,

*t*is the time,

*is the body force, and d/d*F

*t*is the material derivative d/d

*t*=

*∂*/

*∂t*+

*·∇. The Cauchy stress tensor*υ

*in Eq. (*σ

*p*and extra stress tensor

*T*

_{o}:

*refers to the identity matrix.*I

In this paper, we choose the Oldroyd-B model as the constitutive equation, which represents the basic viscoelastic fluid and shows complex rheological behavior of polymer flow. The extra stress tensor _{o} for the Oldroyd-B model is usually decomposed into viscous 2*η*_{s}* D* and polymeric contributions

*:*τ

*λ*

_{0b}is the relaxation time, and the symbol

Moreover, some parameters are introduced, that is, the total viscosity *η* = *η*_{s} + *η*_{p}, the ratio of Newtonian and polymeric viscosity *β*_{0} = *η*_{s}/(*η*_{s} + *η*_{p}), the Reynolds number *Re* = (*ρUL*)/*η*, and the Weissenberg number *We* = (*λ*_{0b}*U*)/*L*, where *U* and *L* represent the characteristic velocity and length, respectively.

Like many previous studies,^{[4,9]} here we also treat the incompressible flows as a slightly compressible fluid by a suitable equation of state, that is^{[10]}

*c*is the speed of sound, and

*ρ*

_{0}is the reference density. It can be shown that the density variation is proportional to the Mach number

*M*(

*M*≡

*V*/

*c*, where

*V*is a typical reference velocity).

^{[4,9]}

In order to accurately simulate the viscoelastic fluid flow, a new corrected particle scheme is proposed. It mainly combines three different improvements, that is, introducing the incompressible condition in the momentum equation to achieve better robustness,^{[26]} applying the corrected symmetric kernel gradient method^{[29]} to the discrete form of the first derivative to improve the numerical accuracy, and employing the shifting particle technique to overcome the tensile instability.^{[25,26,28]}

As is well known, in the traditional SPH method, the particle approximation for a function *f* and its first derivative at particle *i* can be written as^{[4,5]}

*m*and

_{j}*ρ*are the mass and the density of the neighbor particle

_{j}*j*located in the support domain of particle

*i*, respectively,

*W*is the kernel function which usually satisfies some properties,

^{[4]}

*W*=

_{ij}*W*(|

*−*r

_{i}*|,*r

_{j}*h*),

*∂W*/

_{ij}*∂*r =

_{i}*∂W*(|

*−*r

_{i}*|,*r

_{j}*h*)/

*∂*r ,

_{i}*h*being the smoothing length.

The kernel function is related not only to the accuracy but also to the efficiency and the stability of the resulting algorithm. As shown in Ref. [29], the Wendland function can produce more accurate results than the common spline kernel functions. So the following quintic Wendland kernel function is adopted:

*r*= |

*−*r

_{i}*| and the normalization factor*r

_{j}*w*

_{0}is 7/(

*πh*

^{2}.)

According to Eq. (

Considering the discrete gradient scheme

*v*

^{β}is the

*β*-th component of the fluid velocity,

*σ*

^{α β}is the (

*α*,

*β*)-th component of the total stress tensor, and

*x*is the component of spatial coordinate.

^{β}The particle discretized scheme of total stress equation (

The TSPH discretized formulas (^{[15]} As is well known, the TSPH formula for the viscous flow suffers low numerical accuracy and instability, especially for the non-uniform distributed particles and boundary particles (see Refs. [4], [11], and [16]). Moreover, the influence of elastic stress on the numerical simulations of viscoelastic fluid flow is also significant, which easily leads to the instability of numerical schemes for the viscoelastic flow. Therefore, in order to obtain a particle method that possesses higher accuracy and better stability, we introduce three improvements^{[27,28]} into the TSPH formulas (

Firstly, introducing the incompressible condition into momentum equation (

^{[13]}the following discretized scheme of the momentum equation (

*=*r

_{ij}*−*r

_{i}*.*r

_{j}Secondly, the corrected symmetric kernel gradient scheme^{[27]} and the Taylor expansion concept^{[18]} are also applied to the first-order derivative and the elastic stress part of Eq. (

It should be noted that the discretized form Eq. (

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

Finally, a shifting particle technique,^{[28]} which is applied only to the viscous fluid now, is extended to the viscoelastic fluid to obtain a uniform particle distribution and further improve the tensile instability. Next, we describe the shifting particle technique in detail.

After each time-step, the position of particle *i* is shifted according to the following formula:

*ξ*is a constant that ranges from 0.01 to 0.1,

*U*

_{max}is the maximum particle velocity or the characteristic velocity of flow, Δ

*t*is the time-step, and the shifting vector is defined as

*M*is the number of neighboring particles of the particle

_{i}*i*,

*r*= |

_{ij}*|, and*r

_{ij}*i*. The shifting vector

*i*are uniformly distributed.

For simplicity, the method that combines the three different improvements is named the SPH_CS_SP method.

In the SPH simulations, the boundary treatment is an important problem. Next, we will describe the boundary treatment used in this paper.

For the solid boundary, the approach of two types of virtual particles^{[27]} is extended to simulating the viscoelastic fluid flows in this paper. The first type of particle is placed on the solid wall, namely “wall particle”. The pressure and the components of the velocity of the wall particle *i* are calculated according to the following formula:

*j*represents the indices of the fluid particles within the influence domain of the boundary particle

*i*.

The second type of particle, called “dummy particle”, is placed just outside the solid wall, in a way which is similar to how the wall particle is located. In general, three layers of dummy particles are chosen. The pressure, velocity and extra stress tensor on the dummy particles are computed by the following formula:

*represents the vector of variables*S

*d*and

_{D}*d*are the normal distances of the dummy particles

_{F}*D*and the fluid particle

*F*to solid wall, respectively. To specify the values of

*, the interpolation formula (*S

_{F}^{[15]}

In order to solve the system of ordinary differential equations (^{[27]} which has second-order accuracy and better stability. The predictor step consists of a Eulerian explicit evaluation of all quantities for each particle

*represents the vector of the unknown variables*K

_{i}*denotes the vector of the right-hand side of each of Eqs. (*Γ

_{i}As is well known, the time-step and space-step must satisfy the well-known Courant–Friedrichs–Lewy (CFL) condition for ensuring the numerical stability. According to Refs. [14] and [30], we may choose the following stability condition:

*υ*

_{0}=

*η*/

*ρ*

_{0}is the kinematic viscosity.

In this section, the Poiseuille flow of Oldroyd-B fluid is simulated by SPH_CS_CP method to test the capacity for the proposed SPH_CS method to solve the transient viscoelastic fluid, and the numerical results are compared with those obtained by TSPH method and the analytical solutions. It should be noted that the tensile instability does not exist in the simulations of the Poiseuille flow, and the shifting particle technique is not considered in this subsection.

Figure *u* and the stress *τ ^{xy}* profiles at different times for the Poiseuille flow of Oldroyd-B fluid. Here the uniformed particle distribution is used. The employed non-dimensional parameters are as follows: two parallels are located at

*y*= 0 and

*y*= 1, the density

*ρ*

_{0}= 1, the characteristic velocity

*β*

_{0}= 0.3, the kinematic viscosity

*υ*

_{0}= 2, the body force

*F*= 1,

*Re*≈ 0.03, the relaxation time,

*λ*

_{ob}= 0.2, 4 which correspond to

*We*= 0.012, 0.25, respectively, the sound speed

*c*= 10

*U*

_{max}, the time-step Δ

*t*= 10

^{−4}, the smoothing length

*h*= 2.5

*d*

_{0}(

*d*

_{0}is the initial distance between two particles), the number of fluid particles is 31 × 61.

Table *L*_{1}-error norms of the velocity profiles in Fig. *L*_{1}-error norm is defined as

*N*denotes the particle number along the

_{y}*y*axis,

*ũ*and

_{j}*u*are the numerical and the analytical solution at position

_{j}*x*, respectively.

_{j}From Fig.

Figure *We* and *Re*. In Fig. *υ*_{0} = 0.25, *F* = 0.5, *Re* = 1, *λ*_{ob} = 0.04, and *We* = 0.01. From Fig. *t* = 0.6 (see Fig.

The numerical results in Figs.

In this subsection, the 2D Taylor–Green flow^{[28]} of the Newtonian fluid is simulated to show the necessity of the shifting particle technique.

The Taylor–Green flow is a periodic array of decaying vortices in the plane, and its velocity components are given

*U*= 1 m/s, the computational periodic domain [0,1] × [0,1], the kinematic viscosity

*υ*

_{0}= 0.01 m

^{2}/s, Reynolds number

*Re*=

*UL*/

*υ*

_{0}= 100, and the decay rate of the velocity field

*b*= −8

*π*

^{2}/

*Re*.

Figure

Figure *x*/*L* = 0.0 with *Re* = 100, *t* = 0.2 s. The SPH_CS_SP results of velocity component *u*(*y*) are much closer to the analytical solutions than the other numerical results, and small oscillations appear in the TSPH results.

The numerical results show that the proposed SPH_CS_SP method indeed has higher accuracy and better stability than the TSPH, TSPH_SP and SPH_CS method, and it is necessary to combine the shifting particles technique with the corrected symmetric kernel gradient scheme.

The polymer fluid flow around a periodic array of cylinders^{[30–32]} is usually regarded as a challenging test for the particle method. Here, we use it to test the reliability of the proposed SPH_CS_SP method to simulate the viscoelastic fluid.

Figure *F* is an effective body force, *r*_{c} is the cylinder radius, *L*_{H} is the half height of two parallel channels, *L*_{c} is the distance between two neighbor cylinders, and *x*_{c} is the distance from the positive direction of *x* axis to the center of reference cylinder (see Fig. *Ro*=*x*_{c}/*r*_{c}, *L _{r}* =

*L*

_{c}/

*r*

_{c},

*Re*= (

*ρUr*

_{c})/

*η*and

*We*= (

*λ*

_{ob}

*U*)/

*r*

_{c},

*U*′ =

*u*/10

^{−km},

*τ*′

^{αβ}=

*τ*

^{αβ}/10

^{−km},

*P*′ = p/10

^{−kn}(where

*km*and

*kn*are the integers of scientific computing,

*km*= 4,

*kn*= 3 in these simulations), the first normal stress difference

In order to test the reliability of the SPH_CS_SP method of simulating the flow past cylinders, the same example of Newtonian fluid as that in Ref. [9] is first simulated by the SPH_CS_SP method. All the physical parameters are the same as those in Ref. [9], i.e., *L*_{H} = 0.05 m, *r*_{c} = 0.02 m, *F* = 5 × 10^{−5} m·s^{−2}, *ρ* = 10^{3} kg·m^{−3}, *η* = 0.1 Pa·s, *U* = 1.5 × 10^{−4} m·s^{−1}, *c* = 0.01 m·s^{−1}, and *Re* = 0.03; the fluid particles are uniformly distributed in the whole domain, and the particle number *N _{y}* = 81 along the

*y*axis. The time step Δ

*t*= 5 × 10

^{−3}s, smoothing length

*h*= 2.7

*d*

_{0}.

Figure *U* and the pressure at stable state along different paths (see Fig.

Figure *U* and pressure obtained by the SPH_CS_SP method for the flow at the stable state under different boundary conditions (i.e., the flow is confined or not confined in a channel which corresponds to the Newman boundary condition or non-slip wall boundary condition, respectively). Comparing Fig.

In order to demonstrate the capacity of the SPH_CS_SP approach for simulating the viscoelastic flow past a periodic array of cylinders, the Oldroyd-B fluid flow past a periodic array of cylinders is simulated and the particle distribution at stable state is shown in Fig. *L _{r}* = 6,

*L*

_{H}= 0.04 m,

*U*= 1.2 × 10

^{−4}m·s

^{−1},

*c*= 0.02 m·s

^{−1},

*Re*= 0.024,

*We*= 0.9.

*N*= 71, and

_{y}*h*= 2.7

*d*

_{0}. The other parameters are the same as those in Fig.

Moreover, the following drag^{[32]} for the viscoelastic fluid flow past a cylinder is introduced, which is

*C*=

_{D}*F*/(

_{D}*ηU*).

It can be seen that the so-called tensile instability appears near the cylinder in the SPH results (Fig.

From Figs.

According to Ref. [17], the numerical convergence of SPH scheme not only depends on the particle numbers, but also relies on the smoothing length. Next, the simple numerical convergence is investigated by increasing the particle number *N _{y}* with different values of smoothing length

*h*.

Figure *N _{y}* at

*L*= 6,

_{r}*We*= 0.06, and

*L*represents the path over the cylinder centerline with the path being at 45° with respect to the horizontal line. It can be found that

_{xy}*U*,

*τ*′

^{xy}, and

Figure *Ro* = 2 with different values of smoothing length *h* = 2.0*d*_{0}, 2.4*d*_{0}, 2.8*d*_{0}, *L _{r}* = 6,

*We*= 0.8,

*β*

_{0}= 0.2. It can be found that the numerical results have good convergences at

*h*≥ 2.0

*d*

_{0}in our numerical simulations. Meanwhile, the tensile instability at

*h*= 2.4

*d*

_{0}is more serious than that at

*h*= 2.8

*d*

_{0}. The analysis has shown good convergence for the proposed method at

*h*≥ 2.4

*d*

_{0}.

Figures *L _{r}* = 6 reaches a stable value with increasing

*N*. The values of drag coefficient for the Newtonian fluid are nicely consistent with those in Refs. [31] and [32], the values of drag coefficient for the Oldroyd-B fluid are smaller than those for the Newtonian fluid. The change of the drag coefficient with

_{y}*L*= 6 is more complex than that with

_{r}*L*= 2.5 with increasing Weissenberg number (Fig.

_{r}In order to demonstrate the ability of the proposed particle method to capture the complex free surface in the true flow, the 3D viscoelastic jet bulking problem of filling process in a rectangle container is simulated by the SPH_CS_SP method with the MPI parallelization technique. In this section, the instability phenomenon of jet coiling based on the Oldroyd-B fluid is mainly simulated and predicted.

The geometry of 3D jet bulking is shown in Fig. *r*_{c} and the entry speed at the nozzle is *U*. The length, width and height of the rectangular container denote *L _{x}*,

*L*and

_{y}*L*, respectively.

_{z}*Re*= 2

*Uρr*

_{c}/

*η*is related to the viscosity of the fluid, and high value of

*Re*corresponds to fluid with low viscosity.

*We*=

*λ*

_{0b}

*U*/(2

*r*

_{c}) is related to the viscoelasticity of the fluid.

Figure *w*-velocity contour, and the jet coiling phenomenon of the Newtonian fluid is also shown for comparison. The parameters are as follows: the jet width *D* = 5 mm(*D* = 2*r*_{c}), *L _{x}* =

*L*= 5 cm,

_{y}*L*= 6 cm,

_{z}*U*= − 1 m·s

^{−1},

*L*/

_{z}*D*= 12, the reference density

*ρ*

_{0}= 1.1 × 10

^{3}k·gm

^{−3}, the gravitational acceleration

*g*= − 9.81 m·s

_{z}^{−2}, kinematic viscosity

*υ*

_{0}= 0.02 m

^{2}·s

^{−1}, speed of sound

*c*= 8 m·s

^{−1},

*β*

_{0}= 0.4,

*λ*

_{0b}= 0.005 s (

*We*= 1),

*λ*

_{0b}= 0.5 s (

*We*= 100),

*Re*= 0.25, and the time-step Δ

*t*= 10

^{−6}s.

According to previous research,^{[34,35]} the phenomenon of jet coiling based on the Newtonian fluid will occur, if *Re* < 0.5 and *L _{z}*/

*D*> 10 or

*Re*< 1.2 and

*L*/

_{z}*D*> 7 are satisfied. From Fig.

*Re*= 0.25 and

*L*/

_{z}*D*= 12, the reason can be explained as follows. On the one hand, the effect of the viscous force on flow is much bigger than that of inertia force after the fluid has touched the solid plate, thereby leading to the slow flow of the fluid. On the other hand, the fluid is continuously injected from the nozzle. The results of instability phenomenon of jet coiling based on the Oldroyd-B fluid basically are consistent with those obtained by FDM with MAC method in Refs. [30] and [35]. Also, the absolute value of

*w*-velocity near the front of jet fluid is bigger than that of the inlet velocity before the jet fluid reaches the bottom of rectangle container (

*t*= 0.027 s). Thus the required time for Oldroyd-B fluid to arrive at the bottom of the rectangular container is shorter than that for Newtonian fluid, and the jet coiling phenomenon of the Oldroyd-B fluid lags behind that of Newtonian fluid due to the shear thinning behavior of the Oldroyd-B fluid, which reduces the viscosity when subjected to shear strain and is the unique characteristic of the non-Newtonian fluid. The

*w*-velocity of the front of the jet fluid increases then decreases with time going by. Comparing Fig.

In order to clearly show the instability phenomenon of jet coiling, figure *We* = 100. It can be seen that the jet coiling phenomenon indeed occurs which is more obvious than that in Fig.

Next, the influences of macroscopic rheological parameters on the fluid deformation of filling process are predicted and discussed to further illustrate the complex phenomenon in the 3D viscoelastic filling process.

*We*

Observing Fig. *We* has an important effect on the deformation process of viscoelastic jet fluid, and the phenomenon of jet coiling may occur for higher *We* at lower Reynolds number (*We*=100 in Fig.

To further demonstrate the influence of *We* on the deformation process of jet coiling, we investigate the fluid deformations and the jet lengths of the jet buckling for the Oldroyd–B fluid with different Weissenberg numbers at *Re*=0.6. Here the jet length is defined as the distance from the jet front to the nozzle. The employed parameters are as follows: *D* = 6 mm, *L _{x}* =

*L*= 6 cm,

_{y}*L*= 5 cm (

_{z}*L*/

_{z}*D*= 8.3),

*U*= −0.5 m·s

^{−1}, Δ

*t*= 5 × 10

^{−6}s,

*υ*

_{0}=

*η*/

*ρ*

_{0}= 0.005 m

^{2}·s

^{−1}(

*Re*= 0.6),

*λ*

_{0b}= 0.006 s (

*We*= 1),

*λ*

_{0b}= 0.03s (

*We*= 5),

*λ*

_{0b}= 0.6s (

*We*= 100), and the other parameters are the same as those in Fig.

From Fig. *We*=1 at *t* = 0.5s, but it does not occur when *We* changes from 1 to 100 due to the effect of elastic stress and shear thinning. The reason is that moreover, the forefront of the jet fluid is warped as the jet coiling happens, which may be obviously observed from Fig. *t* = 0.09 s). The warped distortion of jet fluid may be another reason for generating the instability phenomenon of jet coiling, which was not mentioned in previous studies.^{[30,34,35]}

Figure *We* indicates that the relaxation time of the fluid increases, now the shear thinning character of the fluid instead of the elasticity dominates the flow in a short time before the fluid retracts. The conclusion is consistent with the remarks in Refs. [30] and [35].

*Re*

Another macroscopic parameter which has an important influence on the deformation of jet fluid is the Reynolds number. Comparing Fig.

In order to further show the effect of the Reynolds number on the viscoelastic flow, figure *Re*=1, 5, and 10. The parameters are as follows: *D* = 8 mm, *L _{x}* =

*L*= 8 cm,

_{y}*L*= 6 cm (

_{z}*L*/

_{z}*D*= 7.5),

*U*= –0.5 m·s

^{−1},

*λ*

_{0b}= 0.006 s (

*We*= 1),

*υ*

_{0}= 0.004 m

^{2}·s

^{−1}(

*Re*= 1),

*υ*

_{0}= 0.0008 m

^{2}·s

^{−1}(

*Re*= 5),

*υ*

_{0}= 0.0004 m

^{2}·s

^{−1}(

*Re*= 10), and the other parameters are the same as those in Fig.

From Fig. *Re*=1 (Figs. *Re* ≥ 1. As the Reynolds number increases up to 5 or 10, the concave phenomenon of forefront fluid appears as the jet fluid impacts on the container bottom (Figs.

In this work, the flow behavior of 3D jet buckling problem based on the Oldroyd-B model is simulated and predicted by a corrected particle scheme (SPH_CS_SP) method. The accuracy of the SPH_CS_SP method and the necessity of the shifting particle technique are first verified by two benchmark problems, that is, the Poiseuille flow and Taylor–Green flow, respectively. Then, the validity of the proposed SPH_CS_SP method to simulate the viscoelastic free surface is shown by the polymeric fluid flow around a periodic array of cylinders. Finally, the present method with MPI parallelization technique is further used to simulate the 3D jet coiling problem of filling process. The influences of macroscopic parameters on the instability phenomena are also investigated. Some conclusions can be drawn from the numerical results below.

**Reference**