Alexander Bernath

**Band 88**

A. Bernath

Numerical prediction of curing and PID of composite structures

# **Numerical prediction of curing and process-induced distortion of composite structures**

Alexander Bernath

**Numerical prediction of curing and process-induced distortion of composite structures**

### **Karlsruher Schriftenreihe Fahrzeugsystemtechnik Band 88**

Herausgeber

### **FAST Institut für Fahrzeugsystemtechnik** Prof. Dr. rer. nat. Frank Gauterin Prof. Dr.-Ing. Marcus Geimer Prof. Dr.-Ing. Peter Gratzfeld Prof. Dr.-Ing. Frank Henning

Das Institut für Fahrzeugsystemtechnik besteht aus den Institutsteilen Bahnsystemtechnik, Fahrzeugtechnik, Leichtbautechnologie und Mobile Arbeitsmaschinen.

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buchs.

# **Numerical prediction of curing and process-induced distortion of composite structures**

by Alexander Bernath

Karlsruher Institut für Technologie Institut für Fahrzeugsystemtechnik

Numerical prediction of curing and process-induced distortion of composite structures

Zur Erlangung des akademischen Grades eines Doktors der Ingenieurwissenschaften von der KIT-Fakultät für Maschinenbau des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Dipl.-Ing. Alexander Bernath

Tag der mündlichen Prüfung: 3. Dezember 2019 Referent: Prof. Dr.-Ing. Frank Henning Korreferent: Prof. Dr. Anthony Straatman

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2021 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 1869-6058 ISBN 978-3-7315-1063-5 DOI 10.5445/KSP/1000125453

# **Vorwort**

Die vorliegende Arbeit entstand im Rahmen meiner Tätigkeit als wissenschaftlicher Mitarbeiter am Institut für Fahrzeugsystemtechnik (FAST) des Karlsruher Instituts für Technologie (KIT).

Mein besonderer Dank gilt Prof. Dr.-Ing. Frank Henning für die Übernahme des Hauptreferats, sowie für die kontinuierliche Unterstützung und das entgegengebrachte Vertrauen. Des Weiteren bedanke ich mich bei Prof. Dr. Anthony Straatman für die Übernahme des Koreferats und für die wertvollen Anmerkungen zur vorliegenden Arbeit. Zudem gilt mein Dank Prof. Dr. Anton Möslang für die Übernahme des Prüfungsvorsitzes.

Bei meinen Kollegen vom FAST möchte ich mich für die anregenden Diskussionen und Ratschläge, die uneingeschränkte Hilfsbereitschaft sowie das angenehme Arbeitsklima bedanken. Besonders hervorherben möchte ich an dieser Stelle Dr.-Ing. Luise Kärger aufgrund ihrer wertvollen Beiträge zur vorliegenden Arbeit und dazugehörigen Publikationen. Meinem langjährigen Bürokollegen Siegfried Galkin möchte ich für die entspannte und freundschaftliche Atmosphäre und die anregenden Diskussionen danken.

Darüber hinaus gebührt mein Dank den Kollegen vom Fraunhofer Institut für Chemische Technologie (ICT) sowie den Kollegen der KIT-Institute Institut für Technische Chemie und Polymerchemie (ITCP), Institut für Angewandte Materialien (IAM) und Institut für Produktionstechnik (wbk) für die umfangreiche Unterstützung bei den verschiedensten Charakterisierungs- und Validierungsversuchen. Dies gilt in gleichem Maße für Wibke Exner und Fabian Groh vom Institut für Faserverbundleichtbau und Adaptronik des Deutschen Zentrums für Luft- und Raumfahrt (DLR) in Braunschweig. Außerdem möchte ich den von mir betreuten Abschlussarbeitern und Hilfswissenschaftlern für das Engagement und die gute Zusammenarbeit danken.

Viele der in dieser Arbeit vorgestellten Inhalte wurden im Rahmen des Schwerpunktprogramms 1712 sowie im Rahmen des Projekts "SMiLE" erarbeitet. Für die finanzielle Unterstützung möchte ich mich bei der Deutschen Forschungsgemeinschaft (DFG) sowie beim Bundesministerium für Bildung und Forschung (BMBF) bedanken.

Nicht zuletzt bedanke ich mich von ganzem Herzen bei meiner Familie - insbesondere bei meiner Frau Christiane - und meinen Freunden für die kompromisslose Unterstützung, den Rückhalt und die Geduld.

Karlsruhe, im Dezember 2019 *Alexander Bernath*

# **Abstract**

Fiber-reinforced materials offer a huge potential for lightweight design of loadbearing structures. However, high-volume production of such parts is still a challenge in terms of cost efficiency and competitiveness. This is particularly true for the automotive sector because of its enormous cost pressure. Besides high material cost of carbon fiber, cycle times are much longer compared to conventional metal based manufacturing processes. A promising strategy to shorten the cycle time is the combined use of fast-curing resins and advanced injection processes like Pressure-Controlled Resin Transfer Molding (PC-RTM). However, finding suitable process parameters is often difficult without extensive practical knowledge. Moreover, when using fast-curing resins, the shortened cycle time makes this even more difficult as the fault tolerance decreases. In addition, it is difficult to monitor and analyze the progress of individual process steps as most of them take place hidden in a closed mold. This is where numerical process simulation can help to provide a better insight and to analyze underlying mechanisms. In this study, the curing process of the resin is investigated with regard to its influence on mold filling and process-induced distortion (PID).

Curing starts to influence the process at the time at which resin and hardener are mixed and injected into the mold. From this point on, the viscosity of the resin depends not only on temperature but also on cure degree. As this interrelations affect the progression of the mold filling, proper and accurate experimental characterization is vital to gain reliable results from simulation. However, the time slot, which is relevant to mold filling and therefore is required to be captured by the measurement, depends on the reactivity of the resin. As is shown in this study, fast-curing resins challenge established measurement techniques. Automation of specimen preparation and insertion into the rheometer setup is needed to reduce the lack of viscosity data in the beginning of the measurement to a minimum. Otherwise, the accuracy of the viscosity model and the mold filling simulation will be impaired. The influence on the evolution and level of the cavity pressure is demonstrated using mold filling simulations of a complex car floor structure. The results show a substantial influence on cavity pressure throughout the course of the filling process. Hence, when simulating process variants like PC-RTM, which rely heavily on the predicted pressure, accurate viscosity characterization is vital.

After mold filling has been completed, the resin further cures and transforms from liquid to solid after reaching gelation. From this moment on, residual strains from chemical shrinkage or temperature change will induce stress. As a result, the manufactured part will show process-induced distortion after demolding, which impairs the dimensional stability and can lead to increased scrap volume. Predicting the resulting PID using process simulation can help to prevent this as suitable countermeasures can be derived without the need for expensive trials. In this study, the underlying mechanisms and their interaction with boundary conditions and design variables of the process are examined. Moreover, elastic and viscoelastic material models are derived and implemented in the form of user subroutines. Validation of the results is based on angled brackets and, as a more practical example, a complex car floor structure. Good agreement between experiment and simulation is achieved in cases in which demolding of the part did not require applying noteworthy forces. If this was unavoidable, scatter of measured values greatly increased and significant deviations between experiment and simulation were observed. Moreover, the sensitivity of the specimens to demolding forces is found to depend on laminate layup. Thus, the demolding process needs to be investigated in more detail in the future, in order to gain more knowledge on this specific process step and be able to further increase the reliability of simulation results.

# **Zusammenfassung**

Faser-Verbund-Werkstoffe ermöglichen Leichtbautragstrukturen, welche mit konventionellen Werkstoffen nicht umsetzbar sind. Deren Großserienfertigung ist allerdings auch heute noch eine große Herausforderung. Aufgrund des hohen Kostendrucks gilt dies insbesondere für den Automobilbau. Hinderlich wirken zum einen die hohen Materialkosten bei Verwendung von Kohlenstofffasern, und zum anderen die im Vergleich zu Stahlwerkstoffen lange Zykluszeit in der Bauteilfertigung. Eine Lösung für Letzteres stellen hochreaktive, schnellhärtende Matrixwerkstoffe in Verbindung mit modernen Harzinjektionsprozessen dar. Durch die schnelle Aushärtung können Zykluszeiten reduziert und die Wirtschaftlichkeit des Prozesses gesteigert werden. Die Beschleunigung der Prozessabläufe erschwert allerdings die Auslegung des Prozesses, da z.B. das mögliche Zeitfenster für die Realisierung der Formfüllung deutlich verkürzt wird. Die Findung einer optimalen Prozessführung erfordert unter diesen Bedingungen ein großes Erfahrungswissen. Dessen Generierung wird durch das Ablaufen des Prozesses innerhalb des geschlossenen Werkzeugs, und damit im Verborgenen, gehemmt. An dieser Stelle greift die Prozesssimlation an, welche einen detaillierten Einblick in die Abläufe der einzelnen Prozessschritte ermöglicht. In der vorliegenden Arbeit wird daher der Aushärteprozess des Matrixwerkstoffs mithilfe von Strömungs- und Struktursimulationen begleitet. Neben der Entwicklung von geeigneten Materialmodellen, wird auch die Ermittlung von notwendigen Kennwerten vorgestellt.

Der Einfluss der Aushärtung beginnt mit dem Vermischen der Harzkomponenten am Anguss des Injektionswerkzeugs. Von diesem Zeitpunkt an ergibt sich die für das Fließverhalten der Matrix charakteristische Viskosität aus der wirkenden Temperatur und dem erreichten Aushärtegrad. Die akkurate experimentelle Charakterisierung dieses wichtigen Zusammenhangs ist eine Grundvoraussetzung für aussagekräftige Simulationsergebnisse. Die hohe Reaktivität der verwendeten schnellhärtenden Harze fordert allerdings etablierte Methoden heraus. Durch deren Trägheit ist eine zuverlässige Messung der Viskosität während des für die Formfüllung wichtigen Zeitfensters erschwert. Um dieses Problem zu lösen, wird eine automatisierte Lösung zur Herstellung sowie Einbringung von Proben in ein Rheometer entwickelt. Hierdurch kann ein Großteil des prozessrelevanten Zeitfensters in der anschließenden Messung erfasst werden. Da diese Daten für die Parametrierung von Viskositätsmodellen herangezogen werden, ergibt sich ein großes Fehlerpotential, wenn unzureichende Informationen zum Verlauf der Viskosität während der Dauer der Formfüllung vorliegen. Dies wird in der vorliegenden Arbeit mittels Formfüllsimulationen einer komplexen Unterbodenstruktur demonstriert, welche einen signifikanten Einfluss auf den Verlauf und die Höhe des Kavitätsdrucks aufzeigen. Diese Erkenntnis ist daher insbesondere für die Simulation von druckgeregelten Prozessen, wie z.B. Pressure-Controlled Resin Transfer Molding (PC-RTM), von großer Wichtigkeit.

Im weiteren Verlauf der Aushärtung findet ein Übergang der Matrix von einem Fluid zu einem Feststoff statt. Dieser Zeitpunkt wird durch den Gelpunkt definiert. Ist dieser überschritten, ist der Werkstoff in der Lage Kräfte zu ertragen. Daher führen prozessinduzierte Dehnungen, aufgrund von chemischer Schwindung oder Temperaturänderung, ab diesem Zeitpunkt zu Eigenspannungen. Nach Entformung des Bauteils führen diese schließlich zu Gestaltsabweichungen, welche die Maßhaltigkeit negativ beeinflussen und somit zu erhöhtem Ausschuss führen können. Die Vorhersage des Bauteilverzugs mithilfe von Simulation ist daher eine wichtige Voraussetzung, um dies zukünftig, bspw. durch eine Anpassung der Werkzeuggeometrie, zu vermeiden. In der vorliegenden Arbeit werden die zugrundeliegenden Mechanismen und deren Wechselwirkung mit Prozessrandbedingungen analysiert. Für die numerische Berechnung des Verzugs werden zudem elastische sowie viskoelastische Materialmodelle parametriert und in Form von Subroutinen implementiert. Die Ergebnisse werden anhand von L-Winkel-Proben experimentell validiert. Zusätzlich wird auch hier wieder die komplexe Unterbodenstruktur als praxisrelevantes Anwendungsbeispiel herangezogen. Eine gute Übereinstimmung zwischen Versuch und Simulation wird dabei immer dann erzielt, wenn die Entformung der Bauteile ohne das Aufbringen von nennenswerten Kräften möglich ist. Ist dies nicht vermeidbar, ergeben sich zum Teil große Abweichungen, welche zudem vom Laminataufbau abhängen. Hier gilt es zukünftig anzusetzen, um sowohl die Prognosegüte als auch die Zuverlässigkeit der Prozesssimulation weiter zu erhöhen.

# **Contents**




# **Nomenclature**

### **Acronyms**



#### **Greek Symbols**



#### **Latin Symbols**





# **1 Introduction**

Lightweight design is a key element in maximizing the energy efficiency of various kinds of vehicles. High performance composites already make up 50 % of the weight of the structure of modern civil aircrafts [1, 2]. In aerospace industry, composite materials are even more important and extensively used. In the automotive sector, however, the situation is quite different, although this industry is under strong pressure to constantly reduce emissions of their vehicles [3]. Fuchs et al. [4] analyzed the importance and use of different materials in light vehicle designs with a focus on the US automotive industry. On average, more than half of the vehicle weight (63 %) stems from steel and iron parts. Although the use of polymer materials increased significantly during the last decades, the materials used were usually non-reinforced or short-fiber reinforced polymers applied to non-structural parts. The use of more advanced composites is found to be rather limited, which is mostly due to material costs. This is especially problematic as the annual production volume increases, leading to a decrease in contribution of fixed costs to final part costs. Consequently, material costs remain as the main cost driver. While this application limit is even more pushed by platform-based concepts, there is still a considerable potential for the application of advanced composite materials in automotive components, especially with respect to hybrid body-in-whites [4].

Regarding lightweight design, the allowed cost per saved mass is what separates the automotive industry from the aerospace industry. Hence, the future use of advanced composites in automotive structures depends on the availability of cost-efficient manufacturing processes as well as reduced material costs. The former can be achieved by using more recently developed processes like Pressure-Controlled Resin Transfer Molding (PC-RTM) or Wet Compression Molding (WCM), both of which allow to significantly reduce cycle times. However, as cycle times are reduced, the error tolerance of the process also decreases. Consequently, a very good and precise knowledge of what happens during processing and which factors influence the progress is vital for the design of a robust and cost-efficient manufacturing process. Besides experimentally gained knowledge, numerical prediction of the course of the process can be of great help [5]. In this thesis, this is investigated with a focus on cross-linking of the polymer matrix during processing and its consequences with regard to increase of viscosity during mold-filling but also process-induced distortion (PID) due to chemical and thermal strains. The evolution of viscosity is a key quantity in the design of injection processes since it influences the progression of the flow front during mold filling as well as the pressure needed to maintain a specific constant mass flow. The accurate prediction of the latter is important for derivation of optimal process parameters of pressure-controlled processes like PC-RTM. Process-induced distortion, on the other hand, impairs the dimensional accuracy of the manufactured component, which can increase assembly costs as well as scrap volume. As a result, the cost-efficiency of the manufacturing process is compromised.

## **1.1 Motivation**

Recently developed Liquid Composite Molding (LCM) processes like PC-RTM and WCM enable large-batch production of high-performance composite structures. However, these processes are not yet fully understood and determination of an optimal process design is therefore impeded. Moreover, LCM processes can involve complex molds and tooling concepts which makes experimental trial and error methods very expensive and, thus, infeasible. As a result, process parameters are often defined based on experience and the limited amount of published data. A promising and more general approach is the investigation of relevant mechanisms by using numerical simulations, which is pursued in this thesis. The common theme throughout the following chapters is the cross-linking

reaction of the resin which influences the process and produced parts in many ways as depicted in Figure 1.1.

**Figure 1.1:** RTM process cycle

Curing of the resin is initiated as soon as a sufficiently high temperature is reached. This typically happens shortly after injection of the resin into the mold. Hence, the first influence of cross-linking occurs during mold-filling in terms of an increase in viscosity which makes it increasingly difficult to impregnate the initially dry fiber textiles of the preform. Moreover, if a specific level of viscosity is exceeded, fiber textiles may be carried away by the resin flow and the probability of imperfect impregnation of rovings will increase dramatically. To avoid this scenario, the molding concept needs to be designed such that it allows for fast filling of the cavity with a high quality of fiber impregnation and by applying the lowest possible pressure to prevent fiber washing and foam core collapsing. Fulfillment of these requirements is hardly possible based solely on experience. Instead, mold-filling simulations are increasingly used to virtually optimize the process beforehand [6–10].

A crucial step in process simulation is proper characterization of the material's behavior. With regard to mold-filling simulation, this includes primarily characterization of reaction kinetics and viscosity evolution. Furthermore, in order to achieve the lowest possible cycle times, fast-curing resins need to be used and,

thus, characterized. However, conventional measurement techniques are often unable to follow the rapid reaction of these materials. As a consequence, more suitable techniques which help to capture as much information as possible need to be developed. Moreover, the relevant time frame in mold-filling begins with injection of the resin and ends either when the cavity has been completely filled or as soon as a critical viscosity has been exceeded. Hence, it is vital to obtain viscosity data right from the beginning of the reaction and therefore involves rethinking of specimen preparation processes in order to avoid curing of the resin prior to data acquisition.

Curing of the resin is accompanied by the evolution of process-induced strain due to chemical shrinkage and thermal contraction. While the former is a direct result of the cross-linking reaction, the latter is caused by the necessity of using a high enough curing temperature that ensures rapid curing but also a sufficiently high glass transition temperature. However, not only residual strain and corresponding stress evolve during curing as the mechanical properties of the resin undergo substantial changes too. While cross-linking reactions take place, the resin transitions from a viscous liquid to a rubbery solid. Afterwards, as the produced part is cooled down to room temperature, the resin again changes its material state to a glassy solid. The latter involves an increase in stiffness by several magnitudes. This complex interrelations challenge the process of choosing a suitable material model but also the interpretation of simulation results. Hence, a thoughtful modeling strategy is needed which enables to analyze the underlying mechanisms based on analysis results.

After demolding of the produced part, the process-induced stress leads to process-induced distortion (PID) as the part is no longer held in place by the mold. However, the subsequent assembly process requires the produced part to fulfill specific tolerances as otherwise problems will emerge due to badly fitting parts. Strategies for compensation of PID, e.g. by shimming or by forcing the part to take on the originally intended shape, are either too expensive in an automotive context or involve the application of unnecessary loads which can lead to pre-damage. Hence, the only reasonable way to fulfill predefined tolerances is by compensating for PID by altering the part's geometry such that it takes on the desired shape only after demolding. This again can be achieved by using accurate numerical simulations [5].

Two types of PID are distinguished: spring-in and warpage. Both types are a result of stress gradients in through-thickness direction of the part. Spring-in refers to an unwanted change of angle in curved areas of a part. Due to the curvature, process-induced strain in through-thickness direction induces a stress gradient which eventually leads to spring-in. Hence, spring-in is only present in curved sections. Warpage, however, can also emerge in flat areas and is caused by various sources of gradients like asymmetric layup, inhomogeneous temperature and thus cure degree distribution or varying mold boundary conditions, to name a few. With the aid of numerical simulation, the influence of these factors can be quantified and predicted, which in a next step, enables to identify possible strategies for minimizing PID.

A typical Resin Transfer Molding (RTM) process cycle as shown in Figure 1.1 begins with loading of the mold with a preform which has been stacked and assembled in a previous process step. During this draping step, the fiber textiles are deformed such that they take on the form of the final part as close as possible. This deformation results in a reorientation of the fibers but is also often accompanied by changes in local fiber volume fraction. Both preform properties can have significant impact on both mold-filling and PID. Similarly, the final state of PID is a result of process-induced stresses which to some extent cause spring-in and warpage. The remaining portion, however, causes a pre-loading of the structure and can therefore lead to premature damage and a reduced structural performance of the produced part. Hence, each individual process step depends on results of previous process steps but also provides vital information for subsequent process steps. This chaining of individual process steps needs to be replicated by numerical methods as shown in Figure 1.2, in order to maximize accuracy and robustness during the whole course of the process. Moreover, such a CAE-chain enables to optimize the whole process chain in a holistic way [11].

**Figure 1.2:** Virtual process chain combining design, process and structural simulation [11]

# **1.2 Thesis Aim**

The aim of this thesis is to provide a better understanding of the influence of cross-linking of the resin on mold-filling and dimensional accuracy of parts manufactured with liquid composite molding processes like RTM or WCM. During mold-filling the relevant material property which is affected by curing is the viscosity of the resin. Here, focus is put on the impact of curing on the evolution of injection pressure during the course of the process, which is investigated numerically. Furthermore, the experimental technique used for viscosity characterization is modified in order to obtain reliable viscosity data for model fitting also in case of fast-curing resins.

The dimensional accuracy is investigated based on spring-in of angled brackets and general process-induced distortion of a large composite structure. Focus is put on the influence of the generally viscoelastic material behavior of the resin. Furthermore, different modeling approaches are compared with regard to spring-in prediction accuracy.

## **1.3 Thesis Outline**

The thesis begins in Chapter 2 with an overview on materials and processes typically used in automotive composite manufacturing. Furthermore, details regarding the materials and processes used within this thesis are provided. In the following Chapter 3, characterization and modeling of the reaction kinetics of the resins are presented. These models are then used as a basis for describing the chemo-rheology of select resins in Chapter 4. Since one of these is of fastcuring type, a new specimen preparation process has been developed, which is detailed in that chapter as well. The subsequent Chapter 5 deals with RTM mold filling simulation which is showcased using a large composite structure. Here, the influence of measurement uncertainties during viscosity characterization on the evolution of injection pressure is investigated. Beginning with Chapter 6, focus is shifted to the chemo-thermo-mechanical behavior of the resin. Hence, in that chapter, numerical models for describing the material behavior using viscoelastic, path-dependent and linear-elastic approaches are introduced. Starting with isotropic material behavior, the models are subsequently extended for use with transversely isotropic materials. Important effects like the influence of viscoelasticity during processing and the determination of the point of vitrification are also addressed in this chapter. In Chapter 7, process-induced distortion is experimentally and numerically investigated. This is carried out using spring-in characterization based on angled brackets, a specimen geometry that is also often used in literature. Furthermore, the complex composite structure that has been investigated before in terms of RTM mold filling is used here again to showcase the phenomenon of process-induced distortion on a component level. The experimental efforts are accompanied by numerical simulations which make use of the kinetic models as well as the chemo-thermo-mechanical models introduced in Chapters 3 and 6, respectively. The thesis ends with a conclusion in Chapter 8 and an outlook in Chapter 9.

# **2 Materials and Processes**

One of the major advantages of composites is the fact that the constituents can be selected from a variety of different types of fibers and resins. This enables to tailor the resulting composite behavior to what is really needed with regard to cost-efficient manufacturing and part performance. If for example short cycle times are needed, a fast-curing resin can be used. Another example is temperature stability, which can be improved by using specialized polymer matrix systems that offer a high glass transition temperature.

In this chapter, relevant fiber and matrix materials as well as suitable manufacturing processes are introduced. Furthermore, the materials and processes used in this thesis are presented. Since only few process types are suitable for cost-effective large-batch production of high performance composite structures, the discussion will concentrate on these.

# **2.1 Fiber reinforcements**

This study focuses on large-batch production of high performance composite structures. Thus, the most relevant fiber types are carbon fiber and glass fiber. While carbon fibers have the advantage of much higher mass-specific mechanical properties, glass fibers are significantly cheaper. Since material costs can have a huge impact on actual component costs, it is vital to choose a suitable fiber type and to design the part in such a way that its performance benefits from the chosen fiber type as much as possible. This fact renders pure substitution of metal by carbon fiber inadequate and ineffective in almost every case.

Besides the fiber type, the way in which these fibers are arranged to form rovings and eventually the final fabric also influences both part performance and processing ability. Individual rovings can be arranged next to and parallel to each other to form a unidirectional (UD) ply of non-crimp fabric (NCF). Multiple layers of a unidirectional ply can be stacked and stitched together to form a multiaxial NCF. Another possibility is given by weaving of rovings oriented 90° to each other. Different weave patterns can be used to tailor important properties of the resulting fabric like draping ability or permeability to the needs of the chosen composite manufacturing process. Figure 2.1 shows fabric types commonly used in designs of composite structures.

**Figure 2.1:** Common types of weave patterns used in composite manufacturing [12]

One major difference of the individual weave pattern types is in the number of crossing points. For a plain weave fabric, this number is much higher than that of a satin type fabric. Each crossing of rovings is accompanied by an out-of-plane waviness, commonly referred to as undulation. Another source of undulation is draping of fiber textiles which can cause both in-plane and out-of-plane fiber waviness. Irrespective of its source, fiber waviness impairs mechanical properties of the manufactured parts [13, 14]. Otherwise, each crossing point facilitates the flow of resin in through-thickness direction which can be beneficial to the quality of the manufactured part.


**Table 2.1:** Properties of the fiber textiles (NCF) used in this thesis

In this thesis, high-performance composites are investigated. Hence, the used fiber textiles are of NCF type which limits fiber waviness to that caused by draping. Table 2.1 gives an overview of the textiles used. The carbon fiber content has been determined by manually separating carbon fibers from the rest of the fabric. The fiber volume content of specimens or parts manufactured with these fabrics is calculated based on this value.

## **2.2 Matrix materials**

The polymer matrix dominates the composite behavior in many ways since it is responsible for load transmission and introduction into the individual fibers. Moreover, the mechanical behavior perpendicular to the fiber direction is directly related to that of the resin. If, for instance, the composite is heated up to a temperature near or above the glass transition temperature T<sup>g</sup> of the resin system, the mechanical behavior of the composite as a whole will deteriorate significantly.

The large variety of polymeric materials is also reflected by the choices of different resin systems used for composite manufacturing. To name a few, these can differ widely in terms of viscosity, reactivity, temperature stability and cost. A good compromise, and hence often applied in prepreg and LCM processes in the past, are epoxy based resins. Today, the same materials are still widely used due to the fact that their characteristics are well analyzed and understood. This is especially valuable in the aerospace industry, in which material certification processes are often lengthy and broad information on the material behavior is needed, including aspects like flame retardation or electric conductivity. Hence, switching to another resin system, e.g. to enable faster processing, is an expensive and lengthy process. This does not apply to the automotive industry, which is why non-epoxy based resins have also come to the fore more recently. In this case, the goal is to find a suitable matrix material which still fulfills mechanical and thermal requirements, but at the same time allows for an increased cost-efficiency of the manufacturing process. The latter can be achieved by reducing cycle time and material cost. Suitable candidates for such a material are polyurethane based resins (PUR), one of which is also investigated in this thesis.



\* Cannot be published due to non-disclosure agreement

Table 2.2 lists the different resins used within this thesis. One slowly curing epoxy resin (SCER1) is used as a baseline for comparison with two fast-curing resin systems, of which one is also epoxy based (FCER3) while the other is polyurethane based (FCPUR). Since part of the work has been carried out within research projects with participation of industry, the exact names of the resins and manufacturers are confidential with the exception of the SCER1 resin. The naming scheme has been chosen similar to that used in Bernath et al. [15] in order to avoid misunderstandings when comparing results.


**Table 2.3:** Manufacturer recommended cure cycle (MRCC) of used resins

The cure cycle applied to a resin is often prescribed by the manufacturer. Table 2.3 lists the manufacturer recommended cure cycles (MRCC) of the resins used in this thesis. The advantage of fast-curing resins regarding cost-efficiency of the process becomes evident by comparing the recommended cure times with that of the slow curing resin. However, as noted in Chapter 1, a resin which cures within 2 min is also very demanding in terms of optimal process design but also accurate material characterization.

### **2.3 Processes**

Many different types of processes exist for manufacturing of composite structures, each having different advantages but also disadvantages. With regard to large-batch production, however, the number of suitable processes diminishes considerably. The main reason for this is cycle time, which needs to be as short as possible in order to maximize the economic efficiency of the process. Only if this is achieved, a composite design can compete and finally outrun alternative solutions like metal based designs. However, composite manufacturing processes are less well understood due to the fact that fiber-reinforced polymers represent a relatively young class of materials.

Another requirement of processes suitable for large-batch production is a high potential for automation of individual process steps which helps to further reduce cycle times. Possible candidates of processes which can keep-up with these requirements are RTM and WCM. Both will be discussed in more detail in the following.

## **2.3.1 Resin Transfer Molding (RTM)**

Figure 2.2 illustrates the individual process steps of a typical RTM process cycle. A large portion of the process is dedicated to cutting, stacking and draping of the dry fiber textiles, which eventually form a preform. Afterwards, the preform is transferred to the RTM mold and injection of the resin is initiated. During moldfilling, the resin starts to cross-link and is allowed to further cure afterwards. In a next step, the cured part is demolded and final trimming is applied.

**Figure 2.2:** Typical RTM process cycle [3]

The above given description of the RTM process should be interpreted as the bare minimum of what is needed and important at a first glance. Depending on part and process design, the given scheme needs to be adopted accordingly. In case of very large and complex shaped structures, draping can be difficult and a subdivision of the geometry may be required. In this case, another step is needed for assembly of the individual sub-preforms to form the final preform. While such a strategy offers a huge amount of possible subdivision solutions, finding an optimal one is difficult and should be accompanied by draping simulations. However, even after identification of a suitable sub-preforming concept, production of each individual sub-preform also depends on many factors. Depending on the geometry, fiber wrinkling or gaping between adjacent rovings may be induced by unsuitable draping concepts. These effects can be minimized using a proper draping setup which involves systematic application of tensile forces, resulting in membrane tension [16]. If multiple layers are draped simultaneously, effects like interply friction are observed [17]. Overall, optimization of draping is complex but numerical simulations can help to better understand the involved mechanisms and to finally identify an optimal strategy [18].

Conventional RTM processing is not suitable for use with fast-curing resins and, thus, is only of limited use for large-batch production. Consequently, the process has been modified in several ways to meet the requirements of the automotive industry. These include high-pressure variants like HP-RTM as well as variable gap variants like CRTM or combinations of these types. Rosenberg et al. [19] studied pressure evolution and distribution during HP-RTM and HP-CRTM processes. In case of HP-RTM, the cavity pressure increases rapidly during resin injection. Contrary to this, high cavity pressures were only found at the end of the process after resin injection stopped and mold closure has been completed when using a HP-CRTM process. While using a high injection pressure is able to speed-up the mold-filling step, the likelihood of pressure induced problems like fiber-washing or foam core collapse is increasing. Furthermore, the investment in larger and more powerful injection units and associated equipment is significant. These disadvantages motivated the development of controlled gap variants, like PC-RTM, which proved to be of similar efficiency without the need for high cavity pressures [20]. Moreover, varying the gap during mold-filling allows to increase the permeability of the preform, reducing the filling time without the need to increase the injection pressure. Consequently, these variants are

especially suitable for manufacturing of composite designs containing any form of foamed cores that might collapse when exceeding a specific pressure level.

When using variable gap process variants, cavity pressure depends heavily on the applied closing profile and the viscosity of the resin. The latter is strongly related to the chosen resin and the reactivity of the polymer. Hence, cure cycle and closing profile must be defined according to the curing characteristic of the resin to enable fast and reliable mold-filling [19, 20]. Again, numerical simulation is vital to find a suitable processing setup, reducing the cycle time as much as possible [21].

## **2.3.2 Wet Compression Molding (WCM)**

During RTM processing, a substantial part of the cycle time is spend on draping and preforming. Furthermore, due to the separate manufacturing of the preform, additional handling steps are needed. Moreover, a significant portion of cycle time is needed for injection of the required amount of resin. The point in time after which this is finished then defines the duration needed for curing of the matrix since the resin material which entered the mold last exhibits the lowest cure degree. Wet compression molding (WCM) aims at avoiding lengthy resin injection and, if applicable, separate preform manufacturing. Thus, two WCM process variants can be differentiated as shown in Figure 2.3.

**Figure 2.3:** Two possible variants of a WCM process cycle [22]

Both process variants involve preparation of a dry fabric stack already cut to fit the dimensions of the final part. Afterwards, forming of the stack is carried out either prior to or during mold closure. The former involves preform preparation similar to that of common RTM processes. The latter, however, allows to significantly shorten the cycle time as separate preforming and associated handling steps are no longer required [23]. Furthermore, resin is applied onto either the preformed or non-preformed stack in both variants prior to introduction of both into the mold. This way, lengthy resin injection is avoided and mold-filling is initiated and controlled by the movement of the upper mold. Hence, the draping step of both variants differs in whether the resin acts as a lubricant (wet or viscous draping) or not (dry draping).

When combining draping and mold-filling within one single process step, both steps interact with each other. This results in a rather complex interaction of properties of resin and fiber textile with process design [24]. Furthermore, the draping ability of the fabric stack may be limited due to the fact that the complete stack is deformed at once. Furthermore, possibilities for the application of forces during draping is limited compared to separate preforming. Otherwise, draping of the stack in wet state can facilitate deformation and sliding of adjacent fabric layers since the resin acts like a lubricant [25]. This effect can even affect and lower the needed press force as found by Fels et al. [23]. As a result, local draping is affected by the flow front propagation and vice versa, which makes it difficult to evaluate and predict [26]. Furthermore, intermediate trimming of the preform is not possible if draping and mold-filling is done simultaneously. Hence, the stack needs to be cut such that it takes on the geometry of the final part after draping.

The disadvantages described in the previous section lead to application limits of the WCM process variant with combined draping and mold-filling. Hence, compared to traditional RTM, WCM is often applied to less complex shaped parts [27]. However, as process simulation advances, this limitation might vanish in the near future.

# **3 Reaction Kinetics of Resins**

The cross-linking process of a reactive polymer dominates its mechanical and thermal properties as it transitions from liquid to solid. In liquid state, the resin exhibits low viscosity, which is beneficial for application in processes that rely on mold filling, like RTM. During the cross-linking process, the material properties further evolve until the resin reaches its final cure state. In this state, the mechanical properties of typically applied resins are sufficiently high to be used for manufacturing of high performance composite structures. The characteristic of the curing process of a resin therefore has a huge impact on the cost-efficiency of the manufacturing process since it dominates the cycle time.

Besides its influence on the cost-efficiency of the whole process, curing also has a significant effect on many of the substeps carried out during manufacturing of a composite structure using LCM processes. As soon as the resin is injected into the hot mold, it inevitably starts to react, leading to a time and temperature dependent increase in viscosity. During subsequent curing, the resin undergoes significant morphological changes, which lead to process-induced strains due to chemical shrinkage. During cooling down of the part after demolding, additional thermal strains occur, whose magnitude is depending on the process temperature which in turn is defined by the cure cycle. As a consequence, these chemical and thermal strains lead to a distortion of the manufactured component which is strongly related to the chosen cure cycle. Moreover, the latter also determines the cure degree that is achieved at the end of the process and defines many properties of the final structure, like glass transition temperature or fracture toughness. Being able to numerically predict the curing process therefore enables to analyze and possibly optimize many aspects of composite manufacturing.

In this chapter, the reaction kinetics of resins used in RTM processes is introduced and characterization as well as mathematical modeling of this phenomenon are presented. The resulting kinetic models are used in subsequent chapters where the cure degree is required as an input parameter.

# **3.1 Introduction and review of related work**

Figure 3.1 shows a simplified schematic of the curing process of a typical thermoset polymer such as an epoxy-diamine system. Cross-linking starts with individual and initially unconnected monomers. As curing advances, these monomers form chains which grow and eventually form chains of several monomers. At some point, each chain and monomer is part of the network and the polymer gels. This is accompanied by a transition from a viscous liquid to an elastic solid and thus marks the starting point of stress formation during polymer processing. However, cross-linking is not completed at gelation. Instead, the chains further cross-link with each other until no further reaction is possible and a dense network is formed. In this state, the polymer reaches its highest possible mechanical properties and glass transition temperature.

The temporal evolution of the cure degree as well as characteristic quantities like time or cure degree at gelation are material specific. Hence, the type of kinetic model must be chosen based on the material's curing characteristic that is to be modeled. Since a multitude of materials exists, each having its own characteristic, many different models have been developed in the past. A comprehensive overview is given by Halley and Mackay [29]. In general, kinetic models can be classified in empirical and mechanistic models. The latter aims at describing the reaction based on the underlying reaction mechanisms and therefore involves comprehensive characterization and analysis of involved chemical processes. Such a detailed study of a material is usually too expensive and often also disallowed by the supplier. As a consequence, the majority of the kinetic models used for process simulations are of empirical type. This class of models does not require deep insights into the reaction mechanisms and relies

**Figure 3.1:** Simplified schematic of thermoset cure (from Menczel and Prime [28]): Initial monomer configuration (a), linear growth and branching (b), gelation (c) and fully cured thermoset (d)

on fitting empirical equations to measurement data. The challenge in this case is the appropriate choice of empirical model and the identification of suitable model parameters. An additional drawback of empirical models is the fact that a parameter set identified for a specific material is only valid for this particular resin system and cannot be used to describe other materials, even in case of similar chemical structure [30]. Furthermore, due to the different nature of isothermal and non-isothermal curing, different kinetic models may be needed [31].

The cross-linking process can be characterized using different experimental techniques which differ in the measured material property that is related to the chemical reaction. The majority of studies uses Differential Scanning Calorimetry (DSC) for this purpose and a comprehensive review of related studies is given by Yousefi et al. [32]. DSC captures the exothermal heat of reaction which is released during cross-linking, and, due to the widespread use, the procedure and data evalution is well described in the literature. Alternative characterization methods include Raman Spectroscopy, Dielectric Analysis (DEA), Fourier-Transform Infrared Spectroscopy (FTIR) and Near-Infrared Spectroscopy (NIR). Although these techniques capture much more information and details of the cross-linking process, data evaluation is more complex compared to DSC [33, 34]. Furthermore, cure evolution as gathered from DEA and Raman spectroscopy can differ from DSC data depending on the temperature [35]. Consequently, DSC is used in this study for characterization of the reaction kinetics.

LCM processes are usually carried out under isothermal temperature conditions in order to minimize cycle times. This type of processing provokes effects like vitrification if T<sup>g</sup> catches up with the curing temperature. This effect is typical for isothermal processing at temperatures below the ultimate glass transition temperature Tg,<sup>∞</sup> of the resin. Since vitrification is common for isothermal LCM processes, using isothermal DSC measurements for kinetic modeling seems to be reasonable. However, this has proved difficult due to the fact that the resin immediately starts to react when using typical process temperatures. As a result, part of the reaction happens while the DSC analyzer is still adjusting initial parameters like cell temperature and therefore does not capture any data points. Subsequent evaluation of the results is therefore difficult and prone to errors. Non-isothermal measurements are more tolerant in this regard and easier to implement since the reaction starts as soon as enough thermal energy has been introduced into the system which can be comfortably controlled by using a sufficiently low starting temperature. As a result, some researchers avoid isothermal measurements and use solely non-isothermal data for the prediction of isothermal curing processes. Spoelstra et al. [36] investigate the chemorheology of a filled epoxy system and find isothermal DSC measurements to be inadequate in case of fast curing thermosets. They even claim to be able to predict reaction kinetics for any temperature profile using their fitted model. This conclusion is only valid as long as there is no vitrification occuring, which is not necessarily the case in isothermal processes. This is also indicated by the results of Hsieh and Su [37], who achieve accurate model predictions for isothermal curing using solely non-isothermal DSC data. As soon as vitrification occurs, the prediction quality is considerably impaired, which is confirmed by the results of Atarsia and Boukhili [38] and González-Romero and Casillas [39]. The latter found isothermal measurements to be even more usable for parameter estimation compared to non-isothermal, although a higher number of measurements is required.

