Thermoplastic pultrusion process: Modeling and optimal conditions for fibers impregnation

Pultrusion is a crucial method for continuous production of fiber-reinforced composites. It was developed several years ago for thermosetting polymer matrices, but the challenge is now to extend it to thermoplastic matrices, with a much higher viscosity. In this paper, we propose an analysis of the parameters influencing fiber impregnation in the conditions of this process. A semi-analytical one-dimensional axisymmetric model based on Darcy’s law and Stokes equations is developed to predict the impregnation profile inside the fibrous phase in the case of a natural impregnation governed by capillary forces. Thanks to a dimensional analysis, it enables to quantify the influence of all the material and testing parameters under the assumption of a slow variation of the geometry of the impregnation die. The cases of a straight and conical dies are discussed. For the first case, an exhaustive numerical study enables to define the optimal processing conditions for a perfect impregnation. Those results are shown to be useful tools for finding an optimal pulling velocity and die length for a given fluid/fibers pair. For the second case, we show that section reductions do not improve impregnation.


Introduction
Over the last 20 years, the investigation of thermoplastic pultrusion process has become more and more prominent in order to manufacture fiber-reinforced composites on a continuous basis and with an advanced quality. Thermosetting resins were the most commonly used matrices as they provide an easy adhesion to reinforcements and an efficicient impregnation due to their low viscosity and good wettability. Nevertheless, they tend to be replaced by thermoplastic resins owing to their advantages for toughness, impact resistance, repairability and processing speed. 1 Most of the past efforts have been directed towards the modelling of pultrusion with thermosetting polymeric matrices such as polyester and epoxy. [2][3][4] Their low-impact strenght is indeed a limiting factor for their spreading in several industrial areas.
The development of pultruded profiles has been critical to industry's ability to design stronger thermoplastic composite materials that compete against thermosets in several applications, mainly aerospace and automobile. One of the major concerns in the production of thermoplastic composites is to figure out how the resin infiltrates the network of reinforcing fibers while achieving the correct shape and preserving an excellent impact resistance. This can be explained by the high viscosity of thermoplastic melts, which is the main difficulty of the cost-effective processing of composites. 5 In the literature, several models dealing with thermoplastic resin impregnation were developed. [6][7][8][9][10][11] However, most of the past efforts remain experimentally unverified. 12 Moreover, the majority of the cited references were only interested in resin impregnation and no complete model is yet available for thermoplastic pultrusion process. Another area of investigation which has recently received attention is the powderimpregnated composites 13,14 and commingled yarns. 15 Thus, a variety of models have been employed to quantify the degree of impregnation. However, there is a lack of work dealing with the direct impregnation technique, similar to what was done with thermosetting resins. It is widely admitted that a proper selection of processing parameters can perform wetting, but this is not frequently discussed in the literature, particularly for the case of unidirectional fiber-reinforced composites. Several variables are involved in the impregnation process. The fibre bundles are pulled through the die with an imposed force or velocity. For a given geometrical configuration, if we increase the pulling velocity, the contact time decreases and impregnation can become poor. Conversely, if one decreases the pulling velocity, the process efficiency can become too low and will have no economical interest. As a consequence, a thorough understanding of the pultrusion physics is crucial for a rigorous prediction of impregnation in order to reduce the time and the cost of fabrication.
The objective of this paper is to analyze the impregnation process within the pultrusion die in a simple and efficient way. In this regard, section One-dimensional model provides a detailed theoretical analysis of the process leading to a one-dimensional (1D) model for an axisymmetric pultrusion case. The next section, after some physical considerations, presents an optimization study based on the use of dimensionless parameters, which enables a practical direct dimensioning rule for pultrusion devices. Finally, the case of a conical die is analyzed.

One-dimensional model
We consider the pultrusion problem illustrated in Figure 1. The process consists in pulling the fibers, initially placed on creels, through a die where the polymer is injected under the atmospheric pressure, P atm .W e distinguish the liquid region (') and the fibrous region (f), which is assumed to be a straight cylinder of radius R i made of aligned fibers. Fibers are pulled through the die at a constant speed V and move along the axial direction e z . During the process, the liquid matrix infiltrates the fibrous medium due to the liquid pressure and the capillary pressure, 16 leading to a full or partial impregnation at the end of the die. The fluid is assumed incompressible and Newtonian with a viscosity . Since the viscosity is assumed to be high and the velocities are rather small in such a process, Reynolds number can be considered as very low so that inertia effects can be neglected in the momentum balance of the fluid. Consequently, the fluid motion in the region (') is governed by Stokes' equation: where P and U are the pressure and the velocity of the fluid inside the liquid region ('), respectively. At the macroscopic observation scale, the fibrous medium ( f) is considered as continuous, porous and homogeneous. It is characterized by its permeability tensor which will be discussed further.

