numerical framework

Polymer xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Polymer journal homepage: www.elsevier.com/locate/polymer

Hydroentanglement of polymer nonwovens. 1: Experimental and theoretical/numerical framework Gen Lia, Christopher Staszela, Alexander L. Yarina,∗, Behnam Pourdeyhimib a b

Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, 842 W. Taylor St., Chicago, IL, 60607-7022, United States 3427 The Nonwovens Institute, Box 8301, North Carolina State University, Raleigh, NC, 27695-8301, United States

HIGHLIGHTS

model of hydroentanglement of polymer fibers was proposed. • AAntheoretical analytical solution of filtration flow inside the nonwovens mat was incorporated. • Numerical simulations of polymer fibers evolution with and without fiber-fiber and self-interaction were performed. • A model hydroentanglement experiment developed. • ABSTRACT

Hydroentanglement is a versatile and important process used to form highly entangled nonwoven materials comprised of polymer fibers. It is a key element of polymer processing in the world, with the nonwovens market worth of tens of billions of the US dollars. No fundamental theoretical models of hydroentanglement were available so far, even though such modeling is required to facilitate the entirely empirical efforts to optimize the process. No model experiments on hydroentanglement aiming at its underlying physics were conducted either. These challenging problems are in focus in the present work. A model experiment is conducted and a quasi-one-dimensional model of individual polymer fibers is proposed and implemented for a number of fibers subjected to water jets impacting as a curtain normally to the nonwoven surface and undergoing filtration motion in the inter-fiber pores. The nonwoven is assumed to be located on a moving substrate with water suction through it. The model allows for a number of two-dimensional water jets and a number of suction ports in the substrate. The quasi-one-dimensional model allows for the three-dimensional motion of individual viscoelastic polymer fibers and the fiber-fiber interaction, which makes them non-self-intersecting. Then, the governing quasi-one-dimensional equations of fiber motion are discretized and solved numerically by using a discrete element method with mesh refinement of the Lagrangian mesh during the time marching. Using the developed code, several basic model problems were solved. Specifically, the behavior of several viscoelastic polymer fibers under the action of water jets with and without fiber-fiber interaction is described. Accordingly, a quantitative measure of the fiber-fiber entanglement is introduced based on the observed morphologies in the accompanying second part of this work.

1. Introduction With the development of the nonwoven technology, the traditional materials have been replaced in numerous industrial fields by novel materials with better performance, higher production rates and lower costs. For example, The National Aeronautics and Space Administration (NASA) has developed a program to use nonwoven polyester fabrics as reinforcement in polyethylene terephthalate (PET) films [1]. Hydroentanglement process is one of the key elements in production of nonwoven materials. Hydroentanglement employs high velocity water jets to entangle and bond the loose polymer fibers in a nonwoven laydown. The subsequent mechanical performance and physical integrity of hydroentangled fibrous nonwovens strongly relies on the fiber-tofiber friction [2], since the friction forces between the entangled fibers

∗

are enormous due to the high contact areas in the knots. The thesis of Xiang [3] describes the production process in which loosely entangled meltblown polymer fibers form an integrated nonwoven web being subjected to hydroentanglement. It should be emphasized that most of the research works dealing with hydroentanglement are purely empirical and a comprehensive theoretical model which will account for the fiber-fiber interaction is unavailable [3,4]. A 2-D model has been developed in Ref. [5] to account for the effect of the specific energy of water jets and the corresponding drag force imposed on the fibers. The disadvantages of this model stem from the fact that it could not correlate the production rate with the data on the resulting nonwoven thickness, as reported in Refs. [3,4]. A 3-D model is required to correctly predict fiber behavior. One of the popular experimental methods used to reconstruct the 3-D microstructure of

Corresponding author. E-mail address: [email protected] (A.L. Yarin).

https://doi.org/10.1016/j.polymer.2018.11.059 Received 14 September 2018; Received in revised form 28 October 2018; Accepted 3 November 2018 0032-3861/ © 2018 Elsevier Ltd. All rights reserved.

Please cite this article as: Li, G., Polymer, https://doi.org/10.1016/j.polymer.2018.11.059

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

hydroentangled nonwovens is the Digital Volumetric Imaging (DVI). DVI results were used to predict the nonwoven permeability, as well as to visualize the internal structure [6,7]. Recently, DVI was used to analyze the internal fiber-fiber entanglement and arrangement in twolayer nonwovens of different colors [8]. The images revealed that the fibers from the top layer penetrate into the bottom layer, which elucidates fiber transport from the surface to the bulk and intermingling of the fibers. Another method which is widely used to investigate polymer nonwovens is the Scanning Electron Microscopy (SEM). Most of the studies employing SEM focus on the top layer of nonwovens (the surface fibers), which allows one to evaluate the in-plane probability density function of fiber orientation [9,10]. However, to investigate the internal structure of nonwovens, the cross-sectional images are required. Such results for hydroentangled nonwovens are scarce, because it is difficult to achieve such imaging without breaking the hydroentangled structure, and no methods which can reveal the evolution of nonwoven structure in time during hydroentanglement are available. The idea of considering water jets penetrating through nonwovens during hydroentanglement was discussed in Ref. [11]. Solid polymer materials, including fibers, are viscoelastic [12], which should be accounted for in modeling of their behavior during hydroentanglement. The governing equations of the dynamics of individual viscoelastic fibers in nowovens are basically the same as in the dynamics of free liquid polymer jets moving in air. These are the quasi-one-dimensional equations introduced in Refs. [13,14]. In the quasi-one-dimensional approximation, the fibers can be considered as a slender body, which can be characterized by the cross-section radius distribution. The quasione-dimensional equations are widely used in the numerical works on free liquid jets, electrospinning, as well as the solution- and meltblowing. The numerical prediction of the formation of the beads-on-astring structure in capillary viscoelastic jets can be found in Refs. [13,15]. The numerical modeling of the bending instability in electrospinning based on the quasi-one-dimensional equations in the momentless approximation was done for the first time in Ref. [16] (also cf. [17,18]). Modeling of meltblowing, also based on these equations, was done in Refs. [19–21]. In all these works the quasi-one-dimensional equations were supplemented with a viscoelastic rheological constitutive equation. The present work aims at the model experiment elucidating the physical aspects of hydroentanglement, and at the development of a comprehensive model of hydroentanglement of viscoelastic polymer nonwovens based on the quasi-one-dimensional equations for the individual fibers and description of filtrating water jets in the framework of the complex analysis. This model ramifies into a robust computational approach to hydroentanglement, which allows for prediction of fibers evolution and the fiber-fiber interactions and entanglement in three dimensions in the framework of an optimized discrete element Lagrangian method. A quantitative measure of the fiber-fiber entanglement is also introduced and implemented and is discussed in detail in the second part of this paper [36]. The experimental setup is outlined in section 2. The theoretical model of fiber dynamics and its numerical implementation are described in section 3. In particular, the filtration water flow velocity fields generated by the impinging water jets are discussed in sub-section 3.1. The quasi-one-dimensional model of viscoelastic fibers is detailed in sub-section 3.2. The preliminary results are presented in section 4. In addition, an optimized numerical method is introduced and its results are presented in Appendix B. Conclusions are drawn in section 5.

