Efﬁcient analysis of dense ﬁber reinforcement using a reduced embedded formulation

In this paper we alleviate limitations of the conventional embedded reinforcement formulation for applications with dense ﬁber contents. We demonstrate that by condensing the ﬁber degrees of freedom during the assembly stage, the condition number of the resultant stiffness matrix is effectively reduced and the use of iterative solvers is facilitated even without preconditioning. Numerical benchmarks consisting of very large numbers of high aspect ratio discrete ﬁbers are performed, and major advantages are reported in terms of the computational efﬁciency of the solution method. We apply the solution method to a set of examples in the modeling of the ﬁber reinforced composites, speciﬁcally the estimation of the homogenized mechanical properties of the discontinuous ﬁber composites and modeling of the compression molding process. In the latter case, a specimen containing more than 2 million discrete ﬁbers is analyzed.


Introduction
The use of computational homogenization and particularly finite element based approaches is common in the study of fiber-reinforced composite materials [4,20,22,30,35]. These methods are desired due to their accurate estimation of the load bearing capacity and precise geometrical representation of non-homogeneously dispersed reinforcements. In addition they allow incorporating the interface imperfection or resistance phenomenon [12] which is important in the study of small size inclusions with high surface to volume ratio [10,31].
Nevertheless, the computational study of composites with high aspect ratio fibers is crucial in a wide range of applications, from carbon nanofibers (CNFs) and carbon nanotubes (CNTs) [6,16] to biological tissues of cytoskeleton filament networks [7,25]. The slender geometry of reinforcements however causes major challenges in generating a finite element mesh of sufficient quality. Even though RVEs with dense distributions of moderate aspect ratio fibers can be generated [23,26], extremely fine finite element descritiza-B M. Goudarzi goudarzi.mohsen@gmail.com 1 University of Twente, Enschede, The Netherlands tions are unavoidable when dealing with high aspect ratio fillers [20]. Adopting simplifying physical and geometrical assumptions, embedded reinforcement techniques [2,3,9,17] alleviate computational costs by eliminating the conformal mesh generation task. Generally, these methods assume a geometrical reduction of solid finite element fibers to equivalent line objects undergoing axial deformations only. Flexural reinforced composites can be modeled in a conceptually similar manner by introducing the bending stiffness of the fibers [24,29]. Ninić et al. [24] proposed an embedded model in which piles are discretized using beam elements, where a frictional contact formulation is used in order to simulate pile-soil interactions. In mentioned models, simplifications regarding the state of continuity of displacement fields are generally adopted [13], which allow the mesh conformity restrictions to be lifted.
To ensure load transfer between fibers and the matrix, a set of constraint equations is enforced in an embedded approach, usually through a penalty or interface type formulation [14] or frictional contact conditions [24]. The coupling scheme handles imperfection by allowing relative displacements (slip) between the fiber and matrix through an arbitrary traction-separation law. As a result, the size of the final system of equations is equal to the sum of the fiber and matrix degrees of freedom (DOFs). Therefore, computing the solu-tion for the resulting system of equations which becomes very large for dense fiber volume fractions is a major computational step when using the embedded reinforcement methods. An equally important issue is the badly scaled system due to the augmented structure of the stiffness matrix. This usually requires the use of expensive direct numerical solvers [14]. Iterative solvers can also be used, but mainly in combination with expensive preconditioning strategies.
Condensation [15,34] is generally used to reduce the size of the system of equations by eliminating a subset of the original DOFs, particularly in the context of the hp-FEM [5,21] where condensation of the internal (bubble) DOFs facilitates the solution of equations. As reported in [33], static condensation is equivalent to partial orthogonalization of basis functions or ILU preconditioning in the context of hp-FEM, where all can result to the same matrices. It is further confirmed that static condensation has positive effects on the conditioning of the system equations, and calculation of the Schur-complement of the full system acts as a natural preconditioner [33].
In this paper, the effectiveness of the condensation of the fiber DOFs is explored in an embedded reinforcement formulation. A reduced solution method is presented on the basis of the embedded formulation previously proposed in [14], which allows imperfections at the fiber-matrix interfaces and is generalized for applications with viscous (fluid) matrix material behavior in Sect. 2. A feature of the fiber model is the total independence of the fiber discretization from that of the matrix, which makes it a suitable candidate for problems with arbitrary shaped fibers (e.g. curved fibers) or situations involving continuously changing geometrical configurations. In addition, condensation of the fiber DOFs is performed during the assembly stage, keeping the size of the system of equations unchanged and equal to that of the intact matrix irrespective of the number of added fibers. The solution algorithm is described with a detailed study on the conditioning of the resultant stiffness matrices including the effect of preconditioning in Sect. 3. Numerical examples are given in Sect. 4 for applications which involve very large numbers of discrete fibers reaching up to around 2.5 million fibers. The effectiveness of iterative numerical solvers along with various preconditioning strategies is also examined throughout the paper.