Flow in the liquid region
Let U U r , U , U z ðÞ denote the liquid velocity inside the liquid region ('). We consider the cylindrical coordinate frame. Due to the axisymmetry of the problem, the dependence on can be ignored and it can be stated that U ¼ 0. Besides, the variations of the external die radius Rz ðÞare assumed to be small with respect to the length L. Owing to such assumption, the radial velocity U r is assumed to be negligible compared to the axial component U z . Subsequently, P is assumed uniform in the thickness of the flow. Stokes' equation can therefore be reduced to the following unidimensional form: The boundary conditions for the above equations being U z Rz ðÞ ðÞ ¼ 0 and U z R i ðÞ ¼ V, the integration of equation (2) leads to: Figure 1. Illustration of impregnation mechanism within the die.
Note that U z depends on r which ranges between R i and Rz ðÞ , this expression being valid whatever Rz ðÞ . This velocity profile is the overlay of a Poiseuille flow and a shear flow due to pulling, similar to the form discussed in Agassant et al. 17 Equation (3) can then be integrated to obtain the flow rate Q along z: For simplification, equation (4) can be written as: Flow in the fibrous medium A dimensional analysis shows that, for an isotropic permeability, the ratio v z =v r inside the fibrous medium is of the order of magnitude of " ¼ R i =L, which is generally very small in pultrusion problems. The theoretical anisotropy ratio between axial and transverse permeability ratio being usually about 3 for composites, 18,19 the kinematic assumption v z % 0i s relevant here.
Let v ¼ vr , z ðÞ e r denote the fluid velocity inside fibers. The flow through the porous medium is governed by Darcy's law: [19][20][21][22] vr , z ðÞ ¼ À k @p @r ð6Þ where pr , z ðÞ is the pressure inside the fibrous medium and k is the transverse permeability. The flow rate balance inside the fibrous medium then brings: Combining (6) and (7) then leads to the following relation: where r f denotes the radial flow front position within the fibrous structure. Due to the capillary effect, the pressure at r ¼ r f is assumed to be equal to the pressure difference ÁP ¼ P atm À P c . This can be approved by Liu and Hwang, 23 when the capillary effect is dominant with respect to the local pressure. The capillary pressure is defined here 21,22 as the suction force divided by the interface liquid-air, which is positive when impregnation is spontaneous. Theoretical value of P c can be computed using Young-Laplace law 24 when the surface tension and the contact angle are known. Thus, we get the expression of the pressure as a function of z: noting that vr f , z ÀÁ ¼ @r f @t . Focusing on the stationary regime filling pattern, we take advantage from the relation z ¼ Vt. Time and axial positions are linked and it is possible to follow the pressure time history undergone by any cross-section of the fiber tow by knowing the pressure axial profile Pz ðÞ . The time derivation in equation (9) can therefore be replaced by a spatial derivation, @ @t ¼ V @ @z , which brings: Governing equations To obtain a closed form of the problem, we now consider the flow rate balance in the fluid zone ('), illustrated by Figure 2(b), which leads to: where q f 4 0 is the amount of fluid per unit length absorbed by the fiber bundle. Reminding that the velocity vR i ðÞ is negative along e r , it can be expressed as follow: Combining equations (8), (11) and (12) yields the flow rate equation: As a result, using equation (5), we deduce a coupled differential equation system describing the pressure profile and the flow front position, respectively: These two equations enable to evaluate both the pressure and the degree of impregnation along the die. The problem thus consists in solving this system of equations and studying the influence of process parameters and die geometry on impregnation, which requires numerical simulations. Note here that equation (14 2 ) can be solved in terms of y ¼ r 2 f which avoids numerical singularities when r f tends to zero. This equation remains valid till fibers are not fully impregnated. Afterwards, the flow front velocity must be set to zero, which requires a careful numerical treatment.