unit (Sunjoe SPX3000) and a suction unit (Shop-Vac 2036000). The encasing chamber was used to house the inlet water jet issued from the pressure unit and the water outlet with vacuum suction directly below (cf. Fig. 1). The nonwoven samples were loaded onto a highly porous metal mesh fixture with a supporting base at a distance of 15 mm from the inlet water jet. In distinction from any pilot experimental platform, the present laboratory setup incorporated only a single water jet issued at the velocity of about 100 m/s from a 1 mm diameter orifice. That means instead of forming a fully hydroentangled nonwoven, only localized hydroentangled lines were aimed here, which is quite sufficient to elucidate the underlying phenomena. In the experiments, hydroentanglement was conducted using two separate polyethylene terephthalate (PET) nonwovens. The first nonwoven sample was a PET nonwoven with the basis weight of 45 GSM (grams per square meter) and the thickness of approximately 4.0 mm. The fiber diameter in this nonwoven was 12 μm. From here on this nonwoven will be referred to as the fluffy nonwoven. This nonwoven was white. The second nonwoven which was also a PET nonwoven was pre-dyed black for visualization purposes. This nonwoven had the basis weight of 265 GSM and the thickness of approximately 3.0 mm. This nonwoven was comprised of fibers with the diameter of 40 μm and was significantly denser then the fluffy nonwoven. Accordingly, the latter nonwoven is referred to as the dense nonwoven. The fluffy nonwoven was placed on top of the dense nonwoven. That means the fluffy nonwoven was facing the water jet, and was pushed by that toward the underlying dense nonwoven in the experimental fixture (cf. Fig. 1). This sandwich-like structure was then prewetted by soaking it with water, while being loaded in the experimental setup with the pressure unit turned off. After that, the pressure unit was turned on and the sandwich-like structure was moved perpendicular to the water jet. As a result of the water jet pressure, the fluffy nonwoven was hydroentangled with the underlying dense nonwoven along the line traced by the water jet. After hydroentanglement, the samples were removed and dried in air. Fig. 2 shows an example of the hydroentangled nonwovens stitched to one another by the impinged water jet tracing a line. The results shown in Fig. 2 reveal that the fibers from the white upper nonwoven (initially at the top) penetrated through the underlying dense black nonwoven due to the traction imposed by the water jet. As a result, the white fibers were hydroentangled with the underlying black nonwoven. It should be emphasized that even though the dense nonwoven was significantly less permeable, the 12 μm and the 40 μm fibers from both nonwovens were entrained by water filtration through the nonwovens and mutually hydroentangled. 3. Theoretical model and its numerical implementation 3.1. Water flow field created by the impinging jets To model the dynamics of polymer fibers in hydroentanglement process, one needs to consider the hydrodynamic drag force acting on the fibers in nonwovens from the impinging and penetrating water jets. A nonwoven mat is delivered by a moving screen which has a velocity Vs and it takes time ts=(b-a)/Vs for the nonwoven material elements to pass under the jet, where a and b are the coordinates of the jet boundaries as shown in Fig. 3b. On the other hand, the characteristic time of jet penetration into nonwoven is tp = h0/V, where h0 is the initial mat thickness and V is the jet impingement velocity. It should be VS . For the realistic dimensions of the water jets emphasized that V and the nonwoven mat thicknesses resembling those in the industry [11,22], the jet diameter (b-a) = 0.127 mm and the initial nonwoven thickness h0 = 0.711 mm. Additionally, in the case of two water jets, the distance between their centerlines is taken as 0.635 mm. The velocity of the water jets and the moving screen can be taken as V = 100 m/s and Vs = 0.75 m/s, respectively. Then, the characteristic time of the water jet penetration into nonwovens tp and the time ts

2. Experimental setup In order to elucidate the underlying physical aspects of hydroentanglement and facilitate the accompanying numerical modeling, a lab-scale hydroentanglement experimental setup shown in Fig. 1 was designed and built. The setup incorporated an encasing chamber, an electric pressure 2

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 1. Schematic of the lab-scale experimental setup.

required for a nonwoven material element to pass under the jet are evaluated as tp = 7.1×10 6 s and t s = 1.7×10 4 s , which shows that tp t s . In other words, the water penetration process can be treated as a quasi-stationary one irrespective of the screen motion. Also, since the water-jet curtain in hydroentanglement is as wide as the substrate belt, i.e. it is much wider than the nonwoven mat thickness, the water jets can be treated approximately as a two-dimensional flow. The two-dimensional velocity fields created by water penetration into a porous nonwoven mat can be described by Darcy's law [23–25].

V=

k p µ

potential is an analytic function of the complex variable z = x+iy, with the Cartesian coordinates x and y being shown in Fig. 3, and i being the imaginary unit, i.e. = (z ) , and thus the conjugate velocity V (z ) is found as

V (z ) =

+i

iv

(3)

Given the boundary conditions of the problem imposed on the ycomponent of velocity v, cf. Fig. 3 (where the velocity and coordinates are rendered dimensionless), it is convenient to use the following analytic function

(1)

(4)

m ( z ) = i V ( z ) = v+ i u