Weak form of the system of equations
The general weak form of the two-phase reinforced system is where ∇ s is the symmetric gradient operator, and the total volume = m ∪ f is subdivided into matrix and fiber parts. The virtual displacement vector and stress tensor are shown by δu m and σ m for the matrix material, and similarly δu f and σ f are defined over the fiber region. The last term in (1) is the virtual mechanical work done by the interface tractions t c as a result of the virtual displacement slip vector δw = δ(u f − u m ) across the two sides of the interface.
For problems involving incompressibility as in the case of a viscous fluid matrix material, a stable u-p formulation [36] is used, where matrix stresses are subdivided into deviatoric (σ dev m ) and hydrostatic (σ p m = −pI) components and p represents the hydrostatic pressure. The weak form (1) is rewritten as: which is accompanied by the incompressibility condition: where ε v = − p K is the volumetric strain and K is the bulk modulus.

Discretized form
Adopting the standard finite element technique, the composite is discretized into matrix and fiber DOFs. The discretized forms are provided for the matrix, fiber and interface separately.

Matrix
The matrix displacement field u m is approximated as where N u contains elemental shape functions, and d m is the nodal displacement vector defined over the matrix material. For a linear elastic matrix material, the stresses are defined as with D m the matrix of elastic constants expressed according to the Hooke's law, and B u representing the conventional strain-displacement relationship. Following the standard procedures applied to the first term in (1), the stiffness matrix is obtained as For a viscous matrix material, the augmented pressure field is discretized in a similar way to displacement components as where N p contains the element shape functions vector and p is the nodal pressure vector. The stress components are written as, where η is the matrix viscosity andḋ m represents the nodal velocity vector. For the considered three-dimensional setting, where m = 1 1 1 0 0 0 T , Following the standard procedures applied to the first and second terms of the weak forms (2) and (3), we have: where the stiffness terms are the conventional terms arising in a u-p formulation [36].

Fiber
Fibers are assumed to be linearly elastic, and are discretized using conventional truss elements. The additional stiffness term due to the fiber segment k as shown in Fig. 1 is: where E m and E k f are the matrix and fiber Young's moduli, respectively, A k is the cross sectional area (A k = π d 2 f /4 with d f the fiber diameter), and L k is the fiber element length. Similar to Ref. [14], the effective stiffness (E k f − E m ) used in (12) cancels the extra matrix contributions at regions where high aspect ratio fibers are present. In 3-D, the fiber orientation matrix is with where x a and x b indicate the coordinates of the fiber element nodes for the fiber segment k as shown in Fig. 1.

Fiber-matrix interface
Zero-thickness interface elements are included between matrix and fiber, where the line interface elements are assumed to be conformal to the fiber discretization. The local interface slip vector w is discretized as: where R is the rotation matrix from the global to the local coordinate system. For the fiber segment highlighted in Fig. 1, the interface shape functions matrix and corresponding displacement vector are written as and respectively, where I is the identity matrix of size three, N f a and N f b are the one-dimensional Lagrange shape functions at nodes a and b defined local to the fiber, and the shape function N m i is evaluated at any point along the fiber segment, and is defined for node i of the parent matrix element, with n being the number of nodes in an element. It should be mentioned that in the present formulation and as shown in Fig. 1, the fiber DOFs can be arbitrarily located with respect to the matrix, and not essentially at the intersections of the fiber and parent elements. Non-straight fibers or situations where matrix/fiber geometry changes during the simulations are therefore handled more efficiently, without complications associated with finding the intersection points as in the traditional approaches [8]. In order to evaluate interface contributions, a nodal integration technique is utilized, where the corresponding parent element and the local coordinates need to be identified for the evaluation of the interface shape functions matrix (16) at each fiber (integration) point. Therefore the formulation is general and each end point of a fiber segment can be located in an arbitrary matrix element. Due to relative sliding between the fiber and matrix and for situations where matrix/fiber geometry changes during the simulations, this procedure is repeated at each step of the simulation, and is further enhanced using efficient algorithms, e.g. kd-tree data structures which are widely applied for the fast determination of neighboring points [11].
Knowing the interface slip vector, the interface traction vector for a mechanical system is defined as: where in practice any interfacial rule can be used. Here, in order to keep the presentation simple, we assume a linear interface law where is defined as the linear interface constitutive matrix. The constant K bn represents the stiffness of the interface in directions normal to the fiber lines, and K bt represents the stiffness of the interface tangential to the fiber axis.
For a viscous matrix material, the interface tractions are defined in a similar way but are instead written as a function of the slip rate vector: where the interface constitutive matrix D be = ∂ t c ∂ẇ is defined similar to the mechanical system.
Finally and upon substituting the interface slip and traction vectors into the last term of (1) or (2), the interface contribution related to the kth interface/fiber element becomes with C k f = π d f the circumference of the fiber, L k the length of fiber-matrix interface element, and s the local coordinate along the interface element. The interface stiffness contribution (20) can be expanded as: being subdivided into matrix and fiber counterparts and ensuring coupling between the matrix and fiber displacements.

