Review of the Numerical Modeling of Compression Molding of Sheet Molding Compound

: A review of the numerical modeling of the compression molding of the sheet molding compound (SMC) is presented. The focus of this review is the practical difﬁculties of modeling cases with high ﬁber content, an area in which there is relatively little documented work. In these cases, the prediction of the ﬂows become intricate due to several reasons, mainly the complex rheology of the compound and large temperature gradients, but also the orientation of ﬁbers and the micromechanics of the interactions between the ﬂuid and the ﬁbers play major roles. The details of this during moldings are discussed. Special attention is given to the impact on viscosity from the high ﬁber volume fraction, and the various models for this. One additional area of interest is the modeling of the ﬁber orientation.


Introduction
Fiber-reinforced polymer (FRP) composites are becoming increasingly attractive in areas where properties such as high strength, stiffness, and crash resistance are sought for as well as low weight, for instance in the automotive industry.
Glass or carbon fibers are often used as reinforcing materials. These fibers may be arranged in tows and weaved or stitched together to form fabrics. Different types of resins are used as matrices, including thermosets such as polyester, vinyl ester and epoxy, and thermoplastics such as nylon.
In the current review, the flow during compression molding of the sheet molding compound (SMC) is in focus. The main advantage with the compression molding of SMC is that the process is relatively fast as compared to many other processes to manufacture FRP composites and that the tooling is relatively cheap compared to metal forming [1]. Raw SMC material is also comparatively cheap compared to continuous fiber prepregs. SMC sheets are produced by spreading chopped fibers on a resin, which in most cases is a thermoset resin. The resin is then pressed between rollers to form sheets of SMC material. At room temperature, these sheets have a rubber-like appearance. A consequence of the random orientation of the chopped fibers in SMC is that the mechanical properties are not as good in comparison with continuous aligned fibers. However, if the fibers can be aligned in desired directions during the molding process, the mechanical properties can be improved extensively. For the process discussed here, pre-impregnated sheets of fibers and matrix are cut into parts of an appropriate size and placed inside a heated mold so that they cover an appropriate part of the tool, this can be seen in Figure 1. The mold is then closed so that the charge starts to flow and form the desired shape. The compression is then held for an appropriate time for the resin to be cured and the finished part can then be removed from the mold.  One major issue with SMC parts, however, is that the quality of the finished parts is very sensitive to variations in the process parameters [2], which means that the properties may be difficult to predict with quality and trust. The chopped fibers that initially are orientated in, more or less, random directions may orientate during flow and the final part may get anisotropic mechanical properties as a result.
During the compression, the charge will flow inside the mold. The exact flow pattern will be determined by: • The properties of the prefabricated charge • The properties of the molten charge • The geometry of the mold • The placement of the charge • The temperature of the charge and the mold • The closing speed of the mold • The tilting of the tools as they are closed Hence numerical modeling of the molding process is complicated; the rheology of the molten compound is complex, there are significant temperature gradients, as well as three-dimensional effects in the flow front, and the geometries of the molds are often quite complex. There exist several commercial software packages, such as Autodesk Moldflow, Moldex3D and 3DTimon, that can be used to model the process. It should, however, be noted that all of the above were initially developed for injection molding, and there is not a significant amount of literature on the validation of these platforms with regards to compression molding of SMC [3].
Noteworthy early work on SMC was done by Barone and Caulk; including heat conduction in SMC [4,5], models of the flow [6,7], and experimental visualization of the flow [8]. Similar experiments have since been conducted by Olsson et al. [9]. An important conclusion drawn by Olsson et al. [9] is that a higher mold closing speed resulted in a more homogeneous flow, while a lower closing speed resulted in more mixing between the layers of the charge. Important early numerical work was performed by Tucker et al. [10] and Lee et al. [11], who developed a model for the flow of the charge based on Hele-Shaw flow (see for instance [12]). Significant work has also been done by Advani [13][14][15][16], especially with regards to the modeling of fiber orientation, and by Dumont and Le Corre, including both numerical modeling of fiber orientation [17,18] and experimental models of the rheology of SMC [19][20][21][22][23]. Some more recent efforts include Li et al. [3], Oter et al. [24], Song et al. [25], and Motaghi [26]. Several relatively new studies have also concerned the modeling of the entire process chain [27,28], concerning both process modeling [29,30] and mechanical modeling [31,32].
The focus of this review will be on the flow during the molding, meaning that aspects such as preparing of SMC or curing will not be considered.
The first step is a more detailed description of the compression molding process, Section 2. Section 3 then contains a more general description of the rheology during the molding, followed by a more in-depth discussion regarding viscosity modeling and fiber orientation in Sections 4 and 5 respectively. In Section 6 some additional newer initiatives involving interesting alternative methods for numerical simulation of the process are briefly discussed. Finally, Section 7 provides some concluding remarks.

