The enhanced properties of continuous carbon fiber-reinforced polymers, such as high stiffness, high strength, and a good stiffness-to-weight ratio, can be used in additive manufactured composite structures (Gebhardt et al., 2019). Various surveys and scientific publications on the performance of these highly optimized structures can be found based on the available printing technology from companies such as Markforged, Anisoprint, 2021, Arevo, 2021, and 9 T Labs to manufacture thermoplastic composite parts (Gregory and Gonzdz 2016; Blok et al., 2018; Goh et al., 2018; Iragi et al., 2019; Borowski et al., 2021).

These companies use a modified version of the FFF method, which combines polymer and pre-impregnated carbon fiber printing. Because of its low cost and ease of use, the FFF method is one of the most widely used processes in additive manufacturing. Pre-impregnated fiber-reinforced materials are used in continuous fiber additive manufacturing, similarly to conventional composite manufacturing methods such as automated fiber placement, automated tape laying, or custom fiber placement. The disadvantage of any additive manufacturing method is that, as the print speed increases, the level of porosity increases, and the requirements for a load-bearing part are no longer fulfilled. Therefore, a consolidation step was required. To address this problem, various approaches such as heat treatment during printing (Bhandari et al., 2019) and annealing or compaction during post-processing can be used (Meng et al., 2019; Zhang et al., 2020). Compaction is the subject of two developments. Using additional rollers to integrate a compaction unit into the printing process. This is the state-of-the-art automated tape-laying method. However, the part can be consolidated using a press and solid mold. The second approach also allows for the addition of more complexity to the part. Eguemann et al. (2013) combined several 3D printed preforms and metallic inserts into a single part using consolidation. The Airbus Helicopters former Eurocopter Germany GmbH application case is a machined steel hinge from rotorcraft model EC 135, weighing 135 g and mounted on a carbon fiber-reinforced polymer door structure. Eguemann et al. (2013) demonstrated the load case, and various composite-related manufacturing processes, such as injection molding, compression molding, and hand lamination, were evaluated and compared to the steel baseline.

This study demonstrates a novel method for simulating the additive fusion (consolidation) process of a complex-shaped composite part to investigate the dependence of the final mechanical properties and residual stresses on the process conditions. The method is presented by the authors as a sequential thermo-mechanical coupled, transient implicit analysis in ANSYS R2022 based on user material subroutines (Grieder et al., 2022). The phase transition behavior of the polymer from solid to molten and back to solid was modeled. Temperature, fiber volume content, porosity, and crystallization all influence the engineering properties of composite materials. The validation of the model for process-induced deformations and porosity in an aerospace part, namely, the helicopter hinge and a bracket to assemble the helicopter door, is the focus of this study. Figure 1 depicts the composite parts under consideration and the manufacturing workflow.

Wijskamp (2005) presented a sequentially coupled thermo-mechanical transient simulation approach that couples relative crystallization with phase transitions from the molten to the solid phase. It is possible to accurately calculate process-induced deformations by making the stiffness during a mechanical simulation dependent on the phase change by using relative crystallization as a status variable. Wijskamp (2005) assumed that when a material’s relative crystallization exceeds 50%, it starts to behave as a solid. This transition point is known as the solidification point, and it is critical for calculating process-induced deformations because after a material passes this point, thermal loads caused by the cooling process can cause deformations.

This approach has also been used in the literature by researchers such as Wijskamp (2005) and Brauner et al. (2014) to analyze process-induced deformation during the thermoforming process of composites. The Nakamura model is used in both models (Nakamura et al., 1972; Nakamura et al., 1973). The cooling phase crystallization kinetics are demonstrated by ∂ θ ∂ t = n ∙ K T ∙ 1 − θ − ln 1 − θ n − 1 n (1) where T ,   θ , and K T represent the temperature, degree of crystallization, and Avrami coefficient, respectively. K T = ln 2 t 1 / 2 n 1 n (2) where t 1 / 2 represents the half-crystallization time: t 1 / 2 = a ∙ ∂ T ∂ t b + c (3)

Based on the experimental data, the parameters n , a , b , c were defined by fitting the crystallization vs temperature dependencies for various cooling rates.