System of equations and reduced form
Inserting the discretized fields into the weak form (1), the elemental system of equations for the elastic matrix is: where d m and d f are the vectors of the matrix and fiber DOFs at the element level, respectively, and the stiffness contributions are defined in the previous section. For the Newtonian viscous matrix, the equations are written in a rate form and assuming a Lagrangian framework, the geometry is updated at each time increment. Thus the following set of equations is solved at each step: where t is the time increment, and the field variables are updated in each time increment, with d m , d f and p being the vectors of matrix, fiber and pressure increments.
Global system Irrespective of the adopted material model for the matrix and assuming that the composite system is composed of n f fibers, the global stiffness matrix and dis-placement/force vectors can be written as: where U M is the global matrix displacement vector (comprising both displacement and pressure DOFs in the case of a u-p formulation) and U F i (i = 1, . . . , n f ) contains the DOFs attributed to each fiber separately. The off-diagonal terms are due to the fiber-matrix interface and ensure coupling between them in view of (21). As the fibers do not interact, no coupling terms arises between them.
Reduced system In order to efficiently solve the global stiffness matrix in (24), a condensation technique is applied. A subset of DOFs (here chosen as the matrix DOFs in vector U M ) is designated as the master DOFs for which the system of equations need to be solved. The DOFs attributed to each fiber are then treated as the internal or slave DOFs, and are eliminated from the system of equations during the assembly procedure without constructing a full size stiffness matrix. This selection ensures that the size of the system which needs to be passed to a numerical solver is noticeably reduced to that of the intact matrix. Condensation is performed by expressing the internal DOFs as a function of master DOFs: where by substituting (25) into the first row of (24), the following reduced system results: where the reduced stiffness matrix is expressed as which is also known as the Schur complement of the augmented stiffness matrix block (24). Simultaneously in absence of fiber forces in (24), the reduced force vector is As the fiber contributions are independent, the set of DOFs attributed to each fiber is treated as an individual set of slave DOFs. Therefore the computation of inverse and arithmetic operations in (26) remains quite cheap, as the size of the fiber stiffness matrices is usually relatively small, and the evaluation of their inverses is straightforward. It should however be mentioned that . . , n f ) must be rank sufficient, which is the case in view of the appropriate coupling terms arising between fibers and matrix through the interface contributions (21).
The matrix DOFs can be directly solved by: and subsequently the fiber DOFs are: In practice and due to the small size of the internal systems, a direct numerical solver can be used for solving the fiber DOFs. In addition and as the systems are independently solved, an efficient parallel implementation can be applied. This should however be noticed that for some practical cases, e.g. evaluation of the homogenized mechanical properties of an RVE discussed in Sect. 4.1, fiber displacements (U F i ) are not explicitly required as all the homogenized properties can be extracted by having the values of the forces along the edges of the matrix.

Reduced versus full system
In the first series of numerical examples, the fiber-matrix system is studied using the mechanical formulation as described in Sect. 2. Various characteristics of the full and reduced systems are compared.