Phenomenological Description of the Process
As briefly described above, during traditional compression molding of SMC, sheets of a rubber-like compound manufactured from a thermosetting resin, chopped fiber bundles, and possible some filler are loaded into a heated mold, the mold is closed and the charge will fill the mold. During this process, the compound is subjected to a mechanical load and heat. This will result in a pressure gradient and a temperature gradient within the compound. The temperature gradient will cause a gradient in viscosity while the pressure gradient will compress the compound (that is air pockets and voids) and initiates motion. Due to the large temperature and viscosity gradients, the flow will not be even, especially not at the start of the pressing. Instead, the compound located closest to the tooling will move first. This is shown in [33] for cases where stacks of 5 circular sheets with a diameter of 100 mm were placed between two parallel heated plates, after which the upper plate was set in motion towards the lower plate at a constant speed.
Different speeds and temperatures of the plates were used and the uneven flow was captured by video recordings of the flow front as shown in Figure 2 for a relatively low speed of 2 mm/s, an lower plate temperature of 135 • C and an upper plate temperature of 165 • C. It should be noted that the compound starts to flow near the lower plate as it had been subjected to heat for a longer time. This phenomenon was first observed by Barone [8] from partial moldings and later confirmed by molding using sheets of varying color by Olsson [9], where it is obvious that the outer layers end up at the flow front. Hence, the outer layers do not stay as outer layers which results in a complicated three-dimensional flow in the front of the SMC during pressing. From Figure 2 it can also be observed how the fibers with an averaged length of 26 mm seem to dominate the flow and there may be air entrapment at the flow front. Another observation in [33] is that due to the high temperatures the resin may start to boil in areas with low pressures creating voids. For an open mold as the one presented in [33], this will happen at the flow front, while for a closed mold it may still happen that the pressure from the compression does not reach all areas within the mold. This may, for instance, be due to the geometry of the mold, partly cured resin, or that the fibers carry the load.
(a) (b) Figure 2. Flow front progression for a relatively low speed of 2 mm/s, a lower plate temperature of 135 • C and an upper plate temperature of 165 • C [33]. The charge bulk of fibers and resin can be seen as the gray area in the left image (a). It can be seen that the flow starts closer to the mold walls before the bulk starts to move (b). These are previously unpublished photos.
In addition to boiling there are often also voids in the compound; air may be entrapped between the layers of the compound during lay-up and air may be entrapped at the flow front during compression, see Figure 2. Larger amounts of voids may influence the pressure build-up, the viscosity distribution and the flow of the compound. The transport of voids during compression was studied by Sentis [34] using 3D synchrotron X-ray for compounds with weight fiber contents of 29% and 50%. It was shown that the initial volume void contents of 1.5%-5.5% and 23% for the respective cases were reduced to nearly zero during the compression. Two mechanisms for the removal of the voids were observed: a transport out of the compound as described in detail in [35,36], and a dissolution into the resin as also discussed in [37,38] for liquid molding processes. To exemplify by using model materials and Particle Image Velocimetry Westerberg [36] show that voids may move faster than the suspension towards the fluid flow front during pressing.
The geometry of the mold is another important aspect of the compression molding of SMC. To exemplify the manufactured parts are often designed with stiffeners such as ribs of different shapes; to make the compound with relatively long fiber bundles flow into these can be challenging.

