Particle Simulation for Predicting Effective Properties of Short Fiber Reinforced Composites

We present a method for evaluating elastic properties of a composite material produced by molding a resin filled with short elastic fibers. A flow of the filled resin is simulated numerically using a mesh-free method. After that, assuming that spatial distribution and orientation of fibers are not significantly changed during polymerization, effective elastic moduli of the composite material are evaluated. The developed micro-mechanical mathematical modelling of effective moduli is aimed to molding process optimization which results in product quality improvement.


Introduction
The subject of the study is a composite material reinforced by short fibers and produced by molding and subsequent polymerization.The goal of developed simulation is to predict elastic moduli of the solid composite material.The simulation includes a flow model of the filled resin, a technique for the solid composite effective moduli computation, and numerical implementation.
In Sec. 2, we propose an advanced mesh-free method for a numerical simulation of a fiber-filled resin flow.The output of this simulation includes fiber volume ratios and orientations for every part of the flow.
Fibers have a complex spatial structure which is continually varying due to the action of the resin flow, which significantly complicates mathematical modeling of the reinforced fluid dynamics.There are two main approaches for studying this dynamics.The first approach [Dinh and Armstrong, 1984;Folgar and Tucker III, 1984;Ranganathan and Advani, 1991;Tucker III, 1991;Iso et al., 1996] consists in applying an averaged model.For example, an empirical model was considered in Chung and Kwon [1995].This model used orientation tensors first introduced in Advani and Tucker III [1987].Within the first approach, many authors used the Folgar-Tucker equation and its modifications for computation of macro-mechanical parameters of the reinforced fluid.In Fan et al. [1999] a hypothesis that fibers intersection is a probabilistic event characterized by a normally distributed random variable was adopted.
All computational methods based on the fiber direction distribution function are similar in their relative computational simplicity and ability to evaluate mean stresses of the fiber-fluid mixture, elastic moduli and fiber directions.Nevertheless in this paper we use the second approach utilizing direct simulation of the individual fiber trajectories which originates from Yamamoto and Matsuoka [1996] and Yamamoto and Matsuoka [1999].It is supposed that this approach can describe the behavior of fibers more adequately to the problem of determining effective elastic properties of a composite material.The particle-based technique allows computing motion of any individual fiber.It is capable to describe heterogeneous distribution of fibers, interaction between individual fibers, fiber breaking and jamming.Therefore, the developed simulation can be used for accurate prediction of material properties.From this point of view, our method is a logical extension of such papers as Advani and Tucker III [1987], Chung and Kwon [1995], Yamamoto and Matsuoka [1996], Fan et al. [1999] and Yamamoto and Matsuoka [1999].We not only use mathematical modelling for prediction of fiber spatial distribution at various time moments but basing on this distribution we also evaluate certain properties of the produced composite thus providing a complete combined approach.
A combined approach for predicting elastic and thermal parameters was developed, e.g., in Gupta and Wang [1993].The authors used Halpin-Tsai equations [Halpin and Kardos, 1976] for effective moduli computation.Fiber orientation state required for this computation was represented in terms of a second-order orientation tensor which in turn was calculated on the base of evolution equation.The latter equation employed fluid velocity field and Advani and Tucker closure formula for fourth-order orientation tensor; however, fiber distribution influence on velocity field was not discussed.A similar fact holds for a commercial software Moldex3D.On contrary, the approach developed in this paper provides strong coupling between fluid velocity field and density and orientation of fibers.
Thus, the main result of the paper is an advanced mesh-free method for direct numerical simulation of the flow of a viscous fluid filled with short inclusions.
An example of simulation of a two-dimensional flow by means of a mesh-free particle method was presented, e.g., in Yashiro et al. [2012].
Mesh-free methods became known after the paper by Gingold and Monaghan [1977] where astrophysical processes were simulated using the method called smoothed particle hydrodynamics (SPH).The methodology of SPH was extended later to simulating fluid and gas flows, large deformations of solids, mixing of heterogeneous materials, destruction and propagation of cracks and various other processes [Liu, 2009].In this paper, we apply the mesh-free method to simulation of large deformations of filled fluid with a viscous liquid matrix and high (above 10%) volume ratio of elastic inclusions.The idea is in approximation of the motion equations by a system of differential equations describing the motion of all particles.The forces acting on particles are computed from the Navier-Stokes equations.
The developed method is different from a typical weakly compressible SPH (WCSPH) method as it uses a more accurate field smoothing procedure.An attractive feature of the method is the fact that it belongs to the Lagrangian type, and at the same time it copes with deformations larger than those that are typical for total Lagrangian approach in the Finite Element Methods.Until recently, mesh-free methods and, in particular, SPH-type methods lost to Lagrange or Euler-Lagrange finite element methods with efficient remeshing technique in accuracy of results and computational cost in simulations of large deformations.The situation changed with the advent of readily available parallel computational hardware such as graphic accelerators.Therefore for filled resin molding simulation we developed a parallel implementation of the introduced mesh-free method.This implementation can be run on CPU or GPU, and in the latter case one obtains a significant performance improvement.
In Sec. 3, we describe a method capable of evaluating effective elastic moduli of a solid composite material with accuracy adequate to experimental data dispersion.We consider both fibers with a uniform spatial distribution and fiber with distributions given by an arbitrary probability density.In the first case, closed-form formulas for effective elastic moduli are well known.Having compared the values obtained by these formulas with experimental data (Sec.4), we came to conclusion that these formulas provide an accuracy that is comparable to the experimental data dispersion.
In order to compute effective moduli of elasticity in the case of arbitrary fiber distribution, we use two distribution functions, namely, spatial distribution of fiber volume ratio and spatial distribution of fiber direction.Both distribution functions are used in discretized form.
In Sec. 4, we give numerical examples of applying the developed simulation method and compare simulation results with theoretical solutions or experimental data.This comparison validates the method.Thus, we propose a complete combined approach allowing one to determine elastic moduli of a solid composite material filled with short fibers.The proposed approach consists of two stages.First, we simulate the resin-with-fibers flow and determine spatial distributions of fiber volume ratio and fiber orientation at the end of the molding process.The simulation is based on a mesh-free technique.At the second stage, we compute the effective (averaged) moduli of elasticity of the polymerized material.The concentration and orientation of fibers, as well as the moduli of elasticity, are the output values that may be used for the molding process optimization.