Condition number assessment
A well conditioned system of equations is generally favored as the performance of numerical iterative solvers is directly proportional to the value of condition number. Therefore a test on the condition number of the reduced and full systems is performed in this section. A unit size elastic cube (Fig. 2) is considered with all the displacement components constrained at the left face (x = 0), and with homogeneous horizontal displacements of 0.05L imposed along the right face (x = L). Unless specified differently, the matrix material is discretized into a grid of 11 × 11 × 11 H8 elements. Matrix Young's modulus and Poisson's ratio are assumed to be E m = 1 MPa and ν m = 0.2, respectively. Analysis is performed considering various values of the fiber Young's modulus, with E f /E m = 10, 100, 1000. The interface constant K bt (= K bn ) varies between 1 and 10 5 N/mm 3 to cover a wide range of weak and stiff interface conditions. Single fiber The simplest case of one fiber with length 0.6L and diameter 0.05L is considered first. The fiber is discretized into 20 equally sized line elements, and is positioned along the center axis of the cube aligned with the global x axis (Fig. 2). To facilitate the numerical calculation of the condition number, the matrix material is discretized into a coarse grid of 5 × 5 × 5 H8 elements. The conventional definition of the condition number as the ratio of the largest to the smallest eigenvalues is evaluated by using the Matlab cond function, which determines the expansion of eigenvalues of the stiffness matrix. The values are reported in Fig. 3 for both the augmented system and the reduced version. Various ratios of the interface constant (K bt ) and fiber Young's modulus (E f ) are adopted. It is observed that a single fiber can considerably increase the condition number of augmented or full system. This however is not the case when the reduced model is used, for which the attributed behavior is very close to that of the intact matrix (with no fiber). This is more visible when the fiber and interface constants are less stiff. In general, by increasing the fiber stiffness and the interface constant, an increased condition number is observed for both models, while in any case the condition number of the reduced model is considerably smaller than that of the full system. Multiple fibers We now study the conditioning behavior of the composite by increasing the number of fibers. Distributions with up to around 10,000 randomly oriented fibers of the same length l f = 0.2L are considered, where each fiber is subdivided into 5 equally sized segments. It is assumed that the fibers are uniformly distributed. As shown in Fig. 4, a behavior similar to that of the single fiber case is observed for the multiple fibers problem, where the estimated condition number of the reduced system is considerably smaller than that of the full system, and is in a range comparable to what is estimated for the intact matrix. Worth mentioning is that for this setup and for both systems, the estimated condition number does not considerably alter by increasing the number of fibers.

Sparsity pattern of the reduced system
The sparsity pattern of the stiffness matrix is another important factor that influences the performance of iterative solvers. The number of non-zero entries and corresponding sparsity patterns of the reduced system for various numbers of fibers are shown in Figs. 5 and 6, respectively. An interesting feature of the reduced model is that the number of non-zero entries or equivalently the sparsity pattern converges to a cer- tain state. This is due to the fact that when a certain number of fibers is reached, all the elements of the bulk mesh become involved in the condensation of fiber DOFs. Such behavior of the reduced system is advantageous on the properties of the applied numerical solvers, with the performance of a sparse solver being usually influenced by the number of non-zero entries of the stiffness matrix. Finally it should be mentioned that the reported values were obtained for a fixed length of fibers, and a different sparsity or number of non-zero entries is expected when the length of the fibers is changed.

Iterative solver
As shown in the previous section, the reduced system of equations is better conditioned compared than the full system. This is an important feature of the reduced model when it comes to the use of iterative numerical solvers, with their convergence rate being directly proportional to the condition number. The generalized minimal residual method (GMRES) is adopted here to solve the resultant stiffness matrices for the multiple fiber example, and the results are reported in Fig. 7 for both the reduced and full systems. As shown, the rate of convergence is considerably slower for the full system than for the reduced one. This is clearly noticeable even for the case where the small number of 500 fibers is added to the system. In line with previous observations, performance of the iterative numerical solvers is therefore considerably better when applied to the reduced system and remains almost unchanged irrespective of the number of fibers. Upon increasing the number of fibers, no considerable performance changes of the iterative GMRES solver are observed. The effect of the interface stiffness parameter (K bt ) is shown in Fig. 8, indicating that by increasing the parameter, an slower convergence rate for the iterative solver is obtained. Again this behavior is attributed to an increased condition number of the reduced model in these cases as previously shown in Fig. 4a.