Rheological Measurements of SMC
The rheology of the SMC during the molding process is complicated for a number of reasons, among them are: • High fiber content • The relatively long fibers • The orientation of the fibers • The fibers are flexible Most of these points (1-4) will be discussed in Sections 4 and 5. Displacement of air and material variations can mostly be considered as things that should be avoided rather than measured or modeled. Anisotropic viscosity is mainly a result of the fiber orientation, this will be discussed briefly in Section 4. Integrated "black-box" values of rheological properties of the SMC may, however, also be captured from experiments, although a significant issue is that ordinary rheometers usually cannot be used due to the long fibers. Some early experimental studies were conducted by Lee et al. [39]. The viscous properties of glass fiber SMC were studied for varying temperatures and average fiber lengths. A simple model for the viscosity of SMC between two parallel plates was suggested with L as the length of the plates, H the separation between them, V f being the volume fraction of fibers and η the basic viscosity of the medium. The authors however stated that this model would probably not be valid for more complex charges and molds. Vahlund and Gebart [40] used a modified press that allowed measurement of rheology in large tools; by using a set of pressure sensors the pressure in the charge was measured for a range of closing speeds ranging from 3 to 30 mm/s and a similarly wide range of press forces. Le Corre et al. [21] also developed a rheometer in a larger scale to be able to properly measure the properties, the results of this will be discussed in Section 4.
A rheological aspect which has more recently been suggested by experiments by Ferré-Sentis et al. [34,41] is that considering the SMC charge as incompressible might be an unrealistic assumption. This is likely connected to air both inside the sheets and between them. Hohberg et al. [29] did similar experiments as Ferré-Sentis and developed a model for how a compressible SMC affects its rheological properties, with reasonably good agreement with experimental results.

Models for the Viscosity
A vast amount of research has been carried out on the relations between solid fraction in suspension and viscosity, see [42][43][44][45] for instance. To exemplify Einstein [46,47] and later Batchelor and co-authors in a series of papers [48][49][50][51] developed models from first principles that are useable for low fractions of monodispersed particles. Einsteins model, for instance, may be written as while the Batchelors model, assuming a Brownian suspension in a Newtonian liquid, can be written as In both these expressions, φ is the volume fraction of particles and it is assumed that the suspension is dilute. For larger solid fractions empirical models have been derived, for example from Ferrini et al. [52].
Work on the extensional viscosity of semi-dilute suspensions of elongated particles was also presented by Batchelor [53], Dihn and Armstrong [54], and Acrivos and Shaqfeh [55].
However, these models become problematic when dealing with prolonged particles. One reason is that the number of possible contacts between particles generally increases with the aspect ratio of particles [56,57] and that orientation and torque of the particles become important factors [58,59].
Several semi-empirical relations for how various factors affect the viscosity of the SMC charge have also been suggested. It was found by Le Corre et al. [21] that the viscosity mainly depends on three factors: strain rate, temperature and fiber volume fraction. The dependence on the strain rate can be described using a Carreau viscosity model as established by Lee et al. [39] where λ is a time constant,γ is the shear rate, and n is a power-law index. Dependence on fiber volume fraction and temperature was also studied experimentally by Le Corre et al. [21] and Dumont et al. [19,20]. They suggested that the fiber fraction dependence could be described as where the subscript 0 indicates initial conditions, η Paste is the temperature-dependent viscosity, f f is the fiber volume fraction,˙ is the strain rate and n is a power law index. The temperature dependence was described as in which η 0 is the initial viscosity of the charge and T is the temperature. A model based on Equation (6) was implemented numerically in Marjavaara et al. [60] and expanded with Equations (4) and (5) in Kluge et al. [61]. In both these works, the goal was to numerically determine the pressure field, and those models were then compared to experimental data with relatively good agreement. A subject that has been discussed is whether the viscosity model should take into account differences between the bulk of the flow and the layers closest to the mold walls. Kluge et al. [61], for instance, evaluated a number of different models with different formulations both with the same equation for viscosity through the charge and with a separate equation for the viscosity of the skin closest to the walls and for the cases studied it was concluded that a one equation model was sufficient to predict the pressure in the mold.
More recently a model was suggested by Bertóti and Böhlke [62], in which the viscosity of the flow was shown to be affected by the orientation of the fibers. A viscosity tensor was introduced that takes into account the fiber orientation, the shape of the charge and the material properties of the charge. The fiber orientation was modeled using the orientation model suggested by Advani and Tucker [13]. This new viscosity model does not consider temperature, fiber volume fraction or curing; however, the authors outline the procedure to include this.
Favaloro et al. [63] also studied the connection between viscosity and fiber orientation and derived an expression for a term referred to as Informed Isotropic Viscosity where η m is the viscosity of the matrix, κ 23 is a factor connecting shear to viscosity, R η is a term describing the anisotropy of the viscosity (that will be equal to 1 for an isotropic viscosity), andγ is the magnitude of the strain rate. A and D r are the fiber orientation tensor and a tensor connected to the evolution of fiber orientation respectively, both of these terms will be discussed further in Section 5. Tseng and Favaloro [64] implemented Equation (7) for injection molding procedures and obtained good agreement with experimental results.

