A model for the carbon fiber breakage along the twin-screw extruder

Carbon fiber-reinforced plastics (CFRP) are widely used in the plastics indus-try. The production of CFRP is mostly carried out on the twin-screw extruder. Many different process variables change the process ability and the final fiber length in the product. In the literature, some models for fiber length decrease along the extruder with glass fibers can be found. The effects of the processing conditions on the final length for glass fiber are determined. In comparison, the incorporation behavior of the carbon fibers are not sufficiently analyzed. The purpose of this paper is to predict the fiber breakage along the screw for carbon fibers. First of all, the buckling case must be determined. The basis of the modeling is a fragmentation equation, for which a break probability and a break rate are defined. A comparison with experimental data shows satisfac-tory results.


| INTRODUCTION
The constant improvement of mechanical and physical properties of products is driving the demand for carbon fiber-reinforced plastics. The reinforcement of polymers with carbon fibers improves the stability and the stiffness. The mechanical properties are determined by the incorporation behavior of the carbon fibers in the twin screw extruder. The fiber diameter, the fiber content and the final fiber length also impact the properties. [1] However, the effects of the processing conditions on the final length in the compound are not always readily apparent and are mostly based on the experience of the compounder. A carbon fiber breakage model for the modular twin-screw extruder is therefore of particular significance for the processing industry. In the literature, many studies can be found for the fiber breakage and the influence of the processing conditions on the fiber length in the compound. Nevertheless, most studies for the fiber breakage use glass fiber as the prediction base.
Bumm et al. [2] and Shon et al. [3] propose a kinetic model for the glass fiber breakage where the fiber breakage usually expressed by the following equation: The average fiber length L is dependent on a kinetic constant k f and a length L ∞ at which breakage no longer occurs. Bumm et al. define the L ∞ by comparing the compressive force and the buckling force according to the Euler model. [4] However, the kinetic constant can only be obtained via experimental studies. The kinetic parameter varies depending on the material and the screw configurations. [2] Shon et al. do not specifically define the L ∞ , which is essential for the prediction. Ville et al., [5] based on, [3] introduce a new approach, where the fiber length degradation is described based on the specific mechanical energy (SME): By using the SME as an influencing factor, Inceoglu et al. validated the modified Shon-Liu-White model with experimental results as described in Reference [6]. The description of fiber length in the final product is accurately predicted, while the fiber breakage along the screw cannot be forecast with certainty. In general, the kinetic models do not take into account that the fiber breakage slows down for shorter fiber lengths. Furthermore, the buckling of the fibers is not considered. [6] The kinetic model in Ref. [5] was adapted for flexible vegetal fibers by Berzin et al. The exponential law for the fiber length and diameter is introduced in Ref. [7] The fiber length is dependent on the ultimate length L ∞ , the cumulative strain Γ and a kinetic constant k LΓ . Since the lignocellulose fibers also separate into individual fibers, a similar equation is established for the diameter. [7] Hernandez et al. [8] developed a pseudo-analytical model which describes the forces leading to fiber damage during flow. Investigations were carried out for l/d ratios between 50 and 300, fiber motions along its axis, perpendicular to its axis and in a shear flow at a − 45 degreeangle. The buckling criterion was modified by a factor of 2 to correct the maximum force at the center of the fiber. Furthermore, the solutions were compared to numerical results performed with the boundary element method (BEM). [8] One of the most extensive research projects on the twin-screw extruder was carried out by Durin et al. [9] In this model, a fiber position in the flow using a Jeffery equation is used. After determining the buckling parameter for each position, the breaking potential is calculated. With the help of a fragmentation matrix, a new length distribution is calculated. For the model, the buckling parameter B u has to be defined in terms of the flow direction. A similar model for fiber length reduction in injection-molded composites was developed by Phelps et al. [10] The model is based on a conservation equation for total fiber length combined with a breakage rate due to buckling under hydrodynamic forces. Furthermore, the model was implemented in a mold filling simulation. [10] From the state of the art, different needs for research on modeling on the twin-screw extruder are derived: 1. Research need 1 (RN1): The assumptions and simplifications used here are based on the processing and knowledge of glass fibers. The transfer of the model assumptions for glass fiber breakage to carbon fibers is a significant step in order to use the present models for other fiber types. An adjustment of the assumptions has to be made in case of non-compliance. 2. Research need 2 (RN2): The geometry of the twinscrew extruder is special and offers different approaches for modeling. In particular, the modular design and process areas (gusset area, gap area and channel area) require special consideration during modeling. Especially their inclusion during the fiber length calculation is essential. 3. Research need 3 (RN3): Up to now, shear rates have only been used to calculate the fiber length for selected parts of the twin-screw extruder. [6,9] calculate the shear rate against the channel height. However, 3D simulations have shown in preliminary investigations that the other geometric areas also have a considerable influence. The aim should be to calculate all shear rates and influencing factors such as viscosity as well as the linkage with RN2. 4. Research need 4 (RN4): The influence of parameter changes on the fiber length for carbon fibers has not been investigated. In particular, the adjustment of the final fiber length during the change is to be analyzed. Furthermore, the determination of the effects of the individual factors on the final length is particularly important.
The main aim of the present study is the development of an analytical model for the carbon fiber breakage for twin-screw extrusion. The focus is on the transmission of the assumptions of glass fibers to carbon fibers. Therefore, the buckling case for the fibers has to be determined. Furthermore, the tensile stretch of the fiber will also be calculated in the breakage probability.