Preconditioning of the reduced system
As shown in the previous section, iterative solvers can be applied for the reduced model without additional treatments. This however only holds true for cases where the intact matrix system is well conditioned, and generally is not the case when a u-p formulation is adopted. For these cases the iterative numerical solvers can still be used provided that suitable preconditioning of the stiffness matrix is applied. A preconditioning matrix M is calculated such that approximates the stiffness matrix K , and its inverse is multiplied to the system of equations (28) from right or left during an iterative algorithm. Among various preconditioners in the literature, the incomplete LU (ILU) method [28] is adopted, which approximates the LU factorization of the stiffness matrix. The resultant matrix has sparsity pattern similar to the original stiffness matrix, and commonly serves as a preconditioner Relative errors versus the number of iterations for the full and reduced systems without preconditioning using the GMRES iterative solver. Solver performance is compared for different numbers of fibers, where E f /E m = 100 and K bt = 100 N/mm 3 . Solid lines: reduced system, dashed lines: full system Fig. 8 Relative errors versus the number of iterations for the reduced system without preconditioning using the GMRES iterative solver. Solver performance is compared for different values of K bt , where E f /E m = 100 and a total of 10,000 fibers are added for the GMRES solver. The simplest choice for M can be constructed such that it is much more sparse than the original matrix by fill-reducing the matrix's unknowns (ILU0, i.e. allowing no level of fill-in [28]). More accurate LU factorizations can be achieved if fill-in is allowed. The level of fill-in is controlled by a non-negative scalar variable defined as drop tolerance, where following the definitions used in Matlab ilu function, all entries which are smaller in magnitude than the drop tolerance are dropped in constructing the preconditioner matrix M. Three levels of fill-in are considered here, where ILU1, ILU2 and ILU3 represent the cases where drop tolerance is equal to 1e−3, 1e−6 and 1e−11, respectively.
The stiffness matrices are premultiplied by the preconditioning matrix, and the estimated condition numbers are compared for the full and reduced systems in Fig. 9. Different preconditioning levels as discussed above are used. The condition number is effectively reduced for both systems, obviously especially when a smaller threshold is used. For the case of ILU3, a highly accurate factorization is performed and the condition number of the resultant matrix is very close to unity. This should however be noticed that for the full system and in case many fibers are involved, the complexity of the procedure for calculating the preconditioning matrix is significant as the size of stiffness matrix increases substantially. The case of multiple fibers will be further studied in the upcoming sections by using the reduced model.

Computational modeling of dense fiber composites
In the next set of numerical examples, the reduced model is applied for situations with dense distributions of discrete fibers. Handling such dense fiber distributions is computationally prohibitive if the full model is used. The numerical simulations start with evaluating the effective mechanical properties of a reinforced box, where the behavior of iterative solvers is also tested. The second numerical example deals with a viscous material in which geometrical features change during the simulations. This example can be realized as a step towards modeling the compression molding process of reinforced polymers [27], which entails incorporating a very large number of discrete fibers.

Homogenized mechanical properties
The numerical algorithm is applied in evaluating the homogenized mechanical properties of the composites with dense fiber distributions. The elastic cube as described in Sect. 3.1 is considered. For extracting the homogenized mechanical properties, the numerical procedure as explained in reference [14, Appendix A] is used, where periodic boundary conditions attributed to various components of the strain tensor are prescribed along the edges of the simulation box. Because the homogenized properties of the elastic composite are of interest, the average stresses can be directly obtained through evaluating the forces along the edges of the simulation box. Fiber distributions are generated using the random sequential adsorption technique [18], with fibers length and diameter set equal to 0.2L and 0.004L, respectively, and each fiber is subdivided into 10 equal sized segments. Fibers are oriented along the x axis and volume fractions up to ν f = 40% are considered (leading to around 138,000 fibers). It is worth mentioning that using a standard formulation in [14], the number of fibers was limited to around 22,000. As discussed earlier in Sect. 2.2, the current reduced formulation can be readily applied for arbitrary fiber geometries as well, and therefore two more fiber curvatures as shown in Fig. 10 are also examined. The homogenized mechanical properties are reported for two different values of the interface constant and various fiber volume fractions in Fig. 11. The homogenized Young's moduli E c x are reported for different mesh densities to confirm the convergence of the numerical results. A good agreement is observed with the reference Voigt estimate with the perfect interface condition (straight fibers) in Fig. 11a, which further confirms that the simulations using the material properties adopted in the present study are not affected by the intraply shear locking reported in [32]. For the smaller interface constant (K bt = 500 N/mm 3 ), the homogenized Young's modulus is considerably lower than the rule of mixture prediction (Voigt estimate). A similar reduction in Young's modulus is observed for the curved fiber distribution. A more detailed study of the fiber curvature effect is shown in Fig. 12 in terms of various mechanical properties, where an increased transverse Young's modulus E c y and shear modulus G c xy are observed. The Young's modulus is however lower for straight fibers along the transverse direction and agrees well with the lower bound estimated by the inverse rule of mixture (Reuss estimate).