Dimensional analysis
In order to efficiently introduce equations (14) in finite element simulations, dimensionless parameters were adopted: Introducing also the dimensionless geometrical function Ã 1 and Ã 2 defined as equation set (14) becomes: This process shows that fluid flow is governed by four dimensionless parameters, P Ã c , 1 , 2 and 3 , the three latter being, respectively: By a simple analysis, assuming that the effect of the pressure gradient is weak in equation (3), 1 can be associated to the tensile stress inside the fiber tow. One can show that: where Ã is equal to Ã ¼ =P atm , denoting the tensile stress undergone by the fiber tow. In practice, 1 should therefore be less than Ã max ln 3 ðÞ , Ã max being the maximum normalized tensile stress. This parameter directly provides an upper bound for 1 depending on the geometrical ratio 3 . Nevertheless, in what follows, it was checked that when taking a realistic value, the maximum tensile stress was never reached in the tested configurations of the optimization study (section Optimization study).
On the other hand, 2 is expressed in terms of the medium properties and process parameters. The duration of fibers transition through the die, t p , is equal to L=V. By computing the solution of equation (10) for a constant fluid pressure Pz ðÞ¼P atm and a constant external radius R, one obtains the characteristic impregnation time t i as: Consequently, 2 can be interpreted as: It is clear that t p should be at least of the order of t i or more for a complete impregnation. This criterion can be useful in optimizing the state of impregnation for given material and die dimensions. 2 is linked to the efficiency of the impregnation process and can be used as a pre-dimensioning variable. This last equation also shows the strong importance of the capillary pressure on this coefficient. Actually, four dimensionless parameters are required.
As it will be shown in the following sections, such a simple analysis holds only in the case of a weak coupling between flow and impregnation. In general, the optimal value of 2 depends on the other dimensionless numbers and no simple result can be drawn without numerical calculations.

Straight die model
We first consider the case of a straight die model and illustrate the influence of some process and material parameters on the impregnation phenomenon. Simulations are performed with the FEM software Comsol Multiphysics 4.2 a. Governing equations (16) can be simplified here, taking into account that Ã 1 z Ã ðÞ and Ã 2 z Ã ðÞare constant, denoted as Ã 1,0 and Ã 2,0 , respectively: where the variable yz Ã ðÞ ¼r f z Ã ðÞ =R i ÀÁ 2 was introduced, and where A reference parameters set is first defined for illustration and further comparisons. We plot in Figure 3. the flow front and the pressure profiles with respect to z Ã for the following material and processing parameters: k ¼ 10 À10 m 2 , ¼ 30 Pa.s, P c ¼ 2:10 3 Pa, Inlet and outlet pressures are imposed at the atmospheric pressure, which is close to most practical industrial situations. This will be the case for all the cases presented below. This is therefore applicable to situations where the liquid polymer is simply brought into contact with fibers without additional overpressure, and where the impregnation zone is separated from the cooling zone. In the presented example, impregnation is completed at a 15% of the die length. At the entrance of the die, the pressure decreases, then we remark a little change of the pressure slope when the radial flow front location tends to zero. In such a case, a strong coupling between equations (21 1 ) and (21 2 ) can therefore be observed. Figure 4 is given here as an illustration of some of the parameters of major influence, namely the viscosity and the capillary pressure. Starting from the reference parameters given above, the fluid viscosity was varied from 10 Pa.s to 1000 Pa.s and we plot in Figure 4(a) the radial flow front position as a function of z Ã . It comes out that for low viscosities, the fluid front gets rapidly to zero. Conversely, for viscosities higher than 100 Pa.s, impregnation becomes difficult for the chosen processing conditions. In this case, it would be interesting to decrease for example the pulling speed.
Another set of simulations was carried out for different values of the capillary pressure P c (see Figure 4(b)). As expected, the higher the capillary pressure, the faster the impregnation. We can notice that in this case, when P c is lower than 10 2 Pa, impregnation almost stops for certain duration and is continued thereafter. This is in agreement with experimental results of Li et al. 24 In some particular cases, especially when the outer radius of the die is small, the coupling between die flow and Darcy flow in the fibre becomes important. In such cases, the pressure drop in the die becomes so close to the capillary pressure that the driving force for impregnation becomes low. Impregnation is therefore slowed down for a while, but as the liquid gets closer to the outlet, the pressure increases and impregnation can start again.