An important material property, which is closely related to the cure degree, is the glass transition temperature Tg. This quantity defines the thermal behavior of the material and manufactured parts. Moreover, as is shown by O'Brien et al. [40], T<sup>g</sup> has profound influence on the relaxation times and small characterization errors can therefore lead to significantly impaired simulation results. However, accurate measurement of T<sup>g</sup> is difficult due to its dependency on the used temperature rate [41]. Furthermore, T<sup>g</sup> should be measured during cooling from a state of equilibrium which is only possible if the material is fully cured. If a partial cured specimen is heated above its Tg, the value of T<sup>g</sup> is altered due to further curing. As a consequence, T<sup>g</sup> for intermediate cure states is often measured during heating within a cyclic DSC run. These measurements combine isothermal and nonisothermal temperature profiles. While the isothermal phase is used to achieve a specific cure degree, the non-isothermal part reveals the corresponding T<sup>g</sup> as well as the remaining reaction enthalpy. An effect that is often observed during such measurements is physical aging [42]. This phenomenon occurs if the resin vitrifies during curing and is exposed to this condition for long times. Physical aging manifests itself in the form of an endothermic peak that superimposes glass transition at the beginning of the remaining exothermic reaction peak. Depending on the severity, evaluation of both T<sup>g</sup> and remaining reaction enthalpy can be significantly hampered. Karkanas et al. [42] propose a rapid quench of the specimen right after the end of the endothermic physical aging peak and reheating afterwards. This technique was also succesfully applied by Wisanrakkit and Gillham [43]. However, in case of fast-curing resins, the temperature range between the end of the physical aging peak and the beginning of the reaction peak can become very small, making it in most cases impossible to eliminate

the physical aging peak without substantially altering the cure degree [44]. Vitrification during cyclic DSC runs should therefore be avoided.

#### Fast curing resins

Regarding fast-curing resins, characterization of the reaction kinetics is a difficult task and demands for tailored measurement procedures. In many cases, these resins already significantly cure before the measurement starts, resulting in an incomplete capture of the curing process [45, 46]. Furthermore, due to the possibility of different chemical reactions acting, depending on the curing temperature, characterization using a much lower temperature and extrapolating to higher temperatures using a mathematical model is prone to errors [47]. One way to reduce the amount of cure prior to the start of the measurement is speeding up the specimen preparation process and the introduction of the sample into the measurement setup, which was presented by the author in Bernath et al. [48]. Unfortunately, this automated process cannot compensate for the temperature equalization procedure inherent to DSC measurements. However, it can be taken into account and compensated for during parameter fitting of the kinetic model as demonstrated in Bernath et al. [44]. Javdanitehran et al. [49] presented an iterative fitting approach which estimates the amount of reaction enthalpy that is not captured by the DSC during heat up prior isothermal measurement. Validation of the fitted kinetic model is carried out by predicting the cure degree of various non-isothermal temperature programs and comparing these predictions with actual measurements. This approach is similar to what was published by the author in Bernath et al. [44]. However, the method presented in Bernath et al. [44] additionally integrates non-isothermal as well as cyclic DSC runs into the fitting process. Thus, the amount of experimental data used during fitting is much larger.

The problem of uncaptured curing at the beginning of isothermal measurements at elevated isothermal temperatures was also takled by Prime et al. [46]. They characterized the cross-linking process of a fast-curing urethane resin at lower temperatures and subsequently extracted a master curve from these measurements using a reduced time approach. However, this approach is not suitable if the resin vitrifies during the process since the reduced time approach is not applicable in this case.

# **3.2 Experimental characterization of reaction kinetics**

Reaction kinetics are often characterized using Differential Scanning Calorimetry (DSC). Compared to more recent analyzing methods like Dieletric Analysis (DEA), DSC offers the advantage of a simple evaluation and interpretation of the measured data. Because of this, DSC was also applied in this study.

### Differential Scanning Calorimetry (DSC)

As cross-linking advances, more and more energy is released due to bonding of reactants. A DSC analyzer is able to measure this energy by monitoring the heat flow between a sample and a well known reference material. Each material is given in an aluminum pan and both are located in the DSC cell, as is schematically shown in Figure 3.2.

Each aluminum pan is located on top of sensors which allow to accurately measure temperature and heat-flux. Furthermore, the cell can be purged using an inert gas like nitrogen, which prevents the specimen and reference from reacting with gases enclosed within the DSC furnace. The DSC cell is designed to allow for precisely controlled heating and cooling of its contents.

DSC measurements can be of the isothermal or non-isothermal type. The latter is much easier to realize since the starting temperature can be chosen so low that no reaction can take place. This enables the assumption of zero degree of cure at the start of the measurement, which is not possible in case of isothermal DSC due to time consuming temperature equalization processes before the actual measurement starts [44, 45]. The drawback of non-isothermal DSC, however, is that the process in reality is carried out under isothermal conditions. Therefore, effects like premature vitrification will occur at some point during the process and

**Figure 3.2:** Principle of DSC analysis

will substantially slow down the cure rate due to the fact that after vitrification, the reaction has to take place in the glassy material state [39, 44, 50, 51]. This effect is not present in non-isothermal measurements. Therefore, using solely non-isothermal data for model parametrization will lead to incorrect isothermal cure predictions, regardless of whether the model is in principle capable of describing the effect of premature vitrification [44]. For RTM process simulation, it is therefore vital to have both, measurement data which shows influence of vitrification as well as a suitable kinetic model that is capable of taking into account vitrification.

Figure 3.3 exemplarily shows the heat flow of a reacting epoxy resin using a constant heat rate of 5 K min−1. The offset in the beginning is caused by the glass transition of the uncured polymer, Tg,0. Due to the very low starting temperature, the specimen is in a glassy state prior to heating and after exceeding Tg,0, the material transitions to rubbery state. Shortly after that, the reaction starts as is visible from the change in slope which eventually forms the characteristic reaction peak.

**Figure 3.3:** Heat flow evolution during cure of epoxy resin from non-isothermal DSC [44]

The reaction peak is what describes the temporal evolution of the cross-linking process. The cure degree α can be evaluated from such data by integration of the heat flow signal over time. However, the measurement data as shown in Figure 3.3 does not only contain heat flow from the reaction but also from constantly heating the specimen, which depends on the specific heat capacity of the material. Therefore, prior to integration, the reaction peak needs to be isolated from the measured heat flow by constructing a baseline. For the hypothetical case that no heat is released during cross-linking, aside from the glass transition at the beginning, the resulting curve would be mostly flat with a smooth step that results from a change in specific heat capacity due to cross-linking. Since it is not possible to measure this baseline using the standard DSC technique, it has to be constructed based on the characteristic of the measured heat flow. A sigmoidal baseline is often used, especially if there is a significant difference in height before and after the reaction peak. Otherwise, a straight line is sufficient (cf. Figure 3.3).

The cure degree describes how much of the cross-linking has already happened. This is equivalent to how much of the specific total reaction enthalpy Δh has already been released by the specimen during a non-isothermal DSC run. The specific reaction enthalpy h is given by the area enclosed by the specific heat flow curve q and the baseline b:

$$h\left(t\right) = \int\_{t\_0}^{t} q\left(\hat{t}\right) - b\left(\hat{t}\right)\,\mathrm{d}\hat{t}.\tag{3.1}$$

Here, t<sup>0</sup> denotes the time at which the reaction starts. The total reaction enthalpy Δh is obtained from this equation by integrating over the whole peak. The temporal evolution of the cure degree α (t) is then given by

$$\alpha\left(t\right) = \frac{1}{\Delta h} \cdot h\left(t\right) = \frac{1}{\Delta h} \cdot \int\_{t\_0}^{t} q\left(\hat{t}\right) - b\left(\hat{t}\right) \,\mathrm{d}\hat{t} \tag{3.2}$$

and the corresponding cure rate by

$$\frac{d\alpha}{dt} = \frac{q\left(t\right) - b\left(t\right)}{\Delta h}.\tag{3.3}$$

The evolution of the cure degree based on the data given in Figure 3.3 is shown in Figure 3.4. Since the underlying DSC measurement is of the non-isothermal type, no vitrification occurs and the material fully cures by reaching a cure degree of α = 1.

Isothermal DSC measurements can be evaluated the same way as is described for non-isothermal DSC above. However, since the temperature at the beginning of an isothermal run is already high enough to enable cross-linking, there is no well defined reaction peak to evaluate. Furthermore, in case of moderate to highly reactive resins, a significant amount of curing can already take place during heating of the DSC cell up to the desired measurement temperature and subsequent temperature equalization before the actual measurement starts. The end of an isothermal measurement is characterized by a vanishing cure rate which is either caused by reaching a fully cured state or by a transition from rubbery to

**Figure 3.4:** Cure degree evolution evaluated from non-isothermal DSC data [44]

glassy state (vitrification). This therefore depends on the interaction of reaction temperature T and glass transition temperature Tg. If T <Tg, the material vitrifies to a glassy state and otherwise remains in rubbery state. However, since it is a smooth rather than a sudden transition, vitrification may already start after a critical value of T − T<sup>g</sup> has been reached. During isothermal measurement, the reaction temperature remains constant whereas T<sup>g</sup> constantly increases due to the ongoing cross-linking. The dependency of T<sup>g</sup> on the cure degree can be characterized by using cyclic DSC measurements, being basically a combination of isothermal and non-isothermal DSC. The idea is to cure the resin isothermally for a specific period of time and than stop the reaction by rapidly cooling down to a much lower temperature. Afterwards, the specimen is again heated using a constant heating rate in order to reach a fully cured state. During this second stage, a glass transition can be observed similar to Tg,<sup>0</sup> of a purely non-isothermal run. This time, however, it represents the T<sup>g</sup> of a partially cured resin where the actual cure state depends on the chosen isothermal temperature and dwell time. Consequently, many different cure states can be realized by varying these values

which enables the characterization of T<sup>g</sup> as a function of the cure degree. Since the first isothermal stage of the measurement cannot be evaluated accurately due to premature curing before the actual measurement, the cure degree which belongs to the measured T<sup>g</sup> is evaluated from the remaining reaction enthalphy ΔhR, that is released during the non-isothermal stage:

$$\alpha = \frac{\Delta h - \Delta h\_{\rm R}}{\Delta h} = 1 - \frac{\Delta h\_{\rm R}}{\Delta h}.\tag{3.4}$$

Evaluation of Δh<sup>R</sup> requires that the reaction peak can be clearly separated from other temperature driven effects which can occur during a DSC measurement. An example of such an effect, which can make data evaluation difficult to impossible, is enthalpy relaxation due to physical aging [52]. This effect usually occurs when reheating a resin sample that vitrified during a measurement and was kept in this state for long periods of time. Furthermore, it depends on the cure degree as reported by Montserrat [53]. This problem can be solved by interrupting the non-isothermal part of the measurement right after the endothermic physical aging effect has been relieved. Afterwards, the sample is rapidly cooled down and heated up again as originally intended showing no influence of physical aging. This solution, however, is only applicable to rather slow curing resins. If the reactivity is too high, it is impossible to relieve physical aging effects without altering the cure state of the specimen. It is therefore beneficial to use multiple temperature and dwell time combinations rather than a single temperature and many different dwell times to avoid physical aging during cyclic DSC measurements. Moreover, when evaluating T<sup>g</sup> from heating experiments, the applied scanning rate should be chosen as small as possible to avoid shifting due to thermal lag and other rate dependent influences as shown by Schawe [41].

### **3.2.1 SCER1 epoxy resin**

The following section describes the characterization of the cure kinetics of SCER1 resin. Some aspects of this have already been addressed by the author in Bernath et al. [44].


**Table 3.1:** Parameters of isothermal and non-isothermal DSC measurements (SCER1)

Table 3.1 contains the temperatures and heating rates of isothermal and nonisothermal DSC measurements which were used for characterization of SCER1 resin. DSC measurements were carried out at the Institute for Chemical Technology and Polymer Chemistry (ITCP) at the Karlsruhe Institute of Technology (KIT) in Karlsruhe, Germany. A DSC Q200 analyzer from TA Instruments (New Castle, DE, USA) and nitrogen atmosphere were used.

#### Isothermal DSC

Figure 3.5 shows the cure rate of isothermal DSC runs of SCER1 resin. The highest cure rate is observed at the very beginning of each measurement, which is what makes the evaluation of such measurements prone to errors due to premature curing as described in Section 3.2. Moreover, the maximum cure rate increases significantly when increasing the cure temperature. Hence, the problem of unrecorded premature curing will also intensify with temperature. In contrast, a lowered cure temperature leads to a decrease in cure rate. While this lowers the risk of unrecorded curing, the cross-linking reaction takes much longer.

The evolution of the cure degree is calculated by applying Equation 3.2 to the cure rate data. The result is shown in Figure 3.6. In accordance to the cure rate data given in Figure 3.5, the cure degree raises faster at the beginning as the isothermal temperature of the measurement is increased. However, the achieved maximum cure degree of each temperature exhibits an unexpected behavior, showing lower levels of cure degree in case of higher curing temperatures. As an example, the calculated final cure degree of the resin achieved using a cure temperature of 120 ◦C is much lower than that of the 60 ◦C run. This is caused by premature curing which occurs before the start of the measurement and therefore

is not recorded by the DSC analyzer. This issue must be taken into account during model fitting.

**Figure 3.5:** Cure rate evolution of SCER1 resin for different isothermal temperatures

#### Non-isothermal DSC

Figure 3.7 shows the cure rate of SCER1 resin for all investigated heating rates. Since the cure rate depends strongly on temperature, the maximum cure rate is higher for increased heating rates due to the fact that higher temperatures are reached in a shorter period of time. When comparing the temperature ranges during which the reaction takes place for each individual heating rate, it is found that the onset temperatures of the reaction peak are more or less equal whereas the endset temperature is shifted to higher values. As a result, the cross-linking process is spread over a much wider temperature range for increased heating rates, leading to a shift of the peak temperature towards higher values. This effect consequently limits the applicable maximum heating rate since the resin might suffer from thermal degradation if too high temperatures are reached, affecting both material properties and measurement data.

**Figure 3.6:** Cure degree evolution of SCER1 resin for different isothermal temperatures

**Table 3.2:** Total reaction enthalpy from non-isothermal DSC measurements (SCER1)


The evolution of the cure degree is calculated from the cure rate data by applying Equation 3.2. The corresponding curves are given in Figure 3.8. The calculated total reaction enthalpies Δh of each heat rate are given in Table 3.2. From these values, a mean <sup>Δ</sup><sup>h</sup> of 499.8 <sup>±</sup> 16.6 J g−1 is found for SCER1 resin. The shift of the cross-linking process to higher temperature ranges in case of high heating rates leads to a similar shift of the cure degree evolution. Similarly, slow heating rates shift the cross-linking process to low temperatures. Depending on the reactivity of the resin, this can lead to temporary vitrification which is observed in form of a temporary decrease in cure rate in case of the slowest heating rate of 1 K min−1 between about 90 ◦C and 130 ◦C. The reason for this is the increase in T<sup>g</sup> as cross-linking advances. As soon as T<sup>g</sup> comes close to

**Figure 3.7:** Cure rate evolution of SCER1 resin for different heating rates

the current temperature of the DSC cell, the cure rate is slowed down due to vitrification which leads to a flattening of the cure degree evolution. This effect is closely related to the cure degree dependency of the glass transition temperature, which is investigated in the following section.


**Table 3.3:** Cure cycles used for cyclic DSC measurements of SCER1 resin

\* Vitrification is observed during isothermal curing.

**Figure 3.8:** Cure degree evolution of SCER1 resin for different heating rates

# Evolution of *<sup>T</sup>***g** during cross-linking

The cure dependency of T<sup>g</sup> is analyzed using cyclic DSC measurements as introduced at the beginning of Section 3.2. The cure cycles used for cyclic DSC measurements are given in Table 3.3. As reported by the author in Bernath et al. [44], vitrification is observed in case of SCER1 resin during isothermal curing for the longest hold time of each curing temperature. As a result, the glass transition at the beginning of the non-isothermal stage is superimposed by enthalpy relaxation due to physical aging. Consequently, evaluation of T<sup>g</sup> as well as the remaining reaction enthalpy Δh<sup>R</sup> may be impaired [44]. To demonstrate this issue, Figure 3.9 exemplarily shows the non-isothermal part of the cyclic DSC runs conducted using an isothermal curing temperature of 80 ◦C. The glass transition, which is visible at the beginning of the heating process, is shifted to higher temperatures as the dwell time tiso is increased. Afterwards, the reaction peak of the remaining cross-linking process begins to develop. The problem of physical aging is most clearly seen in the result of the longest dwell time of 30 min, which shows a strong endothermic effect that superimposes both, glass transition and begin of the reaction peak. Hence, evaluation of the remaining reaction enthalpy is prone to errors since identification of the beginning of the remaining cross-linking reaction is hardly possible. Contrary to this, the glass transition still shows the typical step-like character and therefore can still be evaluated, although its value might be affected by a shift due to phyiscal aging. Since it is favorable to have T<sup>g</sup> data of vitrified material available during kinetic model fitting, cyclic measurements affected by vitrification are evaluated using the strategy proposed in Bernath et al. [44], which makes use of a Tg-model fitted using measurements showing no influence of vitrification in order to reverse calculate the cure degree of measurements affected by vitrification.

**Figure 3.9:** Non-isothermal part of cyclic DSC measurements of SCER1 resin cured using different hold times at an isothermal temperature of Tiso = 80 ◦C

Figure 3.10 shows the dependency of T<sup>g</sup> on the cure degree of SCER1 resin. Data points sharing the same color all belong to the same isothermal temperature and only differ by the applied dwell time. Measurements affected by vitrification are not included in this plot since their evaluation involves the use of a fitted Tg-model, which is therefore shown in Section 3.5.1. In general, the data show an increasing steepness towards high cure degrees that make T<sup>g</sup> more sensitive to small changes in cure degree in this region.

**Figure 3.10:** Evolution of <sup>T</sup>g of SCER1 resin during cross-linking from cyclic DSC measurements

The values of T<sup>g</sup> of uncured and fully cured resin, Tg,<sup>0</sup> and Tg,∞, respectively, are a by-product of non-isothermal DSC runs. The former is found at the very beginning of a non-isothermal measurement if a sufficiently low enough temperature is chosen as the start temperature. Prior to any chemical reaction, a step in the DSC signal is observed, caused by the glass transition of the uncured resin. In case of resin SCER1, this happens at a temperature of −30.64 ◦C. The value of Tg,<sup>∞</sup> is found by reheating a specimen that has been fully cured in a previous non-isothermal run, revealing an ultimate Tg,<sup>∞</sup> of resin SCER1 of 132.85 ◦C.

### **3.2.2 FCPUR polyurethane resin**

The following section describes the characterization of the cure kinetics of the fast-curing polyurethane (FCPUR) resin. Table 3.4 contains the temperatures and heating rates of isothermal and non-isothermal DSC measurements. Measurement data for FCPUR resin was provided by the manufacturer as part of a collaboration within the project SMiLE [54]. Due to the high reactivity of the resin, isothermal DSC measurements were not feasible.

**Table 3.4:** Parameters of isothermal and non-isothermal DSC measurements (FCPUR)


\* Infeasible due to high reactivity of the resin.

#### Non-isothermal DSC

The cure rate of non-isothermal DSC runs using different heating rates of FCPUR resin are shown in Figure 3.11. Two measurements have been conducted for each heating ramp of which the second one is represented by dashed lines. The much higher reactivity of this resin becomes obvious by comparing 10 K min−1 curves of Figures 3.11 and 3.7. While the cure rate of resin SCER1 vanishes at roughly 220 ◦C, this happens already at about 130 ◦C in case of resin FCPUR. Furthermore, the peak value of the cure rate is also much higher for the latter. As a result, isothermal measurements are very difficult given the fact that this was already the case for the much slower resin SCER1.

**Table 3.5:** Total reaction enthalpy from non-isothermal DSC measurements (FCPUR)


**Figure 3.11:** Cure rate evolution of FCPUR resin for different heating rates

The resin shows a significant heating rate dependency of the total reaction enthalpy, which was the original reason for the repetition of the measurements. The corresponding values are given in Table 3.5. The results of the first run exhibit a strong heat rate dependency leading to significantly different enthalpies of each individual heat rate. The second run confirmed these results with the exception of the 10 K min−1 measurement which exhibits a value close to the one of the 20 K min−1 run. A possible reason for this behavior may be a temperature driven change in reaction mechanism as investigated by Karkanas et al. [42]. They propose plotting the normalized reaction rate over the fractional conversion to reveal any abnormal behavior, which is shown in Figure 3.12. This reveals a possible change in reaction mechanism of the resin prior to a fractional conversion of 0.6. However, the effect is much less pronounced as for instance found by Karkanas et al. [42] in case of the F934 resin which shows a much more distinctive change in reaction rate.

Due to the varying reaction enthalpy of resin FCPUR, the evolution of the cure degree from non-isothermal DSC data is calculated using the individual corresponding total reaction enthalpy rather than a mean value of all conducted

**Figure 3.12:** Evaluation of change in reaction mechanism of FCPUR resin

measurements. This is based on the assumption that irrespective of what reaction mechanism is acting, the resin eventually reaches full cure. The resulting evolution of the cure degree for the investigated heating rates is given by Figure 3.13.

**Table 3.6:** Cure cycles used for cyclic DSC measurements of FCPUR resin


\* Vitrification is observed during isothermal curing.

# Evolution of *<sup>T</sup>***g** during cross-linking

Characterization of the cure dependency of T<sup>g</sup> is carried out using a similar approach as is shown above for resin SCER1. The cyclic DSC measurements

**Figure 3.13:** Cure degree evolution of FCPUR resin for different heating rates

were again conducted by the resin's manufacturer and the results were given to the author as part of the collaboration within the project SMiLE.

The cure cycles used for cyclic DSC measurements are given in Table 3.6. In contrast to the amount of cyclic measurements of SCER1 resin, only a limited set of measurements is available. The heatflow evolution of the non-isothermal part of the conducted cyclic DSC runs at 80 ◦C is shown in Figure 3.14. The increase in T<sup>g</sup> with increasing isothermal hold times is clearly visible at the beginning of each curve. The opposite applies to the area of the peak of the remaining crosslinking reaction. The measurement with the longest hold time of 3600 s shows a characteristic typical for an influence of physical aging which is also observed in cyclic DSC runs of SCER1 resin. Compared to the other measurements, it is much harder to separate glass transition and begin of reaction peak from each other. Hence, similar to resin SCER1, the cure degree associated with this measurement is reversely calculated using the Tg-model fitted to the remaining data points.

**Figure 3.14:** Non-isothermal part of cyclic DSC measurements of FCPUR resin cured using different hold times at an isothermal temperature of Tiso = 80 ◦C

Figure 3.15 shows the evolution of T<sup>g</sup> of resin FCPUR with advancing cure degree. Since the start temperature of the non-isothermal runs was chosen too high, the T<sup>g</sup> of the uncured resin, Tg,<sup>0</sup> is not present in the corresponding heatflow curves. Hence, Tg,<sup>0</sup> will be interpreted as a fitting parameter and quantified during model parametrization. This does not pose a problem since a loss in accuracy of the T<sup>g</sup> evolution prior to gelation is negligible. The opposite applies to the value of the fully cured material, Tg,∞. However, evaluation of this property proved to be difficult since the distinctive step in heatflow is not present in the corresponding DSC curve. Instead, this value is derived from Thermomechanical Analysis (TMA) measurements as presented in the following section.

**Figure 3.15:** Evolution of <sup>T</sup>g of FCPUR resin during cross-linking from cyclic DSC measurements

# **3.3 Determination of glass transition temperature**

Characterization of T<sup>g</sup> can be achieved using different techniques. However, only some of them allow to measure the T<sup>g</sup> without any influence of frequency or strain rate, as is typically the case when using Dynamic Mechanical Analysis (DMA) or rheometry. Thermomechanical Analysis (TMA) is based on measuring the dimensional change of a specimen induced by temperature variations, thus neither frequency nor strain rate needs to be superimposed [28]. Hence, TMA measurements, which were originally conducted to measure the coefficient of thermal expansion of the resins, are also evaluated with respect to T<sup>g</sup> in this section. The measurements were carried out at the Institute of Composite Structures and Adaptive Systems at the German Aerospace Center (DLR) in Braunschweig, Germany.

The principle of extracting T<sup>g</sup> from TMA data is based on the change in thermal expansion as the temperature exceeds or undercuts Tg. Moreover, thermal strains increase significantly as the resin transitions from glassy to rubbery state. Hence, the temperature at which this change is observed is attributed to Tg. Figure 3.16 shows the evolution of thermal strain during heating with a constant rate of 3 K min−1, measured using a TMA analyzer. The linear section up to a temperature of about 95 ◦C is dominated by the glassy thermal expansion coefficient of the resin. Above that temperature, the thermal strain shows again a linear region with a much steeper slope, caused by the significantly higher thermal expansion of the material in rubbery state. Hence, the T<sup>g</sup> is defined by the temperature at the inflection point, which is evaluated by constructing tangents to the linear regions above and below Tg. From five identical measurement runs, a Tg,<sup>∞</sup> of 95.77 ± 0.70 ◦C is found for resin FCPUR. The abnormal behavior of the material at a temperature of 105 ◦C to 110 ◦C and above is caused by outgassing and blistering of the resin, which interferes with the strain measurement.

**Figure 3.16:** TMA measurements showing the thermal strain of five repeated trials of fully cured FCPUR resin during heating with a constant rate of 3 K min−1

The evaluation procedure shown above for resin FCPUR is also applied to other resins, of which some were part of the study in Bernath et al. [15]. However, in the mentioned study, values given in corresponding technical data sheets have been used as Tg,<sup>∞</sup> of the resins. The error introduced by relying on these manufacturer provided values must be compensated for using a model parameter as described in the mentioned study. In this thesis, however, more effort is put in correct characterization of Tg,<sup>∞</sup> to avoid the use of calibration parameters in the models used for predicting process-induced distortion. This way, the resulting numerical method is more flexible and robust.

The specimens used for TMA measurements have been cured using the MRCC beforehand. Hence, during a first heating run, the T<sup>g</sup> associated with that cure cycle is measured. As soon as the resin switches to rubbery state, cross-linking continues until either full cure is achieved or the material vitrifies again. The latter is prevented by using a rather slow temperature rate as well as a high end temperature of each heating step. As a result, heating of the same specimen for a second time reveals the ultimate glass transition temperature Tg,<sup>∞</sup> of the resin. This is shown together with the individual cure temperatures as defined by the MRCC in Figure 3.17. Interestingly, when applying the MRCC, all resins but SCER1 show a T<sup>g</sup> that is very close to the cure temperature. Moreover, Tg,<sup>∞</sup> of resins FCER1 and SCER1 is found to be in the same range as T<sup>g</sup> (Tc). This is caused by an already very high value of T<sup>g</sup> after curing with the MRCC, which implies that these resins already reached full cure prior to TMA measurements.

All fast-curing resins included in the presented TMA study show a T<sup>g</sup> (Tc) which is very similar to the cure temperature defined by the MRCC. Hence, although no TMA measurements have been conducted for resin FCER3, it is assumed that this fast-curing resin behaves similarly.

# **3.4 Determination of the cure degree at gelation**

The cure degree at gelationαgel determines at which point the material can sustain any stresses during the course of the process. Hence, it defines the beginning of the development of relevant process-induced strain, eventually leading to residual

**Figure 3.17:** Glass transition temperature <sup>T</sup>g of resins cured using the MRCC and measured using TMA technique. During a first heat-up, the <sup>T</sup>g associated with the MRCC is measured, while during a second run, the ultimate value, <sup>T</sup>g,∞, is revealed.

stress. Different methods exist for characterization of this important property. Often, gelation is quantified using an oscillatory rheometer measurement and by identifying the period of time, after which the curves of storage and loss modulus cross. This period of time is then correlated to the cure degree evolution using DSC measurements or a fitted kinetic model. However, this method can lead to errors if the cure degree at gelation is rather high and therefore may be superimposed by the onset of vitrification, which has a similar effect on the storage modulus as gelation. Another criteria for gelation can be formulated based on the evolution of the phase angle, measured using different frequencies. The point at which these curves intersect and the phase angle becomes independent of frequency is related to gelation. This method was used for determination of the time of gelation of all resins but FCER3, as detailed in Bernath et al. [15]. The measurements are conducted using a constant heating rate of 15 K min−1 within a temperature range of 25 ◦C to 250 ◦C. The actual cure degree is related to the determined time of gelation using DSC measurements carried out using the very same temperature profile, and by applying Equation 3.2. Rheometer and DSC measurements were both carried out at the Institute of Composite Structures and Adaptive Systems at the German Aerospace Center (DLR) in Braunschweig, Germany.

Resin FCER3 was not part of the study presented in Bernath et al. [15]. Hence, gelation of that resin is characterized using an alternative method which makes use of the fact that gelation marks the onset of stress evolution by monitoring the force caused by the ongoing chemical shrinkage of the resin. Prior to gelation, the resin is able to flow and thus can counteract any stress build-up by viscous flow. After gelation, the polymer network is cross-linked to such a high degree that this is not possible anymore. Hence, the strain due to chemical shrinkage induces stress and thus force as soon as gelation is surpassed. This can be measured using parallel-plate rheometry. During such a measurement, the analyzer keeps the gap between the plates at a constant value. Consequently, the period of time after which the force acting on the plates deviates from zero is related to gelation.

Figure 3.18 shows the measured cure degree at gelation of the three resins used in this thesis as well as that investigated in the related study [15]. Resin FCPUR shows the lowest value of 51.7 %. The epoxy resins exhibit higher cure degrees at gelation, ranging from 54.2 % to 79.45 %. Since the cure degree is calculated using a fitted kinetic model, errors in cure degree prediction are directly transferred to the determination of gelation. Hence, the rather abnormally high value of resin FCER3 might be caused by such a shortcoming of the kinetic model since the measurements on which the provided parameter set founds on were not accessible. Furthermore, the reactivity of the resin has no distinctive influence on the cure degree at gelation.

**Figure 3.18:** Cure degree at gelation of the resins used in this study. For reasons of comparison, resins used in Bernath et al. [15] are also included.

# **3.5 Mathematical modeling of reaction kinetics**

Reaction kinetics models can be categorized in mechanistic and phenomenological models. The former class of models aims at describing the reaction based on the fundamental chemical processes that take place during curing. Hence, they require detailed information of the cross-linking process which is often kept secret by the manufacturers or would require extensive characterization effort. Furthermore, even when modeling the individual steps of a reaction, reuse of a fitted kinetic model for another resin with an equal reaction mechanism is not guaranteed [30]. Therefore, in most cases, the second class of models is applied which usually relies on (semi-)empirical equations. The most simple type is named nth-order models and is based on a function like

$$\frac{d\alpha}{dt} \propto \left(1 - \alpha\right)^{n}.\tag{3.5}$$

In this equation, n denotes the reaction order. An important property of this equation is the fact that at the beginning of the reaction (α = 0), the cure rate is non-zero. Relations of this type can be used to describe the curing process of resins under purely isothermal conditions or in cases where several reactions occur simultaniously [42]. However, the majority of thermoset resins, used e.g. for RTM processing, exhibits a self-accelerating characteristic. Such autocatalytic reactions are better described using an equation proposed by Kamal and Sourour [55]:

$$\frac{d\alpha}{dt} \propto \alpha^m \cdot \left(1 - \alpha\right)^n. \tag{3.6}$$

In this case, the cure rate is zero at the beginning of the reaction and will reach its maximum at some point during cross-linking, where the self-acceleration due to the term α<sup>m</sup> is most pronounced. An often applied model of this type is the Kamal-Malkin kinetic model, which is also referred to as Kamal-Sourour kinetic [55, 56]:

$$\frac{d\alpha}{dt} = \left(K\_1 + K\_2 \cdot \alpha^m\right) \cdot \left(1 - \alpha\right)^n. \tag{3.7}$$

This model has been widely applied to different resin systems in the literature [31, 44, 48, 57–61]. Karkanas et al. [31] further modified this equation by decoupling the reaction orders of both reactions, which can yield better accuracy due to the additional degree of freedom:

$$\frac{d\alpha}{dt} = K\_1 \cdot (1 - \alpha)^{n\_1} + K\_2 \cdot \alpha^m \cdot (1 - \alpha)^{n\_2} \,. \tag{3.8}$$

The reaction rate constant K<sup>1</sup> and K<sup>2</sup> are usually modeled using an Arrhenius type equation:

$$K\_i = A\_i \cdot \exp\left(-\frac{E\_i}{\mathbf{R} \cdot T}\right). \tag{3.9}$$

Here, Ai, E<sup>i</sup> and R denote pre-exponential factor, activation energy and gas constant, respectively.

The above mentioned kinetic models are capable of describing a reaction whose cure rate is only limited by chemical processes, which is referred to as chemically controlled. However, as soon as the material vitrifies, e.g. because T<sup>g</sup> reaches a value similar to the current curing temperature, the reaction becomes limited by the ability of the reactants to move [53]. In this case, the reaction rate is diffusion controlled and therefore strongly depending on temperature and time. At this point, predictions of the cure degree start to diverge from experimental data if vitrification is not taken into account by the kinetic model [42]. The latter can be obtained by introducing a maximum achievable cure degree, which depends on the reaction temperature and limits the cure rate as soon as the cure degree reaches that specific value. Shin and Hahn [62] used such an approach for cure kinetic modeling of AS4/3502 resin. However, limiting the cure rate this way can lead to a sudden stop of the predicted cross-linking process. In reality, a more gradual transition from chemical to diffusion controlled reaction is observed. Karkanas et al. [31] successfully applied a diffusion controlled approach based on a splitting of the reaction rate constants:

$$\frac{1}{K\_i} = \frac{1}{K\_{i,\mathbf{c}}} + \frac{1}{K\_{\mathbf{d}}}.\tag{3.10}$$

This equation switches between a chemically controlled reaction rate Ki,<sup>c</sup> and a diffusion controlled reaction rate Kd. While the former is modeled using Equation 3.9, modeling of K<sup>d</sup> can be done using different approaches. Simon and Gillham [63] as well as Karkanas et al. [31] use an Arrhenius type equation, while Wisanrakkit and Gillham [43] use an equation based on the Williams-Landel-Ferry (WLF) approach. Garschke et al. [64] on the other hand introduce a prefactor into Equation 3.7 which gradually fades out the cure rate. Grindling [65] introduced a model for the diffusion controlled reaction rate which switches between three states depending on the distance between T<sup>g</sup> and the reaction temperature. This model offers great freedom for modeling the transition from chemically to diffusion controlled state and is therefore used in the present thesis and related publications [44, 61].

While the model of Karkanas assumes both reaction rate constants to be affected by diffusion, Grindling's model limits this effect to the rate constant K<sup>2</sup> by replacing it with an effective rate constant Keff:

$$\frac{d\alpha}{dt} = K\_1 \cdot (1 - \alpha)^{n\_1} + K\_{\text{eff}} \cdot \alpha^m \cdot (1 - \alpha)^{n\_2} \,. \tag{3.11}$$

Keff is used to switch between chemically and diffusion controlled reaction using an equation similar to Equation 3.10:

$$\frac{1}{K\_{\text{eff}}} = \frac{1}{K\_2} + \frac{1}{K\_{\text{d}}}.\tag{3.12}$$

The value of the diffusion controlled reaction rate is calculated with respect to the current temperature and Tg:

$$\begin{aligned} T &> (T\_{\rm g} + \Delta T\_{\rm g}) ;\\ K\_{\rm d} &= K\_{\rm d} \left( T = T\_{\rm g} + \Delta T\_{\rm g} \right) \cdot \exp \left( -\frac{E\_{2,\rm d}}{\rm R} \cdot \left( \frac{1}{T} - \frac{1}{T\_{\rm g} + \Delta T\_{\rm g}} \right) \right) & \quad \text{(3.13)} \\\\ T\_{\rm g} &\leq T \leq (T\_{\rm g} + \Delta T\_{\rm g}) ; \end{aligned}$$

$$K\_{\rm d} = K\_{\rm d}^{T = T\_{\rm g}} \cdot \exp \left( \frac{c\_1 \cdot (T - T\_{\rm g})}{c\_2 + T - T\_{\rm g}} \right) \end{aligned} \tag{3.14}$$

T <Tg:

$$K\_{\rm d} = K\_{\rm d}^{T=T\_{\rm g}} \cdot \exp\left(-\frac{E\_{\rm 1,d}}{\rm R} \cdot \left(\frac{1}{T} - \frac{1}{T\_{\rm g}}\right)\right) \tag{3.15}$$

Equations 3.13 and 3.15 are of Arrhenius type whereas Equation 3.14 is based on a WLF approach and enables a smooth transition between glassy and rubbery state using the parameters c<sup>1</sup> and c2.

The activation energies E1,<sup>d</sup> and E2,<sup>d</sup> are given by the following set of equations:

$$E\_{1, \rm d} = T\_{\rm g}^2 \cdot \frac{c\_1}{c\_2} \tag{3.16}$$

$$E\_{2, \rm d} = \frac{c\_1 \cdot c\_2 \cdot \left(T\_\text{g} + \Delta T\_\text{g}\right)^2}{\left(c\_2 + \Delta T\_\text{g}\right)^2} \tag{3.17}$$

Compared to the simpler models like the Kamal-Malkin model, the number of model parameters, which need to be identified, is much higher in this case. Furthermore, a large part of the equations used in the Grindling model depend on Tg. Therefore, a model for the cure dependency of T<sup>g</sup> is needed. This relationship is modeled using the approach originally developed by DiBenedetto [66, 67] and which has been often applied in the literature [31, 63, 68]:

$$\frac{T\_{\rm g} - T\_{\rm g,0}}{T\_{\rm g,\infty} - T\_{\rm g,0}} = \frac{\lambda \cdot \alpha}{1 - (1 - \lambda) \cdot \alpha}.\tag{3.18}$$

This equation can be adopted to experimental measurements of the non-linear relationship between cure degree and glass transition temperature using the model parameter λ.

### **3.5.1 SCER1 epoxy resin**