Solver performance
The performance of the model can be improved by applying a preconditioning strategy to the iterative solver. Calculation of the preconditioning matrix however entails additional computational costs. Here instead of computing the preconditioning matrix for each distribution of fibers, one certain configuration is selected as the reference state, for which the preconditioning matrix is computed and is subsequently used for different fiber distributions. For the ideal performance of the iterative solver, the sparsity pattern of the chosen configuration should well approximate the stiffness matrix for each problem. This procedure is simplified for the reduced system as the sparsity pattern of the stiffness matrix resembles that of the intact matrix material, and as shown in Fig. 6, converges to a certain state upon increasing the number of fibers.
The GMRES solver results are shown in Fig. 13 for different fiber volume fractions. As expected, the worst convergence behavior is for the case where no preconditioning is used. An improved behavior is however achieved by applying preconditioning. In Fig. 13a, the preconditioning matrix is calculated for an intact matrix (without fiber) and used at all fiber volume fractions. Various cases with different level of fill-in as discussed in Sect. 3.4 are considered. Differently in Fig. 13b, the preconditioning matrix is calculated for the most dense fiber distribution of 40%. Both preconditioning strategies lead to significantly reduced number of iterations to gain a converged solution for all fiber contents compared to the case where no preconditioning is used. This improvement in performance is more pronounced when smaller thresholds for the fill-in are allowed. Interestingly, using the preconditioner as in Fig. 13b provides a more pronounced improvement at all fiber contents and even when fill-in is not allowed (ILU0). This is attributed to the fact that the structure of the non-zero entries of the stiffness matrix for the reinforced system is similar at various fiber distributions (Fig. 6), and as already mentioned, converges to a certain state upon increasing the number of fibers. As a result and as shown in Fig. 13b, the iterative solver performs best at the most dense fiber dis-  Fig. 14b by using the most dense fiber distribution with ν f = 40%, the iterative solver convergence is achieved pronouncedly faster than in Fig. 14a for the intact matrix.