Effect of the dimensionless parameters
It seems clear from this analysis that the choice of an optimal set of parameters is cumbersome and makes the design of an efficient process very hard, even in the case of a simple straight geometry. The idea of the following is now to take advantage of the dimensional analysis presented in section One-dimensional model to define the best processing parameters. According to the theoretical approach developed in section Dimensional analysis, the dimensionless parameters 1 , 2 , 3 ðÞ as well as the dimensionless capillary pressure P Ã c control the impregnation process. Figure 5 is given to illustrate the influence of 1 and 3 . Results are obtained with P atm imposed at the entry and exit of the die. The capillary pressure was fixed to 10 4 Pa, its effects having been presented in Figure 4(b). We plot the radial flow front location as a function of z Ã for different values of the dimensionless variables.
As mentioned in section Dimensional analysis, increasing the value of 2 is equivalent to increasing the transition time through the die t p so that it promotes a good impregnation. The influence of the two other parameters is less straightforward.  From Figure 5(a), we notice that the quality of impregnation decreases with 1 . Roughly linked to an increase of the pulling velocity, the effect of 1 is antagonist to the one of 2 so that generally an optimum processing condition can be found. Finally, according to the influence of 3 given in Figure 5(b), we notice that the thick die promotes a good impregnation, which is not an a priori obvious result. The intensity of the coupling of equations (21) is actually higher since Ã 1,0 is very sensitive to 3 . Here again, as illustrated here when 3 ¼ 1.3, an optimal situation can be found by varying this parameter, the other ones being fixed.

Optimization study
We now present a numerical methodology to find in a systematic way parameters 1 , 2 , 3 ðÞ leading to a complete impregnation just at the exit of the die. Using the Comsol solver, a Matlab routine based on a dichotomy procedure was written. The objective function d was defined as depicted in Figure 5(b). When impregnation is complete but occurs before the end of the die, d is set to Àz Ã min , the position of the last non null value of r Ã f . When impregnation is incomplete, d is set to þr Ã min , the minimum value of r Ã f . A signed function was thus obtained. The optimization was carried out on the single parameter 2 by fixing the values of 1 and 3 . The developed routine returns a 2 which leads to a complete and optimal impregnation in the die.
Optimization results are presented in Figure 6(a) an (b) for two values of P Ã c , 0.1 and 0.02, that correspond to situations of good and poor wetting, respectively. 1 was varied between 10 À3 and 1 and 3 ranged between 1.2 and 2.5. For a given value of 3 , points under the curve lead to a partial impregnation, whereas points above the curve lead to a too early impregnation. Values of 2 strongly depend on P Ã c , so we plotted the curves in terms of 2 Â P Ã c , which represents the quantity t p =t i . Depending on the values of 1 and 3 , one checks that this ratio is generally far from being 1, as one would have expected. This is due to the non-trivial coupling between equations (21 1 ) and (21 2 ). It tends to 1 only when 1 is small and for high values of the geometrical ratio 3 . These situations correspond to a weak coupling between the die flow and the Darcy flow. For other situations, the control of processing parameters should be considered carefully because the transition time must be chosen much higher (higher values of 2 ). The effect of the capillary pressure is also confirmed to be crucial as it can be checked by comparing Figure 6(a) and (b).
To further illustrate the practical interest of this optimization study, we then present a simple method to determine a die length L and pulling velocity V for a practical situation. For a given cross-section of the die, i.e. given R i and R e , 3 is known and we immediately get the corresponding optimal values of 1 , 2 ðÞ from data of Figure 6. Then, through definitions of equation (17), we deduce the corresponding optimal velocity V and length L by the two following relations: In Figure 7, we plotted curve sets representing the optimal values of V in terms of the corresponding optimal value of L. They are given in the case of two viscosities that correspond to a low-viscosity thermoplastic polymer above its melting temperature (Figure 7(a)) and to a thermosetting resin at room temperature (Figure 7(b)). They are directly deduced from the results of Figure 6 by fixing some of the physical parameters. We chose a good wetting situation, P Ã c ¼ 0.1, a typical value of the permeability k ¼ 10 À10 m 2 and an inner radius R i ¼ 2:5 mm. The resulting curve sets for different geometries defined by the value of 3 have a similar shape for both viscosities. The maximum pulling velocity V is as expected an increasing function of the die length L, but a maximum is reached for high values of L. This is due to the increase in the coupling between Darcy flow and die flow, which results in a decrease of the fluid pressure, so a decrease of the impregnation speed. In industrial applications, one is interested in reaching a maximum pulling velocity, of the order of 1 m.min À1 . Regarding Figure 7(a), it is clear that for thermoplastic polymers, even with a rather low viscosity, a good impregnation requires huge dies more than 1 m long, which seems unrealistic. To the contrary, for a thermosetting resin, high speeds can be reached with die lengths of the order of some centimeters, which make such a process easy to develop.