Simulation of the Fiber Filled Resin Flow
In this section, we present a mesh-free method developed for determination of the density and orientation of fibers appearing in the process of molding.
In general, the developed method belongs to WCSPH family.However, an essential distinction of the developed method from the classic SPH is a substantially different approximation of the velocity field.The classic SPH uses approximations obtained by quadrature formulas applied to convolutions of the considered fields with the smoothing functions.This leads to the fact that the simulation of a viscous flow by the SPH method in regions with a low concentration of particles results in "parasitic" viscous friction.The source of this friction is inaccuracy of the applied quadrature formulas: classic variants of SPH use one-point quadratures.The method used in this paper belongs to the class of so-called partition of unity methods [Liu, 2009].Our method utilizes the well-known Shepard approximation [Shepard, 1968] to the velocity field, which provides an adequate behavior with a small computational cost.At the same time, the incompressibility condition is treated in the same way as in the classic SPH.
The continuous molding model which we use has the following properties.The filled resin is assumed to be a viscous incompressible fluid.The elastic inclusions are simulated as one-dimensional elastic objects possessing extension-compression and bending stiffnesses.
We describe the numerical implementation of the discrete model of the molding process and start with the discretization of the Navier-Stokes and incompressibility equations used for description of the fluid motion.To discretize the incompressibility equation (2.1b), we used the approach typical to weakly compressible particle methods: no incompressibility was forced, instead the special function was used which penalized large deviations from the incompressible flow.
Following the traditional ideology of particle methods, we represent a continuous medium by a finite set of particles.The ith particle is represented by a ball B i of radius ε with the center at some point r i , possessing density ρ i , viscosity µ i , velocity υ i , and also smoothing functions ϕ i (r) and ψ i (r).The ball B i is an "effective neighborhood of the ith particle", i.e., we simulate only the interaction of particles positioned from each other at the distance not exceeding ε.The smoothing functions ϕ i are constructed by shifting and stretching the standard function ϕ: The functions ψ i form Shepard's partition of unity [Shepard, 1968] constructed on the base of the functions ϕ i , i.e., . (2.3) We set the function ϕ equal to the spherically symmetric polynomial of the fourth degree relative to the distance from the origin: The set of particles is associated with the two fields.The first one is a field of smoothed number of particles and the second field is a field of smoothed velocities The field of smoothed number of particles is proportional to the actual density of particles in the discrete medium, the field of smoothed velocities has the meaning of mean velocity.Now, we can formulate an approximation to the Navier-Stokes equation (2.1a) in terms of the introduced variables: (2.7) where P is the function used to provide approximate incompressibility of the discrete medium, i.e., approximate fulfillment of Eq. (2.1b).In the same way as in WCSPH family of methods, penalty function P is chosen as P (n) = K n, where K > 0 is a penalty parameter.The optimal value for K depends on the discretization parameter ε and the time step.We adjusted the value of K empirically utilizing the following approach.We considered Taylor expansion of the state equation p(ρ) of the fluid in terms of (ρ − ρ 0 )/ρ 0 , and selected the penalty parameter K close to the product Cρ 0 , where C was a coefficient in the linear approximation of this expansion.
Following the task of developing a molding simulation, we assumed that the viscosities and densities of particles are the same, i.e., µ i = µ and ρ i = ρ.Clearly, this restriction is not immanent for the developed particle method.
Next, we explain approximations of boundary conditions that were used in numerical examples described below.Suppose the whole boundary Σ consists of two parts.The part Σ 2 is where the pressure p 0 is given.For the boundary part Σ 1 no-flux or no-slip conditions are posed.
Boundary conditions at rigid boundaries can be interpreted in particle methods as collisions between particles and the boundary.Here, we describe the particular model of these collisions used in our simulations (which corresponds to nonelastic collision of classic rigid-body mechanics).Let υ τ be the projection of particle velocity onto the plane tangent to the boundary at the point where the collision happens (i.e.υ τ is the tangent velocity component), and let υ n be the normal component.In this notation, the velocity of particle is υ τ + υ n .The velocity after the collision is set to e τ υ τ − e n υ n , where e τ and e n are parameters.Setting e τ = e n = 0 corresponds to the standard no-slip condition and setting e τ = e n = 1 corresponds to the standard no-flux condition.
For Σ 1 we impose no-slip condition simply setting e τ = e n = 0.As far as Σ 2 is concerned, in the numerical tests described below this boundary is the surface of a piston creating pressure p 0 .This approach is quite general and can be used to simulate various high-pressure injection processes by placing virtual pistons in each injector.The piston is simulated by a moving planar surface.Its velocity at each time step is selected in such a way that the total variation of the particles momentum caused by collisions with the piston results in pressure p 0 .Such a selection of velocity can be made by a multitude of methods.We use a classic binary partition method: starting with the lower and the upper bounds for the velocity value we gradually improve these bounds by evaluating the total particle momentum variation caused by collisions with the piston for the piston velocity equal to the half-sum of the bounds.This approach allows us to simulate collisions with the piston in the same way as we simulate the condition on boundary Σ 1 but with parameters e τ = e n = 1.These conditions are applicable both to fluid and inclusion particles.Now, we proceed to the description of the model of an elastic fiber.An inclusion is modeled by a set of N point masses m i linked with each other by onedimensional elastic elements of length l and rigidity k (the total length of inclusion equals L = (N −1)l).Elastic elements simulate the rigidity of inclusion in extensioncompression.In addition, we simulate the rigidity in bending using a simple linear penalty as follows.If a, b and c are the radius vectors of three adjacent point masses of the fiber, then the force acts on the middle particle c, and halves of this force with opposite directions act on the other two particles a and b.It is important to note that point masses are the part of fluid particles so the no-slip condition between viscous matrix and elastic is satisfied automatically.
If we multiply discrete Navier-Stokes Eq. (2.7) by the volume V of a fluid particle, then the right-hand side of this equation gives forces acting on the particles.If we additionally take into account the masses of inclusions in the left-hand side of the equation and the forces caused by elastic expansion-compression and bending deformations in the right-hand side, then we get the equation (2.9) Here, m i equals the mass of the segment of a rigid inclusion if the ith particle belongs to the inclusion, and equals zero if the ith particle corresponds to fluid.
The force F i is the elastic force acting on the ith particle from neighboring mass points in inclusion.Introducing the notation f i = F i /V, below we use the following form of this equation: (2.10) After full expansion and projection onto the kth axis this equation takes the form where the partial derivatives ∂ k ϕ are calculated analytically using the explicit form (2.4) of the smoothing function ϕ.Equations (2.10) and (2.11) describe the motion of any particle regardless of its type.These equations were derived using spatial discretization and form the system of ordinary differential equations (ODEs).
Certainly, there are a lot of solution methods for systems of ODEs.We utilize the fourth-order Runge-Kutta method which proved to be sufficiently efficient according to our tests.Using notation u i = (r i , υ i ), i = 1, 2, . . ., N, and u = (u 1 , u 2 , . ..) we rewrite Eqs.(2.10) together with the equations ṙi = υ i in the form u = A(u) (2.12) and discretize this system in time by the classic fourth-order Runge-Kutta scheme (RK4): (2.13) Numerical experiments showed that for sufficiently small time step τ there is no significant difference between RK4 and popular schemes of lower order [Butcher, 2003] such as RK2, leapfrog, explicit and semi-implicit Euler methods, etc.Nevertheless, in comparison with schemes of lower order, RK4 allows one to take the time step 10-20 times greater and obtain the same accuracy without getting an unstable behavior.
Computation of the right-hand side A(u) of Eq. (2.12) requires explanation.The expression for the right-hand side for the ith particle includes sums over all particles.However, the summation over all possible pairs at each time step is computationally infeasible for a number of particles required for more or less practical simulations: according to our experience, in most cases this number exceeds 100,000, and it is never less than 10,000.On the other hand, it is easy to notice that these sums can be taken only over particles that are located at the distance not exceeding ε from the ith particle because the other summands in the sums are equal to zero.Therefore, an efficient method for the determination of the "neighboring" particles is critical for the computations.
In the simulations of a large number of interacting particles, the most popular neighbor search method is so-called spatial partitioning [Ihmsen et al., 2011].It is described as follows.The whole space is split into cubic cells with sides αε and the cells are enumerated in some way, for example, by triples of Cartesian coordinates on an integer-valued grid or by a z-type enumeration [Ihmsen et al., 2011].Then, access to particles is performed using a data structure that associates the set of particles located into a cell with the number of this cell.In our simulation, we enumerate the cells in natural order (lexicographically with respect to their coordinates) by consecutive natural numbers, and store references to sets of neighbors in an array.The parameter α was adjusted by test computations.If α is large, cells contain too many particles, and if α is small, the search for neighbors tests too many cells.Our test identified that the best result is attained for the values of α from 0.5 to 2.
In conclusion of this section, we describe the computations that are performed at each step in time.We start with construction of the data structure for efficient search of "neighbors" based on positions of particles.After that we compute the pressure and the friction forces, i.e., the right-hand side of Eq. (2.7).We also compute other forces, i.e., the term f i in Eq. (2.10), and a 1 , a 2 , a 3 , a 4 in the RK4 scheme as well as the "mean velocity" (a 1 + 2a 2 + 2a 3 + a 4 )/6.Then we correct modelling artefacts, e.g., too large velocities appearing occasionally (these velocities are renormalized to the prescribed maximum magnitude retaining their direction), move the particles to updated positions and update their velocities.Next, we process collisions with the boundary.Finally, we perform an additional correction of modelling artefacts, e.g., coinciding centers of particles.The utilization of spatial partitioning allows dividing the computation routine into almost independent threads that are executed in parallel.
Output parameters of the flow simulation are fiber concentration and orientation in the domain of flow.Assuming that the fibers are cylinders of length L with the cross-section area S, we can conclude that the mean volume of a fiber per one particle is LS/N, where N is the number of particles in a single fiber.Taking into account the relation L = (N − 1)l, where l is the distance between the adjacent particles of the fiber, we get the following formula for the fibers volume ratio in the representative volume element (RVE) Ω: Here, m is the number of fiber particles (mass points) in RVE and n is the total number of particles in this volume.
The model of elastic fiber used here allows the particles of the same fiber to be placed at some angle to each other.Therefore, the contributions of fiber segments into the effective elastic properties are taken into account separately, i.e., sum (3.5) is taken over all segments of all fibers inside the RVE.If r a and r b are the position vectors of the segment ends, then the angles ϕ ab , ψ ab can be computed by the formulas where e i are the basis vectors such that C ijkl used for computation of C IJKL in Eq. (3.5) are the components of the elastic moduli tensors in that basis.The set of these angles determines output orientation of the fibers.