Modeling of Fiber Orientation
Most of the work relating to the modeling of fiber orientation in the compression molding of SMC is based on the work of Jeffery [65]. Jeffery extended Einstein's work [46,47] and derived expressions for the motion of a single ellipsoidal particle immersed in a fluid by solving the equations of motion.
Folgar and Tucker [66] adapted Jeffery's work specifically for fibers in concentrated suspensions. This was done by adding a term, called a rotary diffusion term, to account for the interaction between the fibers. A problem with their initial model, however, is that it involves the probability distribution function (PDF) of the angles of the fibers, which is a more realistic case results in a prohibitively expensive computational load.
Advani and Tucker [13] suggested a solution to this problem by representing the PDF as a set of tensors, which results in an expression for the motion of a single fibeṙ and the equation for the evolution of the second-order orientation tensor, where ω is the vorticity tensor,γ is the strain rate tensor, a ij is the second-order orientation tensor, λ is a parameter describing the shape of the particle, p is a unit vector describing the orientation of a fiber with p i being functions of the fiber angles, ψ is the PDF, ∝ψ ∝p is the gradient operator on the surface of the unit sphere created by p, δ ij is the Kronecker delta, and α is a constant connected to whether planar or three-dimensional orientation is calculated. D r is the previously mentioned rotary diffusion term, which in the case suggested by Folgar and Tucker [66] depends on the shear rate and the constant C I , which is called the interaction coefficient. The value of C I must usually be determined by experiments and by letting C I vary it is possible to adapt Equation (8) to many different cases. The approach in [13] results in equations that do not have a closed form as an orientation tensor of the higher-order will appear in every term describing the change of the orientation tensor. A closure approximation, i.e., a way to approximate the equations in a closed form, is thus required; this can be done in a variety of ways [15].
One important difference between the Advani model (Equation (8)) and experimental results is that the fibers align more slowly than the model predicts [67,68]. Several amendments to this model have therefore been suggested. Wang et al. [67] suggested the Reduced Strain Closure (RSC) model, in which the diffusion term is multiplied by a constant term κ, and the closure approximation is modified by a term depending on κ and the current fiber angles. In the case where κ is equal to 1, Equation (8) is recovered. An alternative approach was suggested by Phelps and Tucker [69], which they labeled the ARD-RSC (Anisotropic Rotary Diffusion-Reduced Strain Closure). They suggested that Hand's tensor [70] could be used to connect the diffusion closer to the fiber orientation. To do this five experimental parameters are required. According to Tseng et al. [71], this may result in issues with the numerical stability of the obtained solutions. They thus formulated their own model, called iARD-RPR (Improved Anisotropic Rotary Diffusion-Retardant Principal Rate), which has four constants and fewer terms overall than the ARD-RSC model; however, a critical issue with this formulation is that it is dependent on the coordinate frame used. Tseng et al corrected this issue in a later paper [72]. Worth noting is that parameters in all these models might not be readily accessible in more realistic cases.
A common theme among these fiber orientation models is that they are phenomenological, which means that they might be very good for specific processes to which they have been tuned, but not necessarily good for a more general case. Another common problem is that they are based on an assumption of Newtonian flow which as previously mentioned is a somewhat unrealistic assumption. Dumont et al. [17] numerically studied fiber evolution using both Jeffery and Folgar-Tucker with various different closure approximations, as well as a more direct method for directly modeling the fiber-fiber interactions. They found that the model that best modeled the fiber motion in a more general case is the Jeffery equation; it was also stated that the tested phenomenological models did not perform well in more general cases.
Another detail is that the fibers tend to be regarded either as rigid ellipsoids or as rigid rods, while in reality they can be bent, both during the manufacturing of the sheets and because of the flow during molding. One method has been suggested is the Beam-Rod model [73][74][75], where each fiber is modeled as two rigid segments connected by a bead. This formulation is then used to extend the Folgar-Tucker approach. Another similar possible method for addressing this issue is Direct Fiber Simulation suggested by Schmid [76] and Lindström [77] where the fibers are also modeled using the rod and bead formulation. The difference is that this method does not involve the Folgar-Tucker based method involving the PDF.
An alternative model has been suggested [78,79] where the fiber network is modeled as a deformable porous media. There is some work done regarding this with regards to injection molding processes (see for instance [80]), but to the authors best knowledge this has not been successfully applied to compression molding of discontinuous fibers.

