Micromechanical Modeling of Anisotropy and Strain Rate Dependence of Short-Fiber-Reinforced Thermoplastics

: The anisotropy and strain rate dependence of the mechanical response of short-ﬁber-reinforced thermoplastics was studied using a straightforward micromechanical ﬁnite element analysis of representative volume elements (RVEs). RVEs are created based on the ﬁber orientation tensor, which quantiﬁes the processing-induced ﬁber orientation distribution. The matrix is described by a strain rate-dependent constitutive model (the Eindhoven glassy polymer (EGP) model), which accurately captures the intrinsic response of amorphous polymers. The micromechanical results indicate that the inﬂuence of strain rate and that of the loading direction on the yield stress are multiplicatively decouplable, which conﬁrms previous experimental observations. Moreover, it is demonstrated that the yield stress, to a good approximation, can be directly linked to the ﬁber orientation in the direction of loading. This leads to a new relation that uniquely links the rate dependence of the yield stress to the ﬁber orientation in loading direction.


Introduction
Due to their high specific strength, short-fiber-reinforced thermoplastics are currently used in a multitude of load-bearing applications where structural integrity and long-term reliability are key requirements, such as components used in automotive or aeronautic applications [1]. Hence, the ability to accurately predict the mechanical performance of such systems is of the utmost importance. Unfortunately, however, the accurate prediction of the performance of short-fiber-reinforced thermoplastic components still faces significant challenges.
One of these challenges is related to the fact that the yield strength of short-fiberreinforced thermoplastics is generally strongly anisotropic due to processing-induced fiber orientation. During injection molding, the fibers are oriented during the flow of the molten fluid into the mold [2,3]. Due to spatial variation of the flow conditions throughout the mold, the material will display a marked variation in its anisotropic response on different positions of the molded product [4]. This intrinsic coupling between flow conditions and resulting mechanical properties makes the prediction of the mechanical performance of the molded part challenging. An additional complexity is that experimental observations clearly give evidence for a significant strain rate dependence of the yield strength on fiber-reinforced polymer composites [5][6][7][8]. Since the reinforcing fibers generally remain elastic until fracture [9], this dependence on strain rate is intrinsically related to the strain rate dependence of the mechanical response of the thermoplastic matrix, which is directly related to the influence of applied stress on chain mobility [10,11].
To adequately predict the strength of products based on fiber-reinforced thermoplastics, both the processing-induced fiber orientation as well as the strain rate dependence of the matrix must be considered. The difficulty in finding the optimal shape of a load-bearing component lies in the fact that changes in shape (even small ones) directly influence local flow conditions that subsequently alter the resulting fiber orientation distribution and the related anisotropic, time-dependent response. The inability to accurately predict the local anisotropic, rate-dependent mechanical properties directly from the fiber orientation distribution is the main cause that reliability of such products can only be warranted after full-scale testing-a time-consuming and costly procedure which hampers optimization.
Results from a recent study in our group, however, have opened a new, promising possibility to directly link fiber orientation to the anisotropic time-dependent response [12,13]. In an extensive experimental effort, it was demonstrated that the rate-dependent, anisotropic (yield) strength of short and long fiber-reinforced thermoplastics generally displays a multiplicative decoupling of the effects of strain rate and loading angle [12,13]. Additionally, it was shown that, at a fixed strain rate, the yield strength appears to be uniquely coupled to the average fiber orientation in the direction of testing [14]. This observation, in combination with the decomposition of the influence of strain rate and that of fiber orientation, could lead to a novel, fast method to accurately predict the anisotropic rate-dependent mechanical response for an arbitrary fiber orientation distribution.
The objective of the current manuscript was to further explore the experimentally observed decoupling of the orientation and strain rate dependence by using a micromechanical approach. Micromechanical modeling, based on the spatial fiber orientation distribution and the microscopic matrix behavior, is a method to intrinsically resolve the yield mechanism of composite materials. Compared to mean-field composite stiffness models, such as the Mori-Tanaka model [15,16] and pseudo-grain decomposition followed by two-step homogenization [17][18][19][20] of composite structures, full field models, such as that based on the finite element approach, are able to evaluate micro-level stress and the effects of interaction between adjacent inclusions on the overall properties. Micromechanical finite element analyses of representative volume elements of short-fiber-reinforced composites [21][22][23] showed good agreement between experimental and simulation results under some conditions, but the strain rate-dependent behavior of the composite was not taken into consideration.
A proper rate-dependent constitutive model for polymers is essential to predict the yield strength of the short-fiber-reinforced thermoplastics. A few researchers used a strain rate-dependent constitutive model and extended its utilization for polymer-based composites [24,25]. However, the constitutive model was developed based on the viscoplastic deformation of metals, and lacks the deformation kinetics of the polymers. A landmark in glassy polymer modeling was the pioneering work from Haward and Thackray [26], who proposed a model that considers the contributions of inter-molecular interactions and the molecular network to the mechanical properties. For a glassy polymer, the inter-molecular interactions, as a visco-elastic contribution, determine the moderate strain behavior including yield and strain softening. The molecular network, as a rubber-elastic contribution, governs the strain hardening behavior at large strains. Based on this pioneering work, further major developments were made by Boyce et al. (BPA model) [27][28][29], Buckley et al. (OGR model) [30,31] and our group in Eindhoven (EGP model) [32,33].
In this article, a fine-scale micromechanical finite element model was developed based on the rate-dependent Eindhoven glassy polymer (EGP) model for the thermoplastic matrix of short-fiber-reinforced composites and the processing-induced orientation of elastic fibers. This micromechanical model was used to investigate the anisotropic yield strength at the different strain rates and fiber orientation distributions. The experimentally observed multiplicative decoupling was proven in this study. Furthermore, a unique structureproperty relation between the fiber orientation distribution and the anisotropic mechanical properties was found.