where p is pressure, k is the mat permeability, and µ is the viscosity of water. Equation (1) shows that the filtration velocity V is generated by the potential = (k /µ ) p . Due to the continuity equation, this potential inevitably satisfies the Laplace equation, 2 = 0 . To solve the Laplace equation for the potential, as usual in twodimensional cases, a complex potential , comprised of the scalar poand the function related to via the Cauchy-Riemann tential conditions [26] is introduced as

=

d =u dz

rather than V (z ) . Then, since both v and u are harmonic functions, m(z) is found using the Palatini formula [27,28].

m (z ) =

1 2

v

coth y=0

(t

z) 2

dt +

i 2

v

tanh y=1

(t

z) 2

dt

(5)

where the t is a real dummy variable reckoned along the x-axis. In the particular case of two water jets, with the boundary condition imposed as in Fig. 3, Eq. (5) yields

(2)

Due to the Cauchy-Riemann conditions, the above complex

Fig. 2. Images of the hydroentangled fluffy nonwoven (white) on the dense nonwoven (black). The water jet traced two lines. Accordingly, some white fibers, which were initially at the top became visible on the bottom side of the underlying black material. 3

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 3. Schematic of water penetration into nonwoven mat: (a) a single water jet, (b) two water jets. In the sketch the jets impinge from below and suction is from above the mat.

+ i

=

a

i 2

m (z ) =

1 coth

b

i 2

c

{ln

v1 tanh

c

(1 z ) dt 2 (t

z) 2

+

b

1 coth

a

(t

z) 2

Newton (the momentum balance) should be used with all the forces involved included. The momentum balance should be formulated for an individual infinitesimally small element of an individual fiber in the quasi-one-dimensional approximation [13,20]. The relevant forces imposed on an infinitesimally small element of an individual fiber are comprised of the inertial forces, the internal forces in the fiber crosssection of the viscoelastic origin, and the external force associated with hydrodynamic drag imposed by filtrating water brought by the water jets. It should be emphasized that forces associated with the fiber-fiber interaction will be added later on. The quasi-one-dimensional equations describing fiber evolution are partial differential equations [13,20]. However, these equations can be formulated from scratch in a discretized form convenient for numerical simulations [16] using a set of discrete grid elements denoted by i (different from the imaginary unit in section 3.1). The fiber element which connects nodes i and i + 1 is denoted by an additional subscript u and the one which connects nodes i and i 1 denoted by an additional subscript d as schematically shown in Fig. 5. Thus, the lengths of these elements expressed using the Cartesian coordinates X, Y and Z of the nodes are

dt

dt

sinh [ (a + z ) / 2] sinh [ (b sinh [ (b + z ) / 2] sinh [ (a

z ) / 2] z ) / 2]

}+

i v1

ln

{

cosh [ (c z ) / 2] cosh [ (c + z ) / 2]

}

(6)

where the jets inflow the mat with the dimensionless vertical component of velocity being 1, and suction through the upper surface of the mat in Fig. 3b being given as the dimensionless v = v1. The water mass balance inside the nonwoven mat requires the injection being equal to the suction, and thus

2(b

a)v0 = 2c v1

(7)

which yields

v1 =

b

a c

v0

(8)

The velocity field predicted by Eqs. (6) and (8) is depicted in Fig. 4c and d. Also, Fig. 4a and b shows the velocity field predicted for a single water jet sketched in Fig. 3a. 3.2. Modeling of viscoelastic fibers in nonwovens in hydroentanglement 3.2.1. Description of a single fiber without fiber-fiber interactions Polymer fibers are known to be viscoelastic [12]. To model the behavior of viscoelastic fibers in hydroentanglement, the second law of

ui

=

(Xi + 1

di

=

(Xi

Xi )2 + (Yi+ 1

Xi 1)2 + (Yi

Yi )2 + (Zi + 1

Yi 1)2 + (Zi

Zi )2

(9)

Zi 1)2

(10)

It should be emphasized that the coordinates Xi, Yi and Zi are the Lagrangian coordinates of a material element i. Using the rheological constitutive equation of a viscoelastic material [12,29–32] the internal

Fig. 4. Velocity components v [panel (a)] and u [panel (b)] in the cases of a single jet. Water flow field corresponding to penetration of two water jets into a nonwoven mat is presented in panel (c) where the y-component of velocity (v) is shown, and in panel (d), where the x-component of velocity (u) is depicted. 4

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

balance for an individual material node is expressed by Newton's second law as

m

stresses σui and σdi in these elements can be found from the following equations

G µ

ui

(11)

d di 1 d di =G dt di dt

G µ

di

(12)

i: m

j: m

CD V 2

Vf )2a (

ui

+

di)

Xi Y Z +j i +k i t t t

a + adi ai = ui 2

(Ri + 1

adi2

Ri )

di

di

(Ri

Ri 1)

(17)

) + (v

Xi 2 t

i

Xi )

) + (0

Yi 2 t

i

2 adi di (Xi di

)

Zi 2 ai ( di t

+

ui )

(u

i

Xi t

)

Xi 1)

(u

) + (v

Xi 2 t

i

2 aui ui (Yi + 1 ui

) + (0

Yi 2 t

i

2 adi di (Yi di

Yi )

)

Zi 2 ai ( di t

+

ui)

(v

i

Yi t

)

Yi 1)

(19) d2Zi dt 2

= CD

+

(u

i

2 aui ui (Zi + 1 ui

) + (v

Xi 2 t

Zi )

i

) + (0

Yi 2 t

2 adi di (Zi di

)

Zi 2 ai ( di t

+

ui)

(0

Zi t

)

Zi 1)

(20) It should be emphasized that the velocity components of the surrounding water flow in Eqs. (18) - (20), as discussed in sub-section 3.1, involve zero component in the k -direction, since the filtration flow field is assumed to be two-dimensional. Also, denote the velocity components of the material fiber element as

(14)

where V is the absolute local velocity vector of water calculated in section 3.1, Vf is the absolute local velocity of the fiber, and a is the cross-sectional radius of the fiber. In the framework of the adopted discretization

Vf , i = i