The suitability of a specific kinetic model depends on the reaction mechanisms that participate in the cross-linking process. If more than one reaction mechanism is acting, it is in general possible that one of them dominates the others. This usually depends on temperature as investigated by Karkanas et al. [42]. They propose to check for this by plotting the normalized reaction rate over the fractional conversion, as shown for SCER1 resin in Figure 3.19. From this, no change in reaction mechanism is visible. According to Karkanas et al. [42], resins which do not show a change in reaction mechanism can be modeled using rather simple kinetic models. However, this conclusion is only valid in the absence of vitrification since vitrification does affect the cure rate but is not necessarily followed by a change in reaction mechanism. This effect is present in case of the slowest heat rate in Figure 3.19 and is the reason for the deviant behavior at the end of the cross-linking process. Thus, all conclusions drawn from such a plot are only valid if no vitrifaction occurs.

**Figure 3.19:** Evaluation of change in reaction mechanism of SCER1 resin

The fitting process for SCER1 resin is carried out using the DSC measurements that have been discussed in Section 3.2.1. As was demonstrated by the author in [44], using as much information as available during the fitting process is beneficial to the result. Furthermore, the DSC data contains much information about the influence of vitrification on the cross-linking process. Since this effect is also present in case of cure cycles typically applied in RTM processing, it is favorable to use the Grindling model for describing the cross-linking process of this resin. Due to the fact that vitrification depends on the interaction of cure temperature and Tg, modeling of the latter is inevitable. For this purpose, the DiBenedetto equation given by Equation 3.18 is used to fit the experimentally measured T<sup>g</sup> values of different cure degrees.


**Table 3.7:** Parameters used for modeling the cure dependency of <sup>T</sup>g of resin SCER1 based on the DiBenedetto model

Fitting of the DiBenedetto Tg-model involves identification of a single parameter λ, which controls the curvature of the cure degree dependency of Tg. The remaining model parameters, Tg,<sup>0</sup> and Tg,<sup>∞</sup> have been directly derived from DSC runs as described in Section 3.2.1. Table 3.7 summarizes the parameters used for predicting T<sup>g</sup> of resin SCER1 in this thesis. Furthermore, Figure 3.20 shows a comparison of the fitted model and experimentally measured T<sup>g</sup> values. In this plot, filled symbols refer to cyclic DSC measurements affected by vitrification.

**Figure 3.20:** Comparison of fitted <sup>T</sup>g-model and experimentally measured values of resin SCER1. Values gathered from DSC measurements showing influence of physical aging are marked by filled symbols.

As discussed in Section 3.2.1, fitting of the Tg-model is first carried out without taking into account results of measurements affected by physical aging due to the high uncertainty in evaluation of the remaining reaction enthalpy. Afterwards, the cure states belonging to the T<sup>g</sup> values of these measurement are reversely calculated using the fitted Tg-model. In a second iteration, the Tg-model is fitted again, using also vitrification affected DSC data. This procedure enables to utilize them also in fitting of the kinetic model. The importance of these measurements for the subsequent fitting process of the kinetic model arises from the fact that they contain information about the maximum degree of cure that can be achieved using a specific isothermal cure temperature.

**Figure 3.21:** Distance between curing temperature and <sup>T</sup>g (<sup>T</sup> <sup>−</sup> <sup>T</sup>g) during non-isothermal DSC of SCER1 resin for different heating rates

The knowledge of the dependency of T<sup>g</sup> on the cure degree enables further evaluation of the DSC results. The occurrence of temporary vitrification in case of the lowest heat rate of 1 K min−1 becomes evident when plotting the interaction of cure temperature and Tg, denoted as T − Tg, of non-isothermal DSC over the curing temperature. As shown in Figure 3.21, T<sup>g</sup> equals the reaction temperature over a temperature range of about 40 ◦C, leading to temporary vitrification.


**Table 3.8:** Parameters used for modeling the cure degree evolution of resin SCER1 based on the Grindling kinetic model

Fitting of the Grindling kinetic model is carried out using isothermal, nonisothermal and cyclic DSC data. Since many of these measurements contain information about vitrification, the model will be able to take this effect into account as shown in Bernath et al. [44]. Table 3.8 lists the identified model parameters for modeling the cure evolution of resin SCER1 based on the Grindling kinetic model. A comparison of non-isothermal DSC measurements and numerically predicted cure evolutions is given by Figure 3.22. In general, good agreement is found with the exception of the tendency of the parametrized model to overestimate the cure degree towards the end of the curing process in case of all heating rates. In Bernath et al. [44], this model prediction is compared to predictions of the same model fitted using only non-isothermal DSC data. The latter revealed much better accuracy for this type of measurement. However, prediction quality for isothermal curing was dramatically impaired. Hence, the deviation of the model towards the end of the non-isothermal curing is a result of the simultaneous fitting to isothermal, non-isothermal and cyclic DSC runs. The flattening of the curve due to temporary vitrification in case of the lowest heating rate is correctly predicted by the model.

**Figure 3.22:** Comparison of predicted and experimentally determined evolution of cure degree of SCER1 resin during non-isothermal cure

Figure 3.23 shows a comparison of isothermal DSC data and model predictions. Good agreement between experiments and numerical predictions is found, emphasizing the suitability of the Grindling kinetic model for prediction of isothermal curing processes undergoing premature vitrification.

During evaluation of isothermal DSC in Section 3.2.1, a decrease in maximum achievable cure degree when increasing the isothermal cure temperature, demonstrated in Figure 3.6, was observed. This nonphysical problem is solved during model fitting by introducing additional model parameters that represent specimen preparation times of the individual measurements. As a result, the initial cure degree of each isothermal measurement is predicted based on the fitted kinetic model and the corresponding preparation time. This gives the fitting algorithm increased freedom to match experiment and simulation and to circumvent the problem of unknown amounts of uncaptured cross-linking during measurement preparation. However, these additional model parameters basically enable the algorithm to arbitrarily shift isothermal measurement data since the initial cure degree is calculated based on the fitted model. This will most likely yield wrong results as the primary aim still is to accurately predict the experimentally observed cure degree evolution and most importantly, the maximum achievable cure degree of a chosen cure temperature. It is therefore important to also include cyclic DSC data into the fitting process since these measurements contain that information. Moreover, the method used to evaluate the cure degree of cyclic measurements differs from that applied to conventional isothermal measurements in the important fact that the cure degree is calculated based on the remaining reaction enthalpy. Hence, these measurements are not affected by uncaptured cross-linking prior to data acquisition and thus represent a suitable counterpole to the variability introduced by preparation time parameters.

**Figure 3.23:** Comparison of predicted and experimentally determined evolution of cure degree of SCER1 resin during non-isothermal cure

### **3.5.2 FCPUR polyurethane resin**

Mathematical modeling of T<sup>g</sup> and reaction kinetics of resin FCPUR is carried out using the same models as are used for resin SCER1. However, influence of vitrification needs to be accounted for using a different strategy since no isothermal DSC measurements were feasible due to the high reactivity of the resin. Hence, DSC results are complemented by T<sup>g</sup> values of the fully cured state derived from TMA measurements as depicted in Figure 3.16. The same procedure is applied to a specimen which was cured for five minutes at an isothermal temperature of 90 ◦C. From this measurement, a T<sup>g</sup> of 90.58 ± 0.24 ◦C is derived, which is about 5 ◦C below Tg,∞. This additional data point is vital for proper modeling of the reaction kinetics as it is directly related to the isothermal cure cycle of the manufacturing process. Otherwise, the model would tend to predict 100 % cure also for rather low isothermal curing temperatures due to the fact that the remaining DSC data consists of solely non-isothermal measurements [44].


**Table 3.9:** Parameters used for modeling the cure dependency of <sup>T</sup>g of resin FCPUR based on the DiBenedetto model

Figure 3.24 shows the fitted Tg-model and Table 3.9 lists the identified model parameters. As discussed in Section 3.2.2, the value of Tg,<sup>0</sup> is unknown and therefore is included into the parameter identification process as a variable. The agreement between experimentally measured values of T<sup>g</sup> and the fitted model is good, which is most likely due to the flexibility introduced by Tg,0. The figure also contains the vitrification affected cyclic measurement with a holdtime of 3600 s. Similar to resin SCER1, the cure degree of this data point is reversely calculated from the fitted Tg-model which enables to use it as an additional data point in the subsequent fitting of the kinetic model.

**Figure 3.24:** Comparison of fitted <sup>T</sup>g-model and experimentally measured values of resin FCPUR. Values gathered from DSC measurements showing influence of physical aging are marked by filled symbols.


**Table 3.10:** Parameters used for modeling the cure degree evolution of resin FCPUR based on the Grindling kinetic model

The identified parameter set used for predicting the cure evolution of FCPUR based on the Grindling kinetic model is listed in Table 3.10. Cure predictions derived from this model are shown and compared to experimental measurements in Figure 3.25. The accuracy of the model below roughly 55 % of cure is poor, especially in case of the slowest heating rate of 3 K min−1. This deviation may be caused by the change in reaction mechanism as investigated at the beginning of Section 3.2.2. However, this inaccuracy is uncritical since it only affects the cure evolution in the range before gelation. Consequently, the impact on predicted process-induced strain and stress is negligible.

**Figure 3.25:** Comparison of predicted and experimentally determined evolution of cure degree of FCPUR resin during non-isothermal cure

### **3.5.3 FCER3 epoxy resin**

No DSC measurements were available for fitting of resin FCER3. Instead, a set of model parameters was provided for use within the associated research project. The corresponding parameters are listed in Table 3.11. However, details of the experimental characterization technique as well as raw data is unknown to the author.


**Table 3.11:** Parameters used for modeling the evolution of cure degree and <sup>T</sup>g of resin FCER3 based on the models of Grindling and DiBenedetto

## **3.6 Discussion and concluding remarks**

In this chapter, reaction kinetics of the resins used in this thesis has been analyzed and modeled using suitable mathematical models for the evolution of the cure degree and glass transition temperature. Characterization of the cure rate has been carried out using conventional DSC technique. While non-isothermal measurements give reliable results, isothermal measurements are found to be prone to errors as premature curing takes place before the first data point has been acquired. This is especially problematic in case of fast-curing resins due to their high reactivity. In this study, this has been compensated for by making the initial cure degree an additional parameter in model fitting. However, it would be favorable to avoid this source of error in advance by using alternative measurement techniques. An important obstacle in this regard is thermal inertia of the specimen and the DSC cell, which limits the maximum applicable heating rate and prolongs temperature equalization prior to the measurement. This problem has been solved by a variant of DSC called Flash-DSC. Here, the thermal inertia is minimized by reducing the required specimen mass to a minimum, and by maximizing the area of contact between sensor and specimen. This enables much higher temperature rates and may also allow to shorten the duration of temperature equalization.

During cyclic DSC runs, physical aging is observed in case of long dwell times at isothermal conditions. As this phenomenon is a potential source of error, it would be favorable to characterize reaction kinetics without monitoring thermal properties like heat flow from reaction. Hence, a promising alternative to DSC may be in-situ cure characterization using non-destructive NIR spectroscopy. Erdmann et al. [69] compared NIR and DSC and found good agreement as well as no influence of physical aging when investigating pre-cured specimens. A similar approach is presented by Duemichen et al. [70], who also found good agreement between NIR and DSC based model fitting and subsequent cure predictions. Due to the fact that NIR measurements do not rely on thermal effects like DSC does, NIR technique might be especially useful for characterization of isothermal curing of fast curing resins, which challenges DSC regarding uncaptured cure before the measurement starts. However, glass transition and more importantly the simultaneous measurement of cure degree and T<sup>g</sup> is not possible when using solely NIR.

# **4 Chemo-rheology of Resins**

The viscosity of a resin is an important parameter for mold-filling simulations. Since the material undergoes substantial morphological changes during crosslinking, the viscosity strongly depends on cure degree. This is commonly referred to as the chemo-rheology of a resin. Additionally, the temperature affects both viscosity and cure degree evolution.

In this chapter, the chemo-rheological behavior of resins used in RTM processes is introduced and characterization as well as mathematical modeling of this phenomenon is shown. Emphasis is put on the accurate characterization of viscosity under consideration of the special requirements of fast-curing resins.

# **4.1 Introduction and review of related work**

Many studies exist that investigated the chemo-rheology of resins and a comprehensive overview is given by Halley and Mackay [29]. In order to accurately describe the chemo-rheological behavior of a resin, a suitable combination of models for both reaction kinetics as well as viscosity is needed. The former is discussed in Chapter 3 and predicts the evolution of the cross-linking process represented by the cure degree, whereas the latter takes the cure degree as an input parameter for prediction of the viscosity.

The viscosity of a thermoset resin is closely related to the growth of molecular chains and cross-linking. If the resin is exposed to an isothermal temperature high enough to initiate curing, the viscosity increases proportional to the increase in cure degree and shows a dramatic increase as it approaches gelation, after which it is unable to flow. Under non-isothermal conditions, the viscosity evolution exhibits an initial decrease due to the rising temperature and afterwards rises even faster since the cross-linking reaction is accelerated due to increased temperature. This behavior is schematically shown in Figure 4.1. For LCM processing, a combination of isothermal and non-isothermal temperature conditions acts. At the beginning of injection in an RTM process, the resin is much cooler than the mold. Hence, it is heated until temperature equalization is established. Afterwards, the resin is cured under isothermal conditions.

**Figure 4.1:** Evolution of viscosity for isothermal and non-isothermal temperature profiles

The viscosity is one of the key parameters that define the resistance that the resin flow must overcome in order to impregnate the initially dry fiber textiles and to completely fill the mold. An increased resistance can be compensated by using a higher injection pressure but at some point, the force, which the resin flow causes on the fibers, will lead to fiber washing. Furthermore, fiber tow impregnation quality will also suffer, if the viscosity is too high. Therefore, the viscosity needs to stay within specific limits to enable manufacturing of components of good quality. Consequently, the viscosity range a chemo-rheological model must accurately describe is also limited, simplifying characterization as well as identification of model parameters. A viscosity value of about 200 mPa s is often given as a rough guiding value for the upper limit.

#### Experimental characterization

Experimental viscosity characterization can be done using different techniques of which only few are compatible with cross-linking materials. The major requirement that must be met is interchangeable sample holders as otherwise the analyzer is damaged since the cured material adheres to the sample holders. Therefore, parallel-plate rheometry is often applied in case of thermoset resins [48, 61, 64]. The principle of parallel-plate rheometry is shown in Figure 4.2. The resin is loaded in the gap hgap of two plates of similar radius r, mounted in parallel. By applying a constant angular velocity Ω to one of the plates, a shear stress τ is induced. From this, the viscosity can then be calculated by

$$
\eta = \frac{\tau}{\dot{\gamma}} = \frac{K\_{\tau} \cdot M}{K\_{\gamma} \cdot \Omega}. \tag{4.1}
$$

In this equation, γ˙ refers to the shear rate. The shear stress τ is derived from the torque that is required to maintain the constant angular velocity and a geometry dependent factor K<sup>τ</sup> . The strain rate, is calculated from the angular velocity Ω and another geometry dependent factor Kγ. The factors are given by

$$K\_{\tau} = \frac{2}{\pi \cdot r^3}, \qquad K\_{\gamma} = \frac{r}{h\_{\text{gap}}}.\tag{4.2}$$

**Figure 4.2:** Principle of parallel-plate rheometry

#### Numerical modeling

Chemo-rheological models that are commonly used for thermoset resins can be categorized in two major classes: WLF- and Arrhenius-based models. WLF based models assume that a change in cure degree can be expressed in terms of a change in temperature. As an example, the viscosity model introduced by Hesekamp et al. [71] is given by

$$\ln\left(\frac{\eta\left(T,\alpha\right)}{\eta\left(T\_{\text{ref}},\alpha=0\right)}\right) = \frac{b\_1\cdot\left(T-T\_{\text{g},0}\right)}{b\_2+T-T\_{\text{g},0}} - \frac{b\_1\cdot\left(T-T\_{\text{g}}\left(\alpha\right)\right)}{b\_2+T-T\_{\text{g}}\left(\alpha\right)} - \frac{c\_1\cdot\left(T-T\_{\text{ref}}\right)}{c\_2+T-T\_{\text{ref}}}.\tag{4.3}$$

The first two terms on the right hand side of this equation model the influence of cross-linking on the viscosity, for which T<sup>g</sup> is used as reference. The last term is responsible for effects from changes in temperature and utilizes Tref as reference. The model can be fitted to experimental results using the four model parameters b1, b2, c<sup>1</sup> and c2. Other variants of WLF-based models have been used by Karkanas et al. [72] and Mijović and Lee [73].

An often applied Arrhenius based model was introduced by Castro and Macosko [74] and is used in several studies [48, 61, 64, 75, 76]. In this model, the resin viscosity is given by

$$
\eta = \eta\_0 \cdot \left(\frac{\alpha\_{\text{gel}}}{\alpha\_{\text{gel}} - \alpha}\right)^{C\_1 + C\_2 \cdot \alpha} \tag{4.4}
$$

with the initial viscosity of the uncured resin η<sup>0</sup> defined by

$$
\eta\_0 = A\_\eta \cdot \exp\left(\frac{E\_\eta}{\mathbf{R} \cdot T}\right). \tag{4.5}
$$

#### Fast-curing resins

When using fast-curing resins, characterization of the chemo-rheology in principle is affected by the same problems that were outlined in Section 3.1 for characterization of the reaction kinetics. The difference can be found in the time range that is most relevant for accurate simulation of the process. Since most of the mold filling process takes place during the very beginning of the curing process, measurement of the viscosity during this early time range is very important. Contrary to this, the commonly applied preparation process for rheometer measurements involves many manual and time-consuming steps (cf. Figure 4.3). Similar to DSC measurements, the final steps of loading of the sample into the preheated sample holders and gap adjustment take place at elevated temperatures. This is similar to what the resin is exposed to in the first seconds after it entered the hot mold. Having as much as information on the viscosity during this time span is therefore crucial for accurate mold-filling simulations. Moreover, in contrast to characterization of reaction kinetics, the initial value of viscosity cannot be assumed to be zero. It is therefore vital to start the measurement with the least possible delay to get data points as close to the initial viscosity as possible. This way, the error which is introduced by the fitted viscosity model when used in the time region before the first data point can be greatly reduced.

**Figure 4.3:** Comparison of measurement time and process time

# **4.2 Experimental characterization of viscosity**

The viscosity evolution of two resins is investigated in this thesis: SCER1 and FCER3. The latter is used for mold filling simulations in Section 5. Both resins were characterized using a TA Instruments ARES-G2 parallel-plate rheometer. Corresponding measurements were carried out at the Institute for Chemical Technology and Polymer Chemistry (ITCP) at the Karlsruhe Institute of Technology (KIT) in Karlsruhe, Germany.

Since the FCER3 resin is a fast-curing resin, a method has been developed which significantly shortens the time needed for specimen preparation [48, 77]. The manual preparation process can be divided in several steps, as is shown in Figure 4.4.

**Figure 4.4:** Comparison of manual and automated specimen preparation process [48]

Firstly, resin and hardener are dosed with respect to the required ratio. Afterwards, both components are mixed until a homogeneous liquid is obtained. From this, a small amount is put into the preheated sample holders of the rheometer. From this point on, curing starts. However, before the actual measurement can be initiated, the parallel-plate setup needs to be closed to the desired gap. During this step, it must be ensured that the gap is completely filled with the sample as otherwise erroneous values will be measured. Since the sample mass of the manual process varies from sample to sample, this step holds great uncertainty. Furthermore, approaching the gap, which corresponds to the sample mass, is time-consuming and prevents measurement of the viscosity at the very beginning of the curing process. Geissberger et al. [78] reported a time span of 30 s to 60 s for this process to complete. To overcome this problem, the preparation process is automated using a lab-scale mixing and dosing unit. This way, the sample mass can be directly injected in the gap of the parallel-plate setup, as is shown in Figure 4.5. Furthermore, this setup enables reproducible dosing of the mixed resin sample by prescribing a corresponding dosing time. Hence, the specimen size can be calculated beforehand and the automated closure feature offered by the rheometer can be used to rapidly approach the final measurement gap. It is even possible to directly inject into the measurement gap and to avoid any gap adjustments. This, however, depends on the viscosity of the uncured resin since the maximum injection pressure of the dosing and mixing unit is limited. The diameter of the needle must be chosen accordingly and therefore might be bigger than the desired measurement gap.

**Figure 4.5:** Direct injection of the sample into the parallel-plate setup

The measurements presented in this thesis were conducted using a loading gap of 1.5 mm and a measurement gap of 1 mm. Furthermore, aluminum plates with a diameter of 25 mm were used and the automated specimen preparation process has been setup such that it quickly injects the correct amount of resin. An angled needle with an inner diameter of 0.84 mm, an outer diameter of 1.28 mm and a length of 12.7 mm was used to guide the resin to the center of the preheated metal plates. After injection, the measurement procedure is initiated and the rheometer starts with closing the gap until the specified measurement gap has been reached. For this progress, a closing speed of 5 mm s−1 was chosen which results in a theoretical closing time of only 0.1 s. In reality, this process took about 1 s since the rheometer needs some time for acceleration and deceleration. In total, a time of about 5 s to 7 s elapsed between initiating the measurement procedure and actual data acquisition. For the latter, a rotational frequency of 1 Hz was programmed and the measurement was aborted as soon as the viscosity exceeded a value of 2 Pa s.

## **4.2.1 SCER1 epoxy resin**

The viscosity of SCER1 resin is characterized using three different temperatures. Figure 4.6 shows the viscosity increase during isothermal curing at 90 ◦C, 100 ◦C and 120 ◦C. Each curve represents mean values of three individual measurements. The scatter at the beginning of the measurement arises from temperature inhomogeneities. Since the resin is not preheated prior to injection between the sample holders, temperature equalization takes place at the very beginning of the measurement. However, the amount of scatter decreases quickly after start of the measurement as equalization is facilitated by the beginning rotation.

As was discussed at the end of Section 4.1, the mold-filling stage of a typical, commercially used RTM process needs to be completed after about 30 s. This time frame is highlighted in Figure 4.6. By evaluating the results at the beginning and end of this time range, it can be seen that at the beginning, temperature is the most important influencing factor. Towards the end, however, the mechanism changes. At this point in time, cross-linking has already advanced significantly in case of the highest temperature of 120 ◦C, causing a steep increase in viscosity. Hence, after 30 s, cross-linking is the major influencing factor. This emphasizes the importance of capturing as much viscosity data as possible during this time frame in order to increase the accuracy of simulations.

**Figure 4.6:** Viscosity evolution of SCER1 resin for isothermal temperatures of 90 ◦C, 100 ◦C and 120 ◦C

### **4.2.2 FCER3 epoxy resin**

The viscosity of the fast-curing FCER3 resin is investigated for a single isothermal temperature of 110 ◦C, which is applied in the real process and therefore also used for corresponding mold-filling simulations of Chapter 5.

Figure 4.7 shows the evolution of the viscosity during isothermal curing at 110 ◦C. The presented curves represents three individual measurements, which are mean averaged for subsequent model fitting. In comparison to SCER1, FCER3 resin shows a much faster increase in viscosity. Consequently, a greater portion of the viscosity development already happens during the first 30 s. Proper measurement and characterization of this time period is therefore even more important in case of fast-curing resins.

## **4.3 Mathematical modeling of viscosity**

Based on the introduction of viscosity models given in Section 4.1, the Castro-Macosko model is chosen for modeling the viscosity of SCER1 and FCER3 resin. In both cases, model parameters are fitted using a Python script as well as Differential Evolution [79] as supplied by the Python package SciPy [80]. According to Equations 4.4-4.5, parameters which need to be identified are C1, C2, A<sup>η</sup> and Eη. In order to improve prediction quality in the time range that is most important for RTM processes, model fitting is constrained to viscosity values below 1 Pa s.

### **4.3.1 SCER1 epoxy resin**

Figure 4.8 shows a comparison of measured and predicted viscosity evolution of SCER1 resin. Corresponding model parameters are given in Table 4.1. Good agreement of experiment and model is achieved for each isothermal temperature using the Castro-Macosko model.

**Figure 4.8:** Comparison of model prediction and measured viscosity evolution of SCER1 resin for isothermal temperatures of 90 ◦C, 100 ◦C and 120 ◦C

**Table 4.1:** Identified model parameters for SCER1 resin and Castro-Macosko viscosity model


### **4.3.2 FCER3 epoxy resin**

Figure 4.9 shows a comparison of measured and predicted viscosity evolution of FCER3 resin. Since this resin is a fast-curing one, it makes a huge difference whether fitting is carried out using data points from the very beginning of the reaction or after a delay of 30 s which is typical for fully manual specimen preparation processes. In order to show the difference, two scenarios are investigated. In the "automated" scenario, the model fitting process is conducted using all available data points. Within the "manual" scenario, however, a limited set of data points starting from 30 s upwards is used during model fitting. This way, the fully manual preparation process is simulated. The model parameters of both scenarios are given in Table 4.2.

**Figure 4.9:** Comparison of isothermal model predictions based on manual and automated specimen preparation process and measured viscosity evolution of FCER3 resin

Good agreement of experiment and model is achieved for both scenarios. However, due to the lack of data points prior to 30 s, the fitted model of the "manual"


**Table 4.2:** Identified model parameters for FCER3 resin and Castro-Macosko viscosity model

\* Only viscosity data from 30 s upwards is used

scenario underestimates the viscosity in this time range. In case of the "automated" scenario, the viscosity prediction gives accurate results from about 8 s upwards. Consequently, the time range during which no information on viscosity is available is much shorter compared to the "manual" scenario. Moreover, when using the model of the "manual" scenario for simulation of a resin injection process of about 30 s, the viscosity is significantly under-predicted during the whole mold-filling process. This has huge impact on the resulting injection pressure as is demonstrated in Section 5. A good indication of the error is given by the initial viscosity η<sup>0</sup> = η (t = 0), which is extrapolated using the fitted models. The "automated" scenario results in an initial viscosity of η<sup>0</sup> = 7.14 mPa s whereas a much lower value of η<sup>0</sup> = 3.10 mPa s is calculated from the model of the "manual" scenario.

### **4.4 Discussion and concluding remarks**

Viscosity characterization of fast-curing resins challenges established measurement methods. The problem lies in the premature curing of the resin prior to the beginning of the actual data acquisition. This is on the one hand a limitation which arises from the specimen preparation process, more specifically insufficient precision during sample loading. On the other hand, the measurement equipment, e.g. a rheometer, needs some time for adjustments during which the resin is already in contact with the preheated sample holders. As a result, the start of the measurement does not coincide with the initiation of the curing process. Consequently, a time window exists, during which no viscosity data is gathered. However, during mold-filling, this time window has a huge impact on injection pressure prediction as is shown in Chapter 5. The unknown viscosity evolution at the beginning of the process needs to be extrapolated by a suitable viscosity model. This introduces an inaccuracy in viscosity prediction that depends on the period of time for which no experimental data is available. The presented automated specimen preparation process is able to minimize this uncertainty, leading to more accurate model predictions. However, there is still room for improvement. The needle diameter limits the measurement gap in such a way that the sample needs to be injected while the parallel plate setup is set to a larger loading gap. Consequently, after initiation of the measurement, the rheometer first needs to move the sample holders in order to reach the chosen measurement gap. This is time-consuming and prolongs the time window of unrecorded viscosity evolution. Since the viscosity of resin and hardener is temperature dependent, a possible solution for this problem can be achieved by heating of the resin components. This way, the viscosity is lowered and the same applies to the pressure within the dosing and mixing unit. This enables the use of needles with smaller diameters and subsequently an injection right into the measurement gap rather than a loading gap.

# **5 Mold Filling Simulation**

The economic efficiency of the RTM process is dominated by the cycle time which to a large part is defined by the mold filling step. It is therefore crucial to choose optimal process parameters which help to reduce the cycle time by minimizing the filling time and maximizing the cure degree at the end of the mold filling in order to shorten the cure time needed to reach a sufficiently high cure degree and Tg. The dependency of the resin viscosity on the cure degree and the associated risk of incomplete filling due to too high viscosity further complicates identification of optimal process parameters. Therefore, numerical simulation is needed for prediction of the mold-filling process, which is described in this chapter.

## **5.1 Introduction and review of related work**

The RTM process has been extensively investigated in various studies both experimentally and numerically [7, 81, 82]. The physics of flow in porous media is often described using Darcy's law [83], which is given by

$$
\tilde{v} = -\frac{k}{\eta} \nabla p.\tag{5.1}
$$

This equation relates the filter velocity or Darcy velocity v˜ with the pressure gradient ∇p. The flow resistance that acts against this pressure driven flow arises from the fiber textile permeability k and the resin's viscosity η, which

therefore dominate the filling process. The relation between the actual fluid velocity v and the filter velocity v˜ is given by

$$v = \frac{\tilde{v}}{\Psi} = -\frac{1}{\Psi} \frac{k}{\eta} \nabla p = -\frac{1}{1 - \varphi\_{\rm f}} \frac{k}{\eta} \nabla p,\tag{5.2}$$

in which Ψ and ϕ<sup>f</sup> denote porosity and fiber volume fraction, respectively. The difference between these two velocity quantities is caused by the fact that in a porous media, only a fraction of the total volume is available for fluid flow.

The permeability tensor *k* is defined by the properties of the preform and is anisotropic due to the presence of the fibers. During draping, the permeability is altered through changes in fiber orientation and fiber volume fraction [8, 84]. Consequently, these draping effects greatly affect the injection pressure in experiments during which a constant mass flow rate is applied, as reported by Bickerton et al. [85]. In another paper of Bickerton et al. [86] an extensive study is conducted on the effect of mold corners with different radii on injection pressure and flow front progression. They found a strong increase in injection pressure for decreasing radii and increasing fiber volume fraction. However, flow front progression was not significantly altered by these properties. In Bickerton et al. [87], these experiments were further investigated using mold-filling simulations. Compared to the influence from changing the radii of the mold, a much greater influence of non-constant cavity thickness and thus permeability on injection pressure was found. Louis and Huber [84] investigated the draping of a double curved structure and its influence on permeability due to shearing of the fabric. Allthough the permeability was significantly affected by the shear induced increase in fiber volume fraction, the numerical predictions using the results of a draping analysis showed excellent agreement with experiments. Obviously, the accuracy of mold filling simulations does strongly depend on the precise knowledge of the local fiber structure, and thus local permeability. This information can be supplied by draping simulations, which can be a demanding task in case of complex geometry as common in practice. Furthermore, derivation of local fiber volume fraction from draping simulations is often impossible due to the applied model approaches. Therefore, the simulation needs to rely on estimated values or time consuming experimental measurements of this quantity at several positions of the composite structure. The latter is presented by Magagnato and Henning [8] and led to an improved accuracy of mold-filling simulation of a complex shaped component. Allthough not feasible for large structures, analyzing the local fiber structure using Computed Tomography (CT) scans can help to gain a better understanding of localized effects [10]. Another textile related influencing effect is fiber bed compaction which in case of RTM occurs towards the end of closure of the mold. In general, this effect depends on the textile's structure as well as fiber volume fraction and strain rate [88]. Closely related to this is fiber bridging and wrinkling in convex and concave sections, that can have significant effect on mold-filling due to the local formation of areas of high permeability, commonly referred to as runners [86, 89].

#### Influence of cross-linking on mold-filling

Many studies on influencing factors on injection pressure and flow front progression during RTM mold-filling concentrate on textile properties and mold design aspects, as discussed above. The influence of cross-linking of the resin and the accompanied increase in viscosity is often excluded by using substitute fluids that exhibit a constant viscosity which is similar to that of resins in uncured state at process temperature. In a real application scenario, such a low and constant viscosity during mold-filling can only be achieved by using slow curing resins. In high-volume production, however, fast-curing resins are mandatory since they enable much shorter cycle times and therefore higher cost efficiency. From this perspective, the increase of resin viscosity cannot be neglected anymore and needs to be taken into account in mold-filling simulations. In an early study, Castro and Macosko [74] investigated mold-filling and curing during Reaction Injection Molding (RIM) process. Their measurements show a sharp increase in injection pressure as soon as cross-linking starts. Consequently, pressure predictions using non-constant and constant viscosity start to diverge at this point, leading to significantly different pressure levels. Equal behavior is reported for the flow front progression in RTM mold-filling using a constant injection pressure by Deléglise et al. [21]. The prediction quality of the simulation conducted in both studies was significantly improved by using an appropriate viscosity model. Due to the introduced time and temperature dependency of the viscosity, the spatial distribution of the viscosity becomes inhomogeneous. In combination with local fiber textile properties variations, this can lead to mold filling paths that depend on time and temperature as well.

Another influence from cross-linking is overheating due to the exothermic reaction of the resin. This can lead to a thermal degradation of the composite and typically occurs in the center of rather thick laminate stacks due to the low thermal conductivity of the composite in through-thickness direction. However, as is shown by Sorrentino et al. [90], this effect is dramatically reduced if the laminate thickness decreases. With respect to the thin laminates investigated in this study, it is assumed that this effect can be neglected.

#### Numerical approaches used for mold-filling simulations

A comprehensive introduction into numerical methods for fluid flow problems is given by Ferziger and Perić [91]. Many studies, however, rely on a Finite Element (FE) based approach, often used in combination with a Control Volume (CV) approach [7, 8, 21, 76, 85, 92–94]. The latter is needed since the FE method itself is not strictly mass conservative, especially if multiple inlet and outlet regions are present in the model, which is typical for mold-filling problems. An alternative approach is given by the Finite Volume (FV) method. Since this method inherently assures mass conservation, it is widely used for Computational Fluid Dynamics (CFD). However, only limited results are available from literature that investigate the applicability of the FV method to RTM mold-filling simulations. Magagnato et al. [10] presented a study on the influence of the introduction of metallic inserts as load transmission elements into the laminate on mold-filling behavior. An emphasis was put on the formation of dry spots in the vicinity of the inserts. For this purpose, a FV approach based on the open-source CFD software OpenFOAM [95] was applied. Tracking of the individual air and resin phases is achieved using a Volume of Fluid (VoF) approach. Furthermore, they compared their predictions to FE/CV based simulations and validated both using experimental results. From this comparison, the FV based approach shows higher accuracy in prediction of both, flow front progression as well as formation and convection of entrapped air. Furthermore, due to the fact that OpenFOAM is available and licensed as open source, extension of the solvers as well as tailoring them to the special needs of process simulation is possible. Recently, Seuffert et al. [96] adopted an OpenFOAM based solver for mold-filling prediction of a Compression Resin Transfer Molding (CRTM) process. Due to the possibility of adopting solvers and involved models, the OpenFOAM framework [95] was chosen for mold-filling analysis in this thesis as well.

# **5.2 Numerical method used for mold filling simulation**

The numerical solver used in this thesis is based on the OpenFOAM solver *compressibleInterFoam*. This solver allows to calculate non-isothermal compressible two-phase flow of two immiscible fluids. The use of a compressible approach is motivated by the fact that the air phase needs to be treated compressible in order to correctly capture the behavior of entrapped air during the course of the mold filling process. For this thesis, the *compressibleInterFoam* solver was extended in order to take flow resistance due to the fiber textiles as well as the chemo-rheology of the resin into account. The applicability of the developed numerical method is demonstrated in Magagnato et al. [10].

OpenFOAM is based on the FVM approach and therefore involves decomposition of the computational domain into many control volumes. For each control volume, a set of governing equations including the Navier-Stokes-Equations is then solved. These equations are introduced in the following.

#### Continuity equation

The conservation of mass is ensured by the continuity equation, which is given in differential form in case of a compressible fluid by

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0. \tag{5.3}$$

In this equation, <sup>ρ</sup> and *v* correspond to the fluid density and the flow velocity, respectively. The first term on the left-hand side describes the change in density over time within a control volume. The second term represents the mass flow through the boundary faces of the control volume.

#### Momentum equation

For a compressible fluid, conservation of momentum is achieved by using the following equation:

$$\frac{\partial \rho v}{\partial t} + \nabla \cdot (\rho v \mathbf{v}) = -\nabla p + \nabla \tau + \rho \mathbf{g} + \mathbf{f} \tag{5.4}$$

Here, <sup>p</sup> and *τ* denote pressure and viscous stress, respectively. The first term on the left-hand side of this equation describes the change in momentum in a control volume over time while the second term refers to the flow of momentum in and out of the control volume. The first two terms on the right-hand side introduce normal and shear forces that act on the control volume and are driven by fluid pressure and friction. The third term accounts for the effect of gravity on the flow and is often negligible. The source term *f* can be used to take other sources of momentum into account. In this case, this quantity will be used to include the resistance due to flow through a porous media as shown below. For a Newtonian fluid, the viscous stress tensor *τ* is given by

$$
\tau = \eta \left( 2D - \frac{2}{3} \delta \nabla \cdot \mathbf{v} \right). \tag{5.5}
$$

This equation relates the deformation that acts on the fluid with the resulting viscous stress. Hence, the dominating material property which links these quantities is the fluid viscosity <sup>η</sup>. Moreover, <sup>δ</sup> is the Kronecker delta and *D* represents the deformation rate tensor which is given by

$$D = \frac{1}{2} \left( \nabla \boldsymbol{v} + \left( \nabla \boldsymbol{v} \right)^{T} \right). \tag{5.6}$$

#### Energy equation

The total energy E of a system consists of thermal and mechanical contributions and its conservation is given by

$$\frac{\partial \rho E}{\partial t} + \nabla \cdot (\rho v E) = -\nabla \cdot \mathbf{q} + S\_{\mathbf{Q}} + \nabla \cdot (\sigma \cdot \mathbf{v}) + \rho \mathbf{g} \cdot \mathbf{v}. \tag{5.7}$$

On the left-hand side of this equation, the change in total energy over time is captured by the first term while the second term describes the transfer of total energy in and out of the control volume. The first two terms on the right-hand side represent thermal sources of energy. The remaining terms incorporate energy due to mechanical processes like strain energy. Moreover, *<sup>q</sup>*, <sup>S</sup><sup>Q</sup> and *<sup>σ</sup>* denote heat flux due to diffusion, volumetric heat sources and mechanical stress tensor, respectively. The latter can be decomposed into the thermodynamic pressure p and the viscous stress tensor *τ* , given by

$$
\sigma = \tau - pI \tag{5.8}
$$

in which *I* represents the identity tensor. Moreover, the total Energy <sup>E</sup> can be decomposed into a thermal and mechanical part by using the relation

$$E = e + K.\tag{5.9}$$

Here, the thermal portion is given by the specific internal energy e while K refers to the specific kinetic energy. The latter is calculated from the fluid velocity by using the relation

$$K = \frac{1}{2}|\boldsymbol{v}|^2. \tag{5.10}$$

By inserting Equations 5.8-5.10 into Equation 5.7, an alternative formulation of the energy conservation is found