Conical die case
This section focuses on the conical die illustrated in Figure 1. In practice, such reduction is necessary and eventually expected to improve impregnation. The die  . Optimal pulling velocity in terms of the optimal die length deduced from the optimization study for two practical cases: is composed of two cylindrical straight portions linked by a conical section. Material parameters and processing conditions are those used in section Straight die model, but the exit die radius was set to R 1 ¼ 13 mm while the entry radius to R 0 ¼ 26 mm. P atm was imposed at the entrance and the exit of the die. In addition, the beginning of the conical section was set to L=4 and its length was varied in order to test the effect of the angle of the cross-section reduction.
In Figure 8(a), we present a comparison of our results with a CFD calculation performed in Comsol Multiphysics 4.2 a under a 2D axisymmetric assumption. In this case, 1 ¼ 0:015 and 2 was set to 0 in order to stop impregnation. Results match perfectly which validates our flow model inside the liquid zone ('). It is worth noting that the section reduction induces a very slight pressure increase that reaches a maximum about the exit of the tapered section. Similar results were found by Kim et al. 15 when modeling the pultrusion of thermoplastic commingled yarns.
In a second step, as illustrated in Figure 8(b), the die angle was varied and its effects were observed. In this case, we took again 1 ¼ 0:015, but now with 2 ¼ 11 in order to obtain a good impregnation for the big straight die ( 3 ¼ 2:6) and a bad impregnation for the small straight die ( 3 ¼ 1:3). They were respectively denoted Straight 2 and Straight 1. As visible from Figure 8(b), the effect of the section reduction is very small. Impregnation profiles for all the die angles are close to those obtained for the big straight die. Only in the case Conical 3, the smoothest section reduction, a slight improvement is obtained. To the contrary, in the case of a sharp section reduction (Conical 1), impregnation is even worse, mainly because impregnation is made in the small cross section that does not promote a good impregnation. The pressure rise is actually very small for all conical cases, at most 1% as illustrated in Figure 8(a). Finally, we check that a section reduction is not a solution for improving the impregnation in the pultrusion process. As the fluid flow in the die is due to pulling only, there is no pressure rise due to geometrical variations. Impregnation physics is finally well-approached by the simple straight die analysis presented in sections Straight die model, Effect of the dimensionless parameters and Optimization study.

Conclusion
In this paper, we propose a simplified 1D model for the prediction of fiber impregnation in an axisymmetric pultrusion process. This model is valid when radius variations of the die are slow in the axial direction, when the radial velocity in the liquid region can be neglected compared to the pulling velocity and when the anisotropy of the permeability is not too important. The originality of our approach is to propose a strict coupling between a Stokes flow inside the purely liquid zone and a Darcy radial flow in the fibrous zone in the case of an axisymmetric shape. The present analysis leads to a set of two coupled differential equations in z that give respectively the pressure profile and the radial flow front position in the stationary regime.
Second, the dimensional analysis showed that four dimensionless parameters govern the process in the case of a straight die. The analysis of those parameters, especially 2 Â P Ã c ¼ t p =t i , gives an order of magnitude that can help to an a priori design of the die length and pulling velocity. Nevertheless, due to the non-trivial coupling between the two governing equations, no simple result can be drawn analytically. Following this observation, an optimization study was carried out to find the optimal parameters set 1 , 2 , 3 ðÞ for two given capillary pressures, and an example of process design based on those results was presented. The proposed methodology is shown to enable an optimal choice of pulling velocity and die length for a given fluid and desired final shape. Finally, it was also shown that the induced pressure raise in a section reduction was very small and that the effect on the impregnation velocity was also negligible.
From a technological point of view, we can conclude from this work that thermoplastic pultrusion is generally too slow for industrial applications by a natural impregnation of the desired shape. Other techniques such as a higher injection pressure should be employed in order to reach reasonable die lengths and pulling velocities. Nevertheless, this paper presents a guideline for die dimensioning and process parameters estimation. Such a methodology could be extended to many other geometrical situations such as flat rectangular or more complex shapes.

Funding
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.