Viscous matrix material in the compression molding process
The goal of this final example is to utilize the reduced model in situations with dense fiber contents, where the geometrical features continuously change during the simulations. As already pointed out in Sect. 2, the current numerical model is particularly suitable for such situations because finding intersection points between the fibers and matrix mesh is not needed. Instead, the parent matrix element of all fiber nodes should be identified at every time increment. An application is in the production of the thermoplastic composites, where the compression molding process of fiber reinforced polymer melt is simulated [27]. As shown in Fig. 15, a viscous melt comprising of a large number of high aspect ratio fibers is squeezed under external pressure loading and as a result of the melt flow, the fibers are redistributed. An accurate representation of the fiber reorientation is of importance for identifying the mechanical properties of the resultant composite. The reduced formulation is used, and the incompressible viscous matrix flow is described by the u-p formulation explained in Sect. 2. Here we focus mainly on the basic numerical model, and a more behavioral study is postponed for future contributions.
In the studied problem, the reinforcements are in shape of chopped thermoplastic prepregs or flakes (Fig. 16), where each flake is comprised of parallel high aspect ratio fibers. A cubic sample of size L with a random distribution of same size fibers as shown in Fig. 17 is considered. In generation of the samples, contact between flakes is avoided, and similarly within each flake, fibers do not overlap. Up to around 25,000 randomly distributed flakes are considered, while each flake is comprised of around 100 fibers. As such the number of discrete fibers can reach up to 2.5 million fibers. leading to a volume averaged orientation tensor A ≈ diag (1/3, 1/3, 1/3) revealing isotropy on the macroscopic scale, with the components evaluated as [1,19]: where n f is the number of fibers and p k i (i = 1, 2, 3) are components of the fiber unit vector.
The box is squeezed under a uniform compressive constant velocityu y = 0.01L 1/s from the top face (y = L). Velocity components are constrained at the bottom face (y = 0), and the lateral flows along the x and z directions are allowed by treating the lateral faces of the box as free edges. The problem is solved for a total of 55 time steps with a constant time increment t = 1 s. The state of the viscous matrix is described by a nearly incompressible material represented by a high value of the bulk modulus (K ). The viscosity (η) of polymer matrix material (PEEK) is chosen as 500 Pa s, and the carbon fiber's Young's modulus is E f = 200 GPa. The length and diameter of the fibers are set to 0.2L and 0.0007L, respectively, and the interface constant (K bt = K bn ) is chosen such that (K bt L)/η = 24. The matrix material is discretized using non-structured tetrahedral elements. To ensure numerical stability of the u-p formulation and in order to eliminate the locking in the formulation, the approximation order of the displacement field should be higher than that of the pressure [36]. As such, quadratic tetrahedral elements (T10) are used for the displacement field, and the pressure is discretized using conventional linear tetrahedral (T4) elements. The fibers are discretized using linear one-dimensional Lagrange elements, and each fiber is discretized into 5 segments.
The force displacement curves are shown in Fig. 18a for the most dense distribution of flakes. Results are reported for different mesh densities, and show convergence in the global response of the composite melt. Curves are also reported for various numbers of fibers in Fig. 18b for the finest mesh, where as expected, a stiffer response is observed as a result of adding the fibers. The diagonal components of the orientation tensor A are plotted in Fig. 19 for different fiber distributions. As the model allows fiber reorientation, the initial spatially isotropic distribution of fibers changes to a planar isotropic distribution, with the x and z components of the orientation tensor approaching 0.5, and the y component becoming zero as a result of the applied compression along the y direction. Deformed configurations are also shown in Fig. 20 at different stages of the simulation. It should be mentioned that a re-meshing algorithm was not applied here yet, but potentially becomes inevitable as the deformations are further increased.
An interesting feature of the reduced model which was also discussed in the previous example is that iterative solvers can be used with efficient preconditioning. In this example by considering the u-p formulation with incompressibility, the condition number of the intact matrix is relatively high, and the convergence rate of the iterative solver without preconditioning is very slow. Therefore preconditioning should  be applied similar to what already done for the previous example. Again the intact matrix and the most densely reinforced configuration (2.5 million fibers), both at the initial configuration, are chosen as the reference state for computing the preconditioning matrix, where fill-in threshold is set to 1e−6 (ILU2). As shown in Fig. 21 and in line with previous observations in Sect. 4.1.1, the convergence rate of the iterative solver is significantly improved by applying preconditioning for both reference states. Although the number of iterations is only slightly lower for the case that a preconditioner is computed for the most densely reinforced configuration (solid lines) compared to the intact reference state (dashed lines), the former requires more demanding calculations as the number of non-zero entries is significantly higher. Improvements are most notable at the initial steps, and as the simulations continue and larger displacements incur, the performance of the iterative solver gradually deteriorates. Achievements are however still quite significant for relatively large displacements at step 50. A strategy that can be applied to further improve the performance of iterative solver is to recalculate the preconditioning matrix for the updated configuration in every few steps, and whenever certain displacements increments are reached. Certainly after remeshing to restore proper element shapes after very large displacements, the preconditioning matrix needs to be recalculated.

Conclusions
This paper presents an efficient embedded formulation that is particularly suitable for dense distributions of arbitrary shaped discrete fibers. The proposed algorithm ensures a total independence of the mesh on the fiber distribution from that of the matrix and relies on condensation of the fiber DOFs during the assembly stage, by which major complications of the traditional approaches are bypassed. Through detailed numerical solutions it is shown that the solution algorithm allows efficient computation of the effective mechanical properties and fiber reorientation analysis during the compression molding process without constructing a full-size stiffness matrix. The adopted approach also eliminates the expensive task of solving a full-size stiffness matrix. In addition to condensing the fiber DOFs, the condition number of the stiffness matrix reduces to a range comparable to that of the intact matrix (without fibers). Therefore the approach naturally facilitates the use of iterative numerical solvers, while an efficient preconditioning approach based on the ILU method was shown to further enhance the performance of the GMRES solver. In the proposed preconditioning approach, one certain configuration of the system is selected as the reference state, for which the preconditioning matrix is computed and is subsequently used for different fiber distributions. Another feature of the current numerical framework is the total geometrical independence of the fiber and matrix elements, where calculation of the intersection points between the fiber and matrix elements is not required. This property of the method is particularly advantageous in situations where costs attributed to the calculation of the intersection points are considerable, e.g. the case with non-straight fibers or continuously changing geometrical properties. Both cases have been investigated in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.