$$\frac{\partial \rho e}{\partial t} + \nabla \cdot (\rho v e) + \frac{1}{2} \frac{\partial \rho |v|^2}{\partial t} + \nabla \cdot \left(\frac{1}{2} \rho |v|^2 v\right) + \nabla \cdot (v p)$$

$$= -\nabla \cdot (\alpha\_{\text{eff}} \nabla e) + S\_{\text{Q}} + \nabla \cdot (\tau \cdot v) + \rho \mathbf{g} \cdot v. \tag{5.11}$$

In this equation, the heat flux *q* has been replaced by the relation

$$
\hat{q} = -\alpha\_{\text{eff}} \nabla e,\tag{5.12}
$$

in which αeff denotes the effective thermal diffusivity of the mixture.

Equation 5.11 can be further modified into a temperature based formulation by assuming perfect gas behavior of the fluid. In this case, the specific heat at constant volume c<sup>v</sup> and the specific heat at constant pressure c<sup>p</sup> can be assumed to be independent of temperature and pressure. Hence, the specific internal energy can be calculated by

$$e = c\_\mathbf{v} T,\tag{5.13}$$

which results in a temperature based form of the energy equation:

$$\frac{\partial \rho c\_{\mathrm{v}}T}{\partial t} + \nabla \cdot (\rho v c\_{\mathrm{v}} T) + \frac{1}{2} \frac{\partial \rho |v|^2}{\partial t} + \nabla \cdot \left(\frac{1}{2} \rho |v|^2 v\right) + \nabla \cdot (v p)$$

$$= -\nabla \cdot \left(\alpha\_{\mathrm{eff}} c\_{\mathrm{v}} \nabla T\right) + S\_{\mathrm{Q}} + \nabla \cdot (\boldsymbol{\tau} \cdot \boldsymbol{v}) + \rho \mathbf{g} \cdot \boldsymbol{v}.\tag{5.14}$$

Implementations of this equation in OpenFOAM solvers usually neglect contributions from the mechanical sources ∇ · (*τ* · *v*) and <sup>ρ</sup>**g** · *v*.

#### Interface capturing method

The multiphase solvers provided by OpenFOAM use a Volume of Fluid (VoF) based approach for capturing and maintaining a sharp interface between resin and air phase. This technique was introduced by Hirt and Nichols [97] and is explained in more detail in Ubbink [98] and Albadawi et al. [99].

In a VoF approach, a scalar phase fraction field Θ (*x*, t) is introduced, in which cells completely filled with air have a value of Θ=0 while those completely filled with fluid (resin) possess a value of Θ=1:

$$\Theta\left(x,t\right) = \begin{cases} 1, & \text{element contains only resin} \\ 0, & \text{element contains only air} \end{cases} \tag{5.15}$$

However, from a numerical point of view, it is difficult to maintain a sharp interface between resin and air phase. Due to numerical diffusion, elements near to the current position of the flow front will also show values in the range between 0 and 1. Since this smearing effect will increase with each time step, interface compression techniques have been developed. In OpenFOAM solvers, this is achieved using a technique called *Multidimensional Universal Limiter with Explicit Solution* (MULES). A detailed explanation of this approach is given in Deshpande et al. [100].

#### Equation of state

The above introduced set of governing equations for continuity, momentum and energy still represents an underdetermined system of equations. Hence, an additional equation of state is needed that relates thermodynamic properties of the involved fluids.

The equation of state for an ideal gas is given by

$$pV = n\mathbf{R}T = m\mathbf{R}\_{\mathbb{S}}T,\tag{5.16}$$

in which n, R and R<sup>S</sup> denote amount of substance in moles, gas constant and specific gas constant, respectively. The latter is defined as the ratio of the gas constant R and molar mass M:

$$\mathbf{R}\_{\mathbf{S}} = \frac{\mathbf{R}}{\mathbf{M}}.\tag{5.17}$$

Rearranging enables to derive a different form of equation of state that enables to calculate the density based on temperature and pressure:

$$
\rho\_\mathrm{a} = \frac{p}{\mathrm{R}\_\mathrm{a} T}.\tag{5.18}
$$

This relation is derived from the ideal gas law and assumes a constant heat capacity of the gas which does not change with pressure or temperature. This approach is known as *perfect gas* and is often used in computational fluid dynamics of gases. Hence, it is used in this thesis for the air phase and the index *a* is added to the variable names in order to distinguish between the properties of each phase. The specific gas constant of air is calculated as

$$\mathbf{R\_a} = \frac{\mathbf{R}}{\mathbf{M}} = \frac{8.314 \, 46 \, \text{J} \, \text{mol}^{-1} \, \text{K}^{-1}}{2.896 \, 44 \times 10^{-2} \, \text{kg} \, \text{mol}^{-1}} = 287.06 \, \text{J} \, \text{kg}^{-1} \, \text{K}^{-1}.\tag{5.19}$$

For liquids, a different approach is used which takes their weak compressibility into account:

$$
\rho\_{\mathbf{r}} = \frac{p}{\mathbf{R}\_{\mathbf{r}}T} + \rho\_{\mathbf{r},0}.\tag{5.20}
$$

This approach is known as *perfect fluid* [101] and is used for the resin phase in this thesis. Corresponding properties are therefore marked with index r. It is common to assume a high value for the specific gas constant in order to realize a weak compressibility. Here, a value of R<sup>r</sup> = 3000 J kg−1 K−1 and an initial density of ρr,<sup>0</sup> = 1150 kg m−3 are used.

Properties of the mixture are calculated using the phase fraction Θ introduced in Equation 5.15. For density and viscosity, this is given by

$$
\rho = \Theta \rho\_{\rm r} + (1 - \Theta) \,\rho\_{\rm a} \tag{5.21}
$$

and

$$
\eta = \Theta \eta\_{\rm r} \left( \alpha, T \right) + \left( 1 - \Theta \right) \eta\_{\rm a}. \tag{5.22}
$$

As indicated in this equation, the resin viscosity η<sup>r</sup> is modeled cure and temperature dependent, as introduced in Chapter 4.

#### Cure degree equation

The cure degree is a property of the moving resin and therefore needs to be transported with the resin flow. At the same time, the cure degree increases depending on the local temperature and time step size. This is achieved by using a transport equation for the cure degree which is given by

$$\frac{\partial \rho \alpha}{\partial t} + \nabla \cdot (\rho v \alpha) = \rho \dot{\alpha}. \tag{5.23}$$

The first term on the left-hand side corresponds to the change in cure degree within the control volume. The second term is responsible for transportation of the cure data along with the resin flow. The cure rate as calculated by the reaction kinetics model contributes to this equation in form of the source term on the right-hand side of the equation.

#### Resistance due to textile permeability

In order to take the flow resistance due to the fiber textile into account, Darcy's law given by Equation 5.2 can be used. The implementation in OpenFOAM, however, is based on an extension of Darcy's law that was introduced by Forchheimer [102]:

$$
\nabla p = -\frac{\eta}{k}\tilde{v} - \frac{1}{2}\rho F|\tilde{v}|\tilde{v}.\tag{5.24}
$$

The first term on the right-hand side of this equation corresponds to Darcy's law. The pressure drop calculated by Darcy's law is based solely on viscous drag. Hence, the extension introduced by Forchheimer and given by the second term on the right-hand side allows to also take a contribution from turbulence into account. This is controlled by the Forchheimer coefficient *F*. Whether this is actually required depends on flow velocity *v*, permeability *k* and resin viscosity η. As introduced by Ward [103], a modified Reynolds number can be derived from these parameters:

$$\mathrm{Re}\_{\mathbf{k}} = \frac{v\_i k\_i^{1/2}}{\nu} = \frac{v\_i k\_i^{1/2}}{\eta/\rho}. \tag{5.25}$$

Ward [103] also showed experimentally that Darcy's law is valid if Re<sup>k</sup> < 1 and that Forchheimer's law is valid for a value range of 0.1 ≤ Re<sup>k</sup> ≤ 100. Based on typical values for resin viscosity at the beginning of the process (η = 0.01 Pa s), textile permeability in fiber direction (k<sup>11</sup> = 2.72 <sup>×</sup> <sup>10</sup>−11 <sup>m</sup>2) and resin density (ρ = 1150 kg m−3), the critical value of Re<sup>k</sup> = 1 will be surpassed for a flow velocity of 1.67 m s−1. In RTM processing, much lower values are common and a contribution from turbulence can therefore be neglected by setting the Forchheimer coefficient to *F* = 0. Hence, the flow resistance induced by the fiber textile is incorporated in the momentum equation by using a source term given by

$$f = \nabla p\_{\text{DF}} = -\frac{\eta}{k\left(x\right)}\tilde{v} = -\frac{\eta}{k\left(x\right)}\Psi\left(x\right)v. \tag{5.26}$$

In this equation, permeability tensor *k* and porosity <sup>Ψ</sup> are defined as dependent on the location *x* within the computational domain. This is needed in order to take local changes in fiber orientation and fiber volume content into account. Moreover, the relation between the filter velocity *v***˜** and the fluid velocity *v*, given by Equation 5.2, is used. In this context, multiplication with the porosity ensures that the flow resistance due to the fibers is only applied to the relevant resin portion of the fluid volume.

# **5.3 Mold filling simulation of car floor structure**

In this section, the mold-filling simulation of a car floor structure is presented as a practical application example. As shown in Figure 5.1, the component is large and of complex geometry. In order to facilitate manufacturing of corresponding preforms, the component was decomposed into several subpreforms which were assembled right before resin injection. This is indicated by zones of different colors in Figures 5.2 and 5.3. The motivation for this subdivision of the structure is facilitation of the draping process and the possibility to add local reinforcements in form of unidirectional patches. The latter for example was applied to the dark purple regions in order to improve side-impact crash behavior of the floor structure. Subdivision of the part in several subpreforms, however, requires zones where neighboring subpreforms overlap each other in order to be able to transfer loads. This is shown in detail in Figure 5.3. Due to tolerances in preform dimensions as well as positioning, overlapping zones need to provide some extra space to compensate for this inaccuracy. This is needed to avoid the occurrence of high fiber volume contents and thus low permeability. As a consequence, however, this extra space leads to linear regions of low fiber volume content, known as resin runners. These can have a huge influence on the mold-filling process and therefore are taken into account in this study as presented in Figure 5.3.

The model consists of 1.86 million 3D elements and is subdivided into submeshes according to the subpreform decomposition depicted in Figures 5.2 and 5.3. This enables individual assignment of permeability properties according to the varying layup of each zone. More details regarding the composite layup are given in Appendix A in Figure A.1. Since thickness gradients are neglected in this study, the permeability properties are assumed to be homogeneously distributed within each individual zone. Furthermore, the stacks consist of variable amounts of layers of unidirectional as well as biaxial non-crimp fabric. From these stacks, homogenized permeability values are derived using the approach proposed by Bruschke and Advani [81]. The permeability of both textile types, on which this

**Figure 5.1:** Model geometry and dimensions of the car floor structure

homogenization is based, is given by Table 5.1. The used measurement setup is described in more detail in Magagnato and Henning [104].

**Table 5.1:** Properties of the carbon fiber textiles used in mold-filling simulations for a fiber volume fraction of 50 %. More details of the fiber textiles are given by Table 2.1.


The component is surrounded by a small area of high fiber volume fraction, which serves as fiber clamping. Within the numerical model, a fiber volume fraction of 90 % is defined in this area in order to take it into account during mold-filling. Runners due to sub-preform tolerances (cf. Figure 5.3) are modeled by assigning a fiber volume fraction of 30 % to the corresponding elements. For all remaining elements, a fiber volume fraction of 50 % is defined.

The boundary conditions applied during mold-filling simulation are indicated in Figure 5.2. At the inlet patch, a constant mass flow rate is prescribed whereas

**Figure 5.2:** Model sectioning and boundary conditions of car floor structure used for mold-filling analysis

**Figure 5.3:** Detailed view on an overlapping zone of two subpreforms and the thereby produced resin runners

a zero-gradient is forced for the pressure. A mass flow rate of 140 g s−1 was applied at the inlet gate in the experimental study. However, this value cannot be transferred directly to the inlet boundary condition of the simulation. The reason for this is the fact that the model introduced in Section 5.2 represents a two-fluid approach although a third phase given by the fibers exists in reality. Thus, each control volume will be completely filled by either air or resin in the simulation. In reality, however, the amount of remaining space depends on the porosity that is present in the cavity. As a consequence, the amount of resin required to fill an element in the simulation is unphysically high. Thus, when enforcing a constant mass flow rate of resin at the inlet, the flow rate needs to be scaled according to the porosity. As the mass flow rate is strongly connected to the fluid velocity, an equation similar to Equation 5.2 can be applied:

$$
\dot{m} = \frac{\dot{m}\_{\rm R}}{\Psi} = \frac{\dot{m}\_{\rm R}}{1 - \varphi\_{\rm f}}.\tag{5.27}
$$

In this equation, m˙ <sup>R</sup> refers to the resin mass flow that would be applied in an experiment. In this case, the fiber volume content is 50 % for almost the entire model. Thus, the mass flow rate needs to be doubled. However, due to the use of a symmetry plane, the specified mass flow rate acts only on half of the area of the original inlet patch (cf. Figure 5.2). This in turn leads to a halving of the flow rate. Hence, a mass flow rate of 140 g s−1 is applied to the inlet patch.

At the outlet regions, a constant pressure of 1 bar and a zero-gradient for the velocity field is assigned. With respect to computation times, only half of the symmetric part is considered during simulation by defining a symmetry plane and boundary condition. All remaining boundary faces are assigned to a fixed wall boundary patch for which a zero-gradient in pressure as well as a velocity of zero in all directions is forced. Since the process is isothermal, a constant temperature of 110 ◦C is defined throughout the computational domain. Furthermore, the resin enters the mold with a cure degree of zero. Evolution of cure and viscosity is defined by the models introduced in Chapters 3 and 4.

During viscosity model parametrization in Section 4.3.2, two scenarios have been investigated in order to show the importance of gathering viscosity data as close to the start of injection as possible. In order to assess the consequences of this fact with respect to mold-filling during an RTM process, these two scenarios, namely "automated" and "manual", are also investigated in this section. Furthermore, a scenario "constant" is added that represents mold-filling using a constant viscosity of η = 14.32 mPa s. In this study, the constant viscosity is defined by the mean value of the viscosity evolution of the scenario "automated", evaluated between 0 s to 28 s. The upper limit of this time range was chosen according to the time needed to completely fill the mold using a constant mass flow rate of 140 g s−1. The viscosity evolution of these three scenarios is shown in Figure 5.4.

**Figure 5.4:** Comparison of viscosity evolution of resin FCER3 during mold-filling for three scenarios

Figure 5.5 shows the pressure distribution after 20 s of injection using a mass flow rate of 140 g s−1. Position and shape of the flow front is comparable for the three investigated scenarios. The maximum pressure, however, differs significantly. This is most obvious when comparing the constant viscosity scenario to the scenario "manual".

Figure 5.6 shows the cure degree distribution after 28 s of injection using a mass flow rate of 140 g s−1. At this point in time, the mold is almost completely filled. Only few amounts of trapped air in the vicinity of the fiber clamping area remain. Furthermore, a cure degree difference of roughly 30 % develops between the freshly mixed resin at the inlet and the resin that has been in the mold the longest. As is shown in Figure 5.6, the latter must not inevitably be located near

b) manual specimen preparation -

c) automated specimen preparation -

**Figure 5.5:** Comparison of cavity pressure after 20 s of resin injection for the three investigated viscosity characterization scenarios

a) constant viscosity -

b) manual specimen preparation -

c) automated specimen preparation -

**Figure 5.6:** Comparison of cure degree distribution after 28 s of resin injection for the three investigated viscosity characterization scenarios

the outlet region as it would be expected for a simpler component geometry. Instead, the resin distribution shows island formation which is caused by local fiber textile properties, runners near to sub-preform overlapping zones as well as large gradients in viscosity.

The influence of viscosity gradients on the formation of areas of high cure degree and thus viscosity is evident in a comparison of the three investigated scenarios. This effect is strongest if the manual specimen preparation process is used for viscosity characterization since the viscosity in this case is underestimated by the model, as is shown in Section 4.3.2. Consequently, this scenario shows the biggest increase in viscosity during the first 30 s of the process and hence also the strongest tendency of island formation. On the contrary, in the scenario of constant viscosity, the gradient in viscosity obviously is zero and island formation is the weakest of all. Moreover, due to the absence of a gradient in viscosity, the resulting inhomogeneity in cure degree distribution can be attributed to the influence of local fiber textile properties and runners. Especially the latter act as additional line gates that provide a potential shortcut for the resin due to the much higher permeability in this areas. This is also the reason for the large blue areas of low cure degree in Figure 5.6. As soon as the resin reaches the begin of the overlapping zones located nearest to the inlet, it races along these regions of high permeability, which are indicated in red in Figure 5.2. Because of this additional but usually unintended runners, the mold-filling behavior switches from one typical for point gate injection to one that is more typical for line gate injection.

Island formation as well as the huge gradient in viscosity in case of the manual scenario is also evident from the spatial distribution of the resin's viscosity as given by Figure 5.7. This has significant influence on the evolution of the injection pressure that is needed to maintain the specified constant mass flow rate. To investigate this further, Figure 5.4 and Figure 5.8 show the evolution of resin viscosity and injection pressure over mold-filling time for the three scenarios. Although the viscosity of both non-constant viscosity scenarios surpasses that of the constant viscosity scenario after about two thirds of the mold-filling time, the injection pressure in the constant viscosity scenario remains at a much higher

b) manual specimen preparation -

c) automated specimen preparation -

**Figure 5.7:** Comparison of viscosity distribution after 28 s of resin injection for the three investigated viscosity characterization scenarios

value throughout the process. Similarly, the injection pressure of the manual scenario remains significantly below that of the automated scenario, although they exhibit roughly the same viscosity at the end of mold-filling.

**Figure 5.8:** Comparison of injection pressure evolution of resin FCER3 during mold-filling for three scenarios

## **5.4 Discussion and concluding remarks**

In this chapter, mold-filling simulations of a car floor structure is presented and discussed. The influence of material characterization is analyzed using viscosity model parameter sets which have been determined in Chapter 4 using either all available viscosity data or only data points acquired from 30 s on. This way, the influence of uncaptured viscosity information due to premature curing on flow front progression and pressure evolution is analyzed. Furthermore, a third scenario with a constant viscosity has been added, which represents an often applied approach in which chemo-rheology is completely neglected and a mean viscosity value is used instead.

Mold-filling simulations were carried out for all three scenarios. If the model fitting process lacks viscosity data prior to 30 s (scenario "manual"), the initial viscosity is significantly underestimated. This leads to an increased gradient in viscosity when using this parameter set for mold-filling simulation and, thus, affects flow front progression. However, the increased difference between the resin viscosity at the beginning and end of the mold-filling stage also affects the pressure evolution. The viscosity of the resin as predicted by the scenario "manual" is consistently lower than the prediction of scenario "automated". Hence, the injection pressure which is needed to maintain a constant mass flow rate is also lower in case of scenario "manual". Contrary to this underprediction of the pressure, using a constant viscosity leads to an overprediction. Hence, when using mold-filling simulation for prediction of flow front progression of variable gap processes, accurate viscosity data is vital and should be gathered from the very beginning of the reaction on. Otherwise, pressure calculation is inaccurate throughout the whole course of mold-filling, which is especially critical when simulating pressure controlled processes like PC-RTM. The proposed automated specimen preparation and injection process helps to reduce this error significantly.

# **6 Chemo-thermo-mechanics of Resins**

The material behavior of a resin during a typical process cycle can be very complex. This is mainly due to the fact that the resin undergoes significant morphological changes during this time. These changes affect many material parameters like stiffness and relaxation times but also glass transition temperature. The latter evolves during cross-linking and is therefore closely linked to the achieved cure degree. Hence, the temperature dependency of the material behavior permanently changes during curing. At temperatures above Tg, the resin exhibits rubbery properties and transitions to a glassy state if T<sup>g</sup> is undercut. This transition zone is complex to model but at the same time is very important for accurate process simulation since the curing process inevitably ends within this zone. The reason for this is that once the T<sup>g</sup> has reached a value close to the actual cure temperature, the cure rate vanishes as the resin vitrifies and T<sup>g</sup> stops increasing. This means that at the end of the curing process and just before demolding, the resin exhibits a distinct viscoelastic behavior.

In this chapter, three different modeling approaches are considered: viscoelastic, path-dependent and *Cure hardening instantaneously linear elastic* (CHILE). Before detailing of the individual models, the general viscoelastic material behavior is investigated with respect to influence from processing using representative rheological models. Furthermore, the transition from rubbery to glassy state during isothermal processing is investigated with a focus on the frequency that acts during curing. Mathematical modeling is first detailed based on isotropic material behavior and subsequently adopted for describing transversely isotropic materials.

# **6.1 Introduction and review of related work**

Developing a suitable material model, which is able to capture the most dominant influences, demands for a good understanding of the process and its peculiarities. The aim of this chapter is to provide mathematical models which can be used to predict process-induced distortions of composite structures manufactured within a closed mold. It therefore makes sense to analyze the origin of these distortions first in order to focus on the most important influencing factors during model development.

## **6.1.1 Sources of process-induced residual strains**

Process-induced distortions (PID) are caused by residual stresses which are a direct consequence of what happens to the composite material during manufacturing. Main drivers of this phenomenon are strains due to chemical shrinkage and thermal expansion or contraction. Figure 6.1 schematically shows these strains during a typical process cycle. At the beginning of the process, the barely preheated liquid resin is injected into the hot mold and quickly adapts to the chosen cure temperature Tc. During this period, the resin expands with a coefficient of thermal expansion (CTE) that is significantly higher than that of the resin at the end of the process. This is due to the fact that the mold temperature is much higher than the glass transition temperature of the resin. Hence, the CTE of the resin in rubbery state is active, which is much higher than that of the glassy state. However, this has no effect on the residual stresses at the end of the process since the resin reaches the temperature of the mold quickly and well ahead of gelation. Thus, strains that arise during this period cannot persist due to the dominant viscous behavior of the resin.

At some point during curing, the point of gelation is reached and strains start to contribute to the build-up of stresses. Since the curing process is carried out under isothermal conditions, only strains due to chemical shrinkage arise during this process step. At the same time, the resin undergoes substantial morphological

**Figure 6.1:** Schematic of process-induced strains during RTM cycle

changes as the cross-linking process advances. This is followed by a change in stiffness of multiple magnitudes which directly affects the residual stress level. After the cure rate has vanished, the resin reaches a cure degree that is associated with the applied cure temperature. This happens as soon as the T<sup>g</sup> surpasses T<sup>c</sup> and vitrification hinders any further progress in cross-linking reactions and, thus, cure degree. Afterwards, demolding is carried out and the structure starts to cool down to room temperature. During this cooling stage, the resin quickly finishes its transition to the glassy state. Therefore, during the remaining cooling process, the much lower glassy CTE is active which is indicated in Figure 6.1 by the smaller slope of the volume change compared to that of the heating during mold filling. This time, however, the resulting thermal strains directly contribute to residual stresses due to the dominant elastic behavior of the resin in glassy state.

### **6.1.2 Mathematical modeling**

The instantaneous stiffness of a thermoset resin is a result of the interaction of intrinsic properties like degree of cure or relaxation times, and extrinsic properties given by process temperature and applied strain rates or frequency. Consequently, mathematical modeling of this complex material behavior is difficult and followed by a great effort in experimental material characterization. Another aspect is computational efficiency that greatly depends on the numerical complexity of the involved material models. In general, simple models are based on linear elastic approaches whereas complex models take the viscoelastic properties of the polymer into account.

#### Viscoelastic models

Viscoelastic models offer great flexibility in modeling polymer material behavior over a wide range of temperature and frequency. Furthermore, they allow the material to relax or creep, which are phenomena often observed for many types of polymers, especially if these materials are used at elevated temperatures near to the glass transition temperature. The latter is typically the case in liquid composite molding processes like RTM. A great introduction in polymer viscoelasticity can be found in Bergström [105].

Utilizing the Boltzmann superposition principle, the stress-strain relation of a viscoelastic material in integral form is given by