Micromechanical Modeling Approach
In this section, details of the micromechanical model for evaluating the anisotropic and rate-dependent behavior of short-fiber-reinforced thermoplastics are discussed. An E-glass short-fiber-reinforced polycarbonate material was considered. The anisotropic behavior was captured through finite element analyses (FEA) of representative volume elements (RVEs) at different loading angles. The rate-dependent behavior of the thermoplastic matrix was described with the Eindhoven glassy polymer (EGP) model. The fibers are modeled as isotropic elastic with a Young's modulus of 7.2 GPa and a Poisson's ratio of 0.22. Fibers and matrix are assumed to be perfectly bonded. A schematic work flow of the micromechanical model is shown in Figure 1.

Representative Volume Elements
To quantify the fiber orientation distribution, the fiber orientation tensor [34] was introduced. The second order fiber orientation tensor p is defined as the average of the dyadic product of orientation vectors for a number of m fibers [34]: with − → q i the orientation vector of an individual fiber and where the trace of the fiber orientation tensor p is equal to 1. For investigating the anisotropic behavior of a microstructure, the fiber orientation tensor is rotated by different angles relative to the global coordinate system. The rotated fiber orientation tensor p is given by where R is a rotation tensor. To conduct the finite element simulations of the E-glass short-fiber-reinforced polycarbonate, representative volume elements (RVEs) corresponding to different rotation angles which are generated using software package Digimat-FE 2018.1 [35], ensuring a geometric periodicity of the RVEs. This means that if a fiber passes one boundary, the remaining part of that fiber continues at the opposite side of the RVE. Consequently, if these RVEs are periodically stacked, a geometrically continuous structure is obtained. Due to computational limitations, an RVE with a low weight fraction and aspect ratio was considered in this paper. A total of 18 fibers with a weight fraction of 10% (volume fraction 5%), diameter of 0.02 mm and aspect ratio of 12 were distributed according to the fiber orientation tensor p . Based on this geometric information, the RVE size was determined at 0.3 × 0.3 × 0.3 mm 3 . Consistent results were obtained for the RVE size of 0.3 × 0.3 × 0.3 mm 3 and prove the RVE size to be sufficiently large enough. The fibers were assumed to be cylindrical without intersection or contact. In addition, the fibers were separated as far as possible to avoid clustering effects through tuning the nearest neighboring distance. The representative volume elements are finally meshed with 80 × 80 × 80 linear hexahedral solid elements.
In this article, a transversely isotropic E-glass short-fiber-reinforced polycarbonate was considered with the third direction as the symmetry axis. Therefore, the diagonal components of the 3 by 3 fiber orientation tensor must satisfy p 11 = p 22 = (1 − p 33 )/2, and the non-diagonal components are 0. Three transversely isotropic fiber orientation distributions, with p 33 = 0.94, p 33 = 0.82 and p 33 = 0.55 were, respectively, introduced to represent three microstructures. Each microstructure is used to study the corresponding anisotropic behavior. Because of the limited number of fibers in a representative volume element, the desired transversely isotropic fiber orientation tensor p is approximated with a tolerance of 0.05 in each of its components.
Rotated microstructures are created based on the rotated fiber orientation tensors p , as shown in Figure 2. In order to simplify the description of the fiber orientation, a fiber orientation factor p σσ is extracted from the rotated fiber orientation tensor p and it represents the fiber orientation along the loading direction. When p σσ = 1, all fibers are aligned with the loading direction and when p σσ = 0, all fibers are perpendicular to the loading direction. In the results section, the orientation of the fiber is indicated by p σσ . At each angle, the fiber orientation factor p σσ corresponding to loading in the z direction is illustrated.