Computation of Effective Moduli of Elasticity
We consider composite materials produced from high temperature modified phenolic resins [Bulgakov et al., 2013] and carbon fibers with a typical fiber lengthto-diameter ratio (l/d) equal to 500.The input data consist of the fiber volume ratio γ, Young's moduli of the fiber E f and the matrix E M , the Poisson's ratio of the matrix v M , the ratio k = l/d.Here, we have taken into account that the fiber Poisson's ratio v f does not affect effective moduli of the composite significantly.
A review of methods for computing effective moduli can be found in Thorvaldsen [2011].The history of this topic starts with the paper [Cox, 1952].The method we use assumes that the process of averaging can be divided into two stages (see, e.g., Christensen [1979]).First, we consider a material with parallel fibers.In this case, formulas for an approximate evaluation of five elasticity constants of transversalisotropic material can be derived using the Eshelby tensor [Eshelby, 1957;Christensen, 1979;Hashin and Rosen, 1964;Russel, 1973].
The second stage consists in spatial averaging using the effective moduli for parallel fibers that have been already obtained as the result of the first stage.For this purpose, the systems of parallel fibers with all possible directions are considered and the integrals are computed.Here, C IJKL (ϕ, θ) are the components of the tensor of elasticity moduli for the transversal-isotropic averaged medium associated with parallel fibers directed along the axis of transversal isotropy.The primed coordinate system is determined by the axis of transversal isotropy and is given by two Euler's angles ϕ, θ relative to the basic (not primed) system of coordinates.Equation (3.1) corresponds to the case of uniform spatial reinforcing when Euler's angles determining the directions of fibers vary within the ranges 0 ≤ ϕ ≤ 2π, 0 ≤ θ ≤ π.The components C IJKL are transformed according to the tensor transformation rule where Equation (3.1) can be used to derive the closed-form formulas for effective moduli of composites with uniformly oriented fibers (see, e.g., Christensen [1979]).However, uniform distribution of fibers in a material formed by means of molding is usually violated.In this case, in order to evaluate effective moduli of elasticity, we need to know the probability density ψ(ϕ, θ) of the fiber direction distribution over the angles ϕ, θ: the effective moduli can be obtained by computing the integral This integral is the mathematical expectation of the effective properties C IJKL with respect to all possible orientations of fibers.In many cases, including ours, this integral can be evaluated only numerically.For example, in case of a discrete set of fibers it is more efficient to compute the discrete mathematical expectation In this formula, n is the total number of fibers whose moduli are to be averaged together with resin moduli and ϕ i , θ i are the orientation angles of the ith fiber.