Similarly, melting behavior can be represented by a phase change from 1 (solid state) to 0 (molten state) according to Wijskamp (2005), Greco and Maffezzoli (2003), Brenken et al. (2018). ∂ θ ∂ T = k m b ∙ e − k m b T − T m 1 + d − 1 e − k m b T − T m d 1 − d (4)

The values of the coefficients n , a , b , c from Eqs 1–3 and k m b and d from Eq. 4 for this study are presented in the following section.

Using the relative crystallization and the definition of the solidification point, the Young’s modulus of the polymer was discovered to be temperature dependent during the process. This approach is considered the standard for thermoset-related processes that use the degree of cure. The resin modulus can be linear or nonlinear in its dependence (Wijskamp 2005; Brauner et al., 2011). The polymer modulus can also be used to determine the engineering properties of a unidirectional composite ply. Schürmann (2007) and Halpin and Kardos (1976) provided different homogenization approaches for composite materials, and Brauner (2013) provided an overview. The process-induced deformation and residual stress can be determined using this incremental elastic approach (Zobeiry 2006). These methods are capable of accurately simulating process-induced deformation stresses. However, residual stresses could be overestimated because stress relaxation was not considered. However, the effect of relaxation is irrelevant for thermoplastic composites that are additively manufactured because Young’s modulus is very low compared to that of a thermoset resin.

The transition between the solid and molten material states used in this study is nonstepwise and nonlinear, which improves numerical robustness. This transition was calculated in Grieder et al. (2022) using dynamic mechanical analysis (DMA) and rheometer measurements. This is illustrated in Supplementary Figure S1.

The transverse isotropic mechanical properties were evaluated in this study using the homogenization approach from Halpin-Tsai (Halpin and Kardos 1976; Schürmann 2007). The elastic modulus in the fiber direction E 1 is defined as follows (Halpin and Kardos 1976): E 1 = E f , 1 φ + E m 1 − φ (5)

The elastic moduli in the transverse and thickness directions are defined as (Schürmann 2007) E 2 = E 3 = E m 1 + ζ 1 ∙ φ E f , 2 − E m / E f , 2 + ζ 1 ∙ E m 1 − E f , 2 − E m / E f , 2 + ζ 1 ∙ E m φ (6)

In this case, ζ 1 = 3.3 represents the model parameter fit to the experimental data for carbon fiber-reinforced polyamide 12 (PA12-CF) to provide a modulus value close to the measurement results in comparison to the present study (Grieder et al., 2022). v 12 = v 13 = v 23 = v f φ + v m 1 − φ (7)

Poisson’s ratio is defined as (Grieder et al., 2022).

Shear moduli are defined as follows (Halpin and Kardos 1976; Schürmann 2007): G m = E m 2 1 + v m (8) G 12 = G 13 = G m 1 + ζ 2 ∙ φ ∙ G f , 12 − G m G f , 12 + ζ 2 ∙ G m   1 − φ ∙ G f , 12 − G m G f , 12 + ζ 2 ∙ G m   (9) G 23 = E 2 2 1 + v 23 (10)

In this case, ζ 2 = 1.675 represents the model parameter fitted to the experimental data for PA12-CF to obtain a modulus value that is close to the measurement results (Grieder et al., 2022).

In the solid state, the mixing rules (5–10) are applied to the PA12-CF filament, whereas in the molten state, the material is considered isotropic, with material properties corresponding to pure polyamide 12 (PA12).

The mixing rules for orthotropic coefficients of thermal expansion (CTE) and crystallization shrinkage could be introduced to the model. The coefficients of thermal expansion in the transverse and longitudinal directions are calculated using the following equations (Schürmann 2007): α 1 = S 1 α m ⋅ E m ⋅ 1 − φ + α f , 1 ⋅ E f , 1 ⋅ φ E m ⋅ 1 − φ + E f , 1 ⋅ φ (11) α 2 = S 2 φ ⋅ α f , 2 + 1 − φ ⋅ α m (12) α 3 = S 3 α 2 (13)

