unified approach . . . . . . . .

6 Some examples of components made using QSP ® (a) Aircraft window frame from SHERLOC QSP ® project [12] (b) A low cost composite bumper beam (c) A composite suspension arm [11] . . . . . . . . . . . . . . . . . . . . . . . 1.7 A schematic representation of sub-steps in QSP ® production line . . . . . . . . 1.8 Reproduced Fig. 3 from [14] demonstrating the overview of the literature in the field of simulation of composite forming process upto the year 2018 . . . 1.9 Experimental evidence of existence of resin-rich layer (a) A resin-rich layer at the interface of a [0/90] laminate [15] (b) Resin-rich layers at each interface of a specimen made of IMA-M21 [16]

Introduction

A composite is a type of material that consists of a blend of at least two constituents. The benefit of a composite is that, it can leverage individual benefits of each constituent material to optimize certain desirable functions and reduce the undesirable properties. Even though composites such as wood can be naturally found, man-made engineering composites are developed with the objective to optimize performance, efficiency and production cost.

The usage of composites parts in modern day industry has been increasing significantly mainly due to their high strength to weight ratio which makes them good candidates for structural components. This was one of the main driving factors in the increasing usage of composite parts in the aeronautics domain where even a small weight reduction results in a significant increase in the efficiency of an aircraft (Fig. 1.1a). With the development of various new manufacturing processes that are suitable for mass production, composites have started replacing various components in the automobile sector (Fig. 1.1b). Aerospace domain has started taking advantage of using composites in various parts. Just recently in February of 2021, NASA's Perseverance rover landed on Mars. It has carried the helicopter Ingenuity which has its blades made up of composites (Fig. 1.1c). Numerous other examples of applications of composites in various industries can be found in [1].

One form of engineering polymer composites contains at least two main ingredients: fibres and an organic matrix (Fig. 1.2a). Fibres are generally made up of glass, carbon, Kevlar or vegetal whereas the matrix can be a thermoplastic or a thermoset polymer. Typical fibres have a diameter of the order of several micrometers but have high tensile stiffness and strength. Thus, they are basically the load carrying members in a composite and provide strength. On the other hand, the polymer matrix acts as a binding agent for fibres. It provides a resistance in compression and against chemical agents. Moreover, it provides a good surface finish increasing the durability of the composite.

Thermosets are monomers with low viscosity which get converted into a three dimensional cross linked structure using heat in a chemical process called as polymerization (also referred [4] as curing). Even though a lower viscosity (before curing) can be helpful in the processing of thermoset composites, thermoset polymers cannot be remolded by reheating. On the other hand, thermoplastics contain long chains of polymers which are held together by weak secondary bonds. Because of this, thermoplastics can be reheated and remolded with some degree of degradation after each reheating-remolding cycle. Therefore, recyclability is one of the important advantages that thermoplastic composite possess which has increased their demand in the domains that involve mass production. Despite this advantage, thermoplastic polymers pose some challenges in their processing because of the high viscosity and high processing temperatures.

Fibre architecture of composites plays a key role both in the design of the composite parts as well as in the processes that can be used for their manufacturing. Some examples of fibre architectures such as unidirectional, braided, woven, non-crimp fabrics, chopped strand mat and continuous filament mat are shown in Fig. 1.2b. An overview of various manufacturing processes for both thermoset and thermoplastic composites can be found in [7]. The focus of this work however is on a specific type of manufacturing process: thermoplastic composite forming process.

Forming process for metals has been in practice for many years. The core idea behind forming is to start with a two dimensional sheet of material and then with the application of heat and/or pressure deform it to form a three dimensional shape. Forming process for composites has been derived from the same principle as that of metal forming except here, one uses a two dimensional stack of composite plies instead of metal sheets. Since it is difficult to impregnate thermoplastic polymer in a final three dimensional complex shape due to the high viscosity of the polymer, thermoplastic composite forming process often uses prepregs which are basically sheets containing fibres and a pre-impregnated thermoplastic (sometimes called organosheet).

Forming of ThermoPlastic Composites (TPC) is a process which is well suited for mass production due to its short cycle time of the order of one minute, ease of recycling of thermoplastics and its ability to create components with complex shapes. An overview of various aspects of composite forming can be found in [8,9]. A typical thermoplastic composite forming process begins with a heating stage (Fig. 1.3a) where the thermoplastic prepreg stack is heated above the melting temperature of the polymer. It is then transferred to the mold where it undergoes a forming stage (Fig. 1.3b) in which the stack deforms to take the shape of the mold under the application of force from the punch. To be more comprehensive, one could also refer to the diaphragm forming process where the punch is actually a membrane. Thus, in general the punch could be made up of a stiff material or a more compliant material. This is then followed by a consolidation stage (Fig. 1.3c) where the stack undergoes a through-thickness compressive loading due to the applied pressure on the punch. This stage serves the purpose of healing interfaces between plies and removal of porosities in the stack which is absolutely essential from the point of view of the quality and performance of the component [10]. Finally, the component is cooled down and is then removed from the mold in the demolding stage (Fig. 1.3d).

Quilted Stratum Process: QSP ®

A modified version of the standard thermoplastic composite forming process called "Quilted Stratum Process" (QSP ® ) [11] has been developed in 2012 by CETIM, Centrale Nantes and other partners. The main idea behind this process is to use UD/woven thermoplastic prepreg patches instead of using uniformly shaped prepreg stack as is the case with standard The main advantage of QSP ® over the classical thermostamping process is the use of nearrectangular prepreg patches. A demonstrative example is shown in Fig. 1.5. Due to the possibility to use smaller, near-rectangular patches in QSP ® , one can use the material more efficiently and produce more components compared to the classical composite forming process. This reduction of material wastage not only decreases the production cost but also reduces the environmental impact.

The core motivation behind QSP ® is to start with the basic raw materials and then produce a final product with maximum amount of automation and design freedom. Thus, despite the complexity involving precise patch cutting and stacking, the total time taken for QSP ® is still of the order of about a minute. All the sub-steps of QSP ® are structured in the form of a production line. This makes the process highly desirable for industries that need mass production. The choice of the patch material, its shape and position are determined strategically based on the loading to be experienced by the component. This therefore results in an optimum performance of the part and also a further reduction of material cost. Besides these main advantages, other highlights of QSP ® consist of its ability to produce multi-material, multithickness components and the possibility to perform one-shot assemblies.

Due to the numerous advantages of QSP ® , the applications for QSP ® are growing rapidly. The first stage of QSP ® consists of a continuous pultrusion process generating a cost effective in-line supply of pultruded UD tapes. These are then fed into the patch cutting system, where the prepreg patches of desired shapes and sizes are cut as per the design specifications. The cutting process occurs at a speed of about 600 mm/s. Also, based on the requirements, patches of plies obtained from other material suppliers (UD/Woven) are cut. After this, a two dimensional stack is created with the help of a robotic arm which places the patches at pre-specified locations. During this, the prepreg patches are spot welded using local heating. After this stage, the stack undergoes the standard process of heating, forming, consolidation as shown in Fig. 1.3. These stages are often coupled with certain one-shot operations such as adding metal inserts into the component [11] and an overmolding stage before demolding in order to obtain a net-shape final part. Overmolding can also have a structural and functional role to play in the overall design.

Challenges in the numerical simulation of QSP ®

Numerical simulation plays a key role in a design of any manufacturing process by removing the need to do experimental try-and-error tests. It therefore serves the purpose of shortening the process design cycle both in terms of cost and time. Also, with the availability of realistic simulation tools, process parameters can be optimized virtually, using a computer. Thus, a good numerical simulation needs to be as close as possible to the real physics of the process. At the same time, it should be fast enough to have low impact on overall process design time. Thus, one needs to find a good balance between speed and accuracy while developing a numerical solution.

Driven by the anisotropic nature of fabric coupled with its interaction with the polymer in case of prepregs, numerical simulation of a standard composite forming process itself poses several challenges. An overview about the modeling of composite sheet forming can be found in [13] and more recently in [14]. It is evident from Fig. 1.8 that a major focus of the research in composite forming simulation has been on the modeling of in-plane shear and tension in order to predict the changes in the fibre orientation during forming and the prediction of wrinkle formation. In the case of numerical process simulation of QSP ® , it has some inherent additional challenges compared to a simulation of a standard thermostamping process. They arise due to the unique nature of the process and need consideration. They are as listed as follows:

• In case of standard composite forming process, the blank holders mainly serve the purpose of keeping the plies together during forming and also help in controlling the formation of wrinkles. However, the use of prepreg patches in case of QSP ® makes it impossible to use blank holders.

• In general, a thermoplastic composite forming process involves high temperatures.

In addition to that, QSP ® involves an overmolding process at the end. Thus, during forming and a majority of consolidation phase of QSP ® , the polymer is in melt state due to the high temperature of the process necessary to perform the overmolding process later (Fig. 1.10).

• A resin-rich layer exists at the ply-ply interface (Fig. 1.9) which plays a key role in the deformation of plies [17,18], interply slip. In case of QSP ® , it can result in a possible large sliding and/or rotation of prepreg patches [19]. • QSP ® has the capability to create multi-material and multi-thickness components. This however poses a challenges in terms of design of the stack (and patch positions) for QSP ® . Challenges associated with the design phase of QSP ® known as the "Quilted Stratum Design" (QSD ® ) are discussed in [20].

• The failure and damage mechanisms for the components manufactured using QSP ® need to be studied in order to design the prepreg stack for QSP ® [21].

Mechanisms of interest

In order to design the actual component to be manufactured with QSP ® , one needs to design the two dimensional prepreg stack which after QSP ® would produce the final component. Therefore, a sufficiently good understanding of the mechanisms and physics occurring during the process is needed to obtain the desired final component.

A typical configuration of the prepreg stack for QSP ® can be represented by different layers as shown in Fig. 1.11. A stack shown here consists of three plies. Each of these plies can be further divided into two types of zones. The first type of zone is a composite layer which consists of fibres (UD/Woven) along with a preimpregnated polymer. The second type is a resin-rich layer which predominantly consists of the polymer.

In general the mechanisms discussed below can be classified as intra-ply or inter-ply, inplane or out-of-plane and more importantly as occurring during forming or consolidation.

A distinction between forming and consolidation can be made by an existence of a compressive normal stress and/or reduction of thickness. Also to be noted is the fact that these mechanisms are not necessarily independent of each other and some of them can occur simultaneously. • In-plane loading: This intra-ply phenomenon is the response of a composite layer under the in-plane membrane type of loading. Therefore, this is important in both forming and consolidation phase. This can be further classified into tension/compression/shear (Fig. 1.12a and Fig. 1.12b). More discussion can be found in [8,9].

• Out of plane bending: This is an intra-ply phenomenon which mainly occurs during the forming phase. However, due to the ply-ply interaction, the bending effects are transferred from one ply to the next ply in the stack (Fig. 1.12c). More details about its modeling and characterization can be found in [23][24][25][26].

• Intra/inter-ply slip: This deformation mode occurs due to transverse shear (Fig. 1.12d) [22,27]. An experimental evidence of this mechanism is demonstrated in Fig. 1.13. Since QSP ® typically involves a stack made up of thin plies, the interply slip plays a bigger role compared to the intra-ply slip. This is further discussed in detail in Chapter 2.

• Sliding and rotation of prepreg patches: This is an interply phenomenon which falls under the broad category of interply-slip (Fig. 1.14a and Fig. 1.14b). In QSP ® however, the melt state of the polymer at the interply interface, use of prepreg patches and the inability to use blank holders can result in large sliding (maximum patch sliding reported by [19] is 24.5 mm) and rotation of prepreg patches.

• Interply adhesion: This is an interply phenomenon which mainly occurs during the forming phase. This is due to the tack of melt polymer at interface which provides • Squeeze flow: This is an intra-ply phenomenon which occurs during the consolidation phase. In this mechanism, due to the consolidation pressure, the thickness of the composite layer reduces, creating a resin movement within the ply which in turn convects the fibres with them (provided that there is a room for the material flow in the plane) and modifies the fibre orientations [17,18]. This is schematically represented in Fig. 1.15. Evidence of this mechanism can be observed in Fig. 1.16 in case of an experiment and also on a real component manufactured using QSP ® . This deformation mode is discussed in more detail in Chapter 4.

• Non uniform consolidation/squeeze: This mainly occurs due to the non-uniform distribution of the consolidation pressure (and subsequently non-uniform squeeze flow) which can deform both the composite layers as well as resin-rich layers (Fig. 1.17a). This deformation can affect the final thickness distribution of the composite layer subsequently affecting the final material configuration (Fig. 1.17b) and can result in wrinkle formations during consolidation [28]. More discussion about this deformation

(b) (a)
Figure 1.17: (a) Schematic representation of a non-uniform consolidation/squeeze (b) Experimental evidence of non-uniform consolidation/squeeze [16] mode can be found in Chapter 3-5. Apart from these primary mechanisms, there are other phenomena which could occur during the consolidation. These phenomena can be clearly distinguished from the other consolidation mechanisms based on the fact that these involve a relative motion between fibres and polymer whereas in case of a squeeze flow; both are treated as a single continuum object generally modeled as anisotropic viscous fluid [29].

• In-plane resin bleed: This phenomenon can be defined by exiting of the polymer from the ply (in the plane of ply) due to the applied consolidation pressure (Fig. 1.18 a and Fig. 1.18 b). The resin can bleed both from resin-rich layer as well as from the impregnated fibrous layer. The direction and amount of resin bleed is dependent on the viscosity of the polymer, permeability of impregnated fibrous layer and the configuration of the ply. An experimental evidence of a longitudinal resin-bleed in case of quasi-UD glass fabric [UDT400P] 6 saturated with silicone oils during different stages of the consolidation process can be found in [30].

• Transverse resin percolation: This phenomenon can be defined as the resin movement out of the ply (normal to the plane of ply) due to the applied consolidation pressure. It was hypothesized that this mechanism is probably responsible for the 1.3. Objectives of the thesis formation of the resin-rich layers between plies [31,32].

However, considering the fact that QSP ® and other similar composite forming processes use a highly viscous polymer along with high fibre volume fractions, the simulation of in-plane resin bleed and transverse resin percolation are postulated to be of secondary importance.

In summary, for the standard thermoplastic composites, mechanisms such as in-plane shear and tension are relatively well established in the literature both in terms of simulation as well as their characterization. On the other hand, mechanisms such as interply adhesion, transverse squeeze flow, simulation of consolidation are not yet standardized even though their influence on the final deformed shape is evident.

As for the process simulation of QSP ® , it shares many features with the standard composite thermostamping process but at the same time there are some additional challenges both in terms of actual forming process as well as its numerical modeling. Thus, in order to have a full-scale industrial solution for numerical process simulation of QSP ® , the mechanisms of interply adhesion occurring during forming and transverse squeeze flow occurring during consolidation need to be incorporated into the overall simulation framework.

Objectives of the thesis

The objective of the thesis originates from the goal of the numerical process simulation of forming and consolidation phases of QSP ® . The goal of the numerical process simulation of QSP ® is to predict the following quantities that are of interest from the design point of view:

• Final positions and orientations of prepreg plies and patches

• Deformations of plies

• Fibre orientations within plies

To the best of our knowledge, currently there is no commercial software solution that is able to tackle this problem of performing the numerical process simulation of both forming and consolidation phases of QSP ® in a unified manner using a shell element based approach.

Altair Radioss ™ code is a highly suitable candidate for the full scale simulation of forming and consolidation phases of QSP ® . Altair Radioss ™ is on top of the industrial FE codes which are used for crash simulation as well as for fluid-structure interaction or stamping process simulation. Due to the in-built rich database of constitutive material models, element formulations, various contact algorithms, scalability, speed and its integration with other Altair products such as Hypermesh; it serves as a good tool to build upon. Also, implementation in the industrial code of Altair Radioss ™ is a request made by our industrial partner CETIM.

The overall goal of this thesis is to develop the capabilities of the numerical process simulation of forming and consolidation phases of QSP ® in the commercial code of Altair Radioss ™ . More specifically, the developments in this thesis are subdivided into the following tasks:

• Forming phase:

-Development of a semi-empirical contact mechanism to model interply adhesion in the contacts.

• Consolidation phase:

-Development and implementation of an extension of a classical shell element with selective mass scaling in order to have a capability of a transverse normal stress (σ zz ).

-Development of a constitutive model for nearly incompressible melt polymer that is compatible with this modified shell.

• Unified simulation of forming and consolidation phases:

-Combining the developments in forming and consolidation phases in order to perform a full scale numerical simulation of QSP ® in a unified manner.

Outline of the thesis

The developments done in this thesis are subdivided into several chapters as follows:

Chapter 2 is focused on the forming phase of QSP ® and discusses the importance of modeling interply adhesion. A novel semi-empirical contact mechanism has been developed to model this phenomenon that creates a finite adhesive strength for ply-ply interaction during forming. Further, an experimental procedure to characterize this contact mechanism has been proposed. The chapter concludes with a full-scale industrial validation case where the final position, orientation of plies as well as fibre orientations at several locations were compared with the experimental observations.

Chapter 3 is dedicated to the consolidation phase of QSP ® . It discusses the development of an extension of the hourglass-free full-integration QBAT shell element in Altair Radioss ™ which is in turn is based on the Q4γ24 shell element discussed in [33,34] in order to have the capability of transverse normal stress for the shell. This is achieved by augmenting the QBAT element with additional nodal degrees of freedom. This is then followed by several developments in the contact mechanism to create compatibility with this modified shell element. The chapter concludes with a development of a selective mass scaling technique for this shell which makes it possible to use a much higher time-step for the simulation.

Chapter 4 focuses on the constitutive modeling of nearly incompressible melt polymer which is an essential ingredient in modeling the mechanisms during consolidation such as transverse squeeze flow and its influence on the changes of the fibre orientations in the prepreg. This chapter contains a discussion on the development of an elasto-plastic constitutive model that is compatible with the modified shell element. For a physical comparison with respect to the consolidation process, a practical and economic characterization method is proposed. Finally, the consolidation test results on a [0 • /90 • ] 2 are compared with the simulation results.

Chapter 5 combines the developments from all the previous chapters in order to do a fullscale numerical simulation of QSP ® consisting of both forming and consolidation phases in a unified manner using a shell element based approach. The chapter concludes with a simulation on semi-industrial and industrial models where the results of the simulation are compared with the experimental observations.

Finally, Chapter 6 provides a conclusion of the work done in this thesis along with a discussion about possible scope of future work and perspectives.

Chapter 2

Numerical modeling of interply adhesion in composite forming of viscous discontinuous thermoplastic prepregs

Abstract

This chapter focuses on the development of a semi-empirical contact mechanism that has been developed to model the phenomenon of interply adhesion during the forming phase of QSP ® . An experimental procedure to characterize this contact mechanism is also discussed. Finally, a validation on a full-scale industrial test case is performed where both the ply positions and fibre orientations are compared with those obtained from experiments.

Introduction

Interply adhesion arising due to polymer tack plays an important role in QSP ® due to the presence of resin-rich layer at the ply-ply interface, inability to use blank holders for prepreg patches and high temperature of forming process where the polymer is in melt state. Without modeling interply adhesion in the numerical simulation of their forming, plies delaminate unrealistically without any resistance which in turn results in incorrect final positions of the prepreg patches.

Interply Adhesion and QSP ®

Interply adhesion arises from the phenomenon known as Polymer Tack which occurs due to the presence of molten polymer layer at the ply-ply interface. This phenomenon provides a finite resistance to the interply delamination. The importance of modeling interply adhesion is even more for a process like QSP ® due to the absence of blank holder which would otherwise help in keeping the plies together. It was observed that in the numerical simulation of QSP ® , if interply tack is not modeled, the stack of plies unrealistically delaminates without any resistance (Fig. 2.16a, discussed in detail in Sec. 2.6.1). The phrase "unrealistic delamination" has been used here since, our industrial partner CETIM did not observe such delamination while the preform is formed. This in turn results in incorrect prediction of the final positions of the prepreg patches.

Nomenclature wise, the cohesive and adhesive failure can be clearly defined when the system consists of a substrate and an adhesive with a clearly distinguishable boundary between the two. However, in the case of interface between composite plies, it is a little ambiguous due to the presence of polymer in a continuous manner. At macro level, it can be called as interply adhesive failure (between plies) even though on the micro level, it is a cohesive failure (more details are given in Sec. 2.4.1).

Polymer Tack

The term tackiness is mainly defined in reference to adhesives in liquid state. When one tries to transversely separate two solids connected by an adhesive in melt state, tackiness is the resistance offered by the adhesive to the act of normal (orthogonal) separation [1]. The general theory of tackiness for polymers is discussed in [2,3]. Even though adhesion can originate from various mechanisms such as, mechanical, diffusion, electrostatic and physical/chemical absorption, in regard to the polymer adhesives, the diffusion theory is more relevant [4]. The core idea is that a polymer macromolecule attached to one surface gets diffused into the second surface, creating adhesion between two surfaces.

As one tries to separate the adhesive joint, various mechanisms take place at the interfaces. A detailed experimental analysis was conducted in [5,6] to study these mechanisms for Pressure Sensitive Adhesives (PSA). A typical stress-strain curve obtained in [5] is shown in Fig. 2.1. Two main parameters quantifying the adhesion were considered, adhesive strength (maximum adhesive stress) and adhesion energy. The debonding process of an adhesive joint consists of several stages, that differs with the debonding rate applied to the sample, as well detailed by Poivet et al. [7] for silicon oils. The rupture can be due to the polymer flow and fingering (low debonding rate) or to cavitation (high debonding rate). The experimental evidence of both has been shown in Fig. 2.9b and Fig. 2.9c. In both cases, there is a maximal adhesive measurable force that the sample can handle.

In general, the polymer tack is characterized by measuring the adhesive force response during the debonding phase. A measurement method for tack of polymers was proposed in [8] which involved a layer of polymer deposited on a rigid plate. The flat end of a cylindrical shaped probe made of stainless steel was first brought in contact with this polymer layer. After that, as the probe was pulled in opposite direction (at a controlled debonding rate), the force response was recorded. An ASTM standard for this type of test was proposed in [9] known as the 'Probe Tack Test'. This technique was used to study the low and high debonding velocity regimes of adhesive failure [7]. It was also used in [10] to measure the force-separation relationship for M21 epoxy resin which was used further to study the local buckling of slit tape for Automated Fibre Placement technique. More recently, this was used by [11] to study the dependence of tack on various parameters for Automated Tape Laying (ATL) process.

A peel-type test was also proposed to measure the tack of prepregs for Automated Tape Laying process in [12,13]. This test was employed to study the dependence of tack on different parameters in [14]. A floating roller peel test was employed in [15] to measure the polymer tack against the level of cure (Fig. 2.2). Measurement using peel-type technique is more challenging since it involves correct characterization of other material parameters of prepreg such as bending rigidity which will have a strong influence on the debonding force response.

Several studies have been carried out to measure the adhesion for thermosets (especially on epoxy based) prepregs. The effect of aging and temperature on tack was studied in [16].

A viscoelastic model which requires 4 material parameters was proposed in [17] which was then used to study the tack of epoxy/T-300 carbon fibre reinforced prepregs. For the case of UD-CF/epoxy based system, [18] used Design of Experiments (DOE) to study the effect of impregnation temperature and pressure on tack. The influence of contact time, contact force, temperature, debonding rate and aging on tack was studied in [19]. The adhesion force variation with respect to displacement for the prepreg was compared with the behavior of a pure polymer. It was discussed that the behavior difference is seen due to the structural effects in case of prepregs, which mainly arise due to the presence of fibres and non-uniform thickness of polymer. Various studies based on DOE were conducted to study the dependency of prepreg tack on various parameters [20,21].

When it comes to tack of thermoplastic melt prepregs, there is no standard test reported in the literature. As a consequence, the interply adhesion created during the QSP ® forming phase is not taken into account in QSP ® numerical simulations of thermoplastic prepregs. It is therefore of great interest to propose a novel approach that encompasses numerical and experimental studies to evaluate the influence of the interply adhesion mechanisms on the final spatial positions and orientations of prepreg plies in a full-scale industrial part.

Objectives and content of the study

The main objectives of the study are fourfold:

• Propose a simple semi-empirical model that reproduces the interply adhesion between two molten prepreg plies that plays role during the forming phase of the QSP ® process.

• Quantify experimentally the adhesive strength of a molten polymer to fuel the interply adhesion model.

• Implement the proposed model coupled with experimental characterization in a fullscale numerical simulation of QSP ® to compute the prepreg patch spatial positions and orientations.

• Carry out an extensive comparison between the model's predictions and measurements on a representative industrial part.

After a general presentation of the full-scale simulation strategy, a novel semi-empirical model for prepreg interply adhesion is detailed. Its intentional simplicity requires the characterization of a single parameter which is named adhesive strength. Therefore, in a second part, the experimental setup and methodology developed to quantify the adhesive strength of a molten polymer is presented. Thus, the novel model combined with the experimental data are both integrated in the full-scale numerical simulation of QSP ® using Altair Radioss ™ . The macroscopic numerical results concerning the plies spatial positions and orientations are compared to the experimental results obtained on a full-scale industrial part.

Numerical modeling of interply adhesion

Numerical simulation of thermoplastic composite forming process plays a key role in a design of the prepreg stack to be used for forming. A good simulation of a process can avoid the experimental trial and error tests. It therefore has a potential of shortening the process design cycle both in terms of cost and time. Also, with the availability of realistic simulation tools, process parameters can be optimized virtually, using a computer. Thus, a good full scale numerical simulation should be able to capture the physics of the process. At the same time, it should be fast enough to have low impact on overall design time.

For the numerical simulation of composite forming process, a wide range of approaches exist. A detailed review of these approaches can be found in [22]. Despite the inherent multi-scale nature of the composites, a choice has been made in this work to remain within the scale of macro-level. This choice is made because of the complexity of the components involved in QSP ® and computational performance. The main objective of the numerical simulation of QSP ® consists of predicting the final position of the discontinuous plies and the fibre orientations in individual plies.

The focus of this work is on the forming stage of QSP ® . In simulation, each individual ply is modeled with shell elements at the mid-surface of the ply. The main idea in this work is to replace the resin-rich layer at the interface with an equivalent contact law. This modeling choice has been motivated by the following factors:

• The explicit modeling of resin-rich layer results in very expensive computation [23] making it challenging to be used in full scale industrial models.

• Creating such an interface law instead of explicit modeling of the polymer would benefit in capturing the possible large sliding/rotations of prepreg patches.

• During the forming stage, it can be considered that there is no significant consolidation pressure. Thus, the mechanisms such as bleeding/percolation/squeeze-flow do not occur, keeping the resin-rich layer thickness nearly the same. A numerical model dedicated towards the consolidation phase of the composite forming process is discussed in [24].

There have been several different approaches in the literature aimed at numerical modeling of interply adhesion. One possibility of using cohesive zone elements was discussed in [25]. However, the possible large sliding/rotations make this approach not suitable for QSP ® simulation. Interface adhesion limited by a deactivation distance was considered in [26].

In their work, an educated guess has been used to assign the values for the adhesion strength (0.1 MPa) and the deactivation distance (2.0 mm) where the ply thickness was 0.14 mm. This was done due to the lack of availability of a standardized test for characterization of adhesion for thermoplastic melt prepregs.

Dörr et al. [27,28] modeled interply adhesion with adhesive stiffness and the interply adhesive forces are applied as long as the normal separation between the plies is less than a certain value. The calibration method and the values used for adhesive strength and the deactivation distance are not specified in this work.

Based on the literature review, an interface mechanism providing adhesive strength while permitting large sliding/rotations of plies does not exist. Also, a clear and precise method to characterize such adhesive interface models is needed. Thus, a penalty based, semiempirical contact model for interply adhesion has been developed and implemented in the industrial finite element code of Altair Radioss ™ which is discussed in the next section.

Contact mechanism for interply adhesion in normal direction

Before getting into the details of the development of the contact mechanism capable of adhesion, it is essential to state the requirements of such a mechanism for it to be suitable for QSP ® .

• The model should require the least number of parameters which need characterization. This would make it more suitable for an industry. At the same time, it should have sufficient parameters to capture the physics of the mechanism as closely as possible.

• The model should allow large sliding/rotation of plies while retaining adhesion which is essential for QSP ® .

• The model should provide a finite adhesive strength and an ability to delaminate if the normal stress exceeds this value. It should be able to automatically recreate the adhesive bonds if the separated plies come in contact again in future.

• From a point of view of computational performance, the mechanism should be scalable and should work with both Shared Memory Processors (SMP) and domain decomposition with Message Passing Interface (MPI).

As discussed before, the adhesive strength in general depends on many parameters which include, contact time, contact pressure during contact time, debonding rate, temperature, surface characteristics, thickness and nature of the adhesive. However, it is discussed in [29] that as the contact time and contact pressure increase, the adhesive strength becomes almost independent of these parameters. A similar observation was made for debonding rate. When it increases, the adhesive strength asymptotically approaches a constant value. In another work, it was discussed that the effect of surface roughness diminishes as the contact forces in bonding become high [30]. Based on these observations, a contact algorithm characterized solely by the adhesive strength (σ adh ) and linearly increasing force-separation relationship is proposed here.

The modeling strategy for QSP ® in Altair Radioss ™ is to model each individual ply in a prepreg stack (Fig. 2.3a) with shell elements created at the mid-surface of the solid and to define a standard penalty based node-to-surface contact between the plies (Fig. 2.3b). As an example, a zoomed in view of the actual component model in Altair Radioss ™ consisting of 9 plies is shown in Fig. 2.3c. In the penalty based contact algorithm defined for the interaction between plies (Fig. 2.5a), a penetration resistive spring is added whenever the secondary node enters the penetration zone. The thickness of penetration zone is given by,

Gap = L = 0.5 × (t main + t secondary ) (2.1)
where t main and t secondary are the thicknesses of main surface and secondary node respectively. The representation of the case without penetration and the case with penetration is shown in Fig. 2.4. This standard penalty based contact formulation is modified to have an additional zone of thickness L which extends beyond the penetration zone, called as adhesion zone (Fig. 2.5b).

The core idea is that, as soon as the secondary node touches/enters into the penetration zone, an adhesive spring is created. This spring has no effect in the penetration zone. Thus, the mechanism behaves the same as that of the standard contact algorithm inside penetration zone. However, as the secondary node tries to move away from the main segment and comes in the adhesion zone, adhesive spring (whose base is connected at the interface of adhesion and penetration zone) exerts adhesive forces pulling main and secondary pair towards each other. A simple linear force-displacement relationship is used which is given by Eq. 2.2,

F adhesion = σ adh A L (L -P ) (2.2)
Where, σ adh is the adhesive strength, A is the contact area, P is the normal distance measured with reference to the outer boundary of the adhesion zone. As the node is further pulled and it reaches the outer boundary of adhesion zone, adhesion spring ruptures representing delamination. The adhesive spring is recreated in case the secondary node enters the penetration zone again.

Mechanism verification with a simple test case

Before using this adhesion mechanism on full scale industrial scenarios, it is tested on a simpler simulation scenario consisting of two flat plates each of thickness 0.5 mm. The plates are modeled with 4-node shell elements created at their respective mid-surfaces (Fig. 2.6a).

The bottom plate (Main) is fixed in all degrees of freedom whereas the top plate (Secondary) is subjected to imposed time varying displacement in the normal (z) direction such that the normal interply distance is given by the curve shown in Fig. 2.6b. The corresponding responses of contact force and contact energy obtained from the simulation are shown in Fig. 2.7. Note that the values of adhesive strength used in this plot are not real values. In reality, the adhesive stiffness is several orders less compared to the penalty stiffness. Also it is to be noted that this contact mechanism does not require the main and secondary nodes to be co-incident, the contact is of type node-to-surface where the penetration on the secondary node is calculated based on its projection onto the main surface details of which can be found in [31]. During the path A-B, the plates experience penetration, generating a penetration resistive penalty force which tries to separate the main segment and the secondary node. During this, since the secondary node has entered the penetration zone, adhesion spring is created which creates adhesive forces as the secondary node enters into adhesion zone (after point B). During path B-C-D-E-F, this adhesion spring undergoes cycles of stretching and relaxation. The adhesive force which tries to pull main and secondary towards one another, possesses an opposite sign as that of penetration resistive forces which can be observed in Fig. 2.7. At point G, the adhesion spring is fully stretched (exhibiting maximum adhesion force) and beyond that, it breaks. It is important to note that after point H, even though the top plate enters adhesion zone, the adhesive spring is not recreated until it reaches the penetration zone again (point I). Path J-K demonstrates the second cycle of increasing adhesion and a rupture at the end. The dissipation resulting from each rupture of the adhesion spring is stored irreversibly in the contact energy which gets accumulated with every cycle of adhesion-rupture.

Interply contact mechanism in tangential direction