Periodic Boundary Condition
Periodic boundary conditions were applied to the periodic RVEs to homogenize the response of the heterogeneous materials [36]. The periodic boundary conditions ensure that the deformed shapes of two opposing faces remain identical. For this purpose, the boundary of the representative volume element is split into six boundary faces, here denoted as Γ 4873 , Γ 1562 , Γ 1584 , Γ 2673 , Γ 5678 and Γ 1234 (see Figure 3). The periodicity can then be formulated in terms of the displacements u of the nodes in the boundary faces and corner nodes: The RVE is loaded with a uniaxial strain at a macroscopically constant strain rateε by prescribing the displacement of corner node 5 as where L 0 is the initial length of the RVE in the z direction. Additionally, the displacement of corner node 1 is suppressed in the x, y and z direction, the displacement of corner node 5 is suppressed in the x direction, the displacement of corner node 2 is suppressed in the z direction and the displacement of corner node 4 is suppressed in the y direction. The macroscopic first Piola-Kirchhoff stress tensor P M is obtained as where − → f e i is the reaction force and − → x 0 i is the position vector of the corner nodes i = 1, 2, 4, 5 in the reference state. V 0 is the volume of the RVE in the reference state. The corresponding macroscopic Cauchy stress σ M is then given by where F M is the macroscopic deformation gradient, which is obtained as with − → n 0 the outward normal of the boundary Γ 0 in the initial state.

The Strain Rate-Dependent EGP Model
The rate dependence of the E-glass short-fiber-reinforced polycarbonate stems from the rate-dependent polycarbonate matrix. The EGP model can accurately describe the strain rate-dependent visco-elastic regime, yield stress, strain softening and hardening of amorphous polymers and is used here to describe the polycarbonate matrix. The small strain behavior, including yield and strain softening, is determined by the inter-molecular interactions as a visco-elastic contribution. The larger strain behavior with the strain hardening is governed by the molecular network as a rubber-elastic contribution. Based on these two contributions, the total Cauchy stress σ is split into an elasto-viscoplastic driving stress σ s and an elastic hardening stress σ r : The elastic and plastic parts of the total deformation are defined through a multiplicative decomposition of the deformation gradient tensor F: where the subscripts e and p refer to the elastic and plastic parts of the deformation, respectively. For polycarbonate, two relaxation processes α and β exist. However, since the material is tested at room temperature and the strain rates are lower than 10 0 s −1 , the β process, which dominates at higher strain rates, is ignored [37]. Therefore, in order to save valuable computation time, the driving stress σ s only represents the α process: where subscript k denotes the individual modes; n is the total number of modes; κ is the constant bulk modulus; J = det(F) = det(F e ) is the volume change ratio; I is the secondorder unit tensor. G k and B d e k are, respectively, the shear modulus and the deviatoric part of the isochoric elastic left Cauchy-Green deformation tensor, which is defined as The relation between the driving stress and the plastic rate of deformation is described by using a non-Newtonian flow rule: where a scalar viscosity η depends on the absolute temperature T, equivalent stress τ, pressure p and the state parameter S: where η 0 k is the initial viscosity associated with mode k, H is the activation enthalpy, R is the universal gas constant, and µ is the pressure dependence. τ 0 is the characteristic shear stress, defined as where V is the shear equivalent activation volume associated with Boltzmann's constant k B . The stress and temperature dependency are based on the Eyring theory. The equivalent stress and pressure p are defined as The state parameter S, which incorporates the description of the intrinsic strain softening in a glassy polymer, can be written as a function of the equivalent plastic strain γ p : where the initial state parameter S a represents the thermo-mechanical (aging) history of the material. Intrinsic strain softening is described by a modified version of the Carreau-Yasuda function R γ : with fit parameters r 0 , r 1 and r 2 .
The elastic hardening stress is obtained using the neo-Hookean expression: where G r is the strain hardening modulus and B d = J − 2 3 F · F T . All micromechanical simulations are conducted in Marc.Mentat 2014.0.0. The EGP model is implemented in a user-defined subroutine HYPELA2. The material parameters for the polycarbonate matrix are adopted from a previous study [33] (see Table A1 Appendix A). In Table A2 Appendix A, the reference spectrum for the polycarbonate matrix is shown with 17 modes.