Here, α m , α f , 1 and α f , 2 represent the longitudinal and transverse coefficients of thermal expansion for the matrix and the fiber, respectively. To account for possible changes in the CTE after reaching the glass transition temperature, the correction coefficients S 1,2,3 are introduced into the model. The coefficient values used in this study are presented in the following section.

Shrinkage is a crucial consolidation process-induced effect that has a significant impact on the final warpage of a part. The shrinkage of the PA12 composite was shown to be temperature and printed lay-up dependent (Negi and Rajesh 2016; Li et al., 2020). The shrinkage of composite materials is defined via the following mixing rules (Wijskamp 2005; Schürmann 2007): β 1 = 1 − φ ∙ E m ∙ β / φ ∙ E f , 1 + 1 − φ E m (14) β 2 = β 3 = 1 − φ ∙ 1 + v m ∙ β − v 13 ∙ φ + v m ∙ 1 − φ β ∙ 1 − φ ∙ E m φ ∙ E f , 1 + 1 − φ ∙ E m (15) where β represents the total volumetric shrinkage of PA12. Because the value of volumetric shrinkage was not measured within the scope of the present study, it was based on the average value of 3% (Benedetti et al., 2019; Li et al., 2020).

For the investigated additive-manufactured composite materials, porosity exhibits the greatest impact. Blok et al. investigated porosity during FFF manufacturing for SFR and CFR (Blok et al., 2018). The void content exhibits a significant negative impact on the composite properties. A high void ratio reduces strength, durability, and fatigue resistance, increases moisture absorption, and causes composite properties to vary within one part (Cable 1991; Ghiorse 1993; Mahoor et al., 2018). Therefore, void content should be minimized in bear-loading applications.

Zingraff et al. (2005) and Zingraff (2004) addressed void-process-induced formation in the composite processing literature. Voids may form in autoclave processing when the resin shrinks during thermoplastic crystallization (Eom et al., 2001). Grieder et al. (2022) demonstrated that the level of porosity of additively manufactured composites can reach 20–30% for simple geometry. Leterrier and G’sell, (1994) demonstrated that free heating of the composite specimen could result in an increase of up to 60% in the initial void content. This indicates the need for a consolidation step, which must be incorporated into the process simulation method.

