Le paramètre α 3 utilisé pour décrire le comportement compressible enrichit la formulation incompressible du modèle Dumont proposé en 2003. Le comportement compressible de SMC a été caractérisé dans le cadre de la collaboration avec le laboratoire 3SR à Grenoble-France. Toutes les équations liées au comportement compressible utilisé dans ce travail sont une preuve des communications exprimées le long de ce doctorat afin de proposer un modèle rhéologique SMC compressible macroscopique.

Introduction

Since the decades of the 50s automotive parts has been created using fiber reinforced composites, bringing a whole novel technology for the industry. Aluminum semi-structural parts presented in cars have been replaced by their counterpart made of composites. While replacing for reinforced composites those parts made of aluminum, a reduction in the production cost and the weight of the final part is granted. Such market has gained ground along the years and we find nowadays many producers of automotive components with those reinforced composites. The main challenge of this technique yields on guaranteeing good mechanical behavior while respecting the quality in the design of complexes geometries. The amount of fiber within the composite is proportional to the mechanical exigence that the part undergoes. However, fibers induces complex mechanical flows complicating the process to obtain a suitable part.

Another fact that motivates the use of fiber-reinforced composites parts is the environmental regulations restricting automotive pollution. European laws demand the decreasing of CO 2 emissions emitted by cars for the year of 2020. By looking at the emissions emitted by car in 2012 (Figure 1) the CO 2 level permitted per weight is going to be reduced by a fix amount of 15g/km. For an average car weight of 1400Kg the emissions records in the range of 135g/km. Thus, the reduction of 15g/km implies a reduction of 11% in emissions forcing the producers to improve their fabrication process and provide better solutions.

There are different strategies to reduce CO 2 emissions in a car (Fig. 2). From the mechanical efficiency of the car, one gets the improvement on the motor efficiency, driven chain or the reduction of tire resistance against the road. Any of these solutions are directly translated to a reduction of power to complete 100 km. Specifically for the purposes of this work, another straightforward manner of decreasing emissions consists in lightening the car. Producers of automotive parts are needed to provide with an equivalent automotive component with a remarkably cut down of the weight. For that, specific steel parts will be replaced for its counterparts made of fiber reinforced composites.

So far semi-structural automotive parts made of composites are presented in the market. Nevertheless, in order to fulfill the European regulations the producers must extend their approaches to structural parts. Sheet Molding Compounds materials (SMC) reinforced with glass fiber arises as one possible composite used for this purpose.

As pictured in Figure 3 engineers has already study the feasibility of modifying the fabrication process of structural and semi-structural auto parts gaining up to 110Kg in weight. This implies a decreasing of 10gr of CO2 per 100km. Nevertheless, fiber reinforced composites parts call for challenging mechanics and behavior since fibers reinforcements come together with non homogeneous properties.

The main goal of this work is to model the industrial process when creating a car part of any kind shown in the previous figure. Assisting the engineering task of predicting the mechanical response of the molding process to ensure the quality of the part. The pieces in which we can save the most of weight are submitted to structural loading and engineers need to guarantee the good mechanical resistance of the reinforced part. The fabrication of the automotive parts are done by compression molding process explained in the next section.

Compression molding cycle

Compression molding is a well known manufacturing technique for composites parts. This process has undergone a remarkable growth since the 1950s during the development of Sheet Molding Compounds (SMC) for the automotive industry.

In compression molding process, the amount of material required is placed between the upper and lower molds. During compression, the mold cavity at approximately 150 o C is filled by the flowing of the SMC plies until a chemical reaction occurs. After forming, the final shape is released from the press and the thermal equilibrium is obtained by the As presented by [1] the compression process is composed of four basic phases.

1. Precharge, preparation and placement: A stack composed between three to six layers of molding material (about 7 mm of thickness) is placed in a mold at temperature of 150 o C. This amount of material is weighted before placed in the mold. The position, geometry and composition influence the part quality, affecting the fiber orientation evolution and the appearances of defects during molding. Normally, this precharge has a specific geometry, covers a specific volume of the mold and is denoted as preform, designed to minimize the lost of quality during the process.

2. Mold closure: After the preform is being placed inside the mold, the upper mold part moves down to touch the surface of the preformed piece, this step from the deposition to the contact of the punch against the preform last about 10 s. From this point, the compression begins with a speed in the range of 5-10 mm/s. The preform fills and releases the air entrapped through the shear edges of the mold.

The force required to move the punch during deformation reached a machine limit, decreasing the velocity of the punch and making the compression slower and slower. Finally the mold cavity if filled and the materials is retain between the punch and the mold.

3. Curing: While keeping the mold closed, the molding pressure is maintained for a determinate period of time. This step is performed while the inhibitors of the chemical paste are consumed and the chemical reaction consolidates the piece. The compressed material is then considered as rigid. Curing time depends on the resin mixture formulation, sheet thickness and mold temperature.

Figure 3: View of potential parts to be replaced by fiber reinforced composites in order to gain up to 110Kg of weight saving. In white the parts conceding the exterior bodywork, the gray stands for structural parts and the black for semi-structural parts. (Image courtesy of Plastic Omnium)

4. Part release: After consolidation, the part is removed from the mold walls by means of air ejectors pins to ensure the best quality after molding. The part cools down outside the mold by natural convection and the geometry variations of the part due to thermal contraction reaches its equilibrium. The mold cavity is prepared and sets for the next molding cycle.

The process cycle lasts from one to three minutes depending on the part thickness. For the automotive industry the challenge is to produce one part in one minute. The relative residence time of each stage can be visualized in Fig. 6. As observed, the curing time is the longest stage of the process. For the material, the temperature ranges from 130 o C to 160 o C under a molding pressure of the order of 10 MPa as depicted in Figure 5.

Some of the compression molding advantages can be addressed [1]: High reproducibility of the molded part, wide variety of mechanical properties by controlling the fiber content on the part, interior and exterior surface already finished, it allows as well the fabrication of complex geometries having ribs, curvatures, holes, inserts, etc. Moreover, many components can be consolidated in a single part without secondary assembly steps. SMC parts (thermosets) exhibit better dimensional stability than thermoplastic composites. Thermal expansion coefficient of the SMC can be set to be compatibles to steel or aluminum ones.

Figure 4: The compression Molding Process [1] Figure 5: Temperature and pression during a complete cycle [1] Compression molding pressure and temperature For the automotive industry, the first application of glass fiber reinforced polymer is found in the front panel of the GM Corvette developed in 1953 [2]. Nowadays, the most used form for compression molding is found in SMC materials. Their first applications are found in electrical and industrial products such as electrical fixtures, control boxes, light fittings tool boxes and machine guards [1]. But it was until the early 1970s that the production of exterior body parts leans towards SMCs as molding component. At that time starting by producing grille opening panels and hoods. It was also in the mid-1970s that high strength SMCs were then introduced. Applied to produce semi-structural automotive components, SMCs are found in bumper beams, road wheels, cross-members, tailgates and exterior body panels, namely, More than 1 million tailgates produced only in 2013, 5 million until end of 2014 produced by the enterprise Plastic Omnium for the Range Rover Evoque. This tailgate (See Figure 8 ) incorporates a composite design seen first on the Peugeot 508 SW. The innovative technology on this vehicle features an inner panel made of SMC with a polypropylene thermoplastic outer panel.

SMC as molding material

Sheet molding compound material used for compression molding represents the optimal ratio due to their low cost (cheapest composite per unit mass, about two times more expensive than steel per unit mass) and mold efficiency [1]. SMC is a composite material (in a sheet form) made of basically two components: resin and fiber (See Figure 9). The thermoset resin is mixed with reinforcing fibers adopting a single form denoted as compound. The manufacturing of SMC is illustrated in Figure 11.

Chopping fibers falls on the prepared resin within a thin layer covered by a plastic film made of polyethylene. Such resin contains initiators and fillers during the adopted mixture and entraps the fiber embedding the mixture. A second film, containing the prepared resin is placed on top, thus covering the fibers content. During the process air is entrapped between the fibers and the resin. This sandwich passes by compaction rolls wetting the fibers with the constituents. In Figure 12 is seen the reinforcements (fibers) to the paste resin mixture. The width of the SMC sheet is limited to the SMC manufacturing machine ranging from 0.61 to 1.52 m, whereas the thickness is commonly 3 mm. The resin is slightly cured to an intermediate stage guaranteeing the handling of the SMC plies for rool-up and shipment.

Due to the volatile nature of styrene, SMC is refrigerated since its shelf-life depends on the storage temperature. Extra additives such as Fillers (calcium carbonate (CaCO 3 )) are added in SMCs to reduce the materials cost and minimize the volumetric shrinkage of the resin.

A standard SMC weight content ranges from 30 to 50% of fiber (25 to 75 mm long, frequently E-glass fiber 25.4 mm). Approximately 25% of resin and 25∼45% filler (calcium carbonate, alumina or clay). For structural applications the fiber content ascend from 50 to 70% in weight corresponding to 40 to 55% in volume. In the latter, SMC materials behaves as compressible materials for the first stages of the compression due to the entrapped air between fiber bundles. Materials: Two industrial formulations of SMC with two fiber contents prepared by MCR -Plastic Omnium (Tournon-sur-Rhône, France) are subject of study in this work. A tomography is pictured in Figure 10. The first SMC formulation, which was denoted by the indicator (UL), is typically used to produce lightweight automotive semi-structural parts. Its behavior is mainly incompressible and corresponds to usual models found in the literature. It consists of a polyester-based paste reinforced with 29 wt% of flat glass fiber bundles with a 25-mm length, approximately 0.05 mm in height and 0.5 mm in width. The second SMC formulation is used to produce structural parts. It is denoted by the indicator (HP) and is made of a vinylester-based paste reinforced with 50 wt% of glass fiber bundles having similar dimensions than those of the (UL) formulation. Its behavior is compressible since the porosity content of the preform ranges between 20∼25%. This behavior has not been yet described and it represents a whole section of modeling of this work.

Compressibility in concentrated SMC

Modeling SMC materials implies dealing with a multiphase problem mainly governed by fiber and paste interaction. The fiber concentration plays a role in the resistance of the final part as well as on the complexity of its modeling. Typically, depending on the mechanical loading of the fabricated part, higher fiber concentration are needed. However, as pictured in Figure 3.2 for Ultra Light (UL) parts the mixture can be considered incompressible, whereas for high performance (HP) application, the entrapped air in the composite is translated into a 25% of porosity, triggering the need of compressible models.

An unified model

The compression molding of fiber reinforced composites undergoing a chemical reaction represents an example of an unsteady, non-isothermal three dimensional flow of compressible, viscoplastic fluids. One of the characteristics of the process is the coupling of the heating and the flow into the mold cavity. Whereas the mechanics gradients are mainly in the plane of the structure the thermal variation are more important within the thickness of the sample. In order to get a full understanding of the process a 3D modeling Figure 11: SMC manufacturing process [1] is required.

During the filling, when the stack of SMC flow into the cavity, we can consider the flow to be incompressible, however, the fiber network might contain entrapped air inducing a compressible stage during the compression. While the air is expelled from the composite, the thickness of the material is reduced without flowing. When the porosity are lower then 3% the material start to flow. The contact between the mold and sample creates a mechanical condition defining the nature of the compression. If slip condition is presented, the flow is purely extensional whereas at the presence of any friction the shearing becomes predominant. The process induces both normal and shear stresses. The fiber network induces an orthotropic behavior in the planar coordinates. As well the flow shearing attempts to re-orient the fiber towards the flow direction increasing the local viscosity of the material. During this stage, the thermal effects are only presented in the surface of the sample, since the diffusion in the thickness is slow comparing to the filling stage. At the end of the filling, there is an increase of the pressure and the thermal cycle finishes the process.

During the holding stage the thermal effect induces dilatation of the SMC contained by the mold and punch cavity. Such dilation is normally controlled from the conception process, however it induces thermal stresses. Whereas the temperature increasing reduces the composite viscosity the pressure reduces. However, once the ignition time for the kinetics triggers the reaction, there is an increment of the pressure since polymerization (solidification) of the part begins. Polymerization contracts the material and a strong shrinkage is presented reducing the dimensional geometry of the final part.

When the pressure of the piece seems to be stable, the molding part is considered and assumed to be fully solidified. The rigid part is ejected from the mold. The process simulation and objective of this thesis can be spitted into two main task: -A thermal model of the material includes the kinetics representative of SMC materials, its coupling with thermal properties and chemical reaction and its implication on the mechanical behavior.

Previous works

The main aim of numerical simulation for compression molding process goes further than analyzing the processing stage, but predicting the end-use condition of the molded part. In the work of [3], a simulation of a flatten sample assuming an isotropic mechanical model was used to estimate the filling process. Complemented by the thermal simulation, such approach has presented the first advantages of the simulation of the compression molding of SMC materials. Even-though such example remain a 2D approach it fitted experimental data collected by the authors. Other authors has mainly focused on only one stage of the process. Either the filling or the curing. The mechanical modeling has been studied with the work of Barone et al [4] and later improved with the study of Kotsikos et al [5]. The work of Le Corre [6] concluding that SMC materials behaves anisotropically, gave rise to an orthotropic model developed by Dumont et al in [7]. Regarding thermal approaches, first studies were accomplished by Lee in [8] and by Maazouza [9] using a model based on the radical polymerization mechanisms. The fully 3D modeling of SMC materials however has not been yet tackle because of numerical limitations. Being able to have reasonable computational time remains an open challenged. Techniques of mesh adaptation to reduce computational cost is trending nowadays being that the main advantages on which our work is based with respect to other studies of SMC.

So far, to our knowledge, a fully 3D thermo-mechanical model for dealing with SMC materials in compression molding are not yet proposed. Mostly they take into account

Context, objectives and outline

Plastic Omnium Auto exterior has interest in having a numerical tool to assist the large engineering campaigns, when deciding the designing of an auto part. On that purpose, this general project has been created. This project belongs to an assembly of three projects to study the thermo-chemomechanical behavior exhibited by Sheet molding Compounds (SMC). One of the project is in charge of the characterization of the specified SMC materials provided by Plastic Omnium. This project will formulate the rheology equations and formulation of the material. Here denoted as the 3SR laboratory. The second project is in charge of the thermo-chemical characterization of the SMC as well as the coupling with the mechanical part. They are in charge of the changes exhibited by the material when the reaction takes place. We denoted this team as LTN laboratory. Last but not least, the third project is in charge of all the numerical developments required to perform numerical simulations using the software Rem3D. Before this work, REM3D was not counting with a Compression module to perform compression molding simulations. Additionally the tools available to perform such simulations where mainly:

• a mechanical solver, determining the velocity and the pressure for isotropic materials; • a temperature solver to determine the temperature in the composite.

In general, the material could only be considered as incompressible, thermically and mechanically isotropic as well as no-coupling of the fiber orientation into the mechanical solver. Thermo-kinetical coupling to the mechanical solver due to dilatation or contraction was approximated. Validation of the models were not performed to fiber reinforced composites for molding compression. For that reason, knowing the geometry and the compression parameters, the purpose of our work focuses on:

General objectives

To enrich the software Rem3D to study molding compression processes, dedicated to Sheet Molding Compounds (SMC) for ultralight and high performance SMC:

• position of the flow front at each instant of the compression stage.

• thermodynamical characterization of the materials at each instant of the molding process determining its pressure, velocity, stress tensor, density, temperature, conversion degree.

• prediction of the mechanical force need to deform the piece.

• detection of starting time of the reaction under industrial parts.

• shrinkage and expansion during heating/cooling inside the mold.

Specific objectives

In order to achieve and fulfill the general objective, the scope of the work has been divided into specific tasks: Firstly, material compressibility must be taken into account when dealing with high fiber concentration composites. This represents a more realistic study of the geometry variation of the final form. For this geometrical variations, a corresponding computation of the density needs to be performed. A suitable compressible mechanical solver should be provided in order to describe such evolution of the porosity inside the composite material. This solver has to deal with both compressible and incompressible behavior of SMC. At the same time, a smooth transition between a compressible to an incompressible mechanical solver. Thus, enabling the continuous computation of the molding process. At this point, an important results is the pressure trace on the molding part, a detailed information of what happens inside the mold.

Secondly, because of the reactive nature of the reinforced composite a proper computation of the temperature profile along the process is needed. Even more during the reaction that solidifies the part. Predicting precisely when it starts and how long it last provides an important information of the molding process.

Thirdly, the compression molding itself is characterized for having particular macro calculation such as the compression force. Then, by using this forces most of the producers of auto parts tun the compression speed in order to not exceed an upper bound. Inside our numerical framework this represents as well a non-trivial task. Imposition of the boundary condition to immersed bodies also needs to be studied to avoid numerically perturbations of the velocity profile during the simulation. Another highlighted point is the conservative properties of level set approaches. For geometries of high aspect ratio the choice of mesh size also represents a challenge since the computational cost should be kept bounded. A strategy dealing with all theses aspects is also oriented along this work. Recalling the specific objectives we get:

• Implement a monophase rheology model for SMC behaviors.

• Enrich the monophase model to take into account the porosity exhibited in fiber reinforced composites with high fiber content.

• Implement the kinetics formulation for the case of SMC.

• Propose and implement a fully thermo-chemo-mechanical model for SMC .

• set a strategy to compute the compression force within the immersion approach

• Implement the piloting of the compression process imposing the force .

• Propose a methodology to take into account the friction against the mold walls.

• Determine the fiber orientation for compression molding processes.

• automotic selection of time step and numerical parameters to ensure numerical convergence.

We attempt to deal with all these issues in a single formulation enabling the continuous computation of the molding part. This work is divided into 4 Chapters. After this introduction, we will described the governing equations of the anisotropic compressible reactive flows given by the continuum mechanics for viscous fluids. Then, we introduce our numerical framework introducing the conservative problems encountered during this work and describing the techniques used for the solution. There, a solution is provided for the problem of the mass lost due to the inherent non-conservative properties of level set methods.

Later, the mechanical and themo-kinetical model used for both SMC studied in this work is given. There, the implementation of an anisotropic unified solver is provided. A solver enabling the continuous computation during porosity closure. Additionally a solver that by acting on a structural tensor is able to take into account the fiber network orientation. The solver is also enriched with the thermal coupling adding the thermal dilatation and kinetical shrinkage. The thermal solver is presented and the coupling with the kinetics model is provided.

Summary

Fiber reinforced composites materials as SMC (Sheet Molding Compounds) are used for manufacturing automotive parts. These materials are made of two main components: a mixture of viscous paste denoted as matrix, reinforced by fiber bundles widely-spread covered by a second layer forming a sandwich. Its numerical modeling implies a multiphase material, exhibiting a strong anisotropy since the fiber network is presented in a plane. In Addition, the reinforced materials in raw state is considered viscous and once squeezed in a hot mold undergoes a chemical reaction transforming its viscous phases into a rigid body. Thus, the modeling of SMC materials implies a thermo-mechanical problem coupling thermal effects to the mechanical response and phase transformation.

A fully thermo-mechanical model is proposed taking into considerations all the observe mechanism related in the literature of SMC materials. We use the governing equations to elaborate a model considering the anisotropy of the material, its viscous dependency on the strain rate, its thermal behavior and the description of the chemical reaction transforming the viscous raw state into a consolidated part. This model coordinates the physics and is intended to be used to predict stress evolution along the compression process, temperature evolution during the heating of the part and finally the prediction/description of the chemical reaction.

We proceed to coupled all those mechanism of interaction using continuum mechanics. We supposed a homogeneous compressible anisotropic material and we describe its motion with the velocity and pressure profile. We make used of the heat equation coupled to the kinetical model to describe the thermo-kinetical behavior of SMC.

Normally fully three dimensional simulations (3D) required a high computational cost respect to planar assumptions (2D) or well decoupled systems. However, those approximations are not able to predict the thermo-mechanical response of this materials. The model here presented is used on a numerical platform able to discretize the computational domain smartly by using mesh adaptations techniques. The model is implemented in a Finite Element Library and is based on a Eulerian approach using level set methods.

Résumé en Français

Matériaux composites renforcés de fibres comme SMC (Sheet Moulding Compounds) sont utilisés pour la fabrication de pièces automobiles. Ces matériaux sont constitués de deux composants principaux: un mélange de pâte visqueuse dénommée matrice, renforcée par des faisceaux de fibres largement répandus recouverts d'une seconde couche formant un sandwich. Sa modélisation numérique implique un matériau multi-phase, présentant une forte anisotropie puisque le réseau de fibres est présenté dans un plan. En outre, les matériaux renforcés á l'état brut sont considérés visqueux et une fois pressés dans un moule chaud subit une réaction chimique transformant ses phases visqueuses en un corps rigide. Ainsi, la modélisation des matériaux SMC implique un problème thermomécanique de couplage des effets thermiques á la réponse mécanique et á la transformation de phase.

Un modèle entièrement thermo-mécanique est proposé en prenant en considération tout le mécanisme d'observation lié dans la littérature de matériaux SMC. Nous utilisons les équations gouvernantes pour élaborer un modèle en tenant compte de l'anisotropie du matériau, de sa dépendance visqueuse á la vitesse de déformation, de son comportement thermique et de la description de la réaction chimique transformant l'état brut visqueux en une partie consolidée. Ce modèle coordonne la physique et est destiné á être utilisé pour prédire l'évolution des contraintes le long du processus de compression, l'évolution de la température pendant le chauffage de la pièce et enfin la prédiction / description de la réaction chimique.

Nous procédons á coupler tous ces mécanismes d'interaction en utilisant la mécanique du milieux continu. Nous supposons un matériau homogène compressible anisotrope et nous décrivons son mouvement avec le profil de vitesse et de pression. Nous utilisons l'équation de la chaleur couplée au modèle cinétique pour décrire le comportement thermo-cinétique de SMC.

Normalement, les simulations à trois dimensions (3D) nécessitaient un coût de calcul élevé en fonction des hypothèses planaires (2D) ou des systèmes bien découplés. Cependant, ces approximations ne sont pas capables de prédire la réponse thermo-mécanique de ces matériaux. Le modèle présenté ici est utilisé sur une plate-forme numérique capable de discrétiser le domaine de calcul intelligemment en utilisant des techniques d'adaptation de maillage. Le modèle est implémenté dans une bibliothèque d'éléments finis et est basé sur une approche eulérienne utilisant des méthodes level-set.

Introduction

Modeling SMC materials during compression molding implies the study of non-isothermal unsteady compressible flows. From one side, the fiber network induces an anisotropic mechanism, more specifically, an orthotropic behavior normal to the fiber plane. From the other side, thermal evolution modifies the mechanical response and triggers a chemical reaction until the part is consolidated. After filling the cavity, the coupled thermomechanical mechanisms enable the thermal expansion and chemical shrinkage, potential candidates of the geometry differences between the mold cavity and the final part.

In this chapter, we establish the thermo-mechanical model solved within SMC compression molding simulations and describe the numerical methods employed. The stack of SMC plies is considered as a compressible fluid. Firstly, the hot temperature induces a fluidification (a softening), the material fills the mold cavity and then (or simultaneously) the hot temperature triggers a chemical reaction consolidating the part.

We consider that the anisotropy of the SMC is due to the fiber network and we present hereafter the governing equations describing all the observed behavior of SMC materials during processing. In the first part of this chapter, the equations issued from continuum mechanics on compressible materials are introduced, considering temperature and chemical evolution and its coupling with the mechanical response. The mass and momentum conservation lead to a system having as variables the velocity and the pressure. Neglecting inertial effects, as usual in polymer flows, we get a compressible version of the Stokes problem. The Mixed Finite Element framework is applied to solve such equations. We integrate in the formulation the thermal effect by means of a general definition of strain rate contributions leading to the thermo-mechanical coupling system.

We notice that the density variations within the material are related to porosity closure in the early compression stage. Those density variations of the material are also related to its dilatation and shrinkage coefficients. We present a splitting technique to compute the density evolution. For that, it seems natural to assign the evolution of the porosity volume fraction to a relative density and the thermal expansion and chemical shrinkage to the evolution of a dense state density. In order to couple the resolution of the density equation to the resolution of the mechanical problem, we relate a partial pressure to the trace of the visco-plastic strain rate. This variable change enables the coupling of the two equations, in which a compressible factor evolves during the porosity closure. We present then the mechanical formulation taking into account the porosity dynamics and thermo-kinetic evolution.

The energy equation is introduced supposing also a compressible behavior of the material and the kinetics of the chemical reaction. We generalized the heat equation by considering an anisotropic behavior of the thermal conductivity, the thermal dilatation, the pressure-volume work and the reaction enthalpy.

After defining the whole set of equations of our thermo-mechanical model, the numerical framework is introduced. There, we recall the governing equations by detailing its discrete resolution in our numerical software using the Finite Element Method. The stabilization of the mechanical problem by introducing a bubble space and the stabilization of the thermo-kinetical problem with the Streamline Upwind Petrov-Galerkin (SUPG) strategy are briefly explained. The mesh algorithm to adapt the discrete space is explained and is coupled to a time-adaptation strategy. This adaptive time-step bounds the temporal evolution in order to control the non-linearity of the thermo-kinetical resolution. Finally, at the end of the chapter, the full thermo-mechanical algorithm is presented. The chronological resolution order is provided while picturing the non-linear behavior of the full system.

Momentum Equation

In a general form, the mechanical equilibrium of any material is given by the resolution of the conservative equation. We can then write the equilibrium of the momentum [1] in the form:

∇ • σ = ρ dv dt -f (1.1)
In eq. (1.1), σ is the Cauchy stress tensor, ρ the density of the material, v the velocity and f the volumetric forces. For the case of viscous materials, we neglect the inertia effects by ignoring the term dv dt .

Mass conservation

Additionally, the mass conservation states that any variation on the density is subject to the total variation of the fluxes in the neighboring, here described as the divergence of the flux as:

∇ • v = - 1 ρ dρ dt . (1.2)
In eq. (1.2) the density variation can be expressed by thermodynamics state functions. One can describe it by a function ρ = f (p, e) where p is the pressure and e the internal energy. Then,

∇ • v + χ p dp dt = χ e de dt .
(1.

3)

The expression given by eq. ( 1.3) has been used in other works [2] dealing with compressible flows. Mainly with gases, where the variation of the density can be directly related to energy changes. Furthermore, the energy variation for SMC materials depends on the temperature and on the degree of conversion of the chemical reaction e = f (T, α). The density is then written as ρ = f (p, T, α), where T is the temperature of the composites, α the conversion degree of the polymerization reaction and p the pressure. We can then relate the density variation with respect to other thermo-mechanical properties to obtain:

- 1 ρ dρ dt = -1 ρ ∂ρ ∂T χt dT dt + -1 ρ ∂ρ ∂α χα dα dt + -1 ρ ∂ρ ∂p -χp dp dt (1.4) - 1 ρ dρ dt = χ t dT dt + χ α dα dt -χ p dp dt (1.5)
Then, replacing the expressions provided by eq. (1.4) and eq. (1.5) into the mass conservation, we get:

∇ • v + χ p dp dt = χ t dT dt + χ α dα dt (1.6)
For high concentration SMC a volume variation is observed, where the sample thickness is reduced while the projected surface is unchanged. In this stage, the SMC UL voids are dynamically expelled from the core towards the exterior, whereas for the SMC HP the pores are closed without translation. This porosity closure, called densification, reduces the volume of the SMC up to 25% (HP) and 3% (UL). The densification step increases the pressure of the material since it becomes more and more incompressible. Once the porosity concentration is low enough, a squeeze flow motion is observed. In the studied materials, only high concentrated SMCs evidence a compressible step.

The expression given by eq. (1.6) is a general definition of the mass conservation. However, from the characterization perspective, it results more appropriate to couple directly in the stress tensor the compressible behavior by using a term proportional to the volume variation. Then, an equivalent formulation for the mass conservation is retrieved when connecting the porosity closure to the density. Indeed, the SMC materials or well the densification for porosity closure is a not reversible process. The dependence of compressible behavior proportional to dp/dt implies that the decrease in pressure translates a recovery of the porosity level. For instance, by increasing the pressure, the volume of the material decreases, whereas by decreasing the pressure the negative sign of the factor dp/dt will induce an increment of volume. However, for densification problems, when the material is released and the non-longer forces acts on it, the volume remains unchangeable. This, motivates the reasoning on a compressible factor proportional to the pressure instead of its temporal derivative. As discussed later on this chapter, the proposed model for compressibility of SMC due to porosity evolution is more suitable to be written in the form:

∇ • v + βp = χ t dT dt + χ α dα dt , ( 1.7)
the β factor takes into account the amount of porosity remaining in the material and is associated to the volumetric deformation ε v as explained later. During this compaction the density of the material is modified. By recalling the physical aspects of our thermomechanical model, we remark that the density is modified during the porosity closure and by the thermal evolution. At last, the curing of the part also implies volume variations. The compression step lasts between 3 to 5 seconds while the thermo-kinetical cycle around 1 to 2 minutes. Mostly, the variation of the density corresponding to pore closure acts at the beginning of the process. Such phenomena may be somehow decoupled from the other sources of compressibility such as thermal and kinetical evolution. For that reason, the density is expressed as the multiplication of two contributions. The relative density, ρ r , associated to the viscoplastic strain rate and the dense state density, ρ d , associated to the volume variations in the dense material due to thermal expansion and chemical shrinkage. The density is then written as:

ρ = ρ r ρ dense (T, α) (1.8)
The use of the velocity field obtained from the resolution fo the mechanical problem to obtain the evolution of the relative density has been used in [3]. A similar approach is found in the definition of the Green equivalent strain rate studied as well by Abouaf [4], [5]. Since the nature of our SMC densification problem relies on a mechanism similar to powder compaction, we prefer to use this approach. According to M. Bellet (private communication, July 2016) such approach can be extended to take into account the volume evolution from thermal variations. As in his case, for the study of thermomechanical modeling of spark-plasma sintering of metal powder and using the definition in eq. (1.8), the evolution of the density is given by:

1 ρ dρ dt = 1 ρ r dρ r dt + 1 ρ d dρ d dt (1.9)
For the numerical viewpoint, both contributions of eq. (1.9) need to be computed. In the following, we show how to calculate the evolution of both relative and dense state densities.

Strain rate tensor for thermo-mechanical modeling

Firstly, let us address the strain rate deformation tensor ε as the summation of the viscoplastic strain rate tensor D, the thermal strain rate tensor εth and the chemical strain rate tensor εch :

ε = D + εth + εch = 1 2 ∇v + ∇ T v (1.10)
The expression given in eq. (1.10) defines the total strain rate, which is the one provided by the symmetric part of the velocity gradient. The thermal strain rate tensor εth depends on the expansion tensor χ t and on the temperature evolution. Similarly, the chemical strain rate tensor εch depends on the chemical shrinkage tensor χ α and on the reticulation process evolution. Tensors are used in order to take into account the anisotropic behavior of the thermal dilatation and chemical shrinkage on SMC. As detailed later on this work (??), the thermal and kinetical anisotrope strain rates read:

εth = χ t dT dt , (1.11) εch = χ α dα dt , (1.12)
where χ t and χ α are the tensorial representations of such anisotropic behaviors.

Describing the total contribution in the motion of the SMC as the contribution of the viscoplastic, thermal and chemical strain rates is inspired from principles of thermomechanical modeling as presented in [6] for solid deformation. Here, we extend the approach to highly viscous fluid materials. Without thermo-chemical effects, the mechanical response of viscous materials is related directly to the velocity symmetric gradient (normally, D = ε). Nevertheless, when accounting for other sources of deformation in our material, namely due to thermal and chemical evolution, the viscous strain rate tensor takes the form of eq. (1.13):

D = ε -εth -εch . (1.13)
Similarly, we can associate the trace of this tensor to the velocity divergence through eq. (1.14):

tr(D) = ∇ • v -tr( εth ) -tr( εch ).
(1.14)

One remarks that the viscoplastic strain rate tensor D becomes ε for the cases where thermal or kinetical deformations are neglected. Its trace is thus directly related to the divergence of the velocity (tr(D) = ∇ • v).

When dealing with thermo-mechanical problems, only the viscoplastic strain rate tensor D generates viscoplastic stresses. That is the reason why the dilatation or contraction changes must be suppressed as suggested by eq. (1.13).

Evolution of the relative and dense state density

By integrating the expression of the divergence of the velocity given in eq. (1.14) and the expression of the density evolution given by eq. (1.9) in the mass conservation equation of eq. (1.2), the equation ruling the density evolution is:

tr(D) + tr( εth ) + tr( εch ) = - 1 ρ r dρ r dt - 1 ρ d dρ d dt (1.15)
Furthermore, we considered that the thermal dilatation and chemical shrinkage of the paste simply result on the variation of the dense density ρ d and do not affect the relative density ρ r . Allowing us to split the relation of eq. (1.15) into the relative and dense state densities, produces two separated governing equations. The relative density is ruled by the trace of the viscoplastic strain rate, eq. (1.16):

tr(D) = - 1 ρ r dρ r dt (1.16)
and the dense state density is driven by the trace of the thermal and chemical evolution eq. (1.17):

tr( εth ) + tr( εch ) = - 1 ρ d dρ d dt (1.17)
The relative density can then be rewritten by using eq. (1.14), obtaining:

1 ρ r dρ r dt + ∇ • v -tr( εth ) -tr( εch ) = 0 (1.18)
The expression given by eq. (1.18) is used then to compute the evolution of the relative density during the thermo-mechanical process.

From the mechanical point of view, the evolution of the relative density can be also related to the variation of the porosity fraction φ p inside the material. Another way to express the relative density is by providing it as function of φ p :

ρ r = 1 -φ p (1.19)
This approach is preferred since all other mechanical properties will be then related to the porosity evolution. A phenomenological expression to update the porosity fraction φ p during the compression is provided in Chapter 3. Finally, the evolution of the density in this work is computed by:

ρ = ρ r ρ d (1.20) ρ r = 1 -φ p (ε v ) (1.21) 1 ρ d dρ d dt = tr( εth ) + tr( εch ) (1.22)
and the evolution of the porosity φ p is given by a phenomenological relation of the viscoplastic volumetric deformation ε v computed using the expression:

dε v dt = tr(D) = ∇ • v -tr( εth ) -tr( εch ) (1.23)
One remarks that either solving the relative density equation in eq. (1.18) or solving the relative density by a function of the porosity fraction eq. 1.21 (depending as well as the volumetric deformation governed by eq. (1.23) ), implies the resolution of a differential equation taking as entry the term ∇ • v -tr( εth ) -tr( εch ). However, solving the porosity fraction connects the updating of the fiber fraction φ f evolution during the densification step.

Mass conservation and mechanical problem

In mechanical modeling, the Cauchy stress tensor can be always described by its deviatoric s and volumetric p parts. In order to present the equation governing the mixed problem in terms of the velocity-pressure variables, as usual for a viscoplastic material, a particular splitting is proposed. For SMC compressible materials, we can write its stress tensor as follows:

σ = s -pI with s = dev(σ) and p = - tr(σ) 3 (1.24)
However, in order to describe the porosity evolution inside the SMC by the term tr(D) and to have a clear transition between the compressible behavior of SMC towards its dense state, we prefer to use the following splitting:

σ = ŝ -pI, (1.25)
where ŝ stands for the anisotropic behavior law doctored in [7]:

ŝ = α 0 η eq D + α 1 (M : D)M + 1 2 α 2 (D • M + M • D) (1.26)
There, coefficient α i are determined by experimental test and the structural tensor M stands for the construction of the orthotropic model knowing the normal to the fiber plane in the SMC. Then, p is related to the trace of the viscoplastic strain rate and to a non-dimensional parameter α 3 representing the evolution of the porosity inside the material. The parameter α 3 increases to infinity and tr(D) tends to zero during the porosity closure. The choose of the construction of the term α 3 tr(D) responds to characterization demands since the rheology for compressible SMC are directly related to this term. A more detailed construction is given in Chapter 3. In here, we define the partial pressure to be equal to the multiplication of this terms since it would imply a more stable numerical resolution. This partial pressure is then of the form: p = -α 3 η tr(D), (1.27) This splitting implies a system (D, p) of the following form:

σ(D, p) = ŝ(D) -pI, (1.28) tr(D) = -1 ηα 3 p. (1.29)
We replace the definition of the viscoplastic strain rate given by eq. (1.14) and we introduce the compressibility factor β.

∇ • v -tr( εth ) -tr( εth ) tr(D) = -β 1 ηα 3 p (1.30)
The final set of equations is obtained by replacing the viscoplatic strain rate given by eq. (1.13)

∇ • ŝ( ε -ε th -ε ch ) -∇p = 0 (1.31)
and introducing a compressibility factor β within the mass equation, eq. (1.30), written in its (v, p) formulation:

∇ • v + β p -tr( εth ) -tr( εch ) = 0 (1.32)
The system given by eqs. (1.31-1.32) represents the coupled resolution of the mechanical problem for SMC materials taking into account porosity closure, thermal expansion and kinetical shrinkage. Please notice, that by ignoring all volumetric variations eq. (1.32) corresponds to the incompressible (divergence free ) case solving the Stokes problem. There, if isotropic behavior (α 1 = α 2 = 0), and incompressibility tr(D) = 0, p equals the pressure p. For the general case, a relation between the apparent pressure p and the pressure p is given by:

p = p - tr(ŝ) 3 = p - α 0 η eq 3 (tr(D) + (α 1 + α 2 ) M : D) (1.33)
The advantages of this formulation is the straightforward relation with compressible and incompressible cases. Enabling us to write an unified numerical procedure to fully describe the compressible-incompressible transition. Typically, β decreases in the formulation while porosity gets closed. The formulation then becomes as the incompressible one according to eq. (1.30) .

Energy equation

The last conservative equation is the energy transfer, given by the first thermodynamic law:

dE = δQ -δW, (1.34)
which represents the state variation of the system due to the interactions with the neighboring. In eq. (1.34) E, Q and W are the extensive properties (mass dependent) and account for the internal energy U , the heat entering the system Q and the work delivered by the material W [8]. Those properties depend on intensive properties such as the temperature T , the entropy s , the internal energy e, the pressure p, the velocity v. By replacing the extensive properties to their corresponding intensive one and by taking their differential, we can write the most general energy equation also known as the heat equation for any generic material:

E = Ω ρ e dΩ, ( 1.35)
The energy E represents the energy associated to the control volume. We can show [9] that the energy variation (without assuming constant density) is ruled by eq. (1.36):

dE dt = Ω ρ de dt dΩ, (1.36)
using the definition in eq. (1.34) and the expression of the internal energy eq. (1.36) enables to write the energy equation [1]:

Ω ρ de dt dΩ = - Γ q • ndΓ + Ω ωdΩ (1.37)
After using the Green-Ostrogradski theorem, or the divergence theorem, and using the local form of the energy equation, we obtain:

ρ de dt = -∇ • q + ω (1.38)
In eq. (1.38), q denotes the heat flux inside the material and ω the internal work acting as the dissipation term, since it represents the irreversible work done by the viscous forces. In order to express the energy equation in state variables, we generally need to deploy the definitions of all these terms in variables recognized by our model. Those variables are the temperature, the pressure and the velocity [1] (a full development is given in Chapter A ), the left term can be then written:

ρ de dt = ρ dh dt -χ v t T dp dt -p∇ • v, (1.39)
where we kept the enthalpy definition to be general. Here h = h T + h α , containing the energy subject to thermal variations and the energy subject to chemical reactions. The term χ v t yields for the volumetric dilatation factor (scalar) of our material. Now, the right term of eq. (1.38) is divided in two contributions. First, the internal work, which is defined as:

ω = σ : ε = -p∇ • v + s : ε (1.40)
Then, by using Fourier's law, the heat flux is given by the diffusion of heat along the material (it is expressed here for anisotropic materials), using the tensor of conductivity k that accounts for the anisotropic thermal behavior exhibited by the fiber network. Applying Fourier's law, we relate the heat flux to the thermal gradient through:

q = k∇T (1.41)
In brief, the general expression of the energy equation for compressible reactive anisotropic materials may be written in two different forms. Using the Cauchy stress tensor as:

ρ dh dt -χ v t T dp dt -p∇ • v = ∇ (k∇T ) + σ : ε (1.42)
or using the deviatoric stress tensor as:

ρ dh dt -χ v t T dp dt = ∇ (k∇T ) + s : ε. (1.43)
In the study to be performed in this work, we will prefer the second version i.e. eq. (1.43) using the deviatoric stress tensor for the internal work. please notice that for the incompressible case both expressions remain the same. Moreover, we develop the enthalpy in its thermal and chemical contribution. The first described by the temperature T and the heat capacity at constant pressure c p by means of the term dh T dT = c p . The second depending on the reaction degree α, whose contribution is assigned by the term dh α dα = ∆H α . thus, the energy equation for compressible reactive materials states:

ρ c p dT dt -ρ∆H α dα dt -T χ v t dp dt = ∇ (k∇T ) + s : ε (1.44)

Summary of conservative equations for compressible reactive flows

The model describing the thermo-mechanical behavior of SMC materials is summarized in this section. The problem statements require the computation of the velocity v, apparent pressure p, the volumetric deformation ε v , relative density ρ r , dense density ρ d , temperature T and conversion degree or curing state α. The thermo-mechanical model contains mechanical properties depending on the fiber volume fraction φ f . The evolution of such parameter is given by a direct relation of the porosity volume fraction. In general,

φ f = f (φ p (ε v ))
, for that reason, the equations governing the thermo-mechanical model presented in this section is given by the following system:

∇ • ŝ -∇p = 0 (1.45) ∇ • v + β p -tr( εth ) -tr( εch ) = 0 (1.46) dε v dt = ∇ • v -tr( εth ) -tr( εch ) (1.47) 1 ρ r dρ r dt + ∇ • v -tr( εth ) -tr( εth ) = 0 or ρ r = 1 -φ p (1.48) 1 ρ d dρ d dt + tr( εth ) + tr( εth ) = 0 (1.49) ρ c p dT dt -ρ∆H α dα dt -T χ v t dp dt = ∇ (k∇T ) + s : ∇v (1.50) dα dt = F (α, T ) (1.51)
Those equations need to be enriched by defining the partial deviatoric stress tensor ŝ as function of the viscous strain rate tensor D and the kinetic model for the computation of dα dt as well as, dilatation and shrinkage coefficients defining the strain rate tensor εth and εch .

This full thermo-kinetic-rheological model accounts for velocity-stress relation and also enables the calculation of free dilatation (or contraction) inducing a velocity profile at zero stress thanks to the definition of our viscous strain rate tensor D. Additionally it takes into account the thermal and the kinetical contribution in all the conservative equations.

After having presented all the equations according to our thermo-mechanical model, we proceed to present the numerical framework used in this study. We introduce the weak formulation of the aforementioned system as well as its stabilization.

Discretization and numerical framework

The numerical framework of this work is given by the Finite Element Method [10], assembled with a mesh adaptation tool, in an Eulerian framework. The main variables are the velocity (v), the pressure (p), the temperature (T ) and the conversion degree (α). The problem to be solved is then composed of four equations: First, we write the weak form for the momentum equation by using the splitting given in the later section. In here, p and ŝ are kept as primitive variables. The second equation corresponds to the weak form of the volumetric part of the constitutive equations, taking into account the thermal and kinetics, as well as the description of the porosity closure by the term (β) (functional variable of the volumetric deformation). The third equation states for the energy balance having the anisotropic conductivity as entry and also the heat source coming from the reaction rate. As will be described in Chapter 3 the viscous stress tensor in its anisotropic form is given by the general form:

σ = ŝ(D) -pI, (1.52) with, ŝ(D) = 2ηD + T 4 : D + T 2 D + D T 2 (1.53)
By using the definition of the strain rate tensor1 and following the procedure described in [11], we can write the system as:

D = ε -εth -εch (1.54)
Expressions for the thermal and kinetic strain rate are provided in ??. Notice that here, any assumption regarding the nature of this tensor has been established. Furthermore, the trace of the viscoplastic strain rate tensor D can be associated to the porosity evolution of the material as a pressure proportional factor β. Details on the numerical construction of the porosity evolution is given in Chapter 3. Knowing that the viscous nature of the stress is linear in its formulation, we write a relation similar to the one found in [12], where the viscous stress is only related to the viscous strain rate:

ŝ(D) = ŝ( ε) -ŝ( εth ) -ŝ( εch ) (1.55)
Conservation equations for the compressible reactive SMC in our numerical framework, where inertial forces were neglected and planar isotropic consideration for thermal and mechanical properties were done, states as:

∇ • (2η ε + T 4 : ε + T 2 ε + ε T 2 ) -∇p = ∇ • ŝ( εth ) + ŝ( εch ) (1.56) ∇ • v + β p = tr( εth ) + tr( εth ) (1.57) ∂ρ d ∂t + v • ∇ρ d -ρ d tr( εth ) + tr( εth ) = 0 (1.58) ∂ε v ∂t + v • ∇ε v = ∇ • v -tr( εth ) -tr( εth ) (1.59) ρ c p ∂T ∂t + v • ∇T -∇ (k∇T ) -T χ v t dp dt = ρ∆H α F (α, T ) + s : ∇v (1.60) ∂α ∂t + v • ∇α = F (α, T ) (1.61)
Note that the velocity is present in the resolution of the thermal problem for convective purposes and also in the solution of the reaction degree α. Furthermore, the density variation during the thermo-kinetical evolution also needs to be described. We split in this work the coupling among them as follows:

• velocity (v) -apparent pressure (p) are coupled strongly;

• temperature (T ) and conversion degree (α) are coupled by means of a point fixed resolution;

• density (ρ) is then solved separately by an advection equation

In order to study the finite elements method used to solve the velocity/pressure problem, as well as the temperature/conversion degree one, we have written the weak formulation of the mechanical and the thermal problem in the following section.

Mixed variational formulation and discretization -mechanical problem

Let us introduce the L 2 (Ω) and H 1 (Ω), the classical Sobolev spaces (Lebesque and Hilbert). Let µ and Q be Hilbert spaces and P a Lebesque space, P = L 2 (Ω), with

L 2 (Ω) = v : Ω -→ R; Ω v 2 < ∞ (1.62)
For the mechanical problem, we state Q ⊂ P, Q dense in P , such as q P ≤ q Q . We establish the variational form of eq. (1.56): find (v, p)

∈ µ × Q such that, ∀(w, p) ∈ µ × Q, for (f, g) ∈ µ × Q                      Ω 2η ε(v) : ∇w + Ω T 4 : ε(v) : ∇w + Ω T 2 ε(v) : ∇w + Ω ε(v)T 2 : ∇w - Ω p∇ • w = - Ω f w + Ω s( εth + εch ) : ∇w - Ω p∇ • v - Ω p β p = - Ω p tr( εth ) + tr( εth ) (1.63)
being a a bounded bilinear form on µ × µ, and b a bounded bilinear form on µ × p and d a bounded bilinear form on Q × P . Associating to these forms, we solve the problem:

for (f, g) ∈ µ × Q , find (v, p) ∈ v × Q such that ∀(w, p) ∈ µ × P the system satisfies: a(v, w) + b(p, w) = f, w b(v, p) + d(p, p) = g, p (1.64)
Some conditions are required to the bilinear form in order to guarantee the well representation of the problem [13]. First, a must be coercive on µ, then b must satisfy the inf-sup condition [14] on µ × P and d need to be bounded from below in the P norm [15], [16]:

∃ α > 0 such that a(w, w) ≥ α w 2 µ , ∀w ∈ µ (1.65) ∃ β > 0 such that inf p∈Q sup w∈µ b(w, p) v µ p P ≥ β, w ∈ µ, p ∈ Q (1.66) ∃ γ > 0 such that d(p, p) ≥ -γ p 2 P , ∀p ∈ Q (1.67)
Choosing γ small enough (even negative), it can be shown [17], [2] that the solution (v, p) exists and is unique.

Let Ω h be a discrete space composed of simplex K such as:

Ω h = K∈T h (Ω) K (1.68)
In isotropic mesh cases, h indicates the approximation accuracy of the subspace. This related to the mesh spacing and the diameter of the elements by the relation:

h = max K∈T h (Ω) diam(K) (1.69)
For anisotropic meshes, the indicator of accuracy can be defined by a global error depending on local mesh sizes [18] and related to the physical variable solutions. A more general treatment for the study of the error in anisotropic meshes are subject of study in [19], [20], [21].

Then, the projection operator Π h from the continuous U onto the discrete U h space satisfies:

U → U h u → Π h u = arg min ( u -u h ) (1.70)
Having a dimensional space represented by an Eulerian mesh conformed by simplexes K, we search for the solution of the discrete problem composed of functional spaces V h and

Q h of finite dimensions such that the solution (v, p) ∈ ν × Q is close to (v h , ph ) ∈ ν h × Q h .
This can be expressed in its bilinear form in the variational system:

a(v h , w h ) + b(p h , w h ) = f, w h b(v h , p h ) + d(p h , p h ) = g, p h (1.71)
Furthermore, we assume that, in the discrete problem, a,b and d satisfy the same conditions than in the continuous one. Furthermore, it is proven in [2] the existence and unicity of solution to the discrete problem.

Let d v be the dimensions of V h , d p the dimension of the P h and let us choose {φ m } m=1,...,dv a basis of the space V h and {λ m } m=1,...,dv basis of P h . We write the approximated velocity v h and the approximated pressure p h on theses basis through:

v h = V m φ m (1.72) ph = Pm λ m (1.73)
which represent the discrete interpolation of the solution, (v, p). Assuming our test functions (w h , ph ) to be the approximate interpolation functions (φ m , λ m ), our variational problem may be written in the equivalent matrix form:

A vv B vp Bpv D pp V P = F G (1.74) where V ∈ R, V = (V 1 , .
.., V dv ) t is the velocity solution vector, and P ∈ R dp , P = P1 , ..., Pdv t is the apparent pressure solution vector (in this case d p = 1). The implicit contribution of the variational problem is given by the formulation of the viscous stress (A vv and B vp ) and the velocity-pressure contribution of mass equation ( Bpv and D pp ). As depicted in ??, the trace of the thermal dilatation tensor as well as the kinetical shrinkage is related to the velocity, temperature and the reaction degree variables, acting on the mass equation by means of the expressions:

tr( εth ) = χ v t ∂T ∂t + v • χ t ∇T (1.75) tr( εch ) = χ v α ∂α ∂t + v • χ α ∇α (1.76)
where is introduced the tensor of thermal dilatation χ t built with the thermal dilatation term in each principal direction. The tensor of chemical shrinkage, χ α , is also defined. These tensors account for the anisotropy of the thermal and chemical volume variations.

The volumetric variation coefficients, here denoted with the upper index " v ", is obtained taking the trace of the thermal and chemical dilatation/contraction tensors. The tensorial representation established the directional tendency of the material to undergo geometrical variations due to thermal or chemical variations. The trace however, stands for the total volume variations and is the term found to be inserted in the mass conservation equation. A detailed explanation of this derivation is given in ??. Additionally, by deploying the total temporal derivative in its two part (local and convective), the convective part depending on the velocity v is used as implicit term in our velocity/pressure problem. Finally, the mass equation is written:

∇ • v -v • χ t ∇T -v • χ α ∇α + βq implicit = χ v t ∂T ∂t + χ v α ∂α ∂t explicit (1.77)
Writing the equations of our thermo-mechanical modeling into the corresponding matrix contribution, we get from the following left hand side members:

                               A ∈ M dv,dv (R) , A pm = Ω h 2η ε(v h ) : ∇φ m + T 4 : ε(v h ) : ∇φ m + T 2 ε(v h ) : ∇φ m + ε(v h ) T 2 : ∇φ m + ∇p h • φ m B ∈ M dv,dp (R) , B mp = Ω h -λ m ∇ • v h + λ m v h • χ t ∇T + λ m v h • χ α ∇α D ∈ M dp,dp (R) , D mp = Ω h β q h λ m (1.78)
and on the right hand side the terms:

       F ∈ R dv , F vm = Ω h -f • φ m -σ extra : ∇φ m G ∈ R dp , G pm = Ω h -λ m χ v t ∂T ∂t -λ m χ v α ∂α ∂t (1.79)
being σ extra an extra stress contribution as suggested by eq.section 1.2 as function of the thermal and chemical strain rate tensor ( εth and εch ). This represents the coupling with the thermal and chemical evolution.

For sake of simplicity, we present the formulation in a more compact development. The assembled matrix of the thermo-mechanical problem stands:

A vv = -s(v) : ∇w dΩ (1.80) B vp = -p∇ • w dΩ (1.81) Bpv = -(q∇ • v + d • v q) dΩ (1.82) D pp = β p q dΩ (1.83)
please notice that different to the Stokes problem for incompressible flows, the matrix Bpv does not correspond to the transpose of B vp .

Stabilization of the mechanical problem

We need now to choose independent spaces (V h and P h ) that satisfy the same conditions imposed to the discrete problem [22,23]. In this work, we choose subspaces (V h , P h ) given by the usual name MINI-element P 1 + /P 1, as explained in [23]. This choice enables to have linear continuous velocity and pressure fields on Ω h . The velocity is given by a combination of a linear part and a piece-wise, designated usually as bubble function. The finite element space V h is written as

V h = v h ⊕ b h , where v h = w h ∈ C 0 (Ω h ) d : w h | K ∈ P 1 (K) d (1.84)
being P 1(K) the space of polynomials of degree inferior or equal to one. The bubble must verify the conditions given by eq.1.65. The bubble function vanishes at the boundary of K and is continuous inside the element. The bubble is also defined in K as a polynomial acting on each of the three sub-triangles in 2D and four sub-tetrahedral in 3D, named also the pyramidal version of the bubble function. The discrete space in which is associated states:

b h = b h ∈ C 0 (Ω h ) d : b h | ∂K = 0 and b h | K i = 0 ∈ P 1 (K i ) d , i = 1, ..., D (1.85)
where D is the topological dimension (3 nodes in 2D or 4 nodes in 3D) and (K i ) , i = 1, ..., D is a decomposition of K into D subsimplexes. Those sub-simplexes have a common vertex in the barycenter of K. Finite element spaces for the pressure are defined:

P h = q h ∈ C 0 (Ω h ) : q h | K ∈ P 1 (K) (1.86)
The dimensions of both sub-spaces are:

dim(V h ) = d × (N n + N e ) dim(P h ) = N n (1.87)
being N e the number of elements and N n the nodes in the mesh T h (Ω). By adding the bubble, the global system to solve is now given by:

   A vv 0 B vp 0 A bb B bp Bpv Bpb D pp       V l V b P    =    F l F b G    (1.88)
where V l ∈ R d×Nn is the nodal velocity vector and V b ∈ R d×Ne represents the barycenter velocity vector. P ∈ R Nn is the pressure vector. For numerical construction, we address here some properties of the bubble functions, the first obtained by the Gauss divergence theorem denoted as the transmissible property of the bubble and the second due to the orthogonality of the bubble function respect to C. Both detailed in [24] :

K ph ∇ • b h = - K ∇p • b h (1.89) K C : ∇b h = 0, ∀ tensor C constant in K (1.90)
In order to reduce our system, we can use a classical technique by condensing the bubble function, such that:

A bb V b + B bp P = F b → V b = A -1 bb (F b -B bp P ) (1.91)
this provides a mixed velocity-pressure formulation having as unknowns the nodal velocities and pressures. The final system is written:

A vv B vp Bpv C vb + D pp V P = F l F p (1.92)
with,

C vb = -Bpb A -1 bb B bp and F p = G -Bpb A -1 bb F b (1.93)
In what concerns the implementation of this formulation in Cimlib [25,26,27], we have integrated in the library a new Stokes Anisotrope Compressible solver as a heritage of an available solver (StokesMiniTestBulle) available in the numerical library of CEMEF. The new terms added are underlined.

• for the linear subspace v h :

                               A ∈ M dv,dv (R) , A pm = Ω h 2η ε(v) : ∇φ m + T 4 : ε(v) : ∇φ m + T 2 ε(v) : ∇φ m + ε(v)T 2 : ∇φ m -∇p • φ m B ∈ M dv,dp (R) , B mp = Ω h -λ m ∇ • v + λ m v • χ t ∇T + λ m v • χ α ∇α D ∈ M dp,dp (R) , D mp = Ω h β p λ m (1.94)
• and for the bubble b h :

                   A ∈ M dv,dv (R) , A pm = Ω h 2η ε(v b ) : ∇φ m + T 4 : ε(v b ) : ∇φ b m + T 2 ε(v b ) : ∇φ b m + ε(v b )T 2 : ∇φ b m -∇p • φ b m B ∈ M dv,dp (R) , B mp = Ω h -λ m ∇ • v b + λ m v b • χ t ∇T + λ m v b • χ α ∇α (1.95)

Computation of local matrices -momentum equation

Term T 4 : ε: To introduce the first anisotropic term due the fiber network, we first simplified its integral by noticing that:

Ω (T 4 : ε) ∇v = Ω T 4 : ∇v h : ∇w h (1.96)
In the discrete space, the expression is written:

K d ijkl T ijkl V m ∂φ l ∂x k ∂φ j ∂x i (1.97) Term T 2 ε + ε T 2 :
The second anisotropic term is also simplified by deriving:

Ω (T 2 ε + ε T 2 ) : ∇v = Ω T 2 ∇v h : ∇w h + (T 2 ∇v h ) t : ∇w h (1.98)
and its discrete contribution is written:

K d ijk T ik V m ∂φ k ∂x j ∂φ j ∂x i + K d ijk T ki V m ∂φ j ∂x k ∂φ j ∂x i (1.99)

Computation of local matrices -mass equation

Term v • d: The term to take into account the implicit part of the volumete variation states:

Ω (v • d) q = K d i V m φ i m d i λ m (1.100)
Term βp: The term taking into account the compressible part resulting:

Ω β p q = K β ph q h = K β Pm λ m λ m (1.101)
Remark: The integration space for the mass term implemented in this work are proportional to the multiplication of the test function without derivatives (term v q and term p q). The space assigned to the bilinear form for this case is of second order (multiplication of two piece-linear functions). The Gauss points used in the integration corresponds after multiplication to a second order polynomial expression.

Resolution and stabilization of the thermo-kinetical problem

The heat equation is formulated as the particular case of an unsteady convective-diffusionreaction equation. Classical Galerkin methods fail and generate numerical oscillations of the solution, when the convection term becomes important. Same oscillations problems are founded in diffusion problems during the treatment of thermal shocks. These nonphysical oscillations for the classical Galerkin methods come from the discretization of the first order spatial derivative in the convective term, overcoming others diffusion or reaction contributions. To avoid these instabilities, several techniques are described in literature, among them, the SUPG (Streamline Upwind Petrov Galerkin ) [28], the SCPG (Shock Capturing Petrov-Galerkin) [29] and the RFB (Residual Free Bubble) [30]. All theses methods attempt to stabilize the numerical oscillation by adding an extra diffusion during numerical resolution. In this work, we use the RFB for diffusive dominant problem, whereas the SUPG method for the convective dominant cases. Let us use the spaces defined in section 1.2.1 For the thermo-kinetical problem, the Galerkin variational formulation is given ,in its weak form, as follows: find u ∈ H 1 0 (Ω) and find a ∈ H 1 0 (Ω) such that:

           Ω ρ c p dT dt u - Ω ρ∆H α dα dt u - Ω T χ v t dp dt u - K ∇ (k∇T ) u = K s : ∇v u Ω dα dt a = Ω F (α, T ) a (1.102)
Instead of solving the coupled system, we use a fixed point method and we split both resolutions. For the thermal resolution, the reaction represents a source input denoted as q. So, for the thermal problem, we solve for u ∈ H 1 0 (Ω) such that

Ω ρ c p dT dt u - Ω T χ v t dp dt u - K ∇ (k∇T ) u = K s : ∇v u + K q u (1.103)
and for the conversion degree, α resolution, we look, for a ∈ H 1 0 (Ω) to find α such that:

Ω dα dt a = Ω F (α, T ) a (1.104)
The stabilization of the thermal problem is given by defining test functions u, such that:

u = u h + τ v • ∇u h (1.105)
As noticed in the previous expressions, the term τ represents a numerical variable in order to stabilize the method. Basically, by adding a diffusive term to the scheme, we guarantee a stable solution. The new term stabilizes the convection scheme, since it adds a diffusion in the direction of the velocity v. Many works established that the stabilization term τ can be related to the terms defining the convection-diffusion-reaction equation by the following relation [31]:

τ K = 4k h 2 K + 2 |v| h K + |σ | -1 (1.106)
and σ is by defined using all the terms proportional to T [32]:

σ = ρ c p ∆t + χ v t dp dt (1.107)
h K is the characteristic dimension of the triangle in the streamline direction [31], [33] computed by the form [34]:

h K = 2 k α a i a ∂N α ∂x i -1 (1.108)
The stabilization of the reaction problem is given by defining the test function a as:

a = a h + τ v • ∇a h (1.109)
In a similar ways, as has been done for the thermal problem, the stabilization term is given by taking only the convective term of the stabilization term defined in eq. (1.106), obtaining:

τ K = 2 |v| h K (1.110)
We coupled this stabilization of the numerical resolution with a strategy for adapting the time step. We define our interval ∆t i such that the gradient of the temperature T and the conversion degree α are bounded to guarantee an accurate description during the chemical reaction. The next section deals with the description of strategy for space and time adaption.

Space and time adaptation

This section explains the strategy to improve the mesh discretization with a mesh adaptation technique based on a posteriori edge-error estimation. In Figure 1.1 an example of mesh adaptation is given in the computational domain, whereas in Figure 1.2 the contour of the mesh located at the iso-value zero of the level is presented. The mesh is anisotropic and follow the level set gradients of the embedded bodies. Such mesh is built following a metric construction. The metrics concept is introduced connecting it with the topology of the space. We coupled our resolution with this mesh adaptation strategy in order to reduce the computational cost and increase the accuracy. Additionally, due to the non linear behavior of our thermal and mechanical properties, the variations per increment must be controlled. The needs of time adaptation comes from the strong sharp behavior of the kinetic reaction in the composite. The formulation of the Stokes problem is not time dependent, however, the compressibility factor depends on a volumetric deformation which is time dependent. The heat equation and conversion degree evolution are strongly related to a time discretization. For that reason, a strategy for time adaptation is adopted. The temporal variation per time step is computed taken both its local variations and the convective part. The latter is determined using the velocity and edge gradient and by coupling the temporal strategy to the metric definition. The next time step is then obtained by bounding the projected variation.

Mesh adaptation on a edge-based error prediction

Let a mesh field in the space be described by the junction of edges of elements K. One edge contains two nodes defined by the subindex i, j with the value of the function Figure 1.1: Example of mesh adaptation along the level set fields. We visualize the different domains, namely, the punch and mold in red, the SMC in green, the air in blue.

U i = u h (x i ) and coordinates X i and being any variable F ij = F j -F i . The metric is defined, in dimension d and per node i containing K(i) surrounding elements as:

M i = |K(i)| d   j∈K(i) X ij ⊗ X ij   -1
.

(1.111) Let G i be the gradient of a function U i defined directly at the node i of the mesh. Generally this gradient may be obtained by solving a minimization problem of the piecewise constant gradient on elements (for the case of P1 elements). These gradients are computed on the edges and moved as solution on the nodes through argmin

G G • X ij -U ij 2
or by using the metric through the expression:

G i = d M i • U i .
(1.112) Figure 1.2: Mesh adaptation contour for 3D embedded body. Notice the mesh size along the curvatures of the geometry.

Coupez [21] shows that each measure of the interpolation error on the edge e ij is equal to the difference of the gradient through:

e i,j = G ij • X ij (1.113)
The idea of the remeshing step is to minimize the error e ij in a simplex K, under the constraint of a fixed number of nodes (equivalent to a fixed number of edges, A). The mesh is optimized when the local error on each edge is the same. In order to get this homogeneous error, some nodes need to be moved, deleted or created. A concentration of nodes relies in regions of important gradients. Coupez [21] proposed to modify the length of the edges, obtaining the same result. Firstly, a normalized global error along all the edges e edge is defined depending on the number of edges A given as:

e edge = i =j |K(ij)| e ij A d (1.114)
where K(ij) is the patch of elements sharing the edge X ij , e edge provides an estimation of the minimal average error that the current mesh, limited by the edges A, will have once the mesh is well adapted. Consequently, comparing the e edge with the local error e ij a stretching edge factor s ij is obtained. Notice that the factor 1 d appears as dependent of the space dimension d2 :

s ij = e edge e ij 1 d
(1.115)

Then the new length distribution vector Xij is computed using:

Xij = s ij X ij (1.116) X i X j Xi Xj s ij X ij
The stretching factor s ij not only modifies the edge length but also rotates a single element since the stretching factors are not the same for all the edges of a simplex. Coupez [21] showed that the new error on the edge is e new ij = s ij 2 e ij . In addition, as a consequence of repeating the process iteratively, a mesh in which the error is homogeneous and the stretching factor s ij tends to 1 everywhere is obtained. Visually, we found finer meshes in the direction of higher gradients and coarser meshes in stable domains.

Local Remesh

Generally, the remeshing technique is adopted in the whole computational domain. In order to optimize this task and reduce the computational cost, a more suitable strategy is to adjust only a particular subset of an existent low-quality grid. The local remeshing step removes the low-quality elements modifying its topology until a certain quality is reached. The concept of quality of the element is here extended not only to a geometrical factor, i.e stretched or isotropic elements. Instead, we use the junctions of function U i (a vectorial field) and relate its error to the element quality. For our case, these function U i contains a level set function, the velocity field, the temperature and the reaction degree. By adopting this strategy, we only modify the elements below a given quality.

Lagrangian-Eulerian adaptive mesh strategy

In our particular case, the mobile body (punch) has a defined velocity profile. Its associated level set moves with a fixed velocity. After the mesh is well adapted, the level set is displaced in the next time step causing the changing of the quality of the elements in its vicinity. The punch moves with a known velocity. If the refined mesh related to this level set is also displaced in time, the good quality of the mesh is kept, and the remeshing step can be avoided. For that reason, we treat the mobile body as a Lagrangian space (see Figure 2 By mixing the mesh motion formulation, a transition zone between Lagrangian Adaptive mesh and Eulerian space is represented.

Since our computational domain remains the whole picture, a convection velocity u for the transport of variables needs to defined. The velocity to transport the fields in our formulation states:

u = v -v mesh (1.117)
The mesh velocity is the velocity associated to the mesh motion during the simulation. For nodes that does not move this velocity is zero. For nodes moving at the imposed velocity such as the punch (v mesh = v BC ). For the Lagrangian zone v mesh = v and for the Eulerian part v mesh = 0. Please notice that u = 0 in the Lagrangian mesh, meaning not transport is needed and u = v in the Eulerian, which is the general formulation of our framework.

Time adaptation based on a bounded increment gradient

General idea

The strategy is based on the principle of bounding the increment variation of the variable U . Knowing its temporal derivative, U , we look for:

∆U t = U t+1 -U t = U ∆t ≤ ∆ MAX (1.118)
We search for a time step ∆t such that the function U i does not vary more than the bounded value ∆ MAX .

Strategy in our numerical framework

For a given scenario in time t, computed using a time step ∆t -, we compute now the temporal variation G t i of the variable U , as:

G t i = U i -U i-1 ∆t + v • dM i • U i dU dt = ∂U ∂t + v • ∇U (1.119)
where we use the estimation of the spatial gradient of U by projecting the metric M , defined in eq. (1.111), in the direction of the variable U , as used in the mesh adaptation strategy. Please notice that dM i • U i = G i according to eq. (1.112) . The predictive variation is given by the projection of the gradient in the current time step:

∆ t i = G t i ∆t - (1.120)
Similarly to the space adaptation procedure, we compare the temporal variation with respect to a given bounded value ∆ temp and we determine the temporal stretching factor s t as:

s t = MAX ∆ t i ∆ temp -1 (1.121)
Finally, this stretching factor gives rise to the next time step ∆t with the form:

∆t = s t ∆t - (1.122)
Notice that the stretching factor is 1 when the discrete evolution equals the criterion (∆ t i = ∆ temp ).

Resolution and optimization

The resolution of the linear system is performed using the PETSC (Portable, Extensible Toolkit for Scientific Computation) library, through a preconditioned iterative method.

• Transport of the levelSet function + Correction for mass conservation (See Chapter 2)

Computing of the position of each phase by adopting a convective scheme for the fluid and a rigid motion for the punch and the mold.

• Mesh adaptation

min {e ij } = L 2 : e ij {v, T, α, φ, δ φ } • Mechanical problem (non-linear -fixed point )        ∇ • (2η ε + T 4 : ε + T 2 ε + ε T 2 ) -∇p = σ extra ( εth , εch ) ∇ • v + β p -tr( εth ) -tr( εth ) = 0 -------------- update η, T 2 , T 4 ⇒ (v t+1 , p t+1 )
• Evolution of the viscoplastic volumetric deformation in the material

     ∂ε v ∂t + v • ∇ε v = β p ρ r = 1 -φ p (ε v )
enabling the actualization of the porosity fraction φ p , the local fiber fraction φ f and viscosity consistency η f .

• Thermo-Kinetical problem (non-linear -fixed point )

                     ρ r ρ d c p dT dt -∇ (k∇T ) = ρ∆H α F (α, T ) + T χ T dp dt + s : ∇v dα dt = F (α, T ) 1 ρ d dρ d dt = χ v T dT dt + χ v α dα dt update κ, C p , F (α, T ) ⇒ (T t+1 , α t+1 , ρ t+1 d ) • Calculation of next time-step ∆t bound G t i = G t i {T, α, φ}

Conclusions

A full thermo-mechanical model was introduced to compute the reactive, non-isothermal 3D compressible flows, based on continuum mechanics principles. Compressibility and reaction induce new unknowns with respect to the traditional isothermal incompressible flows found in literature when dealing with SMC materials. The porosity on the material decreases during the compression of the piece and is always a non-reversible process.

From that, its modeling does not depend on the sign of pressure evolution. Thermal and chemical evolution however, does modifies the density in a reversible way. Its modeling, as the one found in literature for other materials, depends on the derivative of the temperature and conversion degree. Integrating a model accounting for all these phenomenological evolution has not been established before for the study of SMC materials.

The full thermo-mechanical coupling, integrates the general contribution of diverse source of strain rates deformations into a momentum a mass conservation equations. We prove the a full thermo-mechanical system may required the resolution of 6 equations having as variables the velocity, the pressure, the volumetric deformation, the temperature, the reaction degree and the dense state density.

In addition, this chapter discussed the integration of such anisotropic behavior of SMC in our numerical framework. There, we present the Finite Element Method approach used in this work. We presented punctually all the contributions of the anisotropic compressible Stokes solver we implemented, with respect to the standard Stokes resolution.

The next chapter focus on the description of the level set method and its strategy when modeling multy-body objects. Some drawbacks of the methodology are highlighted and the solution for the problem of mass lost on dynamics immersed objects is given. Additionally, it is shown the potential of our numerical tools to deal with complicated 3D geometries in compression molding simulations.

Summary

In this chapter, the use of level-set method to describe surface under deformation is considered. The methodology of level-set approach as a signed distance function being mostly narrow-banded is introduced. Particularly, the strategy used in this work to deal with compression molding simulation is explained detailing the advantage and requirements of the implicit description of surfaces. Mainly the punch, the mold and the pre-form are described by the implicit level set function in a whole domain (immersion framework).

The narrow-banded level set used generally to improve and simplify the convection algorithm to describe the motion of the object as well as to support a mesh adaptation strategy is revised. Here, it is used the linear hyperbolic tangent band and the convection is supported on a local booster (reinitialization of normal velocity to move the iso-values for keeping metrics properties).

Consequently, first a study on consevative properties of this method is carried out concluding that during compression the gradient of the velocity within the thickness can induces volume/mass lost of the preform due to the discrete scheme of the level set advection step. This volume lost is a known drawback of the level set method. For the case of a simple disk compression, during this phd was recorded a volume lost of 6%, whereas for industrial applications the lost was recorded up to 30%. A study of the impact of the mesh adaptation and the time step discretization is performed concluding that the mesh adaptation algorithm causes higher lost than using a fix isotropic mesh and that smaller time step for the convection also translates higher volume lost. These results can be understood since the mesh adaptation algorithm has a step of transport of variables and the updated configuration might not be conservative. The time step discretization show that if the time step is smaller, the lost per iteration is smaller. However, since to complete the same simulation more iterations are used, the lost at the end of simulation are higher. Knowing that this topic differs form the main objective of this work, a engineering solutions was proposed. Based on a Newton-Raphson algorithm a corrective procedure is used and detailed to recover the volume/mass lost per iteration along the computation fo industrial parts. An exmaple using industrial geometry proves the decrement from 30% down to 0.1%.

The last section of the chapter introduces a immersed technique to compute the compression force during deformation of the part. A problem for immersion domain

where the standard formula force equal stress times surface can be not obvious since the surface are not explicit. The strategy is based on the viscous power dissipation principle and proved to reduce oscillations and to provide with the same solution that standards highly mesh dependent solutions.

Résumé en Français

Dans ce chapitre, on considère l'utilisation de la méthode de level-set pour décrire la surface sous déformation. La méthodologie de l'approche de niveau comme une fonction de distance signée étant principalement à bande étroite est introduite. En particulier, la stratégie utilisée dans ce travail pour traiter la simulation de moulage par compression explique l'avantage et les exigences de la description implicite des surfaces. Principalement, le poinçon, le moule et la préforme sont décrits par la fonction implicite de réglage de niveau dans un domaine entier (cadre d'immersion). Le jeu de niveaux à bande étroite utilisé généralement pour améliorer et simplifier l'algorithme de convection qui décrire le mouvement de l'objet ainsi que pour soutenir une stratégie d'adaptation de maillage est révisé. Ici, on utilise la bande tangente hyperbolique linéaire et la convection est supportée sur un booster local (réinitialisation de la vitesse normale pour déplacer les iso-valeurs en gardant les propriétés des distances).