(u

2 aui ui (Xi + 1 ui

= CD +

(13)

Vf (V

= CD

d2Yi dt 2

k: m

Fd V2R A/2

ui

ui

(18)

where ρ is water density, Fd is the drag force, VR is local water velocity relative to the fiber element and A is the median area. For the fiber element i , the drag force can be expressed as

Fd =

d2Xi dt 2

+

where t is time, G and µ are the elastic modulus and dynamic viscosity, respectively. The drag force imposed on fiber elements by the surrounding water is found employing the drag coefficient CD

CD =

aui2

where the Ri is the position vector of the element. Equation (17) expresses the balance of the inertial force acting on a fiber element with the internal viscoelastic forces and the hydrodynamic drag. Substituting Eqs. (14) and (15) into Eq. (17), and projecting the latter onto the Cartesian axes, one obtains the following three equations

Fig. 5. Schematic of polymer fiber discretized as a series of material elements serving as nodes in the Lagrangian algorithm.

d ui 1 d ui =G dt ui dt

d 2Ri = Fd + dt 2

(15) (16)

where i, j and k are the unit vectors of the Cartesian axes X, Y and Z, respectively. The subscripts ui and di mean up from node i (i.e. the bond connecting nodes i and i+1), or down from node i (i.e. the bond connecting nodes i-1 and i). If the fiber-fiber interactions are not accounted for, the momentum

Xi = Vxi t

(21)

Yi = Vyi t

(22)

Zi = Vzi t

(23)

In summary, the overall modeling algorithm is shown schematically in Fig. 6 where TD, CD, MD denote the thickness-direction, the crossdirection and the machine-direction, respectively.

Fig. 6. Schematic of the overall modeling algorithm for hydroentanglement.

5

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

3.2.2. Modeling of the fiber-fiber interactions The Lennard-Jones (L-J) potential [33] VLJ is used in the present case to model the fiber-fiber interactions 12

VLJ = 4

6

r

(24)

r

where the is the depth of the potential wall, is the distance where the potential is zero, and r is the distance between two material elements which might either belong to two different fibers or to the approaching elements of the same fiber during self-entanglement. The first term on the right-hand side in Eq. (24) with the higher order corresponds to the repulsive force and the second one - to the attractive force. In the fiber-fiber interactions the most important role will be played by the repulsion at short distances, while the attraction at longer distances will be typically negligibly small compared to the other forces involved. Accordingly, only the term for the repulsive force will be accounted for in the projections of the momentum balance equation (18) - (20). From the L-J potential field, the repulsion force between two fibers is found as

FLJ =

VLJ ds =

s

s

12

4

r

ds

Fig. 7. Schematic of two fibers interacting with each other. Blue arrows show the local hydrodynamic traction acting on the fibers. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

4.1. Prediction of fiber evolution under the effect of water jets without the fiber-fiber interaction First of all, the dimensions of the water jets and their velocity need to be specified. In the case of two water jets shown in Fig. 3b which is used as a generic case in this work, the jet and suction dimensions were set as

(25)

where ds is the length differential along a fiber. In the discrete-node formulation discussed in sub-section 3.2.1, the projections of the L-J force between two fibers with elements denoted as n and m are expressed via (Xm Xn ) , etc. as N

i: Fni = 48

12

(Xm

Xn)( (Xm + 1 ( (Xm

m=0

Xm) 2 + (Ym + 1

Xn) 2 + (Ym

Ym)2 + (Zm + 1

a = 0.08, b = 0.12, c = 0.1, h 0 = 0.25

It should be emphasized that the transformation factor for the length scale (the initial mat thickness h 0 used in Fig. 3 and the Palatini formula in section 3) is now needed for the initial fiber element length L0 used as the length scale in the calculations of the fiber dynamics here. The simulation results discussed here begin with the example of two water jets and a single fiber located initially as a straight line along the Y-axis. It should be emphasized that the X, Y, Z-axes in the following part of the work differ from the coordinate system used in section 3. The x-axis in section 3 corresponds to the Y-axis here and the y-axis in section 3 corresponds to the Z-axis here. Fig. 8 depicts the results of the simulations for this case. The predicted fiber positions corresponding to time moments from t = 0 to t = 0.01 with the interval of t = 0.002 are shown. The fiber is stretched by the drag forces imposed by the powerful water jets and moves to the top (rear) side of the mat. Two bumps in its shape appear due to the action of the two water jets. The following dimensionless parameters were used: B = 4 , Fve = 10 , Q = 0 (cf. Appendix B) and the number of the material elements N = 100. The initial positions of all the material elements in the fiber were set as

Zm)2 )

14

Yn) 2 + (Zm

Zn) 2 )

(26) N

j: Fnj = 48

12

(Ym

Yn )( (Xm + 1 ( (Xm

m=0

Xm

)2

+ (Ym + 1

Xn )2 + (Ym

Ym

)2

+ (Zm + 1

Yn )2 + (Zm

Zn )2 )

Zm

)2 )

14

(27) N

k: Fnk = 48

12 m=0

(Zm

Zn )( (Xm + 1 ( (Xm

Xm )2 + (Ym + 1

Xn )2 + (Ym

Ym)2 + (Zm + 1

Yn)2 + (Zm

Zm)2 )

14

Zn)2 )

(28) where primes denote elements of ‘another’ fiber approaching a given one. Schematic of two fibers interacting with each other is shown in Fig. 7. In the case of self-looping also discussed below, primes denote ‘another’ element (n) of the same fiber approaching a given one (m), and the summation in Eqs. (26)-(28) excludes m = n, i.e., the summation in Eqs. (26)-(28) is done for m n . Equations (26)-(28) describe the pseudo-L-J force associated with a fiber which has N material nodes. These equations are added to the right-hand sides of Eqs. (18)-(20) to account for the fiber-fiber interactions.

Y =

The dimensionless form of the equations is given in Appendix A. The equations (18) - (23) are solved numerically to predict the fiber evolution during hydroentanglement. The initial conditions at t = 0 are imposed as (29)

¯i = 0 V

(30)

=0

(31)

i

0.25 + 0.005i, X = 0, Z = 0.1

(33)

where i is the element (node) number, 0 i N . Fig. 9 depicts the simulation results for the cases with different initial distances between the fiber and the two water jets. Comparing with the results depicted in Fig. 8, the fiber shapes depicted in Fig. 9a, b and c reveal the same tendency to develop two bumps originating from the two water jets from t = 0 to t = 0.01 irrespective of the initial position of the fibers relative to the jets. It should be emphasized that the fibers will be stuck at the substrate surface h 0 = 0.25 at the latter time moments if the simulations would be continued. In Figs. 8 and 9a, the bumped shapes are still not visible close to the beginning of the process at t = 0.002 and the depth of the valley between the two bumps is still very small. When the initial distance of a fiber from the water jet origin decreases, the earlier bump formation is facilitated as the results showed in Fig. 9b and c reveal. The fiber self-looping between the two water jets illustrated in Fig. 10 can be one of the key phenomena characteristic of hydroentanglement. It is driven by the structure of the water flow field between the two jets, with the horizontal flow components impinging onto each other at the borderline where the stagnation zone emerges.

4. Results and discussion

= 0.5

(32)

while the initial configurations of fibers expressed through the coordinates of their elements are discussed below in each case separately. In the numerical implementation the fourth-order Runge-Kutta method was adopted for time marching [34]. 6

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 10. The evolution of the two bumps on the fiber under the action of two water jets. (Red-t = 0 , Green-t = 0.002 , Blue-t = 0.004 , Cyan-blue-t = 0.006, Yellow-t = 0.008, Pink-t = 0.01). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 8. Evolution of a single fiber under the action of two water jets from t = 0 to t = 0.01 (Red-t = 0 , Green-t = 0.002 , Blue-t = 0.004 , Cyan-blue-t = 0.006, Yellow-t = 0.008, Pink-t = 0.01). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

This two-bump structure is very similar to the so-called ‘U’-shape fibers which are bent due to the water impact. Such ‘U’-shape fibers were observed by SEM of hydroentanged fabrics, as reported in Ref. [37]. The initial fiber shape can be also set to be non-straight, for example being harmonically perturbed (with 1 i N )

X =

0.25 + 0.005i,

Y = 0.005 sin(32 X ),

Z = 0.025

(34)

Fig. 11 shows that the water flow structure dictated by the two water jets is the dominant factor, and the fiber rapidly ‘forgets’ its initial harmonically perturbed shape and develops the characteristic doublebump configuration similar to that of the initially straight fibers in Figs. 9 and 10. Another possible initial fiber configuration could be the following one (with 1 i N )

Fig. 11. Fiber evolution starting from a harmonically perturbed case. (Redt = 0 , Green-t = 0.002 , Blue-t = 0.004 , Cyan-blue-t = 0.006, Yellow-t = 0.008, Pink-t = 0.01). It should be emphasized that there is no self-intersection of the fiber here. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 9. Fiber evolution starting from four different initial locations shown in panels: (a) z = 0.1, (b) z = 0.05, (c) z = 0.025, (d) z = 0.0125, (Red-t = 0 , Greent = 0.002 , Blue-t = 0.004 , Cyan-blue-t = 0.006, Yellow-t = 0.008, Pink-t = 0.01). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 7

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 12. Fiber evolution beginning from the initial configuration given by Eq. (35). (Red-t = 0 , Green-t = 0.002 , Blue-t = 0.004 , Cyan-blue-t = 0.006, Yellow-t = 0.008, Pink-t = 0.01). It should be emphasized that there is no self-intersection of the fiber here. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Y =

0.25 + 0.005i, X = Y , Z = 0.025

(35)

Fig. 12 depicts the corresponding results on the fiber evolution. The formation of the double-bumped structure is still clearly visible, which still cannot result in a complete self-intersection due to the structure of the flow field imposed by the water jets. It should be emphasized that in the results shown in Figs. 8–12 and hereinafter the two ends of the fibers are not fixed in space but free to move. However, their actual mobility is very limited and slow compared to the elements located against the water jet inlets, because the fiber ends are located in the stagnation zones of the filtration water flow, and the hydrodynamic drag force acting on them is very small. 4.2. Prediction of fiber evolution under the effect of water jets accounting for the fiber-fiber interactions

Fig. 13. The initial crossbar configuration of two fibers.

By accounting for the pseudo-L-J repulsion as discussed in subsection 3.2.2, the simulation can incorporate the fiber-fiber interactions. In this sub-section, two fibers with the initial crossbar configuration (with 1 i N )

X1 = 0,

X2 =

Y1 =

0.25 + 0.005i,

0.25 + 0.005i,

Y2 = 0,

Z1 = 0.025

(36)

Z2 = 0.04

(37)

are considered The fiber initially located further from the water jets origin is shown in green in Fig. 13. It is set to be fixed and ‘frozen’ serving just as an obstacle for the closer fiber shown in red, which can evolve. Fig. 14 shows the result of the simulation at t = 0.01. The following dimensionless group values were used: B = 4 , Fve = 10 , Q = 5 × 10 7 , and the number of the nodes was N = 100 ; also the time step t¯ = 5 × 10 6 . To illustrate the effect of the fiber-fiber pseudo-L-J repulsion, a comparison of the present results with those from sub-section 4.1 is shown in Fig. 15. The fiber positions are depicted at t = 0.01 and the blue one is the fiber predicted without accounting for the fiber-fiber repulsion. It is seen that the blue line moved through the obstacle fiber shown in green, whereas the red line corresponding to the fiber with repulsion accounted for, formed an entanglement with the obstacte, never crossing it. Fig. 15 also shows that some elements of the fibers can be significantly stretched, as for example those of the fiber shown in red in the right-hand side panel. These elements look like straight lines, because the distances between the nodes of the Lagrangian grid used here become too large. This problem will be eliminated by using a dynamic mesh refinement, as well as by increasing the numbers of nodes, as discussed in Appendix B.

Fig. 14. Red fiber impinging onto a ‘frozen’ green obstacle at t = 0.01. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

In addition, the pseudo-L-J potential is also accounted for to model the self-repulsion of approaching elements of the same fiber. A distance threshold is introduced for different elements to be considered as forming self-looping. Specifically, for each of the discrete elements, the L-J potential-related repulsion from the nearby elements is neglected and only the repulsion resulting from the further located elements is accounted for to achieve non-self-intersection when loops (self-entanglement) are formed. For that aim, in the case of N = 500 , the L-J potential-generated repulsive forces from the nearby thirty elements are excluded. To illustrate the effect of the self-repulsion, the corresponding fiber configurations are depicted at t = 0.016 in Fig. 16. The following

8

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 15. The difference between the cases without and with fiber-fiber interaction at t = 0.01.

Fig. 16. The difference between the cases without and with fiber self-repulsion at t = 0.016. (a) Self-repulsion is not accounted for. (b) Self-repulsion is accounted for.

values of the dimensionless groups were used: B = 4 , Fve = 20 , Q = 5 × 10 7 and the time step t¯ = 1 × 10 6 . In Fig. 16a the fiber shown in red self-looped itself which would be unphysical for real material fibers in a 2-D velocity field as in Fig. 4. By adding the self-repulsion effect as in Fig. 16b, one can see that when distant elements from the same fiber closely approach each other, they repel each other and prevent the unphysical self-looping. In a real physical situation, the ‘frozen’ green fiber of Figs. 13–15 would also start moving during the interactions with another fiber. This situation is explored in Fig. 17, which shows the fiber-fiber intereaction at t = 0.017 when both fibers of Fig. 13 are set free to move. The following values of the dimensionless groups were used: B = 4 , Fve = 20 , Q = 5 × 10 7 and the time step t = 1 × 10 6 . The number of the elements was set as N = 500 . As shown in Fig. 17, when the green fiber is set free, due to the L-J

potential-generated repulsion from the lower red fiber, it will move up in the stagnation area of the water jets. Accordingly, when the two bumps in the red fiber form due to the drag force imposed by the water jets, a depression appears in the interaction area on the green fiber. Moreover, for the red fiber the zoomed-in image in Fig. 17 shows that the self-repulsion forms a complicated loop around the intereaction area in concert with the pseudo-L-J repulsion from the green fiber. In several cases considered above, a significant stretching of the Lagrangian mesh used to describe fibers was observed. Under such circumstances, the numerical results could become unreliable. To resolve this problem, a special dynamic mesh refinement is adopted to control the distances between the nodes (the material elements) following [18]. The detailed process of the dynamic mesh refinement is described in Appendix B.

Fig. 17. The interaction of two ‘free’ fibers accounting for the fiber-fiber repulsion and the self-repulsion. The snapshot at t = 0.017 .

9

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 18. SEM images of the hydroentangled lines. The nonwovens were cut without freezing in liquid nitrogen. Panels (a) to (d) correspond to the magnifications of × 35, × 100, × 300, × 300 (a different location), respectively.

orientated due to its already tightly entangled structure (the yellow line in Fig. 21). The comparison of these experimental results to the predictions of the theory of the multi-fiber entanglement is discussed in detail in the second part of this work [36] where such simulations are conducted. It should be emphasized that the main reason for choosing the two different types of nonwovens in the model experiment is to visualize and distinguish the hydroentanglement of them. To achieve this goal, for a direct observation, we choose different colors, so the entangled strip is clearly seen on the back side of the other material and can be measured accurately. For the SEM imaging, two different clearly recognizable diameters are imperative to distinguish formation of the ‘U’ shaped fibers penetrating and entangling different layers. As for the nonwoven thicknesses, there is no specific reason for choosing them in the model experiment if they are sufficiently thin to be hydroentangled by the water jet employed. Their availability is also a factor.

4.3. Experimental results To observe the micro-structure of the hydroentangled lines, images of the fluffy nonwoven hydroentangled to the dense nonwoven were obtained using SEM Hitachi S-3000. To achieve that, the samples should be cut in liquid nitrogen to prevent local fiber smearing which appear without such freezing (see Fig. 18). As seen in Fig. 18, a cut without using liquid nitrogen will locally smear and disrupt the hydroentangled structure. On the other hand, when cooling in liquid nitrogen is used before cutting, the hydroentangled structure remains intact, as is seen in Figs. 19–22. In order to facilitate distinguishing fibers with different diameters, the original SEM images are processed and parts of the fibers of 12 μm and 40 μm are colored in blue and brown color, respectively, in Figs. 20b and 22b. Figs. 19 and 21 illustrate in detail the effect of hydroentanglement. Note that it is possible to distinguish the fibers from the fluffy nonwoven from those of the dense nonwoven based upon the fiber size. To reinforce the difference in SEM images, some fibers are colored as in Figs. 20b and 22b. Here in the dense nonwoven the fiber size was typically around 40 μm (colored in brown) and is clearly larger in Figs. 20 and 22 than that of 12 μm from the initially fluffier nonwoven (colored in blue). In the areas where water jet impinged and entangled the two nonwovens, an increase in the fiber density is observed (cf. Fig. 19). Here water jet entrained the overlaying 12 μm fibers and pulled them into the underlying denser nonwoven; cf. Figs. 19 and 21. This in turn increased the density of the nonwoven sample locally at the hydroentanglement location. Also, due to the impinging water jet, the fibers in the fluffy nonwoven were reoriented in the direction of the impinging water jet (the green lines in Fig. 21). On the other hand, the fibers in the underlying dense nonwoven mat were not significantly re-

5. Conclusion In this work, a new model was developed to describe and predict the hydroentanglement process. The proposed and implemented numerical model of hydroentanglement is effective and physically sound. Cases without and with the fiber-fiber interaction were explored. Several model problems with a single or two fibers subjected to two water jets were explored in detail. The predicted single fiber evolution revealed a clear tendency to form a two-bump structure under the entrainment caused by the two water jets. This is one of the key phenomena in hydroentanglement. When two fibers interact with each other, the entanglement was predicted, and more complicated local bumped structures revealed. It was also demonstrated how the self-repulsion affects

10

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 19. Samples comprised of two different hydroentangled nonwovens cut in liquid nitrogen and shown at different magnifications: (a) × 35, (b) × 50 , (c) × 100 , (d) × 300 . Here the same location encompassed by the blue circle is progressively magnified from panel (a) to (d). Note that the red arrows point at the 12 μm fibers, and 40 μm fibers are from the initially fluffy and dense nonwovens. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 20. Samples comprised of two different hydroentangled nonwovens cut in liquid nitrogen. (a) Original SEM image in Fig. 19c. (b) Processed SEM image, with a part of the 12 μm fibers and 40 μm fibers colored in blue and brown, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

the process of a single fiber self-looping. The present results for a single fiber and two fibers are fundamental for the realistic numerical modeling of multi-fiber networks in nonwovens tackled in the second part of this work [36].

To investigate peculiarities of the hydroentanglement process experimentally, a lab-scale hydroentanglement setup with a single water jet was designed and constructed, and several experiments were conducted. From the cross-section SEM images of the hydroentangled

11

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

Fig. 21. Jet-impacted area of the hydroentangled nonwoven sample. Here the jet impact direction is sketched by the blue arrows. The location encompassed by the blue circle in panel (a) is magnified in panel (b). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 22. Jet-impacted area of the hydroentangled nonwoven sample. (a) Original SEM image in Fig. 21a. (b) Processed SEM image, with a part of the 12 μm fibers and 40 μm fibers colored in blue and brown, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

samples it was possible to observe how fibers were entrained by the water jet. The experimental results will be further discussed in the second part of this work [36] in comparison with the numerical simulations of multi-fiber networks.

Acknowledgement This work is supported by the Nonwovens Institute, grant No. 17211.

Appendix A. Dimensionless parameters To solve the governing Eqs. (18) - (23) numerically, they are rendered dimensionless as follows. For all the lengths such as the element length ui and the material element positions, such as Xi , the initial length of the fiber L0 is used as the length scale. Time t is rendered dimensionless by the relaxation time µ / G , the stresses - by G , the velocity components - by L0 G /µ . It should be emphasized that the length and velocity scales in the dimensionless expressions for the velocity field in sub-section 3.1 are different from the scales used in the governing equations of the fiber dynamics. Therefore, the corresponding scale transformation is done and implemented in the solver. Taking the x-velocity component of the water jets as an example, this transformation is done as follows

u v0 µ v0 L0 G

u=

(A1)

¯0 = v0/v0 = 1, while where the first multiplier u/v0 is the dimensionless velocity found from Eq. (6) and based on the dimensionless water jet velocity V the scale transformation factor is v0µ / L0 G . The cross-section radius is rendered dimensionless by its initial value a0 . The volume conservation for a material element reads ai2

i

(A2)

= a02 L 0

Denote dimensionless parameters by overbars. Then, Eqs. (18)–(23) and (26) - (28) take the following dimensionless form d ¯ui dt¯

=

d ¯di dt¯

=

[(X¯ i + 1

¯ xi + 1 X¯ i )(V

¯ xi) + (Y¯i + 1 V

¯ yi + 1 Y¯i )(V

¯ yi) + (Z¯i + 1 V

¯ zi + 1 Z¯i )(V

¯ zi)] V

¯ 2ui

[(X¯ i

¯ xi X¯ i 1 )(V

¯ xi 1) + (Y¯i V

¯ yi Y¯i 1 )(V ¯ 2di

¯ yi 1) + (Z¯i V

¯ zi Z¯i 1 )(V

¯zi 1)] V

¯ui

(A3)

¯di

(A4)

12

Polymer xxx (xxxx) xxx–xxx

G. Li et al. d Vxi dt

= B(

+

ui

ui

+ Fve

ui + di

di)

ui

(Xi + 1

Xi )

+Q

= B(

dt

+

ui

ui

+ Fve

ui + di

di)

ui

Yi )

= B(

+

ui

¯ui

+ Fve

ui

(Zi + 1

Zi )

Fve (Zm + 1

N m=0

+Q

di di

(Yi

di

Z n) 2

14

(A5)

Vyi)2 + (Vzi)2 (vi

Vyi)

Yi 1) Yn)

2

Y n) + (Zm

(Zi

Zm)(Zm

X n )2 + (Ym

(Xm

X n)

Z n) 2

Vxi)2 + (vi ¯ di

Vxi)