During solidification, specific pressure and temperature regimes can be used to avoid void formation (Long et al., 1995; Lundstrom and Gebart (1994). Eom et al. (2001) proposed process conditions that yield void-free composites while requiring a minimum amount of pressure (0.2 MPa). Specially designed molds can minimize void formation and process-induced deformation (Eguemann et al., 2013). In this study, a similar approach was developed and implemented in which an advanced mold design was provided to minimize part porosity and warpage.

Experimental and numerical studies on void formation mechanisms in mesostructures and the influence of pores on the mechanical properties of 3D-printed composites have been presented in the literature (Bellehumeur et al., 2004; Wang et al., 2016). Multiscale finite element modeling methods for studying the influence of micro- and mesoporosity and process parameters on the mechanical properties of composites are illustrated in Calneryte et al. (2018), Rodriguez et al. (2000), and Xue et al. (2019). The power law empirical relationship can express stiffness-porosity dependence (Jaroslav 1999). E = E 0 ∙ 1 − ϕ ϕ c S (16)

Here, E represents the effective Young’s modulus of the porous medium, ϕ represents the porosity, E 0 represents Young’s modulus of the fully consolidated solid material, and ϕ c represents the porosity at which the effective Young’s modulus becomes zero. The model parameter S is determined via the porous material structure, pore geometry, and other structural properties. For each material and structure, the value of parameter S must be defined (Wagh et al., 1991; Jaroslav 1999; Grieder et al., 2022).

The power-law porosity dependence (18) was successfully implemented in the finite element model for non-polymer materials with a relatively low initial porosity (Hardin and Beckermann 2007; Morrissey and Nakhla 2018). However, this definition necessitates experimental measurements to tune the stiffness of the porosity function. Moreover, most of the presented approaches only allow for robust numerical simulation for a small range of void content change, failing in the case of rapid porosity evolution from 20 to 30% to almost zero void content.

In this study, a similar approach is presented in which the stiffness is determined using the local porosity content via a power law. Temperature, crystallization degree, and porosity all affect the matrix elastic modulus E m . Therefore, the following empirical dependence was developed (Grieder et al., 2022): E m = E m , s ∙ θ K ∙ 1 − ϕ S + E m , b ∙ 1 − θ K ∙ 1 − ϕ M (17) Here, E m , b represents the molten polymer bulk modulus measured with a pressure evaluation study (Grieder et al., 2022), and E m , s T represents the temperature-dependent elastic modulus in the solid phase measured with DMA (Grieder et al., 2022); K ,   S , and M represent model design parameters that must be fitted to the experimental data depending on the material and composite layup. The values of these parameters relevant to this study are presented in Supplementary Table S1; ϕ represents the void content, which considers the influence of mechanical deformations and air diffusion into the PA12 matrix on porosity.

Composite material properties are defined according to the mixing rules (5)–(15), which determine the proportion of pure fiber to pure matrix properties based on the fiber volume ratio φ : CTE (11)–(13), crystallization shrinkage (14)–(15), and stiffness (17). The relationship between the fiber volume ratio and porosity is not well covered in the literature. However, assuming that the fiber volume ratio is proportional to the void content, Grieder et al. (2022) proposed the simplest linear dependence on porosity: φ = φ 0 1 − ϕ (18) where φ 0 represents the measured fiber volume ratio corresponding to a fully consolidated material with no void content.

Many phenomena occur concurrently during consolidation, including bulk compaction, intimate contact between adjacent plies, interlaminar adhesion, fiber deformation, movement, and molecular diffusion (Koerdt et al., 2022). These phenomena interact in complex ways and are influenced by process parameters such as time, temperature, and pressure. The resin pressure must be evaluated to model void formation. This can be accomplished by simplifications or Darcy’s law when modeling a flow problem. Grieder et al. (2022) developed a simplified porosity model based on the approach of Barari et al. (2020). In the model, porosity ϕ is represented as a dimensionless variable with a value between 0 and 1, where 0 indicates that the material is fully consolidated and one indicates that the volume of the entire considered finite element is empty. The initial porosity value is determined using the material and printing configuration.

According to Barari et al. (2020), the hydrostatic pressure, p h is calculated as follows: p h = P 0 ∙ ϕ 0 ϕ T T 0 1 − V d , ϕ > ϕ min 0 ,   ϕ ≤ ϕ min (19) where P 0 represents the initial (atmospheric) pressure, T 0 represents the initial (room) temperature, ϕ 0 represents the initial porosity, ϕ min = 0.001 represents the minimum considered porosity, and V d represents the relative volume of dissolved air.

The volume of dissolved air is a dimensionless variable with values in the range of 0–1, with 0 indicating that the air occupies 100% of its initial volume and one indicating that the air is completely dissolved in the surrounding matrix. Because it is assumed that air traps have spherical shapes, the volume can be defined as V d = 1 − R b 3 (20) R b represents a dimensionless variable that represents the air bubble radius and exhibits values in the range of 0–1, with 0 indicating that the air bubble is completely dissolved in the resin and one indicating that the bubble radius has returned to 100% of its initial value. Therefore, R 0 = 1 . Eq. 20 is a semi-analytical approach based on dimensionless dissolved volume and bubble radius definitions (Grieder et al., 2022).

The fundamental mechanisms of air bubble formation and reduction in thermoplastic composites remain unknown. Early models incorporated diffusion into the isothermal growth of spherical gas bubbles in a viscous fluid (Moris and Costel 1984). Further studies expanded the model to include many bubbles in highly viscous fluids (Arefmanesh et al., 1990; Arefmanesh 1991). These models were then applied to thermoplastic composites for structural applications (Roychowdhury 1995; Roychowdhury et al., 2001; Yohann et al., 2008; Jin et al., 2020).

Gas diffusion into molten polymer was considered in the literature (Wood and Bader 1994; Zingraff 2004; Zingraff et al., 2005), where the diffusion coefficient for the polymer material was measured. However, the dimensionless diffusion coefficient in the present approach represents not only the diffusion and dissolution processes but also the air evacuation through the void channels along the filament’s printing path. Therefore, the diffusion coefficient D is a model design parameter that cannot be directly measured but must be adjusted for each new material, composite structure, and printing condition (Grieder et al., 2022).

In Grieder et al. (2022), it was assumed that every finite element contains an air bubble that evolves independently of the surrounding medium and other bubbles. The relative bubble radius was calculated as follows (Wood and Bader 1994; Advani and Sozer 2010): R b = R 0 − D ∙ p h − p h , l 0 ∙ t 0.5 (21) where D represents a dimensionless diffusion model parameter and p h , l 0 represents the hydrostatic pressures when melting occurs.

A general description of hydrodynamics is required for polymer melts (Joel et al., 1997). This description considers viscoelasticity, die swelling, and shear thinning. These properties necessitate nontrivial modifications to the hydrodynamic formulation of the fluid, resulting in advanced constitutive equations that significantly increase the complexity of the problem and computational effort. Because of the limitations of the finite-element software ANSYS and the aforementioned issue, this study considers a simplified model of the squeeze flow.

The simplified orthotropic squeeze pressure p s q was only considered in the liquid state and was evaluated using the method demonstrated by Barari et al. (2020): p s q , i = η T 12 ∙ K i ∙ Λ i 2 ∙ ∂ ε m , i ∂ t ,   ϕ > ϕ min 0 ,       ϕ ≤ ϕ min   o r   i = 1 ,   i = 2,3 (22) where η T represents the viscosity of the molten matrix and Λ i represents the directional distance to the nearest void (defined by the printing setup). The voids are represented as channels in the printing direction due to the composite manufacturing method. Therefore, the distance in the printing direction Λ 1 to the closest void is assumed to be zero, implying that there is no squeeze flow in the printing direction. ε m 1 , ε m 2 , ε m 3 are the strains of the matrix material inside the fiber filament are evaluated using Schürmann (2007): ε m , i = ε i l m / l 0 + E m / E f , i ∙ 1 − l m / l 0 ,   i = 1,2,3 (23) where ε 1 , ε 2 , ε 3 represent the longitudinal, transverse, and thickness strains, respectively; E m represents the matrix elastic modulus; E f , i , i = 1,2,3 are the longitudinal, transverse, and thickness elastic moduli of the fiber ( E f , 2 = E f , 3 ), and l m / l 0 is defined as (Schürmann 2007) l m / l 0 = 1 − 4 π φ 0.5 (24) where φ represents the fiber volume ratio. The permeability in the transverse and thickness directions can be defined by simplifying the interlayer region into a rectangular duct as (Zingraff 2004; Advani and Sozer 2010) K 2 = a 2 12 1 − 192 ∙ a w ∙ π 5 ∑ i = 1,3,5 , … 9 tanh i ∙ π ∙ w 2 ∙ a i 5 (25) K 3 = w 2 12 1 − 192 ∙ w a ∙ π 5 ∑ i = 1,3,5 , … 9 tanh i ∙ π ∙ a 2 ∙ w i 5 (26) where a and w represent the height and width of the rectangular duct, respectively, evaluated according to the following: w = w 0 ∙ 1 + ε 2 (27) a = a 0 ∙ 1 + ε 3 (28) where w 0 and a 0 represent the width and height of the printed filament, respectively. The values of these parameters relative to the current study are presented in Grieder et al. (2022).

As a function of the bulk strain ε b u l k , porosity was calculated as follows (Barari et al., 2020): ϕ = ϕ 0 + ε b u l k ∙ 1 − D ∙ p h − p h , l 0 ∙ t 0.5 (29) ε b u l k = 1 + ε 1 ∙ 1 + ε 2 ∙ 1 + ε 3 − 1 (30) where ϕ 0 represents the initial porosity, ε b u l k represents the bulk strain, p h represents the hydrostatic void pressure, and p h , l 0 represents the hydrostatic pressure when melting occurs. D represents a dimensionless diffusion model parameter that considers air dissolution and air evacuation using open-air channels Grieder et al. (2022). The authors’ developed porosity approach is based on (Wood and Bader 1994; Zingraff 2004; Barari et al., 2020) and is described in detail in Grieder et al. (2022). Depending on the internal void hydrostatic pressure, amount of diffused volume, and squeeze flow, void reduction is assumed to occur only in the molten and transitional material states.

The physics of rapid bubble collapse and its impact on the surrounding viscous medium have been discussed in terms of fluid mechanics for liquid media (Flynn 1975; Brujan and Vogel 2006). The bubble collapse process was studied for the solid or molten material state in Jin et al. (2020), where it was demonstrated that violent bubble collapse causes significant changes in material properties. When a porous material is subjected to large deformation, the failure threshold is predicted not directly by the total stress applied but by a certain difference between the total compressive stress and pore pressure. Skempton (1960) investigated the effective stress for volumetric deformation. The pore pressure should be considered in the constitutive equations (Skempton 1960; Barari et al., 2020), but a simple calculation of the pore stress using the ideal gas law results in rapid growth of the internal void pressure during fast void collapse (Grieder et al., 2022), which causes convergence problems when using the finite element method. For a more accurate description of bubble collapse in polymer melt, advanced viscoelastic models are required (Gaudron et al., 2015). Because of the simplicity of the proposed hydrodynamic model, the evaluated void pressure was not considered in the constitutive equations in the current study. However, the pressure indirectly influences the process-induced stresses via stiffness, which is a function of porosity 17) and changes owing to the compaction degree and air diffusion. The stiffness-porosity relationship is controlled by the model design parameters S and M. Therefore, the model can be tuned to compensate for the influence of internal void pressure on the surrounding resin medium (Grieder et al., 2022).