$$
\sigma\left(t\right) = \int\_{-\infty}^{t} E\left(t - \xi\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi \tag{6.1}
$$

in which E and ε denote the relaxation modulus and applied strain, respectively. The integration variable ξ refers to the time at which a specific strain increment has been applied to the material. The lower limit of the integral represents the dependency of the stress evolution on history. Since the resin is a liquid mixture of its individual components at the beginning of the process, no strain history prior to this point in time exists and the integral can be rewritten to

$$
\sigma\left(t\right) = \int\_0^t E\left(t-\xi\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi.\tag{6.2}
$$

Different models are available for modeling of the relaxation modulus and its dependency on cure degree and temperature. Kim and White [106] and O'Brien et al. [40] reported good suitability for modeling epoxy resins using a Prony series given by

$$G\left(t, T, \alpha\right) = G\_{\rm r}\left(\alpha\right) + \sum\_{p=1}^{n} G\_{p}\left(\alpha\right) \exp\left(-\frac{t}{\rho\_{p}\left(T, \alpha\right)}\right). \tag{6.3}$$

Simon et al. [107] used another commonly used formulation of this equation based on weighting factors:

$$G\left(t, T, \alpha\right) = G\_{\rm r}\left(\alpha\right) + \left[G\_{\rm g}\left(\alpha\right) - G\_{\rm r}\left(\alpha\right)\right] \sum\_{p=1}^{n} g\_{p} \exp\left(-\frac{t}{\rho\_{p}\left(T, \alpha\right)}\right). \tag{6.4}$$

In both equations, G<sup>g</sup> and G<sup>r</sup> denote glassy and rubbery modulus, respectively. As shown by other researchers, the glassy modulus can be assumed to be independent of cure degree [40, 106, 108]. The opposite is true for relaxation times ρp, which depend heavily on cure degree and temperature as shown by O'Brien et al. [40] and Hojjati et al. [109]. The principles of time-temperature and timecure-degree superposition are often used to model the shift in relaxation time due to temperature and cure degree changes:

$$
\rho\_p\left(T,\alpha\right) = a\left(T,\alpha\right)\tau\_p\left(T\_{\text{ref}},\alpha\_{\text{ref}}\right).\tag{6.5}
$$

Here, τ<sup>p</sup> refer to the relaxation times which belong to a reference temperature and cure degree, for which a master curve has been generated from measurements. This master curve controls shifting of the relaxation spectra with respect to changes in temperature as well as cure degree and is represented in the equation by the shift-factor a. The physical meaning behind this approach is that an increase in temperature has a similar effect as loading the material during a longer period of time or with a reduced frequency. Time-Temperature-Superposition (TTS) is often modeled using WLF based equations like

$$\log\_{10} a\left(T\right) = \frac{d\_1 \left(T - T\_{\text{ref}}\right)}{d\_2 + T - T\_{\text{ref}}}\tag{6.6}$$

or Arrhenius-type approaches

$$\log\_{10} a\left(T\right) = -\frac{E\_A}{2.303 \cdot \text{R}} \left(\frac{1}{T\_{\text{ref}}} - \frac{1}{T}\right) = -\frac{d\_\text{A}}{2.303} \left(\frac{1}{T\_{\text{ref}}} - \frac{1}{T}\right). \tag{6.7}$$

In these two equations, d1, d<sup>2</sup> and d<sup>A</sup> denote model parameters whereas R and E<sup>A</sup> are ideal gas constant and activation energy. Tref is again the reference temperature, at which the master curve, that is represented by the shift function a (T), has been constructed and parametrized.

Another possibility for modeling the shift factor is given by a generalized Vogel equation which is able to represent WLF and Arrhenius like behavior depending on the chosen parameters [107]:

$$\log\_{10} a\left(T\right) = \frac{d\mathbf{v}}{T - T\_{\infty}} - \frac{d\mathbf{v}}{T\_{\text{ref}} - T\_{\infty}}.\tag{6.8}$$

Using the Vogel equation allows to model the temperature dependency of the relaxation spectrum differently above and below Tg. This may be needed if the material shows a strong change in relaxation characteristic after exceeding or undercutting Tg. As pointed out by Ferry [110], a WLF type equation may be invalid below Tg. Instead, the Arrhenius equation should be used in that temperature range according to Brinson [111].

Time-cure degree superposition can be assumed to follow a similar relationship as time-temperature superposition. Plazek and Chay [108] investigated the creep compliance of an epoxy resin at different cure states. They found that the shape of the shift function is independent of cure degree if the current glass transition temperature T<sup>g</sup> (α) is used as reference. A similar conclusion is drawn by Ferry [110] for changes in molecular weight of a polymer. Hence, Equation 6.8 can be rewritten to [107]

$$\log\_{10} a\left(T, \alpha\right) = \frac{d\_{\rm V}}{T - T\_{\infty}} - \frac{d\_{\rm V}}{T\_{\rm g}\left(\alpha\right) - T\_{\infty}}.\tag{6.9}$$

Figure 6.2 shows WLF and Arrhenius type shift functions for different cure states of a generic resin that exhibits a maximum Tg,<sup>∞</sup> of 115 ◦C. Each individual curve reaches their zero value at a temperature that equals the T<sup>g</sup> of the current cure state. The WLF based model shows a rather broad spreading of the shift factor towards low temperatures and a convergence of the predicted curves for high temperatures, well above Tg. The opposite is true in case of the Arrhenius type equation which emphasizes the different characteristic of both models.

**Figure 6.2:** Comparison of Arrhenius and WLF based shift function formulations and different cure degrees. In both models, <sup>T</sup>g (α) is used as the reference temperature <sup>T</sup>ref.

Proper calibration of the shift factor is vital for correct prediction of relaxation processes. However, in order to achieve this goal, extensive experimental characterization efforts are required, which is one of the reasons why simpler approaches, like path-dependent models, have been developed.

#### Path-dependent models

Compared to a holistic viscoelastic approach, path-dependent models represent a simplification with regard to the modeled relaxation time spectrum and rate dependency. As such, they can be derived from viscoelastic models using the assumption of a relaxation time of zero in rubbery state and infinity in glassy state. Hence, frozen-in stresses are captured by these models but are completely released once the material devitrifies, neglecting any time dependency. Svanberg and Holmberg [112] developed a path-dependent model in which the stress increment is given by

$$\Delta \sigma\_{ij} = \begin{cases} C\_{\text{r},ijkl} \Delta \varepsilon\_{kl} - S\_{ij} \left( t\_{m-1} \right), & T \ge T\_{\text{g}} \left( \alpha \right) \\ C\_{\text{g},ijkl} \Delta \varepsilon\_{kl}, & T < T\_{\text{g}} \left( \alpha \right) \end{cases} \tag{6.10}$$

Here, the history state variable Sij is defined as

$$S\_{ij}\left(t\_m\right) = \begin{cases} 0, & T \ge T\_{\mathbf{g}}\left(\alpha\right) \\ S\_{ij}\left(t\_{m-1}\right) + \left(C\_{\mathbf{g},ijkl} - C\_{\mathbf{r},ijkl}\right)\Delta\varepsilon\_{kl}, & T < T\_{\mathbf{g}}\left(\alpha\right) \end{cases} \tag{6.11}$$

In this equations, Sij contains the frozen-in stresses which arise from the stiffness difference between glassy and rubbery state. For temperatures below T<sup>g</sup> these frozen-in stresses persist and increase in case of non-zero strain increments. As soon as the resin devitrifies through heating above Tg, the frozen-in stresses are completely released. This simplification of the otherwise complex viscoelastic properties of the resin is based on the assumption that both vitrification and devitrification happens rapidly, such that time-dependent relaxation in the transition region can be neglected.

Due to the history term in path-dependent models, they offer a broader application range compared to simple elastic models like CHILE models. While the latter is only applicable for processes without devitrification, path-dependent models are also suitable in processes which include reheating above T<sup>g</sup> like two-step curing cycles. Furthermore, they can be extended by temperature dependency of rubbery and glassy modulus to be able to describe the behavior of Thermo-rheological Complex Materials (TCM) as shown by Ding et al. [68].

#### CHILE models

*Cure hardening instantaneously linear elastic* (CHILE) models are a specialized class of linear elastic models. The stress evolution follows a relationship similar to that of viscoelastic models (cf. Equation 6.2) but lacks the history dependency:

$$\sigma\left(t\right) = \int\_0^t E\left(\xi\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi = \int\_0^t E\left(T, \alpha\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi. \tag{6.12}$$

In this equation, E denotes the instantaneous elastic modulus which depends on temperature T and cure degree α. This modulus for example can be measured using a single rheometer measurement. Hence, the required experimental effort is much lower compared to viscoelastic models.

CHILE models usually assume a monotonic increase of the resin modulus during cure. Johnston [59] developed a CHILE model, in which the stiffness of the resin is modeled by

$$E\left(T,\alpha\right) = \begin{cases} E\_{\rm r}, & T\_{\rm g}\left(\alpha\right) - T < \Delta T\_{c1} \\ E\_{\rm r} + \frac{T\_{\rm g}\left(\alpha\right) - T - \Delta T\_{c1}}{\Delta T\_{c2} - \Delta T\_{c1}} \left(E\_{\rm g} - E\_{\rm r}\right), & \Delta T\_{c1} \le T\_{\rm g}\left(\alpha\right) - T \le \Delta T\_{c2} \\ E\_{\rm g}, & T\_{\rm g}\left(\alpha\right) - T > \Delta T\_{c2} \end{cases} \tag{6.13}$$

In this model, the rubbery modulus E<sup>r</sup> is active as long as T<sup>g</sup> has not exceeded the current temperature T by a value given by the variable ΔTc1. This relationship can also be interpreted the other way around, assuming a fixed value of T<sup>g</sup> and varying T. The stiffness is equal to the glassy modulus if, e.g. during cooling, the temperature T reaches a value that is lower than T<sup>g</sup> − ΔTc2. Between these two material states, a transition is described by the model which depends on the current difference in T and Tg. As noted earlier, CHILE models are not able to account for frozen-in stresses. Hence, reheating to a temperature above T<sup>g</sup> + ΔTc1 does change the instantaneous stiffness but is not followed by a release of frozen-in stresses.

Using a CHILE model to describe the behavior of a viscoelastic material can be prone to errors since the actual material still exhibits viscoelastic properties. Hence, while CHILE models are less demanding in terms of computational modeling and effort, they need to be parametrized carefully such that they are still able to represent the materials behavior as accurately as possible. While most of the complexity of viscoelastic approaches concentrates on modeling the shift factors for both temperature and cure influence, this is all incorporated into the evolution of the modulus in a CHILE approach. Thus, although a single rheometer measurement may be enough to calibrate such a model, this single measurement needs to take into account as much detail of the actual process as possible. It is therefore important to match the frequency of the measurement with that occurring during the course of the process. This is often neglected and a common but arbitrarily chosen frequency of 0.1 Hz or 1 Hz is used instead. As pointed out by Zobeiry et al. [113], this can have significant influence on the predicted stress levels, which is why they introduced a special class of properly calibrated CHILE models, namely *pseudo-viscoelastic* models. In these models, the elastic modulus is modeled according to the actual frequency that acts during the course of the process. Hence, they combine the computational efficiency of linear elastic models with the accurate prediction of time and rate-dependent processes of viscoelastic approaches. However, since comprehensive information on the viscoelastic properties of the resin is needed for proper calibration, the experimental effort increases and can reach similar extends as that needed by viscoelastic models.

#### Influence of curing on material properties

Irrespective of the chosen class of constitutive model, the cure degree dependency of the resin modulus in both rubbery and glassy states needs to be predicted using a suitable model. Gupta et al. [114] investigated the influence of crosslink density on different properties of an epoxy resin. They found that in glassy state, the cross-link density has a minor influence on small-strain properties like stiffness or yield point while it has a bigger impact on large-strain properties like strength, elongation at break and toughness. In rubbery state, the crosslink density has a major influence on both, small and large-strain properties. A similar conclusion is reported by Johnston [59] and Adolf and Martin [115]. With respect to the aim of this study, predicting process-induced distortions, this means we can safely assume the glassy modulus to be independent of cure degree and concentrate on the much more sensitive rubbery modulus instead.

Bogetti and Gillespie [116] introduced a model for the cure dependency of the resin modulus, given by

$$E\_{\mathbf{r}}\left(\alpha\right) = \left(1 - \alpha\_{\rm mod}\right) E\_{\mathbf{r}}^{0} + \alpha\_{\rm mod} E\_{\mathbf{r}}^{\infty} + \gamma \alpha\_{\rm mod} \left(1 - \alpha\_{\rm mod}\right) \left(E\_{\mathbf{r}}^{\infty} - E\_{\mathbf{r}}^{0}\right) \tag{6.14}$$

with

$$
\alpha\_{\rm mod} = \frac{\alpha - \alpha\_{\rm gel}^{\rm mod}}{\alpha\_{\rm diff}^{\rm mod} - \alpha\_{\rm gel}^{\rm mod}}.\tag{6.15}
$$

In these equations, E<sup>0</sup> <sup>r</sup> and E<sup>∞</sup> <sup>r</sup> denote the resin rubbery modulus in fully uncured and fully cured state, respectively. Furthermore, the parameter γ can be used to adjust the balance between stress relaxation and chemical hardening. If, for example, a chosen resin cures at such a rapid rate that the increase in stress exceeds the amount of stress that relaxes within the same time span, the parameter γ can be used to account for this. However, it should be noted that these additional stresses will not be able to relax afterwards or during the remaining part of the curing stage of the process. Hence, γ will often be set to 0 and stress relaxation during curing should be taken into account by using a fully viscoelastic approach instead. However, γ may also be used to tailor the model to resins that show different growth rates of the modulus during curing.

Adolf and Martin [115] developed a model for the resin modulus evolution based on a power law:

$$G\_{\rm r} = G\_{\rm r} \left( \alpha = 1 \right) \cdot \left( \frac{b - b\_{\rm gel}}{1 - b\_{\rm gel}} \right)^{8/3} \,. \tag{6.16}$$

In this equation, b denotes a bond probability which is closely related to the cure degree α. In case of epoxy and polyurethane resins at their stoichiometric ratio, b = α<sup>2</sup> [115]. Moreover, they assume an insignificantly small rubbery modulus prior to gelation. Similar to the parameter γ of Equation 6.14, the exponent of Equation 6.16 can also be used to tailor the model to different resins.

# **6.2 Experimental characterization of chemical shrinkage and thermal expansion**

Numerous material properties need to be quantified in order to achieve meaningful simulation results. While some of these properties are given by technical data sheets provided by the manufacturer of the resin system, others require additional measurements. Properties which are usually not included in technical data sheets are coefficients for thermal expansion and chemical shrinkage as well as the cure degree of the resin, at which gelation occurs. The latter has been characterized in Section 3.4. Determination of the remaining properties is detailed in the following.

#### Chemical shrinkage

Different methods exist for characterization of this property. Chemical shrinkage of resins SCER1, FCPUR and other resins is investigated in Bernath et al. [15] using parallel-plate rheometry as well as video-imaging. Comparison of these two characterization methods in general showed higher effective shrinkage strains when using video-imaging. This is caused by different boundary conditions acting during the measurement. In a parallel-plate rheometer setup, the resin is only allowed to shrink perpendicular to the plate while it is allowed to deform freely during a Video-Imaging (VI) measurement. Hence, the measured value is neither the volumetric nor linear shrinkage of the resin when using rheometry. Consequently, values gathered from rheometer measurements are generally smaller than those measured by video-imaging [15]. Hence, if available, videoimaging measurements are used for determination of chemical shrinkage. This is the case for all resins except FCER3, for which rheometer measurements have been used instead.

Rheometer and VI measurements of all resins except FCER3 were carried out at the Institute of Composite Structures and Adaptive Systems at the German Aerospace Center (DLR) in Braunschweig, Germany. Rheometer measurements of FCER3 resin were conducted at the Fraunhofer Institute for Chemical Technology (ICT) in Pfinztal, Germany.

The Video-Imaging technique is based on optical evaluation of the size of a resin droplet as it shrinks during isothermal curing. For this purpose, the droplet is put on top of a heated cylindrical aluminum socket. Due to surface tension, a nearly spherical droplet is formed. During curing, the specimen's geometry is constantly recorded using a camera and bright backlight. Hence, the assumption of axial symmetry is made as only a 2D cross-section is captured. The data on the droplet's shape is then used to calculate its volume as it decreases with time due to shrinkage. Since this measurement technique is able to capture the whole shrinkage, i.e. also the portion which evolves prior to gelation, the resulting raw shrinkage is much higher than that measured by other techniques like rheometry. However, with respect to residual strains and process-induced distortion, only the portion of shrinkage is relevant, that develops between gelation and full cure. Prior to gelation, the polymeric network behaves like a viscous liquid and therefore is incapable of sustaining any stress. Hence, when using VI data, the effective shrinkage needs to be extracted from raw shrinkage based on the cure degree at gelation using the following equation:

$$
\beta\_{\rm ch} = \beta\_{\rm VI} \cdot \left( \alpha\_{\rm max} \left( T\_{\rm c} \right) - \alpha\_{\rm gel} \right). \tag{6.17}
$$

In this equation αmax (Tc) denotes the maximum achievable cure degree of the resin when cured using the isothermal temperature Tc. If full cure is reached, this quantity becomes one. Moreover, as the cure degree at gelation is needed, shrinkage measurement using VI must be accompanied by separate measurement of αgel, as VI data does not provide any information on gelation. Shrinkage measurement using rheometry does not show this shortcoming as the principle relies on measuring the movement of one of the cylindrical plates to compensate for the tension forces caused by shrinkage. Hence, prior to gelation, no shrinkage is recorded. This enables to measure both shrinkage as well as time of gelation simultaneously. As this technique has been used for characterization of resin FCER3, it is detailed in the following.

When characterizing shrinkage using a rheometer, a small amount of resin is put between two cylindrical metal plates (cf. Figure 4.2). Prior to measurement start, the resulting gap between these plates is measured by the rheometer. During measurement, the rheometer is instructed to only change the gap in case a predefined force limit is exceeded. Hence, during curing and before gelation, no force acts on the plates and therefore no movement is necessary. After gelation, however, chemical shrinkage will lead to tension forces and the rheometer will continuously adjust the gap each time the force limit is exceeded. This way, information on two important properties is gained: shrinkage and time of gelation. As detailed in the related study [15], aluminum plates with a diameter of 40 mm as well as an initial gap of 1 mm are used for all resins but FCER3. Shrinkage of the latter is measured using a plate diameter of 25 mm. The resulting shrinkage is calcluated from the temporal evolution of the gap using the approach presented by Haider et al. [117].

In Figure 6.3, the chemical shrinkage as measured using VI and rheometry is given for the resins used in this study. For reasons of comparison, other resins investigated in Bernath et al. [15] are also included. In this comparison, VI shrinkage is already adjusted according to Equation 6.17 to represent only the relevant portion of volumetric shrinkage. Furthermore, only shrinkage data from rheometry is available for resin FCER3.

**Figure 6.3:** Chemical shrinkage of the resins used in this study and in Bernath et al. [15]

For all resins but FCER2, the shrinkage measured using rheometry is lower than that reported by video-imaging. In case of FCPUR and SCER2, a 20 % to 30 % lower value of shrinkage is found by rheometry, while only a minor difference is found for resin SCER1. Moreover, only in case of resin FCER2, this tendency is reversed. Hence, no distinct relationship between shrinkage measured by VI and rheometry can be derived and used to predict VI shrinkage of resin FCER3. This discrepancy can be explained by several possible reasons. Firstly, an error in the determination of the cure degree at gelation can lead to an underestimation of the relevant shrinkage as measured by VI. In fact, this could explain the reversed tendency in case of FCER2. Another possible source of error may be introduced by Equation 6.17, which implies a linear relationship between cure degree and chemical shrinkage strain. This assumption has been investigated and confirmed for SCER1 resin in Bernath et al. [15] and is supported by other authors as reported by Kravchenko et al. [118]. Garstka et al. [119] report a slightly nonlinear behavior towards full cure for isothermal shrinkage measurements of Hexcel AS4/8552 epoxy prepreg material. However, since they used solely nonisothermal DSC data for fitting of the kinetic model used to predict the cure

evolution during shrinkage measurements, the error might be due to insufficient information on vitrification. This is supported by the fact that measurements conducted at isothermal temperatures of 150 ◦C, 160 ◦C and 170 ◦C all show nonlinear behavior from about 70 % cure on, while that of the highest temperature of 180 ◦C does show much less. Another possible explanation for the differences may be related to the fact that in Garstka et al. [119], a cross-ply stack of prepreg material is investigated, while in Bernath et al. [15], neat resin is characterized. Consequently, different mechanical boundary conditions act during curing and, in case of the cross-ply layup, develop gradually during this process as the bond between fibers and resin becomes constantly stronger. This dynamically changing boundary conditions are also present in shrinkage measurements using rheometry. In this case, in-plane shrinkage strain is hindered by the cylindrical metal plates.

#### Thermal expansion

The coefficients of thermal expansion (CTE) of the individual resins are measured using thermomechanical analysis (TMA), which enables to track properties of a material which are caused by or affected by changes in temperature. Hence, it is often applied for characterization of thermal expansion of polymers as has been done in this study for resins SCER1 and FCPUR, and also for other resins in the related study of Bernath et al. [15]. Corresponding measurements were carried out at the Institute of Composite Structures and Adaptive Systems at the German Aerospace Center (DLR) in Braunschweig, Germany.

The specimens are usually small flat plates or cylinders, cut from bigger thin plates of neat resin that have been cured using the cycles given by Table 2.3. After cutting, the specimens are heated twice from 20 ◦C to 140 ◦C using a constant heat rate of 3 K min−1. The first run aims at identifying the properties of the material while it is in a state similar to that at the end of the curing process in a closed mold process. Due to the non-isothermal nature of the measurement and the chosen high end temperature of 140 ◦C, each heating run ends at a temperature significantly above Tg. Hence, the resin further cures and eventually reaches full cure. Consequently, the properties measured during the second run are that of the fully cured material. Furthermore, since no further curing happens after the temperature surpasses the final Tg,<sup>∞</sup> of the resin, evaluation of the CTE in rubbery state is possible.

**Figure 6.4:** Coefficient of thermal expansion (CTE) of the used resins in glassy and rubbery state. For reasons of comparison, other resins investigated in the related study of Bernath et al. [15] are included. Rubbery CTE of resins marked with ∗ are predicted based on the ratio of rubbery to glassy CTE of the remaining resins.

Figure 6.4 shows the measured CTE of the used resins in rubbery and glassy state. In case of resin FCER3, only the value for glassy state is provided by the technical data sheet. The CTE in rubbery state of this resin but also SCER2 is therefore estimated by using the mean value of the ratio of rubbery to glassy CTE of the remaining resins, which is 2.49 ± 0.17. The relatively small variance of this ratio is remarkable. However, the rubbery CTE is only rarely relevant since in RTM or WCM, cooling is usually started once cross-linking has finished. Consequently, the glassy CTE dominates the build-up of thermal strain during the cooling stage.

## **6.3 Mathematical modeling**

The mathematical modeling used in this thesis is described in the following. During curing, chemical shrinkage and mechanical properties develop. As this process takes place while the material is subject to geometrical boundary conditions caused by the mold, the resulting amount of chemical shrinkage is altered. Moreover, the resin is inherently viscoelastic in nature. Hence, the response to strains from shrinkage, temperature change or demolding depends on the rate, at which these phenomena act. Whether this complex interaction can be represented by a more computational efficient elastic approach, depends on the degree to which viscoelastic effects are pronounced. In this respect, it is especially crucial to correctly determine the point of vitrification, since it is accompanied by significant changes in mechanical but also viscoelastic properties of the resin.

### **6.3.1 Effective shrinkage strain**

The mechanical properties of a resin depend on the advancing cross-linking process. A common model representation of this effect makes use of an array of springs. At gelation, the material for the first time exhibits notable mechanical properties and reacts to applied forces with stress rather than viscous flow. Therefore, the representative model consists of a single spring at the time of gelation. The stiffness increases as the resin cures further, which can be accounted for in the model by introducing more and more springs. This is shown in Figure 6.5 for an idealized material with a cure degree at gelation of αgel = 50 %. In this example, a spring is added to the network every 10 % of cure for the sake of clarity. Furthermore, the black dashed line in Figure 6.5 indicates the unloaded length of the individual springs at which they are introduced in the model.

**Figure 6.5:** Schematic of chemical shrinkage during cure in a closed mold process

Figure 6.5a shows the result of an increase in cure degree from 50 % to 60 %. The spring which represents the stiffness at the time of gelation is lengthened by the incremental chemical shrinkage Δεch caused by the increase in cure degree. After another increase in cure degree of 10 %, a second spring is introduced into the model and both springs are further stretched by chemical shrinkage as shown in Figure 6.5b. This mechanism stops when the maximum cure degree of the applied temperature has been reached. In this example, full cure is assumed and the network development is visualized in Figures 6.5a to 6.5e. As is indicated in Figure 6.5e, the individual springs are not equally stretched but exposed to varying degrees of strain. This is an analogy of the growth and cross-linking of a polymeric network. Each newly established cross-link between molecules is inserted into the existing and distorted network in an undeformed state. Hence, a cross-link which is introduced late in the curing process will experience a much smaller amount of chemical shrinkage as one that is established close to gelation. This fact has an impact on the effective chemical strain that occurs after the constraints of the closed mold are removed as is indicated in Figure 6.5f. In this case, the tension which is present in the springs that were added at an early stage of the curing process leads to a compression of the remaining springs. Consequently, the amount of chemical shrinkage of a resin cured in a constrained state differs from that of the very same resin when it is allowed to deform freely during cure. Furthermore, the balance of tension and compression in this model heavily depends on the evolution of the modulus of the resin. In this study, the model of Adolf and Martin [115], Equation 6.16, is used and the evolution of rubbery shear modulus with cure degree α is given by

$$G\_{\rm r} = G\_{\rm r,\infty} \cdot \left(\frac{\alpha^2 - \alpha\_{\rm gel}^2}{1 - \alpha\_{\rm gel}^2}\right)^{\kappa}.\tag{6.18}$$

Figure 6.6 shows the development of the modulus of a generic resin which gels at a cure degree of 50 % and reaches a maximum rubbery shear modulus of Gr,<sup>∞</sup> of 10 MPa. The parameter κ of Equation 6.18 allows to adjust the rate at which the final shear modulus Gr,<sup>∞</sup> is approached, as is shown in

**Figure 6.6:** Modulus development during cure of a generic resin based on Equation 6.18 for different values of parameter κ

Figure 6.6 for three different values of κ. By increasing this value, a more gentle increase in modulus is obtained. Consequently, springs introduced near to gelation will possess a rather low stiffness since the increments in modulus are small during this period. Hence, the tensile stresses generated by these springs due to shrinkage and thus also the effective chemical shrinkage εeff ch will be lower compared to a model that uses a smaller value of κ. To investigate this further, a one-dimensional mathematical model describing the development of stress in the model introduced above is derived in an incremental form.

Prior to gelation the resin reacts with viscous flow to applied forces, and stress is therefore assumed to be effectively zero during this period. For reasons of simplicity, the time of gelation is referred to as t<sup>0</sup> for the following explanation, hence it follows for the zeroth stress increment

$$
\sigma^{(0)} = 0.\tag{6.19}
$$

After some time, the resin has cured further and therefore gained stiffness but also chemical shrinkage. In the representative model explained above, this equals to a single spring of stiffness k<sup>1</sup> which is exposed to a shrinkage strain Δεch (cf. Figure 6.5a). Mathematically, this can be expressed as

$$
\sigma^{(1)} = \sigma^{(0)} + k^{(1)} \cdot \Delta \varepsilon\_{\text{ch}} = \sigma^{(0)} + k^{(1)} \cdot \left[ \varepsilon\_{\text{ch}}^{(1)} - \varepsilon\_{\text{ch}}^{(0)} \right]. \tag{6.20}
$$

It follows for the stress of the n-th time increment

$$
\sigma^{(n)} = \sigma^{(n-1)} + k^{(n)} \cdot \left[ \varepsilon\_{\rm ch}^{(n)} - \varepsilon\_{\rm ch}^{(n-1)} \right] \tag{6.21}
$$

where the stiffness of the network k(n) equals the sum of the stiffness of the individual springs

$$k^{(n)} = \sum\_{i=1}^{n} k\_i. \tag{6.22}$$

In order to relate this to the evolving macroscopic stiffness of the resin in rubbery state Er, the spring stiffness k<sup>i</sup> is defined as

$$k\_i = E\_\mathbf{r}^{(n)} - E\_\mathbf{r}^{(n-1)}.\tag{6.23}$$

While the discrete spring model is well suited to visualize the underlying mechanism, a more general form is needed for numerical process modeling. Therefore, the temporal evolution of stress is derived from an integral form, which follows from the representative model in case of an infinite number of springs:

$$
\sigma\left(t\right) = \int\_0^t E\_\mathbf{r}\left(\xi\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi. \tag{6.24}
$$

From this, the stress at a time t<sup>m</sup> is given by

$$
\sigma\left(t\_m\right) = \int\_0^{t\_{m-1}} E\_\mathbf{r}\left(\xi\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi + \int\_{t\_{m-1}}^{t\_m} E\_\mathbf{r}\left(\xi\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi
$$

$$
= \left.\sigma\left(t\_{m-1}\right) + \int\_{t\_{m-1}}^{t\_m} E\_\mathbf{r}\left(\xi\right) \frac{\mathrm{d}\varepsilon\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi\right.\tag{6.25}
$$

Assuming linear variation in time of both stiffness E<sup>r</sup> and strain εch, an incremental equation similar to Equation 6.21 is found:

$$
\sigma\left(t\_m\right) = \sigma\left(t\_{m-1}\right) + E\_\mathbf{r}\left(t\_{m-1/2}\right) \cdot \Delta\varepsilon\_{\text{ch}}\left(t\_m\right) \tag{6.26}
$$

where E<sup>r</sup> t<sup>m</sup>−1/<sup>2</sup> denotes the midpoint value of the stiffness and is given by

$$E\_{\mathbf{r}}\left(t\_{m-1/2}\right) = \frac{1}{2} \cdot \left[E\_{\mathbf{r}}\left(t\_m\right) - E\_{\mathbf{r}}\left(t\_{m-1}\right)\right].\tag{6.27}$$

Since cross-linking is the driving force behind shrinkage stress, the time dependency of stiffness and strain in Equation 6.26 is exchanged with cure-dependent values. Furthermore, the shear modulus model given by Equation 6.18 is used in combination with the relationship

$$E\_{\mathbf{r}} = 2G\_{\mathbf{r}} \cdot \left(1 + \nu\_{\mathbf{r}}\right),\tag{6.28}$$

finally leading to

$$
\sigma\left(t\_m\right) = \sigma\left(t\_{m-1}\right) + E\_\mathbf{r}\left(\alpha\_{m-1/2}\right) \cdot \Delta\varepsilon\_{\text{ch}}\left(\alpha\_m\right). \tag{6.29}
$$

Here, α<sup>m</sup>−1/<sup>2</sup> denotes the midpoint value of the cure degree. Moreover, a Poisson's ratio of ν<sup>r</sup> = 0.5 is used in Equation 6.28, which is a valid assumption for a resin in rubbery state. Furthermore, the shrinkage strain increment Δεch is given by

$$\Delta \varepsilon\_{\text{ch}} \left( \alpha\_m \right) = \begin{cases} \frac{\beta\_{\text{ch}}}{1 - \alpha\_{\text{gel}}} \cdot \left[ \alpha \left( t\_m \right) - \alpha \left( t\_{m-1} \right) \right], & \alpha \ge \alpha\_{\text{gel}} \\ 0, & \alpha < \alpha\_{\text{gel}}. \end{cases} \tag{6.30}$$

Distinction of these two cases is needed to ensure zero stress build-up before gelation. Consequently, the effective shrinkage coefficient βch needs to be rescaled to the full range of the cure degree α = [0, 1]. This is achieved by dividing the shrinkage coefficient by (1 − αgel).

**Figure 6.7:** Evolution of shrinkage stress during cure of a generic resin for different values of parameter κ

Figure 6.7 shows the stress evolution of the generic resin introduced above for different values of parameter κ, calculated using Equation 6.29. As is expected, the characteristic of the cure dependency of the modulus has significant impact on the amount of stress generated during cross-linking. Due to the overall higher stiffness in case of p = 1, the shrinkage stress rises faster and eventually reaches a higher value.

The predicted stress evolution can now be used to calculate the effective shrinkage strain

$$
\varepsilon\_{\rm ch}^{\rm eff}(t) = \frac{\sigma\left(t\right)}{E\_{\rm r}\left(\alpha\left(t\right)\right)},\tag{6.31}
$$

which represents the actual relevant shrinkage strain with respect to processinduced distortions. Figure 6.8 shows the effective shrinkage strain again for three different values of parameter κ. By increasing the value of κ, the gently increasing modulus shifts most of the shrinkage-induced stress towards the end of the curing process. Hence, the tensile stress introduced right after gelation is much lower and thus also the restoring force of the network, leading to only a small amount of effective shrinkage strain.

**Figure 6.8:** Evolution of shrinkage strain during cure of a generic resin for different values of parameter κ

### **6.3.2 Viscoelasticity during processing**

In the above discussion, the material is assumed to behave purely elastic. This assumption is based on two facts that are typically applicable to RTM processes. Firstly, most of the time during curing, the cure temperature is significantly larger than the current Tg. Hence, the material is in rubbery state, in which relaxation times are short, especially when compared to the duration of the curing process. Secondly, the volumetric strain due to chemical shrinkage is small and typically in a range of 1 % to 4 % (cf. Figure 6.4). Hence, strain rates caused by chemical shrinkage are very small and the material therefore reacts purely elastic with a stiffness that is represented by the cure-dependent rubbery modulus G<sup>r</sup> (α). The mechanism is illustrated using a *Standard linear solid model* in Figure 6.9(a) and (b). This representative viscoelastic model consists of a spring of rubbery stiffness G<sup>r</sup> (α) and a Maxwell element, which are put in parallel to each other. Moreover, a strain sink for taking into account the chemical shrinkage strain εch (α) acts in a serial manner. The stiffness of the spring in the viscoelastic branch is defined by the difference of the glassy and rubbery modulus. This way, the stiffness of the representative model equals that of the glassy material in case of low temperatures or high frequencies, at which the damper locks. For high temperatures or low frequencies, the stiffness reduces to that of the elastic branch. Intermediate states between these two extremes depend on the relaxation time τ , which is given by the ratio of damper viscosity η<sup>D</sup> and spring stiffness of the Maxwell element:

$$\tau = \frac{\eta\_{\rm D}}{E\_{\rm g} - E\_{\rm r}(\alpha)} \iff \eta\_{\rm D} = \tau \cdot \left[E\_{\rm g} - E\_{\rm r}(\alpha)\right]. \tag{6.32}$$

The viscosity of the damper decreases if the relaxation time is lowered e.g. due to temperature increase. The stress generated by the damper is related to the viscosity by the applied strain rate as given by

$$
\sigma\_{\rm D} = \eta\_{\rm D} \frac{\mathrm{d}\varepsilon}{\mathrm{d}t} = \eta\_{\rm D} \dot{\varepsilon}. \tag{6.33}
$$

During cure, the temperature remains well above T<sup>g</sup> for most of the time. Hence, relaxation time and thus damper viscosity are low. Additionally, the strain rate as defined by the chemical shrinkage is also small. As a result, a negligible amount of stress is build-up in the Maxwell element during cure, while the elastic branch is subjected to strains causing tensile stress. Since shrinkage strain and stiffness are both cure-dependent, the phenomenon discussed in Section 6.3.1 applies here and the single spring of the elastic branch represents the multiple-spring assembly shown in Figure 6.5.

The viscoelastic branch, represented by the Maxwell element in Figure 6.9, starts to contribute to the overall stress development as soon as the stress increase exceeds the stress relaxation in the same time span. Since the viscosity of the damper increases either due to an advancement in cure degree or a reduction in temperature, this is potentially possible at the time of demolding. Here, T<sup>g</sup> has reached or even outrun the reaction temperature and a rather high cure degree is achieved. However, the viscoelastic contribution to stress and thus stiffness increase still depends on strain rate. If the strain rate is high, e.g. by assuming an instantaneous demolding of the part, the damper effectively hinders relaxation of the stretched spring of the elastic branch. Instead, the spring of the Maxwell element gets compressed. This is shown in Figure 6.9(c).

The complex process of demolding is often neglected during process design and associated research. Hence, the assumption of a quasi-instantaneous demolding is made in this study. In reality, the strain rate may differ largely throughout a composite structure since the same applies to the effectiveness of release agent, friction between mold and part as well as spatial distribution of demolding forces. Nevertheless, by using a viscoelastic material model, the simulation in theory is able to incorporate such uncertainties. The situation is more difficult in case of elastic or path-dependent models which do not explicitly model strain rate dependencies. Right after (instantaneous) demolding, a viscoelastic model shows the idealized state given in Figure 6.9(c). Since the Maxwell element prevents the whole model from reaching the state of full shrinkage strain, corresponding process-induced distortions will also be much smaller. Whether both shrinkage strain and process-induced distortions reach their ultimate value depends on

**Figure 6.9:** One-dimensional representative model of viscoelastic closed-mold polymer curing and demolding

what happens in the further course of the process. In reality, cooling of the part down to room temperature is either done in a controlled or uncontrolled fashion. Thus, the cooling rate, which decides for how long further relaxation can happen in the Maxwell element and what state of process-induced distortion is reached, may differ largely. If the cooling rate is high, a great portion of shrinkage might be frozen in the polymeric structure and may only show up if the part is reheated to a temperature close to or above Tg. In case of a more gradual cooling, the stresses in the elastic and viscoelastic branch may fully relax, which is shown in Figure 6.9(d). Elastic and path-dependent models are not able to take into account such effects and careful setup is needed to obtain accurate results as otherwise the accuracy of numerically predicted process-induced distortions will suffer significantly.

It should be noted that the above discussion implies isotropic material properties. The structural behavior will of course differ from that in case of anisotropic material behavior. However, the underlying mechanism still remains the same and is most evident in directions whose properties are dominated by the matrix, e.g. perpendicular to the fibers of a unidirectional ply. In the following section, model development will start with a fully viscoelastic approach from which simpler elastic models will be derived by applying the above discussed assumptions.

### **6.3.3 Isotropic viscoelastic model**

In the isotropic viscoelastic model, stress and strain is decomposed into deviatoric and volumetric parts, as proposed by Zobeiry [120]:

$$
\sigma = \sigma\_{\rm dev} + \sigma\_{\rm vol} = 2G\varepsilon\_{\rm dev} + 3K\varepsilon\_{\rm vol}.\tag{6.34}
$$

In this equation, the deviatoric and volumetric strain tensors are defined as

$$
\varepsilon\_{\rm dev} \quad = \quad \varepsilon - \frac{1}{3} \operatorname{tr} \left( \varepsilon \right) I \tag{6.35}
$$

$$
\varepsilon\_{\rm vol} \quad = \quad \frac{1}{3} \operatorname{tr} \left( \varepsilon \right) I,\tag{6.36}
$$

131

where I denotes the identity matrix. In order to simplify the following derivation of the constitutive equations, the prefactors in Equation 6.34 are included in the strain tensors as

$$
\hat{\varepsilon}\_{\text{dev}} = 2 \cdot \varepsilon\_{\text{dev}} \tag{6.37}
$$

$$
\hat{\varepsilon}\_{\text{vol}} \quad = \quad 3 \cdot \varepsilon\_{\text{vol}}.\tag{6.38}
$$

The deviatoric stress and strain are linked by the shear modulus of the resin G, which is modeled using a Prony series as proposed by Simon et al. [107]:

$$\begin{aligned} G\left(t, T, \alpha\right) &= \left. G\_{\mathbf{r}}\left(\alpha\right) + \left[G\_{\mathbf{g}} - G\_{\mathbf{r}}\left(\alpha\right)\right] \sum\_{p=1}^{n} g\_{p} \exp\left(-\frac{t}{\rho\_{p}\left(T, \alpha\right)}\right) \\ &= \left. G\_{\mathbf{r}}\left(\alpha\right) + \Delta G\left(\alpha\right) \sum\_{p=1}^{n} g\_{p} \exp\left(-\frac{t}{\rho\_{p}\left(T, \alpha\right)}\right) .\end{aligned} \tag{6.39}$$

This approach leads to a shear modulus that switches between a cure-dependent rubbery modulus G<sup>r</sup> (α) and a constant glassy modulus Gg. The independence of G<sup>g</sup> from cure degree has been reported by several studies [40, 106, 108]. The characteristic of the switching is defined by the Prony series using weight factors g<sup>p</sup> as well as relaxation times ρp. The latter is given by

$$
\rho\_p\left(T,\alpha\right) = a\left(T,\alpha\right)\tau\_p,\tag{6.40}
$$

where a (T,α)represents the principle of time-temperature-superposition (TTS), which shifts the reference relaxation times τ<sup>p</sup> in case of a change in temperature or cure degree. Here, the approach used by Simon et al. [107] is used,

$$\log\_{10} a\left(T, \alpha\right) = \frac{d\_V}{T - T\_{\infty}} - \frac{d\_V}{T\_{\text{g}}\left(\alpha\right) - T\_{\infty}},\tag{6.41}$$

in which the parameter d<sup>V</sup> enables to fit the model to the actual material behavior. The cure dependency of the rubbery resin modulus G<sup>r</sup> (α) is described using the model proposed by Adolf and Martin [115] in the form given by Equation 6.18. The same approach is also used for modeling the bulk modulus which links volumetric stress and strain:

$$\begin{aligned} K\left(t, T, \alpha\right) &= \quad K\_{\mathbf{r}}\left(\alpha\right) + \left[K\_{\mathbf{g}} - K\_{\mathbf{r}}\left(\alpha\right)\right] \sum\_{p=1}^{n} k\_{p} \exp\left(-\frac{t}{\rho\_{p}\left(T, \alpha\right)}\right) \\ &= \quad K\_{\mathbf{r}}\left(\alpha\right) + \Delta K\left(\alpha\right) \sum\_{p=1}^{n} k\_{p} \exp\left(-\frac{t}{\rho\_{p}\left(T, \alpha\right)}\right). \end{aligned} \tag{6.42}$$

The temporal evolution of deviatoric and volumetric stress is calculated by numerically solving time integrals given by

$$
\sigma\_{\rm dev}\left(t\right) = \int\_0^t G\left(t-\xi, T, \alpha\right) \frac{\mathrm{d}\hat{\varepsilon}\_{\rm dev}\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi \tag{6.43}
$$

and

$$
\sigma\_{\rm vol} \left( t \right) = \int\_0^t K \left( t - \xi, T, \alpha \right) \frac{\mathrm{d}\hat{\varepsilon}\_{\rm vol} \left( \xi \right)}{\mathrm{d}\xi} \mathrm{d}\xi. \tag{6.44}
$$

Prior to implementation of these equations, they need to be transformed into a suitable formulation. This procedure is shown in detail here for the deviatoric stress given by Equation 6.43 but is also applied the same way to the volumetric stress. Inserting Equation 6.39 into Equation 6.43 yields

$$\sigma\_{\rm dev}\left(t\right) = \int\_{0}^{t} \left[ G\_{\rm r}\left(\alpha\right) + \sum\_{p=1}^{n} \Delta G\left(\alpha\right) g\_{p} \exp\left(-\frac{t-\xi}{\rho\_{p}\left(T,\alpha\right)}\right) \right] \frac{\mathrm{d}\hat{\varepsilon}\_{\rm dev}\left(\xi\right)}{\mathrm{d}\xi} d\xi. \tag{6.45}$$

This equation is first decomposed in individual components for sake of simplicity:

$$
\sigma\_{\rm dev} \left( t \right) \quad = \quad \sigma\_{\rm dev}^{\rm el} \left( t \right) + \sum\_{p=1}^{n} \sigma\_{\rm dev,p}^{\rm ve} \left( t \right) \tag{6.46}
$$

$$\sigma\_{\rm dev}^{\rm el}\left(t\right) \;=\int\_{0}^{t} G\_{\rm r}\left(\alpha\right) \frac{\mathrm{d}\hat{\varepsilon}\_{\rm dev}\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi \tag{6.47}$$

$$\sigma\_{\rm dev,p}^{\rm ve}(t) = \int\_0^t \Delta G\left(\alpha\right) g\_p \exp\left(-\frac{t-\xi}{\rho\_p\left(T,\alpha\right)}\right) \frac{\mathrm{d}\hat{\varepsilon}\_{\rm dev}\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi. \tag{6.48}$$

The stress σel dev can be interpreted as the elastic part while σve dev,p represents stresses that arise from viscoelastic properties of the resin and depend on strain rate and history. For implementation in form of an Abaqus compatible user material subroutine (UMAT), Equations 6.46-6.48 need to be transformed into an incremental formulation, which relates the stress at some point in time t<sup>m</sup> to that of the previous point in time <sup>t</sup><sup>m</sup>−<sup>1</sup>. For <sup>σ</sup>el dev we can write

$$\sigma\_{\rm dev}^{\rm el} \left( t\_m \right) = \int\_0^{t\_{m-1}} G\_\mathbf{r} \left( \alpha \right) \frac{\mathrm{d}\hat{\varepsilon}\_{\rm dev} \left( \xi \right)}{\mathrm{d}\xi} \mathrm{d}\xi + \int\_{t\_{m-1}}^{t\_m} G\_\mathbf{r} \left( \alpha \right) \frac{\mathrm{d}\hat{\varepsilon}\_{\rm dev} \left( \xi \right)}{\mathrm{d}\xi} \mathrm{d}\xi \quad (6.49)$$

where the first term on the right hand side equals the stress state at the previous time σel dev (t<sup>m</sup>−<sup>1</sup>), thus

$$
\sigma\_{\rm dev}^{\rm el} \left( t\_m \right) = \sigma\_{\rm dev}^{\rm el} \left( t\_{m-1} \right) + \int\_{t\_{m-1}}^{t\_m} G\_{\rm r} \left( \alpha \right) \frac{\mathrm{d} \hat{\varepsilon}\_{\rm dev} \left( \xi \right)}{\mathrm{d} \xi} \mathrm{d} \xi. \tag{6.50}
$$

The remaining integral term, which represents the stresses which are to be added within the current time step, is approximated by assuming a linear variation in time of the modulus G<sup>r</sup> (α) as well as the deviatoric strain εˆdev (ξ). Hence, Equation 6.50 is rewritten as

$$
\sigma\_{\rm dev}^{\rm el} \left( t\_m \right) = \sigma\_{\rm dev}^{\rm el} \left( t\_{m-1} \right) + G\_{\rm r} \left( \overline{\alpha} \right) \cdot \Delta \hat{\varepsilon}\_{\rm dev} \left( t\_m \right) \tag{6.51}
$$

in which α refers to the midpoint value of the cure degree within the time step Δt = t<sup>m</sup> − t<sup>m</sup>−<sup>1</sup>, given by

$$
\overline{\alpha} = \frac{1}{2} \left[ \alpha\_m - \alpha\_{m-1} \right]. \tag{6.52}
$$

Equation 6.51 shares great similarity with that of the incremental stress evolution during cure of the representative model discussed in Section 6.3.1 (cf. Equation 6.29). Hence, the effective shrinkage strain is correctly taken into account by this model.

The viscoelastic stress of the individual Prony elements σve dev,p is derived in a similar way. Splitting the integral yields

$$\begin{split} \sigma\_{\text{dev},p}^{\text{ve}}(t\_m) &= \int\_0^{t\_m-1} \Delta G(\alpha) \, g\_p \exp\left(-\frac{t\_m-\xi}{\rho\_p(T,\alpha)}\right) \frac{\text{d}\dot{\varepsilon}\_{\text{dev}}(\xi)}{\text{d}\xi} \text{d}\xi \\ &+ \int\_{t\_m-1}^{t\_m} \Delta G(\alpha) \, g\_p \exp\left(-\frac{t\_m-\xi}{\rho\_p(T,\alpha)}\right) \frac{\text{d}\dot{\varepsilon}\_{\text{dev}}(\xi)}{\text{d}\xi} \text{d}\xi \end{split} (6.53)$$

in which the first term on the right hand side is very similar to the history term of σve dev,p:

$$\sigma\_{\rm dev,p}^{\rm ve}\left(t\_{m-1}\right) = \int\_0^{t\_{m-1}} \Delta G\left(\alpha\right) g\_p \exp\left(-\frac{t\_{m-1} - \xi}{\rho\_p\left(T, \alpha\right)}\right) \frac{\mathrm{d}\varepsilon\_{\rm dev}\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi. \tag{6.54}$$

The only difference is found within the exponential term, where t<sup>m</sup>−<sup>1</sup> rather than t<sup>m</sup> must be used. This is solved by multiplying with an exponential term which can then be written as

$$\begin{split} \sigma\_{\mathrm{dev},p}^{\mathrm{ve}}(t\_m) &= \exp\left(-\frac{\Delta t}{\rho\_p(T,\alpha)}\right) \sigma\_{\mathrm{dev},p}^{\mathrm{ve}}(t\_{m-1}) \\ &+ \int\_{t\_m-1}^{t\_m} \Delta G(\alpha) \, g\_p \exp\left(-\frac{t\_m-\xi}{\rho\_p(T,\alpha)}\right) \frac{\mathrm{d}\hat{\varepsilon}\_{\mathrm{dev}}(\xi)}{\mathrm{d}\xi} \mathrm{d}\xi. \end{split} \tag{6.55}$$

Assuming linear variation in time of strain, temperature and cure degree, the remaining integral term can be simplified to

$$\int\_{t\_m - 1}^{t\_m} \Delta G\left(\alpha\right) g\_p \exp\left(-\frac{t\_m - \xi}{\rho\_p\left(T, \alpha\right)}\right) \frac{\mathrm{d}\hat{\varepsilon}\_{\mathrm{dev}}\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi \qquad (6.56)$$

$$= -\Delta G\left(\overline{\alpha}\right) g\_p \frac{\Delta \hat{\varepsilon}\_{\mathrm{dev}}\left(t\_m\right)}{\Delta t} \int\_{t\_m - 1}^{t\_m} \exp\left(-\frac{t\_m - \xi}{\rho\_p\left(\overline{T}, \overline{\alpha}\right)}\right) \mathrm{d}\xi .$$

Solving the integral finally yields

$$\begin{split} \sigma\_{\text{dev},p}^{\text{ve}}\left(t\_{m}\right) &= \ \exp\left(-\frac{\Delta t}{\rho\_{p}\left(T,\alpha\right)}\right)\sigma\_{\text{dev},p}^{\text{ve}}\left(t\_{m-1}\right) \\ &+ \quad \Delta G\left(\overline{\alpha}\right)g\_{p}\frac{\Delta \hat{\varepsilon}\_{\text{dev}}\left(t\_{m}\right)}{\Delta t}\rho\_{p}\left(\overline{T},\overline{\alpha}\right)\left[1-\exp\left(-\frac{\Delta t}{\rho\_{p}\left(\overline{T},\overline{\alpha}\right)}\right)\right]. \end{split} \tag{6.57}$$

Similar formulations can be derived for the volumetric stress components, namely

$$
\sigma\_{\rm vol}^{\rm el} \left( t\_m \right) \quad = \quad \sigma\_{\rm vol}^{\rm el} \left( t\_{m-1} \right) + K\_{\rm r} \left( \overline{\alpha} \right) \cdot \Delta \hat{\varepsilon}\_{\rm vol} \left( t\_m \right) \tag{6.58}
$$

$$
\sigma\_{\rm vol,p}^{\rm ve}(t\_m) \quad = \quad \exp\left(-\frac{\Delta t}{\rho\_p(T,\alpha)}\right) \sigma\_{\rm vol,p}^{\rm ve}(t\_{m-1}) \tag{6.59}
$$

$$+\quad\Delta K\left(\overline{\alpha}\right)k\_{p}\frac{\Delta \hat{\varepsilon}\_{\text{vol}}\left(t\_{m}\right)}{\Delta t}\rho\_{p}\left(\overline{T},\overline{\alpha}\right)\left[1-\exp\left(-\frac{\Delta t}{\rho\_{p}\left(\overline{T},\overline{\alpha}\right)}\right)\right].$$

The strain increment Δε, from which deviatoric and volumetric strain increments are derived, consists of contributions from external forces, chemical shrinkage as well as thermal expansion or contraction, thus

$$
\Delta \varepsilon\_{\!\!\!\!\/} = \ \Delta \varepsilon\_{\!\!\!\/ \rm ext} - \Delta \varepsilon\_{\!\!\!\/\rm \rm th} \tag{6.60}
$$

$$
\Delta \varepsilon\_{\text{th}} = \beta\_{\text{th}} \cdot \Delta T \cdot I \tag{6.61}
$$

$$
\Delta \varepsilon\_{\text{ch}} = \beta\_{\text{ch}} \cdot \Delta \alpha \cdot I \tag{6.62}
$$

where βth, βch and I denote the coefficient of thermal expansion, the coefficient of chemical shrinkage and the identity matrix, respectively.

## **6.3.4 Transversely isotropic viscoelastic model**

The isotropic material model described in the previous section is able to model the material behavior of neat resin. Composites, however, per definition consist of both fiber and matrix and thus exhibit anisotropic material behavior. In the most general case, this involves 81 independent material parameters. Fortunately, this large amount of parameters can be greatly reduced by taking into account symmetries of the material as well as stress and strain tensors, as shown in detail by Barbero [121]. In this thesis, a transversely isotropic model is used for simulations on macroscopic component level which enables to reduce the number of independent parameters to five.

In order to avoid lengthy indexation, Nye notation is used in the following since it is the notation that is also used in Abaqus UMAT programming. The relationship between tensor and Nye notation is given by

$$\begin{bmatrix} \varepsilon\_1\\ \varepsilon\_2\\ \varepsilon\_3\\ \varepsilon\_4\\ \varepsilon\_5\\ \varepsilon\_6 \end{bmatrix} := \begin{bmatrix} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \varepsilon\_{33} \\ \varepsilon\_{12} \\ \varepsilon\_{13} \\ \varepsilon\_{23} \end{bmatrix}. \tag{6.63}$$

With this, the stress-strain relation of a transversely isotropic material is given by

$$
\begin{bmatrix}
\sigma\_1\\ \sigma\_2\\ \sigma\_3\\ \sigma\_4\\ \sigma\_5\\ \sigma\_6
\end{bmatrix} = \begin{bmatrix}
C\_{11} & C\_{12} & C\_{12} & 0 & 0 & 0\\ & C\_{22} & C\_{23} & 0 & 0 & 0\\ & & C\_{22} & 0 & 0 & 0\\ & & & C\_{44} & 0 & 0\\ & & & & C\_{44} & 0\\ & & & & & C\_{66}
\end{bmatrix} \cdot \begin{bmatrix}
\varepsilon\_1\\ \varepsilon\_2\\ \varepsilon\_3\\ \varepsilon\_4\\ \varepsilon\_5\\ \varepsilon\_6
\end{bmatrix} \tag{6.64}
$$

in which C<sup>66</sup> is defined as

$$C\_{66} = \frac{C\_{22} - C\_{23}}{2}.\tag{6.65}$$

In comparison to the isotropic material model described in the previous section, the transversely isotropic model is formulated without splitting of the strain in deviatoric and volumetric parts. Instead, each independent material property is modeled viscoelastic using a Prony series:

$$\begin{aligned} C\_{ij}\left(t, T, \alpha\right) &= \quad C\_{\mathbf{r}, ij}\left(\alpha\right) + \left[C\_{\mathbf{g}, ij} - C\_{\mathbf{r}, ij}\left(\alpha\right)\right] \sum\_{p=1}^{n} c\_p \exp\left(-\frac{t}{\rho\_p\left(T, \alpha\right)}\right) \\ &= \quad C\_{\mathbf{r}, ij}\left(\alpha\right) + \Delta C\_{ij}\left(\alpha\right) \sum\_{p=1}^{n} c\_p \exp\left(-\frac{t}{\rho\_p\left(T, \alpha\right)}\right) . \end{aligned} (6.66)$$

This approach avoids cumbersome and extensive modeling of underlying engineering constants and their viscoelastic characteristic. Moreover, it facilitates homogenization of UD properties based on RVE simulations as only three load cases need to be applied in order to identify the five independent material parameters [122].

Derivation of the constitutive equations follows a similar procedure as presented in the previous section for isotropic material behavior. Hence, the temporal evolution of stress is given in integral form by

$$
\sigma\_i(t) = \int\_0^t C\_{ij} \left( t - \xi \right) \frac{\mathrm{d}\varepsilon\_j \left( \xi \right)}{\mathrm{d}\xi} \mathrm{d}\xi \tag{6.67}
$$

in which <sup>C</sup>ij denotes the individual components of the stiffness matrix *<sup>C</sup>* defined in Equation 6.64. Inserting Equation 6.66 into Equation 6.67 and splitting into elastic and viscoelastic components for sake of clarity yields

$$\sigma\_i(t) \quad = \ \sigma\_i^{\text{el}}(t) + \sum\_{p=1}^n \sigma\_{i,p}^{\text{ve}}(t) \tag{6.68}$$

$$
\sigma\_i^{\text{el}}(t) = \int\_0^t C\_{\text{r},ij} \left( \alpha \right) \frac{\text{d}\varepsilon\_j \left( \xi \right)}{\text{d}\xi} \text{d}\xi \tag{6.69}
$$

$$\sigma\_{i,p}^{\rm w}(t) = \int\_0^t \Delta C\_{ij}\left(\alpha\right) c\_p \exp\left(-\frac{t-\xi}{\rho\_p\left(T,\alpha\right)}\right) \frac{\mathrm{d}\varepsilon\_j\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi \quad (6.70)$$

The cure dependency of the individual components of the rubbery stiffness matrix Cr,ij (α) is modeled by applying the same approach as used in the isotropic viscoelastic case, given by

$$C\_{\rm r,ij} \left( \alpha \right) = C\_{\rm r, \infty, ij} \cdot \left( \frac{\alpha^2 - \alpha\_{\rm gel}^2}{1 - \alpha\_{\rm gel}^2} \right)^{\kappa}, \tag{6.71}$$

in which Cr,∞,ij denotes the ultimate value of the stiffness Cr,ij when the resin is in fully cured state.

By applying the same procedure as in the isotropic case, the elastic and viscoelastic stress at time t<sup>m</sup> is found to be

$$
\sigma\_i^{\rm el} \left( t\_m \right) = \sigma\_i^{\rm el} \left( t\_{m-1} \right) + C\_{\rm r, ij} \left( \overline{\alpha} \right) \cdot \Delta \varepsilon\_j \left( t\_m \right) \tag{6.72}
$$

and

$$
\sigma\_{i,p}^{\text{ve}}\left(t\_m\right) = \exp\left(-\frac{\Delta t}{\rho\_p\left(T,\alpha\right)}\right)\sigma\_{i,p}^{\text{ve}}\left(t\_{m-1}\right) \tag{6.73}
$$

$$
+ \quad \Delta C\_{ij}\left(\overline{\alpha}\right)c\_p \frac{\Delta \varepsilon\_j\left(t\_m\right)}{\Delta t}\rho\_p\left(\overline{T},\overline{\alpha}\right)\left[1-\exp\left(-\frac{\Delta t}{\rho\_p\left(\overline{T},\overline{\alpha}\right)}\right)\right].
$$

The anisotropy of the material also influences chemical and thermal strains which become anisotropic too. This needs to be taken into account when computing the strain increment Δε, which in this case is defined as

$$
\Delta \varepsilon\_{\!\!\!\!\/} = \ \Delta \varepsilon\_{\!\!\!\/ \rm ext} - \Delta \varepsilon\_{\!\!\/\rm th} \tag{6.74}
$$

$$
\Delta \varepsilon\_{\text{th}} = \beta\_{\text{th}} \cdot \Delta T \tag{6.75}
$$

$$
\Delta \varepsilon\_{\text{ch}} = \beta\_{\text{ch}} \cdot \Delta \alpha. \tag{6.76}
$$

Here, βth and βch represent tensors given by

$$
\beta\_{\rm th} = \begin{bmatrix}
\beta\_{\rm th, 11} & 0 & 0 \\
& \beta\_{\rm th, 22} & 0 \\
\text{sym.} & & \beta\_{\rm th, 22}
\end{bmatrix} \tag{6.77}
$$

and

$$
\beta\_{\rm ch} = \begin{bmatrix}
\beta\_{\rm ch,11} & 0 & 0 \\
& \beta\_{\rm ch,22} & 0 \\
\rm sym. & & \beta\_{\rm ch,22}
\end{bmatrix}.
\tag{6.78}
$$

Due to the assumption of transversely isotropic material behavior of unidirectional CFRP, both transverse directions contribute with the same amount of chemical and thermal strain to the overall strain increment. Thus, it is valid to set the out-of-plane component to the value of the in-plane transverse component in Equations 6.77 and 6.78, respectively.

The presented transversely isotropic viscoelastic material model demands for a significant amount of material data. Similar to Kwok and Pellegrino [122], the relaxation times τ<sup>p</sup> are assumed to be equal to that of the isotropic model. This is motivated by the fact that the viscoelastic behavior is caused solely by the matrix material and thus has the same effect on stress relaxation as in the isotropic case. Contrary to the approach of Kwok and Pellegrino [122], the weighting factors of the Prony series' are also assumed to be equal to that of the neat resin. This is motivated by a difference in model formulation of the components of the stiffness matrix. Kwok and Pellegrino [122] use a model given by

$$C\_{ij} = C\_{\rm r,ij} + \sum\_{p=1}^{n} C\_{ij,p} \exp\left(-\frac{t}{\rho\_p}\right) \tag{6.79}$$

in which weighting factors and increase in modulus due to vitrification are merged in a series of moduli Cij,p directly related to the individual Prony elements. The approach shown above makes use of a similar relationship that isolates the Prony coefficients from the glassy modulus:

$$C\_{ij} = C\_{\rm r,ij} + [C\_{\rm g,ij} - C\_{\rm r,ij}] \sum\_{p=1}^{n} c\_p \exp\left(-\frac{t}{\rho\_p}\right). \tag{6.80}$$

This reduces the number of parameters which need to be identified during homogenization to ten, namely rubbery moduli Cr,11, Cr,22, Cr,12, Cr,23, Cr,<sup>44</sup> and glassy moduli Cg,11, Cg,22, Cg,12, Cg,23, Cg,44.

### **6.3.5 Transversely isotropic path-dependent model**

The viscoelastic models presented above enable accurate modeling of the material behavior but require an extensive effort in material characterization. Furthermore, their demand in terms of computation time and memory is much higher compared to elastic models. It would therefore be favorable to use elastic models which are setup such that they mimic the viscoelastic behavior. The motivation for this has been presented in Section 6.3.2 and is based on the fact that the shrinkage strain evolves very slowly during curing. Moreover, cross-linking usually takes place within several minutes at elevated temperatures near to or above Tg. This leads to the conclusion that strain rate due to shrinkage and relaxation times are very small. Hence, the contribution of the viscous part of the Prony series is negligibly small and the resin therefore exhibits elastic rubbery material behavior. The subsequent cooling process, however, has the opposite effect. Due to the decrease in temperature, relaxation times increase rapidly, especially if T<sup>g</sup> is undercut. Consequently, viscous properties increasingly participate in stress development and the material exhibits glassy behavior during a large portion of the cooling process. With the assumption of a fast switch between rubbery and glassy state, we can define a material state γ which controls this switching based on the acting temperature and cure degree:

$$\gamma\left(T,\alpha\right) = \begin{cases} 0, & T\_{\rm g}\left(\alpha\right) - T \le \Delta T\_{\rm vit} \quad \text{(rubberg)}\\ 1, & T\_{\rm g}\left(\alpha\right) - T > \Delta T\_{\rm vit} \quad \text{(glassy)}. \end{cases} \tag{6.81}$$

Here, ΔTvit refers to the difference in T<sup>g</sup> and T, at which the material vitrifies. While this is a material specific quantity, it can also be used to adjust the model and to compensate for uncertainties of measurements but also models, e.g. kinetic models. By introducing a switching variable γ, the damper in the representative *Standard linear solid* model used above (cf. Figure 6.9) effectively is exchanged with a digital switch as indicated in Figure 6.10.

**Figure 6.10:** Schematic representation of a path-dependent model

If the material is in glassy state (γ = 1), the switch is closed and the viscoelastic branch contributes to the overall stress state. Devitrification causes opening of the switch and a subsequent loss of the stress in that branch.

The stress increment is derived in a similar manner as in the viscoelastic case. However, since it is not intended that the path-dependent model is able to take into account any relaxation spectrum, the Prony series of the viscoelastic component can be limited to one element, which is enough to realize the hard switching between rubbery and glassy state as will be shown in the following. The stress increment is again divided in elastic and viscoelastic contributions:

$$
\sigma\_i(t) = \sigma\_i^{\text{el}}(t) + \sigma\_i^{\text{ve}}(t) \tag{6.82}
$$

$$
\sigma\_i^{\rm el}(t) = \int\_0^t C\_{\rm r,ij} \left( \alpha \right) \frac{\mathrm{d}\varepsilon\_j \left( \xi \right)}{\mathrm{d}\xi} \mathrm{d}\xi \tag{6.83}
$$

$$\sigma\_i^{\rm ve}(t) = \int\_0^t \Delta C\_{ij}\left(\alpha\right) \exp\left(-\frac{t-\xi}{\rho\left(T,\alpha\right)}\right) \frac{\mathrm{d}\varepsilon\_j\left(\xi\right)}{\mathrm{d}\xi} \mathrm{d}\xi. \tag{6.84}$$

The elastic stress remains unchanged compared to Equation 6.69 of the fully viscoelastic approach, thus

$$
\sigma\_i^{\rm el} \left( t\_m \right) = \sigma\_i^{\rm el} \left( t\_{m-1} \right) + C\_{\rm r, ij} \left( \overline{\alpha} \right) \cdot \Delta \varepsilon\_j \left( t\_m \right) \,. \tag{6.85}
$$

Since we only differ between fully rubbery and fully glassy state, the viscoelastic stress is derived for two cases of ρ:

$$\rho\left(T,\alpha\right) = \begin{cases} 0, & T\_{\rm g}\left(\alpha\right) - T \le \Delta T\_{\rm vit} \quad \text{(rubberg)}\\ \infty, & T\_{\rm g}\left(\alpha\right) - T > \Delta T\_{\rm vit} \quad \text{(glassy)}. \end{cases} \tag{6.86}$$

Applying these two limiting cases to Equation 6.84 yields

$$\lim\_{\rho \to 0} \sigma\_i^{\text{ve}}(t) = 0 \tag{6.87}$$

and

$$\lim\_{\rho \to \infty} \sigma\_i^{\rm ve}(t) = \int\_0^t \Delta C\_{ij} \left( \alpha \right) \frac{\mathrm{d} \varepsilon\_j \left( \xi \right)}{\mathrm{d} \xi} \mathrm{d} \xi. \tag{6.88}$$

143

Integrating Equation 6.88 in the same manner as demonstrated above in the derivation of the fully viscoelastic model results in

$$\sigma\_i^{\rm ve}(t\_m) = \begin{cases} 0, & T\_{\rm g}\left(\alpha\right) - T \le \Delta T\_{\rm vit} \\ \sigma\_i^{\rm ve}\left(t\_{m-1}\right) + \Delta C\_{ij}\left(\overline{\alpha}\right) \cdot \Delta \varepsilon\_j\left(t\_m\right), & T\_{\rm g}\left(\alpha\right) - T > \Delta T\_{\rm vit} \end{cases} \tag{6.89}$$

The cure dependency of the rubbery stiffness Cr,ij (α) is described using the same model as in the viscoelastic approach (cf. Equation 6.71).

## **6.3.6 Transversely isotropic linear elastic CHILE model**

The ability of path-dependent models to allow for devitrification of the material and, thus, the release of any frozen-in stress, is not needed for simulating a process that does not contain a reheating near to or above Tg. In this case, the material model can be further simplified by incorporating the switching between glassy and rubbery material behavior into the stiffness of the elastic branch. This enables to omit the viscoelastic branch completely and is shown schematically in Figure 6.11.

**Figure 6.11:** Schematic representation of a CHILE model

The temperature and cure dependent stiffness is modeled based on the *Cure hardening instantaneously linear elastic* (CHILE) approach proposed by Johnston [59] given by Equation 6.13. The model is reformulated such that it describes the transition between glassy and rubbery value of each component Cij

of the stiffness tensor. Furthermore, the rubbery modulus is modeled with cure dependency, yielding

$$C\_{ij}\left(T,\alpha\right) = \begin{cases} C\_{\rm r,ij}\left(\alpha\right), & T\_{\rm g}\left(\alpha\right) - T < \Delta T\_{\rm c1} \\ C\_{\rm r,ij}\left(\alpha\right) + \Phi\left(T,\alpha\right)\Delta C\_{ij}\left(\alpha\right), & \Delta T\_{\rm c1} < T\_{\rm g}\left(\alpha\right) - T < \Delta T\_{\rm c2} \\ C\_{\rm g,ij}, & T\_{\rm g}\left(\alpha\right) - T > \Delta T\_{\rm c2} \end{cases} \tag{6.90}$$

with the transition function

$$\Phi\left(T,\alpha\right) = \frac{T\_{\rm g}\left(\alpha\right) - T - \Delta T\_{\rm c1}}{\Delta T\_{\rm c2} - \Delta T\_{\rm c1}}\tag{6.91}$$

and

$$
\Delta C\_{ij} \left( \alpha \right) = C\_{\rm g,ij} - C\_{\rm r,ij} \left( \alpha \right) \,. \tag{6.92}
$$

In this model, the material behavior transitions from glassy to rubbery state, and vice versa, depending on the difference between T<sup>g</sup> and temperature T. The transition is controlled and fitted to the individual material behavior using the model parameters ΔTc1 and ΔTc2.

This approach is applied to all components of the stiffness tensor with the exception of C<sup>11</sup> as the stiffness in fiber direction is dominated by the mechanical properties of the fiber. Hence, the cure dependency of the rubbery modulus Cr,<sup>11</sup> is negligible.

The stress at time t<sup>m</sup> is given by a relationship similar to that of the elastic part in the viscoelastic and path-dependent models:

$$
\sigma\_i\left(t\_m\right) = \sigma\_i\left(t\_{m-1}\right) + C\_{ij}\left(\overline{T}, \overline{\alpha}\right) \cdot \Delta\varepsilon\_j\left(t\_m\right). \tag{6.93}
$$

### **6.3.7 Determination of the point of vitrification**

In the above discussed path-dependent and elastic chemo-thermo-mechanical models, switching between rubbery and glassy state is controlled by the model parameter ΔTvit. Hence, what is modeled frequency- and temperature-dependent in a viscoelastic model, is represented by a single model parameter, representing the viscoelastic behavior of the material. Figure 6.12 shows the influence of frequency on the modulus evolution of resin SCER1 during curing. The transition of the material from initially rubbery state to glassy state is shifted to earlier times as the frequency increases. This behavior is typical for a viscoelastic material and is caused by the rate or frequency dependency of the viscous property of the material. In the representative *Standard linear solid* model used above, this is reflected in the rate dependency of the damper. Increasing the frequency is followed by an increased participation of the spring of the viscoelastic branch in the overall stress development. At the same time, the relaxation time increases due to the fact that T<sup>g</sup> approaches the cure temperature, leading to an increase of the viscosity of the damper. As a result, glass transition is either shifted to earlier or later times, depending on the acting frequency.

**Figure 6.12:** Modulus evolution during isothermal curing of resin SCER1 using a cure temperature of 140 ◦C. After 43 min of isothermal cure, the specimen is cooled using a rate of 5 K min−1. The storage modulus has been measured using an oscillatory rheometer and multiwave technique.

The influence of frequency on modulus evolution is setup in the path-dependent model using the parameter ΔTvit. Increasing ΔTvit prolongs the glass transition which is similar to the behavior of the viscoelastic material at a reduced frequency. In the CHILE model, the same can be achieved by varying parameters ΔTc1 and ΔTc2, which define onset and endset of glass transition.

The effect of shifting glass transition to later times, e.g. in case of a very small frequency or strain rate, ultimately leads to a case in which no vitrification is observed during the cure step. Instead, the material vitrifies during the subsequent cooling step. While this behavior can be realistic, e.g. in case of cooling from a curing temperature well above the ultimate glass transition temperature, it will in most cases raise the question of what frequency should be assumed in case of a common injection process like RTM. In order to estimate this frequency range, a simple analogy between a fixed frequency rheometer experiment and the temporal evolution of the shrinkage strain is drawn.

The shrinkage strain increases monotonically during curing. Thus, we can try to mimic the curve from the point of gelation on and until the reaction rate vanishes with the help of shifted sinus functions. Figure 6.13a shows the shrinkage strain of Hexcel 8551-7 resin [107] during curing at an isothermal temperature of 163 ◦C, which is also the value of Tg,<sup>∞</sup> at the end of the cure step. The figure also contains two sinusoidal functions of different frequencies which represent a lower and an upper bound of frequency. While the higher frequency is a better representation of the fast advancement in cure right after gelation, the lower frequency is better suited to describe the low strain rate at the end of the curing process. The important fact here is that both frequencies are magnitudes lower than those commonly used in material characterization. As is shown in Figure 6.13b, this has dramatic influence on the stiffness evolution. Both identified frequencies prolong glass transition such that it will only be visible during subsequent cooling. The often applied frequency of 1 Hz, however, leads to a totally different curve, showing an early increase of the modulus which eventually reaches its glassy value. If this modulus development is used in an elastic model, the stiffness will be overestimated by a factor of 100.

**Figure 6.13:** Approximation of the loading frequency during an RTM process based on Hexcel 8551-7 resin model and data of Simon et al. [107]

The lower and upper bound of frequency lie in a range of 1 <sup>×</sup> <sup>10</sup>−5 Hz to <sup>1</sup> <sup>×</sup> 10−4 Hz, which is similar to that reported by Zobeiry et al. [113]. However, these values depend on the reactivity of the resin. As is implied by the limits of the X-axis of Figure 6.13, this resin obviously belongs to the class of slow-curing resins. In this theoretical assessment, this can easily be changed by modifying the parameters of the reaction kinetic model such that it transforms into a faster curing resin while all other properties, i.e. relaxation times, remain unchanged. The result is shown in Figure 6.14. Due to the increased reaction rate, the shrinkage strain develops much faster. This leads to significantly higher frequencies identified for both lower and upper bound. As can be seen from the final stiffness values, this indeed leads to a slight increase in resin modulus caused by a partial vitrification, which implies that the increased frequency introduces viscous stresses that persist until the end of the curing stage. However, this increase is still much lower than that observed for the often applied frequency of 1 Hz.

The shown predicted evolution of shrinkage strain and modulus are strictly valid only for the resin for which the actual measurements have been conducted. However, the qualitative behavior should be to some degree similar across different

**Figure 6.14:** Approximation of the loading frequency during an RTM process based on Hexcel 8551-7 resin model and data of Simon et al. [107]. Kinetic model parameters are modified to enable faster curing.

resins with differences in the effect of an increased strain rate during curing, which is an inherent nature of fast-curing resins. But there are still two major uncertainties that hinder derivation of a universal process frequency. Firstly, measuring the resin's modulus at a frequency as low as 1 <sup>×</sup> <sup>10</sup>−3 Hz during curing is impossible since this would mean that a data point is acquired every <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>3</sup> s. Thus, the material would completely cure before even one measurement point has been recorded. Secondly, shrinkage strain is not the only source of process-induced strains. Thermal strains also evolve during subsequent cooling. However, regarding the actual strain rate, the speed and sudden of demolding may have a much higher impact on the acting stiffness of the material. Furthermore, demolding happens at the same temperature as curing. Hence, depending on the sudden of demolding and the subsequently applied cooling rate, viscous stresses arising during demolding may or may not relax before the temperature has been decreased to a level, at which further stress relaxation in a relevant period of time becomes negligibly small. This can make it very difficult to properly calibrate elastic models. Thus, in case of very complex demolding and processing conditions, the use of a viscoelastic approach may be more appropriate or even

mandatory. In Section 7.3.6, this is further evaluated by means of a parameter study of process-induced distortions.

The coefficient of thermal expansion (CTE) is another property that changes significantly when the resin transitions from glassy to rubbery state or vice versa. In this case, however, the transition is caused by changes in free volume [28], and not by changes in frequency or strain rate. Hence, the transition can be modeled such that it switches between both states depending only on the actual temperature and Tg. This can be achieved using the material state γ introduced in Equation 6.81. With this, the coefficient of thermal expansion βth is defined as

$$
\beta\_{\rm th,i} = \beta\_{\rm th,r,i} + \gamma \left( T, \alpha \right) \left( \beta\_{\rm th,g,i} - \beta\_{\rm th,r,i} \right), \tag{6.94}
$$

in which subscripts r and g denote the CTE in rubbery and glassy state, respectively. In the viscoelastic model, the same approach can be applied. Alternatively, if the shift function is modeled using T<sup>g</sup> as the reference temperature (cf. Equation 6.41), it is possible to couple γ and βth directly to the shift function:

$$\gamma\left(T,\alpha\right) = \begin{cases} 0, & a\left(T,\alpha\right) \le 1 \quad \text{(rubberg)}\\ 1, & a\left(T,\alpha\right) > 1 \quad \text{(glassy)} \end{cases} \tag{6.95}$$

This approach makes use of an intrinsic property of viscoelasticity and thus avoids an additional model parameter (ΔTvit). Equation 6.95 is therefore used in this thesis for modeling the material state and thus thermal expansion in the viscoelastic models.

### **6.4 Discussion and concluding remarks**

In this chapter, elastic and viscoelastic chemo-thermo-mechanic models are derived. Furthermore, the viscoelastic behavior of the resin during the course of the manufacturing process is analyzed using a representative rheological model. During curing, the relaxation time is low as the cure temperature is above the evolving glass transition temperature. Moreover, the strain rate induced by shrinkage is low as well, since it develops over a relatively long period of time. Consequently, the material will be in the rubbery state for most of the time while curing. During demolding, the strain rate depends largely on the design of this process step. If the part is abruptly pushed out of the mold, the strain rate will be much higher compared to gradual demolding. Moreover, the spatial distribution of strain rate depends on part geometry and adhesion to the mold surface, but also on the positions of force application. If viscous stress is generated during demolding, the amount of possible relaxation depends on the cooling rate applied afterwards. A slow cooling rate will give the material more time to relax. Moreover, it will remain longer at a temperature close to Tg, which means low relaxation times. As is demonstrated using the representative model, relaxation leads to an increase in resulting process-induced strain, and thus warpage and spring-in of composite structures. Hence, the prediction accuracy of the simulation depends mostly on the accurate description of viscous stress, that possibly emerges during demolding, and relaxation of this stress during subsequent cooling. This can be achieved by using a viscoelastic material model, which takes the time-temperature-cure superposition into account. However, in comparison to elastic approaches like CHILE, the computational effort is much higher when using viscoelastic models. Hence, it is favorable to use more efficient elastic models instead, in which the influence of temperature, strain rate and curing is incorporated into the evolution of the elastic modulus. The latter is especially vital but also complex in the vicinity of glass transition due to the sensitivity to strain rate and frequency.

The requirements for proper modeling of the modulus for use with an elastic model need to be transferred to the experimental characterization. In this context, it is important to emphasize, that when measuring the evolution of the elastic modulus by using a specific constant frequency, the resulting data will only be valid for this specific frequency. This applies especially to any transition from rubbery to glassy state or vice versa and ultimately raises the question of what frequency acts during a typical RTM process. Zobeiry et al. [113] compared stress predictions of an elastic model using different frequencies with that of a

viscoelastic model. Good correlation between the models is found when using a frequency in the range of 0.001 Hz to 0.0001 Hz, but this may depend largely on the used resin due to differences in cure cycle and amount of chemical shrinkage. Johnston [59] gives a frequency range of 0.1 rad s−1 to 0.0001 rad s−1 for typical autoclave processing. However, using such a frequency for characterization of the material is infeasible. For example, a frequency of 0.0001 Hz would result in a time period of 1 <sup>×</sup> 104 s or 2.8 h, which is needed to acquire a single data point. Such time-consuming measurements are costly and at the same time make it impossible to characterize the evolution of the modulus during cure. Hence, characterization of the material using time-temperature-cure superposition, as needed for parametrization of viscoelastic models, is also beneficial for proper calibration of elastic models, as proposed by Zobeiry et al. [113]. When using this strategy, the experimental effort will be similar to what is needed for viscoelastic models, but the computational effort can be significantly reduced. Otherwise, the actual evolution of the modulus will be unknown and assumptions need to be made. Based on the results of the study presented in Section 6.3.7 a valid assumption may be that the strain rate or frequency of a typical RTM or WCM cure cycle is so low that most of the transition from rubbery to glassy state takes place during cooling. However, this assumption may quickly turn invalid as subsequent process steps like demolding or cooling also involve strain rates which might be much higher. This is investigated in the following chapter based on process-induced distortion of angled brackets.

# **7 Spring-in and Warpage of Composite Structures**

During composite manufacturing, process-induced residual stresses arise and consequently lead to process-induced distortion (PID). As with any other manufacturing process, the produced parts must satisfy certain geometrical tolerances. PID can make this difficult and sometimes impossible to achieve. Hence, numerical prediction of resulting distortions can be of great help in order to match tolerances and reduce geometrical deviations already during the design process of the composite structure.

In this chapter, spring-in of L-shaped specimens is investigated both experimentally and numerically. Moreover, distortion of a large automotive composite structure is analyzed and a comparison between the experimentally measured and numerically predicted displacement field is provided. For this purpose, the models introduced in Chapter 6 are applied and evaluated with respect to their individual prediction quality and possible limitations.

# **7.1 Introduction and review of related work**

Geometrical deviations of composite structures can be categorized in two main classes: spring-in and warpage. While the former is caused by a curved geometry, the latter can also occur in flat areas. Hence, the origin of each phenomenon differs. Figure 7.1 schematically shows the mechanism that leads to spring-in. Initially, the geometry of the shown curved section is in an undeformed state. In a closed-mold process like RTM, this geometry is defined by the cavity. At some

point during the course of the process, the parts surface demolds from the mold surface. This can happen during regular demolding or already before as soon as the force induced by shrinkage and thermal contraction reach a critical value, at which adhesion to the mold surface can no longer be sustained. In both cases, process-induced chemical and thermal strains lead to a physical contraction of the matrix and, albeit to a much smaller degree, of the fibers. As a result, the stiff fibers located at the inner surface of the curved section are forced on a bigger radius, which introduces tensile stress in this region. The opposite happens at the outer surface, introducing compressive stress. Consequently, the surface normal vector of the cutting plane at the end of the curved section is rotated inwards.

**Figure 7.1:** Schematic of the spring-in mechanism

In the above explanation of the spring-in effect, only strains in through-thickness direction of the laminate have been considered. However, during processing, residual strains will also develop in in-plane directions. If the layup is perfectly symmetric, the stresses induced by them will cancel out and no deformation is caused. If this is not the case, a stress gradient is induced, leading to out-ofplane deformations, commonly referred to as warpage. Hence, warpage occurs if gradients are present in through-thickness direction of a part and can therefore also superimpose spring-in in curved sections. Figure 7.2 exemplarily shows spring-in and warpage on specimen level.

**Figure 7.2:** Schematic of spring-in and warpage of L-shaped specimens (left) and initially flat composite stripes (right)

The effect of the spring-in mechanism depicted in Figure 7.1 on part or specimen level is illustrated on the left-hand side of Figure 7.2. The enclosed angle θ of the L-shaped part decreases and therefore deviates from the nominal angle θnom. From this, the temporal evolution of the spring-in angle is calculated by

$$
\Delta\theta\left(t\right) = \theta\left(t\right) - \theta\_{\text{nom}}.\tag{7.1}
$$

On the right-hand side of Figure 7.2, warpage of an initially flat composite stripe is shown. These specimens are often used for investigation of the warpage phenomenon using prepreg material, since this allow to cure the specimens without the need of a rigid mold and therefore enables to study the time dependent deformation as it develops during the manufacturing process [123–125]. However, due to the lack of compaction pressure, laminate quality with respect to pores might be impaired. Hence, in this thesis, emphasis is put on evaluating PID in terms of spring-in of parts manufactured using closed-mold processes like RTM. Moreover, automotive structures usually are manufactured using symmetric layups and thus are primarily affected by spring-in and only to a minor extent by warpage.

### **7.1.1 Sources of process-induced residual stresses**

As briefly outlined during material model derivation in Chapter 6, two sources of residual strains exist: chemical shrinkage and thermal expansion or contraction. The contribution of each of these phenomena depends on the individual material properties as well as on the chosen cure cycle. The same is true for potential relaxation of stresses introduced by these sources of residual strains during the course of the process. In addition, capturing residual stresses and strains experimentally is limited and difficult. Instead, evaluation of resulting stress levels is often done using numerical simulation. Based on a viscoelastic analysis of IM6/3100 prepreg material, White and Hahn [123] conclude that contribution of chemical shrinkage to residual stress is less than 4 %. Contrarily, Prasatya et al. [126] show that shrinkage strain contributes up to 30 % of the total residual stress. Plepys and Farris [127] investigated residual stresses due to chemical shrinkage that arise during curing of a three-dimensionally constrained epoxy resin. Due to the volumetric constrained curing, large residual stresses occur, which indicates that the resulting residual stress level may depend on the boundary conditions which act during curing. These are defined by the used molding concept but also by the layup as neighboring plies interact with each other. Hence, it makes more sense to quantify the influence of chemical and thermal residual strains with respect to their contribution to the resulting PID. Nelson and Cairns [128] report a chemical contribution to spring-in of Δθch = 30 % to 60 %. Yoon and Kim [129] investigated spring-in of L-shaped specimens manufactured using [±φ]2<sup>s</sup> layups of varying angles. Numerical analysis revealed a contribution of chemical shrinkage to spring-in of Δθch = 58 % for ply angles of φ = 0°, 15°, 30° and 45°. Interestingly, this is reduced to 47 % in case of a ply angle of φ = 60°, implying a dependence on layup configuration which was also reported by Radford and Rennick [130]. An even higher contribution of chemical strain to spring-in of L-shaped specimens is reported by Bapanapalli and Smith [125]. They investigated several cross-ply laminates of varying fractions of 0° and 90° plies and quantified the influence of chemical strains to 75 %. Comparability of the mentioned studies in general is limited since model approaches used for separating sources of spring-in vary from study to study. Moreover, different materials and thus also cure cycles are used which exhibit different amounts of chemical shrinkage and thermal expansion or contraction. Hence, it is questionable whether comparing residual stress or PID contributions across multiple studies is applicable. Moreover, the viscoelastic nature of the resin in the vicinity of T<sup>g</sup> introduces even more uncertainties. If, for instance, demolding happens suddenly but cooling is carried out using a very slow rate, most of the PID related to shrinkage will only gradually emerge during cooling of the part. Hence, if the deformed state after demolding and prior to cooling would be taken as basis for deriving the state of PID, this would lead to a significant underestimation.

#### Influence of residual stress on part properties

Depending on the magnitude of residual stress generated during the process, mechanical properties of the manufactured part may be impaired. Griffin [131] analyzed the three-dimensional residual stress state due to curing of symmetric cross-ply laminates and reported stress levels of 34 MPa to 50 MPa. A similar range of stress is found by Zhao et al. [132] using a thermo-viscoelastic material model. Moreover, they investigated interaction of residual stresses with damage initiation in unidirectional composites and found little influence on longitudinal failure and great impact on transverse and shear failure. Whether this influence is detrimental or beneficial to the damage initiation and evolution depends on the applied loading conditions. Moreover, the strength of the material may depend on achieved cure degree and mechanical and thermal loading conditions. White and Hahn [133] investigated the transverse strength of IM6/3100 prepreg material and found it to be strongly depending on cure degree. A transverse strength of 23 MPa is reported for the composite at nearly full cure. Hobbiebrunken et al. [134] analyzed the damage evolution in more detail by differing between yield and fracture. Both have been found to be strongly dependent on temperature. Yield stress of the used RTM6 epoxy resin is reduced from 27 MPa to 12 MPa (−54 %) when the temperature is increased from room temperature to the cure temperature of 180 ◦C. Similarly, the fracture stress is reduced from 90 MPa to 41 MPa (−55 %) for the same temperature delta. Furthermore, they report that yielding leads to a redistribution of residual stresses, effectively lowering their magnitude.

### **7.1.2 Part design-related factors influencing PID**

#### Layup and stacking sequence

The stacking sequence of asymmetric layups has significant influence on the resulting PID. This fact is used by various studies to provoke distortion in order to analyze acting residual stresses. Sarrazin et al. [135] investigated initially flat composite strips cured in a vacuum bagging process. Due to the use of asymmetric [0/90] layups, the specimens show distortion which evolves during the course of the process. It is found that the magnitude of distortion depends on how evenly distributed the individual plies are. Hence, increasing PID is reported when transitioning from [0/90]<sup>4</sup> to [02/902]<sup>2</sup> to [04/904]. When using symmetric layups, this effect is reported to be non-existent by Jain et al. [136]. This conclusion is drawn by comparing spring-in angles of specimens with layups of [90/0/90/0]<sup>s</sup> and [0/90/0/90]s. However, their results imply that the assumption of a symmetric layup no longer applies when using non-closed mold processes. This becomes apparent from a comparison of PID of the aforementioned layups manufactured using either male or female mold setup. Due to tool-part interaction, the resulting impact on spring-in shows contrary tendencies.

Another effect which is strongly affected by the stacking sequence is transverse strain locking. In a purely unidirectional layup, the high process-induced strains perpendicular to the fiber direction can develop freely. If plies of other orientations are introduced in this layup, transverse strain essentially is blocked. This is most pronounced in cross-ply layups, in which each adjacent plies exhibit a difference in orientation of 90°. Due to this transverse strain lock, these layups show much higher strains in through-thickness direction since this is the only direction in which the layup is still free to deform. Garstka et al. [119] analyzed throughthe-thickness strains evolving during processing of AS4/8552 prepreg and report through-the-thickness strains of cross-ply laminates twice as high as those of unidirectional layups. A similar conclusion is presented by Rogers et al. [137], who investigated the influence of ply orientations on thermal expansion. Since through-the-thickness strains have tremendous influence on spring-in of curved sections, this effect needs to be taken into account by numerical simulation and may render some modeling approaches like 2D plane strain inappropriate [59].

The effect of transverse strain locking on spring-in has been addressed by several authors. Groh et al. [138] investigated spring-in of L-shaped specimens manufactured using an RTM process. Their results show an increase in spring-in as more and more layers of the initial [(−45/45)2(90/0)]<sup>s</sup> layup are exchanged with (90/0) plies. Interestingly, the resulting [(90/0)3]<sup>s</sup> layup leads to more spring-in compared to a [(−45/45)3]<sup>s</sup> specimen. Since both laminates are crossply, the effective orientation of the layup relative to the curved section of the part has an effect on PID due to differences in stiffness in circumferential direction. This is supported by findings of Kappel [139] and Kappel et al. [140].

In general, specimens consisting of a purely unidirectional layup whose fibers are oriented along the extrusion direction of the curved section (perpendicular to the tangential direction) do not exhibit spring-in in case of RTM. However, some studies report non zero spring-in which can be traced back to either fabrics which are not entirely unidirectional, as used by Jain et al. [136], or the presence of tool-part interaction. The latter is often encountered in prepreg autoclave manufacturing and induces spring-in and warpage by the fact that only one side of the part is in contact with a rigid mold that undergoes thermal strains. Thus, tool-part interaction leads to stress gradients in through-thickness direction of the part, which results in distortion [139]. While no PID is expected for that specific layup, all mentioned studies investigating distortion of such specimens report significantly higher scatter. This may arise from the fact that no fibers are oriented tangential to the curved section. Hence, the material behavior in this direction is completely dominated by the matrix, which, due to the vicinity to T<sup>g</sup> by the time of demolding, is highly viscoelastic in nature. Thus, any bending of the part during demolding will introduce deformations that might not fully recover afterwards and before the temperature has decreased to a level at which no further relaxation is possible.

#### Part geometry

Jain et al. [136] investigated spring-in of L-shaped specimens of different combinations of radii and laminate thicknesses. They found spring-in to be independent of radius and thickness for radius-to-thickness ratios greater than one. Furthermore, they report smaller spring-in if the enclosed angle is increased. The latter is expected since by further increasing the enclosed angle, the part finally becomes flat and, thus, free of spring-in. However, isolating the effect of a change of a specific geometrical feature is often difficult since this also affects other properties of the part. A prominent example is fiber orientation, which depends on the part's geometry. Altering it possibly requires modifications to the draping strategy, influencing local fiber orientation. Other properties which are strongly related to draping of fiber textiles are fiber bridging and wrinkling in convex and concave sections. These can have significant effect on process-induced distortions due to the local formation of neat resin areas. Only limited research has been published on this topic so far [86, 89, 141].

### **7.1.3 Process-related factors influencing PID**

#### Cure cycle

The cure cycle influences the amount of residual strain in two ways. Firstly, the chosen cure temperature and time defines the resulting cure degree and therefore also the amount of chemical shrinkage. Secondly, the amount of thermal strain depends on the chosen cure temperature. Sarrazin et al. [135] and White and Hahn [124] conclude that by using a lower cure temperature, chemical and thermal strains and thus PID can be reduced. However, as noted by White and Hahn [124], this comes at the cost of impaired mechanical properties due to insufficient curing. The latter can be overcome by using an additional postcure step which ensures sufficient cross-linking, as investigated by Svanberg and Holmberg [5]. In this case, the evolution of residual strains also depends on the point in time, at which the change in cure temperature happens. If the temperature is increased after the resin has passed the point of gelation, thermal expansion will be induced that counteract shrinkage strains. Moreover, heating to the postcure temperature happens in rubbery state of the resin, exhibiting a much higher coefficient of thermal expansion compared to that of the glassy state. Hence, the expansive thermal strains that result from heating to the postcure temperature are much higher than the contracting strains that develop during the subsequent cooling to room temperature. Consequently, a reduction in PID is achived using such a cure cycle [5, 124, 135].

#### Cooling rate

Groh et al. [138] compared L-shaped specimens manufactured using RTM process and cooled afterwards using cooling rates of 5 K min−1 and 0.5 K min−1. No significant difference in spring-in is found, which is in contradiction to the findings of White and Hahn [124], who reported a decrease in PID of 12 % when decreasing the cooling rate from 5.6 K min−1 to 0.56 K min−1. Similarly, Sarrazin et al. [135] found a decrease in PID of about 12 % when decreasing the cooling rate from 4.2 K min−1 to 0.6 K min−1. However, this may arise from the fact that in the studies of White and Hahn [124] and Sarrazin et al. [135], strongly asymmetric [0n/90n]<sup>T</sup> layups are used to manufacture initially flat composite stripes. Thus, the specimens should be more sensitive to relaxation than L-shaped specimens, manufactured using a symmetric quasi-isotropic stack, as done by Groh et al. [138]. This is supported by findings of Svanberg and Holmberg [5], who found no evidence of an influence of cooling rate on spring-in of L-shaped specimens, neither when cooled right after curing nor when undergoing an additional postcure step. Another reason for the divergence of experimental findings may be related to annealing or physical aging which is more pronounced if cooling rates are lowered. White and Hahn [124], who reported a decrease in PID when using slower cooling rates also investigated the impact on the mechanical performance of the specimens. Both, transvere modulus as well as transverse strength show an influence and are found to be lower in case of the more gradual cooling of the specimens.

#### Mold constraints

Despite the use of internal or external release agents, the component will adhere to the mold surface, enabling load transfer. Hence, the mold constraints chemical shrinkage strains and, thus, affects the resulting effective amount of shrinkage as discussed in Section 6.3.1. Moreover, if the process makes use of a nonisothermal cure cycle, the mold will exhibit thermal strains which are transferred to the manufactured part. This is typically not the case in RTM processing as curing is carried out under isothermal conditions and cooling of the part takes place after demolding. The latter, however, means that constraints acting on the part change during the course of a typical RTM process, which therefore must be considered in associated simulations. The effect of load transfer between mold and part is usually referred to as *tool-part interaction* or *forced-interaction* and is especially pronounced in processes, which involve one-sided tooling and/or cooling of the part while it is still fixed by the mold. A process which inherently fulfills both of these conditions is prepreg autoclave technique. Hence, several studies exist on this matter and a comprehensive one is presented by Kappel [140], who found relative displacements at the tool-part interface as well as geometrical constraints, imposed by the rigid mold, to be the most important drivers of forced-interaction. Moreover, even small stress gradients in throughthickness direction can induce significant distortion. Albert and Fernlund [142] investigated the influence of various part-design and process parameters on spring-in of angled specimens and found tool surface and material to be important influencing factors. They propose to lower forced-interaction effects by using special release plies in addition to common release agents to decouple mold and part strain as much as possible. Furthermore, the tool material should be chosen according to its CTE, such that thermal strains of tool and part match as close as possible.

### **7.1.4 Modeling approaches used for prediction of PID**

As mentioned at the beginning of Section 7.1, spring-in is mainly driven by process-induced strains in through-thickness direction of the laminate. Hence, the used numerical model must be able to predict these strains as accurate as possible. For this reason, 2D-modeling approaches are often not suitable if not modified for this specific application. In a recent study, Kappel [143] investigated the most-relevant mechanisms leading to distortion. From numerical results it is concluded that a 3D approach should be used to address all relevant mechanisms. Moreover, methods based on the Classical Laminate Theory (CLT) should not be used when through-thickness strains dominate PID development, like it is the case for spring-in deformation. Similarly, Johnston [59] notes that a 2D plane-strain assumption leads to errors and possibly inverses the influence of transverse strain lock in quasi-isotropic (QI) layups. Using either plane-strain or plane-stress assumptions is also infeasible due to the fact that during curing, the strain in through-thickness direction is forced to zero by mold constraints, inducing stress in this laminate direction. After demolding, however, the stress in through-thickness direction may vanish as the mold constraint no longer exists and, thus, strain in through-thickness direction develops and eventually leads to spring-in. An alternative strategy proposed by Spröwitz et al. [144] makes use of a combined 2D/3D approach, in which curved sections are modeled using 3D solid elements and flat regions are represented by 2D shell elements. Another possibility is the use of a pheno-numerical approach as developed by Kappel et al. [145], in which strain data of representative L-shaped specimens is mapped on more complex components. Hence, a much simpler and computationally more efficient simulation model can be used. However, this method may not be suitable in case of a part geometry, which cannot be subdivided in representative geometrical elements or which involves non-constant radii.

The above discussion makes clear that manufacturing of composite structures is always accompanied by residual stresses. Hence, these stresses are present in specimens used for material characterization but also in produced components and consequently affect their performance. Numerical simulation can assist in evaluating the impact of these stresses by taking them into account in the design process and related structural simulations, e.g. within a holistic CAE chain as outlined by Kärger et al. [11].

## **7.2 Homogenization of UD-layer properties**

Homogenization of unidirectional ply properties is carried out using Representative Volume Element (RVE) analysis. Compared to analytical homogenization using micromechanics, this enables to investigate the viscoelastic material behavior at micro-level using the isotropic viscoelastic material model introduced in Section 6.3.3. A hexagonal fiber arrangement is used and a schematic of the resulting RVE geometry is shown in Figure 7.3. The fiber volume fraction ϕ<sup>f</sup> is defined by choosing appropriate values for height h and width w of the unit cell for a given fiber radius rf. Furthermore, periodic boundary conditions are applied at all nodes on outside faces and edges. Both, determination of the needed size of the RVE as well as formulation of periodic boundary conditions, is implemented the same way as proposed and described in more detail in Barbero [121].

**Figure 7.3:** Geometry of representative volume element (RVE) used for homogenization of unidirectional ply properties

The input data needed for homogenization comprise mechanical properties of the resin in glassy and rubbery state as well as orthotropic properties of the used fiber. Properties of the resins are listed in Table 7.1. Resin SCER3 refers to the resin Hexcel 8551-7, which is used in the parameter study in Section 7.3.6. Properties and model parameters of this resin have been taken from Simon et al. [107]. Parameters of the Prony series used in the viscoelastic material model are taken from Simon et al. [107] which were derived for Hexcel 8551-7 resin. Properties of the individual fiber types are given in Table 7.2.


**Table 7.1:** Material parameters of resins used for homogenization of UD-ply properties

**Table 7.2:** Material parameters of fiber types used for homogenization of UD-ply properties


The implementation of the transversely isotropic material models introduced in Sections 6.3.4 - 6.3.6 is based on the individual components of the stiffness tensor Cij , rather than engineering constants. Moreover, in case of the viscoelastic model, the amount of stiffness assigned to a specific relaxation time is defined by weight factors within the Prony series. In general, these factors may have different values for neat resin and a unidirectional ply. This was taken into account during homogenization in a study of Kwok and Pellegrino [122]. According to their results, the weight factors do not vary much between neat resin and unidirectional ply, and are thus assumed to be the same. Hence, the aim of the homogenization is to identify the five independent stiffness tensor components for both, rubbery and glassy material state, namely Cr,11, Cr,22, Cr,12, Cr,23, Cr,<sup>44</sup> and Cg,11, Cg,22, Cg,12, Cg,23, Cg,44. This can be achieved by calculating the stress-strain response of the RVE for three different load cases [122]. In each load case, ε1, ε<sup>2</sup> or ε<sup>4</sup> is set to a small amount of strain while all other strain variables are set to zero. Depending on the material model complexity, two strategies can be used to derive UD-ply parameters from the resulting stress-strain data. Firstly, the strain can be applied during a very short time such that the material behaves predominantly elastic or glassy. Afterwards, the strain is kept constant for a very long period of time such that stress relaxation emerges. This way, the glassy properties can be identified by analyzing the answer of the material to the abruptly applied strain at the beginning of the load case, whereas the rubbery properties can be derived from the stress that remains after relaxation took place. The alternative to this strategy is motivated by the fact that the material can be forced to behave either glassy or rubbery by choosing a temperature far below or above of the current value of T<sup>g</sup> (α). This strategy doubles the number of needed simulation runs but bypasses the necessity of predicting stress relaxation for a very long period of time, which the material needs to finally reach the rubbery plateau. Furthermore, the use of automated fitting routines is facilitated since the stiffness does not change by orders of magnitudes during a single simulation.

The cure degree dependency of the individual stiffness tensor components is assessed by repeating the identification procedure not only for the fully cured material, but also at intermediate states of cure. Figure 7.4 shows the fitted stiffness components for glassy state of a unidirectional ply of HexPly 8551-7 (SCER3/IM7) material. As is expected for a composite material in glassy state, no influence of the cure degree is visible.

**Figure 7.4:** Cure dependency of the fitted stiffness components of a unidirectional ply of HexPly 8551-7 material in glassy state

The cure degree dependency of the rubbery stiffness tensor components is shown in Figure 7.5. Here, the cure degree has a significant influence on the mechanical performance of the unidirectional ply. During curing, all stiffness components but C<sup>11</sup> increase by about four orders of magnitude. For implementation of the material models introduced in Chapter 6, the cure dependency of the stiffness tensor components is modeled using the approach of Adolf and Martin [115], in the form of Equation 6.71. The prediction of this equation when using a value of 8/3 for parameter κ is plotted as dashed lines in Figure 7.5 for each component. The good agreement with the results of the homogenization proves the suitability of the used model with the exception of component C11, which needs to be modeled cure-independent. This is expected since this quantity represents the stiffness of the unidirectional ply in fiber direction, which is dominated by the mechanical properties of the fiber and, thus, should not depend on cross-linking of the resin.

**Figure 7.5:** Cure dependency of the fitted stiffness components of a unidirectional ply of HexPly 8551-7 material in rubbery state

## **7.3 Spring-in of angled brackets**

Validation of the mathematical models used in this thesis is carried out using spring-in angles of L-shaped components. For this purpose, three series of specimens using resins FCPUR, SCER1 and FCER3 have been manufactured.

## **7.3.1 Experimental setups used for specimen manufacturing and PID quantification**

Specimens used in this thesis are manufactured using two different molding and processing concepts: Resin transfer molding (RTM) and wet compression molding (WCM). The former is applied in a similar way as presented by Groh et al. [138], who kindly provided their RTM tool. Hence, good comparability with results of their studies is provided. The WCM concept uses a mold specifically designed for this process. Both manufacturing options are discussed in more detail in the following.

### Specimen manufacturing using RTM

Figure 7.6 shows a schematic of the mold cavity as well as a picture of the upper and lower mold. Similar to the procedure described in Groh et al. [138], preforms are prepared using a flexible membrane and infrared heating, which ensures activation of the binder. Resin injection and curing is carried out under isothermal conditions as defined by the corresponding recommended cure cycle (cf. Table 2.2). The preforms were manufactured at the Institute of Composite Structures and Adaptive Systems at the German Aerospace Center (DLR) in Braunschweig, Germany. Manufacturing of the specimens was conducted at the Fraunhofer Institute for Chemical Technology (ICT) in Pfinztal, Germany. More details regarding the applied processing can be found in the work of Groh et al. [138]. Parts of the study were also published collaboratively in Bernath et al. [15].

Since inlet and outlet region of the mold are placed outside of the actual specimen geometry, the manufactured parts require trimming. This is done using water jet cutting, which ensures no influence of thermal degradation as may be present when using traditional milling. Afterwards, the produced parts have a flange length of roughly 50 mm. Furthermore, the nominal enclosed angle of the specimens is 90° and the inner and outer radius of the L-shaped part is 5 mm and 7 mm, leading to a nominal laminate thickness of 2 mm. The length of the

**Figure 7.6:** Visualization of the RTM setup used to manufacture L-shaped specimens (taken from Groh et al. [138])

specimens was chosen based on the used layup and its processability. When using highly anisotropic laminates, for instance a purely unidirectional layup in which all plies are oriented the same, fiber washing occurs due to the flow of resin. This was counteracted by introducing additional small stripes of fabric in the middle of the mold to increase fiber volume content and friction in order to prevent a deviation in fiber direction. As a consequence, only the outer portions of each specimen can be used for spring-in characterization. Hence, two specimen configurations have been used. Firstly, a *long* one, which possesses a length of about 218 mm and, secondly, a *short* one with a length of 50 mm, obtained by cutting a specimen in four parts of equal length. In the following sections, the used specimen configuration is referred to as *RTM-L* and *RTM-S* to indicate long and short specimen geometry, respectively. Figure 7.7 exemplary shows both specimen types manufactured using the described procedure.

#### Specimen manufacturing using WCM

Manufacturing of L-shaped specimens using wet-compression molding was conducted at the Fraunhofer Institute for Chemical Technology (ICT) in Pfinztal, Germany. Moreover, a mold specifically designed for this process was made. One of the aims during the design process was to demold a part that already has the dimensions of the final specimens and, thus, does not need any trimming. WCM processing is able to fulfill this requirements since no inlet or outlet port

**Figure 7.7:** RTM-L and RTM-S specimens manufactured using RTM mold and process

is needed. Furthermore, due to the short flow paths of the small specimen geometry, fiber clamping is unnecessary. Figure 7.8 shows the final mold design. On the left, the lower mold is depicted while on the lower right, the upper mold is shown. Each part of the mold is tempered using an individual water temperature control unit which enables to set different temperatures for upper and lower mold. This was used here to compensate for small deviations in temperature between these molds due to differences in mass and flow channel design.

In order to achieve a good specimen quality, application of vacuum during WCM processing is mandatory. However, contrary to RTM, resin flow already starts during closing of the mold. Hence, it is vital to apply a sufficient level of vacuum already before the resin gets in contact with the moving upper mold. This functionality is provided by a circumferential lip seal which enables to evacuate the cavity from a remaining mold gap of about 5 mm on. A detail view of the lip seal is shown on the upper right of Figure 7.8. The groove further inside acts as a resin trap and prevents excess resin from reaching the lip seal and the vacuum port.

**Figure 7.8:** Mold used for manuacturing of L-shaped specimens with WCM process

Another factor that influences specimen quality is the ability to control the closing speed. During manufacturing of the first specimens, it was observed that if the mold is closed too fast or in an uncontrolled fashion, the resin flows mostly on top of the preform rather than in through-thickness direction. As a result, large dry spots form and render the produced specimen unusable. A more controlled way of closing the mold was implemented by using a pair of screws which were initially intended as jacking screws for opening of the mold after completion of the process. During preparation of the mold, these screws are screwed in such that they protrude by 5 mm into the gap between upper and lower mold. This way, the gap remains at this level during evacuation of the cavity and afterwards the mold is closed in a controlled way by slowly unscrewing both screws. By using this procedure, good specimen quality is achieved.

Prior to specimen manufacturing, the mold is prepared by applying the mold sealer Mikon 699 MC and release agent Mikon 700 MC, both of which were kindly provided by Münch Chemie International GmbH. Mold sealer was applied only once while the release agent had to be reapplied as soon as a degradation of the release effect was observed during specimen manufacturing. Moreover, due to the design of the mold, demolding of the part is much easier if it sticks to the upper mold. Hence, the release agent of the lower mold was refreshed more often to provoke this behavior.

After preparation of the mold, a preform is manually stacked and then introduced into the lower mold. While the preform adopts to the mold temperature, the correct amount of resin is dosed and mixed. For reasons of reproducibility, a lab-scale dosing and mixing unit is used. The resin is poured onto the curved section of the preform midway along the width of the specimen. Due to the fact that this region is the lowest point of the cavity, the resin accumulates there and is only forced to flow up both flanges during closure of the mold. Hence, after introduction of the resin, the upper mold is put in place and lowered until a gap of 5 mm remains using the jacking screws as described above. At this point, the vacuum valve is opened and the cavity is evacuated while at the same time, the upper mold is successively lowered until full closure of the mold. This step is finalized by fastening of the screws positioned in a circumferential fashion around the mold to ensure complete closure during the subsequent curing stage. After curing, the mold is opened and the part is demolded. The latter is much easier if the part sticks to the upper mold since this enables to remove the specimen without the need to apply any relevant forces. The opposite is true if the part sticks to the lower mold since then, demolding forces cannot be avoided and may vary significantly from part to part. Figure 7.9 exemplary shows a specimen manufactured using the described procedure.

The specimens produced with WCM process have an overall length of 50 mm and a flange length of 50 mm. Furthermore, the nominal enclosed angle of the specimens is 120° and the inner and outer radius of the L-shaped part is 5 mm and 7 mm, leading to a nominal laminate thickness of 2 mm. Hence, the specimens are very similar to that produced using the RTM process described above and differ only in length and nominal enclosed angle.

#### Determination of the spring-in angle of specimens

Determination of the spring-in angle of specimens manufactured using RTM process makes use of a GOM 3D ATOS scanner. Corresponding measurements

**Figure 7.9:** L-shaped specimen manufactured using WCM mold and process

were conducted and provided by AUDI AG as part of the project SMiLE [54]. Each scan yields a 3D representation of the specimen and therefore enables to determine the enclosed angle and subsequently the amount of spring-in. The latter is carried out in a way similar to that presented by Groh et al. [138] using the software GOM Inspect. However, instead of evaluating the enclosed angle based on the normal vectors of the fitted planes, the intersection lines of each fitted plane with a X-Y-plane located at the center of the specimen are used as depicted in Figure 7.10. This procedure is easier to implement in a post-processing script used for evaluation of the numerical results.

Evaluation of the spring-in of specimens manufactured using WCM process was carried out using a different approach since the 3D scanner used for RTM specimens was not available. This approach is based on an optical analysis as depicted in Figure 7.11. The procedure involves taking photos of the specimens placed on a white sheet of paper and in an upright position as shown in Figure 7.11 a). The white background is vital in order to create enough contrast between specimen and background. Furthermore, the angle will be evaluated using the edge along the outer radius. Hence, the camera is positioned slightly angled relative to the part to enable easier detection of the outer edge by the algorithm. The error introduced by this angulation is negligible. In fact, the angle is about 1° which results in an error of the determined spring-in angle of roughly 0.009°.

**Figure 7.10:** Evaluation of the spring-in angle of RTM specimens using GOM Inspect software

The success of automatic edge detection depends largely on how good the edge of interest can be separated from other disturbing elements like the background. Therefore, the raw image data is modified by converting to grayscale as well as adjusting brightness and contrast until the outer edge is as clearly visible as shown in Figure 7.11 b). This step is similar to what is often referred to as binaryization of an image. The resulting image data is then analyzed using an algorithm which makes use of a *Canny edge detector* [146]. This usually yields multiple detected edges of which the longest ones are extracted. These edges are depicted in Figure 7.11 c) by a red and blue line. Based on the slopes m<sup>A</sup> and m<sup>B</sup> of these lines, the enclosed angle of the specimen is calculated. This procedure is applied to both sides of the specimen in order to account for any gradients along the extrusion direction of the profile. Finally, a mean value of front and backside angle is calculated.

**Figure 7.11:** Procedure of optical evaluation of the spring-in angle of WCM specimens

## **7.3.2 Simulation model and boundary conditions**

The numerical model used for PID simulations of L-shaped specimens makes use of the kinetic modeling introduced in Chapter 3 as well as the chemothermo-mechanical material models introduced in Chapter 6. For this purpose, these models have been implemented in form of UMAT subroutines to make them available during solving using the implicit Finite Element Method (FEM) solver *Abaqus Standard*.

As introduced in Section 7.3.1, the applied RTM and WCM process yield different specimen geometries, which differ in length and enclosed angle. Figure 7.12 shows the discretization of an RTM-S specimen which comprises 2016 quadratic solid elements of type C3D20R. The mesh consists of four layers of solid elements in thickness direction as well as six elements in circumferential direction at the curved section of the mesh. These values provide a good compromise of accuracy and computation time, as was evaluated in a preliminary mesh convergence study. Simulations using the RTM-L specimen configuration use a similar mesh which differs only in element count in extrusion direction of the profile. Hence, these specimens are meshed using 8928 C3D20R elements. The mesh used for WCM specimens is equal to that of short RTM specimens (RTM-S) and only differs in enclosed angle.

**Figure 7.12:** FEM mesh used for PID simulations of L-shaped specimens. Marked nodes are used for evaluation of the enclosed angle during post-processing of the simulation data.

The analysis consists of three steps: Curing, demolding and cooling. During the first, all nodes located at the surface of the part are fixed in all directions to simulate the interaction with the rigid mold. During the demolding step, this fixation is removed within a time span of 1 s and the part is allowed to deform freely during the subsequent cooling step. During demolding, the temperature of the part remains unchanged. Hence, it is assumed that the specimen temperature during demolding does not significantly differ from the curing temperature. During subsequent cool down, a constant cooling rate of 10 K min−1 is used and the simulation ends as soon as room temperature (20 ◦C) is reached.

The temporal evolution of the enclosed angle is evaluated using a post-processing script. In a first step, the displacement data of nodes located at the symmetry plane of the specimen (cf. Figure 7.12) are used to derive a linear fit for each flange of the geometry. The slope of each fit is then used to calculate the enclosed angle as well as the corresponding spring-in angle.

### **7.3.3 Specimens manufactured using resin FCPUR**

The FCPUR specimens were manufactured using the RTM mold introduced in Section 7.3.1. Table 7.3 gives an overview on the specimen configurations used in the FCPUR-study. Layup QI1 is the same as was also used by Groh et al. [138] in their study of the spring-in of specimens manufactured using FCER1 (Momentive Epikote 05475) and SCER2 (Hexcel Hexflow RTM6) resin. More details regarding the used fiber textiles are given in Section 2.1.


**Table 7.3:** Layups used for manufacturing of FCPUR specimens with nominal angle of 90° (FVC = Fiber Volume Content in %)

\* Stacking sequence is given from inner to outer radius

Two asymmetric layups have been chosen to analyze the impact on spring-in. For this purpose, two layers of −45°/45° orientation are removed from the QI1 layup and two unidirectional and 0° oriented plies of the same fiber type and areal weight are introduced instead. The asymmetry is generated by either locating the newly added UD plies at the inner or outer radius of the specimen, resulting in layups AS1 and AS2, respectively.

Layups UD0, UD45 and UD90 consist of six layers of equally oriented unidirectional NCF textile. These specimens are chosen in order to analyze the dependency of spring-in on the orientation of the fibers in the curved section of the part. Since spring-in is affected strongly by matrix dominated properties like low stiffness and high chemical and thermal strains, it is favorable to identify the dependency using differently oriented unidirectional stacks, in which no other oriented ply exists which would make interpretation of the results more difficult.

Figure 7.13 compares experimentally measured spring-in angles of layups QI1, UD0, UD45 and UD90 with numerically predicted values. The latter is further separated in predictions of the path-dependent (FEM/PD) and CHILE (FEM/CH) model.

The experimentally measured spring-in of the quasi-isotropic QI1 specimens is much higher than that of the unidirectional specimens. This is mainly caused by transverse strain locking that leads to an increase in through-thickness strains and thus spring-in. This mechanism is not present in purely unidirectional layups since the process-induced strains can develop freely in transverse direction. Furthermore, spring-in decreases as the fiber direction of the unidirectional layups approaches 90°. In theory, a specimen having only fibers oriented along the extrusion direction of the specimen (UD90) should exhibit zero spring-in. In this study, a mean angle of 0.06° is measured. This deviation is most likely caused by demolding forces, to which these specimens are particularly sensitive. This is supported by the amount of scatter of the unidirectional specimens, which increases as the fiber angle approaches 90°.

**Figure 7.13:** Spring-in angles of FCPUR specimens manufactured using different layups

The numerically predicted spring-in values in general show good agreement with the experimental results. Deviations in case of UD90 specimens are reasoned by the fact that the simulation does not take into account forces applied during demolding. Results of the CHILE (CH) and path-dependent (PD) material model are identical, which is expected since no devitrification occurs.

**Table 7.4:** Detailed comparison of laminate stack of symmetric (QI1) and asymmetric (AS1/AS2) layups


IR = inner radius, OR = outer radius

Figure 7.14 compares spring-in of symmetric QI1 specimens with that of the asymmetric layup AS1 and AS2. The added 0°-oriented UD plies at the inner radius lead to a decrease in spring-in. This can be better explained using a more detailed representation of the stacking sequences as shown in Table 7.4. Here, the first ply is situated at the inner radius while the sixth ply is at the outer radius. Compared to the symmetric layup QI1, the mid-plane of the AS1 stack is shifted to the left. Hence, the inner half of the layup consists to two thirds of UD-0° plies while this is only true for one third of the outer half of the layup. The remaining plies are oriented 90° and ±45° and induce high residual strains in tangential direction. In combination with the less residual strains at the inner half of the layup due to the 0°-oriented plies, this leads to a reduction in spring-in. Moreover, the stiffness of the UD-0° plies in tangential direction is much higher than that of 45° or 90°-oriented plies. Hence, the impact of tensile stresses caused by the latter is small and limited to plies three to six. This mechanism applies in the reverse sense to the AS2 layup and is more typical for warpage rather than spring-in. Thus, PID of L-shaped specimens with asymmetric layups is a superposition of spring-in and flange warpage.

The additional flange warpage in case of asymmetric layups AS1 and AS2 leads to a PID state which depends on the width and length of the L-shaped specimen. Moreover, distortion and thus spring-in angle becomes a function of the lengthwise position, as shown in Figure 7.15. This fact makes it difficult to compare simulation result and experiment only based on a single spring-in angle at the center of the specimen. Moreover, extracting spring-in angles is often accompanied by various fitting operations, which may be different depending on the software used for evaluation of experimental and numerical data. Especially in case of superimposed flange warpage, this can lead to errors. Hence, more complex distortion states should be evaluated on component level, as becomes evident when comparing the apparently bad prediction accuracy of AS2 series in Figure 7.14 with the overall good agreement of simulation and experiment in Figure 7.15.

**Figure 7.14:** Spring-in angles of FCPUR specimens manufactured using symmetric and asymmetric layups

### **7.3.4 Specimens manufactured using resin SCER1**

In this section, experimental spring-in data of specimens manufactured using the slow-curing epoxy resin SCER1 is discussed. These specimens are produced using the WCM process as described in Section 7.3.1. Hence, they exhibit a nominal specimen angle and length of 120° and 50 mm, respectively.

Table 7.5 gives an overview of specimens manufactured using SCER1 resin. Layups CP1 and CP2 are of cross-ply type and differ only in their relative orientation to the tangential direction of the specimens' curved area. Layup CP1 is also manufactured using unidirectional UD200 fabric, made of the same PX35 fiber but with a different areal weight of 200 g m−2. Hence, the number of plies also differs in order to keep the resulting fiber volume fraction in a similar range. Quasi-isotropic layup QI1 also is manufactured using both biaxial and unidirectional textiles. However, in this case, the unidirectional fabric exhibits

**Figure 7.15:** Experimentally measured and numerically predicted surface deviation of a AS2 specimen

**Table 7.5:** Layups used for manufacturing of SCER1 specimens with nominal angle of 120° (FVC = Fiber Volume Content in %)


an areal weight of 150 g m−2 and, thus, each unidirectional layer of the biaxial and unidirectional fabric is of similar weight and thickness. Hence, the only difference should be due to sewing pattern and threads. The unidirectional plies enable much more freedom with respect to stacking sequence. This is used to investigate the impact of accumulating 0° or 90° plies at the outside of the laminate, represented by series' QI2-UD and QI3-UD. It would have been consequent to use UD150 fabric also for CP1-UD specimens since that way, each unidirectional layer would be of same weight and thickness, similar to QI1-BX and QI1-UD. However, the fabric and layup in this case was prescribed by the project for which they have been investigated.

**Figure 7.16:** Spring-in angles of SCER1 specimens of different layups manufactured using unidirectional fabric and a cure temperatures of 120 ◦C

The difference in spring-in angle of unidirectional layups QI1-UD, QI2-UD and QI3-UD is shown in Figure 7.16. The mean spring-in angle seems to decrease if 0°-oriented plies are placed at the outside of the stack. A similar trend is evident for scatter of the experimental values. While it is the lowest in case of QI3, it is highest for layup QI2, most likely due to the fact that all 90°-oriented plies are placed at the outside of the stack. This introduces high process-induced strains which, due to the location relative to the neutral fiber, can have significant influence on PID. However, it also means very low stiffness and strength in tangential direction of the curved area. Hence, the increased mean value as well as the high scatter may be caused by demolding forces which are applied in order to demold the part from the mold. These forces significantly differ depending on whether the part remains in the lower concave or in the upper convex mold. In case of the former, bending of the specimen is inevitable and its magnitude depends on the level of adhesion that acts between part and mold. Demolding of the part from the upper mold, however, can be achieved without applying any bending moments thanks to its convex design. Consequently, if forces have to be applied during demolding, an increase in spring-in as well as scatter is likely.

**Figure 7.17:** Spring-in angles of specimens QI1-UD, QI2-UD and QI3-UD, separated according to demolding conditions (demolded from upper or lower mold)

In Figure 7.17, experimentally measured spring-in angles are separated according to the demolding conditions which were observed during production. Specimens that stuck to the upper mold obviously show much lower values of spring-in than those demolded from the lower mold. Moreover, the difference in spring-in angle is the largest in case of QI2-UD layup in which 90°-oriented plies are located at the outside. Interestingly, the much smaller spring-in values of specimens demolded from the upper mold is in much better accordance with the numerically predicted spring-in. It should be noted, however, that these results are based on a very limited specimen count since only 3 specimens were manufactured for each series, of which different amounts of specimens stuck to the upper mold: one, two and three for layups QI1-UD, QI2-UD and QI3-UD, respectively. Hence, for instance, all QI3-UD specimens stuck to the upper mold.

**Figure 7.18:** Spring-in angles of SCER1 specimens with layup CP1-UD manufactured using different cure temperatures

Specimens CP1-UD have been manufactured using two different cure temperatures in order to investigate the impact on spring-in. The results as well as the numerical predictions are given by Figure 7.18. Increasing the cure temperature from 70 ◦C to 120 ◦C triples the measured spring-in angle. Moreover, scatter is also found to increase which is most likely due to the fact that a temperature of 120 ◦C is close to the ultimate glass transition temperature of the SCER1 resin. Hence, the yield strength of the resin might be relatively low, leading to reduced residual stress levels. Contrary to UD specimens, only one of the CP1-UD specimens had to be removed from the lower mold. Hence, scatter is relatively low when compared to that of specimens QI1-UD, QI2-UD and QI3-UD. Furthermore, spring-in of the one specimen demolded from the lower mold is in line with that of the other specimens of that series. Numerically predicted spring-in angles overestimate the experimentally measured angles. However, the absolute error is in a similar range which indicates a systematic error introduced e.g. by material parameters. Furthermore, predicted values are very similar among the used model approaches.

**Figure 7.19:** Spring-in angles of SCER1 specimens of different layups manufactured using biaxial fabric and a cure temperatures of 120 ◦C

Figure 7.19 compares cross-ply (CP) and quasi-isotropic (QI) specimens manufactured using biaxial fabric and a cure temperature of 120 ◦C. The results imply more spring-in if a cross-ply laminate is oriented 45° (CP2-BX) to the tangential direction of the curved area instead of 0° (CP1-BX). This can be explained by the fact that the process-induced strains in tangential direction to the curved section are higher for this layup since, contrary to CP1-BX, there are no fibers oriented along that direction. However, the significant scatter of CP2-BX series undermines this statement. Hence, in Figure 7.20, the measured spring-in angles are again separated according to the individual demolding conditions. All of the CP1-BX specimens remained in the lower mold after opening of the mold. Interestingly, the resulting scatter is relatively low. Within the CP2-BX series, two specimens have been demolded from the lower and upper mold each. While the scatter of the latter again shows a reduced level, it remains relatively large.

**Figure 7.20:** Spring-in angles of SCER1 specimens of different layups manufactured using biaxial fabric, separated according to demolding conditions (demolded from upper or lower mold)

The slightly smaller spring-in angle of QI1-BX series compared to CP1-BX and CP2-BX might be due transverse strain locking. Since all neighboring layers are oriented perpendicular to each other, transverse strain locking should be more pronounced in purely cross-ply layups compared to a quasi-isotropic laminate, in which this condition is only partially the case.

Numerically predicted spring-in values are in general smaller than those measured experimentally. This also applies when only specimens are considered which were demolded from the upper mold and therefore exhibit smaller springin and scatter. Furthermore, predicted spring-in is found to be rather insensitive to changes of the layup.

## **7.3.5 Specimens manufactured using resin FCER3**

The application example of a car floor structure, which is explained in more detail in Section 7.4, is manufactured using resin FCER3. Hence, L-shaped specimens have been produced to enable a preliminary validation of the simulation method before it is applied to the significantly more complex model of the whole composite structure. Similar to SCER1 specimens, FCER3 specimens are manufactured using the WCM process as described in Section 7.3.1. The used layup and fiber textile equals that of QI1-BX specimens of resin SCER1 as given in Table 7.5. Figure 7.21 shows the experimentally determined spring-in angle of three specimens manufactured using cure temperatures of 100 ◦C and 115 ◦C. The spring-in angles predicted by the numerical models agree well with those experimentally determined.

### **7.3.6 Influence of process parameters on spring-in**

When using linear-elastic material models like CHILE or path-dependent models for numerical prediction of process-induced distortions, viscous effects are assumed to be negligible. As discussed in Sections 6.3.2 and 6.3.7, this is based on the assumption that chemical shrinkage develops at such a slow rate that the stiffness of the material at the end of the curing step equals that of the rubbery state. Hence, only during subsequent cooling, the material transitions to the glassy state. This furthermore implies that strain rates caused by demolding of the part do not contribute to a stiffness increase as would be the case when using a viscoelastic model. However, the latter would also allow these stresses to relax

**Figure 7.21:** Spring-in angle of FCER3 specimens with layup QI1-BX manufactured using cure temperatures of 100 ◦C and 115 ◦C

afterwards which is not possible in case of an elastic model. In this section, the counteracting mechanisms of generation of stress due to demolding strains and relaxation of that stress afterwards due to elevated temperature are investigated using a viscoelastic approach. For this purpose, the model and parameter set presented by Simon et al. [107] for Hexcel 8551-7 resin are used in combination with Hexcel IM7 fiber. This material combination is also known as Hexcel HexPly 8551-7.

The study focuses on the spring-in of L-shaped specimens with a nominal enclosed angle of 90° and a fiber volume fraction of 50 %. Table 7.6 lists the process parameters which are altered between the individual simulations. The cure temperature T<sup>c</sup> is varied to show the interaction of T<sup>g</sup> and T<sup>c</sup> and the resulting influence on the ability of the material to relax. In order to ensure that equal levels of cure degree and T<sup>g</sup> are reached for all simulations of the study, the minimum cure temperature is chosen only 5 ◦C below Tg,<sup>∞</sup> = 163 ◦C of the


**Table 7.6:** Parameters varied to study the influence of cure temperature <sup>T</sup>c, demolding time <sup>t</sup>dem and cooling rate hcool on the evolution of spring-in of L-section specimens

resin. The second and third cure temperature setting is chosen equal to Tg,<sup>∞</sup> and 5 ◦C above Tg,∞. Hence, relaxation is expected to be much more pronounced compared to the lower temperature level. The cure time is set to 1 <sup>×</sup> <sup>10</sup><sup>4</sup> s for all simulations, which according to Simon et al. [107] ensures full cure for all chosen cure temperatures.

The varied process parameters are demolding time tdem and cooling rate hcool. Both are varied in logarithmic space to analyze the influence over a broader range. The term *demolding time* refers to the duration of the analysis step, during which the fixed boundary conditions due to the mold is gradually removed. It should be clear that this process step is much more complex in a real demolding process. However, as long as no detailed demolding simulation is available, an assumption must be made when using viscoelastic models. Hence, this quantity is investigated in this sensitivity analysis in order to evaluate its influence on prediction of PID.

The suddenness of demolding influences the occurring strain rates that act while the specimen deforms from its constrained in-mold geometry to its final deformed state. According to the representative model introduced in Section 6.3.2, this should lead to smaller spring-in angles right after demolding when using short demolding times. The opposite is expected for more gradual demolding since it allows for significantly more relaxation to take place and therefore enables frozen-in chemical shrinkage to increasingly participate in the development of PID. Figure 7.22 shows the evolution of the spring-in angle during demolding of a specimen cured at a temperature of 163 ◦C. Due to the broad variation in

**Figure 7.22:** Influence of demolding time on the evolution of the spring-in angle during demolding of L-shaped specimens cured at 163 ◦C

demolding time, the X-axis is shown normalized. The suddenness of demolding has a significant influence on the spring-in angle at the end of the demolding step. As expected, giving the material more room for relaxation leads to an increase in PID. Moreover, the curves seem to converge for <sup>t</sup>dem <sup>&</sup>gt; <sup>1</sup> <sup>×</sup> 103 s indicating a fully relaxed state.

The effect of demolding time on spring-in at the end of demolding implies a significant influence on the predicted PID. However, this finding is quickly put into perspective when analyzing the evolution of spring-in as it continues during cooling. This is shown in Figure 7.23 for a cooling rate of 10 K min−1. The initial value of spring-in at the beginning of the cool down step corresponds to the value at the end of demolding as shown in Figure 7.22. However, what seemed to have a significant impact on PID prediction proves to be much less important. A good example is given by the curves of demolding times 0.1 s, 1.0 s and 10.0 s which show a significant difference in spring-in right after demolding. This difference vanishes quickly at the beginning of the cool down step since relaxation times

**Figure 7.23:** Influence of demolding time on the evolution of the spring-in angle during cooling with 10 K min−1 of L-shaped specimens cured at 163 ◦C

are still rather short due to a temperature close to Tg. This ability to relax stops as soon as the temperature has decreased by about 5 ◦C, indicated by a kinking of the curves and a strictly parallel progression afterwards. Hence, the time span during which relaxation can take place is defined by the applied cooling rate which is investigated next.

Figure 7.24 shows the evolution of the spring-in angle during cooling of a specimen cured at a temperature level of 163 ◦C and demolded within 0.1 s. The lowest value of tdem is chosen to demonstrate the influence of the cooling rate on the predicted PID as pointed out in the previous section. When the cooling rate is decreased, relaxation can take place for longer times before this is stopped by the falling temperature. Hence, the curves of cooling rates of 1 K min−1 and lower converge and eventually lead to similar levels of spring-in. This effectively cancels out the influence of the demolding time. Increasing the cooling rate otherwise hinders relaxation at an earlier point in time, leading to smaller spring-in angles. Although the cooling rate has an influence on the

**Figure 7.24:** Influence of cooling rate on the evolution of the spring-in angle of L-shaped specimens after curing at 163 ◦C and quasi-instantaneous demolding with a duration of tdem = 0.1 s

final state of PID, the primary cause of any difference in spring-in is still the amount of frozen-in chemical shrinkage that is allowed to contribute to PID due to relaxation. Hence, if the specimen is demolded after it was cooled down in-mold, the final state of PID is insensitive to the applied cooling rate as for example reported by Svanberg and Holmberg [5].

The above discussion shows that demolding time and cooling rate interact and possibly cancel out the impact of each other. Hence, the significance of both process properties depends on what values are applied in the real process. In the literature, cooling rates in a range from 0.5 K min−1 to 5 K min−1 are reported and applied to investigate the impact on PID (cf. Section 7.1.3). As shown in Figure 7.23, the final spring-in angle varies between about 1.0° and 1.1° when varying the demolding time and by using a cooling rate of 10 K min−1. If the same is analyzed for a cooling rate of 1 K min−1, the spring-in angle varies only between 1.07° and 1.10°. Hence, the insensitivity of PID to cooling rate reported in the literature, may be reasoned by the fact that the chosen values lay within a region in which the sensitivity is limited. Moreover, if the impact of both demolding time and cooling rate on PID is negligible, the viscoelastic prediction approaches that of an elastic model. To show this, the prediction of the CHILE model is also plotted as black dashed lines in Figures 7.22 to 7.24. During demolding, the prediction of the elastic model equals that of the viscoelastic model if a long demolding time is used. A similar trend is observed during cooling as shown in Figure 7.24. However, in case of very low cooling rates, the viscoelastic model still predicts slightly higher spring-in angles which is due to the fact that thermal strains emerge only gradually during cooling. Hence, at temperatures close to Tg, a portion of that newly created stress relaxes, leading to a small increase in spring-in.

The curing process is usually continued until the resin has reached a sufficiently high cure degree. As a result, relaxation primarily depends on temperature from that moment on. This has a significant impact on the predicted spring-in angle. Figure 7.25 shows the evolution of the spring-in angle during cooling of a specimen cured at a lower temperature level of 158 ◦C and demolded within 0.1 s. As is shown in Section 3.3, it depends on the individual resin whether it is able to reach a T<sup>g</sup> that is significantly higher than the cure temperature. In this case, the model of Simon et al. [107] predicts full cure and, thus, T<sup>g</sup> = Tg,<sup>∞</sup> at the end of the cure step. Consequently, the only difference between Figure 7.25 and 7.24 is the temperature at which demolding is carried out and cooling is started. This slows down relaxation during these steps which becomes evident when comparing both figures. While the cooling rate dependency remains significant, a much slower rate is needed in order to reach the same level of spring-in and to compensate for the decrease in temperature. If the cure temperature is increased, the contrary is the case as shown in Figure 7.26. Here, the part is cured using a temperature of 168 ◦C, which is 5 ◦C above T<sup>g</sup> = Tg,∞. As a result, increased spring-in is predicted even for rather fast cooling rates.

When comparing Figures 7.24 and 7.25, the biggest impact due to a change in cure temperature is not found in the PID prediction of the viscoelastic material model but in that of the elastic CHILE model. The predicted spring-in angle is reduced by 33 % only due to a decrease in temperature of 5 ◦C. This is caused

**Figure 7.25:** Influence of cooling rate on the evolution of the spring-in angle of L-shaped specimens after curing at 158 ◦C and quasi-instantaneous demolding with a duration of tdem = 0.1 s

by the fact that the model is calibrated to switch from rubbery to glassy state as soon as the temperature undercuts Tg. Thus, at the end of the curing step and for all of the subsequent steps, the material model already changed to glassy state in case of the lower cure temperature. Hence, only a small portion of the chemical shrinkage contributes to spring-in. This huge impact on predicted PID emphasizes the importance of accurate material characterization since this effect can be a result of two scenarios: Firstly, the material actually shows a T<sup>g</sup> that exceeds the applied cure temperature or, secondly, the fact that T<sup>g</sup> > T<sup>c</sup> is only due to an error or uncertainty introduced during material characterization or model parameterization. In case of the former, the lower amount of PID should also show up in experiments. In case of the latter, the model needs to be recalibrated each time the cure temperature is altered since the error introduced by inaccurate models usually is non-linear.

**Figure 7.26:** Influence of cooling rate on the evolution of the spring-in angle of L-shaped specimens after curing at 168 ◦C and quasi-instantaneous demolding with a duration of tdem = 0.1 s

## **7.4 Distortion of car floor structure**

The practical relevance of PID is demonstrated using a complex shaped car floor structure. The geometry is similar to that used for mold-filling simulations in Section 5.3 and differs in two aspects. Firstly, the foam cores of the two crossbeams in the middle part of the component are now included in the model. Secondly, the geometry represents the final trimmed shape of the component. Hence, regions used for fiber clamping or sealing are no longer part of the model. The PID of the trimmed part is analyzed using a GOM 3D ATOS scanner which enables to compare measured and predicted PID. Corresponding measurements were conducted and provided by AUDI AG as part of the project SMiLE [54]. As a result, the prediction accuracy of the numerical model can be evaluated based on a practical component geometry.

## **7.4.1 Simulation model and boundary conditions**

Figure 7.27 shows the geometry of the component. Areas which are cut away during trimming are colored grey. The subdivision of the model in regions of different layup equals that used in mold-filling analysis (cf. Figure 5.2). Details of the composite layup are given in Appendix A in Figure A.1.

**Figure 7.27:** Model sectioning of car floor structure used for PID simulation. Areas cut away during trimming are colored grey.

The model consists of 3 394 118 solid elements of type C3D20R, C3D15 and C3D10. A significant portion of these are needed for meshing of the foam cores. Figure 7.28 provides a close-up of the model showing one of the two foam cores as well as two regions of neat resin. The latter are caused by overlapping zones of individual sub-preforms. During mold-filling, these act like additional runners and therefore influence the flow front propagation. In a PID simulation, neat resin regions potentially influence the resulting deformation due to the high residual strains.

**Figure 7.28:** Close-up of the simulation model showing one of the foam cores and regions of neat resin caused by overlapping of sub-preforms

Due to the high element count of the model, the computational effort in terms of processing power and memory demand is substantial. Hence, the material behavior is modeled using the CHILE approach described in Section 6.3.6, which requires the least amount of state variables as well as mathematical operations. Moreover, the accuracy of this model proved to be very good in terms of spring-in prediction as shown in Section 7.3.5. The boundary conditions applied to the car floor structure are similar to that used in PID simulation of L-shaped specimens in Section 7.3.2. The only difference is a symmetry plane which helps to reduce the time and memory needed to solve the numerical model.

### **7.4.2 Experimental validation of predicted PID**

Figure 7.29 shows a comparison of experimentally measured and numerically predicted distortion of the car floor structure. The color refers to the deviation of the produced part or the predicted shape of the component relative to the nominal CAD geometry. The comparability of these results strongly depends on the alignment of the deformed geometry relative to the nominal CAD geometry. This is ensured by importing the predicted deformed geometry into GOM Inspect software and by applying the same referencing as applied to the experimentally measured component.

In the comparison shown in Figure 7.29, the maximum and minimum values of the legend have been adopted in case of the numerical result in order to facilitate comparison. The overall qualitative displacement field of experiment and simulation shares many similarities. A distinctive feature of the deformed shape is the large deviation from nominal geometry in the front region, shown on the left hand side in Figure 7.29. Contrary to the back of the component, the recessed area in the front is designed open, similar to a box structure of which the top face as well as one of the side faces have been removed. Due to the latter, the bottom face is prone to deformations in the region, in which the geometric stiffening of the box structure is missing.

When evaluating and comparing absolute values, the distortion is underpredicted by about 50 %. The error at the front is even higher as the simulation predicts a deviation of 2.07 mm whereas the experiment shows a much higher value of 12.8 mm. Similarly, the predicted PID at the central measurement point is under-predicted by a factor of 14.

The prediction accuracy of the process simulation as demonstrated above is not sufficient for use during virtual part design or distortion optimization. In contrast, the preliminary study using angled brackets showed good accuracy, especially in case of resin FCER3. Consequently, the under-prediction of PID in case of the car floor structure is surprising and in contradiction to the results of the spring-in study. The most significant difference between manufacturing of angled brackets and the car floor structure is the applied cure cycle. While in case of the former, the resin is cured for 15 min, the car floor structure is demolded after only 5 min. Hence, the cure time in case of the car floor structure is much shorter. In Figure 7.30, the result of isothermal shrinkage measurements using a parallel plate rheometer setup, as discussed in Section 6.2, is shown. Four similar shrinkage experiments were conducted and are colored light gray in Figure 7.30. According to these results, the shrinkage reaches a plateau after about 600 s, which means no further cross-linking takes place. At this point in

**Figure 7.29:** Comparison of experimentally measured (top) and numerically predicted (bottom) PID of the car floor structure. The underlying material model has been calibrated using data provided within the SMiLE project. Moreover, the value ranges have been adjusted to show the qualitative agreement more clearly.

time, the original resin volume, enclosed between the cylindrical metal plates, decreased by 1.88 ± 0.12 %. The corresponding cure degree is 92.57 %, when calculated using the parameter set given in Table 3.11. Hence, the volumetric shrinkage after full cure of the resin can be calculated by

$$
\beta\_{\rm ch} = -1.88 \,\%\, \frac{100 \,\%-\alpha\_{\rm gel}}{92.57 \,\%-\alpha\_{\rm gel}} = -2.94 \,\%. \tag{7.2}
$$

Here, αgel denotes the cure degree at gelation, which has been discussed in Section 3.4. With the assumption of a linear relationship between cure degree and chemical shrinkage, the measurements shown in Figure 7.30 can be used to parametrize a kinetic model. Furthermore, these measurements also contain information on the time of gelation and the amount of shrinkage that develops during curing. Hence, such shrinkage measurements enable to characterize the resin in a holistic way by at the same time providing data on cross-linking evolution, gelation and chemical shrinkage. However, the measurement data contains no information on the cure degree at gelation, which must therefore be determined elsewhere.

**Figure 7.30:** Isothermal shrinkage measurements of resin FCER3 at a temperature of 100 ◦C using parallel plate rheometry

The cure degree at gelation of resin FCER3 has been calculated from shrinkage measurements using the kinetic model and parameter set given by Table 3.11. The resulting value of 79.99 % seems too high as the mean value of all other epoxy resins yields 61.83 %. Hence, the prediction accuracy of the original parameter set is evaluated by comparing the predicted shrinkage strain with that of the measurements. The shrinkage strain is calculated using the following incremental approach:

$$
\varepsilon\_{\text{ch}}\left(t\_m\right) = \varepsilon\_{\text{ch}}\left(t\_{m-1}\right) + \beta\_{\text{ch}} \cdot \left(\alpha\left(t\_m\right) - \alpha\left(t\_{m-1}\right)\right). \tag{7.3}
$$

As the resin is in rubbery state during curing (ν<sup>r</sup> = 0.5), the use of the volumetric shrinkage coefficient is valid. Thus, this equation predicts the temporal evolution of the volumetric shrinkage strain. In Figure 7.30, the corresponding curve is labeled *Original fit*. Since the cure degree used to calculate the amount of shrinkage of the fully cured resin has been determined at a time of 600 s (cf. Equation 7.2), the predicted curve coincides with the measurement curves at this point in time. Before and after this time, the prediction deviates from the measurements. This indicates that the accuracy of the parameter set given in Table 3.11 is insufficient for the applied cure cycle. Hence, the shrinkage measurements have been used to gain a more suitable model fit for reaction kinetics and Tg. This is motivated by the fact that chemical shrinkage develops proportionally to the cure degree (cf. Section 6.2). Moreover, two important assumptions are made. Firstly, during isothermal curing, the T<sup>g</sup> may not exceed the reaction temperature Tc. The validity of this assumption has been discussed in Section 3.3. Secondly, the cure degree at gelation equals the mean value of all other epoxy based resins. Furthermore, the Grindling kinetic model is used since premature vitrification occurs. The fitting process is carried out using a Python script as well as Differential Evolution [79] as supplied by the Python package SciPy [80]. Table 7.7 lists the model parameters fitted to the shrinkage experiments.


**Table 7.7:** Parameters used for modeling the evolution of cure degree, <sup>T</sup>g and shrinkage strain of resin FCER3 based on the Grindling kinetic model and DiBenedetto <sup>T</sup>g model

The predicted shrinkage strain evolution using the new parameter set is labeled Shrinkage fit in Figure 7.30. The measured shrinkage is now accurately predicted by the numerical model. While the changes seem only minor in this respect, the impact on cure degree and T<sup>g</sup> evolution is significant as shown in Figure 7.31. At the time of demolding (t = 300 s), a lower cure degree is predicted by the model fitted to shrinkage data. Due to the non-linearity of the relationship between cure degree and T<sup>g</sup> towards full cure, this has even more influence on Tg. While the original parameter set predicts a T<sup>g</sup> of about 100 ◦C, the alternative shrinkage fit yields a T<sup>g</sup> of only 60 ◦C. This value will further increase during cooling depending on the temperature rate. However, cooling will take place while the resin is in rubbery state and therefore exhibits a much higher CTE. Furthermore, the mechanical properties of the composite in directions perpendicular to the fibers are also much lower than in glassy state.

Figure 7.32 compares experimentally measured PID of the car floor structure and numerically predicted PID using the model fitted to shrinkage measurements. While the prediction quality improved, PID is still underestimated by about

**Figure 7.31:** Evolution of cure degree and <sup>T</sup>g of resin FCER3 during isothermal cure at a temperature of 100 ◦C

20 % to 30 %. Moreover, the accordance between experiment and simulation in the front area of the component is even worse. However, this can be explained by a manifoldness of the numerical solution. The deformation of this initially flat portion of the component stems from spring-in of curved sections from which it is surrounded from all sides except on the front side. The spring-in subsequently leads to a buckling of the flat area and, thus, two stable deformation states are possible. Which of these states eventually persist, depends on the temporal evolution and interaction of PID of the surrounding structures, but also on external forces e.g. applied for demolding of the component. Hence, the difference between experiment and simulation in this region may be attributed to demolding forces as these are not taken into account by the simulation but can lead to significant deformations.

Figure 7.33 shows the positions of the individual ejectors used for demolding of the car floor structure. The part is manufactured upside down as depicted in Figure 7.34, which shows the injected part in the lower mold at the time of demolding. Due to the concave geometry, the part shrinks onto the lower mold and remains there also during mold opening. Consequently, the ejectors were

**Figure 7.32:** Comparison of experimentally measured (top) and numerically predicted (bottom) PID of the car floor structure. The underlying material model is calibrated using shrinkage measurements. Moreover, the value ranges have been adjusted to show the qualitative agreement more clearly.

integrated into the lower mold and push the component upwards for demolding. During this process, substantial forces and deformations act on the component as is evident by the bumps visible in Figure 7.34 at the positions of the ejectors. By comparing the positions of the ejectors with the experimentally measured PID at this locations, demolding forces obviously lead to persistent deformations which superimpose PID. Forces applied by Ejector A during demolding probably cause the large deformation in this region and define the direction of the buckling of the flat area in negative Z-axis. Moreover, Ejectors C and D also seem to contribute to the deformation state as the measurement shows additional out-ofplane deformation at the corresponding positions.

**Figure 7.33:** Positions of ejectors used for demolding of the car floor structure

The chronological sequence of demolding of individual regions of the manufactured part varies significantly from one to another. Hence, the demolding process of the specific component shown in Figures 7.29 and 7.32 is reconstructed by analyzing corresponding video footage. However, as no actual measurement of the deformation during this process step is available, only qualitative results can be drawn from this data.

**Figure 7.34:** Car floor structure during demolding. Bumps visible in the flat areas in front and rear of the component are caused by ejection forces.

Allthough all ejectors are activated simultaneously, the front part of the component comes off first while the rear section maintains its position for slightly longer. During this short time span, the flat region between the two foam cores is bent around the Y-axis as schematically depicted in Figure 7.35. This bending possibly leads to the high amount of positive PID in this area. Moreover, it also contributes to the overall part deformation by increasing the deviation from nominal geometry especially at positions furthest away from the bending axis. This may explain why the simulation underestimates PID at this locations. Furthermore, a small section of the fiber clamping at one side in the middle of the part still sticks to the mold and, thus, causes an asymmetric loading of the component.

The rear section of the component comes off next and most of the part is now detached from the mold surface. As it is further elevated by the ejectors, only the short section of the fiber clamping area still hinders complete demolding which means that the force distribution remains asymmetric until the part is completely detached from the mold. The asymmetric loading of the structure may be the reason for the slightly asymmetric distortion state of the part.

**Figure 7.35:** Schematic of the asymmetric bending of the component during demolding

Besides demolding forces, other sources of uncertainty that have not been included in the presented numerical model may also introduce error. As discussed above, the cooling process has a huge influence on the resulting PID if the component is demolded and cooled while the T<sup>g</sup> is still much lower than the curing temperature. Hence, in case of a gradual cooling, cross-linking will advance, leading to more shrinkage but also a further increase in Tg. In case of rapid cooling, only minor advancement in cure degree will occur and T<sup>g</sup> development stagnates. As a result, cooling from reaction temperature to current T<sup>g</sup> will give large residual strains due to the higher CTE of the resin in rubbery state. Consequently, the cooling rate defines how much further cross-linking occurs and what temperature range will be passed while the material is in rubbery state. In the numerical model presented above, this complex process is simplified by assuming a homogeneous temperature distribution throughout the component as well as a constant cooling rate. Both assumptions will only hold true in reality if the manufactured part is cooled within a temperature chamber in a controlled way using a very slow cooling rate. Otherwise, the rate will differ throughout the part depending on air convection and gradients in ambient air condition. Furthermore, the amount of heat exchanged with the environment will decrease as the component cools due to the continuous reduction in temperature difference. Additionally, foam cores may also introduce uncertainty as their ability to conduct and store heat differs significantly from that of the composite.

## **7.5 Discussion and concluding remarks**

In this chapter, process-induced distortion (PID) of composite structures is investigated. For this purpose, spring-in of angled brackets is experimentally analyzed for different types of resins, layups and cure temperatures. This data is then used for validation of numerical predictions using the models introduced in the previous chapters.

Specimens manufactured with FCPUR resin have been produced using symmetric quasi-isotropic and unidirectional layups of different orientations, as well as two asymmetric layups. Of the symmetric layups, the quasi-isotropic layup showed a larger spring-in angle than purely unidirectional specimens. This is caused by the interaction of neighboring plies of differing fiber orientations, which eventually leads to transverse strain locking of the individual plies. As a result, the through-thickness strain of quasi-isotropic specimens is greater than that of unidirectional ones, leading to increased spring-in. Regarding unidirectional layups, the one with fibers oriented tangential to the curved section (UD0) exhibits more spring-in than those with 45° (UD45) or 90° (UD90) fiber orientation. This is caused by the stiff fibers which are forced to a radius different from the initial one due to through-thickness strains. This induces tensile and compressive stresses at the inner and outer radius, respectively. This mechanism is the origin of spring-in and is not present in UD90 specimens because no fiber is oriented in tangential direction. Hence, in theory, spring-in of these specimens should be non-existent. In reality, however, the lack of fibers in tangential direction makes these specimens particularly prone to influences due to demolding forces, which explains the small amount of spring-in measured as well as the increase in scatter. UD45 specimens lie between these two extremes and therefore inherit a relevant amount of spring-in from UD0 specimens but also the increased scatter from UD90 specimens.

Making the layup asymmetric by placing 0°-oriented plies at either the inner or outer radius is followed by a decrease in spring-in compared to the quasiisotropic layup. It should be noted, however, that this comparison is somewhat invalid as two ±45° plies of the QI layup have been replaced by 0°-oriented ones in order to maintain a similar fiber volume content. Nevertheless, the use of an asymmetric layup leads to a superposition of spring-in and warpage and therefore challenges numerical prediction even more than symmetric ones. However, the superimposed warpage also questions the usual way of characterizing spring-in by breaking the resulting displacement field down to a single quantity, the springin angle. As warpage strongly depends on geometry and part dimensions, the same now applies to the spring-in angle. Hence, experimental characterization of PID in case of asymmetric layups should be done by taking the complete displacement field into account.

The spring-in study using SCER1 resin has been conducted with both unidirectional and biaxial fabrics. The use of unidirectional plies enables layups which are not possible when using biaxial fabrics as for instance the 0°-oriented ply of a 0°/90° layer can be placed within the layup irrespective of the 90°-oriented one. This is utilized to modify an initially quasi-isotropic laminate such that either the 0° or 90°-oriented plies are all located at the outside or in the mid-plane of the layup. The resulting spring-in is found to be lower with 0°-oriented plies at the outside and 90°-oriented plies in the center of the laminate. However, the scatter of this series is also much lower and it is found that this is mainly determined by whether the specimen had to be demolded from the upper or lower mold. In case of the former, demolding of the produced part is possible without bending of the specimen. This is not the case when the part sticks to the lower mold, which involves the application of significant forces and, thus, bending of the specimen. When decomposing the specimens according to from which mold they were demolded, it is found that specimens detached from the upper mold show less spring-in and scatter but also very similar spring-in angles, irrespective of the position of the 0° or 90°-oriented plies. Moreover, specimen demolded from the lower mold show a spring-in angle which is more than twice as high as that of specimen demolded from the lower mold if 90°-oriented plies are placed at the outside of the layup. This is similar to the behavior of purely unidirectional and 90°-oriented layups, which also shows increased scatter in the study of FCPUR resin. In case of biaxial layups, this characteristic is much less pronounced as no accumulation of e.g. 90°-oriented unidirectional plies is possible. Instead, each 90°-oriented ply is accompanied by a 0°-oriented one. Hence, the sensitivity to demolding forces is less pronounced.

Numerical prediction of spring-in is accurate for most specimens investigated in this thesis. Exceptions are the asymmetric layups of the FCPUR study as well as several series' of the SCER1 study. In case of the former, superimposed warpage hinders a meaningful comparison based on a single spring-in angle. However, a visual comparison based on the measured and predicted displacement field shows good correlation. In case of the SCER1 study, demolding conditions dominate both derived spring-in angle and associated scatter. Reasonable agreement is found between simulation and experiment in case of unidirectional layups as long as only specimens demolded from the upper mold are considered. For biaxial layups, the difference in spring-in of specimens demolded from either mold is less pronounced. However, scatter of the 45°/−45° cross-ply specimens is significantly larger than that of 0°/90° cross-ply specimens, even when demolded from the upper mold. Moreover, spring-in angles are underestimated by the simulation for this specimens. These results imply two important conclusions. Firstly, as the numerical model used in this thesis does not take demolding into account, prediction accuracy is impaired if noteworthy forces are applied during this process step. This is important since it affects PID in terms of absolute value but also scatter, and cannot be completely avoided in reality. Secondly, the simulation overpredicts spring-in of unidirectional specimens demolded from the lower mold but underpredicts that of biaxial layups, irrespective of the demolding conditions. A possible reason for this may be related to fabric specific properties like sewing thread and pattern as has also been suspected in Bernath et al. [15]. The type of mathematical model used to describe the material behavior of the matrix during the course of the process is found to have no significant influence on spring-in prediction accuracy. This can be explained by the slow rates at which chemical and thermal residual strains emerge during the process. However, the influence of demolding conditions is difficult to quantify. In order to investigate this interaction of spring-in and process specific quantities, results of a numerical study based on the evolution of the spring-in angle during the individual process steps are evaluated. These results show that abrupt demolding induces high strain rates and thus viscous stress, which potentially relaxes at the beginning of the cooling stage. However, this depends heavily on temperature and time, and therefore on the applied cooling rate. Moreover, the prediction of the viscoelastic model coincides with that of an elastic CHILE or path-dependent model for typically applied processing conditions. It should be noted, however, that this finding does not generally apply to every other composite manufacturing process. If for instance a post cure step is used, a CHILE model will most likely give inaccurate results. Moreover, elastic approaches are also insufficient for numerical simulation of the influence of applied demolding forces as the related deformation mechanisms are of viscous or plastic type.

Investigating the process-induced distortion of angled brackets is rather academic. Hence, PID of a large automotive structure is given as an example of a more practice-oriented geometry. Prediction accuracy is found to be poor when using the parameter set provided for use within the associated research project. By recalibrating the model to shrinkage measurement data acquired using a parallel-plate rheometer setup, the accuracy is improved. However, it is suspected that demolding forces also affected the resulting PID of the component, similar to what has been found in the study of angled brackets. Based on an analysis of corresponding video footage, the produced part is bent during demolding. This again underlines the importance of investigating the process of demolding in more detail, e.g. by using 3D displacement measurement during this process step. This would enable to better correlate experiment and simulation and to identify deficiencies of the numerical model, which need to be resolved to be able to take such influence into account. Important aspects in this regard are strain rate and temperature-dependent plasticity as well as the interaction between mold and part in terms of mechanical contact and heat transfer.

# **8 Conclusions**

In this thesis, the influence of cross-linking on both mold-filling and resulting part quality is experimentally and numerically investigated. During mold-filling, the resins viscosity largely depends on the cure degree and therefore has a significant impact on injection and cavity pressure. Being able to accurately predict this important quantity is vital for reliable process design and optimization. The same is true for the prediction of process-induced distortion. During curing, the resin shrinks and induces shrinkage strain which, together with thermal strains during cooling of the produced part, lead to geometrical deviation from the nominal shape. This potentially leads to increased scrap volume or requires expensive compensation by using for instance manual shimming strategies. Consequently, accurate numerical prediction is a key element in cost-effective design of both the composite structure and the associated manufacturing process.

Three resins are investigated in this study, of which two are epoxy based and one is polyurethane based. The polyurethane resin and one of the epoxy resins are fast-curing while the other epoxy resin is slow-curing. This separation is important as fast-curing resins enable short process cycle times but also challenge conventional material characterization technique. In case of the polyurethane resin, the high reactivity of the resin made isothermal DSC measurements impossible as a significant amount of cross-linking already happened before the first data point is acquired by the analyzer. Even in case of the slow curing epoxy resin, isothermal DSC runs missed a relevant portion of the curing reaction at high temperatures. However, isothermal measurement data is vital for proper prediction of the evolution of the cure degree during the process since the process is isothermal too and, thus, vitrification occurs at some point in time. Hence, the method used for fitting of the kinetic model has been modified to be able to account for this uncaptured premature curing by making the cure state at the beginning of the measurement an additional fitting parameter. This adds flexibility to the fitting process but may easily lead to non-physical results. Thus, it is important to also use data from non-isothermal and cyclic DSC measurements in order to constraint the fitting process to feasible solutions.

Viscosity characterization using rheometer measurements is another experimental technique that is prone to errors in case of fast-curing resins. Due to the high reactivity, the mold-filling process needs to be completed in a rather short time. This period of time, however, is similar to that usually needed for specimen preparation and measurement initialization. As a result, only very limited information is gathered during the process-relevant time span. Thus, when using such data for viscosity model fitting and subsequently for mold-filling simulation, a large portion of the prediction will rely on extrapolated viscosity values. As this impairs prediction accuracy of the simulation, a novel specimen preparation and loading technique has been developed in this thesis. For this purpose, a lab-scale dosing and mixing unit is mounted to a rheometer which allows to precisely inject the needed amount of resin directly into the gap of the parallel plate setup. This way, the time span of unrecorded curing prior to the start of the actual measurement has been reduced from about 30 s to only 7 s. The impact of this difference in acquired viscosity data on the result of mold-filling simulations is demonstrated based on a large composite structure. With a filling time of around 30 s, this practical example is well suited for this comparison, for which a viscosity model is fitted using two different sets of input data. The first one simulates the conventional manual measurement preparation process and, thus, is only allowed to use viscosity data from 30 s upwards. This leads to a significant underestimation of the viscosity in the process-relevant time span. The second one makes use of all acquired data points and therefore is able to predict the viscosity evolution with much better accuracy. During the actual filling process, the underestimation of viscosity leads to a substantially lower injection pressure when compared to that predicted using the accurately parametrized model. Moreover, due to the fact that freshly mixed resin is continuously injected into the mold, the consequence of a prediction error during the first seconds of the curing process persistently influences the pressure distribution during the entire duration of mold-filling. Hence, when designing pressure-controlled injection process variants like PC-RTM, accurate viscosity characterization is vital, particularly in case of fast-curing resins.

Numerical simulation of process-induced distortion (PID) is carried out using three different material models: viscoelastic, path-dependent and cure hardening instantaneously linear-elastic (CHILE). Each model has been implemented in a UMAT for use with the implicit Abaqus FEM solver. The coefficients of thermal expansion are measured using Thermomechanical Analysis (TMA). Chemical shrinkage is determined based on Video-Imaging (VI) and parallelplate rheometry measurements. The latter is also used for evaluation of the cure degree at gelation of the resins which is needed to extract the relevant portion of shrinkage from VI measurements. Predicted spring-in angles agree well with most of the experimentally determined angles derived from angled brackets. Deviations are found when using asymmetric layups or when demolding involved the application of non-negligible forces. The former is primarily an evaluation problem as the asymmetric layup induces both spring-in and warpage, which superimpose each other. Hence, the displacement field becomes more complex and dependent on specimen size and geometry. This leads to a spring-in angle which is no longer independent of the position at which it is evaluated. Specimens with asymmetric layups should therefore be evaluated using a more holistic approach and not by breaking PID down to a single spring-in angle. The influence of demolding forces is found to be significant. Although the study has not been specifically designed to investigate this influence with sufficient statistical backup, the difference in spring-in as well as the increase in scatter in cases involving non-negligible forces is striking. Moreover, this effect is even more pronounced if the outer layers of the laminate are oriented along the extrusion direction of the profile. In this case, demolding forces act transverse to the fibers. In this direction, the composites material behavior is dominated by resin properties. Due to the fact that during demolding, the temperature is very close or above the T<sup>g</sup> of the resin, the material is in rubbery state, in which both stiffness and strength are considerably lower than in glassy state. Hence, the composite is more prone to persistent deformations due to demolding forces, which is evident in the presented absolute spring-in values and associated scatter.

The primary reason for consideration of the viscoelastic material model in this thesis is the possibility to evaluate the influence of suddenness of demolding as well as cooling rate on the development of spring-in. In a parameter study, it is found that a more sudden demolding leads to a decrease in spring-in as the involved strain rates increase. This effect is compensated by relaxation in the subsequent cooling step, which therefore depends largely on the applied cooling rate as relaxation times significantly increase once the resin has been cooled to a temperature substantially below Tg. Consequently, similar spring-in is predicted in case of sudden demolding combined with slow cooling, or gradual demolding and rapid cooling. Furthermore, it is shown that in these cases, the evolution of spring-in is similar to that predicted using a CHILE model. However, whether a CHILE model is able to correctly reproduce the behavior of the inherently viscoelastic material, depends strongly on the chosen parameters. Contrary to a viscoelastic model, the influence of strain rate and temperature on material stiffness is not contained in the model. Instead, the stiffness evolution during the process is estimated based on the expected process conditions. It is therefore vital to understand the interaction of material behavior and process variables to be able to derive valid assumptions. This is especially crucial in the vicinity of glass transition, as it involves a substantial change in stiffness.

As a more practical application example of PID simulation, the car floor structure used for mold-filling simulations is also used for prediction of PID. Compared to angled brackets, this geometry represents a much more complex situation as PID mechanisms emerge at many different locations of the structure and thus interact with each other. Comparison with experimentally measured PID revealed a rather poor accuracy of the simulation. While the prediction of the underlying distortion mechanisms agrees well in qualitative terms, absolute PID is significantly underestimated. Given the very good agreement between experiment and simulation based on spring-in of angled bracket specimen, this is unexpected. The main difference between the manufacturing process applied to angled brackets and that applied to the car floor structure is in the cure time, which is much shorter in case of the latter. Moreover, shrinkage measurements show that after this shorter cure time, the curing process is not completed. Hence, it is suspected that the inaccuracy of the simulation is partly caused by an inaccurate prediction of the evolution of the cure degree and T<sup>g</sup> when using the parameter set provided within the associated research project. Therefore, model fitting has been repeated using the aforementioned shrinkage measurements, which have been conducted using a parallel-plate rheometer setup and by applying the very same temperature conditions as also used in manufacturing of the composite structure. These measurements are well suited for model fitting as each individual measurement contains information on time of gelation as well as temporal evolution of shrinkage, and thus cure degree. Using the newly fit model, the prediction accuracy improves significantly but a substantial quantitative error still exists. By analyzing the demolding process, it is found that distinctive features of the component's distortion can be explained by forces which act on the part during demolding. This conclusion is supported by the study of angled brackets, which also revealed such influence on the resulting PID.

# **9 Outlook**

By using an automated specimen preparation and introduction process, the accuracy of viscosity characterization of fast-curing resins improved significantly. However, this may be even more improved by avoiding the need for a loading gap, which differs from that used during the measurement. This can be achieved by heating of the individual resin components, which will reduce viscosity and thus allows to use needles with smaller diameter. This way, moving the sample holders to the desired measurement gap prior to the start of the measurement is no longer needed.

From an economic point of view, the prediction of process-induced distortions of composite structures is a good starting point. However, the actual important and relevant goal is minimizing PID. This can be achieved by numerical optimization for instance by altering the mold geometry (tool compensation strategy) as proposed by Jung et al. [147] or Wucher et al. [148]. An important requirement of numerical optimization is the availability of a computationally efficient simulation model. The methods used and developed in this thesis have been chosen such that the simulation result is not impaired by shortcomings of the underlying numerical approaches. Thus, 3D continuum mechanics and solid elements have been used as these allow accurate consideration of throughthickness strains. These methods, however, are inferior to 2D or 2.5D approaches in terms of computational efficiency, which is particularly important in case of large structures.

In both experimental studies presented in this thesis, demolding forces are suspected to substantially influence the resulting part distortion. This is problematic

as the demolding process is currently neither experimentally nor numerically researched. Moreover, its impact on part quality is often neglected and not taken into account during design of the demolding concept. This is mostly due to the lack of available information and knowledge. Hence, investigating this process step in more detail may enable to reduce PID already simply by choosing a suitable demolding concept which avoids the application of unnecessarily high forces. However, complete avoidance of these forces is impossible. Hence, it may be advantageous to minimize demolding forces and at the same time use them to purposefully counteract PID. Numerical simulation can again be used to identify an optimal strategy and experimental data could be gathered by using optical 3D coordinate measuring already during demolding and not only for characterization of PID of the finished part. This would enable to better understand and differentiate the origins of PID.

# **A Composite layup of car floor structure**

The composite layup of the car floor structure used in mold filling and PID simulations is detailed in Figure A.1. Red line-like areas are regions of neat resin which are caused by overlapping tolerances of subpreform assembly. For each stack, ply orientation and corresponding textile type are given. Details of the individual fiber textiles are given in Table 2.1. Table 5.1 lists permeability values of the textiles BX300 and UD300. As textile UD200 is very similar to UD300 in terms of fiber type and architecture, it is assumed that the permeability is similar too.

**Figure A.1:** Composite layup of car floor structure

# **List of Figures**








# **List of Tables**




Castro-Macosko viscosity model . . . . . . . . . . . . . . . . . . 77

4.2 Identified model parameters for FCER3 resin and

# **Bibliography**


simulation of plain woven fabrics for wet compression molding," in *Proceedings SAMPE Europe Conference & Exhibition 2017 (Vol. 1 of 2)*, Stuttgart, Germany, Nov. 2017, pp. 86–93.


composites," *Materials Science and Engineering: A*, vol. 452-453, pp. 483–498, Apr. 2007.


[148] B. Wucher, P. Martiny, F. Lani, T. Pardoen, C. Bailly, and D. Dumas, "Simulation-driven mold compensation strategy for composites: Experimental validation on a doubly-curved part," *Composites Part A: Applied Science and Manufacturing*, vol. 102, pp. 96–107, Nov. 2017.

### **Karlsruher Schriftenreihe Fahrzeugsystemtechnik (ISSN 1869-6058)**

Herausgeber: FAST Institut für Fahrzeugsystemtechnik


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.









**Band 88** Alexander Bernath **Numerical prediction of curing and process-induced distortion of composite structures.** 2021 ISBN 978-3-7315-1063-5

Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.

## Karlsruher Schriftenreihe Fahrzeugsystemtechnik

Fiber-reinforced materials offer a huge potential for lightweight design of load-bearing structures. However, high-volume production of such parts is still a challenge in terms of cost efficiency and competitiveness. This is particularly true for the automotive sector because of its enormous cost pressure. Besides high material cost of carbon fiber, cycle times are much longer compared to conventional metal based manufacturing processes. A promising strategy to shorten the cycle time is the combined use of fast-curing resins and advanced injection processes like Pressure-Controlled Resin Transfer Molding. However, finding suitable process parameters is often difficult without extensive practical knowledge. Moreover, when using fast-curing resins, the shortened cycle time makes this even more difficult as the fault tolerance decreases. In addition, it is difficult to monitor and analyze the progress of individual process steps as most of them take place hidden in a closed mold. This is where numerical process simulation can help to provide a better insight and to analyze underlying mechanisms. In this study, the curing process of the resin is investigated with regard to its influence on mold filling and process-induced distortion.

A. Bernath

Numerical prediction of curing and PID of composite structures

**Band 88**