For the full scale implementation, a complete representation of the interface mechanisms is needed. Therefore, in addition to the transverse interply adhesion, a model for tangential interaction between the plies is needed. During the forming phase of QSP ® , the plies do not experience a significant consolidation pressure. Thus for the interply contact behavior in tangential direction, a viscous slip law discussed in [32] was chosen over the Stribeck curve based approach discussed in [33]. Therefore, a tangential force is applied which acts in the opposite direction of the relative sliding velocity of the secondary node with respect to the main surface. This tangential force is applied irrespective of the zone where a secondary node is lying. The magnitude of this force is given by Eq. 2.3,

F T = S η A L (V rel ) (2.3)
Where, η being the viscosity of melt thermoplastic polymer, A is the contact area and V rel is the tangential relative speed of the secondary node with respect to its main surface. In many full scale industrial models consisting of multiple plies with different ply architectures and fibre orientations, it is difficult to measure the exact thickness (and variation of thickness) of resin-rich layer. Therefore, a normalized length L is used in Eq. 2.3. However, a scaling factor S is provided which can be used to scale-up/scale-down the tangential viscous forces based on the available information about the resin-rich layer thickness.

Overall, this mechanism satisfies all the needed requirements described in Sec. 2.3.1. However, the mechanism modeled here takes certain computational liberties which were done mainly in order to reduce the characterization needed for modeling the interply adhesion mechanism. They are:

• The post-peak softening behavior of the force-separation behavior is not modeled.

• The thermal effects are not modeled.

• The adhesion zone is assumed to be extending up to a normal distance of L.

The assumption of adhesion zone thickness to be L is done for two reasons. First, it reduces the needed characterization and in general it is difficult to measure the thickness of this resin-rich zone which can vary throughout the specimen as shown in Fig. 2.8 (reproduced from Fig. 4.12 in [34]). Second, if the adhesion zone thickness is very small, it will increase contact stiffness significantly (defined by k adh = σmaxA L ) which needs to be taken into account for the stability of the time-step calculation. Therefore, this would result in a very small time-step (which is inversely proportional to the stiffness) for the explicit solver, increasing the computational cost significantly.

Experimental determination of the interply adhesion strength

From the characterization point of view, the adhesion model proposed here needs only one main parameter, the adhesive strength (σ adh ). After a detailed literature review on polymer rupture, this section describes the developed experimental setup and methodology used to determine the adhesive strength (σ adh ) of melt thermoplastic polymer PA-66.

Literature review

Cohesive or adhesive rupture

A resin-rich layer usually exists at the interface between two plies of prepreg [35] and the rupture occurring when the transverse stress exceeds the transverse maximal strength adhesion σ max may be adhesive (rupture at the interface between the fibrous reinforcement and the polymer) or cohesive (rupture in the resin bulk). Bastien et al. [36] investigated the fusion between two graphite/PEEK prepreg plies at molten state using a thin amorphous polyetherimide (PEI) film. Authors highlighted that for a contact time t c between both components higher than the PEI reptation time τ rept (which sets the self-diffusion rate of the polymer [37]), the transverse rupture occurs in the bulk of the PEI polymer. Experiments conducted on highly entangled polymers lead to similar conclusions [38]. According to molecular dynamics simulations [39], if t c < τ rept , the polymer entanglement density

The experimental part i.e. the sections 2.4.1-2.4.3 were undertaken by Dr. J. Hemmer as a part of the post-doctoral research conducted at Centrale Nantes remains low and leads to chain pull-out (adhesive rupture); in contrast, if t c > τ rept , the polymer entanglement density is maximal and the cohesive rupture is due to chain scission.

For PA-66, a conventional polymer used to make thermoplastic prepreg, τ rept ≈ 10 -5 s (obtained following Doi Tube theory [37] and using molecular data from [40]). Before the forming process, the prepreg plies remain stacked at melt temperature during much longer time. Therefore, the transverse rupture between two prepreg plies is expected to be cohesive, i.e in the bulk of the molten polymer.

Fluid-like or solid-like rupture

The molten polymer rupture depends on the applied strain rate ε [41,42]. A low strain rate induces a fluid-like rupture (ductile) whereas a high strain rate induces a solid-like rupture (elastic). When the product ετ rept is higher than 1, the polymer macromolecules do not have enough time to reorganize following the imposed strain and the induced rupture is elastic [43,44]. Schach et al. [38] consider that the transition between fluid-like and solid-like rupture occurs for ετ rept varying from 0.8 to 7.

During the forming process, the strain rate ε ∈ [1 -1000] s -1 and again, τ rept ≈ 10 -5 s for PA-66. The cohesive rupture is therefore expected to be fluid-like between two prepreg plies at molten state.

Following the conclusions raised by the present literature review, the adhesive strength between two melt thermoplastic prepreg plies will be investigated through the tensile mechanical behavior of the corresponding molten polymer.

Experimental methods and setup

Fig. 2.9 presents the setup developed to measure the transverse tensile mechanical behavior of a molten PA-66 following the methods employed for silicon oils [7]. A Dynamic Mechanical Analysis device (DMA 150+ from Metravib) whose specifications are detailed in Tab. 2.1 is used as a tensile machine where jaws are two parallel steel platens. A PA-66 sample from Solvay (melting temperature T m =280 • C and melt viscosity 70 Pa.s) is placed on the bottom fixed platen, located in a furnace.

The initial thickness of the molten polymer is set by imposing a compressive constant force (-0.5 N) on calibrated shim plates during 1 minute. Once this bonding step achieved, a constant tensile force rate is imposed and both force and displacement are recorded. The setpoints applied during each step are detailed in Tab. 2.2. The sample volume remains constant during the experimental campaign (diameter of 33 mm and thickness of 0.5 mm) and each test is repeated twice at a given force rate.

Experimental results

Qualitative results obtained on molten PA-66 (Fig. 2.9b) are in good agreement with the one obtained on silicon oils (Fig. 2.9c, extracted from [7]). At low debonding velocity, a radially convergent flow is observed (fingering effect) whereas cavitation occurs at high debonding velocity. In QSP ® , the debonding rates of plies are quite high (discussed in detail in Section 2.4.4) which indicates that the mechanism consisting of extensional separation is more important compared to the mechanism creating radially convergent flow. Barroso et al. [41] distinguished the failure (loss of adhesion) happening when the force

Step

Force Temperature Melting 0 N 10 reaches a maximal peak and the rupture (no polymer threads remain between the two platens) occurring after the failure. This force peak is believed to represent the instability and transition from homogeneous stretching to necking [45] and the use of the Considère criterion to predict it for molten polymer and polymer solutions is still discussed [46,47].

For the present study, the maximal force is found to be as a good indicator to quantify σ max .

Even if test repeatability is difficult to ensure, the molten PA-66 tensile behavior significantly depends on the imposed force rate (Fig. 2.10). The transverse extension is measured using ∆h h 0 where ∆h is the change of thickness with referece to the initial thickness h 0 . The sample rupture, not represented here, occurs right after the force peak at high force rate and long after the force peak at low force rate. The force rate significantly impacts the maximal strength adhesion value (σ max in Fig. 2.11a). Equivalent conclusions are drawn for silicon oils [7]. However, the material viscosity (10 2 to 10 3 Pa.s in [7]) does not influence the force peak value. Therefore, the obtained trend curve (Fig. 2.11a) can be used to determine the maximal adhesive strength for a more viscous PA-66 molten polymer.

Choice of adhesive strength value for simulation

As discussed in [19], the non-uniformity of interfacial resin layer thickness in case of prepregs affects the post-peak behavior of the force-separation curve, however, the linear nature of the pre-peak behavior is still retained. Since, the focus of this work is on pre-peak linear behavior, the adhesion tests conducted on pure polymer can be used to model the interply adhesion for thermoplastic prepregs.

As observed from the experiments, the value of the adhesive strength varies with the debonding rate (Tab. 2.3). The debonding rate was calculated using ∆h h 0 ∆t and the adhesive strength was calculated with Fmax S 0 where S 0 is the original contact area of the sample. However, it was observed by [6,29] that at higher debonding rates, the adhesive strength becomes less dependent on the debonding rate. Inspired from this observation, a curve fitting of the form y = a(1exp -bx ) + c was performed (Fig. 2.11a). The parameters obtained from the curve fitting exercise were a = 0.0052, b = 67.86 and c = 0.001 which produced an R 2 = 0.9921. With higher debonding rates, the value of adhesive strength approaches asymptotically to a value of 0.0062 MPa. The value of local debonding rate was calculated for each node, each contact interface after every 1000 cycles in simulation which resulted in a database of about 2 million data points.

A histogram plot (Fig. 2.11b) for these values shows that the majority of the values occur around the debonding rate of 10 0.225 = 1.68 (1/s) for which the adhesive strength can be safely considered to be 0.0062 MPa. This justifies the use of a constant adhesive strength value = 0.0062 MPa in the simulation. The sensitivity analysis of adhesive strength has been performed and is discussed in Sec. 2.6.2.

Implementation aspects of interply adhesion in Altair Radioss ™

Interfaces solve the contact and impact conditions between two parts of a model. There are several types of contacts in Altair Radioss ™ such as node-to-surface with nonlinear contact stiffness (Type-7), node-to-surface with a constant contact stiffness (Type-25), bonded contact (Type-2), edge-to-edge contact (Type-11) etc. each of which has been designed for a specific purpose and use. The focus of this work is on the Type-25 contact which is a node-to-surface penalty type of contact with a constant penalty stiffness.

Global overview of contact algorithm

Fig. 2.12 shows a global flowchart of the Altair Radioss ™ explicit solver. Altair Radioss ™ being a commercial solver, the actual architecture is significantly more complicated. Individual steps in this flowchart are discussed in detail in Sec. 3.5.1 in Chapter 3. In this Chapter, the discussion has been restricted to the contacts and specifically Type-25 contact interface.

The block Radioss Starter reads the input file containing all the details about the model: nodal coordinates, element connectivity, different elements, contacts, loads, boundary conditions and so on. This block was modified in order to read the additional information related to the interply adhesion, creation of adhesion zones and other related data structures. The block corresponding to the Allocation and initialization was appended in order to allocate and initialize adhesion related data structures and incorporate specific treatments such as the one for excessive initial penetrations.

The main modifications have been done in the block that calculates the contact forces. This has been updated in order to calculate the adhesion related contact forces in the normal direction F adhesion contact and the viscous sliding forces in the tangential direction F tangential contact . Fig. 2.13 shows a simplified flowchart of the contact algorithm. The initial steps in a typical contact algorithm consist of contact sorting and then contact searching which are not discussed in this work (Refer to [31] for more details). Then the normal projection of a secondary node on a main segment is calculated (discussed in Sec. 2.5.2). Based on this normal distance and the 'Gap', it is identified if the secondary node is penetrating the main segment or not. If it is penetrating, then there are two types of forces that need to be calculated, (a) normal force (b) tangential force.

In the standard penalty type of contact, the normal force is directly calculated based on the penetration and contact stiffness. In our case however, the normal force could arise from (a) adhesion or (b) penetration. Thus, it needs to be determined if the secondary node resides in the adhesion or penetration zone and if the adhesive spring exists for that node before we can calculate the normal force.

In case of the tangential direction, depending on the type of law used (such as Coulomb friction, Renard law, Darmstad law etc.), the tangential force is calculated. A block has been added to incorporate the tangential viscous sliding law discussed in this Chapter.

Once the normal and tangential forces are calculated, they are applied on the secondary node and the corresponding nodes of the main segment considering the appropriate direc-tions. Finally, the global assembly of the contact force is performed in order to obtain the global contact force vector F contact .

Calculation of the secondary node projection and penetration

In order to calculate a penetration, it is necessary to calculate the normal distance between the main segment and the secondary node. For this, the projection of the secondary node onto the main segment needs to be calculated. This is shown in Fig. 2.14 for the case of a quadrilateral main segment. The objective is to find the corresponding projection P of the secondary node Q on the main segment ' ABCD'. A, B, C and D are the position vectors (in the global coordinate system) of the nodes forming the quadrilateral main segment. Using the isoparameteric coordinates ξ and η and the bilinear shape functions N 1 -N 4 (Refer to Tab. 3.2); a system of equations is obtained as follows,

P = N 1 A + N 2 B + N 3 C + N 4 D (2.4) s = (1 -η) AB + (1 + η) DC (2.5) t = (1 -ξ) DA + (1 + ξ) BC (2.6) n = s × t |s × t| (2.7) P Q = a n (2.8)
where n is the unit normal to the mid-surface at P and a is the projection distance. Solving this system of equations; one obtains the position of point P and the normal distance a.

Then the penetration is calculated using,

Pene = max(0, 2 Gap -a) (2.9)
where 'Gap' is calculated based on the thickness of the main and secondary as per Eq. 2.1 and Fig. 2.5b. With this, the corresponding penetration force or adhesion force is calculated.

Distribution of penalty forces on main segment

Knowing the unit normal at the point of projection, the contact force on the secondary node is applied directly. However, in case of the main segment, it is necessary to distribute the contact force at the projection onto the nodes of the main segment.

The contact forces F 1 , F 2 , F 3 and F 4 are applied on the nodes at A, B, C and D respectively in the opposite direction as that of the secondary node, such that:

F 1 + F 2 + F 3 + F 4 = -F secondary (2.10)
The forces are a function of the position of the point P = (ξ P , η P ) with respect to the nodes of the main segment,

F i = -N i (ξ P , η P )F secondary (2.11)
where N i (ξ P , η P ) represents the bilinear shape function i evaluated at the projection point P = (ξ P , η P ) which signifies that the main segment node closer to the projection point will experience a higher contact force and similarly, the main segment node farthest from the projection point will have the smallest contribution towards the overall contact force.

Full scale numerical simulation of QSP ® using Altair Radioss ™

This section discusses a numerical simulation of QSP ® for an industrial component performed using Altair Radioss ™ .

Model description and setup

The considered industrial model (Fig. 2.15a) has a prepreg stack consisting of 9 plies of different shapes, fibre architectures (UD/woven) and fibre orientations (Fig. 2.15b). This component has been designed by CETIM after doing several initial design iterations. The woven prepreg plies are TEPEX ® dynalite 101-RG600(x) (Supplier: BONDLAMINATES).

The UD ply tapes are Celestran ® CFR-TP PA-66 GF60-02 (Supplier: TICONA).

The punch and die are modeled with rigid elements. Each ply is meshed with 4-node shell elements (created at the mid-surfaces of plies) with hourglass free full integration elements [48]. Since, a prepreg ply consists of both fibres and the polymer, each ply is modeled by overlapping two sets of shell elements with co-incident nodes. One set of shell elements obeys a fabric constitutive model and the other set obeys the polymer constitutive model. The mesh size of 4 mm was used which generated a total 29,653 nodes consisting of 48,404 quadrilateral elements.

Constitutive model for fabric behavior has been implemented in material Law-58 of Altair Radioss ™ which models the fabric as an anisotropic membrane material. On the other hand, the polymer behavior is modeled with a constitutive model corresponding to Mooney-Rivlin material which is a particular case of material Law-42 of Altair Radioss ™ . More details about these constitutive models can be found in [31]. The material parameters used in the simulation are listed in Tab. 2.4 which have been obtained from the experimental and characterization campaigns [START_REF] Dile | Apprehension et simulation de la mise en forme de preformes nappees discretes par thermocompression[END_REF][START_REF] Marroquin | Variabilité dans la simulation de thermoformage[END_REF]. The value of viscosity has been obtained from the specifications provided by the prepreg supplier.

The ply-ply interaction is modeled with a penalty-based contact (Type 25) formulation in Altair Radioss ™ which has been modified to incorporate interply adhesion developed in this work. Also, specially defined rigid to deformable contacts are defined between the die/punch and plies. The interply adhesion parameters used for the simulation are listed in Tab. 2.5.

Woven (Law-58) UD (Law-58) Polymer (Law-42) The simulation was run on a Intel(R) Xeon(R) CPU E3-1545M v5 @ 2.90GHz (x86_64) machine with SPMD=1 and OMP=24 type of parallel environment. The simulation was run with the natural time-step of the model (approximately 3.0E-7 seconds) and in total 113,458 cycles were completed for the total time of 0.034 seconds. The CPU time used in the simulation was found to be 75,350 seconds whereas the actual time taken for the simulation to complete was found to be 3,198 seconds.