| DETERMINATION OF THE BUCKLING CASE
Kloke developed a model for the length distribution of short glass fiber compounds after processing on the twinscrew extruder. [11] The fibers in a polymer flow are charged by a shear stress τ resulting from the shear flow ( Figure 1). [11] The shear stress can cause different loading cases [11] (Figure 2): 1. The force on the cross-sectional surfaces of the fiber is so high that the fibers are crushed under this load. 2. Due to the force, the fiber buckles. In the case of fiber breakage due to Euler buckling, a distinction must be made between plastic and elastic buckling.
The theory of Euler is based on Hooke's law, which means that the buckling formula is only valid above the slenderness limit or below the proportionality limit according to Reference [12]. The proportionality limit indicates the stress above which the proportionality of Hooke's law of stress is no longer valid. Since carbon fibers (also glass fibers) are brittle materials, the progression of their stress-strain curves is approximately linear.
To determine the buckling case, the limiting slenderness λ 0 is calculated as follows [12] : where E: E-modulus; σ: tensile stress.
In addition, the critical fiber length l c must be considered. According to Euler, elastic buckling takes place above this critical length and plastic buckling below it. [11] l c is defined as follows [12] : where A = π Á d 2 À Á 2 ; I: axial moment of inertia; A: crosssectional area; λ 0 : limiting slenderness; d: fiber diameter.
In the following table, the values from the literature and the calculated results for the critical length as well as the limiting slenderness for chopped fibers can be seen. The limiting slenderness and the critical length are calculated using Equation (4) as well as 5 and the values are highlighted in gray in Table 1. Two types of carbon fiber were used for this calculation, high tenacity (HT) and intermediate modulus (IM).
With an assumed diameter of 5 or 7 μm and the calculated limiting slenderness, plastic buckling is only present below a value of 42 or 27 μm. Elastic buckling can therefore be assumed for the fibers. For the elastic bending force line, the length-diameter ratio is inserted into the Eulerian buckling force. [11] For different l/d-ratio, F Euler (equation (6)) [11] is calculated and can be seen in Figure 3.
where d: fiber diameter; l: fiber length; E: E-modulus. In addition to plastic and elastic buckling, another loading case can occur. The compressive force on the cross-sectional surfaces of the fiber is so high that the fibers are crushed under this load. [11] The compressive force F (until breaking) results from [14] : The calculations for the HT and IM fibers are summarized in the following diagram.
The vertical lines indicate the l/d ratio for the calculated limiting slenderness (cf. Table 1). The critical length is divided by the diameter. Below this limit, plastic buckling occurs. Elastic buckling can be assumed for both fiber types from a l/d ratio of approx. 5.9 (HT) or 5.5 (IM). Transferred to a fiber length, plastic buckling would only occur below a fiber length of 42 (HT) or 27 μm (IM). In first experimental results, the end fiber length is significantly higher, so that the case of plastic buckling can be excluded. [11] For both types of carbon fibers, the compressive force for crushing the fibers is not dependent on the length and is around 1.6 (HT) or 1.1 (IM). In case of fiber breakage by crushing, the point of intersection of the two forces (bending force and crushing) must be considered. The value for crushing is at a l/d ratio below 2.5. If the fibers are around this ratio, the compressive force is smaller or equal to the bending force. Only below this fiber length does the compressive force become relevant for fiber breakage. Transferred to a fiber length, crushing would only occur below a fiber length of 14 (HT) or 8.5 μm (IM). For these carbon fibers, crushing can also be excluded. Therefore, a fiber breakage according to the second Euler buckling can be assumed. [11] 3 | THEORETICAL MODEL To describe the fiber length development, two constants have to be determined [15] : The first constant describes the fiber length decrease caused by the twin-screw extruder. This rate includes all mechanical processes that lead to a breakage. Considering a time interval, this rate expresses the probability that a fiber will break in a given time period. The second rate considers the increase of the different fiber classes by breaking the longer fibers down into shorter fiber lengths. This constant describes the rate at which fibers decrease from one length to another. Viewed over a period of time, it describes the probability of fractures during this period. For the description of fiber degradation, a differential equation is needed, which considers the decrease and the increase of fiber length classes. Length intervals K are defined for this purpose: where i: discrete points; l: fiber length. In Figure 4, the intervals along the fiber are visualized. In this case, the first class is K 1 = [0, Δl] and the last class is K n = [(n − 1)Δl, nΔl], with the initial length of each fiber L = nΔl. The fragmentation equation (9) is given by the following equation based on. [10] The number of fibers N(t,i) in length class K i at time t. The breakage rate is later evaluated with the mean fiber length of the corresponding class (cf. Equation (8) For the solution of the fragmentation equation, the fiber breakage rate k and the probability density p l( j) are needed.

| Fiber breakage rate k
The model in this article is based on the concept that a flow force induces fiber buckling and, with the increasing bending load, a fiber breakage occurs. To break a fiber, the bending load must be at least equal to or greater than the load a fiber can withstand. When a fiber bends, the tensile stress causes the fiber to elongate and the compression stress results in a shortening of the fiber. In the middle of the cross section is the neutral axis, which maintains the length during loading, see Figure 5. [16] The bending load can be calculated with the help of the acting bending moment M, the axial moment of inertia I and the maximum axis distance a max [17] : In this equation, a max describes the distance from the neutral axis to the edge. In the case of a round bar, a max corresponds to d/2 (see Figure 5). Inserting the moment of inertia I = πÁd 4 64 according to Reference [12], a max and converting the bending moment changes the formula as follows: where d: fiber diameter; y max : maximum deflection of the fiber; F flow : force acting on the fiber by the flow. In this equation, y max describes the shift of the neutral axis due to the force from the polymer flow (see Figure 5).
For the calculation of the bending moment in Equation (11), the deflection caused by the force (y max ) and the force itself (F flow ) have to be defined. During processing on the twin-screw extruder, the polymer flow creates a load on the fiber and a substitute force F flow from the shear stress of the flow σ flow and the lateral surface A surface can be determined (Equation 12). [11,17] The polymer flow generates a shear stress on the fiber surface resulting from the shear flow. The shear stress acts on the longitudinal direction of the fiber, for which the sheath surface of the fiber is used (see Figure 1).
The stress can be described with the shear rate and the viscosity of the material η and with some transitions the force is calculated as follows [11] : The maximum deflection of the fiber is described by the following formula [17] : Length intervals K along a fiber F I G U R E 5 Deflection of the fiber due to the force F Euler : Force to bend the fiber according to Euler. In addition to the flow force, the force that a fiber can withstand is calculated based on the Euler buckling theory (see Equation (6)). Inserting F Euler and F flow in Equation (14) results in: With Equations (13) and (15), the bending load caused by the polymer flow σ edge is described as follows: To break a fiber, the bending load by the flow (σ edge ) must be at least equal to or greater than the tension a fiber can withstand. The breaking stress of a fiber is defined in Equation (17) as in Reference [12].
where ε: elongation at break; E: E-modulus. Balancing the two tensions, the following equation is obtained: The fiber breakage rate k for the fragmentation equation (Equation (9)) is defined with the help of the Equation (18) including the residence time t l . For a time interval, k expresses the probability that a fiber will break in a given period of time between the considered element sections.
The rate includes all mechanical processes that lead to a break and can be described as follows: The following graph shows the breakage rate for a fixed point in time as a function of the fiber length ( Figure 6).

| Probability density
For solving the fragmentation equation, in addition to the fiber breakage rate k, a probability density p l( j) is needed. For the probability density, a distribution which is symmetric around the expected value is required. Durin et al. and Phelps et al. use a Weibull or a normal distribution which is valid for all real numbers. A standardization is performed in order to adapt the probability density to the fiber breakage probability. In contrast to a Weibull distribution as in Reference [9], beta distributions can be used for local fractional probabilities, which are defined on arbitrarily closed intervals. The beta distribution is a continuous probability distribution over an interval. This allows a determination of any real number in the interval range without a standardization. The length of the interval is given by the respective carbon fiber length. [18] p y where Γ: Gamma function ðΓ z ð Þ = Ð ∞ 0 t z − 1 e − t dt). [18] With this Equation (21), [18] for each length y [0, L], it is possible to specify a beta distribution on [0, y]. The parameter α of the beta distribution has to be defined and is determined by existing data. The density function changes depending on how alpha α is defined. With the help of the experimental data, the parameter is determined and the density function is fitted. It is specified with the Nelder-Meader method and amounts to 2.06. In this case, the variance is 0.78125 and a standard deviation of σ = 22.1%. In this context, for a fiber with an initial F I G U R E 6 Breakage rate for 5 s for a screw speed of 200 1/min length of 4 mm, it is more likely to break at 2 mm with a standard deviation of approximately 884 μm than at the edge points.

