- Email: [email protected]
Efficient bead-chain model for predicting fiber motion during molding of fiber-reinforced thermoplastics
Accepted Manuscript
Efficient Bead-Chain Model for Predicting Fiber Motion during Molding of Fiber-Reinforced Thermoplastics Toshiki Sasayama , Masahide Inagaki PII: DOI: Reference:
S0377-0257(18)30219-2 https://doi.org/10.1016/j.jnnfm.2018.10.008 JNNFM 4057
To appear in:
Journal of Non-Newtonian Fluid Mechanics
Received date: Revised date: Accepted date:
22 June 2018 4 September 2018 18 October 2018
Please cite this article as: Toshiki Sasayama , Masahide Inagaki , Efficient Bead-Chain Model for Predicting Fiber Motion during Molding of Fiber-Reinforced Thermoplastics, Journal of Non-Newtonian Fluid Mechanics (2018), doi: https://doi.org/10.1016/j.jnnfm.2018.10.008
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights An improved model for direct simulation of fibers in viscous flow is developed.
Fiber motions in simple shear flow are compared between the proposed and conventional model.
Fiber motions in injection molding are also simulated.
The improved model well reproduces the results of the conventional one. The model also provides higher computational efficiency than the conventional one.
AC
CE
PT
ED
M
AN US
CR IP T
1
ACCEPTED MANUSCRIPT
Efficient Bead-Chain Model for Predicting Fiber Motion during Molding of Fiber-Reinforced Thermoplastics
a
Toyota Central R&D Labs., Inc.
CR IP T
Toshiki Sasayamaa,*, Masahide Inagakia
* Corresponding Author (T. Sasayama) Toyota Central R&D Labs., Inc., 41-1, Yokomichi, Nagakute, Aichi 480-1192, JAPAN
AN US
Email address: [email protected] Tel: +81-561-71-7487 Abstract
This paper proposes an efficient model for simulating fiber motion in viscous flow,
M
especially during molding of fiber-reinforced thermoplastics. To increase the
ED
computational efficiency, the proposed model discretizes fibers as chains of fewer spheres than in the bead-chain model previously proposed by the authors, which is
PT
called the Simplified Bead-Chain Model (SBCM). In addition, in order to make the
CE
proposed model reproduce the results of the SBCM, correction factors applied to forces
AC
exerted on each sphere are newly derived so that the equation of motion of a fiber for the proposed model can be equivalent to that of the conventional model. Simulations of rigid and flexible fibers in simple shear flow are first performed to compare the fiber motions between the two models. The proposed model reproduces the results of the SBCM closely for cases of both a single fiber and suspensions containing multiple
2
ACCEPTED MANUSCRIPT
fibers. Finally, fiber orientation behavior during injection molding in a simple-shaped plaque is simulated. The results obtained with the proposed model shows good
CR IP T
agreement with those of the SBCM but at lower computational cost.
Keywords
AN US
Direct simulation; Fiber modeling; Bead-chain model;
M
Fiber-reinforced plastics;
AC
CE
PT
ED
Injection molding.
3
ACCEPTED MANUSCRIPT
1. Introduction To predict the mechanical and failure properties of injection- or compression-molded fiber-reinforced thermoplastics, it is essential to accurately predict
CR IP T
the distribution of orientations and lengths of fibers during the molding process.
Especially, the demand for predicting these properties of long-fiber composites is
AN US
increasing. Long fibers generally undergo large deformation due to flow resistance and interactions with other fibers during the molding processes, resulting in more complicated microstructures (e.g., fiber orientation and length) than those of short fiber
M
composites.
ED
Direct fiber simulations are useful for understanding the behavior of fibers during the molding process in detail (e.g., rheological properties [1-6], fiber orientation
PT
[7-9], and fiber breakage [10]). Concerning fiber orientation, in a practical injection
CE
molding simulation, statistical distributions of the fiber orientation are predicted using a
AC
macroscopic model describing the evolution of the fiber orientation tensor [11], in which the effect of fiber-fiber interactions is introduced as a diffusion term with some phenomenological parameters [12-14]. To investigate the validity of the macroscopic model, Mezher et al. [7] employed direct simulations and compared the obtained orientation tensor of fiber suspensions in simple shear flow with that predicted by the
4
ACCEPTED MANUSCRIPT
macroscopic model. As more practical applications, Yashiro et al. [8, 9] employed a moving particle semi-implicit method [15] for predicting behaviors of fiber orientation and accumulation during injection molding. Kuhn et al. [16, 17] simulated the
CR IP T
compression molding process of a ribbed plate to investigate the distribution of the fiber content within the rib.
AN US
Various models have been proposed for discretizing a fiber. Flexible fibers have been modeled as a chain of spheres [18-21], prolate spheroids [22], and rods
[23-25]. In our previous study [26], we developed an improved model, the Simplified
M
Bead-Chain Model (SBCM), using the work of Yamamoto and Matsuoka [18, 19] to
ED
enhance its computational efficiency by simplifying the governing equations. However, the SBCM has the drawback that the number of spheres required to model a fiber is
PT
equal to its aspect ratio. Therefore, the SBCM is computationally expensive for
CE
simulating long-fiber suspensions compared to the previous models [21-25], which use
AC
fewer components (such as rods and prolate spheroids) for discretizing fibers. It is thus desirable to enhance the computational efficiency of the SBCM by reducing the number of components needed to model a fiber to the levels of the models mentioned above. However, the mass and the moment of inertia of the fiber changes as the number of spheres decreases, which results in numerical errors in the fiber motion. Hence, it is also
5
ACCEPTED MANUSCRIPT
essential to establish a new method to decrease these errors. The purpose of this study is to develop a computationally efficient model that reproduces results as accurate as those of the SBCM. To this end, in the proposed model,
CR IP T
a fiber is modeled with fewer spheres than in the previous SBCM and the forces acting on each sphere are corrected by multiplying them with correction factors, which are
AN US
derived so that the equation of motion of a fiber in the proposed model becomes
equivalent to that in the SBCM. Fiber motions in a simple shear flow are then simulated and compared between the proposed model and the SBCM. Finally, simulations in an
PT
2. Numerical method
ED
model to complex flow.
M
injection molding process are carried out to investigate the applicability of the proposed
CE
2.1. Simplified Bead-Chain Model (SBCM)
AC
In the Bead-Chain Model (BCM) [18, 19], a fiber is modeled as chain of
spheres, as illustrated in Fig. 1(a), where the distance between the centers of two adjacent spheres is equal to the fiber diameter (i.e., closely packed). The SBCM solves only the equation of translational motion of each sphere, in contrast to the BCM, which solves both translational and rotational equations of motion. Thus, the governing
6
ACCEPTED MANUSCRIPT
equation of a sphere is given as
m
dvi Fi h FijS Fijb Fijht FinP , dt j 1 j 1 j 1 n 1
(1)
CR IP T
where m is the mass of the sphere, vi is the velocity of sphere i, the index j indicates the spheres connected to sphere i, the index n indicates the spheres that are not connected to sphere i, Fh is the hydrodynamic force, FS is the stretching force, FP is the
AN US
particle-particle interaction force, Fb is the bending force converted from the bending torque, and Fht is the force converted from the hydrodynamic torque. The SBCM
includes these converted forces for expressing the effect of bending and hydrodynamic
M
torques in the equation of translational motion as in Eq. (1), instead of the equation of
ED
rotational motion as in the BCM. Thus, the SBCM does not need to solve the equation
CE
paper [26].
PT
of rotational motion, and it ignores twisting fiber motion, as described in our previous
AC
2.2 Improvement of Simplified Bead-Chain Model To improve the efficiency of the SBCM, the number of spheres used for
modeling a fiber is reduced. In addition, a correction process is developed for eliminating the errors in fiber motions caused by reducing the number of spheres. First, we consider a fiber modeled with fewer spheres than in the SBCM, as 7
ACCEPTED MANUSCRIPT
shown in Fig. 1(b). Here, the fiber is considered to consist of Nseg segments, each of which is composed of two spheres and a massless rod connecting them. In addition, NR is defined as the number of reduced spheres in a segment (Fig. 1). The indexes of
CR IP T
spheres and segments are assigned sequentially, as indicated in Fig. 1(b). The indices at the upper side of the fiber in Fig. 1(b) are for spheres, and those in the lower side are for
AN US
segments.
Next, the correction procedure to reduce the numerical errors is described. The behavior of the fiber shown in Fig. 1(b) differs from that in Fig. 1(a) without any
M
corrections because the mass and moment of inertia are smaller. This paper proposes the
ED
procedure described below to correct the forces acting on spheres and eliminate these errors. For convenience, the fiber in Fig. 1(b) is referred to as a “reduced fiber” in this
PT
paper. Now, we consider a fiber made of closely packed spheres with the same
CE
conformation as the reduced fiber, as shown in Fig. 2. This paper refers to the fiber in
AC
Fig. 2(b) as an “unreduced fiber”. In each segment of the unreduced fiber, NR spheres (illustrated with dotted outlines in Fig. 2(b)) exist aligned. The unreduced fiber is an imaginary fiber used only for deriving correction factors for forces of spheres in the reduced fiber. The correction factors are derived so that the equations of motion for both fibers become equivalent. In constructing the equations of motion of these fibers, the
8
ACCEPTED MANUSCRIPT
following suppositions are made: (a) each segment is rigid, (b) lengths (i.e., NR) of the segments in each fiber are equal, and (c) distributed forces (e.g., hydrodynamic forces) vary linearly over a segment.
CR IP T
The forces acting on spheres are classified into two types, “distributed” and “non-distributed.” The term “non-distributed forces” means that these forces can be
AN US
regarded as acting only on the spheres at both ends of each segment. Let FD and FND be the distributed and the non-distributed forces, respectively. The distributed force FD includes the hydrodynamic force Fh. In contrast, the non-distributed force FND includes
M
the stretching force FS, the bending force Fb, and the force coming from the
ED
hydrodynamic torque Fht. Note that Fht can be approximated as non-distributed forces provided that the fluid velocity field is linear. The particle-particle interaction forces FP
PT
are assumed to be non-distributed forces in this study, although they are sometimes
AC
CE
distributed throughout the segment. In summary, FD and FND are given as
Fi ND FijS Fijb Fijht Fi Pns ,
(2)
Fi D Fi h .
(3)
j 1
j 1
j 1
ns 1
Note that in the fourth term on the right-hand side of Eq. (2), the index of FP changes from n to ns, which denotes the segment index. That is, interaction forces on sphere i are exerted by segment ns (where ns denotes a segment different from the segment that 9
ACCEPTED MANUSCRIPT
includes sphere i), in contrast to the SBCM. Here, the calculation procedure for interaction forces is briefly described; additional details can be found in Appendix A. In the proposed model, interaction forces are calculated between segments instead of
CR IP T
spheres, as in the SBCM. The calculated forces, which are acting on arbitrary points in the segments, are then converted into equivalent forces that are exerted on the two
AN US
spheres at both ends of each segment.
The equations of translational and rotational motions for each fiber are given in Appendix B. Using these equations, the correction factors for FD and FND are derived as
M
follows so that the equations for the reduced and unreduced fibers become equivalent.
ED
First, focusing on the equations of translational motion, each force exerted on sphere i is corrected as
PT
Fi ND T Fi ND ,
Fi D iT Fi D ,
(4)
CE
where F D and F ND are the corrected forces and T and T are the correction
AC
factors for the translational motion of reduced fibers given by
T
N seg 1 , N seg N R 1 1
(5)
N R 2 N seg 1 N seg N R 1 1 2 iT R seg N 1 N 1 N seg N R 1 1
i 1, N
seg
1
2 i N seg
10
.
(6)
ACCEPTED MANUSCRIPT
Next, focusing on the equations of rotational motion, the correction factors R and
R are obtained as follows: I Rf If
,
(7)
CR IP T
R
I Rf N R 2 N seg 1 3 N R 1 1 N R nˆ1 N 1 N 1 N 1 f seg R 6 N 1 N 1 I Rf 1 2 1 I iR f N R N R 2 nˆ i nˆ i 1 N i3 N i3 N i3 6 N R 1 R I 6 N 1 I Rf 1 N R 2 N seg 1 3 N R 1 N R nˆ N seg N 2 N 2 N 2 f R seg I 6 N 1 N 1
where N seg
i 1
2 i N , seg
i N
AN US
seg
N 1 N seg p 1 nˆ p , N seg
1
(9)
M
p 1
(8)
N pnˆ p ,
ED
2
(10)
p 1
AC
CE
PT
N i3 N i0
I Rf I
f
N
1 seg
1
N1,
(11)
i 1
N i0 nˆ p ,
(12)
p 1
N seg 1
3 rqG rqG q 1
N seg
,
3rNGseg 1 rNGseg 1 3 N R 1 rqG rqG 2aN R rqG nˆ q 2a 2 N seg N R N R 1 2 N R 1 q 1
(13) where riG is a vector from the center of gravity of the fiber to the position of sphere i 11
ACCEPTED MANUSCRIPT
and nˆ i is the unit vector from sphere i to sphere i+1. Throughout this correction process, the notation of this unit vector is different from that used in our previous paper [26] to simplify the equations. Specifically, the unit vector from sphere i to i+1 was
CR IP T
originally written as ni i 1 in [26], but now it is defined as nˆ i ni i 1 .
However, it is not easy to fully correct both the translational and rotational
AN US
motions of the fibers at the same time because the correction factors for each motion are different. This study thus employs the following procedure for correcting translational and rotational motions of the reduced fibers. Let niG be a unit vector from the center
M
of gravity of the fiber to the position of sphere i, and then FD and FND are decomposed
ED
into parallel and perpendicular components to niG , respectively. The correction factors for translational and rotational motions are multiplied to the parallel and the
PT
perpendicular components as
(14)
(15)
CE
Fi ND T Fi ND niG niG R Fi ND Fi ND niG niG ,
iT Fi D niG niG iR Fi D Fi D niG niG ,
AC
Fi D
so that at least the rotational motion of fibers is fully corrected. The calculation procedure of the reduced fiber described above can be
summarized as follows: (i) Calculate the correction factors using Eqs. (5)-(13).
12
ACCEPTED MANUSCRIPT
(ii) Calculate the distributed and non-distributed forces (FD and FND) using Eqs. (2) and (3). (iii) Correct the forces calculated in step (ii) using Eqs. (14) and (15).
m
dv i Fi D Fi ND . dt
(16)
AN US
Note that the model is reduced to the SBCM when NR=0.
CR IP T
(iv) Solve the following equation to update the coordinates and velocities of spheres:
3. Results and discussions
M
The effectiveness of the proposed model was investigated by comparing the
ED
fiber behavior results with those calculated with the SBCM. In this study, the number of segments Nseg was employed to describe the discretization of a fiber. A fiber of aspect
PT
ratio rf was discretized with Nseg=rf 1 when applying the SBCM, whereas 1 Nseg
CE
1 was used in the proposed model. In addition, this study assumed that the fiber and
AC
fluid motions are decoupled. The coupling between fibers and fluid is work for the future.
3.1. Fiber motion in simple shear flow Here, we describe investigation of the motion of fibers in simple shear flow.
13
ACCEPTED MANUSCRIPT
The velocity field of the flow at position x was given as V=( z, 0, 0), where the shear rate was set to =100 s-1. The radius of the fiber and the viscosity of the fluid were a=5 μm and η=100 Pa∙s, respectively. The density of fibers ρf was determined so that the
CR IP T
Reynolds number of a sphere (Re=a2 ρf /η) becomes Re=0.1 so that the effect of inertia can be ignored. Fibers were initially oriented parallel to the x axis (i.e., the flow
AN US
direction).
First, periods of rotation of a rigid fiber were compared. Here, particle-particle interaction forces were not calculated. The aspect ratio of the fiber was rf=100, and
M
Young’s modulus E was large enough to not induce bending deformation. Figure 3
ED
shows the relative error in periods of rotation calculated as 100×|T-TSBCM|/TSBCM, where TSBCM and T are the periods of rotation obtained by the SBCM (Nseg=99) and the
PT
proposed model (Nseg<99), respectively. The filled and open symbols indicate the results
CE
with and without, respectively, the correction process described in Section 2.2. When
AC
the correction process is not applied, the errors increase with decreasing number of segments. In contrast, these errors are mostly eliminated when the correction is applied. Figure 4 represents the ratio of computational time as a function of the number of segments. The calculations in Fig. 4 were performed with the same time step (Δt) for each number of segments. The computational efficiency is found to increase as the
14
ACCEPTED MANUSCRIPT
number of segments decreases. In this case, the calculation becomes about 43 times faster when the number of segments is reduced from 99 to 1. Next, the motion of flexible fibers was compared. Here, the aspect ratio and the
CR IP T
Young’s modulus were rf=50 and E=0.5 GPa, respectively. Here, calculating the
effective stiffness Seff (= E / 64 rf4 ) [27] is useful for evaluating the flexibility of a
AN US
fiber suspended in shear flow. In this case, the effective stiffness was about Seff
=3.9×10-4, which is low enough to cause bending deformation [23, 27]. Figure 5(a) shows the time evolutions of the z-coordinate of the sphere at the fiber end calculated
M
with Nseg = 49 (SBCM), 25, and 15. Figure 5(b) shows snapshots of the fibers at three
ED
time steps. As can be seen in Fig. 5(b), the fiber induces S-shaped deformation during rotation, as reported in our previous study [26]. The results match each other well for
PT
flexible fibers regardless of the number of segments. Of course, not even the proposed
CE
model could trace the deformation of flexible fibers if the number of segments is
AC
insufficient to resolve the deformed fiber shape. Finally, fiber suspensions were simulated considering particle-particle
interactions. The analysis region was a 500-μm cube in which the top and bottom planes in the z direction were wall boundaries and a periodic boundary condition was applied in the x and y directions. The simulations were performed for various fiber volume
15
ACCEPTED MANUSCRIPT
fractions Vf and aspect ratios rf. Fiber orientation tensors Aij [11] were calculated using the directional vectors of each segment. The simulations were run until the fiber orientation tensors almost reached a steady state. In all cases, the fiber suspensions are
CR IP T
homogeneously dispersed and show no heterogeneous structures (i.e., flocculation). This is because this study does not consider both the inter-fiber friction and the
AN US
non-straight equilibrium fiber shape, which are the factors controlling the formation of flocculation [23, 27]. Figure 6 shows the component of a fiber orientation tensor Axx at steady state as a function of Vf rf. Figure 7 plots the average number of interactions per
M
fiber at steady state. The number of interactions is defined as the total occurrences of
ED
lubrication and repulsive forces applied to spheres (for the SBCM) or segments (for the proposed model). Here, the product of the fiber volume fraction Vf and the aspect ratio rf
PT
can be used to estimate whether the semi-concentrated or concentrated regime
CE
dominates the suspension [28]. In Figs. 6 and 7, the black-filled symbols indicate the
AC
results obtained by the SBCM, whereas the grey-filled and open symbols are obtained by the proposed model. The proposed model gives results comparable to those of the SBCM in terms of both the orientation tensor (Fig. 6) and the number of interactions (Fig. 7). In addition, the difference in the orientation tensor between these models shown in Fig. 6 becomes smaller with increasing Vf rf. Figure 8 shows the transient 16
ACCEPTED MANUSCRIPT
behaviors of Axx and Ayy for various aspect ratios rf and volume fractions Vf. The horizontal axis in Fig. 8 is the applied shear strain defined as the product of the shear strain rate and time t. For rf=10 and Vf=0.122 (Fig. 8(a)), the initial behavior slightly
CR IP T
differs between the models depending on the number of segments. However, the differences in both the initial and steady-state regions decrease when Vf=0.201
AN US
(Fig. 8(b)). For rf=20 and Vf=0.056 (Fig. 8(c)), the curves mostly match, although Vf rf is less than that in the case of Fig. 8(a). It is also found that the difference between both models gets larger when Vf rf is less than or equal to unity (not show in the figure).
M
These differences in the orientation behavior between the two models that occur at Vf rf
ED
≤ 1 could be due to the translational motion of fibers not being fully corrected, unlike the rotational motion correction described in Section 2.2. However, these differences in
PT
the translational motion of fibers are considered to decrease when fibers experience a
CE
large number of interactions (i.e., in the concentrated regime). It is thus concluded that
AC
the proposed model reproduces the behavior of fiber suspensions simulated by the SBCM when Vf rf is higher than unity. Figure 9 shows the time evolution of Axx and Ayy in the case of flexible fibers with Young’s modulus E=0.1 GPa. The curves calculated with those models show good agreement in both cases of (rf, Vf)=(20, 0.056) and (30, 0.037).
17
ACCEPTED MANUSCRIPT
3.2. Fiber motion in injection molding process The proposed model was applied to the simulation of fiber motion in a
CR IP T
simple-shaped plaque during injection molding to compare the predicted fiber orientations with those of the SBCM.
AN US
First, a mold-filling simulation for the plaque illustrated in Fig. 10 was
performed using the commercial molding software 3D TIMON (Toray Engineering, Tokyo, Japan). The material was polypropylene with 20 wt% long glass fibers (L-2040P,
M
Prime Polymer, Tokyo, Japan). The inlet and mold temperatures were 230°C and 110°C,
ED
respectively. The flow rate at the inlet was 10 cm3/s. No holding pressure was applied in the mold-filling simulation. Figure 11 shows the flow velocity distribution in the
PT
thickness direction at each region (A, B, and C) delineated with dashed squares in Fig.
AC
inlet.
CE
10. The regions A, B, and C were located 15, 55, and 90 mm, respectively, from the
Then, a direct simulation of fibers was carried out using the transient flow field
obtained by the mold-filling simulation. The fiber radius a, length Lf, Young’s modulus E, and volume fraction Vf were a=10 μm, Lf=1 mm, E=70 GPa, and Vf=0.1, respectively. The fiber orientation at the inlet was set to be nearly aligned in the x (flow) direction.
18
ACCEPTED MANUSCRIPT
The number of segments in a fiber was set to Nseg=13 in the proposed model, whereas Nseg=49 in the SBCM. In calculating fiber-fiber interactions, only repulsive forces were considered to reduce the computational cost. Figure 12 shows the distribution of fiber
CR IP T
orientation tensor components Axx and Ayy along the z (thickness) direction in each
region. Figure 13 depicts fibers at region C simulated with the proposed model (where
AN US
only segments within the region of 2×2×1 mm are shown). In Figs. 12(b), 12(c), and 13, one can observe a layered structure, in which fibers in the shell layer are aligned in the flow direction, whereas those in the core layer show nearly in-plane random orientations.
M
The fibers in the core layer would become transversely aligned to the flow direction
ED
when using a mold shape that causes a diverging flow stronger than that simulated in this study. These results confirm that both models could capture the trends of fiber
PT
orientation that are seen in the injection molding process and that the proposed model
CE
works well even in a complicated flow field. Compared to the SBCM, the proposed
AC
model enables us to use a time step about 3.7 times larger, which yields a calculation about 2.9 times faster. Although the degree of speed-up is smaller than that for single-fiber calculations due to the existence of fiber-fiber interactions, which are computationally expensive, the proposed model succeeds in reducing the computational cost, even when considering these interactions.
19
ACCEPTED MANUSCRIPT
4. Conclusions This study proposed an efficient model based on the SBCM by reducing the
CR IP T
number of spheres needed to model a fiber and correcting the forces exerted on each sphere to reduce the computational cost without loss of accuracy compared to the
AN US
SBCM.
The correction process was applied to simulations of fiber motion in simple shear flow to investigate the accuracy of the proposed model in comparison to the
M
SBCM. First, periods of rotation of single rigid fibers and deformation of single flexible
ED
fibers were compared between these models. The proposed model mostly matched the results of the SBCM, in contrast to the large differences that occurred without the new
PT
correction process. In addition, it was found that the computational efficiency increased
CE
as the number of segments decreased. Second, suspensions of both rigid and flexible
AC
fibers were simulated with consideration of fiber-fiber interaction. The fiber orientation tensors of the suspensions predicted by the proposed model and the SBCM showed good agreement provided that the suspensions were in the concentrated regime (i.e., the product of the fiber volume fraction and the aspect ratio Vf rf is greater than unity). Finally, orientation behavior of flexible fibers within a simple-shaped plaque
20
ACCEPTED MANUSCRIPT
during injection molding was compared between the models. Both models predicted orientations with layered structure that are typically observed in injection molding, and the proposed model mostly reproduced the results obtained with the SBCM.
CR IP T
The validity and the effectiveness of the proposed model were confirmed by
the investigations in this study. Use of the proposed model will produce more realistic
AN US
simulations of the motion of long fibers in a highly concentrated state during molding of fiber-reinforced thermoplastics. Our future work will focus on simulating more practical molding processes by considering fiber breakage phenomena and two-way coupling
ED
M
between fibers and fluid.
Appendix A. Calculation procedure for particle-particle interaction forces
PT
This appendix describes the procedure for calculating particle-particle
CE
interaction forces (repulsion and lubrication) in the proposed model. First, segments are
AC
regarded as capsules for detecting collisions between them (Fig. A.1). When a fiber consists of one segment, a capsule bounding volume is set as illustrated in Fig. A.1(a) with dashed lines. In other cases, segments in the fiber have capsule bounding volumes as shown in Fig. A.1(b) with dashed lines. With these assumptions, repulsive and lubrication forces are calculated as described below.
21
ACCEPTED MANUSCRIPT
(i) Repulsive forces Consider two closely located segments indexed as is and ns, as shown in
force acting on the segment is is calculated by x c xisc 1 ns FisRseg D exp G ns 0 0 2a
x c x c is ns c xns xisc
CR IP T
Fig. A.2(a), where the segment is consists of spheres indexed as i and j. The repulsive
c ( xns xisc 2.001a) ,
(A.1)
AN US
c where D0 and G0 are constants and xisc and xns are the coordinates of the closest
points of the two segments is and ns, respectively. The repulsive forces acting on spheres i and j are determined as
respectively.
x j xi and xi and x j are the coordinates of spheres i and j,
PT
where s xisc xi
(A.2)
ED
M
Rep Rseg Fi ns 1 s Fis ns , Rep Rseg F s F j ns is ns
CE
(ii) Lubrication forces
AC
When calculating the lubrication forces between two segments, this study
accounts for lubrication effects caused by imaginary spheres (illustrated in Fig. A.2(b) with dotted lines), as well as the two spheres at both ends of each segment. Let p and q be the indexes of spheres in segments is and ns, respectively. The lubrication force applied to sphere p by sphere q is given by
22
ACCEPTED MANUSCRIPT
FpqLseg
x xp 3 2 vq v p a q 2 x q x p 2a x q x p
xq x p xq x p
(2.001a xq x p 3a) ,
(A.3)
where η is the viscosity of the fluid. When p (or q) is an imaginary sphere, the velocity
CR IP T
vp of the sphere is linearly interpolated using the velocities of the spheres at both ends of the segment. The lubrication forces acting on the spheres (i and j) of the segment is are calculated as
where s p x p xi
AN US
Lseg Fi Lub ns 1 s p Fpq p q , Lub Lseg F s F j ns p pq p q
x j xi .
(A.2)
M
Finally, the particle-particle interaction forces of the spheres (i and j)
ED
exerted by the segment ns are given by
(A.3)
CE
PT
P Rep Lub Fi ns Fi ns Fi ns . P Rep Lub F F F j ns j ns j ns
AC
Appendix B. Equations of fiber motion The equations of translational motion for the center of mass of the reduced
(Fig. 2(a)) and unreduced (Fig. 2(b)) fibers are given by
N
seg
1 m
dv Rf dt
N seg 1
i 1
T
Fi
ND
N seg 1
i 1
T i
Fi D ,
23
(B.1)
ACCEPTED MANUSCRIPT
N N seg
R
1 1 m
dv f dt
N seg 1
F i 1
ND
i
NR 2 D N 1 Fi F1 FNDseg 1 , 2 i 2
N seg 1
R
D
(B.2)
where vRf and v f are the velocity of the center of mass of the reduce and unreduced
CR IP T
fibers, respectively, and Nseg(NR+1)+1 in the left-hand side of Eq. (B.2) is the number of spheres in the unreduced fiber, including the imaginary spheres illustrated in Fig. 2(b) with dotted lines.
AN US
To write the equation of rotational motion for the center of mass of all fibers, we use the body-fixed coordinate system, in which the origin is fixed in the center of mass of the fibers and the axes are along the principal axis of inertia. The equation of
M
rotational motion of the reduced fiber is written as
~ Rf dI~ Rf N seg 1 ~ Rf dω ~ Rf ~ Rf ~ ~ Rf Rf ~ ~ I ω ω I ω ~ ri G R Fi ND iR Fi D , dt dt i 1
ED
(B.3)
PT
~ ~ Rf are the inertia tensor and the angular velocity of the reduced where I Rf and ω
CE
fiber, respectively, and the tilde ~ means that the quantity is written in the body-fixed
AC
coordinate system. Assuming that the fiber maintains the same shape during time step Δt, the second term in the left-hand side of Eq. (B.3) can be ignored, and then the equation becomes
~ Rf N 1 ~ Rf dω ~ ~ ~ Rf I~ Rf ω ~ Rf ~ I ω ri G R Fi ND iR Fi D . dt i 1
seg
The right-hand side of Eq. (B.4) is calculated as 24
(B.4)
ACCEPTED MANUSCRIPT
N seg 1
N~
~ ~ r~i G R Fi ND iR Fi D 2 N R 1 a i 1
2N 1a ~ N
N seg 1 i 1
0 i
1 N
seg
~ ~ N 1 R Fi ND 1
N 1 2 N 1 a ~1 1 ~ ~ ~0 R ~D R N F 2 N 1 a N 1 iR Fi D . N i seg i 1 seg N 1 N 1 i 2 R
R
N
seg
2
1
~ iR FNDseg 1
seg
(B.5)
diagonal and described as
0 ~ Rf I yy 0
0 0 , ~ I zzRf
(B.6)
AN US
~ I Rf
~ I xxRf 0 0
CR IP T
Here, since the coordinates are along the principal axes of inertia, the inertia tensor is
~ ~ ~ where I xxRf , I yyRf , and I zzRf are the principal moments of inertia. Likewise, the
M
equation of rotational motion of the unreduced fiber is derived as
~f N seg N R 11 ~ f dω ~f ~f ~ ~ f ~ I ω I ω ~ ri G Fi ND Fi D , dt i 1
(B.7)
ED
PT
~ where I f is the inertia tensor of the unreduced fiber. The right-hand side of Eq. (B.7) is calculated as
AC
i 1
N 1 1 ~ ~ ~ ~ ~ r~i G Fi ND Fi D 2 N R 1 a N i0 seg N 1 Fi ND N 1 i 1 R R N 2a 3 N 1 ~1 ~ D R~ N F1 N nˆ1 seg 3 N 1
CE
N seg N R 1 1
seg
~ ~ 2 ~ a N 1 1 ~ ~ N R N R 2 nˆ i nˆ i 1 6 N R 1 N i0 seg N 1 Fi D 3 i 2 N 1
seg
N
R
.
2a 3 N R 1 ~ 2 ~ D R~ ˆ N FN seg 1 N nN seg seg 3 N 1 (B.8)
25
ACCEPTED MANUSCRIPT
To make the motions of the reduced and unreduced fibers equivalent, the following relations need to be satisfied:
~ Rf ω ~f . ω
dt
~f dω dt
.
(B.10)
CR IP T
~ Rf dω
(B.9)
Here, we assume that the following holds for the ratios of principal moments of inertia between the reduced and unreduced fibers for each component:
0 ~ Rf ~ f I yy I yy 0
M
where δ is a unit tensor and
I
Rf
m
N seg 1
~r q 1
ED
Rf I f δ, ~f I I zz
AN US
~ ~ I xxRf I xxf ~ ~ 1 I Rf I f 0 0
G q
~ I zzRf
0 0
(B.11)
~ rqG ,
(B.12)
N N ~ ~ I m ~ rqG ~ rqG m ~ rqG 2apnˆ q ~ rqG 2apnˆ q . i 1 q 1 p 1 N seg 1
seg
R
(B.13)
PT
f
CE
The second term in the right-hand side of Eq. (B.13) is for the imaginary spheres
AC
illustrated in Fig. 2(b). In addition, the right-hand side of Eq. (B.13) is rearranged when it appears in Eq. (13). Finally, by substituting Eqs. (B.4) and (B.7) into Eq. (B10) and using
Eqs. (B.9) and (B.11), the resulting relation for deriving the correction factors αR and βR is given as
26
ACCEPTED MANUSCRIPT
I Rf G R ~ ND R ~D ~ r F F i i i i If i 1
N seg 1
N seg N R 1 1
i 1
~ ~ ~ ri G Fi ND Fi D .
(B.14)
CR IP T
References [1] L. H. Switzer III, D. J. Klingenberg, Rheology of sheared flexible fiber suspensions via fiber-level simulations, J. Rheol. 47 (2003) 759-778.
AN US
[2] M. Yamanoi, J. Maia, Analysis of rheological properties of fibre suspensions in a Newtonian fluid by direct fibre simulation. Part 1: Rigid fibre suspensions, J. Nonnewton. Fluid. Mech. 165 (2010) 1055-1063.
M
[3] M. Yamanoi, J. Maia, T. Kwak, Analysis of rheological properties of fibre
ED
suspensions in a Newtonian fluid by direct fibre simulation. Part 2: Flexible fibre
PT
suspensions, J. Nonnewton. Fluid. Mech. 165 (2010) 1064-1071.
CE
[4] J. Andrić, S. B. Lindström, S. Sasic, H. Nilsson, Rheological properties of dilute suspensions of rigid and flexible fibers, J. Nonnewton. Fluid. Mech. 212 (2014) 36-46.
AC
[5] R. Mezher, E. Abisset-Chavanne, J. Férec, G Ausias, F. Chinesta, Direct simulation of concentrated fiber suspensions subjected to bending effects, Modelling Simul. Mater. Sci. Eng. 23 (2015) 055007. [6] C. G. Joung, N. Phan-Thien, X.J. Fan, Viscosity of curved fibers in suspension, J. Nonnewton. Fluid. Mech. 102 (2002) 1-17. 27
ACCEPTED MANUSCRIPT
[7] R. Mezher, M. Perez, A. Scheuer, E. Abisset-Chavanne, F. Chinesta, R. Keunings, Analysis of the Folgar & Tucker model for concentrated fibre suspensions in unconfined and confined shear flows via direct numerical simulation, Composites Part
CR IP T
A. 91 (2016) 388-397.
[8] S. Yashiro, T. Okabe, K. Matsushima, A numerical approach for injection molding of
AN US
short-fiber-reinforced plastics using a particle method, Adv. Compos. Mater. 20 (2011) 503-517.
[9] S. Yashiro, H. Sasaki, Y. Sakaida, Particle simulation for predicting fiber motion in
M
injection molding of short-fiber-reinforced composites, Composites Part A. 43 (2012)
ED
1754-1764.
[10] S. Yamamoto, T. Matsuoka, Dynamic simulation of flow-induced fiber fracture,
PT
Polym. Eng. Sci. 35 (1995) 1022-1030.
CE
[11] S. G. Advani, C. L. Tucker III, The use of tensors to describe and predict fiber
AC
orientation in short fiber composites, J. Rheol. 31 (1987) 751–784. [12] F. Folgar, C. L. Tucker III, Orientation behavior of fibers in concentrated suspensions, J. Reinf. Plast. Compos. 3 (1984) 98-119. [13] J. H. Phelps, C. L. Tucker III, An anisotropic rotary diffusion model fort fiber orientation in short- and long-fiber thermoplastics, J. Nonnewton. Fluid. Mech. 156
28
ACCEPTED MANUSCRIPT
(2009) 165-176. [14] H. C. Tseng, R. Y. Chang, C. H. Hsu, An objective tensor to predict anisotropic fiber orientation in concentrated suspensions, J. Rheol. 60 (2016) 215-224.
incompressible fluid, Nucl. Sci. Eng. 123 (1996) 421-434.
CR IP T
[15] S. Koshizuka, Y. Oka, Moving-particle semi-implicit method for fragmentation of
AN US
[16] C. Kuhn, I. Walter, O. Taeger, T. A. Osswald, Experimental and numerical analysis of fiber matrix separation during compression molding of long fiber reinforced thermoplastics, J. Compos. Sci. 1(1) (2017) 2.
M
[17] C. Kuhn, I. Walter, O. Taeger, T. A. Osswald, Simulative prediction of fiber-matrix
ED
separation in rib filling during compression molding using a direct fiber simulation, J. Compos. Sci. 2(1) (2018) 2.
PT
[18] S. Yamamoto, T. Matsuoka, A method for dynamic simulation of rigid and flexible
CE
fibers in a flow field, J. Chem. Phys. 98 (1993) 644-650.
AC
[19] S. Yamamoto, T. Matsuoka, Dynamic simulation of fiber suspensions in shear flow, J. Chem. Phys. 102 (1995) 2254-2260. [20] C. G. Joung, N. Phan-Thien, X.J. Fan, Direct simulation of flexible fibers, J. Nonnewton. Fluid. Mech. 99 (2001) 1-36. [21] H. F. Guo, B. G. Xu, A 3D numerical model for a flexible fiber motion in
29
ACCEPTED MANUSCRIPT
compressible swirling airflow, Comput. Model. Eng. Sci. 61 (2010) 201-222. [22] R. F. Ross, D. J. Klingenberg, Dynamic simulation of flexible fibers composed of linked rigid bodies, J. Chem. Phys. 106 (1997) 2949-2960.
CR IP T
[23] C. F. Schmid, L. H. Switzer, D. J. Klingenberg, Simulations of fiber flocculation: Effects of fiber properties and interfiber friction, J. Rheol. 44 (2000) 781-809.
AN US
[24] G. Wang, W. Yu, C. Zhou, Optimization of the rod chain model to simulate the
motions of a long flexible fiber in simple shear flows, Eur. J. Mech. Fluids/B. 25 (2006) 337-347.
M
[25] S. B. Lindstöm, T. Uesaka, Simulation of the motion of flexible fibers in viscous
ED
fluid flow, Phys. Fluids. 19 (2007) 113307.
[26] T. Sasayama, M. Inagaki, Simplified bead-chain model for direct fiber simulation
PT
in viscous flow, J. Nonnewton. Fluid. Mech. 250 (2017) 52-58.
CE
[27] L. H. Switzer, D. J. Klingenberg, Flocculation in simulations of sheared fiber
AC
suspensions, Int. J. Multiphas. Flow. 30 (2004) 67-87. [28] T. Papathanasiou, D. Guell, Flow induced alignment in composite materials, Woodhead Publishing, Cambridge, 1997.
30
ACCEPTED MANUSCRIPT
CR IP T
Figure Captions
AN US
Figure 1. Fiber modeled as a chain of spheres using (a) SBCM and (b) proposed model. In the proposed model, the indices at the upper and lower sides of the fiber are for spheres and segments, respectively. NR indicates the number of spheres
AC
CE
PT
ED
M
reduced per segment.
Figure 2. Two fibers considered in deriving the correction factors for (a) reduced fiber and (b) unreduced fiber. In contrast to the reduced fiber, spheres drawn with dotted lines exist and are aligned within each segment in the unreduced fiber. 31
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 3. Relative error for period of rotation of rigid fibers as function of number of
PT
segments. Open and filled symbols indicate that the correction process is or is
AC
CE
not applied, respectively.
32
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 4. Ratio of the calculation time of the SBCM to that of the proposed model as a
AC
CE
PT
function of number of segments.
33
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 5. Motion of flexible fibers when changing the number of segments: (a) time evolution of the z-coordinate of the sphere at the end of the fiber and (b) snapshots of fibers. 34
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
Figure 6. Orientation tensor component Axx at steady state for various fiber volume
CE
fractions Vf and aspect ratios rf. The results of the SBCM corresponds to
AC
Nseg=rf -1.
35
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
Figure 7. Average number of interactions per fiber at steady state for various fiber
AC
CE
volume fractions Vf and aspect ratios rf.
36
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
37
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
ED
Figure 8. Time evolution of orientation tensor components Axx and Ayy of rigid fiber suspensions: (a) rf=10, Vf=0.122; (b) rf=10, Vf=0.201; and (c) rf=20,
AC
CE
PT
Vf=0.056.
38
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
39
ACCEPTED MANUSCRIPT
Figure 9. Time evolution of orientation tensor components Axx and Ayy of flexible fiber
AN US
CR IP T
suspensions: (a) rf=20, Vf=0.056, and (b) rf=30, Vf=0.037.
Figure 10. Dimensions of the molded plaque model. Dashed squares mark regions
AC
CE
PT
ED
M
where fiber orientation tensors are calculated.
Figure 11. Flow velocity distributions obtained along the thickness direction at each 40
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
region during the mold-filling simulation.
41
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
42
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
ED
Figure 12. Fiber orientation tensor components Axx and Ayy along the thickness direction
AC
CE
PT
in (a) region A, (b) region B, and (c) region C.
43
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
Figure 13. Microstructure of fibers in region C. The fibers are modeled with 13
AC
CE
segments.
44
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure A.1. Capsule bounding volumes (dashed lines) of segments for detecting
CE
PT
ED
M
collisions.
AC
Figure A.2. Calculation procedure for repulsive and lubrication forces acting on each segment. In calculating lubrication forces, effects of imaginary spheres (illustrated with dotted lines) are considered.
45
Efficient Bead-Chain Model for Predicting Fiber Motion during Molding of Fiber-Reinforced Thermoplastics Toshiki Sasayama , Masahide Inagaki PII: DOI: Reference:
S0377-0257(18)30219-2 https://doi.org/10.1016/j.jnnfm.2018.10.008 JNNFM 4057
To appear in:
Journal of Non-Newtonian Fluid Mechanics
Received date: Revised date: Accepted date:
22 June 2018 4 September 2018 18 October 2018
Please cite this article as: Toshiki Sasayama , Masahide Inagaki , Efficient Bead-Chain Model for Predicting Fiber Motion during Molding of Fiber-Reinforced Thermoplastics, Journal of Non-Newtonian Fluid Mechanics (2018), doi: https://doi.org/10.1016/j.jnnfm.2018.10.008
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights An improved model for direct simulation of fibers in viscous flow is developed.
Fiber motions in simple shear flow are compared between the proposed and conventional model.
Fiber motions in injection molding are also simulated.
The improved model well reproduces the results of the conventional one. The model also provides higher computational efficiency than the conventional one.
AC
CE
PT
ED
M
AN US
CR IP T
1
ACCEPTED MANUSCRIPT
Efficient Bead-Chain Model for Predicting Fiber Motion during Molding of Fiber-Reinforced Thermoplastics
a
Toyota Central R&D Labs., Inc.
CR IP T
Toshiki Sasayamaa,*, Masahide Inagakia
* Corresponding Author (T. Sasayama) Toyota Central R&D Labs., Inc., 41-1, Yokomichi, Nagakute, Aichi 480-1192, JAPAN
AN US
Email address: [email protected] Tel: +81-561-71-7487 Abstract
This paper proposes an efficient model for simulating fiber motion in viscous flow,
M
especially during molding of fiber-reinforced thermoplastics. To increase the
ED
computational efficiency, the proposed model discretizes fibers as chains of fewer spheres than in the bead-chain model previously proposed by the authors, which is
PT
called the Simplified Bead-Chain Model (SBCM). In addition, in order to make the
CE
proposed model reproduce the results of the SBCM, correction factors applied to forces
AC
exerted on each sphere are newly derived so that the equation of motion of a fiber for the proposed model can be equivalent to that of the conventional model. Simulations of rigid and flexible fibers in simple shear flow are first performed to compare the fiber motions between the two models. The proposed model reproduces the results of the SBCM closely for cases of both a single fiber and suspensions containing multiple
2
ACCEPTED MANUSCRIPT
fibers. Finally, fiber orientation behavior during injection molding in a simple-shaped plaque is simulated. The results obtained with the proposed model shows good
CR IP T
agreement with those of the SBCM but at lower computational cost.
Keywords
AN US
Direct simulation; Fiber modeling; Bead-chain model;
M
Fiber-reinforced plastics;
AC
CE
PT
ED
Injection molding.
3
ACCEPTED MANUSCRIPT
1. Introduction To predict the mechanical and failure properties of injection- or compression-molded fiber-reinforced thermoplastics, it is essential to accurately predict
CR IP T
the distribution of orientations and lengths of fibers during the molding process.
Especially, the demand for predicting these properties of long-fiber composites is
AN US
increasing. Long fibers generally undergo large deformation due to flow resistance and interactions with other fibers during the molding processes, resulting in more complicated microstructures (e.g., fiber orientation and length) than those of short fiber
M
composites.
ED
Direct fiber simulations are useful for understanding the behavior of fibers during the molding process in detail (e.g., rheological properties [1-6], fiber orientation
PT
[7-9], and fiber breakage [10]). Concerning fiber orientation, in a practical injection
CE
molding simulation, statistical distributions of the fiber orientation are predicted using a
AC
macroscopic model describing the evolution of the fiber orientation tensor [11], in which the effect of fiber-fiber interactions is introduced as a diffusion term with some phenomenological parameters [12-14]. To investigate the validity of the macroscopic model, Mezher et al. [7] employed direct simulations and compared the obtained orientation tensor of fiber suspensions in simple shear flow with that predicted by the
4
ACCEPTED MANUSCRIPT
macroscopic model. As more practical applications, Yashiro et al. [8, 9] employed a moving particle semi-implicit method [15] for predicting behaviors of fiber orientation and accumulation during injection molding. Kuhn et al. [16, 17] simulated the
CR IP T
compression molding process of a ribbed plate to investigate the distribution of the fiber content within the rib.
AN US
Various models have been proposed for discretizing a fiber. Flexible fibers have been modeled as a chain of spheres [18-21], prolate spheroids [22], and rods
[23-25]. In our previous study [26], we developed an improved model, the Simplified
M
Bead-Chain Model (SBCM), using the work of Yamamoto and Matsuoka [18, 19] to
ED
enhance its computational efficiency by simplifying the governing equations. However, the SBCM has the drawback that the number of spheres required to model a fiber is
PT
equal to its aspect ratio. Therefore, the SBCM is computationally expensive for
CE
simulating long-fiber suspensions compared to the previous models [21-25], which use
AC
fewer components (such as rods and prolate spheroids) for discretizing fibers. It is thus desirable to enhance the computational efficiency of the SBCM by reducing the number of components needed to model a fiber to the levels of the models mentioned above. However, the mass and the moment of inertia of the fiber changes as the number of spheres decreases, which results in numerical errors in the fiber motion. Hence, it is also
5
ACCEPTED MANUSCRIPT
essential to establish a new method to decrease these errors. The purpose of this study is to develop a computationally efficient model that reproduces results as accurate as those of the SBCM. To this end, in the proposed model,
CR IP T
a fiber is modeled with fewer spheres than in the previous SBCM and the forces acting on each sphere are corrected by multiplying them with correction factors, which are
AN US
derived so that the equation of motion of a fiber in the proposed model becomes
equivalent to that in the SBCM. Fiber motions in a simple shear flow are then simulated and compared between the proposed model and the SBCM. Finally, simulations in an
PT
2. Numerical method
ED
model to complex flow.
M
injection molding process are carried out to investigate the applicability of the proposed
CE
2.1. Simplified Bead-Chain Model (SBCM)
AC
In the Bead-Chain Model (BCM) [18, 19], a fiber is modeled as chain of
spheres, as illustrated in Fig. 1(a), where the distance between the centers of two adjacent spheres is equal to the fiber diameter (i.e., closely packed). The SBCM solves only the equation of translational motion of each sphere, in contrast to the BCM, which solves both translational and rotational equations of motion. Thus, the governing
6
ACCEPTED MANUSCRIPT
equation of a sphere is given as
m
dvi Fi h FijS Fijb Fijht FinP , dt j 1 j 1 j 1 n 1
(1)
CR IP T
where m is the mass of the sphere, vi is the velocity of sphere i, the index j indicates the spheres connected to sphere i, the index n indicates the spheres that are not connected to sphere i, Fh is the hydrodynamic force, FS is the stretching force, FP is the
AN US
particle-particle interaction force, Fb is the bending force converted from the bending torque, and Fht is the force converted from the hydrodynamic torque. The SBCM
includes these converted forces for expressing the effect of bending and hydrodynamic
M
torques in the equation of translational motion as in Eq. (1), instead of the equation of
ED
rotational motion as in the BCM. Thus, the SBCM does not need to solve the equation
CE
paper [26].
PT
of rotational motion, and it ignores twisting fiber motion, as described in our previous
AC
2.2 Improvement of Simplified Bead-Chain Model To improve the efficiency of the SBCM, the number of spheres used for
modeling a fiber is reduced. In addition, a correction process is developed for eliminating the errors in fiber motions caused by reducing the number of spheres. First, we consider a fiber modeled with fewer spheres than in the SBCM, as 7
ACCEPTED MANUSCRIPT
shown in Fig. 1(b). Here, the fiber is considered to consist of Nseg segments, each of which is composed of two spheres and a massless rod connecting them. In addition, NR is defined as the number of reduced spheres in a segment (Fig. 1). The indexes of
CR IP T
spheres and segments are assigned sequentially, as indicated in Fig. 1(b). The indices at the upper side of the fiber in Fig. 1(b) are for spheres, and those in the lower side are for
AN US
segments.
Next, the correction procedure to reduce the numerical errors is described. The behavior of the fiber shown in Fig. 1(b) differs from that in Fig. 1(a) without any
M
corrections because the mass and moment of inertia are smaller. This paper proposes the
ED
procedure described below to correct the forces acting on spheres and eliminate these errors. For convenience, the fiber in Fig. 1(b) is referred to as a “reduced fiber” in this
PT
paper. Now, we consider a fiber made of closely packed spheres with the same
CE
conformation as the reduced fiber, as shown in Fig. 2. This paper refers to the fiber in
AC
Fig. 2(b) as an “unreduced fiber”. In each segment of the unreduced fiber, NR spheres (illustrated with dotted outlines in Fig. 2(b)) exist aligned. The unreduced fiber is an imaginary fiber used only for deriving correction factors for forces of spheres in the reduced fiber. The correction factors are derived so that the equations of motion for both fibers become equivalent. In constructing the equations of motion of these fibers, the
8
ACCEPTED MANUSCRIPT
following suppositions are made: (a) each segment is rigid, (b) lengths (i.e., NR) of the segments in each fiber are equal, and (c) distributed forces (e.g., hydrodynamic forces) vary linearly over a segment.
CR IP T
The forces acting on spheres are classified into two types, “distributed” and “non-distributed.” The term “non-distributed forces” means that these forces can be
AN US
regarded as acting only on the spheres at both ends of each segment. Let FD and FND be the distributed and the non-distributed forces, respectively. The distributed force FD includes the hydrodynamic force Fh. In contrast, the non-distributed force FND includes
M
the stretching force FS, the bending force Fb, and the force coming from the
ED
hydrodynamic torque Fht. Note that Fht can be approximated as non-distributed forces provided that the fluid velocity field is linear. The particle-particle interaction forces FP
PT
are assumed to be non-distributed forces in this study, although they are sometimes
AC
CE
distributed throughout the segment. In summary, FD and FND are given as
Fi ND FijS Fijb Fijht Fi Pns ,
(2)
Fi D Fi h .
(3)
j 1
j 1
j 1
ns 1
Note that in the fourth term on the right-hand side of Eq. (2), the index of FP changes from n to ns, which denotes the segment index. That is, interaction forces on sphere i are exerted by segment ns (where ns denotes a segment different from the segment that 9
ACCEPTED MANUSCRIPT
includes sphere i), in contrast to the SBCM. Here, the calculation procedure for interaction forces is briefly described; additional details can be found in Appendix A. In the proposed model, interaction forces are calculated between segments instead of
CR IP T
spheres, as in the SBCM. The calculated forces, which are acting on arbitrary points in the segments, are then converted into equivalent forces that are exerted on the two
AN US
spheres at both ends of each segment.
The equations of translational and rotational motions for each fiber are given in Appendix B. Using these equations, the correction factors for FD and FND are derived as
M
follows so that the equations for the reduced and unreduced fibers become equivalent.
ED
First, focusing on the equations of translational motion, each force exerted on sphere i is corrected as
PT
Fi ND T Fi ND ,
Fi D iT Fi D ,
(4)
CE
where F D and F ND are the corrected forces and T and T are the correction
AC
factors for the translational motion of reduced fibers given by
T
N seg 1 , N seg N R 1 1
(5)
N R 2 N seg 1 N seg N R 1 1 2 iT R seg N 1 N 1 N seg N R 1 1
i 1, N
seg
1
2 i N seg
10
.
(6)
ACCEPTED MANUSCRIPT
Next, focusing on the equations of rotational motion, the correction factors R and
R are obtained as follows: I Rf If
,
(7)
CR IP T
R
I Rf N R 2 N seg 1 3 N R 1 1 N R nˆ1 N 1 N 1 N 1 f seg R 6 N 1 N 1 I Rf 1 2 1 I iR f N R N R 2 nˆ i nˆ i 1 N i3 N i3 N i3 6 N R 1 R I 6 N 1 I Rf 1 N R 2 N seg 1 3 N R 1 N R nˆ N seg N 2 N 2 N 2 f R seg I 6 N 1 N 1
where N seg
i 1
2 i N , seg
i N
AN US
seg
N 1 N seg p 1 nˆ p , N seg
1
(9)
M
p 1
(8)
N pnˆ p ,
ED
2
(10)
p 1
AC
CE
PT
N i3 N i0
I Rf I
f
N
1 seg
1
N1,
(11)
i 1
N i0 nˆ p ,
(12)
p 1
N seg 1
3 rqG rqG q 1
N seg
,
3rNGseg 1 rNGseg 1 3 N R 1 rqG rqG 2aN R rqG nˆ q 2a 2 N seg N R N R 1 2 N R 1 q 1
(13) where riG is a vector from the center of gravity of the fiber to the position of sphere i 11
ACCEPTED MANUSCRIPT
and nˆ i is the unit vector from sphere i to sphere i+1. Throughout this correction process, the notation of this unit vector is different from that used in our previous paper [26] to simplify the equations. Specifically, the unit vector from sphere i to i+1 was
CR IP T
originally written as ni i 1 in [26], but now it is defined as nˆ i ni i 1 .
However, it is not easy to fully correct both the translational and rotational
AN US
motions of the fibers at the same time because the correction factors for each motion are different. This study thus employs the following procedure for correcting translational and rotational motions of the reduced fibers. Let niG be a unit vector from the center
M
of gravity of the fiber to the position of sphere i, and then FD and FND are decomposed
ED
into parallel and perpendicular components to niG , respectively. The correction factors for translational and rotational motions are multiplied to the parallel and the
PT
perpendicular components as
(14)
(15)
CE
Fi ND T Fi ND niG niG R Fi ND Fi ND niG niG ,
iT Fi D niG niG iR Fi D Fi D niG niG ,
AC
Fi D
so that at least the rotational motion of fibers is fully corrected. The calculation procedure of the reduced fiber described above can be
summarized as follows: (i) Calculate the correction factors using Eqs. (5)-(13).
12
ACCEPTED MANUSCRIPT
(ii) Calculate the distributed and non-distributed forces (FD and FND) using Eqs. (2) and (3). (iii) Correct the forces calculated in step (ii) using Eqs. (14) and (15).
m
dv i Fi D Fi ND . dt
(16)
AN US
Note that the model is reduced to the SBCM when NR=0.
CR IP T
(iv) Solve the following equation to update the coordinates and velocities of spheres:
3. Results and discussions
M
The effectiveness of the proposed model was investigated by comparing the
ED
fiber behavior results with those calculated with the SBCM. In this study, the number of segments Nseg was employed to describe the discretization of a fiber. A fiber of aspect
PT
ratio rf was discretized with Nseg=rf 1 when applying the SBCM, whereas 1 Nseg
CE
1 was used in the proposed model. In addition, this study assumed that the fiber and
AC
fluid motions are decoupled. The coupling between fibers and fluid is work for the future.
3.1. Fiber motion in simple shear flow Here, we describe investigation of the motion of fibers in simple shear flow.
13
ACCEPTED MANUSCRIPT
The velocity field of the flow at position x was given as V=( z, 0, 0), where the shear rate was set to =100 s-1. The radius of the fiber and the viscosity of the fluid were a=5 μm and η=100 Pa∙s, respectively. The density of fibers ρf was determined so that the
CR IP T
Reynolds number of a sphere (Re=a2 ρf /η) becomes Re=0.1 so that the effect of inertia can be ignored. Fibers were initially oriented parallel to the x axis (i.e., the flow
AN US
direction).
First, periods of rotation of a rigid fiber were compared. Here, particle-particle interaction forces were not calculated. The aspect ratio of the fiber was rf=100, and
M
Young’s modulus E was large enough to not induce bending deformation. Figure 3
ED
shows the relative error in periods of rotation calculated as 100×|T-TSBCM|/TSBCM, where TSBCM and T are the periods of rotation obtained by the SBCM (Nseg=99) and the
PT
proposed model (Nseg<99), respectively. The filled and open symbols indicate the results
CE
with and without, respectively, the correction process described in Section 2.2. When
AC
the correction process is not applied, the errors increase with decreasing number of segments. In contrast, these errors are mostly eliminated when the correction is applied. Figure 4 represents the ratio of computational time as a function of the number of segments. The calculations in Fig. 4 were performed with the same time step (Δt) for each number of segments. The computational efficiency is found to increase as the
14
ACCEPTED MANUSCRIPT
number of segments decreases. In this case, the calculation becomes about 43 times faster when the number of segments is reduced from 99 to 1. Next, the motion of flexible fibers was compared. Here, the aspect ratio and the
CR IP T
Young’s modulus were rf=50 and E=0.5 GPa, respectively. Here, calculating the
effective stiffness Seff (= E / 64 rf4 ) [27] is useful for evaluating the flexibility of a
AN US
fiber suspended in shear flow. In this case, the effective stiffness was about Seff
=3.9×10-4, which is low enough to cause bending deformation [23, 27]. Figure 5(a) shows the time evolutions of the z-coordinate of the sphere at the fiber end calculated
M
with Nseg = 49 (SBCM), 25, and 15. Figure 5(b) shows snapshots of the fibers at three
ED
time steps. As can be seen in Fig. 5(b), the fiber induces S-shaped deformation during rotation, as reported in our previous study [26]. The results match each other well for
PT
flexible fibers regardless of the number of segments. Of course, not even the proposed
CE
model could trace the deformation of flexible fibers if the number of segments is
AC
insufficient to resolve the deformed fiber shape. Finally, fiber suspensions were simulated considering particle-particle
interactions. The analysis region was a 500-μm cube in which the top and bottom planes in the z direction were wall boundaries and a periodic boundary condition was applied in the x and y directions. The simulations were performed for various fiber volume
15
ACCEPTED MANUSCRIPT
fractions Vf and aspect ratios rf. Fiber orientation tensors Aij [11] were calculated using the directional vectors of each segment. The simulations were run until the fiber orientation tensors almost reached a steady state. In all cases, the fiber suspensions are
CR IP T
homogeneously dispersed and show no heterogeneous structures (i.e., flocculation). This is because this study does not consider both the inter-fiber friction and the
AN US
non-straight equilibrium fiber shape, which are the factors controlling the formation of flocculation [23, 27]. Figure 6 shows the component of a fiber orientation tensor Axx at steady state as a function of Vf rf. Figure 7 plots the average number of interactions per
M
fiber at steady state. The number of interactions is defined as the total occurrences of
ED
lubrication and repulsive forces applied to spheres (for the SBCM) or segments (for the proposed model). Here, the product of the fiber volume fraction Vf and the aspect ratio rf
PT
can be used to estimate whether the semi-concentrated or concentrated regime
CE
dominates the suspension [28]. In Figs. 6 and 7, the black-filled symbols indicate the
AC
results obtained by the SBCM, whereas the grey-filled and open symbols are obtained by the proposed model. The proposed model gives results comparable to those of the SBCM in terms of both the orientation tensor (Fig. 6) and the number of interactions (Fig. 7). In addition, the difference in the orientation tensor between these models shown in Fig. 6 becomes smaller with increasing Vf rf. Figure 8 shows the transient 16
ACCEPTED MANUSCRIPT
behaviors of Axx and Ayy for various aspect ratios rf and volume fractions Vf. The horizontal axis in Fig. 8 is the applied shear strain defined as the product of the shear strain rate and time t. For rf=10 and Vf=0.122 (Fig. 8(a)), the initial behavior slightly
CR IP T
differs between the models depending on the number of segments. However, the differences in both the initial and steady-state regions decrease when Vf=0.201
AN US
(Fig. 8(b)). For rf=20 and Vf=0.056 (Fig. 8(c)), the curves mostly match, although Vf rf is less than that in the case of Fig. 8(a). It is also found that the difference between both models gets larger when Vf rf is less than or equal to unity (not show in the figure).
M
These differences in the orientation behavior between the two models that occur at Vf rf
ED
≤ 1 could be due to the translational motion of fibers not being fully corrected, unlike the rotational motion correction described in Section 2.2. However, these differences in
PT
the translational motion of fibers are considered to decrease when fibers experience a
CE
large number of interactions (i.e., in the concentrated regime). It is thus concluded that
AC
the proposed model reproduces the behavior of fiber suspensions simulated by the SBCM when Vf rf is higher than unity. Figure 9 shows the time evolution of Axx and Ayy in the case of flexible fibers with Young’s modulus E=0.1 GPa. The curves calculated with those models show good agreement in both cases of (rf, Vf)=(20, 0.056) and (30, 0.037).
17
ACCEPTED MANUSCRIPT
3.2. Fiber motion in injection molding process The proposed model was applied to the simulation of fiber motion in a
CR IP T
simple-shaped plaque during injection molding to compare the predicted fiber orientations with those of the SBCM.
AN US
First, a mold-filling simulation for the plaque illustrated in Fig. 10 was
performed using the commercial molding software 3D TIMON (Toray Engineering, Tokyo, Japan). The material was polypropylene with 20 wt% long glass fibers (L-2040P,
M
Prime Polymer, Tokyo, Japan). The inlet and mold temperatures were 230°C and 110°C,
ED
respectively. The flow rate at the inlet was 10 cm3/s. No holding pressure was applied in the mold-filling simulation. Figure 11 shows the flow velocity distribution in the
PT
thickness direction at each region (A, B, and C) delineated with dashed squares in Fig.
AC
inlet.
CE
10. The regions A, B, and C were located 15, 55, and 90 mm, respectively, from the
Then, a direct simulation of fibers was carried out using the transient flow field
obtained by the mold-filling simulation. The fiber radius a, length Lf, Young’s modulus E, and volume fraction Vf were a=10 μm, Lf=1 mm, E=70 GPa, and Vf=0.1, respectively. The fiber orientation at the inlet was set to be nearly aligned in the x (flow) direction.
18
ACCEPTED MANUSCRIPT
The number of segments in a fiber was set to Nseg=13 in the proposed model, whereas Nseg=49 in the SBCM. In calculating fiber-fiber interactions, only repulsive forces were considered to reduce the computational cost. Figure 12 shows the distribution of fiber
CR IP T
orientation tensor components Axx and Ayy along the z (thickness) direction in each
region. Figure 13 depicts fibers at region C simulated with the proposed model (where
AN US
only segments within the region of 2×2×1 mm are shown). In Figs. 12(b), 12(c), and 13, one can observe a layered structure, in which fibers in the shell layer are aligned in the flow direction, whereas those in the core layer show nearly in-plane random orientations.
M
The fibers in the core layer would become transversely aligned to the flow direction
ED
when using a mold shape that causes a diverging flow stronger than that simulated in this study. These results confirm that both models could capture the trends of fiber
PT
orientation that are seen in the injection molding process and that the proposed model
CE
works well even in a complicated flow field. Compared to the SBCM, the proposed
AC
model enables us to use a time step about 3.7 times larger, which yields a calculation about 2.9 times faster. Although the degree of speed-up is smaller than that for single-fiber calculations due to the existence of fiber-fiber interactions, which are computationally expensive, the proposed model succeeds in reducing the computational cost, even when considering these interactions.
19
ACCEPTED MANUSCRIPT
4. Conclusions This study proposed an efficient model based on the SBCM by reducing the
CR IP T
number of spheres needed to model a fiber and correcting the forces exerted on each sphere to reduce the computational cost without loss of accuracy compared to the
AN US
SBCM.
The correction process was applied to simulations of fiber motion in simple shear flow to investigate the accuracy of the proposed model in comparison to the
M
SBCM. First, periods of rotation of single rigid fibers and deformation of single flexible
ED
fibers were compared between these models. The proposed model mostly matched the results of the SBCM, in contrast to the large differences that occurred without the new
PT
correction process. In addition, it was found that the computational efficiency increased
CE
as the number of segments decreased. Second, suspensions of both rigid and flexible
AC
fibers were simulated with consideration of fiber-fiber interaction. The fiber orientation tensors of the suspensions predicted by the proposed model and the SBCM showed good agreement provided that the suspensions were in the concentrated regime (i.e., the product of the fiber volume fraction and the aspect ratio Vf rf is greater than unity). Finally, orientation behavior of flexible fibers within a simple-shaped plaque
20
ACCEPTED MANUSCRIPT
during injection molding was compared between the models. Both models predicted orientations with layered structure that are typically observed in injection molding, and the proposed model mostly reproduced the results obtained with the SBCM.
CR IP T
The validity and the effectiveness of the proposed model were confirmed by
the investigations in this study. Use of the proposed model will produce more realistic
AN US
simulations of the motion of long fibers in a highly concentrated state during molding of fiber-reinforced thermoplastics. Our future work will focus on simulating more practical molding processes by considering fiber breakage phenomena and two-way coupling
ED
M
between fibers and fluid.
Appendix A. Calculation procedure for particle-particle interaction forces
PT
This appendix describes the procedure for calculating particle-particle
CE
interaction forces (repulsion and lubrication) in the proposed model. First, segments are
AC
regarded as capsules for detecting collisions between them (Fig. A.1). When a fiber consists of one segment, a capsule bounding volume is set as illustrated in Fig. A.1(a) with dashed lines. In other cases, segments in the fiber have capsule bounding volumes as shown in Fig. A.1(b) with dashed lines. With these assumptions, repulsive and lubrication forces are calculated as described below.
21
ACCEPTED MANUSCRIPT
(i) Repulsive forces Consider two closely located segments indexed as is and ns, as shown in
force acting on the segment is is calculated by x c xisc 1 ns FisRseg D exp G ns 0 0 2a
x c x c is ns c xns xisc
CR IP T
Fig. A.2(a), where the segment is consists of spheres indexed as i and j. The repulsive
c ( xns xisc 2.001a) ,
(A.1)
AN US
c where D0 and G0 are constants and xisc and xns are the coordinates of the closest
points of the two segments is and ns, respectively. The repulsive forces acting on spheres i and j are determined as
respectively.
x j xi and xi and x j are the coordinates of spheres i and j,
PT
where s xisc xi
(A.2)
ED
M
Rep Rseg Fi ns 1 s Fis ns , Rep Rseg F s F j ns is ns
CE
(ii) Lubrication forces
AC
When calculating the lubrication forces between two segments, this study
accounts for lubrication effects caused by imaginary spheres (illustrated in Fig. A.2(b) with dotted lines), as well as the two spheres at both ends of each segment. Let p and q be the indexes of spheres in segments is and ns, respectively. The lubrication force applied to sphere p by sphere q is given by
22
ACCEPTED MANUSCRIPT
FpqLseg
x xp 3 2 vq v p a q 2 x q x p 2a x q x p
xq x p xq x p
(2.001a xq x p 3a) ,
(A.3)
where η is the viscosity of the fluid. When p (or q) is an imaginary sphere, the velocity
CR IP T
vp of the sphere is linearly interpolated using the velocities of the spheres at both ends of the segment. The lubrication forces acting on the spheres (i and j) of the segment is are calculated as
where s p x p xi
AN US
Lseg Fi Lub ns 1 s p Fpq p q , Lub Lseg F s F j ns p pq p q
x j xi .
(A.2)
M
Finally, the particle-particle interaction forces of the spheres (i and j)
ED
exerted by the segment ns are given by
(A.3)
CE
PT
P Rep Lub Fi ns Fi ns Fi ns . P Rep Lub F F F j ns j ns j ns
AC
Appendix B. Equations of fiber motion The equations of translational motion for the center of mass of the reduced
(Fig. 2(a)) and unreduced (Fig. 2(b)) fibers are given by
N
seg
1 m
dv Rf dt
N seg 1
i 1
T
Fi
ND
N seg 1
i 1
T i
Fi D ,
23
(B.1)
ACCEPTED MANUSCRIPT
N N seg
R
1 1 m
dv f dt
N seg 1
F i 1
ND
i
NR 2 D N 1 Fi F1 FNDseg 1 , 2 i 2
N seg 1
R
D
(B.2)
where vRf and v f are the velocity of the center of mass of the reduce and unreduced
CR IP T
fibers, respectively, and Nseg(NR+1)+1 in the left-hand side of Eq. (B.2) is the number of spheres in the unreduced fiber, including the imaginary spheres illustrated in Fig. 2(b) with dotted lines.
AN US
To write the equation of rotational motion for the center of mass of all fibers, we use the body-fixed coordinate system, in which the origin is fixed in the center of mass of the fibers and the axes are along the principal axis of inertia. The equation of
M
rotational motion of the reduced fiber is written as
~ Rf dI~ Rf N seg 1 ~ Rf dω ~ Rf ~ Rf ~ ~ Rf Rf ~ ~ I ω ω I ω ~ ri G R Fi ND iR Fi D , dt dt i 1
ED
(B.3)
PT
~ ~ Rf are the inertia tensor and the angular velocity of the reduced where I Rf and ω
CE
fiber, respectively, and the tilde ~ means that the quantity is written in the body-fixed
AC
coordinate system. Assuming that the fiber maintains the same shape during time step Δt, the second term in the left-hand side of Eq. (B.3) can be ignored, and then the equation becomes
~ Rf N 1 ~ Rf dω ~ ~ ~ Rf I~ Rf ω ~ Rf ~ I ω ri G R Fi ND iR Fi D . dt i 1
seg
The right-hand side of Eq. (B.4) is calculated as 24
(B.4)
ACCEPTED MANUSCRIPT
N seg 1
N~
~ ~ r~i G R Fi ND iR Fi D 2 N R 1 a i 1
2N 1a ~ N
N seg 1 i 1
0 i
1 N
seg
~ ~ N 1 R Fi ND 1
N 1 2 N 1 a ~1 1 ~ ~ ~0 R ~D R N F 2 N 1 a N 1 iR Fi D . N i seg i 1 seg N 1 N 1 i 2 R
R
N
seg
2
1
~ iR FNDseg 1
seg
(B.5)
diagonal and described as
0 ~ Rf I yy 0
0 0 , ~ I zzRf
(B.6)
AN US
~ I Rf
~ I xxRf 0 0
CR IP T
Here, since the coordinates are along the principal axes of inertia, the inertia tensor is
~ ~ ~ where I xxRf , I yyRf , and I zzRf are the principal moments of inertia. Likewise, the
M
equation of rotational motion of the unreduced fiber is derived as
~f N seg N R 11 ~ f dω ~f ~f ~ ~ f ~ I ω I ω ~ ri G Fi ND Fi D , dt i 1
(B.7)
ED
PT
~ where I f is the inertia tensor of the unreduced fiber. The right-hand side of Eq. (B.7) is calculated as
AC
i 1
N 1 1 ~ ~ ~ ~ ~ r~i G Fi ND Fi D 2 N R 1 a N i0 seg N 1 Fi ND N 1 i 1 R R N 2a 3 N 1 ~1 ~ D R~ N F1 N nˆ1 seg 3 N 1
CE
N seg N R 1 1
seg
~ ~ 2 ~ a N 1 1 ~ ~ N R N R 2 nˆ i nˆ i 1 6 N R 1 N i0 seg N 1 Fi D 3 i 2 N 1
seg
N
R
.
2a 3 N R 1 ~ 2 ~ D R~ ˆ N FN seg 1 N nN seg seg 3 N 1 (B.8)
25
ACCEPTED MANUSCRIPT
To make the motions of the reduced and unreduced fibers equivalent, the following relations need to be satisfied:
~ Rf ω ~f . ω
dt
~f dω dt
.
(B.10)
CR IP T
~ Rf dω
(B.9)
Here, we assume that the following holds for the ratios of principal moments of inertia between the reduced and unreduced fibers for each component:
0 ~ Rf ~ f I yy I yy 0
M
where δ is a unit tensor and
I
Rf
m
N seg 1
~r q 1
ED
Rf I f δ, ~f I I zz
AN US
~ ~ I xxRf I xxf ~ ~ 1 I Rf I f 0 0
G q
~ I zzRf
0 0
(B.11)
~ rqG ,
(B.12)
N N ~ ~ I m ~ rqG ~ rqG m ~ rqG 2apnˆ q ~ rqG 2apnˆ q . i 1 q 1 p 1 N seg 1
seg
R
(B.13)
PT
f
CE
The second term in the right-hand side of Eq. (B.13) is for the imaginary spheres
AC
illustrated in Fig. 2(b). In addition, the right-hand side of Eq. (B.13) is rearranged when it appears in Eq. (13). Finally, by substituting Eqs. (B.4) and (B.7) into Eq. (B10) and using
Eqs. (B.9) and (B.11), the resulting relation for deriving the correction factors αR and βR is given as
26
ACCEPTED MANUSCRIPT
I Rf G R ~ ND R ~D ~ r F F i i i i If i 1
N seg 1
N seg N R 1 1
i 1
~ ~ ~ ri G Fi ND Fi D .
(B.14)
CR IP T
References [1] L. H. Switzer III, D. J. Klingenberg, Rheology of sheared flexible fiber suspensions via fiber-level simulations, J. Rheol. 47 (2003) 759-778.
AN US
[2] M. Yamanoi, J. Maia, Analysis of rheological properties of fibre suspensions in a Newtonian fluid by direct fibre simulation. Part 1: Rigid fibre suspensions, J. Nonnewton. Fluid. Mech. 165 (2010) 1055-1063.
M
[3] M. Yamanoi, J. Maia, T. Kwak, Analysis of rheological properties of fibre
ED
suspensions in a Newtonian fluid by direct fibre simulation. Part 2: Flexible fibre
PT
suspensions, J. Nonnewton. Fluid. Mech. 165 (2010) 1064-1071.
CE
[4] J. Andrić, S. B. Lindström, S. Sasic, H. Nilsson, Rheological properties of dilute suspensions of rigid and flexible fibers, J. Nonnewton. Fluid. Mech. 212 (2014) 36-46.
AC
[5] R. Mezher, E. Abisset-Chavanne, J. Férec, G Ausias, F. Chinesta, Direct simulation of concentrated fiber suspensions subjected to bending effects, Modelling Simul. Mater. Sci. Eng. 23 (2015) 055007. [6] C. G. Joung, N. Phan-Thien, X.J. Fan, Viscosity of curved fibers in suspension, J. Nonnewton. Fluid. Mech. 102 (2002) 1-17. 27
ACCEPTED MANUSCRIPT
[7] R. Mezher, M. Perez, A. Scheuer, E. Abisset-Chavanne, F. Chinesta, R. Keunings, Analysis of the Folgar & Tucker model for concentrated fibre suspensions in unconfined and confined shear flows via direct numerical simulation, Composites Part
CR IP T
A. 91 (2016) 388-397.
[8] S. Yashiro, T. Okabe, K. Matsushima, A numerical approach for injection molding of
AN US
short-fiber-reinforced plastics using a particle method, Adv. Compos. Mater. 20 (2011) 503-517.
[9] S. Yashiro, H. Sasaki, Y. Sakaida, Particle simulation for predicting fiber motion in
M
injection molding of short-fiber-reinforced composites, Composites Part A. 43 (2012)
ED
1754-1764.
[10] S. Yamamoto, T. Matsuoka, Dynamic simulation of flow-induced fiber fracture,
PT
Polym. Eng. Sci. 35 (1995) 1022-1030.
CE
[11] S. G. Advani, C. L. Tucker III, The use of tensors to describe and predict fiber
AC
orientation in short fiber composites, J. Rheol. 31 (1987) 751–784. [12] F. Folgar, C. L. Tucker III, Orientation behavior of fibers in concentrated suspensions, J. Reinf. Plast. Compos. 3 (1984) 98-119. [13] J. H. Phelps, C. L. Tucker III, An anisotropic rotary diffusion model fort fiber orientation in short- and long-fiber thermoplastics, J. Nonnewton. Fluid. Mech. 156
28
ACCEPTED MANUSCRIPT
(2009) 165-176. [14] H. C. Tseng, R. Y. Chang, C. H. Hsu, An objective tensor to predict anisotropic fiber orientation in concentrated suspensions, J. Rheol. 60 (2016) 215-224.
incompressible fluid, Nucl. Sci. Eng. 123 (1996) 421-434.
CR IP T
[15] S. Koshizuka, Y. Oka, Moving-particle semi-implicit method for fragmentation of
AN US
[16] C. Kuhn, I. Walter, O. Taeger, T. A. Osswald, Experimental and numerical analysis of fiber matrix separation during compression molding of long fiber reinforced thermoplastics, J. Compos. Sci. 1(1) (2017) 2.
M
[17] C. Kuhn, I. Walter, O. Taeger, T. A. Osswald, Simulative prediction of fiber-matrix
ED
separation in rib filling during compression molding using a direct fiber simulation, J. Compos. Sci. 2(1) (2018) 2.
PT
[18] S. Yamamoto, T. Matsuoka, A method for dynamic simulation of rigid and flexible
CE
fibers in a flow field, J. Chem. Phys. 98 (1993) 644-650.
AC
[19] S. Yamamoto, T. Matsuoka, Dynamic simulation of fiber suspensions in shear flow, J. Chem. Phys. 102 (1995) 2254-2260. [20] C. G. Joung, N. Phan-Thien, X.J. Fan, Direct simulation of flexible fibers, J. Nonnewton. Fluid. Mech. 99 (2001) 1-36. [21] H. F. Guo, B. G. Xu, A 3D numerical model for a flexible fiber motion in
29
ACCEPTED MANUSCRIPT
compressible swirling airflow, Comput. Model. Eng. Sci. 61 (2010) 201-222. [22] R. F. Ross, D. J. Klingenberg, Dynamic simulation of flexible fibers composed of linked rigid bodies, J. Chem. Phys. 106 (1997) 2949-2960.
CR IP T
[23] C. F. Schmid, L. H. Switzer, D. J. Klingenberg, Simulations of fiber flocculation: Effects of fiber properties and interfiber friction, J. Rheol. 44 (2000) 781-809.
AN US
[24] G. Wang, W. Yu, C. Zhou, Optimization of the rod chain model to simulate the
motions of a long flexible fiber in simple shear flows, Eur. J. Mech. Fluids/B. 25 (2006) 337-347.
M
[25] S. B. Lindstöm, T. Uesaka, Simulation of the motion of flexible fibers in viscous
ED
fluid flow, Phys. Fluids. 19 (2007) 113307.
[26] T. Sasayama, M. Inagaki, Simplified bead-chain model for direct fiber simulation
PT
in viscous flow, J. Nonnewton. Fluid. Mech. 250 (2017) 52-58.
CE
[27] L. H. Switzer, D. J. Klingenberg, Flocculation in simulations of sheared fiber
AC
suspensions, Int. J. Multiphas. Flow. 30 (2004) 67-87. [28] T. Papathanasiou, D. Guell, Flow induced alignment in composite materials, Woodhead Publishing, Cambridge, 1997.
30
ACCEPTED MANUSCRIPT
CR IP T
Figure Captions
AN US
Figure 1. Fiber modeled as a chain of spheres using (a) SBCM and (b) proposed model. In the proposed model, the indices at the upper and lower sides of the fiber are for spheres and segments, respectively. NR indicates the number of spheres
AC
CE
PT
ED
M
reduced per segment.
Figure 2. Two fibers considered in deriving the correction factors for (a) reduced fiber and (b) unreduced fiber. In contrast to the reduced fiber, spheres drawn with dotted lines exist and are aligned within each segment in the unreduced fiber. 31
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 3. Relative error for period of rotation of rigid fibers as function of number of
PT
segments. Open and filled symbols indicate that the correction process is or is
AC
CE
not applied, respectively.
32
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 4. Ratio of the calculation time of the SBCM to that of the proposed model as a
AC
CE
PT
function of number of segments.
33
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 5. Motion of flexible fibers when changing the number of segments: (a) time evolution of the z-coordinate of the sphere at the end of the fiber and (b) snapshots of fibers. 34
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
Figure 6. Orientation tensor component Axx at steady state for various fiber volume
CE
fractions Vf and aspect ratios rf. The results of the SBCM corresponds to
AC
Nseg=rf -1.
35
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
Figure 7. Average number of interactions per fiber at steady state for various fiber
AC
CE
volume fractions Vf and aspect ratios rf.
36
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
37
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
ED
Figure 8. Time evolution of orientation tensor components Axx and Ayy of rigid fiber suspensions: (a) rf=10, Vf=0.122; (b) rf=10, Vf=0.201; and (c) rf=20,
AC
CE
PT
Vf=0.056.
38
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
39
ACCEPTED MANUSCRIPT
Figure 9. Time evolution of orientation tensor components Axx and Ayy of flexible fiber
AN US
CR IP T
suspensions: (a) rf=20, Vf=0.056, and (b) rf=30, Vf=0.037.
Figure 10. Dimensions of the molded plaque model. Dashed squares mark regions
AC
CE
PT
ED
M
where fiber orientation tensors are calculated.
Figure 11. Flow velocity distributions obtained along the thickness direction at each 40
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
region during the mold-filling simulation.
41
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
42
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
ED
Figure 12. Fiber orientation tensor components Axx and Ayy along the thickness direction
AC
CE
PT
in (a) region A, (b) region B, and (c) region C.
43
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
Figure 13. Microstructure of fibers in region C. The fibers are modeled with 13
AC
CE
segments.
44
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure A.1. Capsule bounding volumes (dashed lines) of segments for detecting
CE
PT
ED
M
collisions.
AC
Figure A.2. Calculation procedure for repulsive and lubrication forces acting on each segment. In calculating lubrication forces, effects of imaginary spheres (illustrated with dotted lines) are considered.
45