Par conséquent, d'abord, une étude sur les propriétés conservatives de cette méthode est réalisée en concluant que pendant la compression, le gradient de la vitesse dans lépaisseur peut induire une perte de volume / masse de la préforme due au schéma discret de l'étape d'advection dans la méthode level-set. Ce volume perdu est un inconvénient connu de la méthode level-set. Pour le cas d'une compression de disque simple, pendant ce fichier a enregistré un volume perdu de 6 %, alors que pour les applications industrielles, la perte a été enregistrée jusqu'Ãă 30 %. Une étude de l'impact de l'adaptation du maillage et de la discrétisation du pas de temps est réalisée en concluant que l'algorithme d'adaptation du maillage entraîne une perte plus élevée que l'utilisation d'un maillage fixe isotrope et que le temps de passage plus petit pour la convection traduit également une perte de volume plus élevée. Ces résultats peuvent être compris puisque l'algorithme d'adaptation au maillage comporte une étape de transport des variables et la configuration mise à jour pourrait ne pas être conservatrice. La discrétisation temporel montre que si l'intervalle de temps est plus petit, la perte par itération est plus petite. Cependant, puisque pour compléter la même simulation, plus d'itérations sont utilisées, les pertes à la fin de la simulation sont plus élevées. Sachant que ce sujet diffère de l'objectif principal de ce travail, des solutions d'ingénierie ont été proposées. Basé sur un algorithme de Newton-Raphson, une procédure corrective est utilisée et détaillée pour récupérer le volume / la masse perdue par itération le long du calcul des pièces industrielles. Un example utilisant la géométrie industrielle prouve le décrément de 30 % à 0.1 %.

La dernière section du chapitre introduit une technique immergée pour calculer la force de compression lors de la déformation de la pièce. Un problème pour le domaine de l'immersion où la formule standard forçe égale constrainte fois surface n'est pas Ãľvidente car la surface n'est pas explicite. La stratégie est basée sur le principe de dissipation de la puissance visqueuse et s'est avérée réduire les oscillations et fournir la même solution que les normes de solutions fortement dépendantes de mailles.

Introduction

The Level set approach is a conceptual framework for the analysis of surfaces and shapes initiated by Osher [1]. Such modeling eases the numerical computations involving curves or surfaces and can be performed on Cartesian grids without needs of space transformation. Within the Eulerian framework, the level set method enables, by solving differential equations, the tracking of topology changes, when modeling time-varying objects. More often, it provides an implicit function which allows the description of the interface, to which we can link other functionalities. For instance, a direct and accurate computation of the surface normal vector and its curvature, in particular, during ongoing simulations.

Punch (phase 2)

Die (phase 3) Composite (phase 1) Air (phase 4)

Rigid Body

Rigid Body

Body Immersion

Figure 2.1: Schematic illustration of the different phases in compression molding that may be consider for the modeling. All phases "immersed" are in a single computation domain. We distinguish the composite (H1), the punch (H2), the die (H3) and the air (H4) by means of implicit functions.

In fact Level set is a framework to describe boundaries. Using Level set in immersion methods extends such framework to a volumic representation of objects. In order to represent volumes using this implicit function, one considers a single computation domain and each object (body) is "immersed". In Figure 2.1 the different necessary objects for compression molding simulation are the composite, the die and the punch. All immersed in the domain. As a result, four phases appear: one per whole body and the fourth one denoted as the air.

In this work, we solve the governing equations for fluid dynamics applied to compression molding simulation in the whole computational domain. All the bodies immersed in the computational domain are analyzed by adopting adequate behavior laws and constitutive equations. Two objects are considered as rigid tools, whereas the composite material and the air are modeled as explained in Chapter 1. However, the motion of the punch and of the die bodies is pre-defined since the velocity is given as input parameter. In our formulation, those bodies (assumed rigid) are modeled by increasing its viscosity in order to have strain rates much lower than the one of SMC and air. The motion of the interfaces are ruled by: In order to improve the discrete description of each phase, the topology (mesh) is refined in the vicinity of the respective interfaces after immersion. We solve the flow physics in the whole domain and the mesh is adapted (Figure 2.2) along a thickness E in the neighborhood of the interfaces. Figure 2.2a shows the anisotropic mesh sizes, highly dense near the interface. There, a thick region containing the minimum mesh size is observed. Usually, the interface between two bodies is located in the center of the aforementioned region or thickness, as illustrated in Figure 2.2b.

3 bodies + air =                  1 Mobile -→ velocity v imposed 2 Fluids (SMC + Air) -→ velocity v obtained by ∇ • σ = 0 1 Static -→ velocity v = 0
A first option for the implicit function computation is to use the level set function, φ, which is the signed distance function and an independent one should be defined per object. However, one level set function enables the differentiation of two bodies. The relation Phases = Level sets + 1 allows to have one level set less than bodies in our computation domain. Particularly for representing 4 phases we required just of three Level sets. When assigning one level set function to a particular phase, this function is defined positive inside the phase, negative at the exterior and zero in the interface as defined eq. (2.1) :

φ =              +d(x, Γ) if inside the body 0 if in Γ -d(x, Γ) if outside the body (2.1)
Once assigned a Level set per body, we have a fully implicit description where all the interfaces are defined by its corresponding level set function. Knowing the velocity of each object, we reconstruct a distance to an interface by solving a PDE. More specifically, to displace an interface, we solve a convection equation ,eq. (2.2), using the velocity resulting from the resolution of the conservation equations:

∂φ ∂t + v • ∇φ = 0, (2.2)
This resolution will need to be coupled to a re-distancing algorithm to re-establish the metric properties of φ, either by solving the eikonal equation [2] or by adopting an unified formulation [3]. Further details will be given in the next section. We mainly need to adapt the mesh because we need a good continuous description of the properties that are normally discontinuous. Once the level set is moved, the spatial discretization must follow the interface to preserve a smooth discrete description of the continuous φ and of the discontinuous properties.

In order to couple this level set to the mesh adaptation algorithm, the level set function φ is smoothly bounded along a characteristic thickness E (here denoted as narrow Level set). In Figure 2.3 the narrow level set function is introduced without applying a smooth transition. One of the advantages of the narrow level set is the limitation of the convection equation to a space near the interfaces. Additionally, connecting the level set gradient to the mesh provides an intrinsic coupling of the discretization with the kinematics as will be depicted later on this chapter. However, the truncation needs to be smooth since the trucanted level set illustrated in Figure 2.3 contains to point where the gradient is infinite. This causes numerical problems when solving eq. (2.2). After the truncation of the level set φ, the identification of the two materials bounding this interface is performed. The inherent properties of these materials such as the density ρ, the viscosity η, etc must be distinguished in the finite space E. We define then the Heaviside function H per phase to equal 1 in one phase and 0 in the other phases. This function is equal to 0.5 at the interface and normally is discontinuous. We need to make it continuous since the function is applied on discrete spaces and discrete functions can causes numerical instabilities during resolution. We attempt to describe its discontinuity by assigning a continous evolution of the function and its gradient along the characteristic thickness e. Normally, the thickness of this transition is given by a space region e such that e < E, in which any property changes smoothly between one body to another. It was shown in [4,5,6] the evolution in the convergence of the method being strongly dependent on this thickness e. Summing up, the thickness E and e define two important characteristic lengths of the approach:

• E -→ defines the smoothing length for the phase function

• e -→ defines the numerical precision through the smoothing of the physical properties.

Coarse Mesh Fine Mesh

Coarse Mesh

φ after convection φ before convection Air ρ 2 , η 2 Composite ρ 1 , η 1 φ φ t+∆t e E Figure 2
.4: Relative position of the interface within the fine mesh zone, before and after the convection scheme. The fine mesh zone of thickness E contains not only the interface after convection, but also the thick space e for numerical precision.

These two thickness are independent and the only condition established so far is that e < E. Figure 2.4, schematically shows the level set convection within the mesh thickness E. The level set φ is displaced to the position φ t+∆t , but the thickness e after convection remains inside the thickness E. If the interface leaves the finner region, the PDE using for the level set convection may cause instabilities in the solution. Such instabilities will determine unaccurate physical properties of the density, viscosity, etc that will cause further problems during the mechanical resolution. By guaranteeing that after convection the thickness e remains inside the fine mesh zone, the good description of the physical parameters before and after convection are guaranteed. Since the thickness e always remains inside the finner volume delimited by E. Eventually, the re-meshing tool is called and the interface is set again inf the middle of thick region E. This represents one of the advantages implicitly contained by using two independent thickness e, E instead of a only one characteristic length.

All the previous points, addressed in this introduction, embed issues when using the level set approach as phase description support. Those issues were translated into a not proper description of the convective scheme within the level set advection step, loss of volume and penetration of the tools into the bodies, due to rigid motion of the punch and die. They were treated in this work and are detailed in the following.

A first issue to be treated concerns the choice of a smooth function φ and the suitable strategy to choose the thickness e. Former works were performed [7] taking a hyperbolic tangent level set supported on a thickness E. Then, in order to compute the Heaviside function, a non-linear expression to the hyperbolic level set function is calculated. The question remains for the Heaviside function H that is supported on a non-linear function, that needs also to be defined.

Another problematic, associated to the discrete resolution of the convective scheme in the level set method, is addressed hereafter. The non-conservative form of the discrete function eq. (2.2) does not guarantee the conservation of the volume for incompressible flows or mass for more complex problems. Many authors have suggested improvements to the level set method by coupling the convection equation to a correction step, recovering the conservative and metrics properties. However, even if those strategies have improved the convection step, the methods has not been well established for complex 3D industrial cases. Further studies should be conducted to justify their robustness.

In our case, the velocity profile of the mold and of the punch are known in our numerical set up, with a defined velocity profile. In immersed methods only the imposition of Dirichlet conditions in the boundaries of the domain guarantee a stable solution. Imposing to the whole body a rigid motion may introduced instabilities int he solution. The modeling of stable rigid motions in embedded geometries is of interest when dealing with compression molding simulations and is discussed later in this chapter.

The steps described in the following of this chapter are hereafter introduced. Firstly the level set approach is formalized and the step assigned to conserve its metric properties is analyzed. Examples of different smoothed level set and Heaviside functions are given, focusing on their connection between the thicknesses E and e. We will see that the introduction of a linear hyperbolic tangent level set as smooth function enables the exact computation of a Heaviside function and of a Dirac function within the thickness e. A study on local re-initialization is conducted, focused on the influence of the re-initialization velocity, when solving the convection-reinitialization step. There, we propose a strategy to improve the computation of the reinitialization of the level set based on the normal velocity gradient.

Secondly, as a straightforward application of the linear hyperbolic tangent level-set, we have developed a conservative strategy to correct the volume loss during the convection of the level set in discrete schemes. Our methodology introduced as handy in a Newton-Raphson algorithm, requires the computation of a Dirac function as the derivative of our assigned Heaviside function. The Dirac function stands within the thickness e, putting in evidence the advantages of the linear hyperbolic tangent level set. We enriched the method to account for mass conservation, when studying compressible flows. We present numerical cases, where the conservation algorithm proof proficiency and we comment on the limitations of the approach.

In the last section, we address the problem of computing the compression force on a "diffused" or implicit boundary" interface approach. The force is the mechanical variable describing basically all the rheological response of the material under deformation. The force is used to validate the numerical approach when comparing it with the force exerted by the machine in experimental tests. Instead of computing the force by a surface integral as used in standard methods, we use a virtual work principle based on the viscous power. Thus, we associate deformation energy to the compression force getting stable and proper values.

Background in Level set methods

Level set methods were introduced to represent implicitly interfaces and free surfaces Γ by defining a scalar field φ(x), eq. ( 2.3), such that:

Γ = {x/φ(x) = 0} .
(

where φ(x) was defined previously. For multi domain body immersed computations, a level set is defined to each body. By making use of this definition, a smooth Heaviside function eq. (2.4), will be built over a characteristic thickness of 2e, as defined by:

H e (φ) =                1 if φ > e 1 2 1 + φ e + 1 π sin( πφ e ) if -e < φ < e 0 if φ < e (2.4)
This enables a continuous definition of the properties transition along the interfaces. This functions will help during the resolution of the mechanical problem to identify the different phases in the domain, assigning as inlet the properties of each material independently. For instance, given a property η to each material, we can define a single space-dependent for the numerical resolution by doing:

η = H SM C η SM C + H mold η mold + H punch η punch + H air η air , ( 2.5)
The last phase (air) is defined as the complementary of the other phases (H air = 1 -

H SM C -H mold -H punch ).
After resolution of the system given in section 1.1.5, the immersed bodies can be identified by H(φ) and the boundary ∂Ω n of each sub-set by the iso-value zero φ n = 0. Moreover, the interface may evolve since a velocity field affects not only the sub-set but also the interfaces. For tracking the evolution of the interface, each sub-set φ is convected.

Convection and Re-initialization

Once our level sets functions defined; for t = 0 as φ 0 , their motion can be described by the advection relation, using the transport velocity field u:

∂φ ∂t + u • ∇φ = 0, (2.6
)

φ(t = 0, x) =φ 0 , ( 2.7)
which generalizes eq. ( 2.3). Convecting the field φ only guarantees its iso-value zero to be correct. For whatever velocity field, solving eq. (2.6) normally disturbs the metric properties of the level set function since the transport velocity u in the direction of the level set gradient causes the iso-contours φ to get closer or farther from the interface. This implies that for any iso-value φ = 0 the convection step will enhance loss of the distance property. For this reason, in order to recover the distance properties, a re-initialization step is performed, by solving the eikonal equation:

∇φ = 1.
(

Some authors [2] solve directly eq. (2.8) when redistancing is necessary, others [8, 9, 10] couple the convection-reinitialization steps and solve a generalized form of the advection equation given by eq. (2.9):

∂φ ∂t + (u + λU) • ∇φ = λ s(φ), (2.9)
where λ is a numerical coefficient, so-called booster, depending on mesh-size and timestep, U is a re-initialization velocity defined by U = s(φ) ∇φ |∇φ| and s(φ) is the sign function (by convention equals 1 inside the object, 0 at his boundary, -1 at the exterior).

The selection of both λ and the time step are based on the Courant-Friendriechs-Levy condition (CFL condition) computed taking the size of the element in the direction of the velocity u.

Narrow band level-set function

In order to solve the convection equation, an accurate description of the velocity field u in the whole domain must be guaranteed. However, the usage of a level set present lots of advantages. Let us bound φ in the space such that only in the vicinity of ∂Ω n it varies from [-E, E].

• by limiting the value of the distance function, a boundary condition may be fixed and set to -E;

• convecting in the whole domain a level set seems inefficient, since the re-initialization step will always be done.

• we can couple the mesh adaptation to this narrow function bounding the fine mesh to the interest zone, by using its gradient to automatically generate mesh sizes.

Thus, one may use a smoothly truncated distance function, keeping the same iso-value zero and switch properly to a constant value when the distance increases. One example is the hyperbolic tangent level set function, eq. (2.10) [7] :

φ E = E tanh φ E ≈ φ
, when E tends to infinite or in the vicinity of the interface (2.10) Normally, the parameter E defines the bounded values of the level set. The space region where the function varies from [-E, E] approximately equals 3E. Another example of truncated level set, is the linear sinus eq. (2.11) [7]:

φ E =                    2E π if φ > E 2E π sin( πφ 2E ) if -E < φ < E - 2E π if φ < -E (2.11)
In this case, the mixture thickness e is independent of the value of E.

Advection of the narrow level set

If the level set function is smoothed, the re-initialization procedure applied can be better used by taking advantage is smoothing. As explained in [3], the re-initialization step using a pseudo-time to apply the iterative algorithm states eq. (2.12):

             ∂φ E ∂τ + λU.∇φ E = λs(φ E )g(φ E ) φ E (τ = 0, x) = φ t-1 (x) φ E (τ, x = ∂Ω) = -E (2.12) recalling that U = s(φ E ) ∇φ E |∇φ E |
, where the solution at any time in the boundary of Ω is set to the bounded value of the smoothed level-set. The function g(φ) is given by eq. (2.13):

g(φ E ) = |∇φ E | (2.13)
Then, the re-initialization function may now be given by eq. (2.14):

               ∂φ E ∂τ + λs(φ E )(|∇φ E | -g(φ E )) = 0 φ E (τ = 0, x) = φ n+1 (x) φ E (τ, x = ∂Ω) = -E (2.14)
where s(φ E ) is the sign function. The factor |∇φ E | -g(φ E ) represents the difference between the discrete gradient and the exact nodal gradient computed using eq. (2.15) as function of the smoothing function. For example:

g(φ E ) =                  1 for no smoothing 1 - φ E E 2 for hyperbolic tangent smoothing 1 - πφ E 2E 2 for sinus smoothing (2.15)
Once |∇φ E | equals g(φ E ) ,the re-initialization algorithm converges. Moreover, we can define a regularized sign function, s(φ E ), for a smooth transition, as follows, eq. (2.16):

s(φ E ) = φ E |φ E | + E (2.16)
where the transition is proportional to the the size of the element in the velocity direction. The thickness of the transition is proportional to the mesh size (E ∼ h).

Meanwhile, the CFL condition, eq. (2.17), in pure advection problems states:

| v| ∆t ∆x ≤ 1 (2.17)
The choice of the λ parameter is based on the condition:

λ ∼ |λU | ∼ h(U ) K ∆t (2.18)
where h(U ) is the size of element K in the direction of U .

After defining the narrow level set, using one of the previous forms, the smoothed Heaviside shape function can be obtained from eq. (2.19) using:

H E = 1 2 1 + φ E E (2.19)
This expression is equal to 1 inside the object, 0 outside and 0.5 in the interface (φ E = 0 in the interface ). However, respect to the discontinuous case this functions contains a smooth transition respect to the standard Heaviside function explained in the introduction.

Approach

In this section, we present our proposal for the smoothing of level set functions. First, we present the linear tangent hyperbolic level set used in this work and we write our Heaviside function, directly as a function of the smoothed level set. Then, we introduce the formulation for the re-initialization technique by adapting such equation to our linear tangent hyperbolic. Finally, we propose the computation of important numerical parameters, like λ, as dependent of the gradient of the velocity in the normal direction. Arising to the discussion of a local parameter instead of a constant value.

Linear hyperbolic tangent

As previously mentioned, the importance of truncating the level-set enables a faster resolution of eq. (2.9). Also, an accurate resolution is given by guaranteeing a smooth solution for the velocity in the vicinity of φ = 0. Using a hyperbolic tangent automatically smoothed allows us to write the re-initialization step, where the gradient of the level set is represented by a function g(φ) containing explicitly φ.

We seek for an level set function such that the thickness e of the mixing layer can be defined in the region in which φ E is linear, guaranteeing:

φ E φ→0 = φ, ( 2.20)
that imposes some conditions between e and E. For ensure this linearity, in this work, we apply a linear hyperbolic tangent function of the distance, eq. ( 2.21):

φ E =                  e + E tanh( φ -e E ) if φ > e ±φ if -e < φ < e -e + E tanh( φ + e E ) if φ < -e (2.

21)

The main contribution of this Level-set construction yields on the exact value of the smooth level set along the thickness of the numerical precision e (notice that φ E = φ in the region -e < φ < e). The level-set is equal to the distance function. The previous expression (sinus and tangent hyperbolic) where approximate values of the distance, while this construction ensure the exact distance values.

A smoothed Heaviside function is employed over a characteristic thickness of 2e using as input the truncated level-set φ E eq. (2.22) :

H(φ E ) =                1 if φ E > e 1 2 1 + φ E e + 1 π sin( πφ E e ) if -e < φ E < e 0 if φ E < -e (2.22)
We extend the descriptive functions by adding a Dirac function to the approach, thus enabling operations in the subset ∂Ω inside the volume Ω. The Dirac function associated to the Heaviside function defined in eq. (2.22) stands: In the vicinity of the interface, the latter parameter can be determined according to eq. (2.24):

δ e (φ E ) =          0 if |φ| > e 1 2e 1 + cos( πφ E e ) if -e < φ < e (2.
h d = max j,i∈K ∇φ • X ij , (2.24)
where K is a mesh element in the tight sub-layer around the interface and X ij = X j -X i is an edge of K.

For the linear hyperbolic tangent the value of the reinitialization function, g(φ), results in eq. (2.25):

g(φ ε ) =        1 - φ ε -e E 2 if |φ| ≤ e 1 if -e < φ < e (2.25)

Re-initialization velocity -A local description

The coefficient λ, mentioned in eq. (2.18), represents a numerical parameter to line spatial and temporal metrics and has been used for a long set of simulations as presented in [10,11,12]. However as its definition shows, the choice of a re-initialization velocity does not depend on a local evolution of the surface. For instance if rigid body motions λ should equal to 0, since the level set has not been deformed. The booster λ should reflect the amount of reinitialization required depending on how much the level set is locally deformed.

Three sources of level set deformation are here described. Firstly, one related to purely flow kinematics. Secondly, another source due to numerical diffusion associated to the numerical discretization of the advection step given by eq. (2.6). Thirdly, the discretization of sharped edges, which tends to get smooth due to the resolution of eq. (2.14), being difficult to overcome without a mesh adaptation technique.

based on flow kinematics

Let us define Ω s a subset of Ω ∈ R d and let ∂Ω s describes the boundary of Ω s such that ∂Ω s ∈ Ω. Let us now divide the boundary ∂Ω s into two sections such that ∂Ω s = ∂Ω 1 ∪ ∂Ω 2 , in which we apply different boundary conditions. Let us define a level set φ s through eq. (2.26):

∂Ω s = {∀φ/φ s = 0} . (2.26)
This implies that our level set function is defined along the domain ∂Ω s . Let us define, for each sub-domain ∂Ω 1 and ∂Ω 2 , a thick vicinity δΩ 1 and δΩ 2 , compact in a thickness [-e, e], such that the velocity condition is given by eq. (2.27)-eq. (2.28)

δΩ 2 ∂Ω 1 Ω s u = v ∈ ∂Ω 1 , (2.27) u • ∇φ = 0 ∈ δΩ 2 .
(2.28)

When solving eq. (2.14), the solution in the subset Ω s needs the convection of φ = 0 along ∂Ω s and also the re-initialization of its vicinity δΩ s . However, solving in ∂Ω 2 does not imply convection of the φ = 0 and it implies also that no re-initialization is needed. For that reason, the use of a locally λ should be instead implemented. Also, in the cases where the normal velocity on φ is more important (and for that the re-initialization is more important) a different λ should be used. This indicates that assigning a local λ proportional of the local normal velocity (u • ∇ φ) seems a better strategy as was also suggested in [9]. In our work, we propose a λ proportional to the gradient of such velocity, through eq. (2.29):

λ ∼ ∇ u • ∇ φ ∇ φ ∆t (2.29)
being ∇ φ the unit-normalized gradient representing the normal to the interface φ. This gradient is narrow banded, it is zero everywhere except in the vicinity of the interface. This concept points out the need of a λ value proportional to the velocity gradient in the normal direction. This expression of λ is proportional to the time step, namely, if a gradient normal to the velocity exist, the iso-values displacement are proportional to the time step and thus its re-initialization velocity. To better illustrate this point, let us define a circle with an initial diameter of 0.1 m immersed in a box of 0.3m × 0.2m. This circle is submitted to a velocity profile of the form:

Horizontal extension of a circle

v =      2 (x -x o ) 3 2 (y -y o ) 2 (2.30)
This velocity profile in eq. (2.30) is not divergence free. According to the proposed local booster, when convecting the level set function some areas required more and some less speed to reinitialize. In this study, the velocity profile is at the origin of an horizontal stretching of the circle, while shrinking it vertically. The local λ is higher in the horizontal direction while lower in some particular sections of the circle. In Figure 2.6, we compare the solutions for the convection of the circle under the defined velocity profile. The level set convection using a constant λ loses its metric properties, visualized by the deformed iso-values. Notice, that for the case with local λ, the iso-values remains correct.

The simulations were performed setting a fixed time step, ∆t = 0.01s. As indicated by eq. (2.29), the gradient of the velocity in the normal direction is too high a decrease of the time step should be then applied. In such manner, a combination of the CFL condition and a bounded local λ can be studied to improve re-initialization strategies of level set methods. However, a deeper study in such field is out of scope from the objective of this work. Further understanding is proposed as a perspective, at the end of this document.

Another important point when using level set methods is the evaluation of their conservative properties. The convective scheme of the approach is not conservative, providing possible mass lost compromising the numerical results.

Effect of mesh adaptation on volume conservation

In this section, we will focus on the conservative properties of the level set method and its sensibility to the mesh discretization in the case of compression molding, by defining the simple compression of a cylinder in 2D and 3D. We solve a standard Newtonian Stokes problem and compare the mass (or volume since density is constant) of the part during the compression. We set the compression a constant time step and a fixed velocity of the punch. the initial thickness of the sample is 6 mm and the initial diameter is of 50 mm. The compression ends at a thickness of 2.6 mm. The initial and final state are pictured in Figure 3.18a. The point to illustrate in this section is the conservation of the volume of the piece during the compression. We define the volume loss as the difference of volume compared to the initial volume divided by the original volume, defined as:

(a) Simulation -t = 0 (b) Simulation -t = t f inal (c) Simulation -t = 0 (d) Simulation -t = t f inal
error % = V o -V (t) V o × 100 (2.31)
Generally, we observe a relative loss of volume related to the non-conservative convectivescheme of the level set formulation and the numerical diffusion due to the re-meshing routine. When re-meshing, the level set field on the old mesh is linearly interpolated on the new mesh, which induces numerical errors related to the interpolation. Even for the isotropic mesh case, (without re-meshing) volume is lost along the deformation process. In Figure 2.8, we plot the volume lost for four different cases. The first case concern the use of an isotropic mesh, without activating the mesh adaptation tool. We notice that along the whole compression the volume lost reaches 0.7%. However, in order to keep the accuracy of the method define by a minimal mesh size h min , the same everywhere, we have fixed the number of nodes required 400k nodes to perform this simulation. Mesh adaptation allow us to adapt the mesh in the regions of stronger gradients (dynamically) more important. Three cases, with different number of mesh nodes, where chosen to study the influence of mesh sizes.

For 8k, 16k and 32k nodes, the volume lost of the same test case is drawn in Figure 2.8. We notice that the finer the mesh discretization the lower the error in the volume during the simulation. Additionally, we notice that the rate of volume loss increases when the thickness of the sample decreases, probably, due to the high gradients of the velocity profile in such narrow space. There a adaptive time step will then improve this resolution.

A more challenging configuration (to be chosen for the industrial case) has been used as reference study and is presented in Figure 2.9. In this case, the punch compress the preform until the mold cavity is fulfilled and this case will be studied during the development of this manuscript. Hereby, we focus on the initial test, without any optimization, where a considerable mass loss was encountered.

Figure 2.10, illustrates the final deformation of the preform, once the punch reaches its final position. The final configuration pictures how the preform (purple) did not fill the mold cavity (The expected filled regions is visualized in green). The plot of the volume evolution during the simulation shows a remarkable loss up to 30%. As visualized in the image, the compression of flat geometries with a ratio length/thickness l/d >> 1 induces high velocity in the flow front while compressing the part due to the incompressibility condition. In Figure 2.11, the accumulative loss error and the relative loss ∆error (loss of the current time step) is presented during the deformation. We notice also the tendency of the local loss to increase as the thickness becomes smaller.

Basically this lost is due to many factors. At first point, the not-divergence free of the velocity profile from the numerical discretization during the finite element resolution. The level set convection equation that by definition is not presented in a conservative form. The mixed velocity profile computed in the domain (compressible for the air, incompressible for the preform) that can cause a semi-compressible profile in the vicinity of the interface. The acceleration of the fluid during the compression at thickness becomes smaller. Another source is then associated to the time step that should be adapted to responds the flow conditions at any time step.

Particularly, we notice that the lost per iteration is lower than 0.1%, implying that the volume loss is mainly due to the accumulation of minor errors, mainly during the convection scheme. This point is further discussed in the next section, as well as, the methods implemented to decrease this volume loss. Even thought, the conservative (or not) approaches seems out of the scope of the frame of this work, we propose in section 2.3 a suitable strategy to deal with this problem on simple and industrial cases.

Volume/mass Conservation on level set methods

Adapted three-dimensional multi-phase descriptions to simulate compression molding processes enables a direct study of complex phenomena, such as thermo-mechanical inter- actions. In such a context, topological changes and phase interactions are difficult to be reproduced by traditional Lagrangian methods. In these cases, immersion methods and more specifically level set methods [13] are a more suitable approach a seen previously. Re-initialization of the distance function is required to maintain the unitary gradient of the level set function, by solving a Hamilton-Jacobi equation as in [14]. To accelerate computations, the convection step can be modified to integrate the re-distancing in this single equation [15,10].

In general level set methods framework are better than other techniques fo rimmersion domain such as the Volume-Of-Fluid (VOF) [16], since in those methods mostly an interface reconstruction from the element volume fraction is required. Despite, the benefits and ease of use of the level set method, a lack of conservative properties inherent to the approach is found, making the task in multiphase calculations flows difficult since the mass/volume of the body immersed is continuously degradated due to diffusive numerical errors. To overcome this problem, several authors [17,18,8,19,20], have adjusted the methodology proposing a conservative formulation, enabling the use of recovery algorithms to correct the level set after the convection step, keeping not only the unitary gradient, but also the conservation of the mass/volume. More recently methodologies, using fast marching methods [21,22], improved also the accuracy on the normal computation from the smooth distance, enhancing the convection and re-initialization step. The drawback of such methods rely on the use of numerical parameters that should be calibrated.

While adjusting the unitary gradient of the level set is needed to preserve the profiles and improve conservation, the excessive use of it deteriorates the solution by compounding errors as stated in [9]. Such effect is explained by the diffusion term presented in the conservative re-initialization strategy where a fictitious viscosity is used to smooth the surfaces [18], [8], [20]. The starting point to minimize such issue consist to only 2.31) in volume as a function of the sample thickness during compression of an industrial piece using the volume recovery strategy re-initialize regions where the unitary gradient has been modified. Sato [23] relates the amount of local re-initialization on purely flow kinematics. However, in the work of Owkes and Desjardins [24], a calibration metric based on the amount of deformed level set was proposed and they included the effect of the numerical diffusion errors. Such work was enriched by the contribution of [9], improving the regions where very small level set deformation was found, proposing a more localized re-initialization strategy. Such strategy is based on the normal velocity of the interface and the amount of deformation along the vicinity of the surface.

In this work, we present a particular extension of traditional level set approaches to multi-phase flows, with an improvement given by a smoothed level set, supported with a linear region in the vicinity of the zero iso-value, and the use of mesh adaptation coupled with an edge-error estimation metric in order to guarantee a precise description of the gradients. Furthermore, a strategy for conserving the volume after the advection step has been also implemented, where we attempt to avoid the use of problem-dependent parameters and we try to propose a more stable algorithm using a Dirac computed from the distance function.

Objective

The mass of a object occupying the region Ω * (t) = {x/φ s (x) > 0}, and having a density ρ s is given by:

M (t) = Ω * ρ s dΩ = Ω ρ s H(φ)dΩ (2.32)
the convection of the mass implies, for any domain Ω * ⊂ Ω:

∂ ∂t Ω ρ s H(φ)dΩ + ∇ • Ω ρ s H(φ)udΩ = 0. (2.33)
It is important to notice that both eq. (2.6) and eq. (2.33) are valid representations of the interface dynamics for the entire domain [25]. We label eq. (2.6) the levet set advection and eq. (2.33) the mass advection. Both equations are respected if the level set function is properly convected.

The question addressed in this section relies on keeping the conservative properties of the level set functions when solving their advection using finite element methods.

It has already been noticed that the finite element solution of the level set advection does not necessarily guarantee a conservation of the mass advection [20]. For that reason, many works have been performed in order to address this issue [18,8].

Background in conservative strategies

Correction of the level set function to guarantee conservation can be performed in many different ways, mainly by solving an extra equation to find the necessary factor to apply to the already convected level set φ.

Conservative methods

Kees [20] introduced a strategy in which a normal convection-reinitialization of the level set is firstly done using eqs. 2.6 and 2.8 or even any other reinitialization method. Then, the convection of a conservative Heaviside function, denoted by Ĥ, is done using the mass advection equation eq. (2.34):

∂ Ĥ ∂t + ∇ • Ĥu =0 (2.34)
Finally, the author compared the computed Heaviside function H(φ) to the conservative Ĥ. Then, they suggested to find the conservative level set φ c associated to Ĥ using eq. (2.35) such that the discrete smoothed Heaviside obtained through the corrected level set φ c eq. (2.36) equals the conservative Heaviside obtained by the solution of eq. (2.34).

H(φ + φ ) -Ĥ = κ∆φ ∇φ • n = 0 on ∂Ω (2.35) φ c =φ + φ (2.36)
Two remarks can be given to eq. (2.35). The parameter κ avoids the level-set displacement to be a constant along ∂Ω, thus allowing a mapping of corrections i.e. the level set is displacement depends on the space. In addition, k diffuses the new level-set correction guaranteing a stable and smooth conservative level set. However, this parameters needs to be calibrated and can induces to much diffusion during the numerical resolution fo the conservative level set.

One other strategy has been introduced by Olsson [18], using directly the discrete smooth Heaviside as the level set eq. (2.37). In their work, the resolution of the mass advection, eq. (2.33), is used with the variable φ, following eq. (2.38):

φ =H(φ) (2.37) ∂φ ∂t + ∇ • (uφ) = 0 (2.38)
Eq. (2.38) presents conservative properties inherent to the mass equation. The metric properties are retrieve, solving a the following equation:

∂φ c ∂τ = ∇ • (∇φ c • n 0 )n 0 -φ c (1 -φ c )n 0 (2.39)
the conservative re-initialization recommended by Olsson and Kreiss solves the redistancing of the level set by using a pseudo time τ and a diffusion term depending on a fictitious viscosity . The re-initialization step, eq. (2.39), is performed, conditioning that the modification of the distance function in the vicinity of φ = 0 does not modify the location of the interface in discrete schemes. This is achieved by guaranteeing that the flux on the interface Γ is zero (∇φ φ=0 • n 0 = 0). This condition is respect if the choice of the normal n 0 in eq. (2.39) [8], chosen at the initial pseudo time τ = 0 eq. (2.40). The upper index " 0 " is referred to the normal at pseudo time τ = 0, more specifically:

n 0 = ∇φ τ =0 |∇φ τ =0 | (2.40)
For simulations were the contact angle in the surface modifies a certain response in the mechanical model, one can follow the work of [26]. Therein is presented the modification of eq. (2.40) by adding another diffusion term in the tangential direction as well as a smoothing technique to compute the normal n 0 . Both strategies hereabove have been successfully applied in several benchmarks and examples [18,8,20], both providing a local correction of φ. However, both strategies demand a viscosity-diffusion value, or κ. The tuning of these values might be problemdependent. The higher the value of such diffusion, the less-sharper geometries will be obtained. Also, they imply solving either an extra differential equation or well diffusing the solution each time increment. At the same time, the authors have not addressed this problems under adaptive anisotropic mesh techniques. In the next section, we present a simple correction strategy using a Newton-Raphson algorithm which does not require the resolution of other differential equation and more important, the tuning of new numerical parameter.

Proposed strategy

In order to apply the algorithm, the initial mass M 0 is firstly computed at t = 0:

M Ω * 0 = Ω ρ 0 H(φ(t = 0))dΩ (2.41)
Consequently, we can define the mass at any instant as:

M Ω * = Ω ρ(φ) H(φ)dΩ (2.42)
Here, the general expression of the density is used since it is time-space dependent in our problems. After convecting the level set function φ, the mass/volume embedded by the iso-value zero might differ. At any t > 0, the mass in eq. ( 2.41) and in eq. (2.42) enclose errors due to the discrete resolution of the convective scheme. Defining f as the difference between the current mass respect to the initial, one gets:

f (φ) = Ω ρ(φ) H(φ) dΩ -M Ω * 0 . (2.43)
For obtaining a φ c that satisfies f = 0 we use an iterative Newton-Raphson algorithm. Then at iteration n+1, φ n+1 is given by:

φ n+1 = φ n - f (φ n ) f (φ n ) (2.44)
being n+1 the new estimation of φ, and φ c = φ n+1 when f (φ n+1 ) = 0. This algorithm is applied after convecting the level set function using eq. (2.6).

The expression given in eq. ( 2.44) requires the computation of f (φ):

f (φ) = df (φ) dφ = d ( Ω ρ(φ) H(φ)dΩ)) dφ (2.45)
Denoting with the derivative with respect to φ and using the property

d dφ g(φ)dΩ = ∂ ∂φ (g(φ)) dΩ (2.46)
we get:

d ( Ω ρ(φ) H(φ)dΩ)) dφ = (ρ (φ) H(φ) + ρ(φ) H (φ)) dΩ (2.47)
Where ρ and H are obtained by taking into account the fact that:

dH(φ) dφ = δ(φ) (2.48)
where δ(φ) is the Dirac function associated to the Heaviside function defined by eq. (2.22). The variation of the density with respect to the levelset φ, finally is given by:

dρ(φ) dφ = dρ(φ) dH(φ) dH(φ) dφ = dρ(φ) dH(φ) δ(φ) (2.49)
we recall that ρ(φ) is defined as:

ρ(φ) = ρ s H(φ) + ρ ∞ (1 -H(φ)) (2.50)
where ρ s is the density of the subset Ω * and ρ ∞ the density of the exterior media surrounded by ∂Ω * . After derivation, with respect to φ we obtain:

dρ dφ = (ρ s -ρ ∞ ) δ(φ) (2.51)
Getting the two contributions together, we have:

d ( Ω ρ(φ) H(φ)dΩ)) dφ = (ρ s -ρ ∞ ) δ(φ) H(φ) + ρ(φ) δ(φ), (2.52)
and taking a common factor, by using the definition of ρ and noticing that:

ρ s H(φ) -ρ ∞ H(φ) + ρ s H(φ) + ρ ∞ -ρ ∞ H(φ) = 2ρ -ρ ∞ , (2.53)
we get:

d ( Ω ρ(φ) H(φ)dΩ)) dφ = (2ρ(φ) -ρ ∞ ) δ(φ) dΩ (2.54)
In this expressions, ρ is the density space function, ρ ∞ is the density of the external bodies {ρ ∞ = ρ(φ = -E)} and δ(φ) is the Dirac function associated to H(φ) defined in eq. (2.23). Finally, the Newton-Raphson algorithm for correcting the level set φ states:

φ n+1 = φ n - ρ(φ) H(φ) dΩ -M 0 (2ρ(φ) -ρ ∞ ) δ(φ) dΩ (2.55)
Three remarks may be outlined and must be taken into account. Remark 1 : Ihe case where ρ(φ) = ρ ∞ = ρ cte , the algorithm represents also the volume conservation.

Remark 2 : No assumption regarding the density ρ or well the Heaviside function H(φ) were performed. However, the choice of δ(φ) needs to fulfill the condition dH(φ) dφ = δ(φ).

Remark 3 : The correction algorithm takes always the initial mass as reference to convergence, avoiding accumulative errors due to numerical discretization. Similarly, the non-convergence of one time increment will not modify the stability, proving robustness to the method.

Applying this algorithm requires a good description of the Heaviside and of the Dirac function in the discrete scheme. For that reason, we use the anisotropic mesh adaptation tool based on an edge error estimation. Such topic has been subject of discussion in, Chapter 1, section 1.3.

Strategy to avoid penetration of recovered-level-set into tools

For compression cases, is we displaced the level set to the space belonging to the punch and die, we get penetration problems. This section address this issue and how we handle it to displace the level set after the correction and avoid such penetration. The methodology described in section 2.3.3 displaced the level set contour φ keeping its mass constant after the convection scheme. The resulting φ c however, might penetrates other bodies in the case of multi-body immersion physics. This section develops an strategy to force level set displacement on region of free domains. Thus, enabling the displacement to be carried out avoiding penetration problems. For our multi-body problems, the discrete space where the air is defined is considered as the free domain. The punch and the mold are considered as rigid bodies and not displacement towards their sub-domain is allowed. In Figure 2.12, φ c np represents a conservative level set which does not penetrates any tool. The displacement, independent of the spatial position and obtained by the Newton algorithm is modified, to define a spatial displacement (not constant) ∆φ free define only is particular regions of the domain. Such displacement is obtained by adding a f free variable which in the air and the composite equals 1 and equals zero in the tools.

Punch

We can use the Heaviside function of the mold and punch to define the tools space, by creating a variable of the free surface f free , which states that eq. (2.56):

f free = 1 -H(φ tools ).
(2.56)

Notice that such function equals 0 if H tools = 1 which means that φ is in contact with the tool and positive at the free surface.

Adding a factor f free allows us to control to where is being displaced the level set function and to avoid penetration problems. Nevertheless, in the recovery algorithm, states that the convection of the function is performed everywhere in a homogeneous way. In fact, when including the tools, in the system, some regions will not displaced, since the surface in contact do not belong to the free zone. Thus, causing the total mass recovery not to be fully compensated. In order to solve this, it is included a surface factor aiming to take this effect into account by allowing larger displacements in the free regions. This factor relates the difference between a bubble in the free space correction or a compression molding simulation. In Figure 2.12 the blue line represents the conservative non-penetrated level set φ c np . As illustrates in this figure, the displacement is larger. To recover the proper amount of mass loss, we use the factor f surface given in eq. (2.57). This factor also represents the ratio between the full surface of the level set and the surface sharing a free displacement.

f surface = δ(φ)dΩ (1 -2H(φ tools )) δ(φ) dΩ (2.57)
For example, if the body shares the half of its surface with the tools such factor equals 2. Meaning that a double displacement needs to be performed in the free surface to compensate the real mass loss. Following this reasoning, the expression given in eq. (2.55) is then modified to include the local free surface factor f free and the surface ratio f surface obtaining:

φ c np = φ -f surface f free ρ(φ) H(φ) dΩ -M 0 (2ρ(φ) -ρ ∞ ) δ(φ) dΩ (2.58)
Expression given in eq. (2.58) recovers the mass loss during the convective scheme of eq. (2.6) including the tools. Now one can define the local ∆φ free displacement by means of eq. (2.59) and by knowing that the ρ ∞ in the free domain equals ρ air , one may write:

∆φ free = f surface f free M 0 -ρ(φ) H(φ) dΩ (2ρ(φ) -ρ air ) δ(φ) dΩ (2.59)

Maximum volume loss per iteration guaranteeing accurate displacement

When solving the recovery algorithm the level set displaced should remains inside the fine mesh zone (FMZ). We shall thus guarantee:

∆M (2ρ(φ) -ρ ∞ ) δ(φ) dΩ ≤ 3E -e 2 (2.62)
which allows us to evaluate the maximum mass loss to be recovered inside a well described meshed zone by eq. (2.63):

∆M max = 3E -e 2 (2ρ(φ) -ρ ∞ ) δ(φ) dΩ (2.63)
In the monophasique incompressible case, we obtain, eq. (2.64):

∆M max ρ = ∆V max = 3E -e 2 δ(φ) dΩ Surface moved (2.64)
Including the non-penetration condition (if tools) we get eq. (2.65):

∆V max = 3E -e 2 f surface δ(φ) dΩ (2.65)

Exact displacement of the hyperbolic tangent level set

The euclidean distance associated to the level set function is defined always parallel to the normal of the interface, given by the direction of the gradient of the level set gradient. Mathematically, a parallel displacement of a function f (x) is obtained by applying f (x + δ).The displacement of the hyperbolic tangent level set φ to obtain the conservative non penetrated solution φ c np is given by eq. ( 2.66) :

φ c np = φ + E tanh δ free E 1 + φ E tanh δ free E , (2.66)
To Sum up: the complete algorithm to recover the mass/volume loss is given by solving eq. ( 2.58) to compute the displacement and then by using eq. ( 2.66) to compute the conservative value of the modified level set function. But we remark that we need to solve eqs. (2.59) and (2.56) in order to take into account for the presence of the tools. In the following section, we present the case of the reference geometry introduced early in this chapter in order to prove the methodology for mass conservation of the method. For sake of simplicity, we considered constant density of the part, getting by consequence the conservation of the volume. The strategy explained in section 2.3 is tested in this section for the case presented in Figure 2.9. The resolution of the advection problem for the level set displacement of the preform induces a volume loss reaching up to 35%, as illustrated in Figure 2. 19. The strategy evoked in the previous section required the computation of a Dirac function. In order to ensure the convergence of the algorithm, this function should be well described by the discretized mesh. The first point addressed is the definition of the fields assigned to the multi-field mesh adaptation. Then, the features of the problem are defined, like the velocity of the punch and the numerical specifications according to our framework. Later, snapshots of the deformation process are presented during the part compression. The evolution of the volume during compression is presented, compared to the original case, proving the enhancements of the aforementioned algorithm.

Application of the conservative strategy in an industrial geometry

Mesh adaptation on multiple fields

The theory of mesh adaptation have been firstly extended to the case of multi-field adaptation in the works of Gruau [27] and Nguyen [28]. In [11] several examples of mesh adaptation respect to the velocity field have been shown, while in [29] adaptation with respect to a level set field has been addressed. In our numerical schemes, it is necessary to adapt the mesh to several fields. A L 2 norm is computed to adapt the mesh with respect to the distribution of the error on each field, by proposing that the equivalent error êij is given by the L 2 norm of all the particular errors e ij associated to each field:

êij = L 2 : e ij χ φ φ E + e , 2 χ δ δ, χ v arg max (| v| , v min ) (2.67)
Figure 2.15: Mesh adaptation during compression, with the contour along the piece.

The factor χ x represents the weight on the error contribution of each field x to overall the metric calculation. In here, we do not intersect the metric of each field but, instead we add the contribution of each field pounded by the coefficient χ in order to obtain a global scalar error associated to each edge. Figure 2.14 shows the initial adapted mesh according to the level set functions of the preform, punch and mold and the Dirac function associated to the preform. By iterative correction, starting from a isotropic box, the mesh adaptation tool reaches the discretization observed in 2.14a. We notice the particular refined thickness along the preform, where a Dirac function has been defined. In 2.14b, the zoom in the corner of the preform shows the two levels of adaptation: one supported on a thickness E following the level set strategy, a second supported on the thickness e of the Dirac function.

Compression Molding

Compression Velocity 1 mm/s Initial Preform Thickness 6.2 mm Table 2.1: Velocity and initial thickness of the preform. Figure 2.15 provides a view of a section of the preform, highlighting the finer regions arising during the simulation. The mesh follows dynamically the compression process, until the part has attained its final deformation. The images show the level set contours wero of the three bodies, even if all bodies are embedded in a single volumic mesh as pictured in Figure 2.14.

After these definitions concerning the mesh adaptation strategy, we proceed to explain the specifications of the problem hereafter presented.