ρ 1.8E-9 t/mm 3 ρ 1.8E-9 t/mm 3 ρ 1.13E-9 t/mm 3 E 1 10,000 MPa E 1 25,000 MPa ν 0.495 E 2 10,000 MPa E 2 10 MPa µ 1 2.0 G 0 0.069 MPa G 0 0 MPa µ 2 -1.0 G T 3.0 MPa G T 3.0 MPa α 1 2.0 α lock 55 • α lock 85 • α 2 -2.0
For a qualitative validation, the ply positions at an intermediate time step were compared for the case without interply adhesion mechanism against the case with interply adhesion mechanism. A resistance free delamination of plies was observed when interply adhesion was not modeled (Fig. 2.16a). On the other hand, with the added interply adhesion model, the prepreg plies stick together (Fig.

Comparison of the final position of a prepreg patch

The vertical displacement of a corner (Point P shown in Fig. 2.18) was measured from the experiments and it was compared with values obtained from the simulation. The experimentally obtained value was 6 mm (Fig. 2.18a). The simulation without the interply adhesion mechanism predicted this value to be 31.7 mm (Fig. 2.18b) which is significantly different from the real value. On the other hand, adding interply adhesion mechanism improved this value to 6.6 mm (Fig. 2.18c) which is much closer to the experimentally observed value. Therefore, it can be concluded that with added interply adhesion contact mechanism, there is a significant improvement in the prediction of the final ply positions.

A sensitivity analysis was performed on the adhesive strength parameter. The final positions of 6 nodes situated at the corners of Ply 9 (Fig. 2.20) were recorded by post-processing the simulation results. The simulations were performed by varying σ adh from its reference value of 0.0062 MPa. The variation was performed from σ adh = 0.0012 MPa (reference value -80%) to σ adh = 0.012 MPa (reference value +80%). Tab. 2.6 shows the variation of the final position (magnitude) of these nodes in mm with respect to the reference case (σ adh = 0.0062 MPa). In general, σ adh can be thought of as an entity that opposes the separation forces experienced during the forming process. Thus, as long as σ adh exceeds the local separation stress experienced by the plies, the final position of plies remains more of less the same, indicated by the almost similar values of relative distance for positive variation of σ adh in Tab. 2.6. On the other hand, if σ adh is less than the local separation stress, there will be local delamination which can eventually create significant variation in the final position of the plies. This indeed was observed as σ adh value was reduced from its reference value, at the value of 0.0012 MPa, the maximum final position variation was found to be 21.63 mm.

Comparison of fibre orientations for individual plies

After forming this component using QSP ® , the fibre orientations at various locations were experimentally measured for each ply. In order to perform these measurements, the poly-σ adh (MPa) 0.0012 0.0037 0.005 0.0056 0.0068 0.0074 0.0087 0.0112 -80% -40% -20% -10% +10% +20% +40% +80% The fibre orientations measured from the experiments at these locations were compared with the values obtained from the simulation. Details of this information can be found in Tab. 2.7 and Tab. 2.8 for fibre orientation measurement from experiments and simulation respectively.

The nomenclature for the tables is as follows: A point (at which fibre orientation is measured) is identified as PiDj where i denotes the ply number and j denotes the family of fibres (1 or 2). The column Loc contains the identifier for the measurement location of plies (Refer to Fig. 2 The results summarized with respect to each face are presented in Tab. 2.9. It can be seen that for Faces 1, 4, 5 and 6, almost all values from simulation are in quite a good agreement with the experimental observations. For Face 2, even though the overall predictions are in good agreement, there are some locations where some difference can be observed. For Face 3, all values are in good agreement with an exception of just one location (Location 3.4 in Tab. 2.7 and Tab. 2.8) where the simulation differs significantly from the experiments. This in turn has increased the average of absolute difference for this face. Overall, the fibre orientations predicted by the simulation were on an average 5.9 • away from the experiments.

However, it is to be noted that in the current simulation only the forming phase of QSP ® is simulated. It is known that the mechanisms occurring during consolidation such as squeeze flow can play a key role on this face as shown in [23]. Thus, it is necessary to develop the capabilities for the simulation of the consolidation phase in order to investigate this further. This is the subject of the developments done in Chapters 3 and 4. This model will be revisited in Chapter 5 in which the developments of this chapter along with those in Chapters 3 and 4 are coupled to build a unified solution for both forming and consolidation

Conclusion

In this work, the importance of modeling interply adhesion for QSP ® has been demonstrated. It is worth to recall the peculiarities of this process that make the usual models not suitable: (a) usage of discontinuous prepreg patches (b) inability to use a blank holder (c) long distance sliding of patches during the forming phase. Thus, a novel semi-empirical contact mechanism involving an adhesive mechanism in the normal direction and a viscous sliding mechanism in the tangential direction has been developed. The normal adhesive mechanism creates a finite adhesive strength for ply-ply contact. At the same time, it allows sliding of plies along with automatic recreation of adhesion if the plies come in contact again. In terms of computational performance, the model has been implemented for shared memory processors (SMP) as well as for domain decomposition with MPI in Altair Radioss ™ showing its scalability.

In terms of characterization, the interply adhesion numerical model proposed here just requires the values of adhesive strength and viscosity of melt polymer. This makes it highly suitable for industrial applications where an extensive characterization campaign is often a challenge in terms of cost and time. A detailed experimental procedure has been developed to quantify the interply adhesive strength of molten thermoplastic polymer. Also, the sensitivity study showed that a very precise determination of the adhesive strength is not required for this model as long as the adhesive strength lies within a range that is above a threshold value.

Finally, it has been shown that, using this novel numerical model of interply adhesion coupled with adapted experimental data, the prediction of final positions of the plies has improved significantly. Also, the fibre orientation predictions are shown to be in good agreement with the experimental observations made on full-scale industrial QSP ® parts. Thus to conclude, this numerical model of interply adhesion has improved the numerical simulation of composite forming of viscous discontinuous thermoplastic prepregs.

Chapter 3

Development of a full-integration shell with transverse normal stress and selective mass scaling (QBATP)

Abstract

During the consolidation process of thermoplastic prepregs, there exists a through-thickness normal (transverse) stress also referred to as 'Pinching Stress'. This implies that the prepregs are in a three-dimensional state of stress. As the classical 5-parameter Mindlin shells are built with a plane stress assumption; they cannot be used for the numerical simulation of such a process. To overcome this, this chapter focuses on developing an extension to the existing hourglass-free full-integration shell element in Altair Radioss ™ (known as QBAT) in order to add a functionality of pinching stress. This is achieved by augmenting the QBAT element with additional nodal degrees of freedom that can simulate the pinching behavior. A selective mass scaling scheme has been added for this element which increases its critical time-step significantly. This element is referred to as 'QBATP' element in this work. With the added capability to capture pinching stress along with a higher time-step; the QBATP element is a suitable candidate to simulate the consolidation process. This element serves as a first step towards a unified full-scale simulation of both forming and consolidation phases in a generic composite forming process as well as a more specific one such as QSP ® .

Introduction

Modeling of consolidation is an important and necessary step in the full numerical process simulation of a thermoplastic composite forming process. This is because in case of QSP ® and also for a typical thermoplastic composite forming process, the forming stage is almost always followed by a consolidation stage. This through-thickness compression happening during the consolidation process is essential to ensure a proper consolidation of multi-layered thermoplastic composites. The consolidation phase serves several other important benefits such as removal of porosities, healing of ply-ply interfaces, obtaining net-shape parts and improvement of the surface quality to name a few.

Numerical modeling of a consolidation process

The through-thickness compressive force that is applied during the consolidation of composite prepregs is the driving force for various mechanisms such as squeeze flow [1], in-plane resin bleed [2] and transverse resin percolation [3,4]. Thus, during consolidation the material is in a three-dimensional state of stress which is an important aspect to model in the numerical simulation of a consolidation process.

In the literature, there have been three main directions for modeling approaches of a consolidation process from the point of view of the finite element used in the simulation: (1) Solid (Hexahedral) element based approach (2) Solid-shell element based approach (3) Shell element based approach. Each of these approaches has some unique advantages and limitations. An overview of each of these approaches is given in the next sections.

Solid based approach

The main idea behind the usage of a solid (Hexahedral) element to model consolidation is that this type of element can readily capture a three-dimensional state of stress without requiring any changes in the existing element formulations. Also, unlike the classical shells, a modification of the constitutive models is not required for solids. This means it is relatively easier to develop more complicated constitutive models in the existing framework of constitutive models. Such a framework for the constitutive models for solids almost always exists in the in-house or the commercial Finite Element codes.

Numerical modeling of a consolidation process conducted by Belnoue et al. in [5] is performed on a relatively simple (flat) geometry consisting of multiple plies. In this work, each individual ply was modeled using a Hexahedral element (C3D8 element of ABAQUS). This work was further continued on a relatively more complex part in [6] where the simulations were carried out in ABAQUS/Standard with the addition of specific user material models.

An Hexahedral element with 27 nodes with a mixed formulation for velocity, pressure (to model incompressibility) and fibre tension (to model the inextensibility in the fibre direction) was used in [1,7]. The focus of this work was to demonstrate the transverse squeeze flow effect and the importance of modeling a resin-rich layer. This was based on the work of [8] to model the squeeze flow behavior of unidirectional prepregs by treating it as an transversely isotropic viscous fluid.

There are however, several important limitations of a solid element. First, the computational cost. In case of explicit solvers, when a solid element is used to model a thin structure; the time-step is usually governed by the thickness. As the consolidation proceeds, it reduces the thickness further; the time-step reduces as well making this approach computationally very expensive. Even in the case of steady-state implicit approaches such as [1], a geometry of size 5 mm × 5 mm × 5 mm when meshed with 24 × 24 × 49 elements (resulting in about 750,000 DoFs to solve) required 3 hours to solve using 16 cores of two Opteron 6328 @ 3.2 GHz and 350 GB of memory. Secondly, the solid elements are too stiff/rigid in bending mode because of transverse shear locking problems thus not suitable for bending without making some changes in their formulation. Thus, if one wants to develop a unified solution for both forming and consolidation phases for models of industrial scale such as in QSP ® , this approach is not viable.

One approach of using shells for forming phase and then later replacing them with the solid elements for consolidation phase is theoretically possible. However this approach has several challenges which are discussed in detail in Chapter 5 making it unusable in real-life industrial level applications.

Solid-shell based approach

The main advantage of solid-shells is that they only have translational DoFs and no rotational DoFs as is the case with shells. This therefore results in a simplified kinematic description of the element. Compared to solids, a solid-shell element is able to overcome the limitation of the shear locking and thus can be used for bending.

A new prismatic solid-shell element with 7 nodes referred to as SB7γ19 was developed in [9][10][11]. The 7 th node is located at the center of the element with only one displacement DoF along the thickness direction and allows for a quadratic interpolation for throughthickness displacement field. This was then coupled with a viscoelastic constitutive model based on Maxwell approach for compaction in [12] in order to simulate the consolidation process. An indicator for measuring a degree of consolidation (degree of intimate contact) was also proposed. This work provided a unified approach for modeling the simulation of forming and consolidation based on a solid-shell element. The effects such as transverse squeeze flow were not modeled as the focus was more on the woven prepregs instead of UD.

A new hexahedral solid-shell element called SB9 was proposed in [13] in order to have the capability of pinching stress. It has an additional 9 th node at the volume center. Unlike the other 8 nodes, this node has only one DoF which is in the thickness direction. With this extra node, a quadratic interpolation of the through-thickness displacement is possible thereby avoiding the Poisson thickness locking problem as well as allowing the use of a full three-dimensional form of the constitutive model. This element has been implemented in a static and implicit solver of the software 'Code Aster' developed by the company EDF.

Another 8 node hexahedral solid-shell element derived from the work of [14] was extended and implemented in ABAQUS [15] for isotropic material and later extended for anisotropic material model in [16]. Even though the through-thickness behavior was not specifically investigated in these works; investigation was more focused on the avoiding various types of element locking for solid-shell elements followed by a comparison with various solid and shell elements from the ABAQUS element library.

A solid-shell element approach overcomes several important limitations of a solid element approach. If this could be coupled with the selective mass scaling, it has the potential to model both the forming and consolidation phases of a typical thermoplastic composite forming process. Even though the elements in this category have shown a great potential in modeling the numerical simulation on relatively simpler components; a demonstration on a full-scale industrial model with many plies and patches (in case of QSP ® ) is not available yet.

Shell based approach

Historically, the shell elements were specifically designed with a purpose to model thin structures using a three-dimensional surface that can be parametrized using two curvilinear coordinates. In general, shells are fast, easier to mesh with and are therefore used extensively in the design and analysis of various real-life structures.

There have been numerous contributions in the developments of shell elements (Refer to [17,18] for more details about the history of research in the field of shell elements). In this section, only a brief overview is given to introduce the main theories for modeling of shells.

The simplest form of classical shells (3 parameter shells) are based on what is known as the Kirchhoff-Love assumptions [19,20] (representation shown in Fig. 3.1a) which makes three assumptions :

1. the shell is in a plane state of stress i.e. σ zz = 0 2. the shell cannot have a transverse shear i.e. γ xz = γ yz = 0 3. no extensible director i.e. zz = 0. However, a usual practice is to calculate zz by solving the equation σ zz = 0 within the constitutive models. This implies that the thickness of the shell technically can increase/decrease but it has to come from the in-plane compression/tension respectively. One cannot directly modify thickness by applying a through-thickness loading.

Since, ( 1) and ( 3) are contradictory conditions; this type of shell requires modification of the constitutive model and a direct three-dimensional constitutive model cannot be used for this type of shell.

The condition related to the transverse shear was relaxed later with an updated model with 5 parameters commonly known as the Reissner-Mindlin theory [21,22] (representation shown in Fig. 3.1b). Thus, with this theory, one can have a nonzero transverse shear γ xz = 0, γ yz = 0. However, the other two assumptions: plane stress (σ zz = 0) and inextensibility of the director ( zz = 0) still remain, which still require the modification of the constitutive model. Some early works relating to the development of three-dimensional shells based on this theory can be found in [23,24]. QBAT element in Altair Radioss ™ is an hourglassfree full-integration 4-node shell element that falls within this category. It is based on the Q4γ24 element discussed in [25] which in turn is based on the MITC4 element discussed in [26,27].

A 6-parameter shell was later proposed in [28] that added a possibility of a stretching of the director vector (representation shown in Fig. 3.1c). This assumed a linear displacement field along thickness resulting in a constant zz . Thus, with this type of shell; it is possible to have σ zz = 0 and zz = 0 and the plane stress assumption is no longer necessary thereby adding a capability of a three-dimensional state of stress for the shell. However, it is not possible to directly use the three-dimensional constitutive model for this type of shell due to the Poisson thickness locking. The option to resolve this locking is either by decoupling the transverse normal stress with the bending strains [29] or to use a quadratic interpolation of displacement along the thickness [30,31]. The first use of this type of shell in composite process modeling community is attributed to [29] where this type of shell was used to qualitatively demonstrate the importance of capturing the transverse normal stress (σ zz ). Its effect on porosities and subsequently the quality of consolidation in various regions was discussed.

A 7-parameter shell was first proposed in [32,33] which uses a quadratic displacement field along thickness direction (and therefore a linear zz and σ zz ). Despite the complexity involved in its formulation; this by far is the most comprehensive shell formulation for which a full three-dimensional constitutive model can be used directly without any modifications.

A further discussions about this shell can be found in [34][35][36].

To summarize, a general comparison of various different shell theories is shown in Tab. 3.1 which highlights their advantages and limitations.

On the choice of the approach

In conclusion, based on the literature review discussed here; the challenge of numerical process modeling of consolidation of thermoplastic prepregs can be addressed with three main approaches. It is evident that Solid element based approach cannot be used due to its computational cost and challenges in its ability to simulate the forming phase of the process. The Solid-shell elements definitely look quite promising but a demonstration on a full scale industrial model with many plies and patches (in case of QSP ® ) is not available yet. With these observations, a Shell element based approach was chosen for this work.

Clearly the 3-parameter and 5-parameter shell cannot be used for consolidation as they lack a three dimensional representation of stress. The 7-parameter shell was not chosen due to its complexity as well as the additional computational cost (compared to a 6 parameter shell). Thus, based on this and based on the elements available in the current element library of Altair Radioss ™ , a modeling choice was made to extend the QBAT element of Altair Radioss ™ into a 6-parameter shell referred to as 'QBATP' element in this work.

Objectives and content of the study

The overall goal of this work is to develop a full-integration shell element with transverse normal stress (QBATP) and selective mass scaling. The objectives of this work are to:

• Extend the formulation of the existing full-integration shell element (QBAT) in Altair Radioss ™ in order to add the capability of transverse normal stress.

• Add a selective mass scaling mechanism for this element in order to increase its critical time-step giving an advantage in the computational performance.

• Carry out validation tests for this element.

• Implement the interaction of QBATP element with a contact interface and perform validation.

• Demonstrate the capabilities of QBATP element on a semi-industrial part designed by CETIM with the aim to propose a unified simulation of both forming and consolidation phases.

Full-integration shell element with transverse stress (QBATP)

This section discusses the development of the QBATP shell element. The name 'QBATP' stands for a 'Quadrilateral BAToz element with Pinching'. Starting with the discussion about the idea behind the pinching and transverse stress; the geometry and kinematics of this shell element are described next. This is then followed by the calculation of various B matrices in order to calculate strain-rates and strains, followed by a discussion about the constitutive model. Finally, the system of dynamic equations and the explicit time integration technique are discussed. The geometry of a general shell element can be represented as shown in Fig. 3.2a where it is represented by a three-dimensional mid-surface x p (ξ, η) representing the points located on the mid-surface. With this, the position of a general point x q (ξ, η, ζ) within the shell can be written as follows,

Concept of transverse normal stress and pinching

x q (ξ, η, ζ) = x p (ξ, η) + z(ζ) n(ξ, η) where z = 1 2 ζh (3.1)
where ξ, η and ζ are the parametric coordinates.

In general, the nodal fiber direction ζ may not correspond to the mid-surface normal n but the QBAT shell element is built with a postulate that they are equivalent as was the case with the Q4γ24 developed in [25]. This remains true for the case with a uniform thickness within an element or in case of a mild-taper conditions as discussed in [18]. Also, it is considered that the mid-surface represented by the vector x p is continuous however the normals could be discontinuous across elements as shown in Fig. 3.2b.

With this description, the displacements for the general point x q are given by,

u q = x t 0 +∆t q -x t 0 q (3.2) = x t 0 +∆t p + ζ 2 h t 0 +∆t n t 0 +∆t -x t 0 p - ζ 2 h t 0 n t 0 (3.3)
Using the first order approximation of Taylor series for h t 0 +∆t and n t 0 +∆t we get,

u q = x t 0 +∆t p -x t 0 p Term 1 + ζ 2 h t 0 (n t 0 +∆t -n t 0 ) Term 2 + ζ 2 (h t 0 +∆t -h t 0 ) n t 0 Term 3 (3.4)
The Term 1 in Eq. 3.4 corresponds to the displacement DoFs of the mid-surface of the shell.

The Term 2 corresponds to the rotation DoFs of the mid-surface. The Term 3 corresponds to the change in the thickness. This phenomenon is neglected in case of the classic shell theories; as the fiber (not to be confused with the fibres that provide reinforcements for the composites) is assumed to be inextensible [37]. This results in the consideration of only first two terms for the classical shells.

However, the Term 3 is not neglected for the 6 or 7 parameter shells. This term is the basis of creating an additional nodal degree of freedom locally referred to as the 'Pinching degree of freedom'. The word 'pinching' refers to the action of applying a force/pressure through-thickness in this work. Taking into account the Term 3 in the shell formulation, it is possible to have a non-zero transverse normal stress σ zz .

Shell element description

As discussed before, the QBATP element developed and discussed in this work is basically an extension of the existing element QBAT in the Altair Radioss ™ element library. The geometry of this shell is shown in Fig. This being a full-integration shell element; 2 × 2 Gauss points on the mid-surface are used for integrations as indicated by the '×' symbol in Fig. 3.3b. The coordinates and integration weights for the Gauss points are given in Tab. 3.3. Within the element, each node (i) has 6 local degrees of freedom; 3 for standard midsurface displacements (U 1i , U 2i , U 3i ), 2 for rotations (θ 1i , θ 2i ) and additionally 1 for the pinching displacement (U p i ) as follows:

Point ξ η Weight 1 -1/ √ 3 -1/ √ 3 1 2 +1/ √ 3 -1/ √ 3 1 3 +1/ √ 3 +1/ √ 3 1 4 -1/ √ 3 +1/ √ 3
[U i ] (6×1) = U 1i U 2i U 3i θ 1i θ 2i U p i T (3.6)
Thus, the nodal velocity vector for the full element is,

[V ] (24×1) = ... V 1i V 2i V 3i ω 1i ω 2i V p i ... T (3.7)
which are defined with respect to the orthonormal coordinate system (t 1i , t 2i and n i ) defined at each node as shown in Fig. 3.3c. Additionally, one must define the covariant basis vectors (a 1 and a 2 ) as shown in Fig. 3.3d, which are not necessarily orthogonal to each other. More discussion about their derivation can be found in Appendix A of this chapter.

The additional nodal degrees of freedom for pinching are defined with respect to the midsurface. The pinching displacement U p corresponds to the displacement of the artificial nodes located at the top and bottom surface of the shell. This being a relative degree of freedom, a sign convention needs to be defined. The increase of thickness is considered positive i.e. U p > 0 whereas a reduction of the thickness is considered to be negative i.e. U p < 0. This is shown in Fig. 3.4. Besides the postulates used for QBAT element discussed before, some additional computational liberties were taken in this work as follows:

• The integrations were performed with respect to an average thickness (h) within an element. This is justified by the assumption of the so-called 'mild-taper' condition described in [18] which signifies that the spatial changes of thickness within the element are not significantly big. Even though this is an approximation; the areas with a high gradient of thickness variation could be simulated with a basic local mesh refinement.

• The influence of the pinching on the transverse shear which was considered in [29] is not modeled in this work. This is mainly because the main idea here is to use a shell component for each individual ply of the prepreg stack and considering the aspect ratio of a single prepreg ply; it is fair to assume that there will not be a significant transverse shear within a thin ply. It is to be noted that only the pinching contribution to the transverse shear is neglected, the classical contributions coming from the bending are retained.

Calculation of strain rates and strain

With the displacements and velocities defined, the next step is to define the 'strain-displacement' or 'strain rate-velocity' matrices commonly known as the B matrices. The strain rates are calculated based on the rate of deformation tensor,

˙ ij = D ij = 1 2
∂v i ∂x j + ∂v j ∂x i = 1 2 L + L T (3.8)
The calculation of strain rates is split into four strain interpolation matrices (1) Membrane

[B m ] matrix for Membrane

The membrane strain-rate is calculated as follows:

  ˙ xx ˙ xx ˙ xy   = [B m ] (3×24) [V ] (24×1) (3.9)
Where the [B m ] matrix is given by:

[B m ] (3×24) = B m1 (3×6) B m2 (3×6) B m3 (3×6) B m4 (3×6) (3.10) B mi (3×6) =    t T 1 N i,x 0 0 0 t T 2 N i,y 0 0 0 t T 1 N i,y + t T 2 N i,x 0 0 0    (3.11) [B b ] matrix for Bending
The curvatures are calculated as follows:

  χxx χyy χxy   = [B b ] (3×24) [V ] (24×1) (3.12)
Where the [B b ] matrix is given by:

[B b ] (3×24) = B b1 (3×6) B b2 (3×6) B b3 (3×6) B b4 (3×6) (3.13) B bi (3×6) =    t T 1 N bc 1i s 1i N i,x 0 t T 2 N bc 2i s 2i N i,y 0 t T 1 N bc 2i + t T 2 N bc 1i s 1i N i,y + s 2i N i,x 0    (3.14)

[B p ] matrix for Pinching

The pinching strain-rate is calculated as follows:

˙ zz = [B p ] (1×24) [V ] (24×1) (3.15)
Where the [B p ] matrix is given by:

[B p ] (1×24) = B p1 (1×6) B p2 (1×6) B p3 (1×6) B p4 (1×6) (3.16) B pi (1×6) = 0 0 0 0 0 2N i h (3.17)

[B ts ] matrix for Transverse shear

While considering the behavior of transverse shear, it is important to consider the phenomenon of 'Transverse Shear Locking'. This occurs mainly because of the inability of the C 0 shell finite elements to reproduce a pure bending behavior. The coupling between the normal and shear strains for linear elements results in a parasitic transverse shear strains in case of pure bending. This therefore, results in numerical stiffening in bending creating a poor convergence.

There are several solutions discussed in the literature to avoid this type of locking. The methods such as reduced integration of the transverse shear as discussed in [18], usage of C 1 elements as discussed in [37], 'Enhanced Assumed Strain' (EAS) method proposed in [38] and the ' Assumed Natural Strain' (ANS) method first proposed in [39] can be used. This idea of ' Assumed Natural Strain' (ANS) method has been discussed in [26] and further used for Q4γ24 element in [25]. In this work, the same approach has been used for the QBATP elements. The method interpolates the transverse shear strain from the values of the covariant components of the transverse shear strains at 4 edge mid-points as follows (Fig. 3.5), With this, the transverse shear strain-rate is calculated as follows:

γ ξ = 1 -η 2 γ A1 ξ + 1 + η 2 γ A2 ξ (3.18) γ η = 1 -ξ 2 γ B1 η + 1 + ξ 2 γ B2 η (3.19)
˙ xz ˙ yz = [C 0 ] T (2×2) [B ts ] (2×24) [V ] (24×1) (3.20)
Where the [C 0 ] is the local basis at z = 0 (Refer to Appendix A) and the [B ts ] matrix is given by:

[B ts ] (2×24) = B 1 ts (2×6) B 2 ts (2×6) B 3 ts (2×6) B 4 ts (2×6) (3.21) [B 1 ts ] = -nA1 T +AS1 T A1 0 -nB1 T +AS1 T B1 0 [B 2 ts ] = +nA1 T +AS2 T A1 0 -nB2 T +AS1 T B2 0 (3.22) [B 3 ts ] = +nA2 T +AS2 T A2 0 +nB2 T +AS2 T B2 0 [B 4 ts ] = -nA2 T +AS1 T A2 0 +nB1 T +AS2 T B1 0 (3.23)

Constitutive model

The incremental form of the elastic constitutive model discussed in [29] was used in this work as a first step. The membrane and bending contributions to the stress which use the corotational framework are as follows:

        ∆σ xx ∆σ yy ∆σ zz ∆σ xy ∆σ yz ∆σ xz         =            E(1-ν) (1+ν)(1-2ν) Eν (1+ν)(1-2ν) Eν (1+ν)(1-2ν) 0 0 0 Eν (1+ν)(1-2ν) E(1-ν) (1+ν)(1-2ν) Eν (1+ν)(1-2ν) 0 0 0 Eν (1+ν)(1-2ν) Eν (1+ν)(1-2ν) E(1-ν) (1+ν)(1-2ν) 0 0 0 0 0 0 E 2(1+ν) 0 0 0 0 0 0 E 2(1+ν) 0 0 0 0 0 0 E 2(1+ν)                    ∆ xx ∆ yy ∆ zz ∆ xy ∆ yz ∆ xz         (3.24) +          E 1-ν 2 Eν 1-ν 2 0 0 0 0 Eν 1-ν 2 E 1-ν 2 0 0 0 0 0 0 0 0 0 0 0 0 E 2(1+ν) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                  z∆χ xx z∆χ yy z∆χ xy 0 0 0         (3.25)
Compared to the classic QBAT shell, the modifications are as follows:

• The membrane stiffness matrix for QBATP element does not use the plane stress assumption. Therefore, it resembles the full three-dimensional stiffness matrix for elasticity.

• QBATP element has an additional stress compared to QBAT shell i.e. the transverse normal stress σ zz .

• σ zz has been decoupled from the bending DoFs in order to avoid pinching locking (or Poisson thickness locking) as discussed in more detail in [29].

In general, the strains and stresses are numerically evaluated at several through-thickness integration points at each Gauss point on the mid-surface and then the stress resultants and stress couples are evaluated by performing a numerical summation. However, in this work, an alternate approach has been followed which involves performing the through-thickness integrations explicitly and then calculating the stress resultants and couples directly. This approach is known as 'global approach' as it does not use any points of integration in thickness direction. This approach can not be considered to be very generic as it requires calculation of stress resultants and couples explicitly; but it is adopted here as the constitutive model used in this work requires a different treatment for membrane versus bending in order to avoid locking.

However, the benefit of this approach is its speed; as the calculations are performed only on the Gauss quadrature points (4 locations) on the mid-surface. On the other hand, in case of a general approach the calculations will have to be performed at each integration point along thickness for every Gauss point on the mid-surface. If we use 5 points of integration along thickness, this would result in a total of 20 locations.

Stress internal force relationship

The internal force vector (F int ) for the element is then calculated as,

F elem int = elem B T σdV (3.26)
This three-dimensional integration is then split into a through-thickness integration and subsequently a surface integration over the element's mid-surface. Using the global approach, the through-thickness integrations result in the explicit calculation of the stress resultants and stress couples. The stress resultants are calculated as follows:

N x = h/2 -h/2 σ xx dz, N y = h/2 -h/2 σ yy dz, N z = h/2 -h/2
σ zz dz (3.27)

N xy = h/2 -h/2 σ xy dz, N yz = h/2 -h/2 σ yz dz, N xz = h/2 -h/2 σ xz dz (3.28)
And the stress couples are calculated as follows:

M x = h/2 -h/2 σ xx zdz, M y = h/2 -h/2 σ yy zdz, M xy = h/2 -h/2 σ xy zdz (3.29)
This results into a F int for the element as follows,

[F int ] (24×1) = ... F 1i F 2i F 3i M 1i M 2i F p i ... T (3.30)
These are then assembled after performing necessary coordinate transformations into separate global vectors for Forces (F int ), Moments (M int ) and Pinching force (F p int ).

External force and pinching pressure

In order to impose an external load for pinching the element, a surface type of loading for pinching was developed. It is referred to as 'Pinching pressure' in this work. This involves application of an equal and opposite surface load on the top and bottom side of the shell (as shown in Fig. 3.6). Based on the sign convention for pinching (Refer to Fig. 3.4), the external pinching force on each node (i) of the element is calculated as,

F p ext (i) = 1 4 (2P ) A n (3.31)
Where A is the area of the element and n is the normal to the mid-surface at the center of the element.

The above formulation makes certain assumptions as follows: • The normal n is calculated at the center of the element and on the mid-surface of the shell. Even though this is not perfect as one should ideally consider the normals on the top and bottom surface of the shell instead of the normals to the mid-surface. This is not a strong assumption in case of thin shells where the gradient of thickness within the element is assumed to be not very significant.

• Dividing the pinching force equally among the 4 nodes of the element is again not perfect especially for a distorted element but it is done here for simplicity.

Dynamic system of equations

The assembly of internal and external forces, moments and pinching forces results in the following global dynamic system of equations. It consists of three sets of equations each corresponding to displacements, rotations and pinching displacements respectively.

M Ü = F ext -F int (3.32) I θp = M ext -M int (3.33) M p Ü p = F p ext -F p int (3.34)
Where M and M p are lumped mass matrices corresponding to the displacements and pinching displacements respectively. I is the diagonal moment of inertia.

An explicit central difference scheme that falls under the broad category of Newark family scheme with α = 1/2 and β = 0 is used in Altair Radioss ™ [40]. With this scheme, the above dynamic system of equations is solved by first calculating the accelerations at t n as:

Ü (t n ) = [M ] -1 (F ext (t n ) -F int (t n )) (3.35) θ(t n ) = [I(t n )] -1 (M ext (t n ) -M int (t n )) (3.36) Ü p (t n ) = [M p (t n )] -1 (F p ext (t n ) -F p int (t n )) (3.37)
Once the accelerations are known; the velocities are calculated as follows:

U (t n+ 1 2 ) = U (t n-1 2 ) + dt 1 + dt 2 2 Ü (t n ) (3.38) θ(t n+ 1 2 ) = θ(t n-1 2 ) + dt 1 + dt 2 2 θ(t n ) (3.39) U p (t n+ 1 2 ) = U p (t n-1 2 ) + dt 1 + dt 2 2 Ü p (t n ) (3.40)
And once the velocities are calculated; the displacements can be calculated as follows:

U (t n+1 ) = U (t n ) + dt 2 U (t n+ 1 2 ) (3.41) θ(t n+1 ) = θ(t n ) + dt 2 θ(t n+ 1 2 ) (3.42) U p (t n+1 ) = U p (t n ) + dt 2 U p (t n+ 1 2 ) (3.43)
Note that this solution scheme is written for a general case with the use of two time-steps dt 1 = t nt n-1 and dt 2 = t n+1t n . This is because in general the time-step is not constant and changes from one cycle to the next.

This completes the discussion about the development of the QBATP shell and now the focus of the next section is to increase the critical time-step of this element in order to improve its computational performance.

Selective mass scaling for QBATP element

In general, the idea behind any of the mass scaling techniques is to use a higher time-step for the simulation as this would mean a lesser computational cost. Mass scaling is in general very effective in the cases where the high frequency effects are not significant and can be neglected such as forming [37]. Starting with a discussion about the stability and timestep, this section then discusses various mass scaling techniques and finally discussion how a specific selective mass scaling technique was developed and used within this work for QBATP element.

Stability and time-step

The benefit of using an explicit scheme is that the explicit time integration is quite robust; if compared to its counterpart -implicit scheme. An implicit scheme can fail to converge and subsequently stop due to the failure of the numerical algorithm which can occur quite often in case of a nonlinear simulation involving contacts.

But the limitation of the explicit time integration is that it is only conditionally stable. For stability, the model should satisfy the Courant condition (also known as the CFL condition). This condition imposes a limit on the maximum time-step that can be used in case of an explicit time integration with a central difference integration scheme for it to be stable. It is given as follows,

∆t = α∆t crit ∆t crit = 2 ω max ≤ min 1≤I≤N el 2 ω I = min 1≤I≤N el L c I c I (3.44)
Where t crit is the critical time-step. The parameter α is a reduction factor, a typical choice of 0.8 ≤ α ≤ 0.98 is used in the numerical simulation [37]. In Altair Radioss ™ , the parameter α can be chosen by the user but conventionally α = 0.9 is used. In the above expression, ω max is the maximum eigenfrequency of the assembled mesh whereas ω I , L c I and c I correspond to the eigenfrequency, characteristic length and the speed of sound in an element I in the mesh. The characteristic length of an element is defined as the smallest geometric dimension of the element. For solids or solid-shells (Fig. 3.7a), it is given as the ratio of the element volume to the highest surface side surface area, which basically is equivalent to the smallest distance between two opposite faces [40]. Thus, this would be of the order of the thickness (h) of the element.

L c solid = Vol Largest Area (3.45)
On the other hand, for the classic shells, as there is no momentum balance solved across the thickness, the characteristic length for the classic shells (Fig. 3.7b) is calculated as,

L 1 = Area max(AC, BD) (3.46) L 2 = min(AB, BC, CD, AD, AC, BD) (3.47) L c shell = min(L1, L2) (3.48)
Thus, for a classic shell, the time-step depends solely on the mid-surface dimensions and it does not depend on the thickness. In case of pinching shell (QBATP), with the inclusion of additional pinching DoFs, the time-step of QBATP element resembles the BRICK8 (solid) element rather than the QBAT element (classic shell element). This is because now with QBATP element, additional momentum balance equations are written across thickness thereby imposing a limit in the CFL condition just like in the case of solids or solid-shells.

When the stack of prepreg is modeled on a ply-by-ply basis, naturally the thickness of an element is much less compared to the in-plane dimensions. This leads to a high transverse normal stiffness compared to the in-plane stiffness and the overall time-step of the model is governed by the highest stiffness (or smallest L c ) which is usually the thickness. This implies that by default (if mass scaling is not used) the time-step of the QBATP element would be similar to the BRICK8 element. The goal of the selective mass scaling technique discussed in this work is to increase the time-step of the QBATP element so that it is closer to the QBAT element time-step and much higher than the BRICK8 element. To the best of our knowledge, four different mass scaling techniques exist in the literature which are discussed next.

Method A: Uniform mass scaling

This is the simplest mass scaling technique that consists of increasing the density of the smallest k elements in order to meet the target time-step which is higher than the critical time-step of these elements. This decreases the speed of sound c in these element as c ∝ 1/ √ ρ. This therefore increases the critical time-step of these k elements. However, as demonstrated in [41] this type of mass scaling affects all deformation modes (and therefore all eigenfrequencies) of these k elements and can result in incorrect deformations. This would also result in a lack of conservation of rigid body translational momentum.

Thus, a better mass scaling technique is the one that can selectively affect only the higher eigenfrequencies and keep the lower eigenfrequencies unchanged (or have a very less effect on them). This category of mass scaling techniques is called "Selective Mass Scaling".

Method B: Selective mass scaling based on stiffness matrix

One such method discussed in [42] corresponds to adding an extra mass to the original mass (M ) using the stiffness matrix (K) to get a new scaled mass M scaled as follows,

M scaled = M + γK (3.49)
Where γ ≥ 0 is a scaling parameter. The benefit of this strategy is that, since K Ürigid = 0, the translational rigid body behavior ( Ürigid ) does not change due to this type of mass scaling. The eigenfrequencies of this new system with selective mass scaling become,

ω scaled i = ω 2 i 1 + γ ω 2 i (3.50)
Where ω i corresponds to the eigenfrequency of i th DoF of the original system. This shows that with this type of mass scaling; the higher eigenfrequencies will be impacted more (depending on the magnitude of γ) than the lower eigenfrequencies. However, there are two limitations of this approach. First, in a normal scenario without mass scaling; it is a standard practice to use a lumped mass matrix (M ) as this avoids the matrix inversion in the system of equations. However, with the added mass γK; the mass matrix is no longer diagonal. Thus, a matrix inversion becomes necessary. Second, as the added mass depends on K, which changes during the simulation in case of a nonlinear analysis, one would have to invert the matrix at every time-step which will be computationally very expensive and not practical.

Method C: Selective mass scaling based on local mass matrix

To overcome the limitations of the Method B, a second method of selective mass scaling was proposed in [42] with the idea that the added mass should not depend on K. Thus, for the case of a BRICK8 element; where the standard (lumped) mass matrix was given as,

M (24×24) =   m 0 0 0 m 0 0 0 m   where m (8×8) = m e 8      1 0 • • • 0 0 1 . . . . . . 0 1      (3.51)
In this case, the idea that was proposed was to use an added mass matrix as follows:

M a (24×24) =   m a 0 0 0 m a 0 0 0 m a   where m a (8×8) = β m e 56      7 -1 • • • -1 -1 7 . . . . . . -1 7      (3.52)
Where β ≥ 0 is the scaling factor and m e is the mass of a given element. The benefit of this method is that it still ensures M a Ürigid = 0 without depending on K. Thus, one can invert the new mass matrix M scaled = M + M a only once.

It was shown in [42] that even though the method with stiffness matrix based mass scaling (Method B) gave better results; this technique of redistribution of added mass is more practical and computationally less expensive. Of course, in case of deletion of elements and mesh refinement; one has to invert the mass matrix M scaled again adding to the computational as well as storage cost. One solution to this limitation was discussed in [43] where it was discussed that an iterative scheme based on conjugate gradient method could be used at each time-step to invert mass matrix M scaled .

Method D: Selective mass scaling based on acceleration filtering

This technique was proposed in [44] which was specifically targeting the selective mass scaling for a thin walled structure made up of solid (Hexahedral) elements. The key idea was to group the nodes of the mesh that are less than a chosen distance (d) from each other. For a thin walled structure this would mean grouping the nodes along the thickness direction. With this one can then define the individual node's acceleration (a i ) as the sum of the group acceleration (a g ) and the relative acceleration of the node with respect to the group acceleration (a r ),

a i = a g + a r where a g = f g m g and a r = f i m i -a g (3.53)
With these, one can recast the equation for a i in order to have a selective mass scaling as follows,

a i = g a g + r a r (3.54)
With g = 1 and r ≤ 1 which would signify that one reduces the speed of sound selectively for the relative degrees of freedom. The biggest advantage of this method is that, there is no need to invert the mass matrix (it remains diagonal) however the method requires the filtering of nodes based on the distance between them and modifying the overall system of equations in the case of solid (Hexahedral) elements. Another limitation of this work is about the choice of the scaling parameter r (and possible changes in its value during the simulation) which is not obvious to calculate in advance.

Selective mass scaling for QBATP element

In case of solid-shell elements, a technique similar to the acceleration filtering can be directly used to perform selective mass scaling along the thickness direction of the element. This was discussed in [45] where this strategy was developed and also a method for estimating the critical time-step for a parallelepiped element was discussed. This is helpful in estimating the scaling parameters for selective mass scaling techniques. This technique was discussed for a single layer solid-shell elements.

In the case of the pinching shell element (QBATP); since the added pinching degrees of freedom are relative to the mid-surface; a similar approach can be used. Thus it is possible to selectively increase the mass corresponding only to the pinching DoFs (M p ) without changing the mass for the mid-surface DoFs (M ) as follows,

M 0 0 M p Ü Üp = F ext -F int F ext p -F int p (3.55)
This mass scaling technique shall impact only the higher eigenfrequencies (dynamic behavior along thickness) of the model without impacting the low eigenfrequencies (dynamic behavior of mid-surface). To summarize, the advantages of this mass scaling approach are:

• This method falls under the category of the selective mass scaling, thereby avoiding the limitation of Method A.

• The time-step depends on the in-plane dimensions of the QBATP element and does not depend on its thickness.

• There is neither a construction of the mass matrix nor its inversion. This overcomes the limitations of the Method B and Method C.

• The method proposed here calculates (and updates) the mass scaling factor for pinching DoFs automatically during the simulation requiring no input from the user. This overcomes the limitation of Method D.

Mass scaling algorithm for QBATP element

Altair Radioss ™ has two possible ways to calculate the time step (1) Element based approach (2) Node based approach. The element based approach consists of calculating the critical time-step for each element as shown in Eq. 3.44, whereas nodal time-step is calculated based on the calculation of the maximum eigenfrequency by solving the global system,

det K * -ω 2 M = 0 (3.56)
Of course, it is computationally and memory wise expensive to (a) calculate the full stiffness matrix K * and (b) solve Eq. 3.56.

Thus, the strategy employed in Altair Radioss ™ for estimating the maximum eigenfrequency (and subsequently calculating the time-step) is to compute the nodal stiffness for each element and assemble it into a global diagonal matrix K. This, along with the nodal mass, is then used to calculate the maximum eigenfrequency and thereby the time-step.

A similar strategy was used for the pinching shell along with the equivalence of mass scaling to the geometric thickness scaling discussed in [46] for a solid-shell element. Thus, the timestep calculations remain the same as in case of the classical shell (QBAT) and the pinching mass is scaled accordingly in order to have the same maximum eigenfrequency as QBAT shell at every time-step.

Algorithm 1 Mass scaling for QBATP element

(i) = F ext p (i) -F int p (i) M p (i)

Implementation aspects

This section discussed the overview of the implemented framework in the industrial code of Altair Radioss ™ for the implementation of the QBATP element with selective mass scaling. This is further divided into three subsections which discuss a global overview, implementation framework for the internal forces and then the implementation of the constitutive models.

Global framework

Fig. 3.8 shows a global flowchart of the Altair Radioss ™ explicit solver. Altair Radioss ™ being a commercial solver, the actual architecture is significantly more complicated. However, this flowchart is presented here just to provide a simplified global overview of the implementations done in this Chapter.

The block Radioss Starter reads the input file containing all the details about the model: nodal coordinates, element connectivity, different elements, contacts, loads, boundary conditions and so on. This block was modified in order to read the additional information related to the pinching shell (QBATP) element and related functionalities. The block Radioss Engine reads the input file for the solver and various parameters such as time-step, start and end time, specific output requests and writing frequency etc. After the Allocation and initialization block, the explicit time-stepping algorithm begins.

It starts with the calculation of the external forces F ext which was appended to incorporate external pinching forces F p ext coming from the pinching pressure load (discussed in Sec. 3.3.6). After this step, the contact forces F contact are calculated. The interaction of contact and pinching in Type-25 contact interface (will be discussed in Sec. 3.7) was added here resulting in a calculation of the pinching contact forces F p contact . This is then followed by the calculation of the internal forces F int . The calculation of internal pinching forces F p int and their interaction with the classical internal forces was implemented here. After this step, the next time step ∆t is evaluated which is where the subroutines related to the selective mass scaling for pinching (discussed in Sec. 3.4) were added. This is then followed by the calculation of the classic accelerations Ü and pinching accelerations Ü p . Since, no kinematic constraints were implemented for the QBATP shell, the block related to the Kinematic constraints block remains as it is. After this, the classic velocities U and displacements U along with their pinching analogs U p and U p are calculated. Depending on the chosen output frequency; outputs are written in the output files. This completes one time-step of the solver after which the solver continues with the subsequent cycles until a final time t final is reached. The calculation begins with a loop over all QBATP elements and various elemental level calculations are performed which involve, but are not limited to, the calculation of the various coordinate systems within the element, calculation of the characteristic lengths and other preliminary calculations related to the constitutive models. This is then followed by a loop over Gauss points of the mid-surface and both B matrix of the classic QBAT element along with the B p matrix for QBATP element are evaluated. With this, the various strains are calculated next. The subroutines specific to QBATP element have been added to calculate the through-thickness normal strain zz .

Internal forces framework for QBATP element

This information is then passed to the constitutive model framework where the stress resultants and stress couples are evaluated as discussed in Sec. 3.3.4. This is then followed by the calculation of the internal forces both classic F int and pinching F p int which are subsequently assembled into a global vector after necessary local to global coordinate transformations. Also, elemental time-steps are calculated and for the purpose of the nodal time-step calculations, the nodal stiffness are calculated and assembled into a vector.

Constitutive models and their framework

The architectural framework for the most commonly used constitutive models in Altair Radioss ™ is divided into two groups, one for solid elements and other for shell elements. However, the current development of the QBATP element cannot fit into either of these frameworks because of the following challenges:

• For QBATP elements, the full three dimensional constitutive model cannot be used directly without any modifications due to locking issues discussed in Sec. 3.3.4. Therefore, direct use of the constitutive model framework of the solid elements is not possible.

• The classic shell elements are based on the plane-stress assumption, thus both their architecture and the constitutive model itself is different compared to the QBATP shells. Therefore, even the shell element framework of constitutive models cannot be used for QBATP elements.

Thus, a new framework was built for this category of shells as shown in Fig. 3.10. The developments discussed in this chapter are limited to the elastic material behavior and therefore, an elastic type of constitutive model (Law-1) has been added in this framework.

In the next chapter, an elasto-plastic constitutive model (Law-91) is discussed which has been added in this same pinching shell framework. This implemented framework serves as a platform for a future development of various constitutive models for this type of shell element. This completes the discussion regarding the development and the implementation of the QBATP shell element with selective mass scaling inside Altair Radioss ™ . The next section discusses various validation tests performed on this new element.

Validation Tests

In this section, several numerical tests are performed in order to validate the QBATP element. Selective mass scaling is enabled by default for all the test cases discussed here onwards. First, several single element tests are performed in order to assess the transverse normal behavior (Test 3.1), coupling of pinching and membrane loading (Test 3.2) and a case involving both pinching and in-plane shear loading (Test 3.3). For these cases, the results are compared with BRICK8 element of Altair Radioss ™ .

After that, several benchmark cases from the literature are compared with the results obtained with the QBATP element (Test 3.4 -Test 3.8). This set of tests also includes a test case comparing the performance of the QBATP element with the BRICK8 element (Test 3.5) demonstrating the advantage of selective mass scaling.

Test 3.1: Pure pinching pressure driven loading Objective

The objective of this single element test is to verify the response of QBATP element with the BRICK8 element of Altair Radioss ™ when it experiences normal transverse (pinching) type of loading. So in this case, σ zz < 0 whereas σ xx ≈ 0 and σ yy ≈ 0.

Setup

A single element of dimensions 10 mm × 10 mm with a thickness of 2 mm is considered. An increasing pinching pressure (through-thickness compression) is applied on this element that starts at 0 MPa and then reaching a maximum magnitude of 40,000 MPa (using a smoothed curve). The Young's modulus (E) was set to 100,000 MPa and a Poisson's ratio (ν) of 0.3 was used. The side view of the element is shown in Fig. 3.11.

Discussion

The test being symmetric in terms of X and Y directions, the displacements fields are identical in both X and Y. The evolution of the displacement field in Y direction is shown in Fig. 3.12 where it showed a good match between the results obtained with QBATP element with that of BRICK8 element. For a local quantitative comparison, the resulting nodal dis- placements of (N3 of QBATP element) in X, Y and Z (pinching) direction were compared (Fig. 3.13) to the corresponding node of BRICK8 element (same X and Y coordinate but on the top surface of the BRICK8 element) and it showed a good correlation between the two elements. It is obvious to expect a good match for the evolution of thickness (Fig. 3.14) once the pinching (Z) displacements (Fig. 3.13c) are known to be the same. The thickness starting at the initial value of 2 mm, reduces to a final value of 1.34 mm which is 33% reduction of thickness. In general for QSP ® , the % reduction of thickness during consolidation phase is about 15-25 %. The through-thickness pinching stress (σ zz ) of QBATP and BRICK8 is compared in Fig. 3.15a which shows a good match. Also, this being a nearly quasi-steady case, the curves match with the applied pressure as well.

The advantage of selective mass scaling becomes evident when one compares the evolution of time-steps for both elements (Fig. 3.15b). As discussed before, as the thickness decreases, the time-step for the hexahedral element keeps falling as it is governed by the thickness which can be seen for the BRICK8 element in this plot.

However, for QBATP element with selective mass scaling, the time-step does not depend on the thickness and is governed by the in-plane dimensions. Thus, we can see a much higher time-step in case of the QBATP shell when compared to the BRICK8 element. An interesting point to note is that as one reduces thickness; the in-plane dimensions increase due to Poisson's ratio which means the time-step actually increases with consolidation (Tab. 3.4). In total, the simulation with hexahedral element took 56,838 cycles to complete whereas the simulation with pinching shell was completed in just 6,913 cycles.

Test 3.2: Coupled pinching pressure and membrane force loading

Objective

The objective of this test is to test a multi-axial type of loading where both pinching pressure and membrane tension act simultaneously. In this case, σ xx = 0, and σ zz = 0.

Setup

The geometry consists of a single element of dimensions 10 mm × 10 mm with a thickness of 2 mm. The material parameters were set to the same values as for the Test 3.1 i.e. E = 100,000 MPa and ν = 0.3. A linearly increasing pinching pressure that starts at 0 MPa and reaches a maximum of 40,000 MPa was applied. Along with this, a linearly increasing nodal force stretching the element in X direction is applied. It reaches a peak value of 50,000 N. The side view of the element is shown in Fig. 3.16.

Discussion

The displacement fields (Fig. 3.17), nodal displacements (Fig. 3.18a) and normal stresses (Fig. 3.18b) showed a good correlation between the BRICK8 element and QBATP element for this test case of simultaneous transverse and membrane loading.

Test 3.3: Coupled pinching pressure and in-plane shear loading

Objective A typical loading scenario that might appear at certain locations during forming is when an element in a state of stress where it experiences both pinching forces as well as in-plane shearing. Thus, the goal of this test is to compare the behavior of the QBATP element with BRICK8 element when it undergoes in-plane shear and pinching simultaneously.

Setup

The geometry consists of a single element of dimensions 10 mm × 10 mm with a thickness of 2 mm. The material parameters were set to the same values as for the Test 3.1 i.e. E = 100,000 MPa and ν = 0.3. A linearly increasing pinching pressure that starts at 0 MPa and reaches a maximum of 40,000 MPa was applied. In order to create in-plane shearing, a displacement of 10 mm was applied along X and Y direction at node N3 whereas node N1 was fixed in both X and Y directions. This is shown in the setup (Fig. 3.

19).

Discussion

For the purpose of visualization, the initial (shown in wireframe mode) and deformed shapes (shown in grey color) for both QBATP and BRICK8 elements are shown in Fig. 3.20. For a quantitative comparison, the nodal displacements are node N2 for QBATP (and corresponding node in BRICK8 element) were compared (Fig. 3.21) and a good match was observed between the two elements. Similarly it was found that both the in-plane shear stress and pinching stress predicted by the QBATP element are in alignment with those obtained from BRICK8 element (Fig. 3.22a). The same observation can be made for the evolution of the thickness (Fig. 3.22b).

Test 3.4: Crushing under pressure

Objective

The goal of this test is to compare BRICK8 elements versus QBATP elements (with selective mass scaling) under uniform pinching pressure and compare the results with the values reported by Soulat et al. [29].

Setup

The geometrical dimensions are 50 mm × 50 mm × 4 mm. The material properties are: E = 120 GPa and ν = 0.3. The setup is shown in Fig. 3.23. This is a quarter model; therefore symmetry conditions are applied along left (Edge AD) and bottom (Edge AB) edges (QBATP) /faces (BRICK8). The geometry is meshed with an element size of 10 mm and for BRICK8 element only one element was used along the thickness direction. A linearly increasing pinching pressure of peak magnitude 10,000 MPa was applied on all elements and then once the maximum value was reached, the pinching pressure was linearly reduced to 0 MPa.

Discussion

The position of node located at point C with respect to point A in X direction (Refer to Fig. 3.23) is tracked during the simulation. This quantity was compared with the solution reported by [29] and also with the simulation using BRICK8 elements. This comparison is shown in Fig. 3.24a. The maximum value of this quantity occurs in the middle of the simulation i.e. when the pinching pressure is at its maximum value. The comparison of these values is listed in Tab. 3.5 which shows a good match. Similarly, thickness obtained from the simulation using QBATP element was compared with the one using BRICK8 element and also with the one reported by [29] (Fig. 3.24b). The minimum thickness reported by [29] was 3.680 mm which was found to be exactly the same obtained using BRICK8 as well as QBATP elements.

Objective

The goal of this test is to compare the performance of the QBATP element with respect to the BRICK8 element in order to check the impact of the selective mass scaling implemented for the QBATP element.

Setup

The setup of this test is qualitatively identical to the Test 3.4. Same material properties, loading and boundary conditions are used in this test. The only change is in the dimensions of the geometry. The dimensions considered for this case are 800 mm × 800 mm × 4 mm and it was meshed with an element size of 10 mm. The reason to have a much bigger plate is to have a large number of elements for benchmarking without changing the aspect ratio of an element. For a fair comparison, only one element along the thickness was used for the simulation with BRICK8 elements. Thus, in total there are 6400 elements in each of the simulation. The zoomed view of the mesh for BRICK8 and QBATP cases is shown in Fig. 3.25. The simulations were run on a single thread configuration on a local personal computer: Intel(R) Xeon(R) CPU E3-1545M v5 @ 2.90GHz (x86_64).

Discussion

The plot of displacement in X direction (mm) is shown in Fig. 3.26. The summary of the performance of both the cases is tabulated in Tab. 3.6. As discussed before, due to the selective mass scaling, the time-step of the QBATP element is higher than that of the BRICK8 element, thereby reducing the total number of cycles needed for the computation. This coupled with the speedup resulting from the global integration approach, the total speedup was found to be 7.7X which demonstrates the benefit of using the QBATP element from the point of view of the computational cost. thickness (to provide some numerical stability) making its actual thickness to be 2 mm (instead of the full thickness of 4 mm). Thus for a fair comparison (as discussed in Sec. 3.4.1), the real speedup here could be conservatively considered as 7.7/2 = 3.85X which is still significant. This is considering the aspect ratio of in-plane length to thickness L/h was 2.5. The speedup would be higher with a higher ratio of (L/h). Also, this speedup can be further increased as there is a scope for optimization in the subroutines written for QBATP element.

Test 3.6: Clamped square plate Objective

This test is performed in order to check if the QBATP element is able to avoid the locking. For this test, two simulations are compared; one with the classic shell (QBAT) and the second with the pinching shell with selective mass scaling (QBATP).

Setup

A geometry consists of a square plate of dimension 1000 mm × 1000 mm × 2 mm that is clamped (i.e. all displacements and rotations are restricted) on all four edges (AB, BC, CD and AD) as shown in Fig. 3.27. The material properties are: E = 200 GPa and ν = 0.3. The geometry is meshed with 16 elements in each direction. A uniform pressure linearly increasing from 0 Pa up to 800 Pa is applied on all the elements. Two comparisons are to be made as follows:

• The pressure vs deflection curve which is to be compared with the literature

• The deflection of the central node is to be compared with the reference value when the pressure reaches 500 Pa.

Discussion

This being a very thin plate and in the absence of a specific pinching load, the QBATP element should be able to reproduce the behavior of the standard shell (QBAT) without demonstrating locking. The magnitude of the displacement when pressure reaches 800 Pa is shown in Fig. 3.28 where it can be seen that the maximum magnitude of displacement

Test 3.7: Large bending of cantilever

Objective

The goal of this test was to check if the locking is avoided for a large deflection case, similar to the test done in [29]. The large bending type of loading can and does occur in a typical forming scenario.

Setup

The setup is shown in Fig. 3.30a. It consists of a cantilever beam of dimensions 400 mm × 20 mm × 20 mm was fixed at one end (edge AB) and a force of 250 N was applied progressively over 1 s on the other end (edge CD). The material properties are: E = 1000 MPa, ν = 0.3.

Discussion

Due to the applied load, the cantilever undergoes large bending as depicted in Fig. 3.30b showing its deformed shape as the loading increases. The final deformed shape with the plot of displacement field in Z direction is shown in Fig. 3.31 which shows a good match between QBAT and QBATP elements. The evolution of displacement of the tip in Z and X direction using the QBATP element is compared with the QBAT element in Fig. 3.32a and Fig. 3.32b. It shows a good match between the two curves demonstrating that the locking has been avoided even for the large deflection scenario for the QBATP element. The comparison of the final deflection (in Z) of the tip was compared with the values reported in the literature [29,47]. The summary of this comparison is tabulated in Tab. 3.8. It demonstrates that the results of the QBATP element are comparable with the one obtained using QBAT element and also the reference solution.

Test 3.8: Non-uniform pinching pressure

Objective

The goal of this test was to compare the behavior of the QBATP elements with the BRICK8 ones when only a part of the mesh experiences the pinching pressure. This case will often occur in a typical consolidation process where there would always be some regions which undergo higher consolidation than other regions creating a non-uniform pinching pressure distribution.

Setup

The setup is shown in Fig. 3.33 for both BRICK8 and QBATP meshes. The geometrical dimensions are 50 mm × 50 mm × 4 mm. The material properties are: E = 120 GPa and ν = 0.3. This is a quarter model; therefore symmetry conditions are applied as discussed in Test 3.4. A linearly increasing pinching pressure of peak magnitude 10,000 MPa was applied but only on the elements highlighted in red color in Fig. 3.33. There is no pinching pressure on other elements.

Discussion

The displacement field in Y (mm) for both BRICK8 mesh (Fig. 3.34a) and QBATP mesh (Fig. 3.34b) show a good match. The deformed shape has been scaled by a factor of 10X for a better visual representation. Pinching displacements (Z) for nodes N1, N2, N3 and QBATP mesh (Fig. 3.35b) which also show a good correlation. The slight difference between the displacements arises due to the difference of mid-surface normals (which are used for applying the pinching pressure load for QBATP elements) and the top-surface normals (used for applying the pressure load for BRICK8 elements). Comparable results were observed for normal stress field σ xx (Fig. 3.36) demonstrating that the QBATP mesh is able to reproduce the results of a BRICK8 mesh even for a non-uniform pinching load case discussed here.

QBATP element and contacts

In order to use the QBATP element for the simulation of consolidation and/or forming; it is necessary to add the capability of pinching within the contact interface. In this work, this functionality has been added for the contact interface Type-25 in Altair Radioss ™ . A typical main segment of type-25 will have the secondary nodes impacting from both the top (in the direction of the normal of the segment) and bottom (opposite sense of the normal) side of the segment.

A system under equilibrium

The free body diagram of such a system under equilibrium is shown in Fig. 3.37a. Imagine that two equal and opposite compressive forces of magnitude |F | are applied on a system consisting of two secondary nodes S 1 , S 2 and a main segment MM'. The equilibrium for secondary nodes is established as the secondary nodes penetrate into the gap thereby generating a contact force |F c | that opposes the applied load and equal in magnitude (Fig. 3.37b and Fig. 3.37d). The main segment MM' experiences a reaction contact force |F c | due to penetration. This contact force is the cause of the pinching force |F p |. Thus, the main segment is under a through-thickness compressive state experiencing a force of |F c | from both sides (Fig. 3.37c).

Now in case of classic shells (without pinching DoFs); these two opposing forces shall cancel each other out as one calculates the summation of the contributions from the top and bottom side of the main segment. However, in case of pinching shells (or for solids) this cannot be done, as these forces would build a through-thickness compressive stress (pinching stress) which in turn would create an in-plane stretching of the segment. This demonstrates the difference between the contacts for the classic shells versus the pinching shells.

|F upper | = |F 1 | = + 1 2 (|F 1 | -|F 2 |) Mode 1 + 1 2 (|F 1 | + |F 2 |) Mode 2
(3.57)

|F lower | = |F 2 | = - 1 2 (|F 1 | -|F 2 |) Mode 1 + 1 2 (|F 1 | + |F 2 |) Mode 2
(3.58) ). This is the underlying concept for adding the contribution of pinching in the interfaces. This idea when generalized, results in the following global assembly of contact forces for a given main segment:

F contact = F upper + F lower (3.59) F p contact = F n upper - F n lower (3.60)
where F n are the forces projected on the local normal vector n to the contact surface.

With respect to the above equations, note that:

• The forces F upper and F lower are written in the global coordinate system resulting in a positive sign in the expression of F contact and a negative sign for F p contact . • Only the normal components of the contact forces are considered for pinching contact force.

• The case consisting of a system under equilibrium i.e. F upper = -F lower is a special case of the above relationship where there would be no net F contact on the segment but a non-zero F p contact . • The case with no pinching would occur when the main segment moves as a rigid body in the normal direction i.e. F upper = F lower for which F p contact would become zero. With the modified contact interface mechanisms that incorporates the additional contributions of the pinching forces; a validation test was performed which is discussed in the next section.

Test 3.9: Validation of pinching and contacts

Objective

The goal of this test was to compare the behavior of the QBATP elements with the BRICK8 elements when the consolidation pressure is transmitted via contact interfaces.

Setup

The setup consists of three flat plates as shown in Fig. 3.39. Three different models were built. One consisting of BRICK8 elements, the next made up of QBAT elements and the third made up of QBATP elements. The main idea here is to compress the plate in the middle by exerting a pressure on the top plate while the bottom plate remains fixed in all degrees of freedom. This qualitatively mimics the type of loading experienced during consolidation where the pressure applied on the punch generates a through-thickness compression of the prepreg stack.

The top and bottom plates are assigned a very stiff material: E = 2100 GPa, ν = 0.3. On the other hand, the middle plate has been given a relatively softer material in order to observe the in-plane stretching. The material properties used for the middle plate are: E = 4 GPa and ν = 0.495. The initial thickness of the middle plate is 4 mm whereas the in-plane dimensions are 37.5 mm × 37.5 mm. It is meshed with 8 × 8 elements.

A linearly increasing pressure is applied on the top-plate with a maximum magnitude of 800 MPa and it is held constant at this value once it reaches this maximum value. The Type-25 contact is given between the top-middle-bottom plies with a fixed contact stiffness of 500,000 N/mm. Coulomb type of friction with CoF = 0.1 was used. With this setup, all three simulations were run and the results were compared.

Discussion

A first necessary check was to compare the normal contact force in the interface. As expected this matched for all three models (Fig. 3.40). Then the in-plane displacement (X) was compared for all three models. Naturally, in the absence of the pinching DoFs, it is not Figure 3.39: Test 3.9: Geometry, setup, mesh and loading for BRICK8, QBAT and QBATP elements (from left to right) possible for the standard QBAT shell to generate in-plane deformations as a consequence of transverse loading. So no in-plane deformations were observed (Fig. 3.41b). On the other hand, the QBATP element (Fig. 3.41c) demonstrated a good correlation with the BRICK8 model (Fig. 3.41a) for the X displacement field. This indicates that the forces transmitted through the thickness via contacts for QBATP element were able to generate a correct inplane deformations. However, a further investigation was conducted in order to compare the gap between the top and bottom plates vs the thickness of the middle plate. Ideally, one expects that the gap between the top and bottom plates should not be very different from the thickness of the middle plate considering a relatively high contact stiffness was used here (the small difference coming because of penetrations). The comparison of this distance between the In case of BRICK8 element, a good match was observed between d and h. Due to lack of pinching DoFs, for QBAT the thickness cannot decrease but the normal contact force would be transmitted via a penetration. For QBATP element, even though the evolution of h was quite similar to that of a BRICK8 model, the evolution of d did not match with the BRICK8. In fact, it was found to match with the QBAT element.

The reason behind this inconsistency is that, the penetration is calculated with respect to the Gap defined in the interface (Fig. 3.43a). The functionality to update the gap with respect to the change in thickness does not exist in the architecture of the Type-25 interface in Altair Radioss ™ . Thus, the penetration (Pene) is always measured with respect to the initial thickness.

For the classic shells, this feature is not required as (a) the thickness reduction is usually not very big (b) they are not designed for the cases with transverse stress. However, in case of consolidation; this functionality is crucial. Thus, a framework was built inside the Type-25 interface to update the contact gap with the change in thickness dynamically. With this, the penetration would be correctly measured (Fig. 3

.43b).

With this framework implemented, the same test was repeated and compared with the results obtained before. This change should not result in a change in the in-plane displacement field as shown in Fig. 3.44a when compared to the BRICK8 (Fig. 3.44a) as well as the case without dynamic gap update (Fig. 3.44b).

The evolution plot comparing the d and h was then revisited with the new framework

Towards a unified simulation of forming and consolidation

Based on the observations made in the numerical validation tests (Test 3.1 -Test 3.9) performed in previous sections, the QBATP element can be considered to be validated. The next step would be to investigate its usage in an industrial or semi-industrial case.

With the developments conducted in this work, it is evident that the QBATP element with selective mass scaling that has been proposed here has a potential to perform numerical simulation of both forming and consolidation phases in a unified manner. This strategy is discussed in detail in Chapter 5. However, with the developments in this Chapter, one can already perform a qualitative comparison between the QBAT element (Classic shell) and the QBATP element (Pinching shell with selective mass scaling).

Forming and consolidation of a semi-industrial model

For the purpose of the demonstration of the unified simulation of forming and consolidation, a semi-industrial model has been chosen. This mold has been developed by CETIM and will be revisited again in subsequent chapters. The geometry consists of a single sheet of dimensions 400 mm × 120 mm and a thickness of 1 mm. The setup used for the numerical simulation is shown in Fig. 3.46a whereas the real mold geometry and setup is shown in Fig. 3.46b. The mesh consists of elements of size 4 mm. The die and the punch are modeled with rigid elements. The punch was imposed a velocity in Y direction and thereby the distance traveled by the punch is shown in Fig. 3.46c. On the other hand, the die is fixed in all DoFs throughout the simulation. The sheet was modeled with a one shell component representing fibres (Law-58 of Altair Radioss ™ ) (UD fibres oriented along X direction) modeled using QBAT elements. This component was overlapped with the pinching shell QBATP developed in this work which was given elastic material properties with E = 100 MPa and ν = 0.4. It is to be noted that these material properties are used here only for demonstrative purpose. The usage of elastic constitutive model is a limitation which will be addressed in Chapter 4. But overall, this represents a very basic attempt at building a modeling strategy for the unified simulation of forming and consolidation phases. Two simulations were conducted: one using the QBAT element and the other with the QBATP element and the results were compared. The final deformed shape at the end of the consolidation phase is shown in Fig. 3.47. After comparing the width measurements at various locations shown for both QBAT and QBATP simulations; it is evident that there is no increase in the width for the case with QBAT element Fig. 3.47a which is expected as it does not have the capability to transfer the pinching loads to the in-plane directions. On the other hand, the width increments can be observed for the QBATP simulation Fig. 3.47b. To justify that there is a through-thickness consolidation force acting on the sheet; a normal reaction force plot is shown for both QBAT and QBATP simulations in Fig. 3.48a and Fig. 3.48b respectively.

The advantage of using the shell with the capability of transverse stress can then be clearly seen in Fig. 3.49 where it can be seen that obviously the σ zz is zero for the QBAT element (Fig. 3.49a). On the other hand, a non-zero negative value (indicative of a compressive σ zz ) was observed for the QBATP simulation (Fig. 3.49b).

The advantage of using a selective mass scaling approach for the QBATP element can be observed in Fig. 3.50 where it can be noticed that the critical time-step of the QBATP shell

Conclusion

To conclude, in this chapter the existing formulation of QBAT element of Altair Radioss ™ was extended in order to have a functionality of transverse normal (pinching) stress. This implies that the resulting shell element (QBATP) can be in a three dimensional state of stress unlike the classical shells which are under plane-stress condition. The element was validated with various tests and further its interaction with a contact interface was developed, implemented and subsequently validated.

The computational performance of the QBATP element was then augmented with the addition of selective mass scaling that affects specifically the higher frequencies without impacting the lower ones. The increase in the performance was demonstrated with an example which showed about 3.8X speedup when compared to the BRICK8 elements of the same mesh size for an aspect ratio of 2.5.

Finally, it has been shown that the QBATP element serves as a first step towards a unified full-scale simulation of both forming and consolidation phases in a generic composite forming process as well as a more specific one such as QSP ® .

There is however one important limitation concerning the usage of an elastic constitutive model which is debatable both for metal as well as thermoplastic composite forming applications. This serves as a motivation for the next chapter where an elasto-plastic constitutive model compatible with QBATP element is developed.

i=1 ∂N i ∂ξ x i a 2 = dx p dη = 4 i=1 ∂N i ∂η x i (3.62)
x i being the nodal coordinate vectors. The normal to the mid-surface which is a unit vector orthogonal to both the covariant vectors then can be defined as,

n = a 1 × a 2 |a 1 × a 2 | (3.63)
This then creates a mapping [F 0 ] for z = 0,

[F 0 ] = a 1 a 2 n (3.64)
Using this, a general basis at a general location z is then given by,

[F z ] = [F 0 ] + z[F n ] (3.65)
Where the [F z ] mapping is,

dx q (3×1) =   dX dY dZ   = [F z ] (3×3)   dξ dη dζ   (3×1) (3.66)
And [F n ] is given by,

[F n ] = n ,ξ n ,η 0 (3.67)
Where the derivatives n ,ξ and n ,η are calculated as,

n ,ξ = dn dξ = 4 i=1 ∂N i ∂ξ n i n ,η = dn dη = 4 i=1 ∂N i ∂η n i (3.68)
n i being the normal vectors at node i of the element. In order to calculate the inverse mapping

[F 0 ] -1 , first a metric tensor [a] is defined, [a] (2×2) = a 11 a 12 a 21 a 22 = a 1 • a 1 a 1 • a 2 a 2 • a 1 a 2 • a 2 (3.69)
Secondly, the contravariant vectors a 1 and a 2 are such that:

a 1 • a 1 = 1, a 2 • a 2 = 1, a 1 • a 2 = 0, a 2 • a 1 = 0 (3.70)
Using this, the inverse of the matrix [F 0 ] becomes,

[F 0 ] -T = a 1 a 2 n (3.71)
Using eq.( 3.70) and using,

a 1 a 2 = a 1 a 2 [a] -1 (3.72)
the contravariant vectors can be calculated as,

a 1 = 1 det[a] (a 22 a 1 -a 12 a 2 ), a 2 = 1 det[a] (-a 21 a 1 + a 11 a 2 ) (3.73)
With this, the inverse mapping can be written as, dξ dη

(2×1) = (a 1 ) T (a 2 ) T (2×3)   dX dY dZ   (3.74)
The relationship between the global coordinates X, Y and Z and the local coordinate system x, y and z is given by the transformation matrix

[Q],   dX dY dZ   = [Q] (3×3)   dx dy dz   (3.75)
where [Q] is given by,

[Q] = t 1 t 2 n (3.76)
Combining eq.(3.74) and eq.( 3.76), one can build a mapping between the isoparametric coordinates (ξ, η) and the local coordinates (x, y) as,

dξ dη = [C 0 ] (2×2) dx dy (3.77)
where,

[C 0 ] = C 0 11 C 0 12 C 0 21 C 0 22 = a 1 • t 1 a 1 • t 2 a 2 • t 1 a 2 • t 2 (3.78)
The local basis at a point x q (in the shell element but not necessarily on the mid-surface of the shell) can be rewritten as,

[F z ] = [F 0 ] + z[F n ] (3.79) = [F 0 ] ([I] + z[b n ]) (3.80)
where,

[b n ] = [F 0 ] -1 [F n ] =   a 1 • n ,ξ a 1 • n ,η 0 a 2 • n ,ξ a 2 • n ,η 0 0 0 0   (3.81)
With this, the inverse mapping can then be written as,

[F z ] -1 = ([F 0 ] ([I] + z[b n ])) -1 (3.82) = 1 µ(z)     1 0 0 0 1 0 0 0 µ   + z   b 22 n -b 12 n 0 -b 21 n b 11 n 0 0 0 0     (3.83)
where

µ(z) = 1 -2 z H + z 2 K, 2H = -(b 22 n + b 11 n ) and J = det[F z ] = µ(z)J 0 .
For the QBAT as well as for QBATP elements, a linearized version of this relationship is used in order to calculate the strain-rates and strains. Using this relation, another quantity [bc] is defined as,

B matrix for Membrane

The membrane strain-rate is calculated as follows:

  ˙ xx ˙ xx ˙ xy   =   t 1 0 0 t 2 t 2 t 1   [C 0 ] T V ,ξ V ,η = [B m ] (3×24) [V ] (24×1) (3.85)
Where the [B m ] matrix is given by:

[B m ] (3×24) = B m1 (3×6) B m2 (3×6) B m3 (3×6) B m4 (3×6) (3.86) B mi (3×6) =    t T 1 N i,x 0 0 0 t T 2 N i,y 0 0 0 t T 1 N i,y + t T 2 N i,x 0 0 0    (3.87)
where N i,x and N i,y are given by,

N i,x = N i,ξ C 0 11 + N i,η C 0 21 (3.88) N i,y = N i,ξ C 0 12 + N i,η C 0 22 (3.89)

B matrix for Bending

The curvatures are calculated as follows:

  χxx χyy χxy   =   t 1 0 0 t 2 t 2 t 1       [bc] T V ,ξ V ,η + [C 0 ] T β ,ξ β ,η     = [B b ] (3×24) [V ] (24×1) (3.90)
Where the [B b ] matrix is given by:

[B b ] (3×24) = B b1 (3×6) B b2 (3×6) B b3 (3×6) B b4 (3×6) (3.91) B bi (3×6) =    t T 1 N bc 1i s 1i N i,x 0 t T 2 N bc 2i s 2i N i,y 0 t T 1 N bc 2i + t T 2 N bc 1i s 1i N i,y + s 2i N i,x 0    (3.92)
where N bc 1i and N bc 2i are given by,

N bc 1i = N i,ξ bc 11 + N i,η bc 21 (3.93) N bc 2i = N i,ξ bc 12 + N i,η bc 22 (3.94)
and s 1i and s 2i are given by,

s 1i = -t 1 • t 2i t 1 • t 1i (3.95) s 2i = -t 2 • t 2i t 2 • t 1i (3.96)

B matrix for Transverse shear

As discussed in the content of this Chapter, the transverse shear strain-rate is calculated as follows:

˙ xz ˙ yz = [C 0 ] T (2×2) [B ts ] (2×24) [V ] (24×1) (3.97)
where the [C 0 ] is the local basis at z = 0 and the [B ts ] matrix is given by:

[B ts ] (2×24) = B 1 ts (2×6) B 2 ts (2×6) B 3 ts (2×6) B 4 ts (2×6) (3.98) [B 1 ts ] = -nA1 T +AS1 T A1 0 -nB1 T +AS1 T B1 0 [B 2 ts ] = +nA1 T +AS2 T A1 0 -nB2 T +AS1 T B2 0 [B 3 ts ] = +nA2 T +AS2 T A2 0 +nB2 T +AS2 T B2 0 [B 4 ts ] = -nA2 T +AS1 T A2 0 +nB1 T +AS2 T B1 0 (3.99)
Transverse shear strain-rates at the edge mid-points are given by (Fig. 3.51), Figure 3.51: QBATP element: Transverse shear evaluation (reproduced from [25])

γξ = 1 -η 2 γA1 ξ + 1 + η 2 γA2 ξ (3.100) γη = 1 -ξ 2 γB1 η + 1 + ξ 2 γB2 η (3.101) γA1 ξ = a 1 • β + n • V ,ξ ξ=0,η=1 γA2 ξ = a 1 • β + n • V ,ξ ξ=0,η=1 (3.102)
γB1 η = a 2 • β + n • V ,η ξ=-1,η=0 γB2 η = a 2 • β + n • V ,η ξ=1,η=0 (3.103)

Introduction

Through-thickness compression is required to ensure proper consolidation of multi-layered thermoplastic composites. Through-thickness compression force applied during the consolidation of composite prepregs is the driving force for various mechanisms such as squeeze flow [1], in-plane resin bleed [2] and transverse resin percolation [3,4]. As a consequence, several defects can arise during consolidation such as in-plane fibre orientation changes, out of plane fibre waviness [5,6], changes in the in-plane dimensions of a prepreg, all of which have an impact on the performance of the manufactured component [7].

During consolidation, the melt thermoplastic polymer plays an important role because of its near incompressible nature and its ability to flow in the melt state. Because of the incompressibility, as one applies a consolidation pressure transversely on the prepreg stack; in-plane deformations occur. As a result of the deformations and the movement of the thermoplastic polymer coupled with its high viscosity; the fibres are dragged along with it. This phenomenon is known in the literature as the "transverse squeeze flow" [8]. This not only affects the in-plane dimensions of the prepreg as a whole but also affects the spatial field of fibre orientations within the prepreg. The influence of the squeeze flow effect becomes even more important in case of thermoplastic unidirectinoal (UD) prepregs since the in-plane transverse behavior is predominantly governed by the melt polymer. A schematic representation of the transverse squeeze flow phenomenon is shown in Fig. 4.1 for the case of a single UD prepreg undergoing consolidation. This is further discussed in the following sections.

Influence of squeeze flow on in-plane deformations

Depending on the amount of squeeze flow, the in-plane dimensions of the prepreg could experience a significant change. Fig. 4.2a shows the final deformed shape for a UD stack of configuration [0 • ] 5 where the initial undeformed shape is indicated by the green colored rectangular outline. Due to the inextensibility of fibres along X direction and the absence of fibres along Y direction; the major deformations occur along the width of the sample (measured along Y). The initial stack width was 60 mm. After undergoing a 25% consolidation, it was found to have become more than 90 mm which is 1.5 times of the original width.

Even on components manufactured on an industrial scale, the impact of squeeze flow can be seen in Fig. 4.2b. It shows a part manufactured using QSP ® that has a UD patch at the top in its stacking configuration. The fibre directions are schematically indicated by yellow colored lines. The widths are measured at several locations by following curvilinear paths orthogonal to the fibres. The initial width of 7.8 cm was found to have increased significantly at several locations by almost 30% as shown in Fig. 4.2b. Therefore, it is evident that the influence of transverse squeeze flow on the in-plane deformations of the prepregs especially for unidirectional plies is not negligible and hence needs consideration in the numerical simulation of consolidation process modeling. Also, forming of thermoplastic prepregs on doubly-curved parts often induces thickness changes in various zones within prepregs. This implies that thicker zones would experience more consolidation compared to other areas as this change of thickness is usually not anticipated in the mold by adapting the cavity thickness.

Influence of squeeze flow on fibre orientations

The changes in the fibre orientations occurring because of the squeeze flow is important as it influences the mechanical strength and performance of the component. The so called "barrel" shaped deformation pattern resulting from the transverse squeeze flow for a UD (Fig. 4.1) and its influence on the fibre orientations has also been observed and discussed in [9].

In order to investigate the squeeze flow phenomenon further, a consolidation test was conducted on a specimen of a configuration [0 • /90 • ] 6 in [1]. The resulting change in the fibre orientations within the prepreg stack can be clearly seen in Fig. 4.2c. For quantifying the changes in fibre orientations, several copper wires of a very small diameter were placed between the plies that formed the stack of configurations [0 • /20 • /0 • ] and [0 • /80 • /0 • ] in [10]. Changes in the orientations of these copper wires were used as a way to investigate the fibre orientations after consolidation. It was found using micro-tomography that changes in the fibre orientation at certain locations could be as high as 21.5 • . Another experimental investigation conducted in [11] also demonstrated the changes in the fibre-orientations occurring during compression of a prepreg stack. Also the influence of adding interply resin-rich layers on the fibre orientations was investigated. For some configurations, the changes in fibre orientations in plies observed in this work was as high as 22.5 • .

Based on these observations, it can be concluded that in order to correctly simulate the underlying physics of the consolidation process; it is essential to consider the contribution of the near incompressible thermoplastic melt polymer and its interaction with fibres resulting from the squeeze flow behavior.

Constitutive modeling of melt polymer during consolidation

As stated in [12], a thermoplastic polymer behaves as a visco-elastic material at room temperature and visco-plastic at higher temperature. In order to model the behavior of the incompressible thermoplastic melt polymer, various attempts have been made previously. These can be subdivided into two broad categories; one type of approach is followed in the research conducted by the composite forming community whereas a different approach is followed by the consolidation and the rheology community.

Constitutive modeling of melt polymer during forming

An isotropic-viscoelastic constitutive model was used in [13] to model the behavior of the polymer where a membrane element was used for its modeling. Using a micro-mechanical approach and without explicitly modeling fibres and matrix separately; [14] used a constitutive model for prepregs which was based on the viscosity of matrix, weave pattern and initial spacing between fibres. This was compared with results obtained from picture frame tests performed on two prepregs: glass fibre with polypropylene matrix and carbon fibre with epoxy matrix.

Polypropylene polymer was modeled using a Maxwell model in [15] and the dependence of viscosity on strain rate was modeled using a power law. This coupled with a finite element model for fabric was then used to perform various numerical simulations on the Double Dome Benchmark in order to study the influence of material parameters on final deformations and punch force. It was observed that simulation test with viscosity (maximum observed shear angle: 31.88 o ) and without viscosity (maximum observed shear angle: 32.08 o ) resulted in a very minute difference in terms of deformation. The variation in punch force for the two scenarios was found to be 5.1%. Based on this observation the simulation performed in [16] did not consider matrix viscosity.

In the study of thermoset prepreg compression molding done in [17], two sets of overlapping membrane elements with coincident nodes were used where one type of element represents fabric and the other type represents polymer. Neglecting viscous effects, an elastic type of constitutive model was used and the simulations were performed by varying the Young's modulus of polymer and the shear rigidity of fabric.

In [18], the thermoplastic prepreg was modeled as an anisotropic hyperelastic model where the strain energy density function was written as the sum of contributions from fibres, polymer and their interaction. The constitutive model used for polymer was an incompressible Mooney-Rivlin model and no viscous effects were modeled. Since it is difficult to conduct tests on the melt matrix alone; they used the data of bias extension test performed on the prepreg itself to characterize the polymer behavior. It was postulated that the polymer can be considered as the major contributor to the prepreg's shear stiffness upto a moderate deformation (a strain of about 1.5%) mainly due to the fibre undulations present in the fabric. This approach was extended in [19], where the polymer was modeled with a neo-hookean constitutive model along with temperature dependence. Even in this work, the viscosity of melt thermoplastic polymer was neglected.

To summarize, in this broad category of modeling the behavior of melt polymer for forming applications, many researchers have demonstrated that the viscous effects of the matrix can be neglected in the simulation without significantly impacting the final results.

Constitutive modeling of prepregs during consolidation

The consolidation behavior of toughened uncured prepregs was investigated in [20]. These types of prepregs exhibit both types of behavior: they tend to experience squeeze flow (like thermoplastics) and also resin bleed (like thermosets). A hyper-viscoelastic model was proposed in order to model their behavior using a UMAT subroutine in ABAQUS. The experimental work for the same can be found in [21]. The compaction behavior of a woven prepreg (carbon fibre 2×2 twill and PA-66 matrix) was simulated in [22] using a viscoelastic constitutive behavior based on a generalized Maxwell approach where a prismatic solidshell element was used to represent a prepreg.

Even though there have been several attempts at modeling the consolidation behavior for a woven prepreg, the literature for the UD prepregs is predominantly focused on the modeling of squeeze flow phenomenon.

Constitutive modeling of UD prepregs and squeeze flow

The research focused on the modeling of transverse squeeze flow relies on the idea that a thermoplastic UD prepreg behaves as a continuum that can be modeled with a Transversely Isotropic Fluid (TIF) model. The origin of this modeling strategy is in the theory of transversely isotropic fluids proposed by Ericksen [23] in 1960 and the theory for the incompressible and inextensible materials proposed by Pipkin [24] and later formalized by Spencer in [25,26]. These two ideas were combined in [27] to model the squeeze flow behavior of unidirectional prepregs by modeling it as an incompressible and inextensible (in the fibre direction a) transversely isotropic viscous fluid. This modeling approach is known as the 'Ideal Fibre Reinforced fluid Model' (IFRM) approach. The finite element analysis using this technique was performed for with plane stress assumption resulting in a 2D scenario in [28]. Simulations on a full three dimensional model were performed in [1,29] along with the explicit modeling of the resin-rich layer. Following the indicial notation, this constitutive model is given by,

σ ij = -P δ ij + T a i a j + 2µ T d ij + 2(µ L -µ T )(a i a k d kj + a j a k d ki ) (4.1)
where P (pressure) and T (tension in the fibre direction) are the terms arising due to the constraints of incompressibility and inextensibility in the fibre direction respectively. d ij is the strain rate tensor and δ ij is Kronecker delta. The material parameters are µ T and µ L which are known as transverse and longitudinal viscosity which need characterization. Computational cost and obtaining the values for transverse and longitudinal viscosities are two main challenges with this IFRM approach. In order to obtain these values, a semiempirical relationship from Christensen [30] was used in [1] whereas η L = η T was used in [31].

On the choice of constitutive modeling strategy

Based on the literature review, various constitutive models have been used to simulate the behavior of melt polymer. However, there appears to be no standardized approach in its modeling probably due to the complexities and challenges involved in its characterization. Keeping this in mind, several modeling choices are made in this work that not only simplify the numerical aspects of the modeling but also reduce the complexities in the characterization procedure. They are discussed below:

• Following the approach used predominantly by composite forming community; viscous effects are not modeled within the prepreg. A plasticity-based approach is chosen in order to emulate the flow of melt polymer.

• The core idea of modeling the squeeze flow based on the IFRM model is to model the prepreg as an equivalent single continuum without differentiating between polymer and fibres. This automatically incorporates the coupling between the two phases.

In this work, this interaction is achieved by overlapping two shell elements (one for polymer and the other for fibres) with shared and coincident nodes.

The modeling strategy of this work is schematically represented in Fig. 4.3 for a general stack consisting of 4 plies. Each ply contains two shell element components with shared coincident nodes. One of these two components represents fibres whereas the other represents polymer. The shell for fibres is modeled using a full integration standard shell element (QBAT) of Altair Radioss ™ and uses an anisotropic material model (Law 58). On the other hand, the shell for polymer is modeled using the full integration pinching shell element (QBATP) developed in the previous chapter and shall use an elasto-plastic constitutive model (referred to as Law 91) that will be discussed in this work. There are several advantages of following this approach over the other approaches previously discussed in the literature. They are as follows:

• Due to the use of the QBATP element, a full three dimensional state of stress can be simulated while still remaining in the shell element framework. Thus, one can model a through-thickness stress which is the main driving force for the consolidation. This is not possible with the standard shells which are built with a plane stress assumption resulting in a 2D model which cannot capture the squeeze flow behavior correctly [1].

• The approach to use hexahedral elements with IFRM approach is computationally very expensive [29] and is difficult to use for industrial full-scale models involving complex geometries and curvatures. This becomes even more challenging in the case where prepreg patches are used in the stack such as QSP ® process.

• By following the plasticity based approach, the challenging task of determining transverse and longitudinal viscosities in case of an IFRM approach is avoided. This eliminates the need to rely on the semi-empirical relations [1] or to make assumptions about the longitudinal and transverse shear viscosities [31].

• Even though this means the characterization of the material parameters for the plasticity model, as will be discussed in Sec. 4.6, this can be achieved directly using flat areas of the actual mold by conducting a consolidation test on a single UD ply. This avoids the need of a dedicated test setup for the material parameter characterization thereby speeding up the overall design workflow in a typical industrial setting.

• This approach serves as a crucial step in moving towards a unified solution methodology for both forming and consolidation without the need to model both of these processes separately. This aspect is discussed in more detail in Chapter 5.

Objectives and content of the study

The overall goal of this work is to basically develop a strategy to simulate the process of consolidation using a pinching shell (QBATP). With reference to the postulates stated above, the objectives for this work are to:

• Propose an elasto-plastic constitutive model compatible with the pinching shell (QBATP).

• Perform a numerical validation for this constitutive model by comparing with solid element (BRICK8) and also with several other benchmarks in the literature.

• Characterize this model using a consolidation test on a thermoplastic UD prepreg.

• Carry out validation of this model using experiments.

An elasto-plastic constitutive model for QBATP element

A standard mechanism for plasticity in shells involves a through-thickness calculation of stress resultants using several integration points, storing plastic strain at each of them. The other approach is to perform this through-thickness integration analytically and calculate stress resultants directly. This second approach is called as 'global plasticity'.

Since only one point (on the mid-surface) is used to represent the through-thickness state of the element, the usage of global plasticity would mean a faster computation and a reduced memory storage. But it also possesses a drawback that it is not possible to simulate a propagation of plasticity through thickness using global plasticity.

In case of a consolidation process involving a stack of prepregs, the inability of global plasticity algorithm to capture a through-thickness spread of plasticity could be an important phenomenon to capture if the full stack is modeled by just one layer of shell elements across thickness. However, the modeling strategy used in this work differs from this approach. Here, each ply in the prepreg stack is modeled separately as a shell layer and an appropriate contact mechanism is added to model the interaction of plies with one another. Therefore, in this approach, even though the through-thickness spread of plasticity cannot be captured for an individual ply, it can be captured from the point of view of the stack.

Global plasticity for a standard shell

One of the earlier attempts of writing the plastic yield criterion for a shell using directly the stress resultants was proposed in 1948 by Ilyushin [32]. It was written for perfectly plastic thin shells. Thus, neither hardening nor transverse shear stresses were considered in their work. The yield criterion used was,

f = N 2 eq + 1 √ 3α NM eq + 1 α 2 M 2 Eq -(σ 0 y ) 2 = 0 (4.2)
Where σ 0 y is the yield stress; N eq , M eq are the equivalent stress resultants expressed in the form of Von-Mises stress and NM eq is the contribution associated with the coupling between membrane and bending, and α is a modeling parameter. Hardening was included later in the work by Shapiro [33]. The limitation of this stress resultant based approach was that the material yields instantaneously along the thickness. However, in case of a bending induced plasticity, the material on the outer layer of shell crosses the yield stress first and then the plastic zone extends through the layers towards the mid-surface when loaded further. Therefore, there exists a hardening effect before the entire section goes into plastic region. To address this, a pseudo-hardening effect was introduced by [34] by adding another parameter γ defined as a function of equivalent plastic curvature. This work was further extended by [35][36][37][38] but is limited to the case of standard shells.

The development of the pinching shell discussed in [39] is limited to the elastic constitutive model and to the best of our knowledge, there has not been any attempt of modeling global plasticity for the pinching shell which was proposed in [39].

Global plasticity for a pinching shell

The elasto-plastic constitutive model developed in this work has 4 material parameters: Young's modulus (E), Poisson's ratio (ν), Yield stress (σ y ) and linear isotropic hardening modulus (HM). This constitutive model (referred to as Law 91 in Altair Radioss ™ material model library) has been developed uniquely for the pinching shell (QBATP). For this shell, the stress resultants are N x , N y , N z , N xy , N yz and N xz whereas the stress couples are M x , M y and M xy . For generalizing one can define,

Ñi = N i h = 1 h h/2 -h/2 σ i dz Mi = M i h 2 = 1 h 2 h/2 -h/2 σ i zdz (4.3)
Thus, the term Ñi can be interpreted as a through-thickness average of stress σ i . This strategy of generalization of stress resultants has been followed in [36] for standard shells and has also been used in the software architecture of Altair Radioss ™ [40].

For the sake of simplicity, the following assumptions are made in the development of this model.

• For the criterion, the coupling term of the membrane and bending part is not considered. For the locations that are undergoing consolidation, it is a justifiable to assume that they would be predominantly under a membrane (non-bending) state of stress rather than a bending state of stress thereby justifying the choice of neglecting the coupling term.

• The contribution of transverse shear is also not considered in the criterion. Thus, the resultants N yz and N xz are always updated assuming elastic behavior. This is not a strong assumption, since for thin shells it has been shown in [38] that the results obtained by considering elastic response for transverse shear showed no significant difference when plasticity was considered for transverse shear.

With these assumptions, the yield criterion for the global plasticity has been modified so as to be usable for the pinching shell. Before considering the criterion for the pinching shell, it is worth to revisit the origin of the global plasticity criterion for standard shells. For the case of J2 (Von-Mises) plasticity for a standard shell, the yield surface is given by the Von Mises equivalent stress given by its standard definition in terms of the principal stresses σ 1 , σ 2 and σ 3 ,

σ 2 v = 3J 2 = (σ 1 -σ 2 ) 2 + (σ 2 -σ 3 ) 2 + (σ 3 -σ 1 ) 2 2 (4.4)
Calculation of the Von Mises stress for the membrane and bending part of the stress separately before adding them would reproduce the uncoupled version of eq. 4.2 with α = 0.25.

Owing to the possible 3D state of stress that can occur in the case of QBATP element, σ zz = 0 is not imposed. Thus, the eigenvalues of the membrane (non-bending) part are calculated considering the 3D state of stress and using the stress resultants,

Algorithm

Using the modified yield criterion for pinching shell, the step by step algorithm is given below:

1. Calculate the deviatoric part of resultant of normal stresses { Ñ n } using,

P n = 1 3 ( Ñ n xx + Ñ n yy + Ñ n zz ) { Ñ n dev } =        Ñ dev x Ñ dev y Ñ dev z        =        Ñ n x + P n Ñ n y + P n Ñ n z + P n        (4.

Test 4.1: Displacement driven membrane loading

Objective

The objective of this test is to verify the elasto-plastic response of QBATP element by comparing it with BRICK8 element of Altair Radioss ™ when it undergoes displacement driven membrane loading. In this case, σ zz ≈ 0 whereas σ xx = 0.

Setup

A single element of dimensions 10 × 10 × 0.5 mm is stretched along X axis using linearly increasing imposed displacement (Fig. 4.4). The Young's modulus (E) was set to 2000 MPa and a Poisson's ratio (ν) of 0.3 was used. The yield stress (σ y ) is set to 100 MPa. Three simulations were performed by varying the hardening modulus (HM) as 0, 500 and 1000 MPa.

From Altair Radioss ™ element library, BRICK8 element with an appropriate elasto-plastic model was used as a reference. The resulting nodal displacements (N3 of QBATP element) in Y and Z (pinching) direction were compared. Also, a comparison was made for the stress in X direction and the evolution of equivalent plastic strain eq p .

Discussion

It can be seen that a good match was observed between the BRICK8 and QBATP element both in pre-yielding as well as post-yielding phases in terms of membrane stress σ xx (Fig. 4.5a) and the equivalent plastic strain (Fig. 4.5b) for all three tests. Also, the comparison of nodal displacement for N3 demonstrated a good match for both Y (Fig. 4.6a) and Z (pinching) direction (Fig. 4.6b).

Setup

The same element from Test 4.1 was used in this case with same material properties. This being a load driven test instead of a displacement driven one; the case with perfect plasticity (HM= 0) couldn't be simulated as it would have resulted in an instability. The values of hardening modulus (HM) used were 200, 500 and 1000 MPa. A linearly increasing pinching pressure that starts at 0 MPa and reaches a maximum value of 200 MPa was applied. This is more than the yield stress (σ y ) of the material which was set at 100 MPa. The side view of the element is shown in Fig. 4.7).

Discussion

The test being symmetric in terms of X and Y directions, the displacements fields are identical in both X and Y. The displacements in membrane directions predicted by the QBATP element perfectly matched with those obtained from BRICK8 element (Fig. 4.8a). Similar observation could be made about the pinching displacements of the QBATP element both before and after yielding (Fig. 4.8b).

This test being similar to a uniaxial (σ xx ≈ 0, σ yy ≈ 0) compression test, the slopes of the linear hardening curves in the pinching stress (σ zz ) vs pinching strain ( zz ) plot (Fig. 4.9a) were found to be equal to E•HM E+HM as expected. Also the predictions of equivalent plastic strain ( eq p ) obtained using QBATP matched correctly with those obtained from BRICK8 (Fig. 4.9b).

Objective

The objective of this test is to test a multi-axial type of loading where both the pinching pressure and membrane tension act simultaneously creating a multi-axial state of stress within the element.

Setup

The geometry of the element was set to 10 × 10 × 1.0 mm and E = 2000 MPa, ν = 0.3, σ y = 100 MPa with HM = 1000 MPa were used as material parameters. A linearly increasing pinching pressure that starts at 0 MPa and reaches a maximum value of 500 MPa was applied along with a linearly increasing nodal force of 200 N. The side view of the element is shown in Fig. 4.10).

Objective

The objective of this test is to perform uniaxial in-plane stretching of a rectangular plate with a circular hole in the center. This test is inspired from a similar test from the literature discussed in [37,38].

Setup

The rectangular plate is of size l = 36 mm, b = 20 mm with a circular perforation of radius r = 5 mm at its center. Different simulations were performed by varying thickness of the plate (h = 0.1, 0.5, 1.0 and 1.5 mm) for both QBAT and QBATP elements. E = 70 MPa, ν = 0.0, σ y = 0.243 MPa with HM = 0.2 MPa were used as material parameters. Only one quarter of the plate was modeled due to symmetry and subsequently appropriate symmetry boundary conditions were applied on edges AB and DE (shown in Fig. 4.14a). An imposed velocity in Y direction was applied along edge CD and reaction force (Y) vs displacement (Y) curves were compared for these simulations (Fig. 4.14b).

Discussion

In this case, since σ zz ≈ 0, it is expected that the QBATP element should give almost identical results as that of QBAT which can be seen in Fig. 4.14b and also with the field of σ yy shown in Fig. 4.15. Besides, the force needed for yielding calculated using F yield = σ y (b -2r)h is also plotted which was found to be matching with the instant at which yielding begins. The objective is to perform the benchmark test popularly known as the Cook's membrane test which is designed to investigate the volumetric locking phenomenon. It has been used in the literature for plane-strain shells [41] as well as solid-shells [42] and hexahedral (solid) elements [43][44][45] with some modifications. Here it is used to demonstrate the influence of performing selective reduced integration on the solution convergence in case of a near incompressible material.

Setup

The geometry consists of a tapered plate of thickness 1 mm (Fig. 4.16a). Edge AD is fixed in all DoFs whereas edge BC is subjected to a force in Y direction. The material properties used are: E= 200 MPa and the Poisson's ratio is varied from a case of typical elasticity (ν=0.3) to the cases of near incompressibility (ν=0.495). Simulations are performed with mapped meshes as: 2×2, 4×4, 6×6, 8×8, 12×12, 16×16, 24×24 and finally 32×32. Two sets of simulations were performed; one using the selective reduced integration (abbreviated here as SRI) and other without SRI.

Discussion

The deformed shape for the case of ν=0.495 is depicted in Fig. 4.16b. To investigate the convergence behavior of the model and formulations; vertical displacement (Y) of point C was plotted against the number of elements per side and the curves of without and with SRI are compared as ν approaches 0.5 (full incompressibility). It can be seen that performing a selective reduced integration (with SRI) didn't improve the convergence behavior very significantly compared to the full-integration of spherical component of stress (without SRI). However, as the material becomes nearly incompressible (ν → 0.5) the convergence behavior of the case without SRI gets worse due to volumetric locking. On the other hand, the convergence of the case with selective reduced integration can be seen to be much better.

Test 4.6: Cylinder under shear loading

Objective

The objective of this test is to reproduce a benchmark test from [38] consisting of a hollow cylindrical object with elastic perfectly plastic material properties undergoing a large shear.

Setup

Discussion

Reaction force in Y direction obtained from the simulation was plotted against the Y displacement of the top edge (Fig. 4.18b). It was found that the response of QBATP element was quite close to the response obtained by [38]. It was also found to be very similar in response to the ABAQUS result involving S4R elements (curve reproduced from [38]).

Objective

The objective is to test QBATP element with the present constitutive model for the benchmark case of a pinched cylinder with linear isotropic hardening. This was proposed in [35,46] and later used in [38].

Setup

The geometry consists of a hollow cylinder with a length of 600 cm, radius of 300 cm and a thickness of 3 cm (Fig. 4.21a). The material properties assigned to the cylinder were: E = 3000 kN/cm 2 , ν = 0.3, σ y = 24.3 kN/cm 2 and the linear isotropic hardening modulus (HM) was set to 300 kN/cm 2 . Both ends of the cylinder are constrained by a rigid diaphragm. Hence, a boundary condition that restricts u x , u y and r z was applied at the ends of the cylinder. Two opposing point forces (radial direction) were applied at locations A and B creating a pinching action. The displacement of point B was measured and plotted against the applied load (Fig. 4.21b).

Discussion

The results obtained using QBATP element with the present constitutive model were similar to the ones obtained in [38] which is evident from the Fig. 4.21b. Some deviation is observed when the displacement exceeds 250 cm, but this could have been caused by the snap-through behavior observed by Dujc in [38].

(a) (b) (c) (d) (e) (f)

Materials

The characterization in this work is done solely on UD plies. The UD ply tapes used in this work consist of Glass fibres impregnated with PA-66 polymer. The specification for the same is given by Celestran ® CFR-TP PA-66 GF60-02 (Supplier: TICONA).

Obtaining Young's modulus and Poisson's ratio

To obtain these parameters, first the bulk modulus of the PA-66 polymer was estimated using the Tait equation which is given by,

V p (T, P ) = V p (T, 0) [1 -C(T )] ln 1 + P B(T ) (4.23)
where V p is the specific volume of polymer, T is the temperature, P is the pressure and other variables (B(T ) and C(T )) are material parameters. The material parameters for PA-66 polymer tabulated in Table I of [47] were used . The parameters are: C(T ) = 0.0894, V p (T, 0) = 0.7657e 0.00066T and B(T ) = 3226e -0.00504E-3T . Further, by assuming a linear equation of state and assuming an isothermal process; the bulk modulus of the polymer (K) was found to be 226.46 MPa.

In future work, these assumptions can be resolved by implementing a specific and temperature dependent equation of state. The incompressibility of melt polymer is the motivating factor to use Poisson's ratio (ν) as 0.495. Using, the bulk modulus and Poisson's ratio, the Young's modulus was calculated using E = 3K(1 -2ν) and it was found to be 6.8 MPa.

Obtaining Yield stress and Hardening modulus

To obtain the remaining two material parameters yield stress (σ y ) and hardening modulus (HM), a UD consolidation test is proposed here in which the UD stack is squeezed in a flat mold at a high temperature equivalent to QSP ® . The final deformed shape obtained from the process is used to characterize these two parameters.

The key idea here is to use the flat sections of the actual mold to characterize the constitutive model avoiding a dedicated setup for characterization. This point makes it highly suitable for industrial use where a dedicated characterization campaign is often difficult to conduct as it can impact the design cycle time and cost significantly.

It is worth mentioning that an approach to use a Digital Image Correlation technique was attempted by marking an initial grid on the specimen so that the local deformations could be studied. However, due to the squeeze flow occurring at such a high temperature, the markings were spread thereby making it impossible to study local deformations with this approach.

The workflow for the characterization procedure is shown in Fig. 4.24. The stack of UD prepregs C undergoes the consolidation in the flat sections of the mold B . A schematic representation of the temperature and loading cycles is shown in A . After the completion of the consolidation experiment, one can measure the final average thickness D and final maximum width F using the deformed shape of the prepreg.

Due to the use of the industrial hydraulic press, this is representative of industrial conditions. Thus neither the accurate compaction force nor the punch displacement can be controlled during the test. This also is a typical situation in the industrial setups where these quantities cannot be controlled or measured very precisely. Thus, the characterization strategy was planned solely based on the final deformed shape after consolidation.

The final deformed shape can be subdivided into two quantities: out of plane deformation (thickness reduction) and in-plane deformation (increase of width and length due to the transverse squeeze flow). These two are the only outputs needed for the characterization method proposed here.

Experimental procedure

The experimental campaign on a stack of UD prepregs with a configuration [0 • ] 5 made up of the material specified in Sec. 4.6.1 and dimensions 60 × 60. The total initial thickness was measured to be 1.55 mm. The experiments were conducted using the Omega mold provided by CETIM which was installed in Centrale Nantes. The setup of the installed mold is shown in Fig. 4.25 where the flat sections of the mold that were used for the experiments are highlighted by blue color. The Scamex machine of 750 kN capacity was used to conduct the experiments. Only one sample was tested for this campaign. Metal spacers of appropriate thicknesses were used as a way to restrict the maximum amount of consolidation.

Determining Yield stress and Hardening modulus

Referring back to Fig. 4.24, the value of the final average thickness D serves as an input to the simulation. Also, the material parameters obtained using Tait equation and incompressibility condition (E and ν) are input to the simulation G . Using these inputs, a set of simulations E were performed by varying the values for σ y and HM for the material H . The final maximum width obtained from the simulation (W sim ) I was compared with the final maximum width obtained from the experiments (W exp ) F in the block J . This difference (W sim -W exp ) served as a quantitative measure in order to perform the simulations with the next guess of σ y and HM H . This feedback loop is indicated by the red colored arrows in Fig. 4.24. A set of simulations consisting of 25 simulations with a range of values for the parameters σ y =[0.05, 0.1, 0.25, 0.5, 1] MPa and HM = [0.05, 0.1, 0.25, 0.5, 1] MPa were performed. Based on the obtained results and comparing the quantity (W sim -W exp ), a further refinement in the material properties was performed. Thus, the set of values for which the difference (W sim -W exp ) was minimum were chosen as the material parameters.

In future, a dedicated optimization framework can be used in order to obtain more precise values of the material parameters and also multiple experimental repetitions should be done in order to identify the sensitivity of the experiments and thereby the variation of the W exp .

Die/Punch (Law- The final deformed shape obtained from the simulation is compared to the experiment in Fig. 4.27. Fibres being nearly inextensible, the displacements along X are very small as indicated by the measurements at various locations. It can be observed that even though some curvature is seen in the simulation depicting the barreling effect, it does not match exactly with the experiment. This difference could have come from the ply-mold interaction and the heat transfer happening at their interface which is not modeled in the simulation. The in-plane deformations observed from the simulations are compared with those obtained from the experiments in Fig. 4.28. The width predictions at the corners (locations A , C and E in Fig. 4.28) are predicted quite well from the simulation (with a difference of less than 1 mm). Also, the deformation patterns at the corners are predicted with a good accuracy. However, the amount of widening predicted by the simulation at the midpoints of sides of the stack (locations B and D Fig. 4.28) is less than what has been observed in the experiments. There are several possible hypotheses that might explain this difference as follows:

• Since the simulation is performed with the assumption of an isothermal process, the ply-mold interaction and specifically the heat transfer happening at their interface is not modeled in the simulation.

• In the current experimental setup, a uniform temperature distribution at all points on the mold surfaces could not be ensured due to the existing design limitations. This can result in a temperature inhomogeneity that can impact the squeeze flow behavior observed in the experiments.

• In this work, the Stribeck theory for modeling the tangential viscous friction [48] is not considered. Therefore, the influence of the normal pressure on the tangential sliding behavior is not taken into account here.

• The physics resulting from the phase changes (crystallization) in polymer is not modeled in this work.

Thus, even though the overall in-plane deformations of the prepreg stack were captured both qualitatively and quantitatively, a further investigations both from the point of view of numerical and experimental work can be helpful in identifying and resolving the differences observed in the results of the simulation and the experiment.

Conclusion

In this work, an elasto-plastic constitutive model based on the global plasticity approach has been developed in order to model the behavior of the thermoplastic polymer during the consolidation process. It is designed to be compatible with the pinching shell element (QBATP) so that together they provide a shell based modeling solution for the consolidation process modeling of thermoplastic prepregs which has not been done before.

Additionally, a methodology for the characterization of this constitutive model has been proposed that uses the flat sections of the actual mold thereby avoiding the need to have a dedicated test setup. This makes it highly suitable for full-scale industrial applications where design cycle time is a very critical factor to be considered.

Finally, the characterized constitutive model was used to simulate a relatively simple yet essential stack configuration of [0 • /90 • ] 2 . Overall, the results obtained from the simulation showed a good correlation with the experimental observations for this prepreg stack after undergoing a consolidation process.

• Filling small gaps between the forming prepregs and the mold-walls/edges.

Therefore, it plays an important role in the overall final quality and performance of the component. Also, consolidation has an influence on the spring-back of components in the demoulding stage indicating the importance of including the numerical modeling of consolidation in the process simulation of thermoplastic composite forming. The state of the art of process modeling of thermoplastic composite forming from an industrial point of view has been presented in [2]. Various features and capabilities of different industrial solutions such as PAM-Form, AniForm, LS-Dyna and ABAQUS were discussed. However, it can be observed that the main focus of various industrial simulation software is related to the development of new capabilities and improving existing ones predominantly for the forming phase of the process.

The research work of Soulat et al. [3] did address the need to have a unified solution using a shell based approach that comprises of both forming and consolidation phases. The pinching shell element developed in their work was implemented in the commercial code of PAM-Form software. However, only the elastic material model was discussed in their work which cannot be used readily for the simulation of prepregs.

On the other hand, the numerical modeling of consolidation process conducted in [4] is performed on a relatively simple (flat) geometry and not on an industrial model. A relatively more complex part involving multiple plies was used in [5] where the simulations were carried out in ABAQUS/Standard with the addition of specific user material models. In both of these works, each individual ply was modeled using a Hexahedral element (C3D8 element of ABAQUS).

A potential approach for providing a unified solution is proposed in [6][7][8][9] with the use of a prismatic solid-shell element. However, to the best of our knowledge, there exists no unique numerical solution that is able to address both the forming and consolidation phases (with or without long distance sliding of prepreg patches) together in a unified manner that is solely based on shell elements. This is the main focus of this Chapter.

Objectives and content of the study

The overall goal of this Chapter is to demonstrate the strategy that combines the use of interply adhesion, QBATP element and elasto-plastic constitutive model for melt polymer to provide an unified solution for both forming and consolidation phases of QSP ® . The main objectives of this chapter are to:

• Discuss the two-step approach where the simulation of the forming and consolidation phases is achieved sequentially.

• Discuss and demonstrate the challenges associated with the two-step approach.

• Discuss the strategy for a unified solution purely based on shells and its advantages over the two-step approach.

• Perform a full-scale numerical simulation on an industrial component developed by CETIM.

Two-step approach

Principle

One of the ways to address the simulation of both forming and consolidation phases is the two-step approach. This two-step procedure is motivated by the fact that the primary mechanisms for forming phase and consolidation phase are quite different. Thus the core idea of a two-step procedure of simulation of thermoplastic composite forming is to perform the simulation of forming phase and consolidation phase separately in a serial manner. The flowchart for this type of procedure is shown in Fig. 5.2.

Step The sub-steps for each of the forming and consolidation phases are given as follows:

• Forming Phase Forming phase setup starts with a setup of a 2D prepreg stack where each ply is modeled as a shell element component with an appropriate contact model between them. The setup phase is complete after assigning other geometrical, material, boundary conditions and the process parameters. After that, the actual simulation for the forming phase begins where the stack is deformed under the forming force exerted by the punch. The stack is deformed until it takes the shape of the mold. The next substep is to Export deformed shape of all elements into a mesh format. This exported mesh serves as the input to the consolidation phase.

• Consolidation Phase

The setup for the consolidation phase begins by importing the deformed shape information generated by the forming phase. In order to simulate a through-thickness compressive behavior, it is essential to add a HEX elements layer between two shell components representing the melt polymer. After this, the setup needs to be modified based on this updated geometry and mesh. This involves redefining the contact algorithms, boundary conditions, process parameters and material parameters. Once this step is complete, the simulation of the consolidation phase can be performed.

Challenges with the two-step Strategy

The main advantage of a two-step approach is that, it can use the existing solvers for forming and consolidation phases as long as one is able to write converters to read and write data in a desired format. This methods works well for academic or relatively simple cases. However, when it comes to industrial or even semi-industrial models (such as the one shown in Fig. 5.3), there are several challenges associated with this two-step approach.

Mesh quality of HEX layer

At the end of the forming stage, the shell components are already in a 3D deformed configuration. Therefore, it is already challenging to introduce a HEX elements layer in between the two shell layers. In an ideal scenario shown in Fig. 5.4a, one can expect a more or less uniform gap between two shell components. However, in reality, certain regions with a very small gap were observed. This is schematically represented in Fig. 5.4b. This is often the case with the industrial components which have more complex geometry and curvatures.

This process of adding an HEX elements layer was attempted for the industrial model from company Dallara shown in Fig. 5.3. Despite using the advanced meshing capabilities of Hypermesh; the mesh quality is not good enough in several regions. One such region is shown in Fig. 5.4c where the green colored surface shown the top shell layer and yellow colored HEX elements are the ones created using the deformed mesh of top and bottom shell components. Clearly, the small gaps resulting from the deformations of top and bottom shell components cause very small thicknesses of the HEX elements at these locations. When the consolidation simulation was performed with this mesh, it resulted in a very small time-step and eventually stopped due to excessive deformations at certain locations.

Possibility of intersections

During forming phase, Type25 contact law is used to model ply-ply interaction. For this contact mechanism, the penalty stiffness is constant (i.e. it does not depend on the penetration) [10]. Thus, with the use of Type25 contact law, it is possible to get intersections which are defined when the mid-surface of a shell crosses the mid-surface of other shell. Despite the mid-surface crossing the other mid-surface, there are penalty forces which are oriented based on the direction from which the penetration/intersection has occurred. This is schematically represented in Fig. 5.5 showing the difference between a penetration and intersections. High Penetrations would result in a very small thicknesses of the HEX elements between the shell components; whereas an intersection would mean that it is not even possible to create the HEX element layer at that location. If the intersecting nodes are far away form the region of interest and are very few in number, these nodes can be moved to remove intersections but keeping high penetrations. However, this is very subjective to where and how many intersections are observed and hence cannot be considered as a robust solution to the problem.

Challenges in automation

Because both high penetrations and intersections need a user intervention to decide about their resolution, it is a challenging task to make the whole process automated. From an industrial point of view, this is a big hurdle where the goal is to provide a one-click solution with maximum amount of automation and minimum human intervention. Also, the adding of an HEX element layer between two deformed shell components is a tedious task whose complexity will increase as the number of plies increase.

Possibility of local consolidation even during forming phase

There is however, one more challenge associated with this two-step approach which is about deciding when the forming phase is complete and when the consolidation phase shall begin.

As it is possible that porosities in some regions are removed during the forming phase indicating a possibility of a local consolidation even during the forming phase. This serves as an evidence that certain regions of the component can begin undergoing consolidation before other regions, therefore a global split of forming and consolidation phase may not be possible.

In conclusion, based on these observed challenges and limitations; a unified approach is more desirable compared to a two-step strategy.

Towards a unified solution for a full scale simulation of QSP ®

The idea behind a unified approach is to build a modeling strategy that can be used for both forming as well as consolidation phases without the need to explicitly separate them.

A schematic representation for comparing the two-step approach with a unified approach is shown in Fig. 5.7. The unified approach in general simplifies the overall workflow and avoids the major limitations of a two-step workflow.

Ingredients of the unified approach

The developments of the previous chapters can be combined to build this unified approach. The global overview of how these individual developments serve as the building blocks of the unified approach is depicted in Fig. 5.8.

To recall from Chapter 4, the modeling strategy is schematically represented in Fig. 5

Benefits of the unified approach

The benefits of the unified approach over the two-step approach are as follows:

• One does not have to create the HEX layer mesh after forming stage at all. If the elements used right from the beginning are capable of providing a through-thickness stress such as QBATP element, they will be able to demonstrate a consolidation behavior. This avoids the limitations that arise due to excessive penetration, intersections as there is no more the need to build a HEX layer in between.

• From user point of view, there is no intervention necessary at the intermediate steps such as import/export of data, HEX mesh creation, resolving mesh quality related issues and so on. This makes it a highly suitable solution from an industrial point of view.

• With the unified approach, one does not have to make a decision to decide when the forming phase ends and consolidation phase begins. Thus, this permits the possibility to have a local consolidation in certain regions even when the stack is still undergoing forming.

With the proposed modeling strategy and the unified approach, one can now perform fullscale numerical simulations of a typical composite forming process and also for specialized processes such as QSP ® .

Revisiting the industrial model from CETIM

In this section, the industrial model from CETIM discussed previously in Chapter 2 is revisited but with the updated modeling strategy that consists of using the QBATP element with the elasto-plastic constitutive model.

Model description and setup

To recall, the considered industrial model (Fig. 5.10a) has a prepreg stack consisting of 9 plies of different shapes, fibre architectures (UD/woven) and fibre orientations (Fig. 5.10b).

The woven prepreg plies are TEPEX ® dynalite 101-RG600(x) (Supplier: BONDLAMINATES).

The UD ply tapes are Celestran ® CFR-TP PA-66 GF60-02 (Supplier: TICONA). Similar to the discussion in the Chapter 2, a mesh size of 4 mm was used which generated a total 29,653 nodes consisting of 48,404 quadrilateral elements. The material parameters and parameters related to the interply adhesion are summarized in Tab. 5.1.

Results and discussion

Overall, with reference to the final deformed shape, the results of simulation with unified approach (Fig. 5.11b) are comparable with the ones obtained from the experiment (Fig. 5.11a). However, now with this new unified approach, one has access to additional information such as through-thickness normal stress. The maximum initial stack thickness is 3.5 mm. Thus, a rough estimate can be made to predict that the consolidation phase will begin when the distance between the punch and die reduces to 3.5 mm (after adjusting for the thicknesses of die and punch themselves).

Based on this prediction, the transverse normal stress for all plies was plotted along with the magnitude of contact force experienced by the punch at various instants during the simulation. Fig. 5.12 is obtained when the punch-die distance is 4.78 mm which is before the beginning of the consolidation as per estimate. As expected, there is neither a noticeable σ zz nor any contact force on the flat portions of the punch. It can be seen that the maximum compressive stress is on the flat section of the top face of the component. Whereas, other regions are contributing less in the sharing of the consolidation pressure. This can be qualitatively cross-verified with the observed contact force plot. Ply-wise through-thickness normal stress is an important information which can be used in the design of the mold and/or stack. This helps in identifying the regions which will have a better consolidation and which regions will have a poorer quality of consolidation. The regions with lower consolidation are likely to have higher porosities and therefore are possible sites of weakness/failure for the component. Thus, based on the information obtained from this type of simulation, one can either modify the mold design or the stack itself, in order to have a better consolidation behavior.

The transverse normal pinching stress (σ zz ) at three locations on Ply-1 was tracked throughout the simulation which are shown in Fig. 5.15a. One location was chosen at corner indicated by marker A , the second on the top face indicated by marker B and third was chosen on the front face indicated by marker C . The evolution of σ zz at these locations is shown in Fig. 5.15b. It can be seen that the consolidation at the corner started much before the location on top and front face. This demonstrates the advantage of using the unified approach where local consolidation occurs even during the forming stage can be simulated. Similar to Chapter 2, the fibre orientations after the simulations using the unified approach were compared at various locations as shown in Fig. 5.16 with those obtained from the experiments. The fibre orientations obtained from the unified simulation approach are listed in Tab. 5.2. To recall, the nomenclature for the table is as follows: A data-point (at which fibre orientation is measured) is identified as PiDj where i denotes the ply number and j denotes the family of fibres (1 or 2). The column Loc contains the identifier for the measurement location of plies (Refer to Fig. 5 It is worth comparing the fibre orientations of the unified approach with both the experimental values as well as the values obtained from the simulation of a purely forming based approach. Fig. 5.17 to Fig. 5.21 show the graphical comparison for each face in order to compare these three data-sets: experiment, pure forming approach, unified approach. Based on the comparisons, following observations can be made:

• Referring to Face-1 data (Fig. 5.17), overall there is no noticeable difference between the forming versus the unified approach. However at data-points 1.4-P9D1 and 1.3-P9D1 (which are the locations on the UD prepreg patches); the predictions of the unified approach are closer to the experimental observations.

• Referring to Face-2 data (Fig. 5.18), again the results of the unified solution are almost similar to purely forming case. But some observable improvement in the prediction can be found at data-points 2.1-P6D1 and 2.1-P5D1. But also there are data-points such as 2.2-P2D1, 2.2-P4D2, 2.4-P2D1 and 2.4-P3D1 (which are data-points on the woven prepregs) where in fact the forming solution is closer to the experimental observation.

• For Face-3 data (Fig. 5.19), significant improvement is observed at location 3.4 where the pure forming solution was quite different from the experimental observations. At this location, even though the unified solution is still not within the bounds of the experimental observations; the predictions are much better compared to the pure forming approach. Also, improvements can be observed at several other data-points on this face, reducing the overall average error on this face to 7.1 • compared to forming approach where it was 9.8 • .

• For Face-4 (Fig. 5.20), with the exception for data-points 4.3-P2D2 and 4.4-P1D2 (which are the locations on the woven prepregs), the unified solution has provided a similar solution as that of the forming solution.

• For faces 5 and 6 data (Fig. 5.21), a significant improvement is observed at data-point 6-P9D1 which is located on the topmost UD patch. The results summarized with respect to each face are presented in Tab. 5.3. The error in this table is calculated by taking an average of the absolute difference between the values from the simulation and the mean values from experimental measurements. Overall when compared for all faces, the fibre orientations predicted by the unified simulation were on an average 5.4 • away from the experiments. This value is similar to the value obtained with a pure forming solution (which was 5.9 • degrees).

To conclude, this section has demonstrated that the fibre orientations are still in alignment with the experimental observations even with the unified approach. This verification is essential since there have been several changes in the modeling approach for the forming approach versus the unified approach such as the usage of QBATP elements with elastoplastic constitutive model and other consequences of the pinching behavior both within the components as well as at the interfaces.

Even though there are less data points available for UD compared to the woven prepregs; as a general trend, it has been observed that the predictions for the locations on the UD prepregs seem to have improved. The locations where the error increased are predominantly located on the woven prepregs, which might indicate a need to further investigate the applicability of the constitutive model for the polymer in case of a woven prepregs. It is possible that there needs some further development and/or modification to the existing constitutive model in order to correctly simulate the woven prepreg behavior.

It is perhaps expected that in this component, there is not much of a squeeze flow to be observed because of mainly two reasons. First, the compressive loading predominantly passed through a zone on the top face (Fig. 5.14). Second, due to the overall stack configuration, there may not be enough room for the squeeze flow to occur. This may not be the case with other components and stacking sequences. An example designed to show the squeeze flow on a semi-industrial component is discussed next.

Numerical simulation of a UD ply

In order to demonstrate the impact of squeeze flow phenomenon on the overall deformation, the unified approach described in this chapter has been applied to a semi-industrial Omega mold provided by CETIM.

Model setup

The geometry consists of a single UD ply of dimensions 400 mm × 120 mm and thickness of 1 mm. The setup simulation is shown in Fig. 5.22a. The fibres are oriented along the X direction. The UD ply was meshed with an element size of 4 mm. The punch undergoes an imposed displacement as shown in Fig. 5.22b first performing the forming stage and subsequently followed by the consolidation phase. The unified simulation starts with the initial position (Fig. 5.23a). As the punch begins to move downwards 1 (in -Y direction), it starts deforming the UD ply (Fig. 5.23b). As the distance between punch and die reduces further 2 -3 , the forming continues (Fig. 5.23c).

Results and discussion

As the gap between die-punch (after adjusting for their thicknesses) approaches the ply thickness 4 ; it can be seen that the UD ply has been deformed to take the shape of the mold (Fig. 5.23d). After this instant, a further downward movement of the punch would As the consolidation begins 5 ; one starts observing the effects of squeeze flow which can be seen by the increase of the width (Fig. 5.23e). As the punch travels even further 6 ; the width increases even more (Fig. 5.23f). In this case, the maximum width increase is about 14.3%.

Consolidation in localized regions

An interesting observation can be made about the regions where squeeze flow phenomenon was observed. It is evident from Fig. 5.23f; that it has occurred only in the flat regions of Omega mold and almost no width increment was observed on the slant regions of the Omega mold. This behavior can be understood when the reaction forces on the die and punch are plotted which are shown in Fig. 5.25. They reveal that due to the geometry/design of the mold, there is almost no contact force in the slant regions and the load-path almost entirely passes through the flat sections of the ply.

In fact, this information can be helpful to the designers as it helps in identifying the regions

Iterative design update process

To show the iterative design process for this demonstrative example, the geometry of punch was locally modified by translating nodes normal to their plane using Hypermesh. The modifications are done based on the observed deformations from simulation. The original design discussed in Sec. 5.6.1 is referred to here as Design-1 shown in Fig. 5.26a. Recalling that the consolidation was found to occur only in zones A , C and E . On the other hand, the zones B and D experienced no noticeable consolidation.

Based on this observation, the nodes on faces B and D of the punch were moved inwards (orthogonally to the face) by a small distance of about 0.6 mm (which was later found to be significant) to arrive at Design-2. The simulation outcomes of Design-2 are shown in Fig. 5.26b. It showed the exact opposite trend as that of Design-1. The consolidation was found to occur in zones B and D whereas the other zones A , C and E showed no significant consolidation. Again, this is evident from the reaction force distribution plot (Fig. 5.26b). This provided an indication that for a good consolidation in all zones, the geometry of the punch on should be somewhere in between Design-1 and Design-2.

Also, an important observation can be made with respect to the fibre orientation distribution for various designs as shown in the top view in Fig. 5.27. By tuning the mold design and thereby tweaking the amount of squeeze flow occurring in each zone, it was possible to maintain the straightness of the fibres in case of Design-3 which was not possible with the initial designs Design-1 and Design-2 where the local barreling type effect made significant changes in the fibre orientations.

Based on this conclusion, Design-3 was created this time by moving the nodes on faces B and D by 0.1 mm with respect to the Design-1. The simulation on this design, showed a relatively uniform distribution of consolidation on all faces as shown in Fig. 5.26c. Note that this iterative design procedure is discussed mainly for demonstrative purpose and it can be continued further to tune the design. This simulation consisted of 6,000 deformable shell elements (total nodes ≈ 30,000 including the nodes on die and punch). In terms of performance, it took ≈ 10 minutes with 4 threads with shared memory processing on a local computer: Intel(R) Xeon(R) CPU E3-1545M v5 @ 2.90GHz (x86_64) with the natural time-step of the model.

Comparison with the experiment

In order to validate the results quantitatively; an experimental test was conducted on the Omega mold discussed in the previous section. The material used was a UD ply stack [0] 4 of dimensions 420 mm × 60 mm and thickness of 0.3025 mm each resulting in a total Based on the comparison, several observations can be made. With the exception of face E , all other faces have comparable results. Considering all faces, the average of absolute errors in the widths calculated from simulation and obtained from the experiments was found to In summary, this validation shows that the simulation using the unified approach are able to simulate the squeeze flow behavior. It is also able to predict the final deformed shape with a reasonable accuracy with the current modeling approach.

Conclusion

In this work, a unified approach to the modeling of both forming and consolidation phases has been proposed. It not only simplifies the overall workflow from a user point of view but also overcomes several major difficulties experienced by the two-step approach.

This approach was used in a full-scale numerical simulation of QSP ® on an industrial part developed by CETIM. It was found that the results obtained from the unified approach are in agreement with the experimental observations when the fibre orientations are various locations were compared. Additionally, with this unified approach, one has access to the additional ply level information such as transverse normal stress σ zz which can serve as a tool to determine the locations with low level of consolidation and also help in the design of the molds and/or the 2D stack itself.

This approach was then used on a semi-industrial model consisting of a UD ply first to qualitatively demonstrate the impact of squeeze flow during consolidation on the final deformed shape as well as fibre orientations. The final deformed shape predicted by the simulation was then quantitatively validated against the experiments and it was observed that the simulation had correctly predicted the final deformed shape with a reasonable accuracy.

Chapter 6

Conclusion and perspectives

Conclusion

The 'Quilted Stratum Process' (QSP ® ) shares several features with the standard composite thermostamping process but at the same time there are some additional challenges both in terms of actual physical process as well as its numerical modeling. These challenges mainly come because of the usage of thermoplastic prepreg patches (both UD and woven), inability to use a blank holder, possible long distance sliding of prepreg patches during forming and the transverse squeeze flow occurring during consolidation especially for UD patches.

Thus, the overall goal of this thesis was to develop the capabilities of the numerical process modeling of forming and consolidation phases of the QSP ® process in the commercial code of Altair Radioss ™ . This was achieved by developing and implementing new functionalities in the code which have their origin in the capturing of different physical phenomena.

To provide a global overview, Chapter 2 focused on the numerical modeling of forming phase of QSP ® whereas Chapter 3 and Chapter 4 focused on the developments related to the numerical simulation of consolidation. Finally, the intention of Chapter 5 was to combine the developments from all the previous chapters and develop a unified simulation approach for both forming and consolidation phases of QSP ® .

In Chapter 2, a novel semi-empirical contact mechanism to model the phenomenon of the interply adhesion is developed as the usual adhesion models cannot be used for QSP ® due to its specific nature. This model requires only the values of adhesive strength and viscosity of melt polymer as inputs which makes it highly suitable for industrial applications where an extensive characterization campaign is often a challenge in terms of cost and time. Using this model of interply adhesion, it is shown that the prediction of final positions of the plies has improved significantly. Also, the fibre orientation predictions are shown to be in good agreement with the experimental observations made on a full-scale industrial QSP ® part designed by CETIM.

In Chapter 3, the existing formulation of a full integration shell element of Altair Radioss ™ (QBAT) is extended in order to have a functionality of transverse normal (pinching) stress. Thus, the resulting shell element (QBATP) can be in a three dimensional state of stress unlike the classical shells which are built with a plane-stress assumption. This makes it a suitable candidate to model the consolidation process where the compressive force applied in the through-thickness direction is the main driving mechanism for deformations. After validating the element itself and its interaction with the contact interface, a further improvement is made in terms of its computational performance. By adding a selective mass scaling that specifically affects the higher frequencies without impacting the lower ones; the critical time-step increased significantly. Overall, it is shown that this QBATP element serves as a first step towards a unified full-scale simulation of both forming and consolidation phases in a generic composite forming process as well as a more specific one such as QSP ® .

Chapter 4 focused on the development of an elasto-plastic constitutive model based on the global plasticity approach in order to model the behavior of the melt thermoplastic polymer during the consolidation process. It is designed to be compatible with the pinching shell element (QBATP). Additionally a methodology for the characterization of this constitutive model is proposed which uses the flat sections of the actual mold thereby avoiding the need to have a dedicated test setup. This makes it highly suitable for full-scale industrial applications. The characterized constitutive model is then used to simulate a relatively simple yet essential stack configuration of [0 • /90 • ] 2 . The results obtained from the simulation showed a good correlation with the experimental observations for this prepreg stack after undergoing a consolidation process.

Finally, Chapter 5 is dedicated to developing a unified approach to the modeling of both forming and consolidation phases of QSP ® by combining the developments done in all the previous chapters. This unified approach not only simplifies the overall workflow from a user point of view but also overcomes several major difficulties experienced by the twostep approach. This approach was used in a full-scale numerical simulation of QSP ® on an industrial part developed by CETIM. It is found that the results obtained from the unified approach are in agreement with the experimental observations when the fibre orientations are various locations are compared. Additionally due to the use of pinching shell (QBATP) in this approach; one has access to some additional ply level information such as transverse normal stress σ zz . This can serve as a tool to determine the locations with low level of consolidation and also help in the design of the molds and/or the 2D stack itself. A semi-industrial model consisting of a UD ply is simulated and it is found that the final deformed shape predicted by the simulation was in good agreement with the experimental observations within a reasonable tolerance.

Perspectives and future work

This work is a step towards a full-scale numerical process simulation of all the phases of the QSP ® process. This work was focused on the phases of forming and consolidation; but there are several other sub-steps such as overmolding and demolding in QSP ® , that need to be considered. Following points could serve as interesting next steps in continuation of the work done in this thesis.

• In terms of physics, effects such as the spring-back are known to impact the final shape of the component. Thus, it would be interesting to investigate if and how this phenomenon could be captured using the QBATP element.

• In terms of material modeling, development of a visco-plastic constitutive model (compatible with QBATP shell) is an essential ingredient for correctly simulating the behavior of the melt polymer. Ideally, it would be interesting to have a single shell element simulating both the fibre and the polymer instead of using two overlapping components as was done in this work. This would require adding anisotropy and the specific anisotropic constitutive models to the QBATP shell.

• In order to have even better representation of the actual process, a coupled thermomechanical simulation approach could be considered. This along with advanced constitutive models that can take into account the viscosity of the polymer, changes of material parameters with temperature and the crystallinity of the thermoplastics could improve the numerical process modeling of QSP ® even further.

• It was found during the experimental campaign that an excessive squeeze flow can literally separate the fibres resulting in some local through-thickness gaps. It could be an interesting next step to build a failure criteria based on the plastic strain in the constitutive model developed in this work to try and predict these local failures.

• In terms of capabilities of the QBATP shell element itself, it can be extended into a 7parameter shell which would make it possible to use all the existing three dimensional constitutive models directly without any modifications even though this will require a significant amount of work. Abstract: 'Quilted Stratum Process' (QSP ® ) uses strategically placed thermoplastic prepreg patches in its stack to form a part. Thus, even though it shares several features with the standard composite thermostamping process; there are some additional challenges that arise due to the usage of prepreg patches, along with the inability to use a blank holder, possible long distance sliding of prepreg patches and the transverse squeeze flow occurring during consolidation especially for UD patches.

A novel semi-empirical contact mechanism to model the interply adhesion is developed in the commercial code of Altair Radioss TM which has improved the prediction of patch positions and fibre orientations significantly.

A new full-integration shell element with the capability of transverse normal stress and selective mass scaling is developed by modifying the shell element in Altair Radioss TM for the numerical modeling of consolidation where the classic shells cannot be used due to their plane-stress assumption.

An elasto-plastic constitutive model for the new shell element is developed in order to model the behavior of the melt thermoplastic polymer during consolidation. A practical method for its characterization is proposed.

Finally, using these developments, a unified approach is developed for the numerical modeling of forming and consolidation. It is validated on an industrial part and results are found to be in good agreement with the experimental observations.

6 . 2

 62 Perspectives and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures 1.1 Applications of composites in various industries (a) Composite parts of a C-5 cargo aircraft [2] (b) Composite parts being used in the automotive industry [3] (c) Composite blades of the Ingenuity helicopter which is part of NASA's Mars mission [4] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 (a) An example of a typical polymer composite [5] (b) Various fibre architectures of polymer composites [6] . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 A schematic representation of various stages in a thermoplastic composite forming process (reproduced from [10]) (a) Heating stage (b) Forming stage (c) Consolidation stage (d) Demolding stage . . . . . . . . . . . . . . . . . . 1.4 (a) Flat two dimensional stack of prepregs before QSP ® (b) Final formed part 1.5 (a) Desired final component (b) Ply-cutting with classical approach (c) Plycutting with QSP ® . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.

4 . 2

 42 . . . . . . . . . . . . . . . . . . . . . . 1.10 A schematic representation of a typical evolution of temperature during QSP ® (Note: the image is not drawn to scale) . . . . . . . . . . . . . . . . . . . . . . 1.11 A schematic representation of the prepreg stack . . . . . . . . . . . . . . . . . 1.12 Classical mechanisms of interest . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 (a) Demonstration of intra-ply slip for [0] 10 laminate [22] (b) Demonstration of predominantly inter-ply slip for a [0, 90, ±45, 90, 0, ±45, 0 2 ] laminate [22] . 1.14 Mechanisms at interface which are specifically more important for QSP ® (a) Sliding of patches (b) Rotation of patches (c) Out-of-plane tension . . . . . . ix 1.15 Schematic representation of a transverse squeeze flow for a single UD ply undergoing consolidation with free lateral edges (a) side view (b) top view . . 1.16 (a) Deformed shape of a [0] 4 laminate after consolidation where initial ply size is shown by green rectangle (the initial fibre orientation is shown by yellow colored lines) (b) Deformation of a UD prepreg patch of a component after QSP ® (the initial fibre orientation is shown by yellow colored lines). Width measurements at various locations demonstrate the deformation of this patch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.17 (a) Schematic representation of a non-uniform consolidation/squeeze (b) Experimental evidence of non-uniform consolidation/squeeze [16] . . . . . . . . 1.18 Schematic representations of (a) In-plane transverse resin bleed (b) In-plane longitudinal resin bleed (c) Transverse (out of plane) resin percolation . . . . 2.1 A typical stress-strain curve obtained from a probe-tack test (reproduced here from the Fig. 6 in [5]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Schematic of the floating roller peel test (reproduced here from Fig. 4 in [15]) 2.3 (a) Schematic representation of prepreg stack (b) Modeling strategy for prepreg stack for QSP ® in Altair Radioss ™ (c) Zoomed in view of the actual component with 9 plies modeled in Altair Radioss ™ . . . . . . . . . . . . . . . . . . . . . . 2.4 Schematic representation of a case with (a) no penetration (b) penetration (denoted by P ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 (a) Standard contact mechanism (b) Modified contact mechanism with transverse adhesion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 (a) Flat plates geometry and setup (b) Imposed normal interply distance . . . 2.7 Interply contact force and contact energy obtained from simulation (a negative force signifies penetration resistive force whereas a positive force signifies adhesive force) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 A microscopic image of a test specimen showing the irregularity of the thickness of the resin-rich zone (reproduced from Fig.4.12 in [34]) . . . . . . . . 2.9 (a) Global view of the setup designed to measure the tensile mechanical behavior of molten PA-66 (b) Experimental evidence of the fingering and cavitation effects occurring respectively at low and high debonding velocities: cross-section views, (c) Top views extracted from [7] . . . . . . . . . . . . . . 2.10 Quantitative results showing the raw force and displacement recorded for several force rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 (a) Curve fitting for adhesive strength experimental values (b) Histogram of local debonding rates from simulation . . . . . . . . . . . . . . . . . . . . . . 2.12 Overview of the explicit solver flowchart (green colored boxes show the implemented and/or modified framework for the developments in this Chapter) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Overview of the contact algorithm flowchart (green colored boxes show the implemented and/or modified framework for the developments in this Chapter) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 2.14 Calculation of the projection of the secondary node on the main segment . . . 2.15 (a) Industrial model geometry and simulation setup (b) Prepreg stack configuration (Ply 1 to 9: from bottom to top in the stack) . . . . . . . . . . . . . 2.16 (a) Without interply adhesion: resistance free delamination (b) With adhesion: no delamination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Comparison of the final position of plies (a) Without interply adhesion model (b) With interply adhesion model . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Vertical displacement of a patch corner (a) Experiment: 6 mm (b) Simulation without interply adhesion: 31.7 mm (c) Simulation with interply adhesion: 6.6 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.19 Comparison of the angle made by top ply with horizontal (a) Experiment: 60 • (b) Simulation without interply adhesion: 27.9 • (c) Simulation with interply adhesion: 64.4 • . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 Position measurement locations on Ply 9 for sensitivity study . . . . . . . . . . 2.21 Fibre orientation comparison information (a) Local coordinate systems (b) Measurement locations (c) Simulation with adhesion: Fibre direction 1 (d) Simulation with adhesion: Fibre direction 2 . . . . . . . . . . . . . . . . . . . 3.1 Representation of various shell theories: (a) Kirchhoff-Love theory (b) Reissner-Mindlin theory (c) 6 or 7 parameter shells . . . . . . . . . . . . . . . . . . . . 3.2 (a) QBAT shell geometry and nodal fiber representation (b) QBAT element representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 QBATP element(a) Geometry (b) Reference element (c) Nodal orthogonal coordinate systems (d) Covariant basis . . . . . . . . . . . . . . . . . . . . . . 3.4 QBATP element pinching DoF sign convention . . . . . . . . . . . . . . . . . . 3.5 QBATP element: Transverse shear evaluation (reproduced from [25]) . . . . . 3.6 F p ext : Development of pinching pressure type of loading . . . . . . . . . . . . 3.7 Characteristic length (a) Solid or Solid-shell element (b) Shell element . . . . 3.8 Overview of the explicit solver flowchart (green colored boxes show the implemented and/or modified framework in this work) . . . . . . . . . . . . . . 3.9 Overview of the F int and F p int calculation (green colored boxes show the implemented framework in this work) . . . . . . . . . . . . . . . . . . . . . . . 3.10 Framework of constitutive models inside Altair Radioss ™ (green colored boxes show the implemented framework in this work) . . . . . . . . . . . . . . . . . 3.11 Test 3.1: Geometry, setup and loading (side view (left) and top view (right)) . 3.12 Test 3.1: Global comparison of Y displacement (mm) field for QBATP (bottom) and BRICK8 element (top) at various time instants (a) t = 0.000 s (b) t = 0.005 s (c) t = 0.010 s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Test 3.1: Comparison of evolution of nodal displacements for Node N3 for QBATP element and BRICK8 element (a) X displacement (mm) (b) Y displacement (mm) (c) Z displacement (mm) . . . . . . . . . . . . . . . . . . . . 3.14 Test 3.1: Comparison of QBATP and BRICK8 element: Evolution of thickness . xi 3.15 Test 3.1: Comparison of QBATP and BRICK8 element (a) Evolution of the through-thickness normal (Pinching) stress (MPa) (b) Evolution of the timestep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Test 3.2: Geometry, setup and loading (side view) . . . . . . . . . . . . . . . . 3.17 Test 3.2: Displacement field comparison between BRICK8 element and QBATP element (a) Displacement X (mm) for BRICK8 (b) Displacement X (mm) for QBATP (c) Displacement Y (mm) for BRICK8 (d) Displacement Y (mm) for QBATP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Test 3.2: Global comparison of QBATP and BRICK8 element (a) Nodal displacements X/Y/Z (mm) at N3 (b) Normal stresses σ xx , σ yy and σ zz (MPa) . . 3.19 Test 3.3: Geometry, setup and loading (side view (left) and top view (right)) . 3.20 Test 3.3: Comparison of deformed shape (a) BRICK8 element (b) QBATP element (The black wireframe shows the initial configuration) . . . . . . . . . 3.21 Test 3.3: Comparison between BRICK8 element and QBATP element nodal displacements at N2(a) Displacement X (mm) (b) Pinching displacement Z (mm) for QBATP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.22 Test 3.3: Comparison of QBATP and BRICK8 element (a) In-plane shear stress (σ xy ) and pinching stress (σ zz ) in MPa (b) Evolution of thickness (mm) . . . . 3.23 Test 3.4: Geometry, setup and loading (top view) . . . . . . . . . . . . . . . . 3.24 Test 3.4: Comparison of solution using QBATP elements, using BRICK8 elements and from Soulat-2006 (a) Relative position in X direction (in mm) of point C with respect to point A (b) Evolution of thickness (in mm) . . . . . . . 3.25 Test 3.5: Zoomed view of the mesh with (a) BRICK8 elements (b) QBATP elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.26 Test 3.5: Comparison of displacement field X (mm) (a) BRICK8 element (b) QBATP element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.27 Test 3.6: Geometry, setup, mesh and loading . . . . . . . . . . . . . . . . . . . 3.28 Test 3.6: Comparison of magnitude of displacement (a) QBAT element (b) QBATP element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.29 Test 3.6: (a) Evolution of Z displacement of the central node for QBAT and QBATP element simulations (b) Applied pressure versus observed deflection of the central node for QBATP element case and referecnes from the literature [23, 29] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.30 Test 3.7: Geometry, setup, mesh and loading (a) 3D view (b) Side view of the deformed shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.31 Test 3.7: Comparison of the Z displacement field for (a) QBAT element (b) QBATP element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.32 Test 3.7: Evolution of displacements at the tip for QBAT and QBATP element (a) Displacement Z (mm) (b) Displacement X (mm) . . . . . . . . . . . . . . 3.33 Test 3.8: Geometry, setup, mesh and loading for (a) BRICK8 and (b) QBAT elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.34 Test 3.8: Deformed shape (10X) and Y displacement (mm) field comparison for (a) BRICK8 elements (b) QBATP elements . . . . . . . . . . . . . . . . . . xii 3.35 Test 3.8: Comparison of the Z displacement (in mm) field for (a) BRICK8 element (b) QBATP element . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.36 Test 3.8: σ xx comparison for (a) BRICK8 elements (b) QBATP elements . . . . 3.37 Free body diagram (FBD) for an interface under equilibrium (a) Representation of the interface (b) FBD of secondary node S 1 (c) FBD of the main segment MM' (d) FBD of the secondary node S 2 . . . . . . . . . . . . . . . . . 3.38 Free body diagram for the main segment under unequal contact forces from top and bottom (a) Representation of the interface (b) Forces on the main segment (c) Contact forces in case of a standard shell (d) Additional contact forces in case of a pinching shell . . . . . . . . . . . . . . . . . . . . . . . . . 3.39 Test 3.9: Geometry, setup, mesh and loading for BRICK8, QBAT and QBATP elements (from left to right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.40 Test 3.9: Evolution of normal contact force (N) . . . . . . . . . . . . . . . . . 3.41 Test 3.9: Comparison of displacement X (mm) for (a) BRICK8 (b) QBAT (c) QBATP elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.42 Test 3.9: Comparison of three elements: BRICK8, QBAT and QBATP with respect to (a) evolution of the gap between the top and bottom plates (b) evolution of thickness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.43 Test 3.9: Representation of the inconsistent behavior and resolution strategy (a) Gap does not get updated with thickness (b) Gap gets updated with thickness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.44 Test 3.9: Comparison of displacement X (mm) for (a) BRICK8 (b) QBATP without dynamic gap update (c) QBATP with dynamic gap update . . . . . . . 3.45 Test 3.9: Comparison of (a) distance reduction and (b) thickness evolution for BRICK8, QBATP without gap dynamic gap update and QBATP with dynamic gap update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.46 Forming and consolidation of a semi-industrial model: (a) Geometry, setup and mesh (b) Actual mold and setup installed at Centrale Nantes (c) Imposed punch displacement (mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.47 Forming and consolidation of a semi-industrial model: Comparison of displacement Z (mm) for (a) QBAT simulation (b) QBATP simulation . . . . . . . 3.48 Forming and consolidation of a semi-industrial model: Demonstration of the overall distribution of the normal reaction force for: (a) QBAT simulation (b) QBATP simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.49 Forming and consolidation of a semi-industrial model: Comparison of though thickness normal stress (local σ zz ) for: (a) QBAT simulation (b) QBATP simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.50 Forming and consolidation of a semi-industrial model: Evolution of the model time-step for QBAT and QBATP simulations . . . . . . . . . . . . . . . . . . . 3.51 QBATP element: Transverse shear evaluation (reproduced from [25]) . . . . . 4.1 Schematic representation of a transverse squeeze flow for a single UD ply undergoing consolidation (a) side view (b) top view . . . . . . . . . . . . . . xiii Impact of transverse squeeze flow: (a) In-plane deformed shape of a [0 • ] 5 laminate after consolidation (initial geometry is shown by green rectangle and the initial fibre orientation is shown by yellow colored lines) (b) Width measurements for a UD prepreg patch of a component manufactured using QSP ® (the fibre orientation is shown by yellow colored lines). (c) Changes in the fibre orientations for a stack [0 • /90 • ] 6 (reproduced Fig. 5 from [1]) . . . . 4.3 Schematic representation of the modeling strategy in this work . . . . . . . . 4.4 Test 4.1: Geometry, setup and loading . . . . . . . . . . . . . . . . . . . . . . 4.5 Test 4.1: Comparison of QBATP and BRICK8 element for various values of hardening modulus (HM: 0, 500 and 1000 MPa) (a) Evolution of σ xx (MPa) (b) Evolution of equivalent plastic strain eq p . . . . . . . . . . . . . . . . . . . 4.6 Test 4.1: Comparison of QBATP and BRICK8 element for various values of hardening modulus (HM: 0, 500 and 1000 MPa) (a) Displacement Y of N3 (mm) (b) Pinching displacement Z of N3 (mm) . . . . . . . . . . . . . . . . . 4.7 Test 4.2: Geometry, setup and loading . . . . . . . . . . . . . . . . . . . . . . 4.8 Test 4.2: Comparison of QBATP and BRICK8 element for various values of hardening modulus (HM: 200, 500 and 1000 MPa) (a) Evolution of Displacement X/Y of N3 (mm) (b) Pinching displacement Z of N3 (mm) . . . . . . . . . 4.9 Test 4.2: Comparison of QBATP and BRICK8 element for various values of hardening modulus (HM: 200, 500 and 1000 MPa) (a) Pinching stress σ zz vs pinching strain zz (b) Evolution of equivalent plastic strain eq p . . . . . . . . 4.10 Test 4.3: Geometry, setup and loading . . . . . . . . . . . . . . . . . . . . . . 4.11 Test 4.3: Comparison of X displacement (mm) field for QBATP (bottom) and BRICK8 element (top) at various timesteps (a) t = 0.002 s (b) t = 0.004 s (c) t = 0.006 s (d) t = 0.008 s (e) t = 0.010 s . . . . . . . . . . . . . . . . . . 4.12 Test 4.3: Comparison of QBATP and BRICK8 element (a) Evolution of Displacements X/Y/Z of N3 (mm) (b) Evolution of normal stresses σ xx , σ yy and σ zz (MPa) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Test 4.3: Comparison of equivalent plastic strain eq p for QBATP and BRICK8 element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Test 4.4: (a) Geometry, setup and loading (b) Comparison of QBAT and QBATP element : Force vs displacement behavior for plates with different thicknesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Test 4.4: Comparison of σ yy for QBAT and QBATP element (h = 0.1 mm) . . . 4.16 Test 4.5: (a) Geometry, mesh (32×32), setup and loading (b) Deformed shape (ν=0.495) (Undeformed shape indicated by the wireframe) . . . . . . 4.17 Test 4.5: (a) Geometry, setup and loading (b) Deformed shape . . . . . . . . . 4.18 Test 4.6: (a) Geometry, mesh (64 x 36 QBATP elements) and loading (b) Comparison of reaction force Y (N) vs displacement Y (mm) for (1) QBATP element (2) Reference solution from [38] (3) ABAQUS reference solution obtained from [38] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.19 Test 4.6: Final deformed shape (a) Isometric view (b) Visualization of shear (c) Equivalent plastic strain field ( eq p ) . . . . . . . . . . . . . . . . . . . . . . xiv 4.20 Test 4.6: Final deformed shape (reproduced from Fig.25 of [38]) (a) ABAQUS reference solution obtained from [38] (b) Reference solution from [38] . . . 4.21 Test 4.7: (a) Geometry, mesh and loading (b) Comparison of load (KN) vs displacement (cm) for node B (1) QBATP element (2) Reference solution from [38] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Test 4.7: Deformed shape comparison (a) QBATP (u y =150 cm) (b) QBATP (u y =200 cm) (c) QBATP (u y =250 cm) (d) Dujc-2012 (u y =150 cm) (e) Dujc-2012 (u y =200 cm) (f) Dujc-2012 (u y =250 cm) . . . . . . . . . . . . . . . . . 4.23 Representation of the characterization procedure for the elasto-plastic constitutive model proposed in this work . . . . . . . . . . . . . . . . . . . . . . . 4.24 Representation of the characterization workflow . . . . . . . . . . . . . . . . . 4.25 Experimental setup for characterization of material parameters . . . . . . . . 4.26 Flat UD [0 • ] 5 case: (a) Initial geometry and mesh (b) Displacement Z (mm) at 5% consolidation (c) Displacement Z (mm) at 10% consolidation (d) Displacement Z (mm) at 15% consolidation (e) Displacement Z (mm) at 20% consolidation (f) Displacement Z (mm) at 25.1% consolidation . . . . . . . . 4.27 Flat UD [0 • ] 5 case: (a) Deformed shape from experiment (b) Deformed shape from simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Flat [0 • /90 • ] 2 case: (a) Measurements on the final deformed shape from experiment (b) Measurements on the final deformed shape from simulation (c) Field of the magnitude of displacement from simulation . . . . . . . . . . . . 5.1 An overview of the literature up to 2018 in the field of simulation of composite forming process (Reproduced Fig.3 from [1]) . . . . . . . . . . . . . . . . 5.2 Flowchart of a two-step procedure for thermoplastic composite forming simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 An example of a typical double curvature industrial model (Courtery of Dallara)166 5.4 Challenges in ensuring a good mesh quality of HEX layer (a) Ideal deformed shape after forming phase (b) Real deformed shape after forming phase (c) Example of small gaps and bad mesh quality for an industrial model . . . . . 5.5 (a) Penetration (b) Intersection . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Micrograph of curved and flat region of a component at the end of the forming stage (reproduced from [3]) . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Schematic representation of two strategies for simulation of composite forming process: a two-step approach versus a unified solution . . . . . . . . . . . 5.8 Ingredients of the unified approach . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Schematic representation of the modeling strategy for unified approach . . . . 5.10 (a) Industrial model geometry and simulation setup (b) Prepreg stack configuration (Ply 1 to 9: from bottom to top in the stack) . . . . . . . . . . . . . 5.11 Comparison of final deformed shape (a) Experiment (b) Simulation with a unified approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Before consolidation phase begins (a) Plot of through-thickness normal stress σ zz (MPa) for all plies (b) Magnitude of contact force (N) on punch . . . . . xv 5.13 During the consolidation phase (a) Plot of transverse normal stress σ zz (MPa) for all plies (b) Magnitude of contact force (N) on punch . . . . . . . . . . . 5.14 Final stage after the consolidation phase (a) Plot of transverse normal stress σ zz (MPa) for all plies (b) Magnitude of contact force (N) on punch . . . . . 5.15 Through-thickness pinching stress σ zz (MPa) during simulation (a) Locations on Ply-1 where the stress is tracked (b) Evolution of σ zz (MPa) during simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Fibre orientation comparison information (a) Local coordinate systems (b) Measurement locations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.17 Comparison of fibre orientations on Face-1 . . . . . . . . . . . . . . . . . . . . 5.18 Comparison of fibre orientations on Face-2 . . . . . . . . . . . . . . . . . . . . 5.19 Comparison of fibre orientations on Face-3 . . . . . . . . . . . . . . . . . . . . 5.20 Comparison of fibre orientations on Face-4 . . . . . . . . . . . . . . . . . . . . 5.21 Comparison of fibre orientations on Face-5 and Face-6 . . . . . . . . . . . . . 5.22 Unified simulation of a UD ply (a) Geometry, setup and mesh (b) Imposed punch displacement (mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.23 Unified simulation of a UD ply: Deformed shape and displacement Z (mm) at various instants (a) Initial position (b)-(c) Forming phase (d) Beginning of consolidation phase (e)-(f) Consolidation phase . . . . . . . . . . . . . . . . . 5.24 Unified simulation of a UD ply: Fibre orientations within the ply at various instants (a) Beginning of consolidation phase (b) During consolidation phase (c) End of consolidation phase . . . . . . . . . . . . . . . . . . . . . . . . . . 5.25 Unified simulation of a UD ply: Reaction force distribution on punch and die during consolidation phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.26 Unified simulation of a UD ply: Displacement Z and Reaction force distribution for: (a) Design-1 (b) Design-2 (c) Design-3 . . . . . . . . . . . . . . . . . 5.27 Unified simulation of a UD ply: Top view showing fibre orientations: (a) Design-1 (b) Design-2 (c) Design-3 . . . . . . . . . . . . . . . . . . . . . . . . 5.28 Unified simulation of a UD ply on Omega mold: Final deformed shape from the experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.29 Unified simulation of a UD ply on Omega mold: Width comparison between experiments and simulation (measurements are in mm) (initial width is shown with the white colored lines) . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi List of Tables 2.1 Specifications of the force, displacement and temperature cells for the DMA Metravib 150+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Experimental protocol steps with the corresponding imposed force and temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Experimentally obtained values of adhesive strength of PA-66 at high temperature for different debonding rates . . . . . . . . . . . . . . . . . . . . . . 2.4 List of parameters for each constitutive model used in simulation where ρ is the density, E 1 and E 2 are the Young's modulii of fibres for warp and weft directions, G 0 is the initial shear modulus, G T is the tangent shear modulus at locking angle α lock . For resin, ν is the Poisson's ratio, µ 1 , µ 2 , α 1 and α 2 are material parameters for the Ogden material law. . . . . . . . . . . . . . . . . . 2.5 List of interply adhesion parameters used in simulation . . . . . . . . . . . . . 2.6 Sensitivity of adhesive strength . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Fibre orientations (in degrees) measured from experiments . . . . . . . . . . 2.8 Fibre orientations (in degrees) measured from simulation . . . . . . . . . . . 2.9 Summary of fibre orientation comparison for each face . . . . . . . . . . . . . 3.1 A broad overview of different shell theories . . . . . . . . . . . . . . . . . . . 3.2 QBATP element: Shape functions and their derivatives . . . . . . . . . . . . . 3.3 QBATP element: Gauss points and weights . . . . . . . . . . . . . . . . . . . . 3.4 Test 3.1: Comparison of time-steps (in seconds) for QBATP and BRICK8 element 3.5 Test 3.4: Comparison of the maximal relative position in X direction (in mm) of point C with respect to point A . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Test 3.5: Performance comparison of BRICK8 vs QBATP element . . . . . . . . 3.7 Test 3.6: Comparison of the deflection (in mm) of the central node . . . . . . 3.8 Test 3.7: Comparison of the Z displacement (in mm) of the tip . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Applications of composites in various industries (a) Composite parts of a C-5 cargo aircraft [2] (b) Composite parts being used in the automotive industry [3] (c) Composite blades of the Ingenuity helicopter which is part of NASA's Mars mission[4]

Figure 1 . 2 :

 12 Figure 1.2: (a) An example of a typical polymer composite [5] (b) Various fibre architectures of polymer composites[6]

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: A schematic representation of various stages in a thermoplastic composite forming process (reproduced from [10]) (a) Heating stage (b) Forming stage (c) Consolidation stage (d) Demolding stage

Figure 1 . 5 :

 15 Figure 1.5: (a) Desired final component (b) Ply-cutting with classical approach (c) Ply-cutting with QSP ®

  One such example is the SHERLOC QSP ® project. Structural HEalth monitoring, manufacturing and Repair for Life management Of Composite fuselage (SHERLOC) is a part of CLEAN SKY 2 initiative. Within this project, CETIM was selected to manufacture demonstrators of aircraft window frame (Fig. 1.6a) using QSP ® . Besides this, QSP ® has demonstrated its effectiveness with several other components such as low cost bumper beam (Fig. 1.6b) and suspension arm (Fig. 1.6c).

Figure 1 . 6 : 1 . 2 . 2

 16122 Figure 1.6: Some examples of components made using QSP ® (a) Aircraft window frame from SHERLOC QSP ® project[12] (b) A low cost composite bumper beam (c) A composite suspension arm[11]

Figure 1 . 7 :

 17 Figure 1.7: A schematic representation of sub-steps in QSP ® production line

Figure 1 . 8 :

 18 Figure 1.8: Reproduced Fig. 3 from [14] demonstrating the overview of the literature in the field of simulation of composite forming process upto the year 2018

Figure 1 . 9 :

 19 Figure1.9: Experimental evidence of existence of resin-rich layer (a) A resin-rich layer at the interface of a [0/90] laminate[15] (b) Resin-rich layers at each interface of a specimen made of IMA-M21[16]

Figure 1 .

 1 Figure 1.10: A schematic representation of a typical evolution of temperature during QSP ® (Note: the image is not drawn to scale)

Figure 1 .Figure 1 . 12 :

 1112 Figure 1.11: A schematic representation of the prepreg stack

Figure 1 .Figure 1 . 14 :

 1114 Figure 1.13: (a) Demonstration of intra-ply slip for [0] 10 laminate [22] (b) Demonstration of predominantly inter-ply slip for a [0, 90, ±45, 90, 0, ±45, 0 2 ] laminate[22]

Figure 1 . 15 :

 115 Figure 1.15: Schematic representation of a transverse squeeze flow for a single UD ply undergoing consolidation with free lateral edges (a) side view (b) top view

Figure 1 . 18 :

 118 Figure 1.18: Schematic representations of (a) In-plane transverse resin bleed (b) In-plane longitudinal resin bleed (c) Transverse (out of plane) resin percolation

Figure 2 . 1 :

 21 Figure 2.1: A typical stress-strain curve obtained from a probe-tack test (reproduced here from the Fig. 6 in [5])

Figure 2 . 2 :

 22 Figure 2.2: Schematic of the floating roller peel test (reproduced here from Fig. 4 in [15])

Figure 2 . 3 :

 23 Figure 2.3: (a) Schematic representation of prepreg stack (b) Modeling strategy for prepreg stack for QSP ® in Altair Radioss ™ (c) Zoomed in view of the actual component with 9 plies modeled in Altair Radioss ™

Figure 2 . 4 :Figure 2 . 5 :

 2425 Figure 2.4: Schematic representation of a case with (a) no penetration (b) penetration (denoted by P )

Figure 2 . 6 :

 26 Figure 2.6: (a) Flat plates geometry and setup (b) Imposed normal interply distance

Figure 2 . 7 :

 27 Figure 2.7: Interply contact force and contact energy obtained from simulation (a negative force signifies penetration resistive force whereas a positive force signifies adhesive force)

Figure 2 . 8 :

 28 Figure 2.8: A microscopic image of a test specimen showing the irregularity of the thickness of the resin-rich zone (reproduced from Fig. 4.12 in [34])

Table 2 . 1 :

 21 Specifications of the force, displacement and temperature cells for the DMA Metravib 150+

Figure 2 . 9 :

 29 Figure 2.9: (a) Global view of the setup designed to measure the tensile mechanical behavior of molten PA-66 (b) Experimental evidence of the fingering and cavitation effects occurring respectively at low and high debonding velocities: cross-section views, (c) Top views extracted from[7]

Figure 2 . 10 :

 210 Figure 2.10: Quantitative results showing the raw force and displacement recorded for several force rates

- 15 -Figure 2 . 11 :

 15211 Figure 2.11: (a) Curve fitting for adhesive strength experimental values (b) Histogram of local debonding rates from simulation

Figure 2 . 13 :

 213 Figure 2.13: Overview of the contact algorithm flowchart (green colored boxes show the implemented and/or modified framework for the developments in this Chapter)

Figure 2 . 14 :

 214 Figure 2.14: Calculation of the projection of the secondary node on the main segment

Figure 2 .

 2 Figure 2.15: (a) Industrial model geometry and simulation setup (b) Prepreg stack configuration (Ply 1 to 9: from bottom to top in the stack)

Figure 2 .Figure 2 . 17 :

 2217 Figure 2.16: (a) Without interply adhesion: resistance free delamination (b) With adhesion: no delamination

Figure 2 . 18 :

 218 Figure 2.18: Vertical displacement of a patch corner (a) Experiment: 6 mm (b) Simulation without interply adhesion: 31.7 mm (c) Simulation with interply adhesion: 6.6 mm

Figure 2 . 19 :

 219 Figure 2.19: Comparison of the angle made by top ply with horizontal (a) Experiment: 60 • (b) Simulation without interply adhesion: 27.9 • (c) Simulation with interply adhesion: 64.4 •

Figure 2 . 20 :

 220 Figure 2.20: Position measurement locations on Ply 9 for sensitivity study

  mer was burned off to unstack the laminate. These measurements were verified at several locations using X-ray tomography scanning technique.The measurement locations are shown in Fig.2.21b along with the local coordinate systems used for each face (Fig.2.21a). For a quadrilateral face, the measurement locations lie on the diagonals and are 25 mm away from the intersection of diagonals. For a triangular face, the measurement location lies on the angle bisector of the top angle and is located at one third distance from the base of the triangle. In terms of nomenclature for Fig.2.21b, a generic location is globally identified with two numbers I.J, where I indicates the face and J indicates the location on Ith face. For a woven fabric prepreg, the orientation of both the fibre families was recorded. This in total resulted in 186 data points. These values are compared with the fibre orientations at the same locations for the simulation with added interply adhesion contact model (σ adh = 0.0062 MPa).

Figure 2 . 21 :

 221 Figure 2.21: Fibre orientation comparison information (a) Local coordinate systems (b) Measurement locations (c) Simulation with adhesion: Fibre direction 1 (d) Simulation with adhesion: Fibre direction 2

3. 7 . 1 A

 71 system under equilibrium . . . . . . . . . . . . . . . . . . . . . . 101 3.7.2 General state of contact and pinching . . . . . . . . . . . . . . . . . 101 3.7.3 Test 3.9: Validation of pinching and contacts . . . . . . . . . . . . . 103 3.8 Towards a unified simulation of forming and consolidation . . . . . . . 108 3.8.1 Forming and consolidation of a semi-industrial model . . . . . . . 108 3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 3.10 Appendix A : Geometrical calculations for QBATP . . . . . . . . . . . . . 116 3.11 Appendix B : Construction of strain-rate matrices for QBATP . . . . . . 118

Figure 3 . 1 :

 31 Figure 3.1: Representation of various shell theories: (a) Kirchhoff-Love theory (b) Reissner-Mindlin theory (c) 6 or 7 parameter shells

Figure 3 . 2 :

 32 Figure 3.2: (a) QBAT shell geometry and nodal fiber representation (b) QBAT element representation

3 . 3 4 i=1N i u i ( 3 . 5 )Figure 3 . 3 :

 3343533 Figure 3.3: QBATP element (a) Geometry (b) Reference element (c) Nodal orthogonal coordinate systems (d) Covariant basis

Figure 3 . 4 :

 34 Figure 3.4: QBATP element pinching DoF sign convention

  [B m ] (2) Curvatures [B b ] (3) Pinching [B p ] and (4) Transverse shear [B ts ]. The details for the individual strain rate matrices can be found in Appendix B of this chapter.

Figure 3 . 5 :

 35 Figure 3.5: QBATP element: Transverse shear evaluation (reproduced from [25])

Figure 3 . 6 :

 36 Figure 3.6: F p ext : Development of pinching pressure type of loading

Figure 3 . 7 :

 37 Figure 3.7: Characteristic length (a) Solid or Solid-shell element (b) Shell element

Fig. 3 .

 3 Fig.3.9 shows a flowchart for internal force calculation for both classic F int and pinching F p int .Loop over elements

Figure 3 . 9 :

 39 Figure 3.9: Overview of the F int and F p int calculation (green colored boxes show the implemented framework in this work)

Figure 3 . 10 :

 310 Figure 3.10: Framework of constitutive models inside Altair Radioss ™ (green colored boxes show the implemented framework in this work)

Figure 3 . 11 :

 311 Figure 3.11: Test 3.1: Geometry, setup and loading (side view (left) and top view (right))

Figure 3 . 12 :

 312 Figure 3.12: Test 3.1: Global comparison of Y displacement (mm) field for QBATP (bottom) and BRICK8 element (top) at various time instants (a) t = 0.000 s (b) t = 0.005 s (c) t = 0.010 s

Figure 3 . 13 :

 313 Figure 3.13: Test 3.1: Comparison of evolution of nodal displacements for Node N3 for QBATP element and BRICK8 element (a) X displacement (mm) (b) Y displacement (mm) (c) Z displacement (mm)

Figure 3 . 14 :Figure 3 . 15 :

 314315 Figure 3.14: Test 3.1: Comparison of QBATP and BRICK8 element: Evolution of thickness

Figure 3 . 16 :

 316 Figure 3.16: Test 3.2: Geometry, setup and loading (side view)

Figure 3 . 17 :

 317 Figure 3.17: Test 3.2: Displacement field comparison between BRICK8 element and QBATP element (a) Displacement X (mm) for BRICK8 (b) Displacement X (mm) for QBATP (c) Displacement Y (mm) for BRICK8 (d) Displacement Y (mm) for QBATP

Figure 3 . 18 :Figure 3 . 19 :Figure 3 . 20 :Figure 3 . 21 :Figure 3 . 22 :

 318319320321322 Figure 3.18: Test 3.2: Global comparison of QBATP and BRICK8 element (a) Nodal displacements X/Y/Z (mm) at N3 (b) Normal stresses σxx, σyy and σzz (MPa)

Figure 3 . 23 :

 323 Figure 3.23: Test 3.4: Geometry, setup and loading (top view)

Figure 3 . 24 :Table 3 . 5 : 3 . 6 . 5 Test 3 . 5 :

 3243536535 Figure 3.24: Test 3.4: Comparison of solution using QBATP elements, using BRICK8 elements and from Soulat-2006 (a) Relative position in X direction (in mm) of point C with respect to point A (b) Evolution of thickness (in mm)

Figure 3 . 25 :

 325 Figure 3.25: Test 3.5: Zoomed view of the mesh with (a) BRICK8 elements (b) QBATP elements

Figure 3 . 26 :

 326 Figure 3.26: Test 3.5: Comparison of displacement field X (mm) (a) BRICK8 element (b) QBATP element

Figure 3 . 27 :Figure 3 . 28 :Fig. 3 .

 3273283 Figure 3.27: Test 3.6: Geometry, setup, mesh and loading

Figure 3 . 30 :

 330 Figure 3.30: Test 3.7: Geometry, setup, mesh and loading (a) 3D view (b) Side view of the deformed shape

Figure 3 . 31 :

 331 Figure 3.31: Test 3.7: Comparison of the Z displacement field for (a) QBAT element (b) QBATP element

Figure 3 . 32 :

 332 Figure 3.32: Test 3.7: Evolution of displacements at the tip for QBAT and QBATP element (a) Displacement Z (mm) (b) Displacement X (mm)

Figure 3 . 33 :

 333 Figure 3.33: Test 3.8: Geometry, setup, mesh and loading for (a) BRICK8 and (b) QBAT elements

Figure 3 . 34 :

 334 Figure 3.34: Test 3.8: Deformed shape (10X) and Y displacement (mm) field comparison for (a) BRICK8 elements (b) QBATP elements

Figure 3 . 35 :

 335 Figure 3.35: Test 3.8: Comparison of the Z displacement (in mm) field for (a) BRICK8 element (b) QBATP element

Figure 3 . 37 : 2 3. 7 . 2

 337272 Figure 3.37: Free body diagram (FBD) for an interface under equilibrium (a) Representation of the interface (b) FBD of secondary node S 1 (c) FBD of the main segment MM' (d) FBD of the secondary node S 2

Figure 3 . 38 :

 338 Figure 3.38: Free body diagram for the main segment under unequal contact forces from top and bottom (a) Representation of the interface (b) Forces on the main segment (c) Contact forces in case of a standard shell (d) Additional contact forces in case of a pinching shellNote that the convention followed here is that |F upper | is oriented downwards (towards main segment) and |F lower | is oriented upwards (towards main segment). This is also depicted in Fig.3.38b. Now when the modes are summed up with the appropriate direction, two forces are obtained: first the classic contact force |F contact | (Fig.3.38c) and additionally a second one for pinching mode |F p contact | (Fig.3.38d). This is the underlying concept for adding the contribution of pinching in the interfaces. This idea when generalized, results in the following global assembly of contact forces for a given main segment:

Figure 3 . 40 :

 340 Figure 3.40: Test 3.9: Evolution of normal contact force (N)

Figure 3 . 41 :

 341 Figure 3.41: Test 3.9: Comparison of displacement X (mm) for (a) BRICK8 (b) QBAT (c) QBATP elements

Figure 3 . 42 :Figure 3 . 43 :

 342343 Figure 3.42: Test 3.9: Comparison of three elements: BRICK8, QBAT and QBATP with respect to (a) evolution of the gap between the top and bottom plates (b) evolution of thickness

Figure 3 . 44 :

 344 Figure 3.44: Test 3.9: Comparison of displacement X (mm) for (a) BRICK8 (b) QBATP without dynamic gap update (c) QBATP with dynamic gap update

Figure 3 . 45 :

 345 Figure 3.45: Test 3.9: Comparison of (a) distance reduction and (b) thickness evolution for BRICK8, QBATP without gap dynamic gap update and QBATP with dynamic gap update

Figure 3 . 46 :

 346 Figure 3.46: Forming and consolidation of a semi-industrial model: (a) Geometry, setup and mesh (b) Actual mold and setup installed at Centrale Nantes (c) Imposed punch displacement (mm)

Figure 3 . 47 :

 347 Figure 3.47: Forming and consolidation of a semi-industrial model: Comparison of displacement Z (mm) for (a) QBAT simulation (b) QBATP simulation

Figure 3 . 48 :Figure 3 . 49 :

 348349 Figure 3.48: Forming and consolidation of a semi-industrial model: Demonstration of the overall distribution of the normal reaction force for: (a) QBAT simulation (b) QBATP simulation

Figure 3 . 50 :

 350 Figure 3.50: Forming and consolidation of a semi-industrial model: Evolution of the model time-step for QBAT and QBATP simulations

[bc] = bc 11 bc 12 bc 21 bc 22 = b 22 n 3 .

 223 11 Appendix B : Construction of strain-rate matrices for QBATP This section shows the details about the construction of various [B] matrices.

Figure 4 . 1 :

 41 Figure 4.1: Schematic representation of a transverse squeeze flow for a single UD ply undergoing consolidation (a) side view (b) top view

Figure 4 . 2 :

 42 Figure 4.2: Impact of transverse squeeze flow: (a) In-plane deformed shape of a [0 • ] 5 laminate after consolidation (initial geometry is shown by green rectangle and the initial fibre orientation is shown by yellow colored lines) (b) Width measurements for a UD prepreg patch of a component manufactured using QSP ® (the fibre orientation is shown by yellow colored lines). (c) Changes in the fibre orientations for a stack [0 • /90 • ] 6 (reproduced Fig. 5 from [1])

Figure 4 . 3 :

 43 Figure 4.3: Schematic representation of the modeling strategy in this work

2 xy( 4 . 5 ) 2 xy

 2452 This makes the definition of the Von-Mises stress as, Standard Shell : Ñ 2 x + Ñ 2 y -Ñx Ñy + 3 Ñ Pinching Shell : Ñ 2 x + Ñ 2 y + Ñ 2 z -Ñx Ñy -Ñy Ñz -Ñz Ñx + 3 Ñ

Figure 4 . 4 :

 44 Figure 4.4: Test 4.1: Geometry, setup and loading

Figure 4 . 5 :Figure 4 . 6 : 4 . 5 . 2 Test 4 . 2 :

 454645242 Figure 4.5: Test 4.1: Comparison of QBATP and BRICK8 element for various values of hardening modulus (HM: 0, 500 and 1000 MPa) (a) Evolution of σxx (MPa) (b) Evolution of equivalent plastic strain eq p

Figure 4 . 7 :

 47 Figure 4.7: Test 4.2: Geometry, setup and loading

Figure 4 . 8 :Figure 4 . 9 : 4 . 5 . 3 Test 4 . 3 :

 484945343 Figure 4.8: Test 4.2: Comparison of QBATP and BRICK8 element for various values of hardening modulus (HM: 200, 500 and 1000 MPa) (a) Evolution of Displacement X/Y of N3 (mm) (b) Pinching displacement Z of N3 (mm)

Figure 4 . 10 :

 410 Figure 4.10: Test 4.3: Geometry, setup and loading

Figure 4 . 11 : 4 . 5 . 4 Test 4 . 4 :

 41145444 Figure 4.11: Test 4.3: Comparison of X displacement (mm) field for QBATP (bottom) and BRICK8 element (top) at various timesteps (a) t = 0.002 s (b) t = 0.004 s (c) t = 0.006 s (d) t = 0.008 s (e) t = 0.010 s

Figure 4 . 14 :

 414 Figure 4.14: Test 4.4: (a) Geometry, setup and loading (b) Comparison of QBAT and QBATP element : Force vs displacement behavior for plates with different thicknesses

Figure 4 . 16 :Figure 4 . 17 :Fig. 4 .

 4164174 Figure 4.16: Test 4.5: (a) Geometry, mesh (32×32), setup and loading (b) Deformed shape (ν=0.495) (Undeformed shape indicated by the wireframe)

  In this test, a hollow cylinder of radius r = 285 mm, height H = 850 mm and thickness h = 5 mm undergoes shear loading (Fig.4.18a). The bottom edge was fixed in all degrees of freedom and an imposed velocity was applied on the top edge in Y direction until a displacement of 100 mm is reached. All other DoFs on this edge were constrained. The material properties were set to: E = 210000 MPa, ν = 0.3, σ y = 240 MPa and hardening modulus (HM) was set to 0 making it a case of perfectly plasticity.

Figure 4 . 18 :

 418 Figure 4.18: Test 4.6: (a) Geometry, mesh (64 x 36 QBATP elements) and loading (b) Comparison of reaction force Y (N) vs displacement Y (mm) for (1) QBATP element (2) Reference solution from[38] (3) ABAQUS reference solution obtained from[38]

Figure 4 . 19 :Figure 4 . 20 : 4 . 5 . 7 Test 4 . 7 :

 41942045747 Figure 4.19: Test 4.6: Final deformed shape (a) Isometric view (b) Visualization of shear (c) Equivalent plastic strain field ( eq p )

Figure 4 . 21 :

 421 Figure 4.21: Test 4.7: (a) Geometry, mesh and loading (b) Comparison of load (KN) vs displacement (cm) for node B (1) QBATP element (2) Reference solution from[38]

Figure 4 . 22 : 4 . 6

 42246 Figure 4.22: Test 4.7: Deformed shape comparison (a) QBATP (uy=150 cm) (b) QBATP (uy=200 cm) (c) QBATP (uy=250 cm) (d) Dujc-2012 (uy=150 cm) (e) Dujc-2012 (uy=200 cm) (f) Dujc-2012 (uy=250 cm)

Figure 4 . 23 :

 423 Figure 4.23: Representation of the characterization procedure for the elasto-plastic constitutive model proposed in this work

Figure 4 . 24 :

 424 Figure 4.24: Representation of the characterization workflow

Figure 4 . 25 :

 425 Figure 4.25: Experimental setup for characterization of material parameters

2 :

 2 List of parameters for each constitutive model used in simulation where ρ is the density, E 1 and E 2 are the Young's modulii of fibres for warp and weft directions, G 0 is the initial shear modulus, G T is the tangent shear modulus at locking angle α lock . E is the Young's modulus, ν is the Poisson's ratio, σy is the yield stress and HM is the hardening modulus. idation of 25.1%. This being a comparatively small model with less contacts and a much shorter distance to travel for the punch, the simulation took about 3-4 minutes to run on a local personal computer. The displacement field along direction Z at various stages of the consolidation is shown in Fig. 4.26b to Fig. 4.26f.

Figure 4 . 26 :

 426 Figure 4.26: Flat UD [0 • ] 5 case: (a) Initial geometry and mesh (b) Displacement Z (mm) at 5% consolidation (c) Displacement Z (mm) at 10% consolidation (d) Displacement Z (mm) at 15% consolidation (e) Displacement Z (mm) at 20% consolidation (f) Displacement Z (mm) at 25.1% consolidation

Figure 4 .

 4 Figure 4.27: Flat UD [0 • ] 5 case: (a) Deformed shape from experiment (b) Deformed shape from simulation

Figure 4 . 28 :

 428 Figure 4.28: Flat [0 • /90 • ] 2 case: (a) Measurements on the final deformed shape from experiment (b) Measurements on the final deformed shape from simulation (c) Field of the magnitude of displacement from simulation

Figure 5 . 1 :

 51 Figure 5.1: An overview of the literature up to 2018 in the field of simulation of composite forming process (Reproduced Fig.3 from [1])

Figure 5 . 2 :

 52 Figure 5.2: Flowchart of a two-step procedure for thermoplastic composite forming simulation

Figure 5 . 3 :

 53 Figure 5.3: An example of a typical double curvature industrial model (Courtery of Dallara)

Figure 5 . 4 :

 54 Figure 5.4: Challenges in ensuring a good mesh quality of HEX layer (a) Ideal deformed shape after forming phase (b) Real deformed shape after forming phase (c) Example of small gaps and bad mesh quality for an industrial model

Figure 5 . 5 :

 55 Figure 5.5: (a) Penetration (b) Intersection

Fig. 5

 5 .6 shows the result given in[3] of the micrograph at different locations of the same component (Fig.5.6b) at the end of the forming stage. It was observed that the porosities in the curved part of the geometry (Fig.5.6a) were removed during the forming phase itself whereas they still exist in the flat vertical region of the geometry (Fig.5.6c).

Figure 5 . 6 :

 56 Figure 5.6: Micrograph of curved and flat region of a component at the end of the forming stage (reproduced from [3])

Figure 5 . 7 :

 57 Figure 5.7: Schematic representation of two strategies for simulation of composite forming process: a two-step approach versus a unified solution

Figure 5 .

 5 Figure 5.10: (a) Industrial model geometry and simulation setup (b) Prepreg stack configuration (Ply 1 to 9: from bottom to top in the stack)

Figure 5 . 11 :

 511 Figure 5.11: Comparison of final deformed shape (a) Experiment (b) Simulation with a unified approach

Figure 5 . 12 :

 512 Figure 5.12: Before consolidation phase begins (a) Plot of through-thickness normal stress σzz (MPa) for all plies (b) Magnitude of contact force (N) on punch

Figure 5 . 13 :Figure 5 . 14 :

 513514 Figure 5.13: During the consolidation phase (a) Plot of transverse normal stress σzz (MPa) for all plies (b) Magnitude of contact force (N) on punch

Figure 5 . 15 :

 515 Figure 5.15: Through-thickness pinching stress σzz (MPa) during simulation (a) Locations on Ply-1 where the stress is tracked (b) Evolution of σzz (MPa) during simulation

Figure 5 . 16 :

 516 Figure 5.16: Fibre orientation comparison information (a) Local coordinate systems (b) Measurement locations

Figure 5 . 17 : 1 Figure 5 . 18 : 2 Figure 5 . 19 : 3 Figure 5 . 20 : 4 Figure 5 . 21 :

 5171518251935204521 Figure 5.17: Comparison of fibre orientations on Face-1

Figure 5 . 22 :

 522 Figure 5.22: Unified simulation of a UD ply (a) Geometry, setup and mesh (b) Imposed punch displacement (mm)

Fig. 5 .

 5 Fig. 5.23a to Fig. 5.23f show the deformed shape of the UD ply and the displacement in the Z direction (along the width of the UD ply) at 6 instants during the unified simulation. Each of the instants corresponds to each marker 1 to 6 shown in Fig. 5.22b. The Z displacement has been chosen here to demonstrate the lateral spreading of UD plies due to squeezing.

Figure 5 . 23 :

 523 Figure 5.23: Unified simulation of a UD ply: Deformed shape and displacement Z (mm) at various instants (a) Initial position (b)-(c) Forming phase (d) Beginning of consolidation phase (e)-(f) Consolidation phase

Figure 5 . 24 :

 524 Figure 5.24: Unified simulation of a UD ply: Fibre orientations within the ply at various instants (a) Beginning of consolidation phase (b) During consolidation phase (c) End of consolidation phase

Figure 5 . 25 :

 525 Figure 5.25: Unified simulation of a UD ply: Reaction force distribution on punch and die during consolidation phase

Figure 5 . 26 :

 526 Figure 5.26: Unified simulation of a UD ply: Displacement Z and Reaction force distribution for: (a) Design-1 (b) Design-2 (c) Design-3

Figure 5 . 27 :

 527 Figure 5.27: Unified simulation of a UD ply: Top view showing fibre orientations: (a) Design-1 (b) Design-2 (c) Design-3

Figure 5 . 28 :

 528 Figure 5.28: Unified simulation of a UD ply on Omega mold: Final deformed shape from the experiments

Figure 5 . 29 :

 529 Figure 5.29: Unified simulation of a UD ply on Omega mold: Width comparison between experiments and simulation (measurements are in mm) (initial width is shown with the white colored lines)

Titre:

  Modélisation numérique unifiée du formage et de la consolidation des composites thermoplastiques avec des patchs préimprégnés Mots clés : simulation numérique, adhésion inter-pli, formage, consolidation, coque, préimprégné Résumé : Le « Quilted Stratum Process » (QSP ® ) est un procédé qui utilise des patchs de préimprégné thermoplastique placés stratégiquement pour former une pièce composite. L'utilisation de patchs préimprégnés introduit des particularités par rapport au formage classique, comme l'impossibilité d'utiliser un cadre serre-flanc, la possibilité de grands déplacements des patchs préimprégnés pendant le formage et la présence d'écoulement de compression pendant la phase de consolidation, en particulier pour les patchs unidirectionnels. Un nouveau modele semi-empirique de contact est développé dans le code commercial d'Altair Radioss TM pour prendre en compte l'adhesion inter-pli. Il a permis une amélioration de la prédiction des positions des patchs et des orientations des fibres dans les plis. Un nouvel élément de coque à intégration complète est développé en modifiant l'élément de coque existant dans Altair Radioss TM base sur l'hypothèse de contraintes planes. Le nouvel élément permet de calculer les contraintes normales dans l'épaisseur. Il est enrichi par une méthode d'ajout artificiel et sélectif de masse pour réduire les temps de calcul. Un modèle constitutif élasto-plastique de la nouvelle coque est développé afin de modéliser le comportement du polymère thermoplastique pendant la phase de la consolidation. Une méthode d'identification du modele est proposée. Enfin, en utilisant l'ensemble de ces développements, une approche unifiée est proposée pour la modélisation numérique du formage et de la consolidation. Elle est validée sur une pièce industrielle. Title : Unified Numerical Modeling of Forming and Consolidation of Thermoplastic Composites with Prepreg Patches Keywords : process simulation, interply adhesion, forming, consolidation, shell element, prepregs

3 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

  List of parameters for each constitutive model used in simulation where ρ is the density, E 1 and E 2 are the Young's modulii of fibres for warp and weft directions, G 0 is the initial shear modulus, G T is the tangent shear modulus at locking angle α lock . E is the Young's modulus, ν is the Poisson's ratio, σ y is the yield stress and HM is the hardening modulus. . . . . . . . . . . . . . . . . 5.1 List of parameters for each constitutive model used in simulation where ρ is the density, E 1 and E 2 are the Young's modulii of fibres for warp and weft directions, G 0 is the initial shear modulus, G T is the tangent shear modulus at locking angle α lock . For the polymer, E is the Young's modulus, ν is the Poisson's ratio, σ y is the yield stress and HM is the hardening modulus. For interply adhesion, η is the viscosity of the melt polymer and σ adh is the adhesive strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Benefits of QSP ® . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Sub-steps in QSP ® . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.3 Challenges in the numerical simulation of QSP ® . . . . . . . . . . . 7 1.2.4 Mechanisms of interest . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.

	4.2
	4.1 Material model parameters for melt PA-66 obtained from characterization . . 152
	xvii

Contents 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Quilted Stratum Process: QSP ® . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1

Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

	2.1.1 Interply Adhesion and QSP ® . . . . . . . . . . . . . . . . . . . . . . 23
	2.1.2 Polymer Tack . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
	2

.2 Objectives and content of the study . . . . . . . . . . . . . . . . . . . . . 25 2.3 Numerical modeling of interply adhesion . . . . . . . . . . . . . . . . . . 26

	2.3.1 Contact mechanism for interply adhesion in normal direction . . . 27
	2.3.2 Mechanism verification with a simple test case . . . . . . . . . . . 29
	2.3.3 Interply contact mechanism in tangential direction . . . . . . . . . 31
	2

.4 Experimental determination of the interply adhesion strength . . . . . 32

	2.4.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
	2.4.2 Experimental methods and setup . . . . . . . . . . . . . . . . . . . 33
	2.4.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . 34
	2.4.4 Choice of adhesive strength value for simulation . . . . . . . . . . 35
	2

.5 Implementation aspects of interply adhesion in Altair Radioss ™ . . . . . 37

	2.5.1 Global overview of contact algorithm . . . . . . . . . . . . . . . . . 37
	2.5.2 Calculation of the secondary node projection and penetration . . . 40
	2.5.3 Distribution of penalty forces on main segment . . . . . . . . . . . 41
	2

.6 Full scale numerical simulation of QSP ® using Altair Radioss ™ . . . . . 41

	2.6.1 Model description and setup . . . . . . . . . . . . . . . . . . . . . . 41
	2.6.2 Comparison of the final position of a prepreg patch . . . . . . . . . 44
	2.6.3 Comparison of fibre orientations for individual plies . . . . . . . . 45
	2

.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Table 2 .

 2

	m

• C/min to T 2: Experimental protocol steps with the corresponding imposed force and temperature

Table 2 .

 2 3: Experimentally obtained values of adhesive strength of PA-66 at high temperature for different debonding rates

	Debonding rate (1/s) Adhesive Strength (MPa)
	2.00E-04	7.72E-04
	6.00E-04	1.40E-03
	5.00E-03	2.31E-03
	1.50E-02	4.28E-03
	4.00E-02	5.75E-03

Table 2 .

 2 4: List of parameters for each constitutive model used in simulation where ρ is the density, E 1 and E 2 are the Young's modulii of fibres for warp and weft directions, G 0 is the initial shear modulus, G T is the tangent shear modulus at locking angle α lock . For resin, ν is the Poisson's ratio, µ 1 , µ 2 , α 1 and α 2 are material parameters for the Ogden material law.

	Parameter	Symbol	Value
	Viscosity of melt polymer (280 • C)	η	275 Pa.s
	Adhesive strength	σ adh	0.0062 MPa
	Viscous scaling factor	S	1.0

Table 2 . 5

 25

: List of interply adhesion parameters used in simulation

Table 2 . 6

 26

: Sensitivity of adhesive strength

Table 2 .

 2 .21)

	Loc	P1D1	P1D2	P2D1	P2D2	P3D1	P3D2	P4D1	P4D2	P5D1	P6D1	P7D1	P8D1	P9D1
	1.2	-1.0	93.7	-3.0	91.5	0.5	92.5	94.0	3.3	98.0	15.7	15.7 106.7 147.0
	1.3	3.7	93.3	-1.5	92.5	1.5	91.5	92.3	1.0	94.7	4.0	4.0 101.3 147.0
	1.4	6.3	92.0	2.0	94.0	0.0	96.5	99.3	-1.7	94.7	5.3	5.3	98.3 140.3
	2.1	3.0	92.0	7.5	93.5	17.5	95.5	92.3	16.3	88.0	2.7
	2.2	1.7	89.3	15.5	88.5	14.0	86.0	80.3	15.7	81.7	1.3
	2.3	0.3	91.0	-10.0	92.5	-10.0	93.5	92.3	-11.7	85.3	1.7
	2.4	7.3	93.0	12.5	91.5	13.5	93.5	94.0	14.3	89.3	2.3
	3.1	11.7 100.7	13.5 102.5	24.0 104.5 104.7	28.3	96.3	16.3
	3.2	14.7	98.3	18.5	91.5	29.5	87.5	90.0	30.3	90.0	16.7
	3.3	14.7	97.3	18.5 100.5	24.5 102.0	96.3	32.7	92.3	17.0
	3.4	11.0 105.7	21.0 107.5	26.5 106.0 103.3	29.3 105.3	17.3
	4.1	167.7	84.0 161.0	85.0 152.0	84.5	87.0 152.0	83.3 158.3 158.3
	4.2	170.0	83.3 163.0	84.0 160.0	81.0	76.3 158.0	75.7 162.7 162.7
	4.3	166.3	78.3 164.0	57.5 154.0	59.0	63.7 152.7	65.3 155.3 155.3
	4.4	165.7	81.0 163.0	78.0 153.0	75.5	76.3 151.7	77.0 155.7 155.7
	5	128.0	72.5 130.5	66.0 135.5	66.5	69.5 136.5	75.0 135.5
	6	71.5 127.0	64.5 120.0	59.0 128.0 129.0	54.5 133.0	62.5	62.5 135.0	23.0

7: Fibre orientations (in degrees) measured from experiments

Table 2 .

 2

	Loc	P1D1	P1D2	P2D1	P2D2	P3D1	P3D2	P4D1	P4D2	P5D1	P6D1	P7D1	P8D1	P9D1
	1.2	1.9 100.0	2.6	99.4	3.1 100.0 100.0	3.3 100.1	6.6	9.1 106.2 147.5
	1.3	4.1	96.4	4.8	96.5	5.6	97.0	97.6	6.7	98.7	10.4	11.3	99.7 134.4
	1.4	3.7	95.1	2.1	95.3	1.4	95.7	95.9	1.4	95.9	-3.8	6.2	97.9 140.3
	2.1	8.8	97.3	3.9	95.8	22.0	97.5	96.5	16.9	97.0	-13.8
	2.2	-2.8	88.0	22.3	85.2	-2.0	84.3	85.0	11.7	84.9	1.2
	2.3	4.1	87.8	3.1	87.7	3.2	88.2	88.4	2.2	88.5	-1.2
	2.4	6.5	92.9	6.1	93.2	5.6	93.3	93.7	4.9	94.5	7.3
	3.1	25.5	96.1	25.4	96.0	25.8	95.8	95.8	26.8	96.3	26.7
	3.2	23.5	99.2	26.1	98.0	26.3	97.3	96.7	26.7	96.1	27.9
	3.3	27.2 102.1	26.5 103.2	27.2 103.3 104.3	27.6 104.3	27.1
	3.4	25.7 137.6	25.6 137.4	26.5 138.0 138.1	26.9 136.9	27.1
	4.1	157.5	79.8 157.7	80.2 154.6	78.8	82.4 155.4	79.7 157.6 157.0
	4.2	167.8	78.5 166.2	82.5 163.0	79.6	76.3 161.5	79.7 161.2 161.3
	4.3	157.9	50.5 158.0	50.7 157.2	49.5	49.1 158.0	48.3 158.7 159.1
	4.4	153.7	72.2 153.5	68.2 153.5	68.1	67.8 153.3	69.9 152.5 154.3
	5	136.3	69.8 135.9	70.2 135.9	70.6	70.9 135.8	70.4 135.9
	6	66.1 132.2	64.7 131.7	64.7 131.8 131.9	64.1 132.0	62.3	61.8 133.4	5.2
				Table 2.8: Fibre orientations (in degrees) measured from simulation
	phase of QSP ® .
			Face Data Points Average of absolute difference (degrees)
			1		39				3.9
			2		40				5.1
			3		40				9.8
			4		44				5.6
			5		10				3.2
			6		13				5.1

9: Summary of fibre orientation comparison for each face

2 Objectives and content of the study . . . . . . . . . . . . . . . . . . . . . 3.3 Full-integration shell element with transverse stress (QBATP)

  Numerical modeling of a consolidation process . . . . . . . . . . . 3.1.2 Solid based approach . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Solid-shell based approach . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Shell based approach . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.5 On the choice of the approach . . . . . . . . . . . . . . . . . . . . . 3.. . . . . . 3.3.1 Concept of transverse normal stress and pinching . . . . . . . . . .

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 3.

3.2 Shell element description . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Calculation of strain rates and strain . . . . . . . . . . . . . . . . . 3.3.4 Constitutive model . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.5 Stress internal force relationship . . . . . . . . . . . . . . . . . . . 3.3.6 External force and pinching pressure . . . . . . . . . . . . . . . . . 3.3.7 Dynamic system of equations . . . . . . . . . . . . . . . . . . . . . 3.4 Selective mass scaling for QBATP element . . . . . . . . . . . . . . . . . 3.4.1 Stability and time-step . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Method A: Uniform mass scaling . . . . . . . . . . . . . . . . . . . 3.4.3 Method B: Selective mass scaling based on stiffness matrix . . . . . 3.4.4 Method C: Selective mass scaling based on local mass matrix . . . 3.4.5 Method D: Selective mass scaling based on acceleration filtering . . 3.4.6 Selective mass scaling for QBATP element . . . . . . . . . . . . . . 3.4.7 Mass scaling algorithm for QBATP element . . . . . . . . . . . . . .

3.5 Implementation aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.5.1 Global framework . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Internal forces framework for QBATP element . . . . . . . . . . . . 3.5.3 Constitutive models and their framework . . . . . . . . . . . . . .

3.6 Validation Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	. .
	3.6.5 Test 3.5: Performance of QBATP element (crushing under pressure)
	3.6.6 Test 3.6: Clamped square plate . . . . . . . . . . . . . . . . . . . .
	3.6.7 Test 3.7: Large bending of cantilever . . . . . . . . . . . . . . . . .
	3.6.8 Test 3.8: Non-uniform pinching pressure . . . . . . . . . . . . . . .
	3.7 QBATP element and contacts . . . . . . . . . . . . . . . . . . . . . . . . .

3.6.1 Test 3.1: Pure pinching pressure driven loading . . . . . . . . . . . 3.6.2 Test 3.2: Coupled pinching pressure and membrane force loading . 3.6.3 Test 3.3: Coupled pinching pressure and in-plane shear loading . . 3.6.4 Test 3.4: Crushing under pressure . . . . . . . . . . . . . . . . .

Table 3 .

 3

		3	No	No	No
	Reissner-Mindlin	5	No	Yes	No
	Simo-Rifai	6	Yes	Yes	Yes (modified)
	Buchter-Ramm	7	Yes	Yes	Yes

1: A broad overview of different shell theories

Table 3 . 2

 32

: QBATP element: Shape functions and their derivatives

1 Table 3 . 3

 133

: QBATP element: Gauss points and weights

Table 3 .

 3

4: Test 3.1: Comparison of time-steps (in seconds) for QBATP and BRICK8 element

Table 3 .

 3 6: Test 3.5: Performance comparison of BRICK8 vs QBATP element

		BRICK8 QBATP
	Elements	6400	6400
	DoFs	39366	39366
	Cycles	30,974	9,201
	Time (s)	995	129
	Speedup	1	7.7

Table 3 .

 3

7: Test 3.6: Comparison of the deflection (in mm) of the central node summarized in Tab. 3.7. It demonstrates a good agreement of the QBATP solution with the reference and also with QBAT solution.

Table 3 .

 3

8: Test 3.7: Comparison of the Z displacement (in mm) of the tip

  With this algorithm of the elasto-plastic constitutive model for the pinching shell, several numerical tests are performed in order to validate the model. First, several single element tests are performed in order to assess both the in-plane (Test 4.1) and out-of-plane (pinching) behavior (Test 4.2). A case involving simultaneous pinching and membrane loading is discussed after that (Test 4.3). The results are compared with the BRICK8 element of Altair Radioss ™ . Then, several benchmark cases from the literature are compared with the results obtained with the developments of this work (Test 4.4 -4.7). Each test is further subdivided into three sub sections which specify the objective of the test, setup of the test finally providing a discussion on the obtained results.

	3. Calculate other trial stress resultants ( Ñ t xy ) and trial stress couples (M t x , M t y , M t xy ) 4.5 Validation tests
	same as in case of the elastic constitutive model.
	4. Calculate P n+1 assuming a linear equation of state (g) and performing a selective
	reduced integration to avoid volumetric locking (K is the bulk modulus).
	P n+1 = g(K, ρ n+1 , ρ 0 )	(4.9)
	5. Calculate the updated yield stress based on the equivalent plastic strain and hardening
	modulus (HM) as,
		σ y = σ 0 y + HM eq p	(4.10)
	6. Calculate the modified yield criterion calculated in eq. 4.6 which can be rewritten
	for the non-bending part using the deviatoric components of the stress resultants as
	follows,
	f = 16 M t eq	2	+ Ñ t eq	2	-σ y	(4.11)
	M 2 eq = M 2 x + M 2 y -Mx My + 3 M 2 xy Ñ 2 eq = 3 2 Ñ dev x 2 + Ñ dev y 2 + Ñ dev z	2	+ 2 Ñxy	2	(4.12) (4.13)
	7. If f ≤ 0, it implies that the material is not undergoing plastic loading signifying ˙ eq p = 0.
	{ Ñ n+1 dev } = { Ñ t dev } { Ñ n+1 } = { Ñ n+1 dev -P n+1 } Ñ n+1 xy = Ñ t xy	(4.14) (4.15) (4.16)
	{ M n+1 } = { M t }	(4.17)
	8. If f > 0, it implies that the material is undergoing plastic loading. Thus, a factor for
	radial return (R) is to be calculated with subsequent updates in the stress resultants,
	stress couples and equivalent plastic strain.
	R =	16 M t eq	σ y 2	eq + Ñ t	2	7) (4.18)
	2. Calculate a trial deviatoric stress resultant { Ñ t { Ñ n+1 } = { Ñ t dev }R -P n+1 Ñ n+1 xy = Ñ t xy R dev } assuming elastic response (G is the (4.19) (4.20) elastic shear modulus). { Ñ t dev } = { Ñ n dev } + 2G{∆ dev } (4.8) (4.21) { M n+1 } = { M t }R ∆ eq p = f /(E + HM) (4.22)

Table 4 .

 4

	1)		UD (Law-58)	Polymer (Law-91)
	ρ 7.85E-9 t/mm 3	ρ	1.8E-9 t/mm 3 ρ 1.13E-9 t/mm 3
	E 210,000 MPa	E 1	25,000 MPa	E 6.8 MPa
	ν 0.3	E 2	0.1 MPa	ν 0.495
		G 0	0 MPa	σ y 0.1 MPa
		G T	3.0 MPa	HM 0.105 MPa
		α lock 85 •

  .9 for Figure5.8: Ingredients of the unified approach a general stack consisting of 4 plies. Each ply contains two shell element components with shared coincident nodes. One of these two components represents fibres whereas the other represents the polymer. The shell for fibres is modeled using a full integration standard shell element (QBAT) of Altair Radioss ™ and uses an anisotropic material model (Law 58). On the other hand, the shell for the polymer is modeled using the QBATP element developed in Chapter 3. The elasto-plastic constitutive model (Law 91) developed in Chapter 4 is assigned to this polymer shell. The ply-ply interaction is modeled with Type 25 contact law with interply adhesion as discussed in Chapter 2 and in[11].

			Unified Approach
	Chapter 2		Chapter 3	Chapter 4
	Contacts		Elements	Materials
	Interply Adhesion		Pinching Shell (QBATP)	Elastoplastic Constitutive Model (Law 91)
	Ply 3		Ply 4
				Fibre :	: QBAT + Law58
	Contact Law	Contact Law	Resin :	: QBATP + Law91
	Ply 2	Contact Law
	Ply 1
	Figure 5.9: Schematic representation of the modeling strategy for unified approach

Table 5 .

 5 1: List of parameters for each constitutive model used in simulation where ρ is the density, E 1 and E 2 are the Young's modulii of fibres for warp and weft directions, G 0 is the initial shear modulus, G T is the tangent shear modulus at locking angle α lock . For the polymer, E is the Young's modulus, ν is the Poisson's ratio, σy is the yield stress and HM is the hardening modulus. For interply adhesion, η is the viscosity of the melt polymer and σ adh is the adhesive strength

	Woven (Law-58)		UD (Law-58)	Polymer (Law-91)	Interply Adhesion
	ρ	1.8E-9 t/mm 3	ρ	1.8E-9 t/mm 3 ρ 1.13E-9 t/mm 3	η	275 Pa.s
	E 1	10,000 MPa	E 1	25,000 MPa	E 6.8 MPa	σ adh 0.0062 MPa
	E 2	10,000 MPa	E 2	0.1 MPa	ν 0.495
	G 0	0.069 MPa	G 0	0 MPa	σ y 0.1 MPa
	G T	3.0 MPa	G T	3.0 MPa	HM 0.105 MPa
	α lock 55 •	α lock 85 •

Table 5 .

 5 .16).

	Loc	P1D1	P1D2	P2D1	P2D2	P3D1	P3D2	P4D1	P4D2	P5D1	P6D1	P7D1	P8D1	P9D1
	1.2	2.4	96.2	2.9	94.4	2.9	94.4	94.0	3.9	95.4	8.6	8.6 101.4 145.3
	1.3	6.8	95.5	3.4	95.0	2.3	95.1	95.5	2.8	94.5	5.4	10.4 100.4 141.9
	1.4	1.1	89.5	1.4	90.4	-1.3	93.1	94.4	-1.5	92.5	6.7	5.2	97.7 139.6
	2.1	4.0	89.8	-1.1	90.5	11.6	89.3	89.8	-0.9	89.2	-0.5
	2.2	-1.1	87.9	3.9	87.9	0.8	88.0	88.6	1.0	87.7	-1.7
	2.3	-2.3	86.2	-0.1	87.4	1.5	87.5	87.9	-0.9	89.4	-0.8
	2.4	4.4	90.4	-0.3	91.9	-1.0	91.6	92.2	0.3	90.2	0.9
	3.1	11.5	91.5	12.5	91.3	12.9	90.7	90.1	13.0	91.0	13.2
	3.2	15.1	90.8	16.7	89.5	16.8	89.9	89.1	17.1	89.9	16.5
	3.3	16.0	94.9	17.0	95.9	16.3	97.2	97.3	16.6 100.2	14.6
	3.4	15.2 119.1	15.9 120.3	12.7 120.1 119.2	15.2 114.8	13.4
	4.1	162.7	83.9 165.4	87.2 160.8	86.1	88.1 160.8	78.9 169.6 166.8
	4.2	169.0	88.1 167.3	83.5 165.0	82.9	81.7 164.1	85.5 164.8 166.3
	4.3	161.6	46.6 161.3	46.6 161.7	46.4	41.9 161.6	42.2 158.3 160.2
	4.4	159.3	60.2 159.2	76.1 158.2	76.3	77.7 158.7	73.3 160.3 158.9
	5	123.2	64.1 124.1	63.0 135.1	64.8	62.9 131.8	69.9 140.3
	6	66.8 126.6	67.0 126.4	67.5 128.0 127.7	67.2 128.8	62.8	63.8 132.3	27.0

2: Fibre orientations (in degrees) measured from solution with the unified approach

Table 5 .

 5 3: Summary of fibre orientation comparison for each face

	Face Data Points Avg. Error Avg. Error
			Forming	Unified
	1	39	3.9	2.7
	2	40	5.1	5.7
	3	40	9.8	7.1
	4	44	5.6	6.6
	5	10	3.2	4.6
	6	13	5.1	3.8

This chapter has been reproduced from: PD Mulye, J Hemmer, L Morançay, C Binetruy, A Leygue, S Comas-Cardona, P Pichon, and D Guillon. Numerical modeling of interply adhesion in composite forming of viscous discontinuous thermoplastic prepregs. Composites Part B: Engineering, 191:107953, 2020

3.7. QBATP element and contactsFigure 3.36: Test 3.8: σxx comparison for (a) BRICK8 elements (b) QBATP elements

Appendix A : Geometrical calculations for QBATP

This section discusses calculations of various geometrical and physical quantities needed for the shell element. The quantities involve the calculations of various co-ordinate systems, evaluation of certain quantities at specific points within the element etc. Calculation of these quantities is essential and is used in the calculation of the B matrices discussed in Appendix B.

Mapping of mid-surface into the iso-parametric space is given by,

Where the covariant basis vectors a 1 and a 2 are given by,

where, β and V are the nodal variables corresponding to the rotational velocities and translational velocities. β is the rotational velocity vector defined as,

The vectors nA1, nA2, nB1 and nB2 in eq.(3.99) are given by,

where n A1 , n A2 , n B1 and n B2 are the unit normal vectors at edge mid-points A1, A2, B1 and B2 respectively. These mid-point unit normals are calculated by taking the average of the normals at the edge end-points. As an example for the point B1 which lies at the midpoint of nodes N1 and N4, the unit normal n B1 is given by,

The quantities AS1 A1 to AS2 Ab in eq.(3.99) are given by,

(3.107)

(3.108)

where, t 1i and t 2i correspond to t 1 and t 2 at node i.

Chapter 4

An elasto-plastic constitutive model for pinching shell (QBATP) element for numerical simulation of consolidation of UD prepregs

Abstract

This chapter focuses on the development of a constitutive model for the nearly incompressible thermoplastic melt polymer which is a key ingredient in the numerical modeling of consolidation process. It becomes even more important in case of unidirectional (UD) thermoplastic prepregs where the mechanisms such as the transverse squeeze flow can impact not only the in-plane dimensions of the prepreg but also the fibre orientations within the prepreg. Thus, based on the global plasticity approach, an elasto-plastic constitutive model compatible with the pinching shell element (QBATP) is developed which can model the behavior of the melt polymer and the three dimensional state of stress experienced during consolidation process. The numerical part of the constitutive model is validated by comparison with the solid element in Altair Radioss ™ and also using several benchmark tests from the literature. For a physical comparison with respect to the consolidation process, a practical and economic characterization method is proposed. Finally, the consolidation test results on a [0 • /90 • ] 2 are compared with the simulations results.

Discussion

The displacements fields (Fig. 4.11), nodal displacements (Fig. 4.12a), normal stresses (Fig. 4.12b) and equivalent plastic strain (Fig. 4.13) showed a good correlation between the QBATP element and BRICK8 element for this test case of simultaneous transverse and membrane loading.

Consolidation process modeling of thermoplastic UD prepregs

In this section, various test cases of consolidation of UD thermoplastic prepregs are discussed. Following the modeling strategy discussed in Sec. 4.2.4, this section aims to demonstrates the capabilities of using the elasto-plastic constitutive model developed in this work with the QBATP element in order to do consolidation process modeling.

This section begins with a brief discussion about the flat [0 • ] 5 case. Even though the characterization of constitutive model was done using this prepreg stack configuration; only the maximum final width of the specimen was used for the characterization. Thus, it is interesting to observe the overall final shape and to compare that with the experimental observations. This is then followed by a validation test on a flat [0 • /90 • ] 2 sample that uses the characterized constitutive model. The simulation predictions are compared with the deformations obtained from the conducted experiment.

Flat [0] 5 Case

For the simulation, a flat UD ply of dimensions 60 mm × 60 mm and a thickness of 1.55 mm was modeled between two flat molds. The punch and die are given elastic material properties. As per the modeling strategy discussed in Sec. 4.2.4, the prepreg ply is modeled as two overlapping shell element components, one modeling the UD fibres and the other modeling the melt polymer. The initial fibre orientation is along X axis as shown in Fig. 4.26a. Overall, the model consists of 3450 4-node shell elements with a total of 10,445 nodes. The summary of the material properties used for this simulation is given in Tab. 4.2.

An imposed displacement is applied on the flat punch which induces a transverse squeeze flow in the prepreg stack whereas the die is fixed in all DoFs. The simulation is run until a final average thickness of 1.161 mm was obtained which corresponds to a total consol-Chapter 5

Numerical modeling of forming and consolidation of thermoplastic composites with prepreg patches using a shell based unified approach

Abstract

The objective of this chapter is to demonstrate how the developments done in the previous chapters (interply adhesion, QBATP shell element and elasto-plastic constitutive model for QBATP shell) when brought together can provide a unified solution for a full-scale numerical simulation of QSP ® consisting of both forming and consolidation phases in a unified manner using a shell based approach.

Introduction

A relatively recent review conducted by Bussetta et al. in 2018 has summarized the state of the art for numerical methods for composite forming simulation [1]. Referring to Fig. 5.1 (reproduced from their work); it can be seen that the main focus of research in composite forming simulation is focused on the mechanisms occurring during the forming phase. However, for QSP ® and also for a typical thermoplastic composite forming process, the forming stage is almost always followed by a consolidation stage.

Consolidation is an integral part of the overall process as it helps in the following aspects:

• Removal of porosities

• Healing of ply-ply interfaces

• Improvement of surface quality

• Controlling the thermo-mechanical transformations of the material at the onset of the cooling phase