| VERIFICATION
To solve the fragmentation equation, a code was written in MATLAB®. For the fracture rate k, the residence time, the viscosity and shear rate values are taken from the simulation program SIGMA (simulation of co-rotating twin-screw extruders). The shear rate is calculated separately in the twin-screw geometric areas channel, gap and intermeshing area, see Figure 7. The differential Equation (9) is solved using the classical Runge-Kutta method. After all iteration steps have been completed, a matrix N is generated for discrete points in time of the fiber length distribution, which can be plotted as a bar chart.
To check the general function of the presented calculation model, the analytical model is first checked with previously measured data. The first calculations and measurements were made for a speed of 200 1/min and a fiber content of 35%. For the investigations, a twin-screw extruder (KraussMaffei Group GmbH, ZE28 BluePower Ultraglide, Germany) was used. For the materials, a polypropylene (Borealis AG; HK060AE, Austria) and chopped fibers with a length of 4000 μm (C.A.R. FiberTec GmbH, Germany) were chosen. The screw and barrel design can be found in Figure 8. The positions for the sampling are marked in the screw design with blue arrows. While the first position is directly after the fiber feeding, the last position is at the screw tip. [19] In Figure 9, the simulated graph and the measurement points are shown. The fiber feeding starts at 24.95 s and the screw tip is reached at 54.25 s. The measurement points are close to the simulated data, and significant deviations cannot be found. The fiber end length of 54.25 s matches the simulated value well. The initial fiber reduction in the first 5 s or so generally corresponds with the simulated values, but slightly higher deviations can be seen in comparison to the other points. The first measured points are located shortly after the fiber entry into the melt. Some experimental investigations regarding homogeneity over fiber content have shown that these points have the lowest values. The sampling position and the relatively long fiber length could be possible reasons for the deviations.

| CONCLUSION
In this article, a model for the carbon fiber breakage along the twin-screw extruder was introduced. For this purpose, a fragmentation equation was set up and solved numerically. A further possibility for describing probabilities and breakage rates was presented. The parameter α of the beta distribution was defined with experimental data for carbon fibers. The buckling case was also determined using the values for carbon fibers. In general, this model can be applied to other fibers by adjusting the beta distribution and reviewing the buckling case. The first experiments show good agreement with the resulting length development. For a validation of the model, more data should be collected to increase the informative value. The time steps for the length calculation in the MATLAB code can also be made smaller to obtain the best possible results, even if the computing time increases significantly.