• Immersion of the three geometries in a box. (1m × 0.28m × 0.4m

Test Conditions

An incompressible SMC material, having an initial geometrical configuration (preform) is compressed to fill the cavity (mold) by the motion of an upper die (punch). Informations concerning velocity of the punch and on the initial thickness of the preform are given in Table 2.1. In Figure 2.16, the discretization is provided, as well as the mesh size and the level set definition parameters used in the simulation.

The preform is deformed until filling the mold cavity. Figure 2.17 shows the error in the mass/volume, enhancing the improvements of the presented algorithm in the industrial piece.

V o -V V o × 100
Figure 2.17: relative error in volume as function of the sample thickness during compression of industrial piece using volume recovery strategy

The punch deforms the material until a thickness of 2.3mm is reached. In the following, only a section of the 3D geometry is used. The reduced part is shown for sake of simplicity. First, different snapshots illustrate the filling of the mold. The evolution of the part is given in Figure 2.18. There, the algorithm of mass conservation was activated and the final volume lost registered was 0.4%. Compared to the lost without mass conservation algorithm of 35%, a clear improvement is obtained by the algorithm proposed in section 2.3.3. The time step for the simulation reaches a minimum of 0.02s. The velocity in the flow front increases with the thickness reduction. In Figure 2.19, the comparison of the full geometry is given. The compression of the piece at the beginning of this work, and the compression of the piece after implementation of the recovery Newton algorithm. Despite the ease of the implemented solution in this work, the improvements visually are remarkably. The algorithm prove to be robust and of easy implementation. A more adapted strategy for future works should be used working directly in the level set convection step.

short-shots comparison

In Figure 4.29, a comparison of the simulated part with partial-filled molded parts test performed in Plastic Omnium are provided. There, a different thickness of the the preform deformation is compared the flow front position. Basically, this information allow us to verify the filling of the part. The four comparisons presented prove the good agreement of the filling of the piece with the numerical simulation.

The force during the compression of the piece remains one of the quantity to compare with the experimental data. In immersion methods, more precisely supported on a level set strategy, the computation of the force is not a evident task. The upcoming section deals with a strategy adopted in this work in order compute the compression force for any given geometry. Test on a reduced volume of industrial geometry.

Stable and robust strategy for the computation of compression force in immersion domains

For immersed body approaches the computation of fields implying a surface integral in 3D and a line integral in 2D, is not a evident task. The representation of implicit surfaces by implicit functions has its advantages that were discussed early in this chapter. However, in order to compute scalar macroscopic variables such as the compression force to deform a body, surface integral appears of the type:

σ n dS, ( 2.68)
being σ n any variable to integrate along a surface S ⊂ V, in here the normal stress. Several authors addressed this problem by assigning a Dirac function supported on the level set function. In [30] , [31] , [32] , [33] , [34], [35] many examples of Dirac functions, strategies to reduce the numerical oscillations are found. In particular, in [36] and [37] the computation of the force in fluid dynamics problems was addressed, by the resolution of a linear system. In [38], the force is determined explicitly by the resolution of a penalization strategy to impose boundary conditions. The force is then computed at the interfaces. Eventhough all these strategies rely on the level set method, or well in immersion body methods, they are susceptible to numerical oscillations. In the level set approach supported with mesh adaptation was proven that the major error of the method yield on the interfaces. In brief, the force is computed in the regions or major uncertainty of the method. For our applications, the force correspond to the macroscopic term more important to be compared to the rheological experiences. In this section, we propose a methodology supported on the principle of virtual work transforming a surface integral into a volumetric integral without making appear the Dirac function.

Viscous Power and Compression Force

Deformation and flow material requires energy. For polymers having high viscosity such energy can be related to viscous dissipation. Let a volume V of material to be deformed by a stress σ on its surface. The rate of work W done on the material is calculated by integrating (n • σ • v) over the surface.

W = S (n • σ • v) ds,
(2.69) n being the outward unit normal on the surface and (n • σ) is the force per unit surface. Multiplication with the velocity gives the rate of work per unit surface. By using Gauss Theorem, this surface integral is equivalent to a volume integral taking its divergence in the volume V contained by all the surface ∂V :

S (n • σ • v) ds = V ∇ • (σ • v) dV, (2.70)
the term inside the integral is replaced by the indentity:

∇ • (σ • v) = σ : ∇v + v • ∇ • σ. (2.71)
The second term vanishes for fluid satisfying Stokes equation (∇ • σ = 0) since inertial and gravity are neglected. Hereby however, we consider the general case where the divergence of the stress is provided by a force per unit volume f such as:

∇ • σ = f . (2.72)
Finally, the relation between the surface integral to the volume integral stands for all the surface embedding the volume V :

S (n • σ • v) ds = V σ : ∇v dV + V f • v dV.
(2.73)

The surface integral acts on three main zones in our multy-phases problem. The contact with the Mold at (v = 0), the contact with the punch at (v = v punch e z ) and air/fluid interface obtained by the surface tension (n•σ = 0) hereby neglected. Rewriting expression eq. (2.73) one gets:

v punch punch (e z • σ) ds fc = V σ : ∇v dV + V f • v dV.
(2.74)

Therefore, for the case of compression molding simulations, the compression force f c can be obtained by the relation:

f c = 1 v punch V σ : ∇v dV + V f • v dV . ( 2

.75)

The main interest on this strategy relies on a direct calculation of n • σds by means of a volume integration supported on a mathematical development instead of a discrete scheme.

In Figure 2.21, the force comparison between the numerical approach and the analytic solution [39] for the slip and no-slip case are compared. The example of the compression of a cylinder at constant speed. Along the deformation, the material required more force to be compressed. This force required energy. This is the principle of this strategy. Notice the good agreement between the numerical prediction and the reference solution. The strategy reduce remarkably the oscillations found in previous work [30,31,32,33,34,35] 2.5 [39] for the slip and no-slip case with numerical predictions of this work.

for the computation of surface integrals.

Additionally when accounting for complex geometries such as the industrial cases, the surface integral demands the projection of the stress tensor on the normal direction. In our methodology, the orientation is given by the direction of the punch velocity. Since the viscous energy computation is a scalar factor, the strategy is robust and non direction dependent.

Conclusions

In this chapter, the numerical framework, namely, the level set method using immersed geometries in a single computation domain was addressed. Additionally, the strategy to account for multi bodies simulations with different physical properties was discussed and the Heaviside side function was defined. We discuss on the level set convection and the reinitialization, in order to preserve the metric properties. For the objectives of this work, three main contribution were addressed in this chapter. We discussed on the computation of a local booster depending on the deformation of the iso-values of the level set. Meaning, the gradient of the normal velocity. This strategy enables a local computation of the booster. We prove by means of the extension of a circle that the local booster proposed enable a wide range of values improving the level set reinitialization step. This point needs a further study for other kind of simulations such as fluid dynamics in other approaches than compression molding. Another point discussed was the conservative properties of standard level set methods. After the literature review, we propose here a solution for the numerical lost of the volume during the convection step. The lost of volume/mass is basically due to the convection step of a non-conservative form of the advection equation and also the interpolation after the re-meshing step. The spatial discretization was proven to be also an indicator of the mass lost. We proposed an strategy based on a Newton-Raphson algorithm that recover the volume progressively based on the principle of free space. We performed a test on a 3D geometry proven the robustness of the algorithm and its easy implementation in any numerical tool For our case, we reduce from 35% of mass lost to 0.4%. Such methodology implies the computation of a Dirac function and can be used for volume or mass conservation since it is the first approach that really accounts for the density variation on its algorithm. The third point addressed here, was the computation of the compression force. Here an strategy based on virtual work principle enable the computation of the force by means of the viscous dissipation. This strategy proof to provide the exact solution for the slip and non-slip case. Additionally, it reduces the oscillations connected to the other methods.

Summary

The study of the rheological behaviour of SMC materials implies the study of the viscous paste joined with the interactions fiber-fiber and fiber-matrix. The presence of fiber induces an anisotropic behavior on the direction normal to the fiber plane. Previous work concluded that the viscous planar isotropic model proposed in the literature conveys the best results. The anisotropy of the materials are more important when increasing the fiber content. So far the contact fiber-matrix have been supposed perfect and then the material mixture was supposed incompressible. However, the higher the fiber content the more difficult to the fiber to coextensively align with the matrix. For the purposes of this work, high fiber content in volume is used on the SMC sample study herein. The lost of contact during the fabrication process of the fiber with the matrix causes the macroscopic behavior of the stack of plies of SMC to behave as a compressible material the first 30% of material deformation. At isothermal conditions, the porosity within the material originally at 30% closes progressively until the material is fully incompressible.

In order to deal with both approaches in this work an unified formulation is proposed, describing the transition compressible/incompressible of the material. This compressible behavior modifies the evolution of the viscosity along the compression as well as other rheological features used to describe the stress during the compression. For the purposes of this work the high performance SMC material (SM-HP) is described by this procedure. The unified methodology enables the description of incompressible SMC as well. Three main cases has been study: the compression in plane strain, simple compression and compaction test.

The parameter α 3 used to described the compressible behavior enrich the incompressible formulation of the Dumont model proposed in 2003. The compressible behavior of SMC has been characterized in the frame of collaboration with the laboratory 3SR in Grenoble-France. All the equations related to the compressible behavior used in this work is a prove of the communications expressed along this Phd in order to proposed a macroscopic compressible SMC rheological model.

Résumé en Français

L'étude du comportement rhéologique des matériaux SMC implique l'étude de la pâte visqueuse jointe aux interactions fibre-fibre et fibre-matrice. La présence de fibres induit un comportement anisotrope sur la direction normale au plan de la fibre. Les travaux antérieurs ont conclu que le modèle isotrope planaire visqueux proposé dans la littérature donne les meilleurs résultats. L'anisotropie des matériaux est plus importante lorsqu'on augmente la teneur en fibres. Jusqu'à présent, la matrice de fibres de contact a été supposée parfaite et le mélange de matière était supposé incompressible. Cependant, plus la teneur en fibres est élevée, plus la fibre est difficile à aligner avec la matrice. Aux fins de ce travail, une teneur élevée en fibres en volume est utilisée dans l'étude de l'échantillon SMC dans le présent document. La perte de contact pendant le processus de fabrication de la fibre avec la matrice provoque le comportement macroscopique de l'empilement de plis de SMC pour se comporter comme matériau compressible les 30 premiers% de déformation matérielle. Dans les conditions isothermes, la porosité dans le matériau à l'origine à 30 %, se ferme progressivement jusqu'à que le matéteriau soit totalement incompressible.

Afin de traiter les deux approches dans ce travail, une formulation unifiée est proposée, décrivant la transition compressible / incompressible du matériel. Ce comporte-ment compressible modifie l'évolution de la viscosité le long de la compression ainsi que d'autres caractéristiques rhéologiques utilisées pour décrire le stress pendant la compression. Aux fins de ce travail, le matériel SMC haute performance (SM-HP) est décrit par cette procédure. La méthodologie unifiée permet également la description de SMC incompressible. Trois cas principaux ont été étudiés: la compression plane, la compression simple et le test de compactage.

Introduction

Modeling fibrous materials requires an accurate description of its anisotropic behavior. The fiber network associated to the thermoset matrix enhances the mechanical properties of the overall composite [1], but also complicates its mechanical response. Nowadays, glass fiber composites processed by compression molding are used as semi-structural parts since they exhibit advantageous mechanical resistance, relative lower weight and are used for the production of car body parts [2]. It is interesting to better understand the fibermatrix interaction when undergoing large deformations, in order to predict the composite behavior in practical applications. To attain this target, improved models considering fiber and fiber/matrix coupling are needed. Many works attempt to study the flow of SMC preform during compression molding [3,4,5,6]. Based on the observations made in those works, their rheological behaviors are addressed but can seem incomplete [7,8,9,10,11,12,13]. In fact, the SMC has been studied numerically using simplified models assuming isotropic pseudo-plastic materials [14,15,16]. The isotropic models therein used are not consistent with the anisotropic microstructure of SMC, since widely spread fibers build a higher-flow-resistance structure in the plane parallel to sheets. The introduction of an anisotropy representation of the stress tensor starts from the 90s [17]. Despite the good agreement with experimental data shown in that work, this model presents some limitations on the stress behavior for different configurations. That has casted doubt on this initial approach and has encouraged further authors to extend this model. After a decade of technological improvement, Dumont and co-workers [18] have proposed a viscous and transverse isotropic model for reinforced composites, fitting the stress levels obtained experimentally. Therein, a rheometer was used to reproduce realistic deformation on SMC sample, enhancing the characterization of such materials under real industrial conditions. Shear-compression tests performed in [19,2] on SMC samples were complemented with the plane strain and simple compression test performed by Dumont [20] on SMC materials for several fiber concentrations. It is found that the SMC materials have shear-thinning behavior which depends also on the fiber concentration. After the anisotropic model, which is the departure point of this work, many works were carried out to compare experimental data to its predictions [21,22,23,24,25]. This model was also able to describe the behavior of other kind of fiber plane materials such as BMC [26]. So far, the viscous isotropic transverse model is the only model which is able to reproduce all kinds of kinematics applied to incompressible SMC materials, namely, the simple compression, the plane strain and the shearing. Therefore in this chapter, we describe the implementation of such rheological behavior in our numerical platform [27] and computations are performed on these three typical flows: simple compression, plane strain and shear. Usually, simulations on SMC are performed considering a single and incompressible phase Figure 3.1. However as explained in the introduction chapter, air is entrapped in the stack of SMC layers as the fiber content increases. This induces a compressible step in the early stages of preform deformation. To take into account for such phenomenon, the model is enriched by porosity evolution term which modifies the diagonal of the stress tensor according to the observations of the porosity evolution during compression experiments [28]. The extended model developed in collaboration with 3SR laboratory (Grenoble), enables the computation of the compressible stress during the porosity closure. A volumetric term evolves during the material deformation turning the compressible material into a fully incompressible one.

In order to deal with both incompressible and compressible behaviors, an unified anisotropic formulation is proposed to perform numerical computations. The implementation of both models in our numerical software is explained in section 3.3 in which both models are assembled in one unified approach. Detailed explanations on all the parameters needed to handle this anisotropic compressible model for SMC materials are given. The characterization of the compressible behavior of high concentrated SMC materials will be addressed in section 3.3.2, the evolution of the fiber volume concentration during the porosity closure is then discussed in section 3.3.3. Then, following a micro-macro extension, it is presented the model for the viscosity of SMC in section 3.3.1. The numerical implementation detailing the coupling with the mass equation is given and its considerations for other anisotropic problems are established. A summary of all the constitutive equations to be solved are addressed and compressible cases are studied. Here-under, computations are mainly made for the compaction case (also called OEDO test) and stress levels are compared from compressible case to the incompressible one. That gives a clear picture of the stress and material evolution during the compressible/incompressible transition. Finally, the classical compression case is studied when the material changes from a densification state (pore closure step) towards a pure squeeze flow motion.

Modeling Incompressible SMC

From their experimental observations [18], Dumont et al. have proposed a macroscopic viscous transverse isotropic model for SMC materials [20] which is considered as continuum media. This constitutive model assumes that there is a plug flow and the fibers spread in an in-plane way. Therefore, fibers are considered to be randomly distributed in the plane along the sheets during the whole deformation. Averaging the fiber directions modifies homogeneously the macroscopic stress, leading to a plane isotropy. A symmetric structural tensor, M , is then built considering the unitary normal vector to the plane of the sheet. The stress tensor σ is written depending on the strain rate tensor D by the form:

σ = α 0 η eq D + α 1 (M : D)M + 1 2 α 2 (D • M + M • D) -p I (3.1)
where α i are rheological parameters depending on the fiber fraction φ f and the orthotropic strain rate invariants; p stands for a partial pressure acting on the material; η eq stands for an equivalent viscosity at a given strain rate. The total pressure is obtained by taking the volumetric contribution of the stress tensor defined by:

p = - 1 3 tr(σ) (3.2)
Taking into account of incompressibility condition tr(D) = 0 and the symmetric form of M and D gives tr(M • D) = tr(D • M ) = M : D. As M is a structural tensor build with orthonormal rotational vector that implies tr(M ) = 1 and the expression of the material pressure is simplified to:

p = p - α 0 η eq 3 (α 1 + α 2 )(M : D) (3.3)
which gives the link between the partial pressure p and the total pressure. Thefore, the arbitrary pressure p assigned in [20] represents a partial contribution of the volumetric stress.

Modeling Compressible SMC

The viscous orthotropic model presented in the previous section 3.1 has been used for volume fiber concentration of 3.5%, 7.1%, 10.8%, 14.7% and 18.8%. For these low fiber concentrations, the anisotropic behavior of the SMC can be assumed incompressible and it is described by the expressions given in [18] and presented in the previous section (Equation (eq. (3.1). However, as automotive applications demands a strong cut down of vehicle weight, the fiber concentration in the SMC reaches up to 24% and 38% in volume.

The fabrication process of the SMC shown in Figure 12 entraps air within the stack of layers. This entrapped air, denoted here as porosity, becomes important when the fiber volume concentration increases. As depicted in Figure 3.2, the porosity captured through the fiber network inherent in the fabrication process, modifies the SMC behavior. Applying a compression to SMC sheet, we notice a compressible step until a critical value from which the volume strain remains constant. The colored tomography pictures show the porosity inside the material before compression and the normal stress profile indicates a modification of the rheological behavior during the compression. For UL SMC (24% of fiber concentration) the porosity reaches 2%, whereas for the HP SMC (38% of fiber concentration) it reaches 30%. This amount of porosity induces a compressible behavior for the SMC. The air is expelled when material deformation occurs. For molding compression, the porosity concentration decreases with deformation. Accordingly, the porosity closure depends on the volumetric deformation and it is independent of the strain rate [29]. Therefore, modeling of SMC for high fiber concentration has to be adapted to take into account of this compressible behavior. With respect to the previous model, it has to be improved thanks to experimental measurements made by D. Ferre-Sentis during his PhD work. By adding a term proportional to the volume variation of the material, here represented as the trace of the viscoplastic strain rate, the stress tensor is written as follow:

σ = α 0 η eq D + α 1 (M : D) M + α 2 2 (M • D + D • M ) + α 3 tr(D)I (3.4)
where the compressibility of the material is taken into account with the term α 3 tr(D)I. This compressible modeling has been found in other works on isotropic materials in which the stress tensor is linked to the velocity divergence [30]. Following the procedure briefly explained in chapter Chapter 1, a different formulation is used in order to facilitate the numerical implementation.

The new contribution to the modeling of SMC materials from the rheological point of view is given by expression eq. (3.4). Note that the anisotropic model presented here enables the computation of SMC compression molding materials as well as problems of fiber reinforced injection simulations.

By comparing the expressions of stress tensor for both incompressible and compressible cases, it follows that the term α 0 α 3 ηtr(D) becomes p once the porosity gets closed. During the densification the α 3 term increases and the volumic variation tr(D) decreases until the incompressible condition tr(D) = 0 is achieved. From a numerical point of view, multiplying two terms having two different behaviors (one tends to infinity whereas the other tends to 0) might cause numerical instabilities. More important, when tr(D) tends to zero the term multiplied by α 3 becomes a partial pressure p which needs to be coupled with the mass conservation equation. In the upcoming section, we propose to define a variable β which allows dealing with both cases and at the same time a direct coupling with the mass conservation equation. Thus, a clear and smooth transition between the compressible case and the incompressible case is recovered.

Unified Compressible-Incompressible Model for SMC

For compression molding simulations, the models presented in eq. (3.1) and eq. (3.4) are needed and we have to deal with a compressible/incompressible transition. In this section, we unify both approaches into an unified formulation which is able to change from one behavior to another by using a single parameter β. We show in the sequel that both compressible and incompressible cases can be written in a compact form already presented in eq. (3.5):

σ = 2ηD + T 4 : D + T 2 D + D T 2 + qI (3.5)
where η represents the equivalent viscosity of composite, T 4 a fourth-order tensor describing the anisotropy of the material; T 2 a tensor also describing the anisotropic contributions of fiber network and q states for all the volumetric contributions presented in the incompressible or compressible models.

In the sequel, the equivalent form for the compressible case is presented, as well as, the coupling of the velocity-pressure problem for the unified formulation.

For the compressible case, we start writing down the expression given in section 3.2, the stress tensor stands:

σ = α 0 η eq D + α 1 (M : D) M + α 2 2 (M • D + D • M ) + α 3 tr(D)I . (3.6)
This expression takes into account the compressibility of material by the term α 3 tr(D). By replacing this term by the variable p, the Cauchy stress can be presented as:

σ = s -pI (3.7) σ = α 0 η eq D + α 1 (M : D) M + α 2 2 (M • D + D • M ) s + α 0 η eq α 3 tr(D)I -pI (3.8) α 0 η eq α 3 tr(D) ∇•v = -p (3.9)
where s is the viscous stress for SMC materials introduced in [18], and p is a measurement of the compressibility (or incompressibility) depending on the factor α 3 tr(D)I. The definition of σ is then changed by this splitting, whereas the relation given for p defines the new velocity divergence condition.

In this point, we introduce the compressibility factor β as a function of the rheological parameter α 3 .

β = 1 α 0 η eq α 3 (3.10)
Which allows to write a relation between the divergence of the velocity with the partial pressure p:

tr(D) + β p = 0 (3.11)
Please notice for α 3 → ∞, β → 0 (=> ∇ • v = 0; incompressible case). This strategy actually unifies the compressible and the incompressible cases thanks to factor β.

Under such construction the system to solve the mechanical problem is given by:

∇ • S -∇p = 0 (3.12) (3.13) ∇ • v + β p = 0 (3.14)
and the total pressure for both compressible or incompressible cases can be then obtained by taking:

p = p - α 0 η eq 3 (tr(D) + (α 1 + α 2 ) M : D) (3.15)
If α 3 → ∞, the fluid becomes incompressible (tr(D) = 0) and p corresponds to the partial pressure due to the incompressibility, whereas the total pressure is

p = p - α 0 η eq 3 (α 1 + α 2 )M : D (3.16)
the Cauchy stress σ is given by eq. (3.7) while the deviatoric stress is given by:

s = σ + pI (3.17)
or well given in term of S:

s = S + α 0 η eq 3 (tr(D) + (α 1 + α 2 ) M : D) I (3.18)
By solving the set of equations described in the system of equations 3.13, we obtain the velocity profile v and a partial pressure p. This construction simplifies the numerical task of building the anisotropic solver. Additionally, we have as well used a more general integration using a fourth-order tensor T 4 . This strategy comes for the idea of generalized the solver up to other anisotropic behavior laws (for instance the Lipscomb Model [31]). The solver build under the frame of this PhD also enable the description of the mechanical motion for injection of polymer reinforced by fibers. The reader is invited to check an extension of this work in Appendix (Chapter F).

Viscosity model for SMC materials

Based on homogeneous rheometry experiments performed on SMC7520-26 and SMC719 formulations, a continuum one-phase rheological model has been proposed in [29]. In order to take into account of the planar micro-structure induced by the fiber bundle network, the nonlinear viscous character of its deformation and the compressible behavior of the part, the composite is seen as a shear-thinning compressible fluid. The viscosity follows a power law and is given by the formula:

η eq = η 0 (φ f -φ c ) 2 D n-1 eq (3.19)
where n is the index of the power law and D eq an equivalent shear rate which takes into account of the anisotropy of our composite. The dependence with respect to fiber concentration is due to the material compressibility (the volumic concentration is directly related to the porosity). Finally, the viscosity of the part increases during the densification by a quadratic expression (See Figure 3.4). Here, φ f stands for the volume fiber concentration of the material and φ c a fitting parameter allowing a better agreement with experimental data. During the compaction of SMC, the fiber volume fraction φ f increases whilst φ p decreases. Then, the viscosity increases since the number of contact fiber-fiber increases. The proposed form for D eq (the equivalent strain rate) is:

D 2 eq = α 0 D : D + α 1 (M : D) 2 + α 2 (D • M) : D + α 3 tr(D) 2 (3.20)
where the viscosity η 00 correspond to the axial viscosity recorded during homogeneous plane strain compression for an axial strain rate of 1 s -1 applied on a SMC sample. The temperature dependency of the viscosity is given by an Arrhenius relation as found in [20] and [26]. For modeling such temperature dependence, we use values found for similar SMCs studied in the literature [24] (See Figure 3.3). The Arrhenius expression states:

η 0 = η 00 e b( 1 T -1 To ) (3.21)
with η 00 the viscosity at the reference temperature and the reference strain rate. The data can be read from Table 3.1 and Table 3.2. The viscosity evolution without chemical reaction stands:

η = η00 e b( 1 T -1 T o ) (φ f -φ c ) 2 D n-1 eq (3.22)
and the evolutions of each terms as function of the temperature and volume fiber fraction are given in Figure 3

φ f (φ f -φ c ) 2
viscosity evolution with the fiber contacts fiber-fiber Contact Figure 3.4: Evolution of the SMC viscosity as a function of fiber-fiber interaction (φ c = 0.02).

The Compressible/Incompressible transition -Evolution of parameter α 3

The parameter α 3 is used to describe the compressible nature of the SMC when it traps air because of the fabrication process. While compressing the SMC sample, expelling the porosity (or densification) occurs in the first stage when there is no motion in the transverse direction. The height decreases while no material flow is observed. When the air is compacted by the deformation, the α 3 term reaches a high number and the materials becomes incompressible. According to [29], the evolution of the constitutive parameter α 3 as function of fiber concentration is proposed on the form:

α 3 = k e -b φ c f + H (3.23)
In Figure 3.5, the evolution of α 3 defined by eq. (3.23) with respect to the fiber concentration is plotted. The volume fiber concentration φ f is a bounded value. For HP SMC, the initial value of the volume fiber fraction equals 28% before the porosity closure. Once the porosity reaches 0%, the incompressible fiber fraction φ dense f equals 38% (additional information are given in Table 3.4). The value α 3 increases meaning that the porosity is expelled. The parameters k, b, c and H are given in Table 3.3. In the mechanical resolution, the evolution of α 3 is taken into account by means of the compressibility term β (see eq. (3.10)) and its variation is also plotted in Figure 3.5. It is positive while the material is compressible and tends to 0 as the porosity decreases and therefore the material becomes incompressible. 3.3. α 3 is the parameter proposed by [29] whereas β is the parameter used for the computations. They are linked by the relation eq. (3.10)

Evolution of porosity and fiber fraction concentration

For the high profile SMC, the initial porosity fraction is around 25% of the volume of the material. During the compression, this value decreases modifying many properties of the material such as viscosity, density and compressibility. By knowing the fiber concentration in the dense state and computing the volume evolution of the sample, the According to the rheological measurements and the phenomenological observations, we can associate the volumetric deformation with the porosity evolution. At macroscopic level and a simple squeeze flow, the volumetric deformation is obtained using:

ε v = ln V V o (3.24)
For our simulations, we describe the volume variation as a local variable. For that the volumetric part of the strain rate tensor is taken and the volumetric variation follows:

dε v dt = tr(D). (3.25)
The porosity is then calculated using an experimental expression eq. (3.26) proposed in [29] and used in previous work [32]:

φ p = - δ 2 (-ε v -ε vo ) -(ε v -ε vo ) 2 + γ 2 + φ po (3.26)
where the fitting parameter δ stands for the slope of the linear transition of the volumetric deformation with the porosity closure; γ is the numerical parameter used to avoid that the volumetric deformation exceeds the critical value ε vcrit . The value ε vo in eq. (3.26) is obtained through the expression:

ε vo = - γ 2 - 2ε vcrit δ 2 4ε vcrit δ (3.27)
Finally the matrix concentration is updated as well as the fiber fraction concentration thanks to this sequence of equations: In Figure 3.6 is plotted the evolution of the porosity, fiber and matrix volume concentration as function of the volumetric deformation ε v . In the fitting presented in eq. (3.26) δ stands for the slope of the curves presented in Figure 3.6 while the material is still compressible. Please notice the volumetric deformation for compression cases are negatives. Once the volumetric deformation reaches the critical value ε vcrit = -0.25 the linear behavior is modified by the convex parameter γ assigned to the concavity of the transition once the incompressible limit is reached. Then, the relative density is obtained:

φ m = φ dense m (1 -φ p ) (3.28) φ f = φ dense f (1 -φ p ) (3.29) or φ f = 1 -φ m -φ p (3.
ρ r = 1 -φ p (3.31)
Notice that the information of φ p enables the direct computation of the relative density and the information of the fiber concentration φ f allows us to update the viscosity of the material η according to our rheological model defined in section 3.3.1.

Additionally, the evolution of the relative density can be obtained by solving eq. (1.18) directly, as suggested in Chapter 1. However, the procedure depicted in this section also enables the updating of the fiber concentration φ f and also the viscosity consistency as required in section 3.3.1.

Summary of mechanical resolution

After having introduced the mechanical behavior of SMC and their constitutive parameters, a summary of set of equations to solve numerically are addressed hereunder. The planar isotropic model for SMC is defined by a structural tensor denoted as M accounting for the normal e 3 of to the fiber network (M = e 3 ⊗ e 3 ). The anisotropic stress tensor is written in its compact form:

σ = 2ηD + T 4 : D + T 2 D + D T 2 + qI (3.32)
where η is the updated viscosity and T 4 and T 2 account the anisotropy exhibits for such materials and they are functions of the structural tensor M . With this formalism, the connection with model described earlier in this chapter is made by taking :

η = α 0 η f (D eq ) n-1 2 Viscosity (3.33) T 4 = 2 α 1 η (M ⊗ M ) 4th order Tensor (3.34) T 2 = α 2 η M Tensor (3.35)
The consistency of the material depends on the local fiber concentration φ f

η f = η 0 (φ f -φ c ) 2 (3.36)
and the equivalent viscoplastic strain rate D eq retrieved accordingly to the expression:

D 2 eq = α 0 D : D + α 1 (M : D) 2 + α 2 (D • M) : D + α 3 (tr(D)) 2 (3.37)
Finally the mechanical problem is given by finding (v, q) such as:

∇ • σ(D, q) = 0 (3.38) tr(D) + β p = 0 (3.39)
For compressible SMCs, the fiber fraction concentration φ f and the relative density ρ r need to be update during the compression molding process. As described in section 3.3.3, those parameters depend on the volumic deformation ε v obtained by solving:

dε v dt = tr(D) (3.40)
That enables the computation of φ p and then φ f , which is used as entry to update the viscosity η and the compressibility β in the mechanical problem.

Numerical cases

In order to check the implementation of the mechanical resolution in our Finite Element Library, several numerical tests were performed. First, the rheological cases for the determination of the orthotropic incompressible viscous model explained in section 3.1are reproduced. There, three cases are presented: simple compression, plain strain compression and shear test. The equivalent stress level and equivalent strain rate are summarized in a single comparative plot. The model prediction is superimposed putting in evidence the validation of the incompressible model. Secondly, the feasibility of the unified model for the case of high concentrated SMC and its compressible/incompressible transition is studied on a compaction case. There, the materials is densified until the incompressible state is reached. The normal and lateral stress are compared to analytical formulas and experimental measurements. Thirdly, the compression of a sample with 38% of fiber concentration is studied, evolving from a compressible case towards a fully incompressible compression. In this latter case, the influence of slip and no-slip boundary conditions are studied. We observe the stress evolution during the compaction and during flow motion. Finally, the plain strain compression of the same material is reproduced and the stress prediction is compared to the ones provided by the compressible model presented in section 3.2. All simulations are made by using the unified model presented in section 3.3.

Standard SMC -incompressible case

Assuming that the compositeis a homogeneous phase (composed of the fiber and the resin paste) as illustrated in Fig. 3.1, the two equivalent components (stress and strain rate) are defined according to eq. (3.41) and eq. (3.42) (deduced from the viscous transverse isotropic model presented in section 3.1):

σ 2 eq = 1 2 (1 + 2H) s : s + (5 + H -6L) (M : s) 2 -2 (1 + 2H -3L) (s • M ) : s (3.41) D 2 eq = α 0 D : D + α 1 (M : D) 2 + α 2 (D • M : D (3.42)
all rheological parameters are constant and the coefficient α i can be expressed as function of H and L (for example, H = 0.446 and L = 15.8 for UL SMC). The tensor s states for the deviatoric stress tensor. The orthogonal direction n represents the normal to fiber network (here, the vertical direction) and it intervenes in the model by means of the tensor M = n ⊗ n. The parameters L, H account for fiber concentration, their determinations and relations are detailed in [18]. The geometries of the three rheological experiments are plotted in Figure 3.7. The mechanical equations solved in this section are presented in section 3.4. Here, the compressibility factor β is set to 0 (incompressible case) and the equation for the volume strain evolution is not necessary. Finally the system solved is:

∇ • σ(D, p) = 0 (3.43) tr(D) = 0 (3.44)
The initial dimensions of studied sample are respectively: for simple compression, a cylinder used having a height of 7.6 mm and a radius of 50.0 mm (Figure 3.18a ); for plane strain, a cube with a height of 7.6 mm, a width of 60 mm and a depth of 40 mm (Figure 3.18b). They are compressed until reaching both a height of 4mm (47% of compression ). For shear computation, we use the same sample as for the plane strain problem.

Experiments presented in [18] gave values of equivalent strain rate D eq and equivalent stress σ eq defined by eq. (3.41) and eq. (3.42) for the three experiments, two different fiber concentrations (3.5% and 14.7%) and three applied strain rates [10 -4 , 10 -2 , 1]s -1 . The composite resulting consistency stated for the two different concentrations η f 3.5% = 0.763 M P as and η f 14.7% = 4.886 M P as Theses experimental values were compared with the analytical solutions and the results of our numerical computations. A total of sixteen simulations were performed with slip boundary conditions. Since the Cauchy tensor is constant in the domain, a space constant stress tensor defines the stress state. The results presented in Figure 3.8 point out a good agreement between numerical computations and the experimental measurements used to characterize the material. Important to mention that the experimental cases hereunder compared satisfy slip wall condition. In fact, such friction-less conditions were reproduced experimentally by lubricating the contact sample/wall. The numerical strategy used to tackle this boundary condition in the framework of immersion methods is detailed in Appendix:Chapter B. We also notice the linear tendency in this log-log plot. The reference value of the stress at D eq = 1s -1 in this curve represents the viscosity η 0 that depends on the fiber concentration φ f which is assumed constant for incompressible SMC. The slope of the line stands for the viscoplastic strain rate dependency index n, here n = 0.44. Clearly, the viscosity of the mixture increases with fiber concentration.

In this section, We prove that for the incompressible model the numerical computations predict equivalent stress correctly. It is worth notice that the equivalent stress tensor eq. (3.41) depends only on anisotropic parameters H and L. As proved in [20], the incompressible model fits very well experimental data and by extension our numerical results as well. In the next section, a numerical test is carried out using the unified model to predict the stress evolution during the compressible step of compression of SMC sample.

High performance SMC -Compressible case

In the following, the compressible case is studied when high performance SMC is modeled. We present the analysis of compressible SMC for three configurations: a purely compaction case, a classical squeeze flow and a plane strain compression.

Compaction of HP SMC -Oedometric Compression

High performance SMC materials exhibits a compressible step during compression molding. The dedicated configuration to check this compressibibility is the compaction set up pictured in Figure 3.9. In this apparatus, the sample is blocked inside a circular cavity and it is deformed until the material becomes incompressible. The cavity surrounded the SMC sample restricts its radial movements. Therefore, only vertical deformation is allowed and the SMC porosity decreases progressively during the compaction process . The internal pressure in sample increases while porosity closes. That translates into an increment of the stress needed to deform the sample. A constant strain rate deformation ḣ/h is imposed and the evolution of the stress during vertical deformation is recorded, When the level of porosity is low enough, the material behaves as an incompressible flow and the stress tends to infinite. In the experimental test, the compression ends when the machine limit compression force is reached. The vertical displacement is directly related to the macroscopic volumetric deformation by means of relation eq. (3.24).

The compaction test (also called densification test) has been performed on two samples of HP SMC. The normal and lateral stress level were recorded. The normal stress is obtained by dividing the force of the machine to the surface of the sample. In this case, the surface of the sample does not change as only axial deformation is allowed. Note that the normal experimental stress presented here is obtained by assuming an uniform stress level within the material. An approximation that might not be truth at the extremities of the sample as shown in the numerical simulation.

THe initial sample has a diameter of 120 mm and a height of 9.5 mm. To get an uniform top surface a pre-compression of 200 N is carried out. This stands for the precharge in compression molding process for industrial parts. The machine is set to respect a constant strain rate compression v/h (v being the punch velocity ) and it stops when the compression force reaches 90 kN . The compression test conditions are described in Table 3.5. The sub-index 1 or 2 are used along the text to identify the experimental conditions. The immersion of the sample mold and punch is pictured in Figure 3.10. As observed in Figure 3.11 the volumetric deformation ε v equals the deformation rate ε. During this densification, the porosity volume fraction decreases and the volume fiber fraction increases as well as the relative density ρ r . For the compaction case, we notice also that the porosity fiber fraction does not follow a linear tendency.

Experimental setup for HP SMC Case

v/h T imp h 0 h f F max 1 0.003 s -1
Hereunder, a spatial description of the porosity closure is pictured in Figure 3.12. During the densification, the porosity evolves depending on its position. The porosity in the center tends to close before the porosity close to the external radius. This nonuniformity is due to boundary conditions at the mold/sample contact. In Figure 3.12a, the sample contains a homogeneous initial porosity fraction of 25%. During the compression the volumetric deformation follows the equation eq. ( 3.40) and the material is densified. Its porosity along the sample is obtained thanks to the relation given in section 3.3.3. In Figure 3.12b and Figure 3.12c, we noticed that the porosity after some deformation reaches 10% in the center and is about 15% at the mold lateral contact. In Figure 3.12d the material porosity is lower than 0.2% and the simulation stops. In Figure 3.11, the volumetric deformation ε v and the volume concentration of the fiber matrix and porosity evolution during the axial deformation is presented. As expected, we observe that the volumetric deformation is linear with respect to vertical deformation until a deformation larger than 0.2. Afterward, a change in the slope of the curve is observed since the porosity begins to be low enough and the material becomes incompressible. The pressure increases during the porosity closure since the material becomes incompressible (See Figure 3.13). According to the model, the pressure departs from a stable value 3.13a and increases during deformation Figure 3.13b and Figure 3.13c. Logically, the lower the porosity inside the material, the lower its capacity of being compressed, the higher the pressure. In 3.13d a snap of the pressure is presented before the simulation stop. In order to avoid numerical problems once the material is fully incompressible, the simulation is stoped when the porosity is lower than 0.1%. Consequently, the velocity profile presented in Figure 3.14 moves from a compressible profile as in Figure 3.14a towards an incompressible pattern visualized in Figure 3.14d characterized by the squeeze flow profile [33,34]. From a linear profile in the vertical direction and zero-radial velocity as depicted in Figure 3.14b, it moves to a standard squeeze flow profile at the incompressible state. The modification of the compressible profile towards the incompressible one occurs at ε = 0.22 (see Figure 3.14c) where there is the compressible/incompressible transition. When the materials becomes incompressible (or almost) the immersion method meets difficulties to impose the velocity inside the punch and the mold.. The pressure increases until a point where the viscosity affected to the mold (hereby chosen as 1000 times higher than viscosity of SMC), is not sufficient to fulfill its role. The strain rate on the mold side reaches the same order of magnitude than for the SMC.

Comparison of stress level with experimental data

According to our model, the normal and lateral stresses during the compaction test (compression OEDO) states (See details in section D.2.3):

σ L = η 0 (φ f -φ c ) 2 α 3 (1 + α 1 + α 2 + α 3 ) n-1 2 D n 33 (3.45) σ n = η 0 (φ f -φ c ) 2 (1 + α 1 + α 2 + α 3 ) n+1 2 D n 33 (3.46)
where η 0 is the viscosity of the HP SMC given in Table 3.1. The volume fiber concentration φ f changes during the simulation and modifies mainly the coefficient α 3 . The constitutive parameters α i are described in Table 3.2 and section 3.3.2. Assuming constant strain rate during the deformation, these analytic stress levels are plotted in Figure 3.15 and Figure 3.16. They are compared with the numerical stresses provided by the numerical simulations and the experimental measurements performed on samples of SMC materials. Theses comparisons are made by looking at the evolution of the stress during the vertical deformation for the two strain rate velocities summarized in Table 3.5. The behavior of the stress is the same for both velocities. There is only a modification of stress level since the model is of purely viscous in nature. The analytic model and the numerical simulations predict the same stress level as pictured in Figure 3. Velocity profile evolution during oedo compression and compressible/incompressible transition: the velocity profile is vertical and linear at the beginng and has a typical squeeze flow patterns at final state.

SMC behavior.

The differences found with respect to experimental data during the densification are mainly due to the difficulty to reproduce equivalent experiments for SMC samples. In fact, it is not surprising to found high dispersion (until 50%) with measurements performed for SMCs. This dispersion brings problem when deciding the reference curve to set the model parameters values, namely, to characterize the mechanical properties of the material. In addition, the model proposed for the evolution of the porosity parameter by means of α 3 term remains an open question. The model proposed here supposes an isotropic closure of the porosity since the orthotropic parameters α 1 and α 2 remain constant. The characterization was performed for many strain rates and different temperatures, here we only presented two cases at ambient temperature. Despite the difference between experimental and numerical values, the model proves to follow the evolution of stress level according to the experimental data. Meaning that previous of this work this stress evolution was not possible. It remains as perspective to try to improve the modeling by assigning a orthotropic evolution of the porosity and by determining a better description of the porosity closure in relation with the given values for the viscosity proposed in [29]. From the numerical point of view, our reference curve is the predictions of our model, and we notice a good agreement implying that our implementation in the software was done correctly. The main improvements of the model will be noticed for squeeze flow and plane strain compression.

The stress level using the numerical simulation also reproduced the model predictions for the second strain rate Figure 3.16. There, the model and the numerical predictions are in agreement, whereas the results with experimental data differs. We notice by the curve shape that the lateral stresses behaves differently during the compaction indicating once again the non-isotropic evolution of the stress during the porosity closure. Deformation ε = -ln h h 0

Contraint [P a]

Case 1 The ratio between the two stresses, the lateral and the axial, is a measure of the compressibility and anisotropy of material. Indeed, according to the model, the initial stresses ratio equals the rheological parameter H = 0.446. During the compaction, the lateral stress increases faster than the axial ones. Thus, when incompressible state arises, such ratio equals the unity. The evolution during the densification is strongly related and ruled by the evolution of parameter α 3 (and therefore by its experimental fitting from experimental data). In Figure 3.17 the evolution of both strain rate compression are compared to the one predicted by the model.

σ L Model σ L 2D σ L 3D σ L Exp (
The compaction case enables the verification of the volumetric deformation with the porosity closure. The important value for the molding process is the amount of deformation before it behaves as an incompressible material. And therefore, it is necessary to develop a progressive strategy to deal with the compressible and/or incompressible case when the material is not avoided to flow in real configurations. The simple compression case is a more representative configuration in relation with the compression molding pro-

Stress ratio [-]

Evolution of stress ratio Model σ L /σ n case 1 σ L /σ n case 2 Figure 3.17: Stress ratio σ rr σ zz evolution during compaction test cess. There, the material is free to flow in the radial/horizontal direction. That case is studied as follow in order to evaluate when the material starts to flow, at which axial deformation and how the stress behaves in such case.

Compressible/Incompressible Transition for a Squeeze Flow

For the simple compression case the numerical configuration is shown in Figure 3.18.

The analytic solutions are provided in section D.2.1. Two wall boundary conditions are studied: slip contact and sticky (no-slip) contacts. In molding compression applications the friction between the mold wall and the sample modifies the mechanical behavior and therefore the force applied on the punch.. For a slip contact, the material is subjected to a purely elongational flow whereas a sticky contact causes that the flow is dominated by the shear. The friction condition is a combination of both cases. The most predominant kinematic in SMC compression molding is an equivalent of the squeeze flow and it is studied by numerous authors [35,36,37,38]. The comparison between analytic solutions and the numerical computations can only be done for the slip boundary condition for which it is possible to exhibit such solution by using continuum mechanics. Therefore, further to such studies, we also present in this section the case for a non-slip boundary condition. There, we point out the differences in the behavior of the porosity closure and stress level for both boundary cases.

In Figure 3.19 and Figure 3.20, the evolution of the compressible parameter β introduced in this work is presented. In Figure 3.19a the initial porosity concentration φ po = 28%. At deformation about ≈ 25% the material reaches the transition state pictured in Figure 3.19c. We notice than also in Figure 3.21 the volumetric deformation presented in Figure 3.24a has a change in its slope and its tendency becomes flatter. Thus part of the vertical deformation is transformed in surface deformation. When the material reaches the incompressible state Figure 3.19c it flows. A further snap is taken at deformation ε ≈ 0.6 where the material is fully incompressible 3.19d since its porosity concentration in lower than 0.1% as shown in Figure 3.21b.

Accordingly for the slip contact, the analytic solution states as shown in section D.2.1 and the variable α 3 is updated knowing that the trace of the viscoplastic strain rate stands in this case as:

σ 33 = η 0 (φ f -φ c ) 2 α n+1 2 0 1 + α 1 + α 2 + α 3 1 + 2α 3
T r(D) = 1 1 + 2α 3 D 33 (3.48)
In Figure 3.22 the stress for the slip condition is compared to the analytic solution. Firstly, it is important to observe the behavior of the stress evolution for this viscous model. At constant strain rate, normally, an incompressible viscous model predicts a constant stress value. Clearly, the shape of this stress curve along deformation does not behave like that. The evolution of the stress along deformation will depend on compressible behavior and friction condition against the wall. The compressible transition is determined by noticing that the slope of the linear behavior of the stress in the early compression stages (deformation < 0.2) changes abruptly. In the previous section, we have shown that the porosity closure of the material increases the internal pressure making it more difficult to compress. The same behavior is also found during these computations. Then, the material flows and the friction against the mold play an important role. For the slip case, once the material becomes incompressible, the viscous stress remains constant. For the no-slip case, the stress will linearly increases during the deformation since more surface is in contact with the mold, requiring more energy to displace the material. The reference curve eq. (3.47) for the slip case is plotted in dashed line inFigure 3. the two boundary conditions are depicted with blue and red lines. In the framework of immersed domain method, the slip boundary condition is recovered by adding a constant viscous layer as explained in [39], the numerical prediction are represented by the blue line. This numerical stress is slightly higher and increases weakly with the vertical deformation. This latter means that the viscosity of the additional has to depend on the local tangential stress to better satisfied the slip condition. It seems that there is a small friction for large deformation. For the no-slip case, there is no analytic solution but the analyze of curve shape can still be made: the stress behaves as for the slip case during the densification; then the stress level increases as the surface in contact with the mold does. These two curves represent the bounded values of stress and all other friction cases remain inside this range. To check the capability of the developed model and its numerical implementation, compression moldings for plane strain compression (ps) were simulated. For that purpose, an initial SMC sample of dimension w o × L o × h 0 = 40 mm × 60 mm × 6.5 mm located in the center of a rectangular mold having a channel length of 80 mm. The charge is compressed until the mold cavity is fulfilled. The deformation follows the kinematics given in Figure 3.23.

Plane Strain Compression

For this case, we study the evolution of the porosity φ p , the volume strain ε v and the plane stresses during the compression. Only, the no-slip boundary condition is addressed as there are numerical difficulties to set slip boundary condition for the plane compression. The viscous layer approach under this configuration induces penetration of the sample into the mold, since it induces a non-zero velocity of the material towards the mold. There is also the problem of vertical wall. Further strategies should be integrated to deal with the slip case during compression in immersed domain method (Nistche method for example). The molding configuration gives rise to two main stresses: a first one along the axis of compression σ 33 , a second one along the vertical plane of mold in contact with the material σ 22 . The analytic solutions for slip boundary condition are determined in Appendix:section D.2.1 and stand as:

σ 33 = η 0 (φ f -φ c ) 2 α n+1 2 0 1 + α 1 + α 2 + α 3 1 + α 3 n+1 2 D n

33

(3.49) The analytic and computed evolution of the stresses during the plane strain compression have been plotted in Figure 3.25. The materials behaves as incompressible after a vertical deformation of 0.4 until its porosity concentration vanishes. We also notice changes on the shape of curves depicting the stress behavior at this point. The stresses (σ 33 and σ 33 ) given by analytic solutions stand for slip boundary condition and remain nearby constant when the material becomes incompressible. As in [26], the analytic solutions and experimental stresses differ due to friction along the mold. From this comparison, we can see the influence of the friction or well the boundary condition during compression with respect to the prediction of stress level. The compression curves during this experimental test are rather related to the configuration for which there is a high friction against the mold.

σ 22 = η 0 (φ f -φ c ) 2 α n+1 2 0 α 3 1 + α 3 1 + α 1 + α 2 + α 3 1 + α 3 n-1 2 D n 33 ( 3

Conclusions

In this chapter, we addressed the problem of compressible SMC materials and its rheological model. Large volume fiber concentration in SMC materials improves mechanical properties of the rigid part once molded. However, as presented in this chapter the fabrication process entraps air among the fiber bundles inducing a compressible behavior of the SMC. Usally, such compressibility is not taken into account since the porosity is less than 1%. For standards SMC, the model proposed in literature by Dumont can be used. This model is the first model assigned to describe SMC as an orthotropic material. Basically different behaviors stand in the thickness and in the plane of fibers. The numerical simulations described the model are proposed for several molding conditions and fit well experimental data.

However, for the SMC studied in this section containing 38% of volume fiber concentration carries 25% of porosity. The compressible behavior was modeled and its evolution is described by a term related to the volumetric changes in the expression of stress tensor. The model was used to predict the mechanical response in a compaction test. A new formulation for the compressible SMC is here presented followed by the work performed during this PhD as part of the collaboration with the laboratory 3SR in Grenoble-France. The compressible model however, did not show a clear transition between the compressible SMC behavior and the incompressible one. We presented in this chapter an unified model. Such model described by the term of a compressibility factor, β, enables a single formulation for both compressible and incompressible SMC and allows a continuous computation. For the incompressible cases, we predict the stress level of the simple compression, plane compression and shear. The numerical model used described the same stress level of the model proposed by Dumont et at. For the compressible case, the pure compaction computations show the ability to get the compressible/incompressible transition. Further, plane strain and axisymetric compressions are addressed for two boundary conditions at the SMC/mold interface: the slip and the no-slip case. We also point out that the porosity closure can be detected by looking at the slope variation of the stress during the compression.

In conclusion, this model enables thanks to the compressible factor β, an easy computations of dense state of SMC. The model here proposed improves the rheological model, since it is based on a stable numerical formulation when the materials becomes incompressible. The model using the fourth order tensor as input also enables simulation for fiber-reinforced simulations in injection molding cases extending its implementation to further fiber reinforced problems. Compare to current models in literature, no models, to our knowledge, have been introduced to deal with compressible and incompressible SMC. More important, this model enrichs the model presented in [20].

Coupling the compressison of SMC with the thermal exchange between the mold stands forward a more complete modeling of SMC. The next chapter deals with the thermal treatment assigned to SMC materials and we compare such predictions under chemical reaction with experimental data. ature of the sample causes the chemicals components of the matrix mixture to react and trigger a chemical reaction denoted as curing. Such mechanism also causes the transformation of the matrix from a viscous paste state towards a consolidated part considered rigid.

In this chapter a themo-kinetical model for studying the temperature evolution during curing as well as the phase transformation is proposed. The model used is the so-called Bailleul model. The materials has been characterized using a device know as the PVT-α. This work has been performed by the laboratory of thermo-kinetical of Nantes -France (LTN). This laboratory has proposed a fully thermo-kinetical model for the two types of SMC study herein. As complement of this collaboration, the laboratory also provided us with experimental data for 2 different configurations in each one of the two materials here studied. We use the model and test it during the heating of a cylinder sample encapsulated in a hot mold and during the compression in a channel.

Two main cases are in this chapter compared. Before those test all the equations related to the thermo-kinetical model are depicted. The chapter focus in detail on describing the heat source injected into the energy equation that is proposed as a function chemical reaction evolution. the variable used to described the curing state is dneoted as α. So we move from a standard temperature finite elements resolution to a coupled temperature-conversion degree resolution (T-α unknowns). The fist test is a cylinder sample blocked where the temperature of the surface increases by 3 degres per minute (slow reaction). The main objective was to make evolve the temperature of the sample controlling the thickness gradient. After about 1500 seconds the material reacts for SMC-UL and after 1800 second for the SMC-HP. The numerical prediction of the same test proves that the model allows us to fit experimental data based on temperature evolution in the core of the sample as well as thickness evolution during the steps of heatingreacting -and cooling. As complement of this study a sensibility case of characterized properties such as chemical energy, conductivity, heat capacity, thermal contact resistance was performed concluding that the variables the modify the most the reaction time and the temperature peak are the thermal contact resistance and the chemical energy of the paste. For that reason it is important to well characterized such properties for further materials.

A second test case was studied where the material is compresses during heating. This case consisted in model a real compression molding experience where the hot wall of a mold compresses the sample until reach certain thickness. There the molds holds until the chemical reaction takes place and the material is fully consolidated. The experimental data recorded force during the test, temperature in the center of the sample and in lateral positions as well as heat fluxes sensor where placed on the mold and punch surfaces. The comparison of the experimental data with the numerical simulations predictions are in general a nice estimation of this work under real industrial conditions. The temperature evolution for SMC-UL is in good agreement with the experimental data. The core surface temperature and the heat fluxes computed numerically proves that the model and the simulation set describes fairly good the experimental test. The force during the compression however is relatively not so well predicted based on the main reason that the friction against the mold is here not precisely modeled. The same procedure is repeated for the SMC-HP but there the chemical reaction describe numerically are in less agreement with the experience. Finally, at the end of the chapter the prediction of the thermo-kinetical model for a industrial geometry is given showing that for a piece used in Plastic Omnium Auto Exterior the reaction time under the condition of mold temperature around 150Celsius is about 45 seconds.

Résumé en Français

L'étude du comportement thermo-cinétique des matériaux SMC fait l'objet d'une étude dans ce chapitre. Les murs de moules chauds chauffent le SMC en augmentant sa température. La température de l'échantillon fait réagir les composants chimiques du mélange de matrice et déclencher une réaction chimique indiquée comme durcissement. Un tel mécanisme provoque également la transformation de la matrice d'un état de pâte visqueux vers une partie consolidée considérée comme rigide.

Dans ce chapitre, on propose un modèle thérmo-cinetique pour étudier l'évolution de la température pendant le durcissement ainsi que la transformation de phase. Le modèle utilisé est le modèle Bailleul. Les matériaux ont été caractérisés à l'aide d'un périphérique connu sous le nom PVTα. Ce travail a été réalisé par le laboratoire de thermo-cinétique de Nantes -France (LTN). Ce laboratoire a proposé un modèle entièrement thermo-cinétique pour les deux types d'étude SMC dans le présent document. En complément de cette collaboration, le laboratoire nous a également fourni des données expérimentales pour 2 configurations différentes dans chacun des deux matériaux étudiés ici. Nous utilisons le modèle et le testons lors du chauffage d'un échantillon de cylindre encapsulé dans un moule chaud et pendant la compression dans un canal.

Deux cas principaux sont comparés dans ce chapitre. Avant ces essais, toutes les équations liées au modèle thermo-cinétique sont représentées. Le chapitre se concentre en détail sur la description de la source de chaleur injectée dans l'équation d'énergie qui est proposée comme une fonction d'évolution de la réaction chimique. La variable utilisée pour décrire l'état de durcissement est calculée comme alpha. Nous passons donc d'une résolution standard d'éléments finis à une résolution de degré de conversion de température couplée (Tα inconnus). Le premier test est un échantillon de cylindre bloqué où la température de la surface augmente de 3 degres par minute (réaction lente). L'objectif principal était de faire évoluer la température de l'échantillon contrôlant le gradient d'épaisseur. Après environ 1500 secondes, le matériau réagit pour SMC-UL et après 1800 secondes pour SMC-HP. La prédiction numérique du même test prouve que le modèle nous permet d'adapter les données expérimentales basées sur l'évolution de la température dans le noyau de l'échantillon ainsi que l'évolution de l'épaisseur pendant les étapes de réaction au réchauffement et de refroidissement. En complément de cette étude, un cas de sensibilité des propriétés caractérisées telles que l'énergie chimique, la conductivité, la capacité calorifique, la résistance au contact thermique a été réalisé en concluant que les variables qui modifient le temps de réaction et le pic de température sont la résistance au contact thermique et l'énergie chimique De la pâte. Pour cette rai-son, il est important de bien caractériser ces propriétés pour d'autres matériaux.

Un deuxième cas de test a été étudié où le matériau se comprime pendant le chauffage. Cette affaire consistait en un modèle d'expérience en moulage par compression réelle où la paroi chaude d'un moule comprime l'échantillon jusqu'à atteindre une certaine épaisseur. Les moules se maintiennent jusqu'à ce que la réaction chimique ait lieu et que le matériau soit complètement consolidé. Les données expérimentales ont enregistré une force pendant le test, la température au centre de l'échantillon et dans les positions latérales ainsi que le capteur de flux de chaleur placé sur les surfaces du moule et du poinçon. La comparaison des données expérimentales avec les prédictions des simulations numériques est en général une bonne estimation de ce travail dans des conditions industrielles réelles. L'évolution de la température pour SMC-UL est en accord avec les données expérimentales. La température de la surface centrale et les flux de chaleur calculés numériquement prouve que le modèle et l'ensemble de simulation décrivent assez bien le test expérimental. Cependant, la force pendant la compression n'est relativement pas bien prédite en fonction de la raison principale pour laquelle le frottement contre le moule n'est ici pas précisément modélisé. La même procédure est répétée pour le SMC-HP, mais la réaction chimique décrivant numériquement est moins conforme à l'expérience. Enfin, à la fin du chapitre, la prédiction du modèle thermo-cinétique pour une géométrie industrielle est donnée montrant que pour une pièce utilisée dans Plastic Omnium Auto Exterior, le temps de réaction sous l'état de la température du moule autour de 150Celsius est d'environ 45 secondes.

Introduction

Consolidation mechanism modifies the viscous nature of the polymer arrangement, by connecting the whole molecule chains, hardening the macroscopic behavior of polymer compound. This process transforms a viscous composite body into a rigid part. The materials studied in this work, once consolidated are used as lighted structural pieces and therefore it is important to model its behavior from the viscous state to the rigid state. The resin in thermoset preform of Sheet Molding Compound (SMC) has the chemical component which at hot temperatures reacts and consolidates the product. The final composite is then a rigid part, containing mineral fillers and long chopped fibers of 25-50 mm length [1]. Due to its rigid mechanical behavior once consolidated, those products are used to manufacture semi-structural parts such as tailgates, front end carriers, bumpers, inside hatchbacks and trunk floor or well other applications mentioned in Chapter 3. During compression molding of composites parts, specifically those used in the automotive industry, two main features have to be analyzed: firstly, the deformations step for which a proper mechanical modeling describing the anisotropic behavior of the SMC was presented in the previous section; secondly, the study of the cross-linking reaction and its influence on its thermal/mechanical properties which is the main subject of this chapter, . Many authors have focused on studying the mechanical behavior during the mold filling and some few have proposed or characterized a thermochemical model adapted to SMC. Barone et al [2] and Kotsikos et al [3] started describing the SMC as viscosplastic deformable body before curing. Le Corre, however, [4] improved the modeling by observing the anisotropic behavior on the mechanical response and by adding the thermal influence on their viscosity. Later, Dumont et al in [5] proposed a model that accounts for the anisotropic behavior exhibited by reinforced materials by means of structural tensors. In [6] a study on the molding process for parts molded at room temperature was provided. However, models taking into account the thermo-kinetic nature of such composites remains a complementary key towards the thermo-mechanical modeling of reinforced parts. First studies were accomplished by Lee in [7] and Maazouza [8] using a model based on the radical polymerization mechanisms. For structural parts, where fiber concentration increases up to 38% in volume, the mechanical model modified by the fiber concentration was study in Chapter 3 and in the following, a thermo-kinetical model coupled with the mechanical response is introduced based on the collaborative work performed in [9].

Within the framework of producing structural parts, updated formulations with stronger curing mechanism are found. However, those formulations lead to different molding behavior with respect to standard SMC products [5]. During industrial SMC compression, a strong thermo-mechanical coupling between thermal and chemical mechanisms defines the final geometry of the molded part. Specifically heating, flow, cross-linked reaction and cooling are active conditions faced during the compression cycle. The increase in fiber concentration within the material decreases the thermal conductivity of the compound, translated into a very local exothermic reaction, inducing temperature and conversion degree gradients through the thickness [10] [11] [12]. Coupling those physics with thermal expansion and chemical shrinkage enhances the modeling of the final geometry and improve the estimation of the residual stresses [13] [14] [15]. This reaction increases the viscosity of the part since the cure of thermoset resins converts liquid monomers into three dimensional networks [16]. As part of this project, the high concentrated SMC materials were characterized [17] and a consistent model is here presented. In this study, we observe not only the thermal evolution of the sample, but we also include resin reaction mechanism. Both, thermal and kinetic evolution are coupled to the mechanical problem by means of the variation of dense density. A first introduction of the thermo-mechanical model was presented in [18] and here we extend the approach to fully coupled simulations.

In this chapter, firstly the fully thermo-kinetical model is introduced following the work in [17]. The form of the cross-linking reaction using Bailleul model is explained, as well as, its connection to all thermal and kinetical properties. The evolution of the density with respect to thermal expansion and chemical shrinkage is presented complementing the full description provided in Chapter 1. Moreover, this model can also deal with the anisotropic volume variation induced by the fiber network. Secondly, we study the temperature evolution in the sample under reaction and thickness evolution in a static PVT-α experiment and comparisons between computations and experimental measurements are made. After validating the model, compression molding simulation under plane strain condition for isotherm setups at 50 o C and 80 o C are studied. The aims is to analyze the viscosity evolution with temperature without reaction and characterized the friction of the SMC sample with the mold. Then, plane strain compression is performed on real industrial conditions at mold temperature of 150 o C while the SMC sample has an initial temperature of 30 o C. The SMC sample is deformed until filling the mold cavity. The sample is hold until reaction is triggered and the temperature, reaction degree and thickness evolution during the molding cycle are compared with experimental data coming from [17].

Thermo-kinetical modeling of SMC

Let us consider a sample of fiber reinforced composite, such as SMC, described by an orthotropic non-Newtonian compressible constitutive law and modeled by the Compressible-Stokes equation presented in Chapter 3. For non-isothermal computations, we extend the modeling by adding the thermal evolution of sample, describing its thermal and kinetical evolution as active variable (T, α). Including a model for the evolution of the conversion degree α, the governing equations for the thermal evolution is given by the compressible heat transfer equation:

ρ c p dT dt -ρ∆H α dα dt -T χ v t dp dt = ∇ (k∇T ) + s : ∇v (4.1) dα dt = F (τ, α, T ) (4.2)
with T the temperature, α the reaction degree, ρ the density of the mixture, s the deviatoric stress tensor, v the velocity, p the pressure and χ v t the volumic thermal expansion coefficient of the mixture. The thermal parameters c p states for the heat capacity of the mixture, κ the conductivity tensor, ∆H α the chemical energy in the composites and F (τ, α, T ) the evolution function of the conversion degree dα/dt. The variable τ states for the spare time that requires the system before reacting, known as induction period.

Expressions for the evolution of thermal properties

The thermal properties of the SMC materials under study depend on the temperature T . Likewise, accounting for the polymerization process undergone by the paste, its properties depends on the conversion degree α. According to the characterization given by the Laboratory of thermo-kinetics in Nantes (Table 4.1) [9], the expressions for the heat capacity c p and the transverse conductivity κ stand:

c p (α, T ) = X f c p,f (T ) + (1 -X f ) [α c p,cured (T ) + (1 -α) c p,uncured (T )] (4.3)
The heat capacity responds to a mixture of the fiber with the paste. Firstly, it depends on the massive fiber fraction concentration X f and its heat capacity c p,f . Then, the contribution of the paste before reaction c p,uncured and in cured state c p,cured .

For the conductivity, the characterization is provided directly with the SMC solving an inverse heat conduction problem [19]. There, the material is studied as a homogeneous mixture. The properties of the axial conductivity (normal to the fiber network) reads:

κ (α, T ) = α κ cured (T ) + (1 -α) κ uncured (T ) (4.4)
The planar orientation of the fiber network may induces an anisotropic behavior of the thermal evolution in SMC materials. Thus, opening up the discussion of considering as well a tensor of conductivity. A further discussion followed by a model to take into account this anisotropic behavior is carried out later in this chapter (see section 4.1.5). Another property evolving during the thermo-kinetical cycle is the density. Since the materials undergoes thermal and chemical changes, thermal dilatation and chemical shrinkage takes place during the molding of the part. Thus, modifying the geometrical aspect of the piece, caused by density variations. Such evolution responds to a continuous variation of the volume ruled by the mass conservation equation. This variation during the thermal and chemical evolution is studied in section 4.1.4. In the following, we introduce the kinetical model employed to describe the curing mechanism of SMC profiles.

Kinetic Evolution of SMC -Bailleul's Model

Cross linking reaction is usually described by a macro-scale variable denoted as reaction degree. Despite polymerization includes many chemical reactions, the empirical models stand the hypothesis of describing the curing behavior by a unique local reaction. Several kinetic models has been reported for thermoset resins and were associated to composite materials. Those models are mainly distinguished as phenomenological [20], [21], [22] and mechanistic [23] [24]. According to the work presented in [17,9], the cure of SMC high profiles parts are modeled by the empirical kinetic equation developed by Bailleul [20]. In such model, the induction period is considered. This corresponds to a spare time while the inhibitors presented in the reactive system are consumed, preventing the initiation of the reaction by neutralizing free radicals. This ignition time is governed by its own kinetic and by the imposed temperature. However, the model do not consider the vitrification phenomenon (or diffusion effect). If the curing temperature occurs in the vicinity of the glass transition temperature, a slowing down of the reaction takes place. This is caused by the decrease of mobility of the polymer chains. Such effect might lead to partial conversion and must be described by a gradual diffusion effect on reaction kinetic during the cure [25]. Fortunately, this effect is more probable to occur for low curing temperatures compare to imposed mold temperature in compression molding process. The kinetic model used to describe the curing of SMC materials stands:

dα dt = W (h(t)) K(T ) G(α), (4.5)
The computation of the reaction rate evolution F (α, T ) eq. (4.2) is usually described by the product of a temperature function K(T ) eq. (4.6), a function depending on the conversion degree G(α) eq. (4.7) and by a boolean function depending on the induction time:

K(T ) = k ref exp -A T ref T -1 (4.6) G(α) = a i α i . ( 4.7)
The induction period function h(T ) (eq. (4.8)) reads:

h(T, t) = t ref -τ (T, t). ( 4.8)
The two contributions of the kinetic model are plotted in Figure 4.2. We notice that the shape of the reaction is supported by the evolution of function G(α). Also, at higher temperatures the reaction is stronger, notice the increment of K(T ) respect to T . The activation function W (h(t)) equals to zero if the ignition period has not been reached yet (h(T, t) > 0) and to 1 once all the inhibitors of the resin have been consumed (h(T, t) ≤ 0). The function τ (T, t) evolves according to

dτ dt = exp -B T ref T -1 (4.9)
The ignition time τ can also be expressed in its differential form, in order to be used in a convection scheme during the compression molding simulations and reads:

∂τ ∂t + v • ∇τ = exp -B T ref T -1 (4.10)
In general, all thermal parameters depend not only on temperature T but also on the conversion degree α. From a numerical perspective, we focus on modeling the the thermo-rhelogical-kinetic coupling of fiber reinforced materials studied by R. Cardinaud [9] Under an isotherm condition, the time to trigger the reaction is obtained by vanishing eq. (4.8) in which the ignition time follows eq. (4.11). That gives a reaction time t on

t on = t ref e -B(T ref /T -1) (4.11)
which plotted in Figure 4.3. From this expression, important information can be obtained: the higher the temperature the faster the reaction is triggered; the high performance profile (HP) is more reactive than the lighter profile (UL). Notice that when a point of the material is at temperatures higher than 150 o C (423 K), the material reacts in less than 1 second. However, the energy released by the UL profile is higher causing that in real conditions the maximum temperature in the system is found for this profile. Under real conditions, the induction time can not be estimated easily since it depends on the thermal history. For compression molding process, the temperature in the sample evolves from 25o The mold hot temperature imposition will then define the thermal history according to the conductivity of the material. The temperature on the surface of the SMC is also affected directly by the thermal contact resistance. Thus, it modifies the determination of the beginning of reaction. In addition, this provides a first notion of thermal gradients within the thickness. Due to this thermal gradient, then it is natural to find that the reaction starts in the surface earlier than in the core. Whereas the mechanical gradients, namely, pressure and stresses, occur along the piece, the thermal evolution is more predominant within the thickness. The justification of fully 3D simulations is then hereby justified.

Kinetic and thermal variations cause local changes in the material density which is translated into an induced velocity. Hereafter, we model the effect of thermal and kinetic variations on the strain rate tensor. The coupling with mechanical problem is introduced and focused on a 3D examples.

Expansion and Shrinkage of SMC

Let consider the thermal strains produced by temperature changes during the heat conduction and the chemical strains during the crosslinking reaction. These strains have inherently a volumic nature (thermal expansion or chemical contraction) and do not cause any shear, placing their displacements only in the main directions of the strain rate tensor. The total strain rates ( ε) are then due to the (additive) contribution of the mechanical strain rate (here considered as viscoplastic) (D), i.e., those produced by the stresses, the thermal strain rate ( εth ) and the chemical strain rate ( εch ):

ε = D + εth + εch . (4.12)
Finally, the anisotropic volumic variations due to thermal and kinetic evolution are modeled by the tensor of dilation and shrinkage. For example, the thermal strain rate is a tensor computed using the dilatation factor χ t and the temporal temperature variation:

εth = χ t dT dt (4.13)
where χ t states, in the case of SMC materials, for the tensor of anisotropic dilatation coefficients defined by

χ t =    χ x t 0 0 0 χ y t 0 0 0 χ z t    (4.14)
where usually χ x,y t < χ z t if the vertical axis is the plane perpendicular to fiber network. Therefore, a thermal variation of dT /dt induces different displacements in the plane of the fiber and in the thickness (for isotropic materials that displacements would be the same).

In the Eulerian framework, we express the temporal derivative by taking in account of the convective part which gives rise to

εth = χ t ∂T ∂t + χ t v • ∇T (4.15)
The trace of this tensor defines a volumic variation which can be translated as a density variation. For example by taking the trace of thermal dilatation tensor, we can write:

tr( εth ) = (χ x t + χ y t + χ z t ) ∂T ∂t + (χ x t + χ y t + χ z t ) v • ∇T (4.16)
by defining the volumetric dilatation coefficient χ v t = χ x t + χ y t + χ z t and by noticing that the second terms can be grouped, we write:

tr( εth ) = χ v t ∂T ∂t + χ v t v • ∇T (4.17)
The same reasoning for the chemical shrinkage gives:

εch = χ α ∂α ∂t + χ α v • ∇α (4.18)
and its trace stands:

tr( εch ) = χ v α ∂α ∂t + χ v α v • ∇α (4.19)
.

By taking tensors defined by dilatation and shrinkage coefficients, we can study the volumetric variation as a directional depending phenomena. For instance, we encounter lower dilatation and contraction terms in the fiber plane. The fiber network induces a mechanism of braking the natural dilatation or contraction of the matrix in our reinforced composite. It is easier to move along the normal direction and more difficult towards the fiber network alignment. From thermo-mechanical point of view, we found in here another source of anisotropy.

Thermal and chemical dilation in any reference system

By knowing the normal to the fiber plane in any point of our material, we can define the anisotropic dilatation and shrinkage tensors. Let be M = n ⊗ n the structural tensor built with the vector normal n to the fiber plane. The tensors are built automatically by using the following expressions:

χ t = χ plane t (I -M) + χ normal t M (4.20) χ α = χ plane α (I -M) + χ normal α M (4.21)
SMC materials under study have different compressible properties under heating than cooling. A more general approach should include thermal expansion coefficients that depends on the sign of the cooling rate. In the scope of this work, we only study the compressibility of SMC under heating.

By assuming n = e 3 the tensor of thermal dilatation coefficient χ t and the tensor of chemical shrinkage coefficients χ α states:

χ t =    χ plane t 0 0 0 χ plane t 0 0 0 χ normal t    (4.22) χ α =    χ plane α 0 0 0 χ plane α 0 0 0 χ normal α    (4.23)
For UL and HP SMCs, the values are presented in Table 4.3. There, we noticed the high ratio (≈ 20) between the normal and transverse coefficients. The strain rate induced in the normal direction are by far the ones inducing geometrical modifications of the piece. During heating or transformation (from viscous state towards consolidated state) the volumic variations induces variation of material density. In the following, we link such volumic variations to the density evolution of the part.

Evolution of the density due to thermal and kinetical evolution

According to the multiplicative decomposition introduced in Chapter density evolution is presented in Chapter 3 and in this section, we focus on the evolution of the dense density mainly due to the thermal and kinetic effects. Consequently, the dense density is a function of the temperature and the degree of reaction and is ruled by the differential equation: Notice that for isotropic materials, if the volumetric thermal dilation coefficient χ v t and the volumetric chemical shrinkage χ v α are constant, the differential equation has an exact solution of the form:

ρ d = ρ 0 e χ v t (T -T ref )+χ v α α (4.25)
where ρ 0 corresponds to a reference density at T ref and in raw state (α = 0) (some values are given in Table 4.4). In our simulations, we solve the differential equation taking into account both convective and local variation of the dense density as:

1 ρ d ∂ρ d ∂t + v • ∇ρ d = χ v t ∂T ∂t + v • χ v t ∇T tr( εth ) + χ v α ∂α ∂t + v • χ v α ∇α tr( εch ) (4.26)
which enhances the formulation given in Chapter 1. The anisotropy on the thermal/chemical volume changes is also present for other thermal properties such as the conductivity κ.

The transverse conductivity faces the fiber network while the planar conductivity acts all along the fiber lengths. That point is analyzed in the next section.

Anisotropic Conductivity in thermal resolution

Due to the fiber network, the thermal properties of the SMC are anisotropic. That intervenes for the conductivity as well as the thermal dilatation and chemical contraction coefficients. In particular for SMC, the properties can differ between its planar and normal direction. In order to take into account such anisotropy, the thermal solver should be as well enriched to account for anisotropic conductivity inputs κ. The conductivity in any framework is expressed as the rotation matrix R of the main system into the reference framework and the reference conductivity tensor κ 123 :

κ = R t κ 123 R (4.27)
and the reference conductivity tensor is built using the structural construction tensor M introduced earlier to model the anisotropic mechanical behavior of SMC. There, two inputs are required: the conductivity in the normal κ n and the planar κ t directions:

κ 123 = κ t (I -M 123 ) + κ n M 123 (4.28)
After having introduced all the details of the thermo-kinetical modeling of SMC profile, we noticed the strong links between variables and the statement of anisotropic thermal behavior inducing as well an anisotropic mechanics. The heat equation and the kinetic models are solved and the density is updated. To account for inhibition time an extra differential equation needs to be solved. In section 4.1.7 a summary of all equations to be solved is given as well as the link to the expression require for the fully thermo-kinetical resolution. At the same time, such summary helps to condense the equations and work as reference for the current work.

Mechanical coupling -Complex viscosity

A complex viscosity evolution can be noticed from experimental test. According to compression test performed to the past of the high profile SMC (HP), the Modified Castro & Macosko model is used to model the crosslinking effect on the viscosity.

η = η 0 (T )   α gel α gel -α A+Bα -1   (4.29)
where η 0 (T = 118 o C) = 7000P a.s, A = 11, B = 1 and α gel = 0.56. These parameter values were identified from oscillatory tests done at isothermal temperature (around 118 o C or 391 K ) and for a fixed frequency of 0.5 Hz. This formulation complements the thermal Arrhenius law usually applied to SMC. This modified Castro & Macosko model was proposed in [9] for the paste of HP SMC. We did not include this viscosity evolution in our modeling, since the material once reacted is in static molding conditions.

The mechanical transition between viscous SMC stresses until consolidated solid mechanics is a requirement for the computation of residual stresses. Connecting a viscous approach to a solid motion is considered as a further topic to be discussed for future work.

Thermal resolution for SMC materials

For the resolution of SMC materials, the heat equation is solved taking into account the non-constant values of the thermal variables and it is coupled to the kinetic equation by solving:

ρ r ρ d c p dT dt -∇ (k∇T ) = ρ r ρ d ∆H α dα dt + T χ T dp dt + +s : ∇v (4.30) dα dt = F (τ, α, T ) (4.31)
The dense density evolution is ruled by the expression:

1 ρ d dρ d dt = χ v t dT dt + χ v α dα dt (4.32)
with χ v t and χ v α being respectively the thermal expansion and chemical shrinkage coefficients. In order to take into account the induction period before the beginning of reaction and to determine this starting time, a convection scheme of the ignition time variable τ is solved:

∂τ ∂t + v • ∇τ = exp -B T ref T -1 (4.33)
This ignition time is compared to the reference time depending on the material and will define the beginning of reaction. For the whole thermo-kinetical problem, we solve four differential equations coupled by a point fixed method.

In the next section, we perform numerical computations of experimental apparatus with the numerical model proposed in this chapter. There, we look at the numerical predictions on temperature evolution, detection of reaction, heat fluxes and others to verify and validate our numerical approach and the thermal characterization of the SMC under study. to its relative low thermal conductivity. The apparent chemical shrinkage has an amplitude of 4.3% and the final volumetric shrinkage is around 2.5%. However the peak of temperature seems underestimated by the numerical computations. The reaction degree evolution during the thermal cycle is also presented.

In section 4.5 the velocity of the part induced by the thermal gradient is depicted. We notice the velocity profile along the sample due to density variations. The SMC level-set is convected with the velocity profile determining implicitly the material new surface. Due to the thermal coupling, a thickness evolution is observed. Additionally, the ignition time resolution enables the prediction of the beginnings of reaction. For UL SMC, the reaction begins in places where τ > 325 s. This value depends on the thermal history. The thickness evolution and the density variations are plotted in Figure 4.10. The experimental measurements and the numerical predictions are in agreement.

Study of HP SMC

Similarly, in Figure 4.11, the temperature, the reaction degree, the thickness evolution and the density during PVT-α apparatus is presented. The HP SMC material induces a weaker reaction (observed by a lower temperature peak compared to UL SMC). The reaction takes longer time and that means that the reaction mechanism is slower than for UL SMC. Thickness and density evolution are shown in Figure 4.12. A sensibility analysis has been performed in order to study the thermal parameters which have the most influence on the the reaction. A set of height simulations were performed by modifying the energy released by the reaction, the conductivity of the mixture, the heat capacity, and the thermal contact resistance by 10%. In this manner, we evaluated their influences on the temperature peak and the beginning of the reaction. In Appendix: E the effect of each parameter on the time evolution of the temperature and the reaction for the PVT-α test is given. For sake of simplicity, we present in Figure 4.13 a summary. We notice the maximum temperature find in the center of the sample along the reaction and the time variation for the beginning of the reaction, with respect to the values given by the first computations earlier in this section. For example, a positive value gives the amount of second that it is necessary to wait before the beginning of the reaction. This value is taken by comparing the time difference when the core of the sample has a reaction degree α = 0.5. We notice that thermal contact resistance (RTC) and the reaction energy are parameters which have the most influence on the reaction modeling than conductivity and heat capacity. They delay the beginning of the reaction and also control the peak in the temperature found in the sample more importantly. This sensibility analysis gives us information about the accuracy needed to characterize the materials. During the simulation of this test, the thermal contact resistance RT C is assumed constant and chosen equal to 5 10 -4 m 2 K W -1 . During the simulation the value of RT C is assumed constant. The thermo-kinetic numerical scheme used also impacts the evolution of temperature and reaction degree. At each time step, we used a fixed point scheme in order to solve all the nonlinearities. The thermo-kinetic problem is a coupled system and the non-updating of properties and variables along the iterations may results in incorrect solutions. In addition the use of adaptive time step improves the efficiency of the numerical resolution by detecting the time of the reaction beginning and describing it properly. The conversion degree value α is very sensible to the numerical scheme. If the temperature-reaction-property coupling is not well satisfied, non-realistic values such as α > 1 can be obtained. We have checked that this fix point scheme is necessary to obtain the good behavior of the variables. In the next section, we present the case of the SMC compression within a channel (plane compression) which corresponds to apparatus built in ??. These computations are done at iso-temperatures of 50 o C (323 K) and 80 o C (353 K) for which there is no reaction. The force applying on the punch to fulfill the mold cavity is compared to experimental measurements. We evaluate the viscosity of the material using the Arrhenius law and compare the slip and no-slip case during the compression of the sample for both material (UL and HP).

Isothermal Cases -Plane Compression

The plane compression test consists on compressing the material while conducting motion along a channel. The SMC is placed in the center of a mold and compressed while the geometry of the channel enables only deformation in one direction, namely, longitudinal direction. In Figure 4.14, it is shown an example of plane compression test. It is seen how the material deforms along a plane. This test has a plane state of deformation along ε 11 and ε 33 , while stresses are mainly in σ 22 and σ 33 direction.

The tested materials are the UL SMC and HP SMC supplied by Plastic Omnium (Saint Julie, France). The produced sheets are made of a polyester based paste matrix containing 24% of glass fiber bundles for the UL SMC and 38% of fiber in volume for the HP SMC. Compression molding experiments were performed on a laboratory hydraulic press (maximum axial force of 9 kN . The press was equipped with a standard rectangular mold able to produce thin (≈ 2 mm) plates with a in-plane (e 1 , e 2 ) surface of 140 × 40 mm 2 . Temperature, pressure and heat flux sensors were located in the upper and

Experimental procedure -case 1

The first experimental procedure was to set plane compression without reaction.

Isothermal compression test were performed at temperatures of 50 o C and 80 o C. The main point was to correlate the viscosity thermal translation given by the Arrhenius expression in 3 with respect to temperature and to compare the experimental force during the test for two boundary conditions: slip and no-slip case. In such manner,it would be possible to get information on the friction between the mold and the SMC.

The mechanical modeling of the three bodies in our numerical approach is given by assigning a high viscosity to the mold and the punch as explained previously in 3. Here, η m = η p = 10 3 η SM C eq . In such ways, the mold and the punch can be considered rigid by comparing to the SMC. The boundary condition for the velocity is only imposed at the boundary of computational domain to avoid numerical singularities. Therefore, the border of the mold is set to v = 0, while it is v x = v y = 0 and v z = v imposed for the punch border. In the experimental test, the initial thickness of the UL SMC sample is 7.3 mm whereas it is 6 mm for the HP SMC. The velocity imposed on the punch v imposed is nonconstant. As explained previously, the test condition attempts to imposed constant strain deformation rate v z /h. However, once the maximum force of the machine is reached the velocity is imposed to keep that a constant force of 9 kN . This implying a change on the slope in the curve given the thickness evolution during the test. In Figure 4.16, the thickness evolution along the test during the isothermal experiments are presented for both materials and both iso-thermal conditions. Please notice that the compression test start somewhere around 9 to 11 seconds. This is due to the spare time between the initialization of the recording data compare to the time the punch start to move and touches the SMC.

In our numerical simulation, we imposed the velocity of the punch to follow up the experimental condition. In that manner, we ensure the same kinematics conditions during the test (velocity and thickness evolution), while comparing the force prediction in our numerical simulations.

In Figure 4.17 and Figure 4.18, the force prediction along the test is compared to the experimental recorded data. For the time of comparison of this plot, the friction of the SMC sample against the mold has not been experimentally measured. For that reason, we proceed to compare the force to the limit cases: slip and no-slip conditions.

At 50 o C, we notice that for the slip and no-slip cases the departure of the force in the experimental data is closer to the slip case. However, during the compression, the experimental force gets close to force computed with the non-slip condition. Before reaching the final thickness and the end of test, the force is located between the slip and no-slip case curves. For isothermal cases the evolution of the viscosity is only related to the kinematics aspects. For the HP case, we notice that the curve of the experimental force is located mostly between the slip and no-slip curves.

At 80 o C, the experimental force for the UL case is closer tot the slip case at the beginning of compression. After 3 seconds, the force is located between the slip and no-slip cases meaning a friction at the SMC/mold interface. The temperature at 80 o C makes the SMC less viscous and a lower compression force induces the SMC motion. Therefore the kinematics is more associated to slip cases. The HP curve indicates that the experimental data and the slip case are basically the same during the expelling of porosity (10-10.5s). After materials becomes incompressible, the force is placed in between the non-slip and slip cases meaning a friction during the test.

These comparison proves the evolution of the friction nature with respect to temperature and flow kinematics. The response of the experimental force compared to slip and no-slip cases gives us information on the friction temperature dependency of SMC against the mold and flow motion (or stress) dependency. A campaign to characterize the SMC friction should be carried out. However, noticing the kinematic conditions of the friction according to [27] the plane compression test is not the best test to characterize such condition. They proposed a circular rheometer, where pressure sensors are located along the radial direction. Thus, along the compression of the sample, the difference on the pressure signal along the radial direction of the sample enables by analytic continuum mechanics equations to obtain an approximation of the friction coefficient without assuming any a priori model. This work remains for perspectives and is encouraged for future work.

Compression, heating and reaction of SMC

In this section, we use the same hydraulic press described in Figure 4.15. The compression strategy remains the same, however in order to take into account the reaction of the material and to simulate the real compression molding process non-isothermal conditions are studied. The mold temperature is set to an usual industrial value and the SMC sample to room temperature. During the compression, the thermal exchange between mold and SMC heats up the part triggering the reaction. We follow the reaction thanks to heat flux sensors and the measurements of temperature along the mold surfaces.

Experimental procedure -case 2

The tests were carried out setting a constant temperature in the mold and punch located at 50 mm from the surface respectively (see Figure 4.15). The imposed temperature is kept constant to 150 o C (423 K) and a vertical kinematic is chosen for the punch. It corresponds to a deformation velocity v/h of 0.3s -1 . The experimental procedure has three main steps. First, the sample of initial dimensions 60 mm × 40 mm × h 0 at room temperature is set in the middle of the hot mold, at that point the material enters in contact with the lower and lateral sides of mold whereas the upper side lefts free. A spare time t r of ≈ 15 s is intentionally imposed in order to simulate the charge positioning in real compression molding. Therefore there are thermal transfers between the mold and the SMC before compression begins. The punch goes down at constant speed until it enters in contact with the SMC sample. After this contact, the punch speed is piloted by keeping the strain rate constant until the maximum force is reached. Then, the SMC deformation is made by applying a constant force. When the material fulfills the mold cavity, no longer vertical deformation is detected. The material remains kept into the mold until the chemical reaction is triggered. The sensor of flux and temperature detect the end of the reaction and the sample is released from the mold. The final thickness h f of the sample is measured. These values as well the experimental conditions are summarized in Table 4.5.

Experimental conditions Material

Modeling of thermal evolution in the mold and punch

Since the temperature imposition is set at 50 mm from the SMC surfaces, the thermal evolution inside the mold and punch is necessary to have an accurate description of thermal transfer. The thermal properties of the mold and punch -made with 40CMD8 steel -are provided in Table 4.6.

Thermal properties of Mold and punch

Property value units Once described the modeling of the mold assigned by its conductivity, the thermal contact resistance during the test is not assumed constant as explained in the sequel.

c p 475 J kg -1 K -1 ρ 7850 kg m -3 κ 36 W m -1 K -1

Thermal contact resistance (RTC) in SMC materials

Concerning the RTC values during the compression test, we take into considerations two situations labeled as RTC-static and RTC-dynamic: the static case during the first 15 s where the SMC sample is deposed into the mold and the dynamic case where the SMC is compressed by the punch. It is natural to understand that the RTC values are not the same in these two cases. For that reason, we model a different RTC(-static) when the contact between the SMC and the mold is done at ambient pressure and a RTC(dynamic) when the punch squeezes the SMC part. Here, we use RTC = 2 10 -3 m 2 .K.W -1 for the static case and RTC = 5 10 -4 m 2 .K.W -1 once the punch touches the part. The change from the static to dynamic case is performed immediately once t = t r . A linear evolution between this two limit values to take into account the evolution of the thermal contact between mold walls and SMC surfaces during the compression step escapes from the scope of this work. It can be described by the evolution of the deformation (or as a function of the thickness reduction) according to the observations in [9].

The experimental data recorded on the sensor located according to Figure 4.20 and Figure 4.21 measure a surface temperature evolution in the mold and punch during the thermal cycle including the local increments due to the chemical reaction. In Figure 4.19 the velocity field and the temperature along a partial section during the compression of the piece is presented. We can notice as well the mesh during the simulation.

Analysis in the Mold and Punch Surfaces

In Figure 4.20 the temperature evolution on the mold surface is compared to the recorded measured in the sensor T2. The temperature in the sensor begins at about 145 o C (418 K) with respect to the imposed value 150 o C (423 K). The resistance of the steel mold imposes 5 degree less in the surface of the SMC due to 50 mm distance between the thermal regulation and the SMC surface. The contact between the cold SMC and the mold surfaces decreases the temperature in the sensor until 138 o C (411 K). The simulation is able to detect such mold cooling since the thermal equation is solved on the Figure 4.20b shows that the numerical evolution of the sample thickness meets the measured data. This is obtained by imposing a velocity directly related to the punch velocity used during the experimental procedure. Finally the conduction in the piece defines the new unsteady state. According to this comparison, the numerical model predicts the temperature evolution correctly. In Figure 4.22, the heat flux measured in the center (T2) and in the lateral position (T1) are compared. The continuous line stands for the experimental data whereas the dashed lines for the numerical predictions Figure 4.22a. The heat flux measured in position T 2 detects the contact of the mold with the sample at ≈ 15 s resulting in a negative heat flux. It means that the heat goes from the mold towards the part. The experimental heat flux in the lateral position T 1 shows the cross-linking process. The reaction moves from the hottest part of the sample towards the center. The heat flux is detected to be higher in the sensor located at position T 2. The shapes of numerical and experimental curves are in relative good agreement. However, the numerical simulation predicts the beginning of the reaction with a little delay. For this matter, according to the sensibility analysis perform in the case of PVT-α test, the value of thermal resistance during the plane compression must be the most determinant factor to fit the experimental data with the numerical prediction.

The total heat flux measured in the center is obtained by the addition of the flux measured in the two sensors T2 and T3. The numerical and experimental total heat fluxes are then compared. We notice a nice agreement between both results. We quantify the delay of reaction for the numerical computations and it is around 4 s for UL SMC. In Figure 4.22b the thickness evolution measured by the LVDT sensor shows a dilatation zone due to the thermal expansion followed by a shrinkage due to the chemical reaction.

In Figure 4.21 the temperature for the HP SMC is compared. Similar to the case of UL SMC, the modeling of the mold and SMC enables the visualization of the mold cooling when contacts is set (after ≈ 15 s) Figure 4.21a. For this case the mold temperature does not present a peak since this material is less reactive because of the presence of more fiber in its volume. The variations of experimental and numerical thicknesses are compared in Figure 4.21b noticing the good agreement which comes from the suitable boundary conditions imposed on the punch. The heat flux presented in Figure 4.23 shows the numerical prediction to anticipate the reaction Figure 4.23a. Additionally the shape on the heat flux curve is narrow for the numerical case. The wider experimental curve seems to measured a smoother reaction. The total heat flux is then compared The measured data recorded during the experimental campaign was compared to the numerical prediction. Measured data at the surface level of the mold and punch. Those data provide information on the thermal evolution outside the SMC. To get information inside the SMC, the numerical data are presented in the next section. The numerical simulations allow us to see the evolution of the cross-linking reaction inside the material.

Analysis of thermal evolution inside the SMC

After comparing the mold and punch experimental data with computational results, we can analyze information of the thermal evolution given by the numerical simulation.

The temperature in the SMC and the reaction evolution during the compression process is presented for the UL profile in Figure 4.24 and for the HP profile in surface in blue. The compression molding process hereby simulated has three main stages as mentioned before. During the static step, where the piece is placed in the mold and is heated up from the bottom, the surface temperature at the depth of 0.5 mm is about 360 K whereas the core departs at 300 K. During the first 15 s, we notice a conduc- tion from the mold surface to the core of the sample. When the second step begins and the compression of the part is carried out, we notice how the conduction from the mold surface to the core of the part occurs faster. Basically, the reduction of the thickness makes the heat to travel faster towards the core having as result an increment on the temperature rate. During the compression step, we notice that the sensor located at the surface captures the beginning of the reaction. At about 45 s, the reaction occurs in the core as notice in Figure 4.24b. The reaction begins at the surface of material but finish first in the core. Figure 4.24b shows the sharp behavior of the reaction in the core of sample. The strong gradient of the reaction creates numerical instabilities that were solved by the use of adaptive time method. The mesh is also adapted with respect to the conversion degree gradient. These two last improvements are not discussed further in this section, however there are directly related to our numerical strategy used to deal with the multiphysics thermo-mechanical problem.

For the high profile SMC, the temperature peak is almost not detected at the surface. The core temperature in Figure 4.25a does not present a high difference with respect to the mold temperature. The reaction evolution in Figure 4.25b proves that for this material the reaction in the core and the surface is less differentiate than in the case for UL SMC. That enables a more homogeneous reaction along the thickness. The high profile SMC seems to have a less sharp thermo-kinetical behavior than lighter profile.

Evolution of velocity profile during plane compression for high performance SMC

The velocity profile evolution during the compression molding process taking into account all the physics involved is summarized in Figure 4.26. At first, the thermal evolution induces the thermal dilatation assigning a velocity profile as pictured in Figure 4.26(a). Notice that the mold and punch velocities are zero. This velocity shows the thermomechanical coupling of the mass equation which modifies the density of the material. In

Heat transfer and reaction in industrial like geometries

In order to prove the the efficiency of our model for complex geometries, we study the thermo-kinetical evolution on a mold which gathers all the industrial features Figure 4.27.

In this section, only the thermal part of the process is computed. The mechanical deformation of the part is not presented in this section, however it was introduced and analyzed in Chapter 3 and a summary of the deformation mechanism is recalled in Figure 4.28. We focus this section on the static heating and the chemical reaction after the mold filling. The SMC thickness is of 2 mm and the mold and punch are supposed to be in the static position touching the preform part as depicted in ??. The mold cavity is then fulfilled and the punch is at rest. The initial conditions and thermal boundary conditions for the thermal cycle are summarized in Table 4.7. For accomplishing this simulation, the use of time adaptation is required in order to describe properly the evolution of the cross-linking reaction. The time step found during the reaction reaches 6 ms. Indeed the heating of the part is done in 720 increments (up to 25 s) and the reaction from [25 s, 45 s] is performed in 1650 increments. In total for modeling 60 s of process, we make 2550 iterations using adaptive time steps.

The numerical prediction under the aforementioned thermal boundary conditions established that a UL SMC is fully reacted after 45 s. This values are in agreement with the industrial process for which the pieces are demanded to be fully processed for a total time below than 60 s.

The simulation describes the temperature evolution in industrial part, where the conduction mechanism heats the piece then triggering the chemical reaction. The simulation predicts a beginning of reaction at about 28.1 s and a full consolidated part at 44 s. The model here presented is numerically stable and can be used for other industrial geometries.

Conclusions

This chapter dealt with the thermo-chemical evolution of SMC materials under compression molding simulations. From the physical and mechanical point of view, the thermal heating of the piece induces changes to the piece far beyond only temperature modifications. In the following, we address the physical mechanisms encountered in the part and how we dealt with that in this work. The hot mold conducts heat towards the SMC part. The standard heat equations for incompressible flows have been extended to account for the compressible behavior of SMC. The temperature evolution induces first a volume change due to thermal dilatation. In this work, we integrated this issue by means of a volumetric term in the mass equation during the resolution of the mechanical problem. This term is integrated implicitly to avoid numerical instabilities. The density is then updated by the resolution of a convection-reaction equation. The chemical reaction of SMC resin does not occur immediately. First, the inhibitors are consumed, translated into a spare time before the beginning of reaction. Therefore, the presented chemical model establishes a dependence on the thermal story of the part. For this reason, we proceed to model an ignition time by solving a convective equation taking into account the spatial thermal history. Once this ignition time τ reached a reference value t r , the reaction is triggered. We model the reaction by a convection equation coupled by a fixed point method to the resolution of heat equation. A set of four equations are then solved to take into account the thermal evolution of the part on thermal and mechanical properties.

Each sections of this chapter are illustrated by examples coming from laboratory or industrial apparatus. The thermo-kinetic model is checked by comparing the temperature profiles with experimental data for both SMC. This comparison proves detection and predictions of the reaction similar as the experimental data. The dilation and contraction of the part coupled to the thermal and kinetic problem is also compared for experimental tests where the SMC sample increases its thickness while being heated and decreases its thickness during the reaction. The model here presented coupled the dilatation terms and predicts the geometrical evolution as the experimental data recorded. For both materials, good agreement was obtained. Further on the modeling, we predict the temperature and heat fluxes evolution during the compression of a sample in plane strain compression. The aim is to estimate the SMC thermal prediction under dynamics conditions. A compression from 6 mm to 2.8 mm is achieved and then the mold is kept closed until reaction occurs. In those cases, the predicted reaction is very similar to the experimental data. The heat fluxes data are compared and a good agreement is found. The modeling of the mold and punch also improves the global numerical simulation of the process.

In general, we conclude that the thermo-kinetic model assigned to the SMC presented in this work describes very well the experimental data and encourages further improvements related to mold/SMC interfaces (RTC and friction conditions). The differences encountered during these comparisons yield on the estimation of the thermal contact resistance, the friction coefficient, the underestimated viscosity coupling with respect to temperature and conversion degree. For further studies, we recommend to look at these variables in order to get a proper fully modeling of the process. In this PhD work, we have addressed the problem of solving compression molding simulation for SMC products, by describing and predicting all thermo-mechanical aspects to model a fully industrial process. To reach this objective, we used an existent Finite Element Library as well as literature models and information representing the latest state of the art on this subject.

After three years of work, a manuscript with four chapters has been written, detailing every macroscopic aspect of SMC products. Literature review on these materials concludes that the more adapted model so far, used is the one declared in [1]. The review also indicated that thermo-mechanical coupling studying SMC materials is not frequent in literature. In fact, only few examples have coupled those mechanisms and none of them adopted appropriate mechanical models accounting for anisotropic behavior due to the concentration of long fibers in the material. Usually authors tend to assign in their viscous material an standard isotropic model such as the power law model. Additionally, coupling mechanical and thermal evolution is the only way to guarantee that the filling of the mold cavity occurs before the triggering of the cross-linking reaction. If the material reacts during the filling, the final part will be fragile. The thermal coupling during and after filling allow us to predict the beginning of the reaction. Clear objectives of this work stated to model a fully thermo-mechanical process of the SMC compression molding on complex geometries found in industrial applications.

Synthesis and conclusion

For accomplishing this task, this work discuss has first presented the mechanical behavior of SMC materials during molding and the kinetic model that describes the cross-linking reaction of the resin (paste). Taking into account for the anisotropy and compressibility are the two main contributions of this work from the mechanical viewpoint. In the following, a description of each aspect is given from both physical and numerical perspectives.

Rheological behavior

When studying the rheological behavior of SMC during the filling of the mold, we have presented a model accounting for the anisotropic behavior of such composites. The model assumed a viscous transverse isotropic behavior supposing a power law for the flow dynamic description. Meaning fibers are widespread along the SMC layers but the fibers are not orientated in the normal direction (usually the axe of compression). The fiber concentration (>15%) induces an homogeneous behavior in the plane but not in the thickness of the stack of SMC layers, translated into anisotropy. In fact, the plane of the fiber is assumed to be isotropic whereas the anisotropy behavior is modeled assuming a different strain rate-stress relation along the normal plane. Numerically, a fourthorder tensor model has been proposed and implemented from the description given in [1]. The Dumont's models for compression molding and the Lipscomb's model for injection molding can be both generalized from our approach.

The model has also been enriched to take into account for the compressible behavior due to porosity closure during the early stages of compression, representing the densification of the material under deformation. When the air is expelled, the material becomes less compressible, until reaching an incompressible state with no-porosity. This process is non-reversible. According to the rheological model proposed to account for this phenomena, a direct relation to the pressure in the mass conservation equation can be found. Thus, the compressibility is assigned to a term proportional to the pressure and a compressible factor. Such factor depends on the volumetric deformation of the product and tends to zero when the materials becomes incompressible. The comparison with experimental data shows the agreement on the shape of the curve during the evolution of the compressible/incompressible stress. An anisotropic compressible solver was implemented in REM3D R that accounts for anisotropic by means of a structural tensor and an evolving compressible behavior by a pressure-proportional term. Behavior which is related to a density evolution of SMC.

Continuum mechanics conservation equations for viscous anisotropic compressible materials are considered. Porosity closure and thermo-kinetical changes have a direct influence in the density evolution of the material. Since the density can be considered as a thermo-dynamic variable, conservation of mass is developed in terms of velocity, pressure, temperature and chemical reaction. Eliminating the direct contribution of the density in the mass equation gives rise to the notion of isothermal compressibility, dilatation, shrinkage. Those parameters are material dependent coefficients that can vary or remain constant.

Mass and momentum equations are solved coupled at a given temperature state of our computation domain. This means that a splitting scheme between the mechanics and the thermal response have been used and implemented, leading to the resolution of a system of equations having as unknowns the velocity and pressure (for a given temperature). Then, the thermal problem is solved.

Numerically, the anisotropic compressible Stokes problem was solved using a Mixed Finite Element method. The computational domain is discretized using tetrahedral elements and an enriched continuous test function (combination of a linear and a piecewise) for the velocity, and a piece-wise linear function for the pressure, temperature and conversion degree fields. The linear system is non-symmetrical and non-linear (due to the viscosity power law relation and the compressible coefficient). Both anisotropic and compressible responses of the implemented solver have been checked with the model predictions.

Heat transfer and chemical reaction -thermo-kinetical model

The chemical reaction of SMC materials is an exothermic reaction ruled by a thermokinetical model accounting for the inhibitors presented in the resin. The SMC materials heat up and after reaching certain thermal state where all inhibitors were consumed, a chemical reaction releasing heat causes the material to consolidates and becomes rigid.

The reaction releases heat increasing the temperature of the part up to a temperature of 50 o C higher than the mold. For a given temperature, the material requires a certain time to react, which depends on the inhibitors presented in the mixture of the resin and on the thermal story of the SMC. The considered thermal model required the resolution of the heat equation coupled to a phenomenological model describing the reaction evolution, enriched with the equation ruling the inhibition time required to consume the inhibitors presented in the SMC mixture. Formulation of the heat equation included the heat coming from the reaction by a source term. The heat released by the reaction was computed directly knowing the density, the chemical energy of the system and the conversion degree evolution given by the Bailleul's model [2], a temperature-reaction dependent model. In order to know when to activate such model, the ignition time needed to be computed. This thermal time depends on the thermal history of the part and once this spare time reaches a reference value the reaction is triggered. Basically, the inhibitors consumption is space and time dependent. To compute the ignition time and to take into account its time-space dependency, a convection equation was solved.

This methodology is consistent with thermo-kinetical model and has proven to predict the reaction of the part in agreement with experimental results.

Numerically, the thermo-kinetical model is solved by an intermediate coupling scheme. The heat and the conversion degree evolution equation given by the Bailleul's model, as well as the convective scheme for the ignition time, are coupled by a fixed point scheme. The chemical reaction lasts few seconds with respect to the thermal cycle. As an example, the reaction of a sample takes approximately 60 seconds to start but lasts around 10 seconds. The reaction profile distribution in time and space is a sharp function. For that reason, a mesh adaptation and time step adaptation strategies were employed and modified for this case.

Thermal and kinetical changes of the piece induce a variation on the local density of the material. For that reason, the density is updated as well, after solving the thermokinetical problem. In the first chapter, we have introduced a technique to use a multi-plicative contribution from the porosity closure and the thermal variations, named relative and dense state densities. The relative density is modified following the volumetric deformation whereas the thermo-kinetical evolution modifies the dense state density, which decreases if the temperature increases meaning an increment on the volume of the piece. When the reaction occurs, the material consolidates, which also increases the density meaning also the part's shrinkage.

Thermal dilatation and chemical shrinkage coefficients are function of the thermodynamic state, mainly depending on the reaction degree and the temperature. These coefficients are also founds to be directional dependent. For that reason, they are written in its tensorial representation. In the fiber plane, dilatation or shrinkage terms differs from the normal direction. The thermo-kinetical strain rate tensors depend on these directions. Meaning that for a given temperature change, the dilatation of the part is higher in the thickness than in the plane.

Finally, a full thermo-mechanical model was here proposed, taking into account all the observed phenomena during SMC compression. It has been proven that it predicts similar results as the one observed in experimental tests and, to our knowledge, is the first model coupling such physical phenomena to be implemented in a numerical tool.

Compression molding applications

At the beginning of this work, REM3D was a software mainly devoted to injection molding applications. Today, the developments and strategies proposed during this PhD allow the software to deal with compression molding problems. In this work, we shown its ability to cover the main stages of the compression molding process: namely, from deposition of the piece, filling of the mold cavity, consolidation of the part and part release. The first result of this work is thus the extension of Rem3D to compression molding problems, from deposition to ejection of the part. Here, we considered also porosity closure, the compression stage, the thermal evolution and the conversion of the viscous material to a rigid part. Determination of physical variables such as pressure, temperature, conversion degree, stresses and reaction time are also possible with our approach.

A splitting technique of the density evolution by porosity closure and by a thermodynamical evolution enables a de-coupling of both contributions, easing the numerical description of the density. Therefore, the information on the relative density can be used to estimate the fully expelling of the air during compression, allowing us to know if the density of the part has evolved for the porosity or for the thermo-kinetical variations, providing important information on what is happening inside the material.

The anisotropy exhibited by the composite is also taken into account by coupling the mechanical solver to a orientation tensor or structural tensor, defined according to the application. The mechanical solver predicts, the stress level, close to the one found in experimental test.

In this work, comparison with the literature and experiments have been performed, validating the thermo-mechanical model and detecting the domains where the approach needs to be improved. To our knowledge, this work has produced the only software that considers fully 3D anisotropic compressible behavior in all stages of the compression molding process simulation. The software is able to provide thermo-mechanical solutions considering the evolution of the temperature and stresses and its coupling.

Perspectives and improvements

At the current stage of development, improvements in the numerical methods, as well as in the physical models can be proposed for future developments. As in the previous case, they are divided in the mechanical, thermo-kinetical and compression molding areas.

For the mechanical resolution...

From the mechanical viewpoint, the stress evolution during the compaction of the SMC sample in the oedometric test shows differences with the model proposed in this project. Basically, these differences may come from the fitting curve used in α 3 or well by the isotropic nature of the stress in the term α 3 tr(D) I. Further models including a transverse isotropic porosity evolution will include two terms acting on the direction of the structural term M . This will allow an orthotropic evolution of the porosity during the densification stage.

A more stable proposition for the equivalent strain rate should be provided. This term is used as input for updating the viscosity and its transition towards the incompressible case has only been considered for the densification. The multiplication of the term α 3 with the trace of the viscoplastic strain tensor tr(D) carries some numerical difficulties when α 3 increases, causing thus numerical instabilities on the computation of the viscosity during the simulation. A formulation implying the integration of the β factor instead seems more logic. However, the dependence of β to the viscosity η and to the equivalent strain rate introduce a non-linear dependency not so easy to explore.

From the numerical perspective, the structural tensor M should be included as a local field. Meaning that for 3D complex geometries, the M tensor needs to be build from the normal to the local fiber plane all along the part. IN this work we assumed that the normal to the fiber plane was always the vertical direction. Consequently, the Folgar and Tucker equation for the evolution of the orientation of fibers may be computed in order to compute the normal of the fiber plane. The principal directions of the orientation tensor will define the fiber plane by noticing the space in which most of the fiber are aligned. Then, a mechanical coupling can be obtained by taking the normal to this plane and by building the tensor M . In that way, a fully 3D anisotropic mechanical problem can be solved for complex geometries.

For the thermo-kinetical resolution... A sensibility study on the boundary conditions will conclude the dependence of the model and of the numerical response to fit experimental data. Some thermal aspects involved during the simulation presented here were not fully discussed on this work. For instance, the thermal contact resistance between the wall and SMC plays a roll on the heat transfer during the compression and holding stages. Additionally, this value can evolve at the beginning of the compression step. In appendix, we have detailed the methodology used to account for the thermal contact resistance and how it can be easily implemented, even in immersed methods. For that, we use a narrow layer of different thermal characteristics close to the contact mold/SMC.

The characterization of the SMC materials on thermal evolution is usually assigned to follow an Arrhenius type law. However, the friction against the mold during the characterization test can induce errors on the parameters associated to such fitting. A fully thermo-mechanical coupling using inverse analysis would produce a more accurate characterization. A joint work on experimental-numerical inverse characterization should be performed, since the model is fully described and the software has been enriched to account for all these physics.

For the compression molding domain... The improvements of the boundary condition imposition play an important roll on the results presented in this work. The mechanical and thermal conditions should be implemented from real industrial conditions or well by real data recorded during experimental tests. The main problem remains on imposing boundary conditions in immersed methods in regions inside the computational domain. Imposing fields inside the domain can cause instabilities on the solution and cause divergence in some cases. A penalization strategy can be implemented for this purpose. The Nitsche method [3] seems to be a good approach. The same methodology for penalization can include a term to account for friction between the SMC and the mold walls. In this way, imposing either slip or no-slip conditions can become more natural when simulating the process. The analysis, characterization and modeling of friction for contact between SMC materials and mold remains in this work an open question as well. The friction coefficient depends on local conditions such as normal stress, relative velocity and temperature. A strategy for its characterization was proposed in [4]. Friction on immersion methods differs from Lagrangian approach since here it has to be supported on a level set approach. Meaning that an implicit function of the surface should carry the friction law and integrate a coupling on the mechanical problem without inducing instabilities. A viscous layer approach was proposed in [5]. Nevertheless for the use of such technique, an imposition of the full immersed body (for the mold and punch) should be considered. This carries instabilities next to the implicit surfaces since the velocity field presents a discontinuity. A more adapted approach should be implemented to deal with this problem.

A level set convection scheme where the mass is conserved. This point has been addressed in this work and a pragmatical solution was used. Works on a conservative level set convection scheme during re-meshing supported on stabilized Finite element methods is a wide open field. The causes of mass lost come from the stabilization of the numerical schemes, the time and mesh discretization, the interpolation error and the interpolation during the re-meshing, as noticed in this work. These are the main points to address when proposing a conservative method. Future work may concern the introduction of a strategy supported on level set approaches to deal with this problem. The complexity of such method implies a fully parametric study to determine the main causes of loss with level set approaches.

The use of a multi-phase approach for the SMC. Here, the fiber and the resin where considered as one single homogeneous phase. The use of a biphasic-model for the fiber and the matrix, including a third phase (the air) should improve the mechanical modeling. The computational cost added to industrial simulations may be compensated by the description of fiber migration. In fact, the compression process induces fiber migration modifying the local fiber concentration. This migration, once the part is consolidated, may affect its mechanical resistance. In order to account for this fiber segregation fiber and matrix should compute a different velocity profile. In a homogeneous approach, as the one presented in this work, the computed velocity represents the velocity of the mixture of fiber-and-matrix. First attempts for a biphasic model were already proposed in [6] and [7] but remained a theoretical model to be integrate in numerical tools.

Furthermore, this work can be improved by assigning an evolution of the viscosity with respect to the conversion degree coefficient during the curing. A modified Castro-Makosco equation is presented in chapter 4 and has been only characterized for the paste and not for SMC materials. Further characterization should conclude if the modified relation remains the same for the composite (paste + fibers).

After consolidation of the part, the SMC is a rigid composite. The behavior of the rigid phase (assumed solid) can be introduced by means of a displacement-based formulation. However, this might lead to numerical problems such as bad matrix conditioning when interacting both solid and viscous phases. In order to have an stable formulation, a multi-phase formulation can be exploited, translating into a fluid-structure problem. Furthermore, an accurate rheological behavior of the solid phase is required to predict the correct residual stress in the part after deformation and equilibrium. An equivalent strain rate can be defined depending on the rheology selected for this layer. In this case, for example if we take the layer to be powerlaw type, we write the equivalent strain rate as: Where α represents the friction coefficient between the two surfaces. Normally, obtained through experimental results. This allows us to make an straightforward link to the approach study here.

The viscous friction approach mentioned in here requires to induce the same stress. From equations B.4 B.5, this relationship is obtained:

α f = K L e m L (B.6)
We recall that the strategy of the relation B.6 vary from the rheology assigned to the boundary layer and the friction law used when characterizing experimentally the phenomenological friction. In any case, the Global definition states: We desire to study the influence of the contact between the two fluids when this contact moves from no-slip to slip. The only parameter to be modified is the viscosity of the boundary layer η L that varies from 100 to 5e -05. In this study we neglected the inertial terms.

τ Layer = τ f riction (B.

Results

For the extreme cases studied in here, we present the velocity field along the domain in Fig. B.3. We visualize how for the no slip case the velocity goes smoothly from zero to one having a constant gradient. Whereas for the slip case, there is a jump in the velocity field from zero to one exactly in the position of boundary layer.

However, this solution do not describe the conditions addressed in this section, since a temperature jump is establish in the interface, we should set:

T C 1 = T C 2 .
For that reason, we add a third body of thickness e between the two previous domains:

Ω 1 Ω 2 0 L 1 L 2 K e T - e T + e e T 1 imposÃľe T 2 imposÃľe
In this way, the new weak formulation stands: making the same conditions applied previously, we get:

K 1 ∇T 1 , ∇w 1 Ω 1 -K 1 ∇T 1 • n 1 , ∇w 1 δΩ 1 = 0 + CL in Ω 1 K e ∇T e ,
T C 1 = T - e et T C 2 = T + e
Which allos us to describe by means of a smooth transition, the temperature jump in the interface between the two sub-domains Ω 1 and Ω 2 . In fact, we get a jump along the thickness e, visualized on the following picture: After having proposed a solution to the problem of the temperature jump. We will then evaluated the value of the conductivity in the thickness associated to the third body in order to guarantee the same heat flux.

Appendix D

Rheological setups for dense and porous SMC -Analytical formula for purely elongational flow D.1 Rheological setups for incompressible case

No slip condition are assumed on the wall, and therefore the stress tensor and the pressure are constant on the whole domain.

D.1.1 Simple compression

Assuming the incompressibility

D 11sc = D 22sc = - 1 2 D 33sc = - 1

D.2 Rheological setups in the compressible case

In the sequel, it is assumed that the equivalent viscosity follows a power law η eq = η 0 D eq D 0 n-1 the parameters η 0 et n depend on fiber concentration. The combination of the two last equations with power law gives

σ 33sc = η 0 α n+1 2 0 1 + α 1 + α 2 + α 3 1 + 2α 3 n+1 2 D 33sc D 33sc D 0 n-1 (D.16)
For α 3 → ∞ we recover the expression of compressible case.

D.2.2 Plane strain compression (ps)

For plane strain compression (there is motion in the 2nd direction) , one has scopic analysis scale it is described. However, it is too small from the rheological point of view. [Ausias,2007]. An homogenisation procedure is used to obtain the equivalent orientation of a representative Volume Element (VER). It is required to used an intermediate scale in which th ehteory of continous mechanics works. This allows to described an tensor orientation description in such VER, explained as follows.

F.0.2 Orientation tensor description -Macro scale

For the population of fibres containing in a VER, a distribution orientation function is used, as the probability to find a fiber aligned to a particular direction. This enables th eformulation fo the macroscopic orientation tensor a (second order tensor) and A (fourth order tensor) described in If we represent the tensor A in that system, most of the components are zero, whereas the others can be expressed as function of (λ 1 , λ 2 , λ 3 ). In fact it is shown in [Cintra & Tucker, 1995] that it can be written as a contracted tensor A nm ⇐⇒ A ijkl . Then, transporting this tensor A nm to the tensor on the main refence system A * nm , we need to define each component A * nm as function of the eigenvalues λ 1 , λ 2 (remember λ 3 = 1 -λ 1 -λ 2 ). for the rest of the terms, we use the polynomial proposition of [Wetzel,1999] base on the coefficients C mm of [VerWeyst,1998].

A * nm = f (C nm , λ 1 , λ 2 ) (F.9)

All the details are found in [Ausias,2007]. It is important to say that this model is the best between all the models described in the literature until the year 2000. However, different models came out due to the complex industrial processes requiring a better description since normal models have fallen in accuracy.

F.0.3 Complex Models for Fiber Orientation Evolution

In the last decade new better adapted models have been proposed to account for fiber orientation prediction. These models correct some acceleration on the orientation evolution that can be found using previous models. They include an anisotropic diffusion coefficient instead of the standard isotropic C I , which describe better the physical interaction between fibers.

R ésum é

Ce travail porte sur la simulation num érique et la mod élisation du comportement thermo-m écanique des mat ériaux composites renforc és par des fibres. Sp écifiquement les mat ériaux SMC (Sheet Moulding Compound) sont utilis és dans le processus de moulage par compression pour construire des pi èces automobiles de haute performance. Ce travail est divis é en quatre chapitres, d écrivant tout d'abord un mod èle thermo-m écanique enti èrement coupl é pour les mat ériaux SMC standards et innovants à haute concentration en fibres (> 25% en volume). Le SMC est trait é comme un m élange incompressible de fibre et de r ésine compl ét é éventuellement par une phase de porosit é compressible.

Son anisotropie est mod élis ée au moyen de tenseurs structurels. La cin étique de r éaction et de consolidation de la pi èce est également mod élis ée et étudi ée. Les donn ées exp érimentales m écaniques et thermiques enregistr ées sur des échantillons de mat ériaux SMC sont compar ées au mod èle et à la solution num érique fournie par ce travail. D'un point de vue num érique, nous utilisons la m éthode des domaines immerg ées o ù chaque phase est distingu ée par une fonction distance sign ée. Nous d écrivons le proc éd é de moulage par compression en proposant une r ésolution compressible anisotrope unifi ée capable de d écrire la transition compressible / incompressible du mat ériau SMC sous d éformation. Cela permet de d écrire la r éponse m écanique du SMC et de pr édire localement la consolidation (durcissement) de la pi èce le long du cycle thermique.

Mots Cl és

Abstract

This work deals with the numerical simulation and modeling of thermomechanical analysis of fiber reinforced composites materials. Specifically for SMC (Sheet Molding Compound) materials that are used in compression molding processes to build automotive high performance parts. The work is divided into four chapters, firstly describing a fully coupled thermo-mechanical model for standard SMC materials and for innovative SMC with high fiber concentration (> 25% in volume). The SMC is treated as an incompressible mixture of fibers and paste complemented by a compressible porosity phase. Its anisotropy is modeled by means of structural tensors. Kinetic of reaction and consolidation of the part is also modeled and studied. Mechanical and thermal experimental data recorded on samples of SMC materials are compared to the model and numerical solution provided in this work. A numerical framework, we use the immersed boundary method and the level set method. We describe the compression molding process by proposing an unified anisotropic compressible resolution able to describe the transition between compressible/incompressible of SMC materials under deformation. We are able to describe the mechanical response of the SMC and to predict locally the consolidation (curing) of the part throughout the thermal cycle.

4. 4

 4 Compression, heating and reaction of SMC . . . . . . . . . . . . . . . . 4.4.1 Analysis in the Mold and Punch Surfaces . . . . . . . . . . . . . . 4.4.2 Analysis of thermal evolution inside the SMC . . . . . . . . . . . 4.4.3 Evolution of velocity profile during plane compression for high performance SMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Heat transfer and reaction in industrial like geometries . . . . . . . . . . 4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Possible ways to reduce CO2 emission.

Figure 2 :

 2 Figure 2: View of potential parts to be replaced by fiber reinforced composites in order to gain up to 110 Kg of weight saving. (Image courtesy of Plastic Omnium)

Figure 6 :Figure 7 :

 67 Figure6: Relative time for each stage of molding compression process[1]

Figure 8 :Figure 9 :

 89 Figure 8: Examples products made of composites

Figure 10 :

 10 Figure 10: Tomography of SMC studied in this work. Left SMC UL -Right SMC HP (High profile). In white the fiber bundles, in black the air entrapped in gray the paste.

Figure 12 :

 12 Figure 12: Photography showing chopped glass fiber bundles falling onto the plastic film covered by paste resin

Figure 13 :

 13 Figure 13: Evolution of the porosity as function of the equivalent deformation. High Performance SMC are compressible whereas Ultra Light one can be consider as incompressible

Figure 14 :

 14 Figure 14: Project launched by Plastic Omnium Auto exterior

26 1. 1 . 2 27 1. 1 . 3 31 1. 1 . 4 32 1. 1 . 5 34 1. 2 35 1. 2 . 1 37 1. 2 . 2 1 . 3 1 . 3 . 1 1 . 4 Time adaptation based on a bounded increment gradient . 48 1 . 4 . 1 1 . 5

 261227133114321534235213722131311414115 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.1.1 Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . . Mass conservation and mechanical problem . . . . . . . . . . . Energy equation . . . . . . . . . . . . . . . . . . . . . . . . . . Summary of conservative equations for compressible reactive flows Discretization and numerical framework . . . . . . . . . . . . Mixed variational formulation and discretization -mechanical problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Stabilization of the mechanical problem . . . . . . . . . . . . . 40 1.2.3 Computation of local matrices -momentum equation . . . . . . 42 1.2.4 Computation of local matrices -mass equation . . . . . . . . . 42 1.2.5 Resolution and stabilization of the thermo-kinetical problem . 42 Space and time adaptation . . . . . . . . . . . . . . . . . . . . 44 Mesh adaptation on a edge-based error prediction . . . . . . . 44 Resolution and optimization . . . . . . . . . . . . . . . . . . . . 49 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Figure 1 . 3 :

 13 Figure 1.3: Strategy for body immersion with an adaptive mesh procedure. The meshes, Lagrangian in the punch and Eulerian in the air, die and preform

2 . 4 2 . 4 . 1 2 . 5 2 . 6

 242412526 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Background in Level set methods . . . . . . . . . . . . . . . . . 2.1.2 Narrow band level-set function . . . . . . . . . . . . . . . . . . 2.2 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Linear hyperbolic tangent . . . . . . . . . . . . . . . . . . . . . 2.2.2 Re-initialization velocity -A local description . . . . . . . . . 2.2.3 Effect of mesh adaptation on volume conservation . . . . . . . 2.3 Volume/mass Conservation on level set methods . . . . . . 2.3.1 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Background in conservative strategies . . . . . . . . . . . . . . 2.3.3 Proposed strategy . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Strategy to avoid penetration of recovered-level-set into tools . 2.3.5 Maximum volume loss per iteration guaranteeing accurate displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.6 Exact displacement of the hyperbolic tangent level set . . . . . Application of the conservative strategy in an industrial geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mesh adaptation on multiple fields . . . . . . . . . . . . . . . . 2.4.2 Test Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 short-shots comparison . . . . . . . . . . . . . . . . . . . . . . . Stable and robust strategy for the computation of compression force in immersion domains . . . . . . . . . . . . . . . . 2.5.1 Viscous Power and Compression Force . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

( a )

 a Visualization of the finner discretization along the interfaces of the body, where we notice a thick region of small mesh sizes. (b) Visualization of the finner discretization in a thickness E. The interface is located in the middle of the thickness.

Figure 2 . 2 :

 22 Figure 2.2: Immersed Body Approach: the composite and the dies are immersed and tracked by a level set approach coupled with mesh adaptation.

Figure 2 . 3 :

 23 Figure 2.3: Representation of level set function φ as a distance and the truncated level set function φ

  Figure 2.5: When using E = 0.15 and e = 0.1 we define level set eq. (2.21) and its respective Heaviside eq. (2.22)

  (a) -constant booster t=0s (b) -constant booster t=1s (c) -constant booster t=2s (d) -local booster t=0s (e) -local booster t=1s (f) -local booster t=2s

Figure 2 . 6 :

 26 Figure 2.6: Deformation of the circle and convection using a constant (a-b-c) and local defined λ (d-e-f).

Figure 2 . 7 :

 27 Figure 2.7: Initial and final configuration of the simple compression case. Final deformation -ln(h f /h o ) = 0.8. (a) and (b) velocity profile in the piece at the initial and final position. (c) and (d) cross section of the compression of the fluid.

Figure 2 . 8 :

 28 Figure 2.8: Influence of mesh adaption on volume conservation. h min = 1e -03mm, ∆t = 1e -02s

  (a) Compression of an Industrial part. Punch Mold Preform Body Immersion (b) Schematic illustration immersion of different phases: punch, composite and die.

Figure 2 . 9 :

 29 Figure 2.9: Industrial case presentation. (a) the three bodies defining the mold, punch and the geometry of the preform. (b) Schematic of the immersion in our numerical tool.

  (a) View of incomplete mold filling during reference compression case due to material loss.

  Volume variation during Simulation for the reference case.

Figure 2 . 10 :

 210 Figure 2.10: Volume loss on an industrial part Compression from 6mm to 2.5mm. Without any optimization, a global 30% of loss is observed in particular by the incomplete filling of the mold cavity.

Figure 2 .

 2 Figure 2.11: relative error eq. (2.31) in volume as a function of the sample thickness during compression of an industrial piece using the volume recovery strategy

Figure 2 . 12 :

 212 Figure 2.12: Schematic illustration of the displacement of the iso-contour zero after applying recovery our algorithm

2 Figure 2 . 13 :

 2213 Figure 2.13: Schematic illustration of the displacement of iso-contour zero inside the Fine Meshing Zone (FMZ), after applying the recovery algorithm

( a )

 a View of the mesh adapted (b) Adaptation respect to Dirac function

Figure 2 . 14 :

 214 Figure 2.14: Mesh view at the boundary of domain. Three bodies immersed in a single numerical domain (a). Close up of mesh adaptation following the Dirac function on the surfaces (b).

Figure 2 . 16 :

 216 Figure 2.16: Initial view of the immersed bodies. Information on numerical discretization, CPU used and level set thickness definition e,E.

Figure 2 . 18 :

 218 Figure 2.18: View of the deformation during the descending of the punch. Velocity profile of the part. The cavity of the mold is filled, while the thickness of the part is reduced. Test on a reduced volume of industrial geometry.

Figure 2 . 19 :

 219 Figure 2.19: Volume lost on industrial piece after deformation of the part. 30% lost observed by the incomplete filling of the mold cavity (a). 0.1% lost after implementation of recovery algorithm using Newton method (b).

Figure 2 .

 2 Figure 2.20: Molding Compression -Comparison of flow front real piece / simulated part. At a thickness of 2.8mm.

Figure 2 . 21 :

 221 Figure 2.21: Compression force for the compression of a cylinder Figure 3.18a. Comparison of analytical case[39] for the slip and no-slip case with numerical predictions of this work.

3 . 3 3 . 4 3 . 5 3 . 6

 33343536 Incompressible SMC . . . . . . . . . . . . . . . . . 99 3.2 Modeling Compressible SMC . . . . . . . . . . . . . . . . . . 99 Unified Compressible-Incompressible Model for SMC . . . 101 3.3.1 Viscosity model for SMC materials . . . . . . . . . . . . . . . . 3.3.2 The Compressible/Incompressible transition -Evolution of parameter α 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Evolution of porosity and fiber fraction concentration . . . . . Summary of mechanical resolution . . . . . . . . . . . . . . . 108 Numerical cases . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3.5.1 Standard SMC -incompressible case . . . . . . . . . . . . . . . 3.5.2 High performance SMC -Compressible case . . . . . . . . . . . 3.5.3 Comparison of stress level with experimental data . . . . . . . 3.5.4 Compressible/Incompressible Transition for a Squeeze Flow . . 3.5.5 Plane Strain Compression . . . . . . . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

Figure 3 . 1 :

 31 Figure 3.1: Matrix and fibers are modeled as a homogeneous mixture

Figure 3 . 2 :

 32 Figure 3.2: Porosity in tomography pictures for a cylinder sample of HP SMC (diameter 50 mm and height 6 mm). Evolution of the stress, volume strain as function of the axial deformation in a simple compression. .

Figure 3 . 3 :

 33 Figure 3.3: Evolution of the reference viscosity as a function of the temperature. Arrhenius law with b = 4900 T R = 296 K.

0. 27 Figure 3 . 5 :

 2735 Figure 3.5: Evolution of compressible state in SMC as function of the volume fiber fraction: (blue line) Dimensionless rheological parameter α 3 ; (green line) numerical parameter β. Parameters for the SMC HP are given in Table3.3. α 3 is the parameter proposed by[29] whereas β is the parameter used for the computations. They are linked by the relation eq.(3.10)

4 :( 1 -

 41 Parameters needed to compute the volume fractions of HP and UL SMCs. new concentrations of fiber, paste and porosity are obtained. Firstly all initial concentrations are calculated knowing the initial porosity state φ po = 0.25 and the dense fiber concentration φ dense f = 0.38 (known by the fabrication process). The initial fiber concentration is then obtained by φ f o = φ dense f φ po ) and the matrix concentration by φ mo = φ dense m (1 -φ po ) or by the mass conservation as φ mo = 1 -φ f o -φ po .

Figure 3 . 6 :

 36 Figure 3.6: Evolution of the porosity, fiber and matrix volume concentration with respect to the volumetric deformation of HP SMC.

Figure 3 . 7 :

 37 Figure 3.7: Sketchs of geometries studied in our numerical applications: simple compression, plane strain and oedo test.

Figure 3 . 9 :

 39 Figure 3.9: Oedo configuration for compressible SMC: the SMC sample is blocked in a ring in oder to allow only vertical deformation.

Figure 3 . 10 :Figure 3 . 11 :

 310311 Figure 3.10: Compaction of SMC sample: the densification occurs due to porosity closure; Visualization of the sample contained by the mold and punch (left); 2D visualization of a middle plane cut showing the pressure at an intermediate time (right).

(a) t 1 : 2 (c) t 3 : 28 Figure 3 . 12 :

 12328312 Figure 3.12: Porosity concentration evolution during oedo compression: compressible/incompressible transition of SMC.

(a) t 1 : 28 Figure 3 . 13 :

 128313 Figure 3.13: Pressure profile evolution during oedo compression and the compressible/incompressible transition: the pressure increase exponentially when the material becomes incompressible.

  15a and Figure 3.15b where the vertical and the lateral stresses during the compaction are plotted. Thus, that validates the numerical implementation of the set of equations modeling the

(a) t 1 :

 1 compressible velocity profile (b) t 2 : velocity profile at ε = 0.2 (c) t 3 : compressible/incompressible transition ε = 0.22 (d) t 4 : incompressible velocity profile -end of simulation -squeeze flow

Figure 3

 3 Figure 3.14:

0 5 • 10 -

 510

  b) Lateral stress during porosity closure.

Figure 3 . 15 :

 315 Figure 3.15: Stress profile evolution during Oedo compression for the case 1. Comparisons between model predictions, simulation results for 2D and 3D cases and experimental data.

0 5 • 10 -Figure 3 . 16 :

 510316 Figure 3.16: Stress profile evolution during Oedo compression for the case 2. Comparisons between model predictions, simulation results for 2D and 3D cases and experimental data.

Figure 3 . 18 :

 318 Figure 3.18: Numerical immersion of the different phases (punch, mold and SMC) in the computational domain: (a) for the simple compression case, (b) for the plane compression.

n+1 2 D n 33 ( 3 Figure 3 . 19 :

 2333319 Figure 3.19: Compression of HP SMC. Porosity concentration in the material at the beginning of compression. Notice the evolution until φ p ≈ 0 when volumetric deformation equals critical state. Compressible/incompressible transition as critical state is reached. Complete squeeze flow behavior at the end (incompressible).

Figure 3 . 20 :

 320 Figure 3.20: Compressible/incompressible transition during porosity closure: evolution of the compressible term β defined in eq. (3.10)

Figure 3 . 21 :

 321 Figure 3.21: Evolution of the viscoplastic volume deformation ε v and the volume concentration in compression test during compressible/incompressible transition of SMC: slip boundary conditions (continuous line); No-slip boundary condition (dashed line).

Figure 3 . 22 :

 322 Figure 3.22: Normal stress during compressible/incompressible squeeze flow. For the analytic solution, the slip case is reproduced by the numerical simulation. The non-slip case reproduces higher normal stress due to the shearing of the material. In both cases, the stress during the compressible step increases linearly with the volumetric deformation.

Figure 3 . 23 :

 323 Figure 3.23: Plane strain: Compression of SMC within a channel.

σFigure 3 . 24 :

 324 Figure 3.24: Compressible/incompressible transition for plane strain compression: (a) evolution of viscoplastic volume strain ε v as a function of vertical deformation; (b) evolutions of stresses for no-slip numerical computations (NS) and slip analytic solutions.

Figure 3 . 25 :

 325 Figure 3.25: Compressible/incompressible transition for plane strain compression: (a) evolution of viscoplastic volume strain ε v as a function of vertical deformation; (b) evolutions of stresses for no-slip numerical computations (NS) and slip analytic solutions.

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: Thermo-kinetic function G(α) for ultra light (UL) and high performance (HP) SMCs.

Figure 4 . 3 :

 43 Figure 4.3: Reaction time at constant temperature: time taken for the material to react according to the temperature.

Figure 4 . 6 :

 46 Figure 4.6: Reaction, heat released and temperature evolution during crosslinking process.

Figure 4 . 7 :

 47 Figure 4.7: (a) Temperature profile in Kelvin and reaction degree (b) along the sample during beginning of cross-linking reaction.

Figure 4 . 8 :

 48 Figure 4.8: (a) Numerical temperature compared with the records of measured temperature in PVT-α test. (b) The numerical simulation predicts the evolution of the cross-linking reaction.

  (a) Cross section: induced velocity in mm/s of SMC sample due to thermal gradient. (b) Density kg/m 3 along the piece for SMC UL due to consolidation process.

Figure 4 . 9 :Figure 4 . 10 :

 49410 Figure 4.9: (a) Velocity profile during reaction. In (b) density along the sample during consolidation.

  Conversion degree evolution according to Bailleul Model.

Figure 4 .

 4 Figure 4.11: (a) Numerical temperature compared with recorded measured temperature in PVT-α test. (b) The numerical simulation predicts the evolution of the cross-linking reaction.

Figure 4 . 12 :

 412 Figure 4.12: (a)Evolution of sample thickness during PvT-α test. Comparison measurements with numerical prediction. (b) Numerical prediction of density evolution during experimental test.

Figure 4 . 13 :

 413 Figure 4.13: Influence of the thermal parameters accuracy on maximum temperature during reaction and on variation on beginning of reaction

Figure 4 . 14 :

 414 Figure 4.14: Plane strain compression velocity. SMC materials flows within a channel. Snaps of deformation for various times

Figure 4 . 15 :

 415 Figure 4.15: Sketch of experimental setup: definition of direction (e 1 e 2 e 3 ) and positions of sensors.

Figure 4 . 16 :

 416 Figure 4.16: Thickness evolution during isothermal test

  Compression stress mold side center position for SMC UL, T=50 o C, h o = 7Compression stress mold side center position for SMC HP, T=50 o C, h o = 5.92mm

Figure 4 . 17 :

 417 Figure 4.17: Evolution of stresses during channel compression for SMC HP and UL at isothermal conditions of 50 o C and 80 o C without reaction.

  (a) Compression stress mold side lateral position for SMC UL, T=80 o C, h o = 7.27mm (b) Compression stress mold side lateral position for SMC HP,T=80 o C, h o = 5.92mm

Figure 4 . 18 :

 418 Figure 4.18: Evolution of stresses during channel compression for SMC HP and UL at isothermal conditions of 50 o C and 80 o C without reaction.

Figure 4 . 19 :

 419 Figure 4.19: Thermo-mechanics of plane strain compression: views of velocity and temperature during the thermo-mechanical simulations

  Evolution of pressure and thickness

Figure 4 . 20 :

 420 Figure 4.20: Comparisons between experimental measurements and simulations of UL SMC for a plane compression: (a) temperature; (b)pressure and thickness.

  Evolution of pressure and thickness

Figure 4 . 21 :

 421 Figure 4.21: Comparisons between experimental measurements and simulations of HP SMC for a plane compression: (a) temperature; (b)pressure and thickness.

  Figure 4.23b. With respect to UL SMC, the LVTD sensor does not record a thermal dilation before the chemical shrinkage for HP SMC. Total heat flux recorded on sensors.

Figure 4 . 22 :

 422 Figure 4.22: Evolution of Heat Flux recorded during experience of channel compression for UL SMC. Comparison with numerical prediction. Continuous line stands for experimental data and dashed lines for numerical results.

  Figure 4.25.We notice the evolution of the temperature in Figure4.24a in two different positions( at the surface and in the middle of SMC sample). The core is highlighted in red and the 0

Figure 4 . 23 :

 423 Figure 4.23: Evolution of Heat Flux recorded during experience of channel compression for HP SMC. Comparison with numerical prediction. Continuous line stands for experimental data and dashed lines for numerical results.

Figure 4 . 24 :

 424 Figure 4.24: Evolution of temperature and reaction recorded during plane compression ofUL SMC.

  Reaction evolution in plane strain compression

Figure 4 . 25 :

 425 Figure 4.25: Evolution of temperature and reaction recorded during plane compression of HP SMC.

Figure 4 .

 4 26(b) the punch compresses the sample and the the velocity profile is linear during the porosity closure. The porosity of the HP SMC basically allows deformation along the thickness. Notice it goes from it maximum (punch speed) towards zero (mold speed). In Figure4.26(c) the compressible/incompressible transition is observed on the velocity profile. By looking at the common velocity profile found in squeeze flow motion (d), the transition velocity profile while the porosity concentration is not 0 but not high enough, is a combination of the profiles pictured in (b) and (d).

Figure 4 . 26 :

 426 Figure 4.26: Velocity profile evolution of HP SMC

Figure 4 . 27 :

 427 Figure 4.27: Compression of industrial part.

Figure 4 .

 4 Figure 4.28: Molding Compression -Dedicated mold

  (a) t = 6.8s (b) t = 17.5s (c) t = 22.5s (d) t = 25.7s

Figure 4 . 30 :

 430 Figure 4.30: Temperature evolution of UL SMC heats by walls maintained at temperature of T = 420 K: the temperature profiles are shown on a cut section before the beginning of chemical reaction.

  (a) t = 28.1 s -Temperature (b) t = 28.1 s -Reaction (c) t = 32.8 s -Temperature (d) t = 32.8 s -Reaction

Figure 4 . 31 :

 431 Figure 4.31: Temperature and Conversion degree evolution snapshots during crosslinking reaction. The temperature in the piece increases due to exothermic reaction.

  (a) t = 40.1 s -Temperature (b) t = 40.1 s -Reaction (c) t = 44.3 s -Temperature (d) t = 44.3 s -Reaction

Figure 4 . 32 :

 432 Figure 4.32: Temperature and conversion degree profiles on cut planes in the last five seconds of the reaction.

1

 1 Synthesis and conclusion . . . . . . . . . . . . . . . . . . . . . 171 5.2 Perspectives and improvements . . . . . . . . . . . . . . . . . 175 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177

4 )

 4 write the shear stress τ xy as function of the consistency of the material (K L ) and the power index m:In literature, friction can be described, for instance, using the Norton friction Law[ref] or Coulomb friction Law. The tangential friction stress exerted to body 1 from body 2 τ 1-2 states: τ 1-2 xy = α f v m r (B.5)

7 )Figure B. 2 :η 1 η 2 > 1

 7221 Figure B.2: Simple Shear test performed on square 2x1 (m) containing two Newtonian fluids with η 1 /η 2 = 4

2 Figure C. 1 :

 21 Figure C.1: Temperature field along the bodies in contact : Contact thermal resistance modeling by the insertion of a third thin body of thickness e.

whereF 2 eq = α 0 1 + α 1 + α 2 + 1 2 D 2 33 1 + α 2 )η eq α 0 1 + α 1 + α 2 F 3ps h L 0 l 0 h 0 ; σ 22ps = F 2ps l 0 h 0 2 33 3 η 3 η

 21122121120233 3sc is the axial force measured on the upper wall. Starting with the model expresses in Chapter 3, one getsDThe equation σ 11sc = 0 gives p = α 0 η eq D 11sc (atmospheric pressure is neglected) and finally the axial stress tensor and the total pressure areσ 33sc = η eq α 0 (1 + α 1 + α 2 ) D 33sc -D 33sc + p (D.2)orσ 33sc = η eq α 0 1 + α 1 + α 2The incompressible model expresses in the first section givesD 2 eq = α 0 (2 + α 1 + α 2 ) DAs in the previous section, the equation σ 11sc = 0 gives p = α 0 η eq D 11ps (atmospheric pressure is neglected). Finally the non null components of stress tensor and the pressure areσ 33ps = η eq α 0 (1 + α 1 + α 2 ) D 33ps -p (D.5) σ 22ps = -p (D.6) p = -1 eq α 0 (α 1 + α 2 ) D 33ps + p (D.7) or σ 33ps = η eq α 0 (2 + α 1 + α 2 ) D 33ps (D.8) σ 22ps = η eq α 0 D 33ps (D.9) p = -1 eq α 0 (3 + α 1 + α 2 ) D 33ps (D.10)

η 0 = 2 . 2 D. 2 . 1 +α 3 (Finally σ 33sc = η eq α 0 1 + α 1 + α 2 3 η 2 eq = α 0 1 + α 1 + α 2

 02221311232112 η 00 (φ -φ c ) Different apparatus proposed here correspond to M = e 3 ⊗ e 3 andD 2 eq = α 0 D 2 33 (1 + α 1 + α 2 ) + D 2 11 + D 2 22 + α 3 (D 11 + D 22 + D 33 ) Simple compression (sc)For a simple compression and perfect slip conditionσ = D 11sc + D 22sc + D 33sc ) σ 11sc = σ 22sc = 0 (the pressure of atmosphere is neglected), one getsD 11sc = D 22sc =eq α 0 1 + α 1 + α 2 + α 3 1 + 2α 2 D 33sc (D.14)and D

2 eq = α 0 1 + α 1 + α 2 Finally σ 33ps = α 0 η eq 1 + α 1 + α 2 33ps = η 0 α n+1 2 0 1 + α 1 + α 2 D. 2 . 3 = α 0 [ 1 + α 1 + α 2 + α 3 ] D 2 33oc(D. 27 )α 1 + α 2 + α 3 α 3 (= σ 22oc = η 0 α n+1 2 0α 3 [ 1 + α 1 + α 2 + α 3 ]σ 33oc = η 0 α n+1 2 0[ 1 +D. 2 . 4 +α 3 (σFigure E. 1 :

 211211221122311222713231123212431 Figure E.1: Influence of the thermal parameters accuracy on maximum temperature during reaction and on variation on beginning of reaction

  [Advani & Tucker,1987]. a = p α ⊗ p α ψ(p α ) dp α (F.2)A = p α ⊗ p α ⊗ p α ⊗ p α ψ(p α ) dp α (F.3)The equivalent orientation evolution equation for the macroscopic tensor is given by: aD -2A : D) (F.4) Later improved by [Folgar & Tucker, 1984], which takes into account the particles interactions (for non -diluted regimes ), by adding a diffusion term. aD -2A : D) + 2C I D eq (I -3a) (F.5) After many experiments the diffusion term C I can be obtained from an empirical relation : C I = 0.0184exp(-0.7148f r) [Bay,1991]. Being f the fiber fraction in volume and r the radius of the cylinder fiber. The relation F.5 has a complex a non linear term A : D. Since the computation of the fourth order tensor, requires many different mathematical treatments, an approximation of the tensor A has to be done, usually called a closure approximations. Closure Approximations Nowadays, the two families of closures the most recent are the natural closure [Verleye & Dupret, 1994] and the orthotropic one [Wetzel & Tucker, 1999]. Both coming from the same basics. The eigenvalues of the tensor a (λ 1 , λ 2 , λ 3 ) are used to compute the the fourth order tensor A. Where the eigenvectors are the same as the tensor a.

A * 11 =

 11 f 11 (λ 1 , λ 2 ) (F.6) A * 22 = f 22 (λ 1 , λ 2 ) (F.7) A * 33 = f 33 (λ 1 , λ 2 ) (F.8)

Figure F. 1 :

 1 Figure F.1: Influence of fiber on mechanical field, difference of a simulation in which fiber orientation does not modify the rheology of the fluid (left), simulation when coupling both effects (right).

Figure F. 2 :

 2 Figure F.2: Closure approximation effect on fiber orientation response.

Figure F. 3 :

 3 Figure F.3: Representation of the fiber in its main direction using tensorial properties

  .3 and Figure 3.4 respectively.

	Parameters for evolution of η		constitutive parameters
	Material HP SMC UL SMC	Material SMC HP SMC UL
	η 0	20M P a s 120M P a s	α 0	1.06	1.06
	φ c	0.02	0.02	α 1	1.36	1.42
	b	3200	4900	α 2	-1.92	-1.97
	T o	296	296	n	0.25	0.35
	Table 3.1: Parameters for viscosity	Table 3.2: Rheological parameters of SMC

Table 3 .

 3

	Parameters for evolution of α 3
	Material SMC HP	SMC UL
	k	6e+05	incompressible
	b	0.5	-
	c	3.2	-
	H	0.446	0.446
	Table 3.3: Parameters needed to computation of α 3 and coming from experimental mea-
	surements [29].
	Parameters for the evolution of φ p , φ f , φ m
	Material SMC HP	SMC UL
	φ dense f	0.38	0.24
	φ dense m	0.62	0.76
	φ po	0.25	0.026
	φ mo	0.47	0.74
	φ f o	0.28	0.234
	δ	1.0	1.0
	γ	0.01	0.01
	ε vcrit	0.25	0.0265

  Numerical predictions for Simple Compression, Plane Strain Compression and Shear. Computations are made for standard SMC having fiber volume fractions φ = 3.5% and φ = 14.7% at room temperature.

		10 1	Model eq. (3.41)	14.7%
			Simple compression
	eq [M P a]	10 0	Plane strain compression Shear	3.5%
	σ
			D eq [s -1 ]
	Figure 3.8: Incompressible Model for SMC: Comparisons between Experimen-
	tal/Analytical/

Table 3 .

 3 5: Experimental conditions for the densification case of HP SMC

			20 o C 9.5 7.2 90 KN
	2	0.03 s -1	20 o C 9.5 7.3 90 KN

in Kelvin for conductivity and in Celsius for rest of propertiesTable 4 .

 4 1: Thermal properties for HP and UL profiles

		Thermal parameters for SMC under study
	parameters	HP SMC	UL SMC
	X f	0.5	0.28
	c p,f	1.247 T + 781.0	1.247 T + 781
	c p,cured	4.77 T + 1045.0	3.63 T + 1095
	c p,uncured	3.43 T + 1246.0	2.64 T + 1386
	κ cured	1.0E -04 T + 0.571	0.25 + 0.0005 T
	κ uncured	5.0E -04 T + 0.497	0.28
	Temperature

Table 4 .

 4 . The values for parameters k ref , A, T ref , a i , t ref , B depend on the resin composition. They are determined by R. Cardinaud and they are reported in Table 4.2. 2: Kinetical parameters for UL and HP SMCs.

	Material HP SMC	UL SMC
	∆H α	76000 J/kg 142560 J/kg
	-	K(T )
	k ref	0.05055	0.08246
	A	46.7518	32.8898
	T ref	393.0K	373K
	-	G(α)
	α 0	7.5590e-04	5.0932e-05
	α 1	3.7372e-01	2.2797e-01
	α 2	-1.1417	-7.7730e-01
	α 3	1.6229	1.2529
	α 4	-1.2695	-1.0881
	α 5	4.1529e-01	3.8455e-01
	-	h(T, t)
	t ref	98.1	325.14
	B	78.76	68.85

Table 4 .

 4 1, the density is split in two terms called the relative density ρ r and the dense density ρ d . The relative 3: Anisotropic Dilatation and shrinkage coefficients for UL and HP SMCs.

	Anisotropy in SMC materials during volume variations
	Material		HP SMC UL SMC
			uncured
	dilatation coefficient in the plane χ plane t	5.6e-5	5.7e-5
	dilatation coefficient through the thickness χ normal t	2.8e-4	2.86e-4
			cure
	dilatation coefficient in the plane χ plane t	3.3e-5	3.8e-5
	dilatation coefficient through the thickness χ normal t shrinkage coefficient in the plane χ plane α shrinkage coefficient through the thickness χ normal α	1.66e-4 -1.8e-3 -3.6e-2	1.89e-4 -2.5e-3 -5.09e-2
	Density references for SMC under study
	Material	HP SMC UL SMC
	ρ(φ p = φ po , T ref , α = 0)	1550	1367
	ρ(φ p = 0, T ref , α = 0)	1940	1400
	ρ cured	2000	1444

Table 4 .

 4 4: Values of reference density for UL and HP SMCs.

Table 4 .

 4 5: Molding conditions for both UL and HP SMCs

		v/h	T imp	h o	h f	t r	F max
	UL	0.3s -1 150 o C 7.34 2.6 15 s 9 kN
	HP	0.3s -1 150 o C 6.37 2.16 15 s 9 kN

Table 4 . 6

 46

: Mold thermal properties

  ∇w e Ωe -K e ∇T e • n e , ∇w e δΩe = 0 + CL in eK 2 ∇T 2 , ∇w 2 Ω 2 -K 2 ∇T 2 • n 2 , ∇w 2 δΩ 2 = 0 + CL in Ω 2

the symmetric gradient tensor is defined as ε = 1

∇v + ∇ t v

In the work of Coupez[21], it was proved that factor 1/d actually is replaced by p/p + 2, which measures a limit condition on the mesh when forcing a fixed number of edges. Here, for sake of simplicity, we present the strategy for a non-limit edge condition. The reader is referred to[21] for details when the number of edges is fixed. The procedure evoked here remains the same as in the original paper.

Summary

The study of the thermo-kinetical behavior of SMC materials is object of study in this chapter. The hot molds walls heats up the SMC increasing its temperature. The temper-

Kinematic reaction of SMC in a PVTα apparatus

The PVT-α mold (Figure 4.4) is a device developed in the kinetic laboratory of Nantes [26] The thermo-kinetic model is validated by comparing experimental recorded measurements with numerical predictions for core temperature and sample thickness variation. The PVT-α apparatus contains a thermocouple located at the center of sample. The thermocouples enable the recording of temperature evolution along the heating and the LVDT sensor measures the thickness variation along the heating and during curing.

Numerical modeling

The PVT-α test was simulated in our numerical platform. The mold, punch and SMC sample are immersed in a computational domain of dimension 100mm×100mm×100mm. In Figure 4.5, the numerical representation of the PVT-α is given. The computational domain contains four phases. A rigid mold with conductivity κ, in which its rigid motion is modeled by a high viscous fluid, here η mold = 10 3 η SM C . The mold embeds the whole sample of SMC from the bottom to the lateral along the circular surface. The punch is modeled by two bodies. A rigid punch with the same characteristics as the rigid mold and a moving punch. The moving punch is used to allow the SMC to dilate with the heating. The thermal conductivity of the moving punch remains the same as the rigid punch, however, we set a lower viscosity η d for this body, such as, η d = 10 -2 η SM C . The SMC is placed between the mold and the punch and the thermal contact resistance is modeled by an additional small layer with a specific conductivity. The details of this modeling is explained in Appendix: C. The heating in during the simulation is obtained by the updating of the temperature on the boundary of the computational domain ∂Ω. We imposed a temperature T (t) in ∂Ω of the form, T (t) = 25 o C + 3 o C/min × t. We imposed v = 0 along ∂Ω. We use a mesh with 4000 nodes and adaptive time step. We perform the simulation for both SMC. The temperature evolution in the center of sample is recorded. For estimating the thickness evolution along the simulation, we follow the vertical coordinate of the level set in the symmetry plane.

In Figure 4.6, we follow through the snapshots the locations where the reaction starts and where it ends. The pictures prove that the reaction occurs first in the surface and evolves progressively whereas in the core it evolves faster.

Study of UL SMC

For the UL SMC, the reaction begins around 1500 s and lasts approximately 80 s. We observe in Figure 4.7 that temperature along the sample increases due to the exothermic reaction. The heat is conducted towards the SMC sample until triggering the chemical reaction translated into the evolution of the conversion degree. The thermal evolution is depicted during two steps. First, the heating part is dominated by only conduction heat transfer. There, several snaps depict the evolution of the temperature along the piece in Figure 4.30 Second, once the reaction is triggered as shown in Figure 4.31, the snaps picture the evolution of the temperature and the conversion degree α in Figure 4.32.

During the heating up of the part from the hot walls, the evolution on temperature is shown in Figure 4.30. The three dimensional piece is cut in by a plane along its largest dimension to show the thermal evolution within the thickness. Snaps at times equals to 6.8 s, 17.5 s, 22.5 s and 25 .7s show the heating time of the part. Notice the temperature range inside the interval [320 K, 420 K].

Secondly, at around 28.1 s the chemical reaction starts to become visible. In Figure 4.31, snaps of the temperature evolution during reaction are complemented by the evolution of reaction degree. The chemical reaction releasing heat induces an increase of part temperature higher than the mold temperature set at 420 K. As pictured in the snapshots of Figure 4.31, the maximum temperature found in the part reached the temperature of 500 K (230 o C). It is shown in three cut planes: two in the transverse direction and one along the largest dimension of the part. The mesh adaptation is also plotted.

We notice that the reaction starts from the hottest point defined as the regions where the mold is covering most of the SMC material. That correspond to the end of the ribs or well the regions with small thicknesses. We notice that the reaction is not homogeneous along the length and is localized also within the thickness. Notice the temperature scales up to 500 K.

Consequently, the heat released warms up the part and the inactive parts start to react in their turn. at about 40 s most of the materials is consolidated in rigid state and

Appendices

Appendix A variation of internal energy for reacting compressible materials

For the modeling of SMC or any compressible materials that undergoes a chemical reaction the internal energy is function of the temperature T , the specific volume υ and the chemical reaction α. We write then following a similar procedure depicted in [1]:

We will expressed these two term as function of the state variables p, T , υ.

-We use now the state equation that links the variation of the specified volume with the pressure and the temperature. and Maxwell relations, we get:

the term ρ de dt from the energy equation using the continuity equation, the relation v = 1/ρ and introducing the dilatation term χ T is then written as:

this represents the heat equation used in this work since the material are not supposed to be incompressibles.

Appendix B Friction Modelling

In this section, we present a numerical strategy to take into account for friction during compression molding, by considering a viscous thin layer between the two bodies in contact, whose objective is to modify the slip contact condition. These first studies on this topic need to be improved for more complex compression molding geometries, but they have provided a good framework to take into account for friction phenomena in Rem3DÂő numerical basis context.

B.1 Relation of viscous layer with Friction theories

Accounting friction using a viscous film is inspired on the boundary layer approach (See On this incompressible layer, the shearing transmitted from Body 1 to Body 2 will cause a jump in the velocity of both bodies. The shear rate γ is expressed according to:

We recall that (v 1 -v 2 ) is exactly the definition of relative velocity v r . The strain rate tensor in the layer writes: The problem in discussion can be solved analytically. Such solution is as presented in eq.B.8, where δ states for the slip condition in the layer:

The factor δ equals infinite or zero for the slip and no-slip cases respectively.

Appendix C Thermal Resistance modelling C.1 Weak formulation

Let us start by writing the weak formulation on each domain. We impose the temperatures T 1 and T 2 in the boundary and we note T C 1 and T C 2 as the contact temperature respectively:

In this way, one gets the weak formulation in the form:

We search for a solution in the whole domain Ω.

We have Ω = Ω 1 ∪ Ω 2 et n 1 = -n 2 , obtaing the following weak formulation:

Analitical solution 1D case

We seek to determine the value of the conductivity K e of the added layer. Denoting T e as the temperature in the layer, we ensure the continuity of the heat flux in the body by the expression:

After computation we get the following expression:

Sensibility analysis on thermal characterization of SMC light profile

A sensibility analysis on thermal properties is explained in the Appendix. The main information is found in the text in Chapter 4. Here; you can find the plots of the evolution of Temperature and reaction according to the sensibility analysis and comparison to the reference value obtained with the characterized values.

Evolution of the fiber orientation and Impact on rheological response

Fibers are widely used to reinforce polymer materials and improved its mechanical properties. However, the introduction of fibers within a polymer or composites will be translated into an anisotropy behavior in mechanical and thermal response of such material. Since bundle-fibers normally differ in mechanical and thermal properties in compared to the composites.

The orientation of the fibers enable us to define characteristic resistance properties of the piece in compression molding after solidification. Normally we would like to have the fibers oriented to the main directions of the mechanical solicitations. Therefore, we need to know and follow all the fiber path during the compression process so as to solve these unknown.

In this section, we introduce the mechanics of the fiber movements and its formal formulation according to the available literature. Special attention is focused on the latest models developed during the last five years.

F.0.1 Fibre orientation description -Micro scale

This forces are the ones provided to the suspensions (fibers) caused by the surrounding fluid. It is divided by translation movement of the particle and rotation movement of the particle. Which defines a really understandable relation of the evolution of the fiber. The fiber will be described by its gravity center. The translation is obtained in the core of the fiber at the send velocity than the medium, and the vorticity of the fiber obtained from the angular velocity is given by the velocity gradient of the medium. [Jeffery,1922] has them proposed the base of the suspension models:

the coefficient λ depends on form ratio of the ellipsoid, in here the geometry of the fiber is cylindrical, so this value is close to 1.

For Industrial application ingeneral, this approach is hardly costly. Modelling thousands of fiber will required extremely power form a machine. This is due to the micro-

Anisotropic Rotary Diffusion Model for Fiber Orientation in short and longfiber Thermoplastics

The Folgar and Trucker model, widely used accounts for fiber-fiber isotropic interaction through the isotropic diffusion coefficient C I . Nevertheless, such model does not match all aspects of composites with long discontinuous charges. Therefore, this models focus on takes into account an Anisotropy Rotary Diffusion (ARD), deduced form kinetics theories. Basically such anisotropy depends itself on the current orientation state. As a second improvement exhibited in this model, they use a Reduced Strain Closure (RSC) which translates to set a parameter in charge of the well prediction of the orientation evolution to fit experimental data.

Bibliography [1] Jean-François Agassant, Pierre Avenas, Jean-Philippe Sergent, Bruno Vergnes, and Michel Vincent. La mise en forme des matières plastiques. Technique & Documentation-Lavoisier, 1996.