Numerical Analysis of Thermoplastically Welded CFRP Structures

Carbon fibre reinforced thermoplastic polymers (CFRTP) are of particular interest to the aerospace industry. The possibility of thermoplastic welding as a joining method makes CFRTP an enabler in the pre-installation of systemic functionalities and cabin elements. This can be achieved by dust-free joining. In this work, thermoplastically welded joining zones are numerically analysed and evaluated for their failure behaviour. Finite element models for the evaluation of the peel strength (L-pull) are defined. In particular, the respective beginning of the damage as well as the damage propagation within the thermoplastic joining zone are of interest in order to identify the critical regions and to derive possibilities for design improvements.

of particular interest because thermoplastic welding can establish fast and cost-effective joining processes. The mathematical description and evaluation of the failure behaviour of CFRTP structures and especially of the thermoplastic joining zones has been the subject of research for several decades. Starting with the relatively simple analytical approaches, such as equivalent stress criterion, maximum principal stress or strain criteria, the interpolation criterion offers for the first time the possibility to evaluate multiaxial loads in a single layer or in a complete laminate, [1]. By Puck [2] a further very accepted failure criterion was introduced, which allows the distinction between fibre and interfibre breakage by evaluating the current load in the fracture process zone. With the introduction of stress intensity factors by Irwin [3], cracks in mechanical structures were characterized for the first time. The differentiation of the different fracture modes based on stress, planar shear and non-planar shear finally led to the definition of three different stress intensity factors. At the same time, this can be seen as a starting point for linear elastic fracture mechanics (LEFM). The disadvantage of all mathematical models from the field of LEFM must be seen in the stress singularities at the respective crack tip. It is known that a plastic zone is formed at the crack tip, which limits the resulting stresses. The description of ductile material behaviour at the crack tip was introduced by Wells [4] through the crack tip opening displacement concept (CTOD). For ductile materials it is assumed that the crack behaviour depends only on the plastic deformations at the crack tip. The concept can be seen as the starting point for Elastic Plastic Fracture Mechanics (EPFM). Dugdale [5] and Barenblatt [6] were the first to implement the idea of further developing the crack tip into a fracture process zone in order to avoid stress singularities. By doing so, cohesive forces act on the fracture process zone until the crack opening has reached a critical value. The cohesive forces then degrade until the crack surfaces are finally stress-free and the crack can spread further. At the same time the concept of the Cohesive Zone Method was born. The description of the formation of discontinuities withhile discrete interface elements was mainly driven by Needleman [7], Xu and Needleman [8,9] and Camacho and Ortiz [10]. For an overview of the numerical implementation of the cohesive zone approach within the framework of the finite element method, this paper mainly gives [9,[11][12][13][14][15][16].
Thermoplastic matrix systems enable the application of alternative joining processes. In the aviation industry, where mechanical riveting processes or duroplastic bonding processes are primarily used today, thermoplastic matrix systems also permit the use of welding technologies or fusion bonding technologies. The following disadvantages must be considered for mechanical riveting and thermosetting bonding processes: While stress concentrations occur at the rivet holes during mechanical riveting, complex surface preparation is usually indispensable for the thermosetting bonding of structural components. In contrast, thermoplastic welding processes, such as electric resistance welding, ultrasonic welding or electromagnetic induction welding, offer promising alternatives to these joining processes, [17] and the references therein. The continuous further development and optimization of thermoplastic welding processes is an important component today, since the demand for thermoplastic composite materials in the aerospace industry is also constantly growing in comparison with metallic and thermoset materials, in order to be able to better withstand the static and fatigue loads in the structure.
In Hoppe [18], a numerical crack simulation of fibre-matrix debonding in single fibre pullout tests was carried out using finite element method and a peridynamic model approach to determine the interface parameters of the model. They have shown, that the traction separation law parameters of Mode II and the free fibre have the most important influence on the pull-out test results in finite element simulation, whereas in peridynamics, the energy release rate as well as the horizon radius were found to be the most important parameters in addition to the free fibre length. The work of Gribanov et al. [19] is dedicated to efficiently implement the Park-Paulino-Roesler cohesive zone model (PPR-model) for an efficient simulation of material deformation, fracture and post fracture behaviour compliant, brittle and ductile materials. They developed a library for modelling polycrystalline materials under various types of loading conditions. The PPR-model therefor is able to distinguish between normal and tangential fracture modes at consistent fracture energy. Fang et al. [20] developed a special dual Langrangian multiplier algorithm to account for cracks with quasibrittle behaviour by avoiding stress oscillations and ill-conditioned systems matrices. To ensure the well posedness of the matrix, the proposed dual Langrange multiplier method uses the area-weighted average and biorthogonality conditions to define the contact constraints. The efficiency of solving the linear equations was shown to be more than 10 times increased. In de Brauer et al. [21], a level-set based method was functioned to simulate the evolution of damage in ductile materials under high velocity impact conditions using ductile fracture in bulk material and interfacial debonding at material interfaces. The level-set method therefor is used to create the crack around the fully damaged zone. Chowdhury and Wu [22] performed a cohesive zone based finite element analysis to examine the effects of processing parameter, i.e. clay nanoparticle volume fraction and aspect ratio. The showed that this approach is feasible to achieve high load-carrying capacity when increasing clay nanoparticle volume fraction, aspect ratio or alignment of the particles. Bian et al. [23] proposed a phase-field based cohesive zone model to account for cracking simulation. This approach is able to not only account for bulk cracks but also for interfacial cracks. Numerical examples were carried out to verify the cohesive model. Furthermore, they have shown that the interfacial mechanical properties were not affected by the bulk material properties. In Bayat et al. [24], a modelling approach for brittle failure of interface in case of linear elastic material is presented by integration of an extrinsic cohesive zone model into the incomplete interior penalty Galerkin variant of the discontinuous Galerkin method. A reduced integration scheme was employed on the boundary terms of the underlying equations to prevent any artificial stiffening within the bulk. Furthermore, by using hanging nodes on the discontinuities and interfaces, the calculation costs could be reduced. Seon et al. [25] introduced a new finite element-based method to simulate the debulking by using pore-pressure cohesive elements (PPCE) at ply interfaces. The PPCE therefor are used to consider entrapped air pockets as well as to model the air flow during debulking in prepreg materials. Experimental investigations on two-ply carbon/ epoxy specimen have proven the suitability of the numerical approach. The work of Russo and Chen [26] focuses on the problem that a very fine finite element mesh is required to resolve the cohesive zone with sufficient accuracy. They proposed a cohesive element-based method by use of slender structural elements for the plies to model delamination in composites. The method was implemented in 2D and validated on various models, i.e. DCB 1 , ENF 2 , MMB 3 and buckling-driven delamination and double-delamination problems. The proposed approach showed excellent agreement with the reference solutions. In standard cohesive zone modelling, the crack paths must be determined beforehand. To overcome this issue, the work of Roth and Kiefer [27] introduces an approach that allows for embedding cyclic cohesive laws into phase-field modelling. The results show that the phase-field approach is able to achieve the same results as the discrete formulation using cohesive element. The calculation time with this approach is up to 28 times higher than with the discrete formulation. However, the cohesive zone can be easily adjusted and changed without the need for extensive remeshing. A new approach to simulate the delamination initiation under cyclic loading conditions is in focus of [28]. They used the hysteresis cohesive zone modelling in conjunction with gradual degradation of the interface properties. The numerical results showed a good agreement with experimental findings taken from literature. Furthermore, they point out that damage initiation is caused by gradual degradation of the interface. Therefor, stiffness reduction occurs before damage initiation as a result of fatigue loading. Li et al. [29] introduced the shifted fracture method. For this, they combined an approximate fracture geometry representation in conjunction with approximate interface conditions. The advantage of his method is that the fracture behaviour is imitated with standard integrals on the approximated fracture surface. Therefor, they achieve accuracy comparable to XFEM/GFEM approaches at much lower computational complexity.
In this work the out-of-plane peel failure as well as the shear failure of thermoplastically welded joints are numerically investigated and evaluated. For this purpose different finite element models are constructed and provided with cohesive contacts at the corresponding joints. Particular emphasis was placed on the qualitatively correct mapping of the beginning of the damage and the propagation of the damage within the joining zones. The Limit of effectiveness is introduced to enable a simple and quick evaluation of the residual strength behaviour of L-pull specimens with filler material compared to L-pull specimens without filler material. The numerical investigations show very reasonable results. This leads the authors to the conclusion that the numerical models described here can be regarded as verified. A final and necessary validation of the models is part of further research tasks.

L-pull Sub-components
The L-pull sub-components are used to determine the maximum peel forces of the welded joint. In order to investigate and evaluate the quality of the thermoplastically welded joints, two different configurations A and B were selected. In configuration A, the roller bars are positioned parallel to the base of the L-angle. In addition, in configuration A, the distance between the roller bars was chosen so that one bar sits directly on the angled base with a roller bar distance of 50 mm (L50 model) and on the base plate 4 completely with a roller bar distance of 80 mm (L80 model) on the other hand. In configuration B 5 , these are perpendicular to the angle base at a distance of 130 mm (L130 model).
Both the L-angle and the corresponding base plate consist of 12 UD layers of LM-PAEK 6 material. This is a thermoplastic material resistant to high temperatures with a single layer thickness of 0.18333 mm. The entire symmetrical laminate structure was formed from two auxiliary laminates, each consisting of 6 individual layers. The stacking sequences of the two laminates were a = [+45/135/90/0/135/+45] and b = [+45/135/0/90/135/+45], respectively. The stacking sequence itself is based on the laminate structure of current and future CFRP shells and structures in aircraft construction.
It is known that the loads on L-pull specimens cause large notch effects in the area of the angular connection to the base plate. This significantly reduces the maximum loads that can be transferred. These notch effects can be reduced and the transmissible forces increased by using an appropriate filler in the area of the radius of the L-angle. For this reason, models of the L-pull sub-components were also defined, which contain a corresponding filler in the radius range, Figs. 1, 2 and 3. The filler itself is made of LM-PAEK

The Material Model
Linear elastic material behaviour is assumed for the base material. We assume that an elastic potential exists for the CFRTP material, so that the main symmetry can be assumed for the stiffness tetrad. Due to the corresponding symmetry properties of the stress tensor and the strain tensor , we also assume the left and right sub-symmetry for the stiffness tetrad. We also assume transversal isotropic material behaviour for a CFRTP single layer. This means that 5 independent material constants exist in the material law, Eq. 2 7 .
The engineering constants given in the table below which were taken from the literature are assumed for the CFRTP base material Table 1. where the following relationships exist between these engineering constants and the stiffness components in the stiffness tetrad  7 The material law is given in Voigt notation, i.e. (1) The density of the material was assumed to be Linear elastic isotropic material behaviour was defined for the filler material. The defined engineering constants are given in Table 2.

The Damage Model
The damage to the material was realized using cohesive zone models. Basically, there are two ways to work with cohesive zones in AbAqus. The direct use of cohesive elements allows the control of interface properties as well as the use of different failure modes. Furthermore, the interface has a finite thickness. A disadvantage of this modelling technique has to be seen in the fact that the complete interface has to be modelled a priori, which can be very complex. Cohesive elements are very thin and stiff. For this reason, they require very small and stable time increments when using an explicit solver. In addition, a mass definition is required for the element definition.
As an alternative to the cohesive elements, it is also possible to work with cohesive surfaces. The thickness of the interface is zero. The cohesive interface can be added to the model at any time by defining a corresponding contact. Disadvantages of this technique are that the user has no control over the interface properties. Furthermore, only one type of failure can be used. On the other hand, there is the advantage that no mass definition is necessary for the material definition of the cohesive surface. Compared to cohesive elements, the use of cohesive surfaces requires a sufficient mesh fineness to adequately map the failure behaviour. Due to the simplicity of the application and the possibility of being able to adequately depict the damage behaviour when the L-pull specimens fail, cohesive contacts are used in this work to define the failure behaviour. The beginning of the damage is evaluated using the Quadratic Traction (QUADS) law predefined in AbAqus where t n , t s and t t denote the actual traction values in normal, transverse as well as shear direction, and t 0 n = 30 MPa, t 0 s = 70 MPa and t 0 t = 70 MPa are the critical tractions in normal, transverse and shear direction. The values were chosen in close consultation with the industrial project partner. This is a stress-based criterion which calculates the stress ratios of the current stress state with the respective strengths. If the sum of the quadratic stress ratios is 1, the start of the damage is assumed.
In order to be able to consider mixed mode loads in the damage evolution, the power law evolution criterion based on energy release rates also predefined in AbAqus is used for this purpose where G I , G II and G III denote the actual energy release rates and G IC = 0, 5 N/mm, G IIC = 2 N/mm and G IIIC = 2 N/mm are the critical energy release rates in mode 1, mode 2 and mode 3, respectively. Again, the values were chosen in close consultation with the industrial project partner. was chosen to equal 1 to account for a linear influence of the energy release rates onto the damage.
The principle sequence of damage propagation using a cohesive zone approach can be seen in Fig. 4. At the beginning of the crack tip, the applied external forces finally induce loads, so-called tractions, which cause further cracking of the crack. The cohesive zone is still able to absorb further loads, which can be seen in a further increase of the traction above the respective separation. If the crack opening has reached a critical value c , the maximum of the load transfer in the cohesive is also reached. From this point on, a further load at the free edges of the cohesive zone finally leads to a successive loss of traction. Thus, this area is also called strain-softening branch. If the separation has finally reached the final value f , then the final failure of the cohesive zone is reached and a complete separation of the former contact partners is present.

Finite Element Models
For the L-pull models the respective UD layers were modelled as a 3D continuum shell. The single layer thickness was 0.1833 mm each. For each UD layer a local coordinate system was introduced to define the fibre angles. The most vulnerable regions of the respective models were determined by preceding sensitivity analyses. I.e. in order to shorten the calculation time and modelling time, numerical preliminary investigations were carried out. Based on these calculations, the damage behaviour of L-pull samples was assessed numerically. These investigations have shown that the maximum strains and the delamination only occur in certain areas. In these regions, cohesive contacts were defined, which allow a damage analysis. In the other areas of the models, the individual layers were linked via tie contacts, Fig. 5. Furthermore, between the gripper elements and the L-pull model, a surface to surface tie contact condition was used. The linear and reduced integrated 8-node element C3D8R was used for meshing. Appropriate stabilization techniques were used to avoid hourglassing effects. By application of the CourAnt-FriedriChs-Levi criterion, the maximum permissible time step is finally determined at where l char denotes the characteristic element length, Δt is the smallest time increment and c ist the wave propagation speed, respectively. The following Table 3 shows the defined element sizes for the L-pull models for the areas of the cohesive contact and the tie contact. The L-pull models are sensitive to the load application (velocity/angular velocity). For this reason, a smooth step was introduced for the load introduction, Fig. 6.
The failure of L-angles is a highly dynamic process that must consider geometric and material non-linearities. For this reason, the explicit AbAqus solver with mass scaling was used in the present work for the solution of the FE models. The stable time steps could be determined from previous investigations. For the L-pull models the stable time step is t=5E-8 s.

Numerical Results
In the following pictures the damage behaviour of the L-pull models with and without filler at different times of the experiment is shown. It can be seen that the filler leads to a significantly improved failure behaviour of the sample. Compared to the model without filler, the interfacial crack starts much later. In addition, the filler prevents a too large notch effect in the specimen, which leads to early failure. If one looks at the reaction forces achieved by both models, it can be seen that the model with filler has about 150% more load capacity. The very positive effect of the filler on the overall load-bearing behaviour has already been qualitatively proven.

L50 Model
In the L50 model a very early marking of the L-angle from the baseplate in the case without filler could be observed, Fig. 7. This means that significant damage to the L-angle radius occurs at a very early stage of loading. Accordingly, the damage continues to spread very rapidly in the further course, so that the final failure of this model can also be observed at a very early stage. The onset of damage for the models with and without filler is 0.008 s. At t=0.01066 s, the model without filler already shows a strong detachment of the L-angle from the base plate. Instead, in the model with filler, only half of the filler has come off. At t=0.012 s, more than half of the bracket in the model without filler has already detached from the base plate, whereas in the model with filler, only the filler has detached a little further. The connection of the L-angle is still completely intact. At t=0.02 s, the bracket has completely detached from the base plate in the model without filler. In the model with filler, more than half of the connection between the L-angle and the base plate is still intact. This behaviour clearly speaks for the use of filler to increase the loadbearing capacity and to reduce possible notch effects. This can also be seen in the course of the reaction force curve.
There is a large residual strength in the model with filler, which is noticeable in a second increase in the reaction force curve. After that, the damage in the model expands so that at about t=0.016 s another local minimum becomes visible in the course of the reaction force. Only then does the reaction force continue to increase until the final failure is recorded  Fig. 8. Furthermore, this model still has about 75% of its maximum load capacity when the model without filler has already reached its final failure.

L80 Model
In the case of the L80 model, a very early scribing of the sample without filler is shown again, Fig. 9. Accordingly, damage propagation caused by further delamination occurs very early. In contrast to the L50 model, the initiation of damage in the L80 model with and without filler is at t=0.02 s. At t=0.03 s, the L-angle has already detached up to half of the base plate in the model without filler. Whereas in the model with filler, only the filler begins to detach from the baseplate. Here, a failure behaviour that is qualitatively comparable to the L50 model can be observed. At t=0.04 s, the angle of the model without filler has already completely detached from the base plate. The connection of the angle to the base plate in the model with filler is still completely intact here, until the load maximum is reached at t=0.051 s. However, a very strong load bearing capacity exists for this model, which is confirmed in the course of the reaction force curve. The roller bar distance of l=80 mm ensures that there is no direct contact between the foot of the L-angle and the roller bar. For this reason, the model without filler fails much earlier in this case, whereas the model with filler has a much later failure due to the much smoother load transfer into the L-angle. It is obvious that the filler improves the load transfer, as well. For this reason, the use of filler is of great advantage when using L-angles. If the model has already reached its final failure without filler, the model with filler only reaches its maximum load bearing capacity here. Overall, the load-bearing capacity of the model with filler is increased by a factor of about 2.5 compared to the model without filler, Fig. 10.

L130 Model
The L130 model shows both without and with filler very early scribing of the samples, Fig. 11. In the further course of the load, the respective delamination spreads rapidly, so that both samples fail almost simultaneously. Overall, the load bearing capacity of the Fig. 13 Comparison of all reaction force curves with and w/o filler. One notices the largest influence of the filler material in case of the L50 model, followed by the L80 model. Less influence of the filler material can be seen in case of the L130 model. The reaction forces are normalized by the maximum reaction force of the L50 model with filler. The diagram also shows the limit of effectiveness. It can be said that all curves are above the limit of effectiveness, which means that a positive effect of using a filler can be seen for all three configurations. If one of the curves were below the limit of effectiveness, this would mean a reduction in the load-bearing capacity of the structure due to the use of a filler model with filler is not so dominant compared to the previous models with L=50 mm and L=80 mm roller bar distance respectively. Rather, it shows that the model with filler has a maximum load capacity increased by a factor of only 1.6 compared to the model without filler, Fig. 12. The authors assume that this property is due to the significantly increased roller bar distance on the one hand and to the 90 • rotated arrangement of the roller bars compared to the L50 models and L80 model on the other hand. With this model, the positive influence of the filler can therefore be regarded as significantly smaller. A general comparison of the reaction force curves with and without filler material for all cases investigated in given in Fig. 13, where the limit of effectiveness is introduced to allow the possibility of a quick comparison and evaluation of the quality of the filler in the different geometric configurations studied. The reaction forces are normalized by the maximum reaction force of the L50 model with filler. A summary and evaluation of the numerical results obtained is given in Table 4 5

Conclusion and Further Work
In this paper three different models were defined to investigate the tear-off behaviour of L-angle samples with and without filler. The L-angle samples consist of 12 layers of carbon fibre reinforced thermoplastic polymer with a quasi-isotropic layer structure. The first two models have roller bar distances of L=50 mm and L=80 mm, respectively, whereby the roller bars were aligned in 0 • direction of a single layer. The third model examined had a roller bar distance of L=130 mm, whereby the orientation of the roller bars in this model was rotated by 90 • . For all cases examined numerically, it could be proven that when using filler material to introduce and transfer the structural loads, a significantly greater load could be transferred than in the case without the filler material. The significantly smoother load introduction into the samples ensures a significantly later initial failure combined with significantly higher residual strength. Accordingly, the use of filler material to connect L-profiles to thermoplastic CFRP base plates is to be preferred to the use of L-angles on base plates without filler material. In case of L50 model, a significantly earlier detachment of the L-angle without filler material compared to the model with filler was observed. The models with L=50 mm and L=80 mm roller bar distance showed in the case with filler a clearly improved load bearing capacity compared to these models without filler. The filler increased the maximum load by a factor of 2.4 and 2.5 respectively. The model with L=130 mm roller bar distance behaved somewhat differently compared to the first two models. Both the model with filler and the model without filler showed a very early marking of the L-angle from the base plate. The model with filler showed a significantly improved load bearing capacity compared to the model without filler. However, with a factor of 1.6 this is significantly lower than with the L50 models and L80 model. In addition, the final failure of both models with L=130 mm roller bar distance was achieved at almost the same time, which is also a difference compared to the L50 models and L80 model. The L80 model showed a significantly earlier detachment of the L-angle without filler material, as well. Compared to the L50 model, the L80 model with filler material showed a softer behaviour due to the greater roller bar distance. Looking at the model with L=130 mm roller bar distance, a simultaneous detachment of the L-angles was observed for the models with and without filler material. Furthermore, delamination within the base plates as well as in the L-angles could be identified. In contrast to the L50 models and L80 models with and without filler material, the models with L=130 mm roller bar distance revealed the lowest reaction forces.
In order to allow a quick and relatively simple assessment of the quality of fillers in relation to the residual forces to be achieved, the size of the limit of effectiveness was introduced in this work. In the present work, this variable clearly shows that the use of filler always has an advantage in terms of the achievable residual forces. The L80 model showed the greatest potential, even if the residual forces achieved in absolute terms turned out to be somewhat lower than in the case of the L50 model.
From a numerical point of view, it can be said that the models and the results obtained are very reasonable and verified. For a final validation of the models, future experimental investigations must be carried out to identify the material parameters of the material used. On the other hand, experimental L-angle tear-off tests must be carried out in order to be able to quantitatively confirm the qualitatively achieved numerical results.