Possible Areas of Future Research
There are several areas of research that are of interest, especially with regards to the issue of fiber orientation since this is arguably the greatest challenge associated with numerical modeling.
The previously mentioned method where the fiber network is modeled as a deformable porous medium might be an interesting alternative to the traditional fiber orientation models. Another method was proposed by Oter et al. [24], where lubrication theory is used to describe the flow, and the fiber orientation is described without using equations based on Jeffery's equation. This has still only been used for a two-dimensional problem, but suggestions were made for generalizing their method to more complex geometries. Yet another was presented by Schommer et al. [81], who described the charge using an Arbitrary Lagrangian-Eulerian formulation. Another possible method is Discrete Element Modeling (DEM) (see for instance [82]). By modeling the fibers as composite particles, i.e., built as sequences of lesser particles, both fiber bending and breakage could be modeled.
An associated issue with any model, suggested above or otherwise, is that models should also preferably be applicable in an industrial context, which imposes limits on the complexity and time consumption. Connected to this is the issue of predicting the mechanical properties of the finalized SMC parts. This will not be discussed in detail in this paper, but the approach of Görthofer et al. [27] allowed the mapping between results of compression molding simulations into structural models with good accuracy regarding predictions of structural properties. Until new models for the fiber orientation have been developed a possible route forward is to use commercial software as those presented in the introduction and calibrate them with experiments.

Concluding Remarks
A review regarding the practical difficulties of numerical modeling of the compression molding of SMC has been presented. Several key issues such as modeling of fiber orientation and modeling of viscosity have been discussed, as well as approaches for dealing with these issues. One conclusion that can be drawn is that there is yet room for significant future work. This is especially true with regard to the numerical modeling of fiber orientation. Funding: This work is a part of the PROSICOMP project, funded by the Swedish Research Innovation Agency (VINNOVA) and the industrial partners.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: SMC Sheet molding compound FRP Fiber-reinforced polymer