This section introduces the material and provides detailed characterizations. The following experimental characterizations were used: differential scanning calorimetry (DSC) was used to determine the crystallization/melting behavior, and thermomechanical analysis (TMA) was used to measure the coefficient of thermal expansion. The computer was used to determine the final shape and the porosity content of the consolidated composite part.

The developed numerical approach was implemented using a user material subroutine in the commercial finite-element software ANSYS (v. 2022 R1). The layup was created using the SpaceClaim plugin Fibrify^® and then mapped to 3D finite elements. Data exchange with the using hierarchical data format five ACP.

The results of the developed numerical approach for samples A and B are presented in this section, as demonstrated in the previous section. The crystallization, process-induced deformations, residual stresses, porosity, and squeeze pressure results of the sequentially coupled thermal and mechanical solutions are presented. The model was validated using CT porosity and process-induced deformation analysis results.

In this study, the composite consolidation process was investigated experimentally and numerically. The initial and final porosities, CTEs, and crystallization and melting kinetics of the pure PA12 and composite PA12-CF materials were determined experimentally. On the basis of these findings, a mathematical model for the consolidation process was developed and numerically solved using the finite element method. By considering the orthotropic composite properties, the proposed finite-element model predicts the final process-induced deformations of the consolidated part. Temperature, crystallization, and porosity all influence these deformations. The proposed model can predict the final composition of a part with high accuracy. The predicted warpage followed the same trend as the observed warpage. However, their absolute values differed. This difference could be caused by various factors, including inaccuracy of the thermal simulation, limitations in the simplified porosity model, defined homogenization mixing rules, imprecision in the initial data provided to the model, and inaccuracy in the CT measurements. The developed approach includes a set of design variables that must be adjusted for each material and composite structure. Furthermore, due to the ACP fiber filament mapping method, the model is extremely sensitive to the finite element mesh. Nonetheless, the developed approach enables rapid simulation of the consolidation process considering multiple physical and structural effects. The final residual stresses and porosity are predicted by the model. This is an important step toward achieving “first time right” composite production. The proposed approach provides a digital model of the consolidation process, which reduces the number of expensive prototyping iterations. The combination of highly accurate 3D-printing and post-printing consolidation, as well as the ANSYS finite-element model and fiber filament layup design, allows for the transition to serial production of additively manufactured continuous fiber composite parts.

The authors intend to use a local initial porosity distribution that depends on the printing process parameters and layup in a future study. Furthermore, novel measurement techniques to investigate porosity and crystallization are required.