Numerical Examples
The first example concerning Poiseuille flow was considered for verification of the developed flow modelling method.The classic Poiseuille solution is known for an infinite cylindrical pipe of radius R with no-slip conditions on its boundary [Sutera and Skalak, 1993].The pressure does not depend on the radial coordinate and varies along the pipe with a constant gradient (falls by ∆p at distance of length ∆x), and the velocity at any cross-section is directed strictly along the pipe.Its magnitude is given by formula where r is the distance to the center of the pipe.For pipes of finite length with a steady laminar flow this formula is valid sufficiently far from the ends of the pipe.
We simulated the flow in a long pipe of radius R, and the ratio ε/R of the discretization parameter to the radius of the pipe was taken equal to 0.14.The simulation was performed for different values of the pressure and viscosity.The ratio υ/υ max was computed in a transverse section of the laminar flow domain.A typical simulation result is illustrated in Fig. 1 which represents the dependence of this ratio on the distance to the center of the pipe along with the theoretical values of the ratio.
One can see that the typical parabolic dependence was observed at least within the range of 0.7R from the center of the pipe.The deviation from a parabolic curve near the boundaries of the pipe is explained by the nonlaminar flow in simulation.The theoretical value of the maximal velocity was within 10% bounds compared to the computed value.So the first example concerns simulation of a pure fluid without fibers.The second example concerns the flow of a viscous resin with short fibers.The aim of this  example is to demonstrate the whole methodology in the case of a problem not accompanied by difficulties related to a complex geometry of the mold, boundary conditions, etc.The problem description is the following: the resin with short inclusions is forced by the pressure p 0 to flow into a rectangular parallelepiped ("mold") of size A × B × C through a cylindrical pipe ("channel") of diameter d.The parallelepiped is centered at the origin.The flow domain is schematically shown in Fig. 2. No-slip boundary conditions are posed on the walls of the channel and mold.The initial distribution of elastic fibers over the volume and fiber direction is uniform.The simulation results in effective moduli of elasticity computed in the domain of the mold.
We carried out simulation I with a small number of particles (18,000 particles in total, 3,000 particles corresponding to fiber; volume ratio of fiber ≈ 5%) and simulation II with a high number of particles (270,000 particles in total, 120,000 particles corresponding to fiber; volume ratio of fiber exceeds 20%).In simulation I, as well as in the problem of Poiseuille flow simulation, the ratio ε/d of the size of particles to the diameter of the channel was taken equal to 0.07.In simulation II, this ratio was taken equal to 0.03.The form had the sizes 7 × 1.5 × 4 relative to the radius of the channel.The distance from the walls of the mold to the center of the channel D was taken equal to the diameter d of the channel.
The flow process in simulation I is presented at different moments in Fig. 3.The resulting distribution of the fibers in the obtained sample was heterogeneous and anisotropic.Fibers concentration near borders was higher than at the center.Fibers orientation was most isotropic under the channel and most anisotropic along the borders of the form.
Computation time in simulation I was about 10 min on the AMD A10 APU graphic processor while simulation II took about 4 h.
If N denotes total number of particles, then in case of a fixed time step theoretical computational cost has asymptotics O(N log N ).One can see that this asymptotics holds for the considered example.If the simulation resolution is increased while the size of the computation domain is retained, the asymptotics becomes worse and takes the form O(N 4/3 log N ) because the time step has to be decreased proportionally to the decrease of the particle size.
We computed the concentration and orientation of fibers and the effective properties of the resulting composite at 10 control points indicated by thick marks in Fig. 2.These points are positioned at the height equal to one third of the parallelepiped height (half of the channel radius) equidistant from the front and the back boundaries.Distance between the neighboring control points in horizontal direction equals half of the channel radius.The control points are enumerated "left-to-right".The third, fourth and fifth control points are located directly under the channel.
Table 1 presents fiber volume ratios at control points in a representative neighborhood of a point computed in simulation II.This table also shows the degree of anisotropy (DoA) i.e., the operator norm Id/2 − 3A/2 , where Id is the identity operator and A is known as the second-order orientation tensor [Advani and Tucker III, 1987] given by the formula Here, n i is the unit vector directed along the ith fiber and π i (x) = n i • x.If the DoA value is close to one, then all the fibers are arranged along the same direction.Distributions close to uniform give the DoA close to zero.The effective Young's modulus E 11 and the shear modulus µ 12 were computed using Eq.(3.5) for the following values of elastic constants of the matrix and inclusions: E f = 230 GPa, v f = 0.3, E M = 5.4 GPa, v M = 0.3.The resulting values of E 11 and µ 12 are presented in the same Table 1.
We additionally carried simulation III that was similar to simulation II.The only difference was in the location of the channel: in simulation III the channel was centered relative to the form, i.e., the center of the channel was directly over the center of the mold.Table 2 presents fiber volume ratios and the DoA for simulation III along with the effective moduli E 11 and µ 12 .The tables show that the resulting distribution essentially depends on the position of the channel.This indicates that the method considered here can be used for optimization of molding parameters.Finally, we considered the real molding process.This molding process was performed on the injection molding machine DE3732R (max.Young's modulus computation adequate to the variance of experimental values, i.e., about 10%. In the second series of tests we also measured the Poisson's ratio.Its values are presented in the first column of Table 6.The second column contains the values of the Poisson's ratio computed by Eq. (3.4).The last column presents the values of the differences relative to experimental values.The maximal value of these differences is again about 10%.
Thus, the test of the accuracy provided by close-form formulas for uniformly oriented fibers indirectly indicates the feasibility of applying Eq. (3.5) to evaluate the effective moduli for a resin with short fibers for the volume ratio up to 30%.
Interestingly, an application of another well-known approach for determination of effective moduli, namely, computation of FEM solution for boundary-value problem formulated in the RVE with specific boundary conditions, resulted in values that were further from the experimentally measured ones: simulation of fibers using one-dimensional finite elements led to a deviation of up to 23%, and simulation of fibers as cylinders with slenderness ratio 100 led to a deviation of 15-20% for fibers concentration of 15-30%.

Conclusion
The simulation technique presented in the paper provides a complete methodology which can be used to predict the fiber distribution in the flow domain and subsequently to compute the elastic moduli of the polymerized material.The described methodology consists of two stages: the flow modelling by a weakly compressible SPH-like particle method and the computation of elastic properties using the fiber orientation and concentration distributions obtained at the first stage.
The primary feature of the developed particle method which distinguishes it from other similar WCSPH methods is the utilization of a physically justified partition of unity smoothing procedure instead of a single-point quadrature.This refinement keeps a computational simplicity of the SPH-family methods which provides an opportunity to divide the computation routine into almost independent parallel threads of execution using the modern parallel particle-interaction-system algorithms (e.g., the spatial hashing by means of parallel count sort which we used).
The presented technique is capable of identifying the domains where fibers get stuck or break, as well as domains where matrix failures of different kinds are observed.The knowledge of detailed fiber distribution additionally gives a way to estimate strength of the material produced by molding.
Along with the simulation method, we describe key points of its software implementation which delegates the task of computing the flow process to a GPU, thus achieving a significant performance improvement compared to the usage of a single CPU due to a parallelization of the computations.
An important fact regarding parallel methods is that nowadays there is a tendency in the hardware engineering to simplify processing units obtaining overall computational performance gain not by increasing clock frequency of a single processor core but, instead, by exploiting the massive parallel networks of very simple and cheap processors interconnected by a high-bandwidth memory bus.This fact makes development and implementation of parallel methods especially important.
The developed method can be utilized for the optimization of the flow input parameters (such as pressure, fiber concentration, etc.) thus increasing the quality of products produced by molding and subsequent polymerization.
The presented comparisons of the numerical results obtained by applying the developed simulation method with theoretical solutions and experimental data prove the validity of the method.

Fig. 1 .
Fig. 1.The velocity ratio υ/υmax for a transverse section of a flow in a cylindrical pipe.

Fig. 2 .
Fig. 2. Schematic representation of the simulation domain and control points.

Fig. 3 .
Fig. 3. Molding process for a resin with inclusions in simulation I.
injection volume −250 cm 3 , max.injection pressure 600 bar) at the Department of Chemical Technologies and New Materials of the Faculty of Chemistry, Lomonosov Moscow State University.Injection mold for standard testing samples with size of 10 × 4 × 80 mm according to ISO 178:2010 was used.The photos of injection mold and resulted samples are shown in Fig. 4. Samples of STN 150 resin with carbon fibers with total volume concentration of fibers 20% were prepared by double screw extruder, pelletized and injected in abovementioned mold.Curing of the resin was performed at 180 • C for 2 h inside the mold.

Table 1 .
Fibers volume ratios, the DoA and the effective moduli E 11 and µ 12 for simulation II.

Table 2 .
Fibers volume ratios, the DoA and the effective moduli E 11 and µ 12 for simulation III.