Xi 1)

Ym )(Ym

(ui

ui di

Vyi)2 + (Vzi)2 (ui

Vxi)2 + (vi

X n ) + (Ym

ui + di

di)

(Xi

Y n)2 + (Zm

2

(Xm

d Vzi dt

Fve (Ym + 1

N m=0

+Q

di

Xm )(Xm

(ui

ui di

(Yi + 1

di

X n )2 + (Ym

(Xm

d Vyi

Fve (Xm + 1

N m=0

Vxi)2 + (vi

(ui

ui di

14

(A6)

Vyi)2 + (Vzi)2 ( Vzi)

Zi 1) Zn )

Y n)2 + (Zm

Z n) 2

14

(A7)

xi = Vxi t

(A8)

yi = Vyi t

(A9)

zi = Vzi t

(A10)

The above equations involve three dimensionless groups

B=

CD a0 L02 m

(A11)

Fve =

a02 µ2 mL0 G

(A12)

Q=

48 µ2 12 mG 2L013

(A13)

Appendix B. Dynamic mesh refinement Briefly, this method consists of three steps described below. First, the distances between the nodes on a fiber are calculated and compared to a preset reference length ref at fixed time intervals tref . If any of the node-to-node distance i > ref , the specific node will be marked with is and the refinement process begins. [0,1] will be set up based on the following equation Second, a curvilinear coordinate i

=

i

1

m

(B1)

m=1

i m=1 m

where the sum is the partial length of the fiber from i = 0 to i = m . The values of i will be the mesh system instead of the subscript index of the nodes i and all the variables will be tabulated as fi = f ( i ) and used in the calculation of cubic spline. The cubic spline is constructed based on Akima's algorithm [35]. Third, a new mesh system can be set up from the marked node is . Due to the refinement, the numbers of nodes will increase at i > is and the new numbers are calculated as

n =s+

n l m=s+1 m

(B2)

lref

Then the new mesh curvilinear coordinate is defined as i

=

i

i

=

i 1

if i +

(B3)

s

1

s

nref

if i

s

(B4)

n l /l m = s + 1 m ref

where nref = , and it is used for the updated new numbers of the nodes at i > is . Finally, the new values f ( i ) can be calculated by spline interpolation that is based on the step two. This dynamic refinement process has already been implemented in the code that was used to simulate the two-fiber case discussed above (cf. Figs. 13–17). To illustrate the efficiency of this process, several additional simulations were done based on different numbers of nodes. The following 13

Polymer xxx (xxxx) xxx–xxx

G. Li et al.

parameters were used: B = 4 , Fve = 1, Q = 5 × 10 11. The simulation results discussed below correspond to t = 0.014 . Fig. B1a depicts the simulation results at the number of nodes N = 100 without mesh refinement. It is clearly seen that near the entanglement location a tremendous stretching of the Lagrangian mesh appears and long sections of the fibers connect only two nodes and thus are straight. Figure B1b reveals that at the larger number of nodes N = 500 (without mesh refinement), the Lagrangian mesh stretching is significantly diminished, albeit some salient points in the fiber shape are still visible. When mesh refinement is introduced at N = 500 , the fiber shape near the entanglement becomes even smoother, i.e. the mesh refinement is definitely helpful.

Fig. B1. Comparison of the simulation results without and with dynamic mesh refinement. (a) N = 100 (without mesh refinement); (b) N = 500 (without mesh N = 999 (with dynamic mesh refinement). refinement); (c) N = 500

The improvement of the results due to the dynamic mesh refinement in Fig. B1 is obvious and most of the straight sections and salient points vanished in Fig. B1c. In the latter panel the number of the nodes increased from N = 500 to N = 999 and would keep on increasing when the fiber would be kept stretching. It should be emphasized that an increase in the number of nodes due to the dynamic mesh refinement results in more timeconsuming simulations, as well as in an additional complexity in the code. So far, this method has been only adopted for the refinement of a single fiber. Moreover, to guarantee the accuracy of the cubic splines, the calculation time step needs to be set small enough.

[19] [20] [21] [22]

References [1] M.A. Said, V. Thomas, J. Plastic Film Sheeting 13 (1997) 212–220. [2] E. Ghassemieh, M. Acar, H.K. Versteeg, Compos. Sci. Technol. 61 (2001) 681–1694. [3] P. Xiang, Numerical Modelling and Experimental Investigation of the Hydroentanglement Process, PhD Thesis North Carolina State University, 2007. [4] P. Xiang, A.V. Kuznetsov, A.M. Seyam, J. Textil. Inst. 100 (2009) 293–304. [5] A.M. Seyam, D.A. Shiffler, H. Zheng, Int. Nonwoven J. 14 (2005) 25–33. [6] S. Jaganathan, H.V. Tafreshi, B. Pourdeyhimi, Chem. Eng. Sci. 63 (2008) 244–252. [7] L.B. Venu, E. Shim, N. Anantharamaiah, B. Pourdeyhimi, Microsc. Microanal. 18 (2012) 1368–1379. [8] L.B. Venu, E. Shim, N. Anantharamaiah, B. Pourdeyhimi, J. Textil. Inst. 108 (2017) 301–313. [9] E. Ghassemieh, M. Acar, H.K. Versteeg, Part L: J. Mater. Des. Appl. 216 (2002) 199–207. [10] E. Ghassemieh, M. Acar, H.K. Versteeg, Part L: J. Mater. Des. Appl. 216 (2002) 211–218. [11] P. Xiang, A.V. Kuznetsov, A.M. Seyam, J. Porous Media 11 (2008) 35–49. [12] R. Lakes, Viscoelastic Materials, Cambridge University Press, Cambridge, 2009. [13] A.L. Yarin, Free Liquid Jets and Films: Hydrodynamics and Rheology, Longman Scientific & Technical and John Wiley & Sons, Harlow, New York, 1993. [14] V.M. Entov, A.L. Yarin, J. Fluid Mech. 140 (1984) 91–111. [15] J. Li, M.A. Fontelos, Phys. Fluids 15 (2003) 922–937. [16] D.H. Reneker, A.L. Yarin, H. Fong, S. Koombhongse, J. Appl. Phys. 87 (2000) 4531–4547. [17] A.L. Yarin, B. Pourdeyhimi, S. Ramakrishna, Fundamentals and Applications of Micro- and Nanofibers, Cambridge University Press, Cambridge, 2014. [18] M. Lauricella, G. Pontrelli, D. Pisignano, S. Succi, J. Comput. Sci. 17 (2016) 325–333.

[23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]

14

A.L. Yarin, S. Sinha-Ray, B. Pourdeyhimi, J. Appl. Phys. 108 (2010) 034913. S. Sinha-Ray, A.L. Yarin, B. Pourdeyhimi, J. Appl. Phys. 108 (2010) 034912. A.L. Yarin, S. Sinha-Ray, B. Pourdeyhimi, Polymer 52 (2011) 2929–2938. A. Begenir, The Role of Orifice Design in Hydroentanglement, MS Thesis North Carolina State University, 2003. A.E. Scheidegger, Physics of Flow through Porous Media, University of Toronto Press, Toronto, 1974. J. Bear, Dynamics of Fluids in Porous Media, Elsevier, New York, 1972. G.I. Barenblatt, V.M. Entov, V.M. Ryzhik, Theory of Fluid Flows through Natural Rocks, Kluwer, Dordrecht, 1990. L.G. Loitsyanskii, Mechanics of Liquids and Gases, Pergamon Press, Oxford, 1966. G. Polya, G. Latta, Complex Variables, John Wiley & Sons, New York, 1974. P. Henrici, Applied and Computational Complex Analysis, John Wiley & Sons, New York, 1993. G. Astarita, G. Marrucci, Principles of Non-Newtonian Fluid Mechanics, McGrawHill, New York, 1974. R.B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids. Fluid Mechanics vol. 1, John Wiley & Sons, New York, 1987. Y.N. Rabotnov, Creep Problems in Structural Members, North-Holland Publishing Company, Amsterdam, 1969. S.A. Theron, A.L. Yarin, E. Zussman, E. Kroll, Polymer 46 (2005) 2889–2899. L. Verlet, Phys. Rev. 159 (1967) 98–103. J.H. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Springer, Heidelberg, 2012. H. Akima, J. ACM 17 (1970) 589–602. G. Li, C. Staszel, A.L. Yarin, B. Pourdeyhimi, Polymer (2018), https://doi.org/10. 1016/j.polymer.2018.11.004 (in this issue). N. Mao, S.J. Russell, Compos. Sci. Technol. 66 (2006) 80–91.