Results and Discussions
Using direct finite element simulations, the homogenized stress-strain responses of RVEs, generated from specific fiber orientation tensors, are obtained. The RVE is meshed with 512,000 linear elements, and the FEA of the RVE is conducted with eight CPUs and 52 GB RAM. The computational running time until the yield point depends on the loading angle, since a finer time increment is required at lower loading angles for better convergence. In general, the computational running time until the yield point is 3 days at lower loading angles and less than 1 day at higher loading angles. To investigate the effect of different RVEs with identical orientation distributions, four RVEs were generated with p 33 = 1 (see Figure 4a). The simulated stress-strain curves of these RVEs are presented in Figure 4b for loading angles of 0 • and 90 • at a strain rate of 10 −1 s −1 . The yield point is obtained as the maximum stress of the stress-strain curve, which is presented as the marker in Figure 4b. The observed variation is clearly small; for the 90 • loading angle, the stress-strain curves for the different RVEs are almost indistinguishable as they are on top of one another. In the 0 • direction, where changes in fiber placement may be expected to have a larger influence, the yield stress variation is still very small, 88.7 MPa ± 1.8 MPa, i.e., a variation of little more than 2%. These small variations between RVEs are a strong indication that the RVE size is large enough and can indeed be regarded as representative.
To investigate the influence of fiber orientation on the mechanical response, further simulations were performed at a strain rate of 10 −1 s −1 for different values of the principal fiber orientation p 33 and various loading angles. The RVEs used were presented in Figure 2. Illustrative examples of the simulated stress-strain curves are shown in Figure 5, where it is clear that the fiber orientation, as well as the loading angle both markedly influence the mechanical response. Upon a decrease in value of p 33 , a gradual decrease in the yield stress and elastic modulus is observed in the 0 • loading direction. A similar change in stress-strain response is observed when the loading angle increases gradually from 0 • to 90 • , where it should be noted that the value of the principal fiber orientation p 33 has almost no influence on the response in the 90 • loading direction. Even in the moderate case studied here, a composite with a fiber loading of 10% weight fraction, a clear orientation-induced anisotropy is observed.   Since the matrix material displays a significant dependence on the applied strain rate, the results of simulations at various strain rates on a representative volume element with p 33 = 0.94 are presented in Figure 6, for various loading angles. Some curves end earlier because their simulations were manually terminated earlier after the yield (maximum) stress was reached. Meanwhile, some simulations indeed stop earlier due to the nonconvergence, especially at lower loading angles. The resulting yield stresses at the different strain rates were determined and are presented in Figure 7a in a semi-log representation of yield stress versus the logarithm of strain rate. Figure 7a visualizes the pronounced influence of loading angle on the value of the yield stress. Additionally, the results seem to indicate that the strain rate dependence of the yield stress, here quantified as the slope of the lines (guides to the eye), also displays a variation with the loading angle. To study this effect in more detail, extra simulations were performed on RVEs with p 33 -values of 0.09, 0.47, 0.71, 0.82, 0.94 and 1.00, for a loading angle of 0 • (p σσ = p 33 ). The resulting yield stress values, at various strain rates, are presented in Figure 7b, where the response of the unfilled matrix is also added. With the larger range in fiber orientation, its influence on the strain rate dependence of the yield stress is clearly more pronounced. The trend becomes clearer when the values of the simulated strain rate dependence are ranked with respect to increasing p σσ , the fiber orientation in the direction of loading. These are presented in Table 1, which shows that, initially, up to a p σσ of 0.5, the slope only changes marginally, and stays more or less constant. At larger fiber orientations, there is a marked increase in the strain rate dependence with increasing p σσ . A similar trend is observed in the yield stress which also shows a marked increase at larger p σσ values.
The observed trends appear in good agreement with experimental observations of glass fiber-reinforced polycarbonate [10,[12][13][14]. In these experimental studies, it was demonstrated that the influence of loading angle/fiber orientation and strain rate dependence are, to a good approximation, factorizable-implying that the influence of the strain rate and that of the loading angle are actually decoupled. This was shown by plotting the yield stress data versus strain rate in a double logarithmic plot. This procedure yielded parallel lines, demonstrating that the influence of fiber orientation and that of strain rate can be separated in a multiplicative manner. This factorizability was not only observed for the low fiber weight fraction used in this study, but also for other fiber-reinforced systems with fiber weight fractions and aspect ratios over the entire range commercially available.  As shown in Figure 8a,b, presenting the double logarithmic versions of the semi-log plots presented in Figure 7a,b, a similar observation can be made. For different fiber orientations, the strain rate dependence of the simulated yield stresses, to a good approximation, is also transformed to a set of parallel lines. Remarkably, the simulated response of the unreinforced matrix seems to fit this trend well, clearly witnessing that the polymer matrix is the origin of the observed rate dependence, and that the effect of fiber-reinforcement and orientation merely shifts the response-a geometrical effect which can be identified by micromechanical simulations. In practical terms, the results presented in Figure 8 imply that the dependence of the yield stress σ y on strain rateε and loading angle θ can be expressed as For a logarithmic stress-axis, the multiplicative decoupling transforms into: log σ y (ε, θ) = log( f (ε)) + log(g(θ)), where f (ε) = σ y (ε, θ re f ) represents the strain rate dependence of the yield stress for a chosen reference state, the loading angle θ re f , which is generally chosen to be 90 • . To obtain a suitable expression for g(θ), the Hill equivalent stress, based on the well-known Hill yield function [38] is used. To simplify the expression, a plane stress state was assumed, in which case: with: where the constants R 11 , R 22 and R 33 are ratios of the yield stress in different principal material orientations to the reference yield stress σ y,re f . The value of R 23 is the ratio of the shear yield stress to the reference shear yield stress τ y,re f = σ y,re f √ 3 . In our specific case, it should be noted that the R-values are all independent of strain rate. Since strain rate and loading angle dependence are decoupled, the latter are determined at an (arbitrarily chosen) strain rate. An example of the application of the Hill equivalent stress approach is shown in Figure 9, which presents the loading angle dependence at a strain rate of 1 × 10 −1 s −1 for a principal fiber orientation p 33 of 0.55, 0.84 and 0.94. The solid lines are fits with Equations (19) and (21), taking the yield stress at 90 • as a reference yield stress. It is clear that the Hill function adequately represents the observed loading angle dependence. The Hill parameters obtained from the fit are presented in Table 2.    Figure 9 clearly demonstrates the strong influence of fiber orientation on local anisotropy; a different fiber orientation distribution clearly reflects marked changes of the load angle dependence. Of course, in a realistic situation, where a complex shaped component is injection-molded, the variation in the flow conditions over the product shape will lead to a pronounced spatial variation of the local fiber orientation throughout the product. Simultaneously, this implies that the anisotropic mechanical response will demonstrate a similar spatial variation. This variation of properties throughout the product, and its complex relation to processing conditions, product geometry, flow and mechanical characteristics of the fiber-reinforced polymer is one of the major challenges encountered in the design of load-bearing fiber-reinforced polymer components today. The usual approach is to address the problem in a two steps; firstly, a simulation of the processing step, in which the designer can obtain information on the local fiber orientation distribution throughout the component. Methods to achieve this were developed over the years [34,39,40], and are readily available in various commercial software packages dedicated to injection molding. The second step is the translation of the local fiber orientation distribution to the resulting, local, anisotropic mechanical response. This step is complex, since it usually involves either a mean-field homogenization such as Mori-Tanaka model or a full-field homogenization through the FEA of RVEs. The mean-field homogenization was able to rapidly provide the overall response and is commonly used for actual product simulations. Meanwhile, the full-field homogenization provides a better accuracy and is used in the current study. From these, the local parameters of a suitable macroscopic constitutive model must be determined, leading to a set of parameters that may vary from element to element.
The results from our simulations, however, also reveal the possibility of an alternative route. From Figure 8a, it became evident that, as mathematically defined in Equation (19), the dependence of the yield stress on strain rate and loading angle are multiplicatively decoupled. This implies that the rate dependence can be characterized on a reference orientation of choice. Figure 8b additionally evidences that the same rate dependence is actually applicable over a great variety of p σσ (p, θ), obtained from different orientation distributions tested in various loading directions. Hence, it appears interesting to investigate how the value of p σσ (p, θ) influences the value of the yield stress at a single strain rate. The result is presented in Figure 10a, for a strain rate of 10 −1 s −1 , and demonstrates that the yield stress is intimately related to the fiber orientation in loading direction p σσ . A similar observation was also experimentally made on a polycarbonate composite system with 20% weight fraction of short glass fibers, although in the experimental case, the range of fiber orientations was more limited (0.3 < p σσ < 0.75) [14]. The observed trend is described well by σ y (ε re f , p σσ ) = c p p q σσ + σ y (ε re f , p σσ,re f ), (23) whereε re f is the reference strain rate chosen (10 −1 s −1 in this case), p σσ,re f is the reference orientation chosen to determine the rate dependence (most obvious choice p σσ,re f = 0), and c p and q are parameters to be determined by fitting. As can be observed in Figure 10a, a satisfactory fit is obtained with c p = 20.8 MPa and q = 3. Equivalent to Equation (19), it is convenient to rewrite Equation (23) to where, following Equation (23): and σ y (ε, p σσ,re f ) can be described using an Eyring expression, in the isothermal case simplified to: which provides the linear relationship between the yield stress and the logarithm of strain rate that is generally observed over large ranges of strain rate and temperature [10]. On a log base 10 scale, the slope of the reference orientation can be expressed as Of course, with increasing fiber orientation, the slope dσ y (ε, p σσ )/dlog(ε) will also gradually increase, and from Equation (24), it can easily be deduced that this can be expressed as with h(p σσ ) defined by Equation (25). With all the parameters required determined from Figure 10a, Figure 10b presents a direct comparison between the measured slopes σ for different p σσ , listed in Table 1, and the trend that is predicted by Equation (28). Overall, the deviations are in the order of 5%, except for the perfectly aligned case. The results in Figure 10 suggest that, with adequate material characterization, expressions can be found that conveniently link the yield stress and its strain rate dependence to the fiber orientation in the direction of loading. Such a relation can be employed to create direct analytical links between the local fiber orientation distribution (as obtained from injection molding simulations) and the parameters required for the characterization of its anisotropic macroscopic response (Hill parameters, rate dependence).

Conclusions
In this study, we employed a straightforward micromechanical modeling approach to study the anisotropy in the strain rate dependence of the mechanical response of a short (glass) fiber-reinforced polycarbonate. The matrix response was captured using the EGPmodel, which was shown earlier to accurately represent the rate-dependent mechanical response of polycarbonate [33], and the glass fibers were modeled as linear elastic solids. RVEs were generated based on several transversely isotropic orientation distributions of short glass fibers with fixed length. The results of the simulations demonstrate that the influence of the loading direction and that of the strain rate can be regarded independently (factorizable), and that they are multiplicatively decoupled. This implies that the rate dependence of the system can be fully characterized in a single (reference) loading angle and the dependence on loading direction at a single strain rate. Moreover, the results also clearly demonstrate that the yield stress, to a good approximation, can be directly linked to the fiber orientation in the direction of loading. This leads to a new relation that uniquely links the rate dependence of the yield stress to the fiber orientation p σσ in the loading direction. We expect that the observed multiplicative decoupling can be the basis of a novel method to rapidly estimate the anisotropic mechanical response of short-fiber-reinforced thermoplastics on the basis of fiber orientation distribution. Such a method would be invaluable with respect to the development of new predictive methodologies that can take processing-induced fiber orientation into account in the very early stages of design of load-bearing components.
For future work, the micromechanical modeling in this study could be extended from an amorphous to a semicrystalline polymer matrix, since a considerable amount of polymer composites are based on semicrystalline polymers. It has been shown that the modular set-up of the EGP model provides a versatile tool that can also be successfully applied to semicrystalline polymers [41].