Numerical Investigation of Residual Stresses in Welded Thermoplastic CFRP Structures

Using thermoplastics as the matrix in carbon fiber-reinforced polymers (CFRP) offers the possibility to make use of welded joints, which results in weight savings compared to conventional joining methods using mechanical fasteners. In this paper, the resulting temperature distribution in the material due to resistance welding is investigated by transient finite element (FE) simulations. To examine the effects on the component structure, a numerical modeling approach is created, which allows determining the residual stresses caused by the welding process. It is shown that the area of the structure, especially near the joining zone, is highly affected by the process, especially in terms of residual stresses. In particular, the stresses perpendicular to the fiber direction show failure relevant values up to a maximum of 221 MPa, which might lead to the formation of microcracks in the matrix. In turn, that is assumed to be critical in terms of the fatigue of welded composite structures. Thus, the suggested modeling approach provides residual stresses that can be used to determine their effects on the strength, structural stability, and fatigue of such composite structures. In a subsequent step, these findings could play an important role in the design process of thermoplastic composite structures.


Introduction
The use of thermoplastic matrix materials in fiber-reinforced polymers is increasingly becoming the focus of applied aerospace research. Compared to epoxy resin systems predominantly used in the industrial environment today, there are far-reaching advantages when using thermoplastics. First of all, the possibility of in situ consolidation while using thermoplastics plays a decisive role, which results in significantly shorter cycle times in the production of such components. Secondly, the recyclability of thermoplastic materials is another significant factor [1]. Both aspects are to be considered concerning the possible production and resource efficiency that can be achieved and thus in turn represent key aspects with regard to environmentally friendly aviation.
Another benefit of thermoplastic materials is their weldability. Compared to conventional joining techniques using mechanical fasteners [2,3], significant weight advantages can be achieved through welding, which in turn represents a concrete action with respect to eco-lightweight design and, as a result of lower fuel consumption, leads to lower CO 2 emissions. Several methods have been investigated on an experimental and numerical level so far [4][5][6][7][8][9][10][11][12]. While in the paper at hand, the resistance welding method is examined and described in detail later, further relevant technologies will be described in the following as well. One of them is the induction welding method. For the induction welding of thermoplastics, the electromagnetic field of an external inductor heats the heating element 2 of 10 in the joining zone. This leads to the melting of the matrix and in consequence, the adherends can be joined under pressure. In [4], a procedure based on numerical methods has been developed, allowing a more precise understanding of the process parameters involved. The joining quality is evaluated based on single-lap-shear tests (SLS tests) and microscopic examination of the failure surfaces [5]. Another well-studied method is the laser welding of thermoplastics, for which a laser as a heat source is used. In [6], a thermal finite element (FE) model of a laser transmission welding process is developed using an unreinforced thermoplastic. This allows determining the size of the heat-affected zone and investigating the resulting internal stresses. Furthermore, a process model under consideration of the composite's microstructure was developed, which enables clarifying the correlation between fiber volume content and welding quality [7]. Familiarization with the literature shows that ultrasonic welding could also play an important role in future applications. For ultrasonic welding, the vibrations generated by an ultrasonic generator are converted into mechanical vibrations via sonotrode. The targeted introduction of these vibrations in the area of the joining zone leads to the melting of the matrix. With the help of experimental methods, it has so far been shown that this variant, which is usually implemented in a spot welding topology, could be a useful alternative to conventional rivets [8]. A process modeling procedure was shown in [9]. However, there still is no valid concept for how the ultrasonic welding process could be used for large continuous joints. According to the current state of research, one of the most promising variants, and therefore one of the best understood, is resistance welding. Here, a resistance element located in the joining zone is heated by current feed so that the matrix is melted. Subsequently, the adherends are joined under pressure as a result of cooling [10]. For this process, it was possible to demonstrate, mainly on an experimental level, the potential of its use in aircraft construction. For example, possible process windows could be identified, and the quality of the joint is determined by optical observation of the failure surfaces and SLS tests. In addition, the first numerical investigations were carried out [11,12].
All these joining processes have in common that a considerable amount of heat is introduced into the material during the process. As a result, different cooling rates within the material, different fiber orientations of the individual layers, and different coefficients of expansion of the fiber and matrix cause stresses to develop without external forces, so-called internal or residual stresses. These occur on different levels, which are briefly explained in the following. First of all, the residual stresses have to be considered on a micromechanical level. These occur primarily due to the different thermal expansion coefficients of fiber and matrix and thus arise at the fiber-matrix interface. The second level, which is given by the macro-mechanical relationships in the laminate, is also characterized by the occurrence of internal stresses. These occur in particular with different fiber orientations within the laminate and thus result in interlaminar shear forces. The third level of residual stress can be named global stresses. These become relevant in the case of relatively thick laminates, in which the high thickness leads to areas of different cooling speeds and thus forming stress gradients [13,14]. For the investigation of residual stresses, several methods are applicable. In [15], the nanoindentation test was used to measure the field distribution of the residual stress in functionally graded ceramic-metal composites. In the following step, a comparison of the results with corresponding finite element method (FEM) simulations showed a qualitative agreement.
The presented work examines the connection between the method of resistance welding and the occurrence of residual stresses on a macro-mechanical level. In particular, the high amount of thermal energy that is introduced into the laminate during the welding process inevitably leads to the development of internal stresses. First, a model of the process of resistance welding using FE-based methods to determine the temperature distributions responsible for the development of the residual stresses will be displayed. This allows the investigation of the relationship between the welding parameters and the thermal load of the structure. Moreover, the process model can be used as a basis for research on process parameter optimization. Furthermore, a numerical method that enables the resulting residual stresses to be determined on a mesoscopic level and evaluated for further investigations is presented thereupon. These results serve for the investigation of the welding's thermomechanical effect on the component concerning preloading and material damage. In addition, the residual stresses represent useful data for future experimental validation and comprehensive structural design.

Materials and Methods
The composite material used for the investigations in this paper is the unidirectional prepreg Aromatic Polymer Composite 2 (APC-2) certified for aircraft structures, consisting of high strength, standard modulus carbon fibers (AS-4) with a nominal fiber volume content of 61%, and the thermoplastic matrix material Polyether Ether Ketone (PEEK). The polymer has semi-crystalline properties with a glass transition temperature of 145 • C and a melting temperature of 390 • C. All calculations are carried out using temperaturedependent material properties taken from [16]. The elastic and thermal properties at room temperature used are given in Table 1. A simplified resistance welding process is modeled using the commercial Finite Element code ABAQUS ® . Two coupon specimens according to ASTM D1002 standard are used as adherends. To reduce the computational time for the modeling procedure, both adherends are assumed to behave similarly, so that only one coupon specimen is modeled with adequate boundary conditions. The schematic resistance welding process based on a typical experimental setup and the corresponding simplification for modeling is shown in Figure 1. resulting residual stresses to be determined on a mesoscopic level and evaluated for further investigations is presented thereupon. These results serve for the investigation of the welding's thermomechanical effect on the component concerning preloading and material damage. In addition, the residual stresses represent useful data for future experimental validation and comprehensive structural design.

Materials and Methods
The composite material used for the investigations in this paper is the unidirectional prepreg Aromatic Polymer Composite 2 (APC-2) certified for aircraft structures, consisting of high strength, standard modulus carbon fibers (AS-4) with a nominal fiber volume content of 61%, and the thermoplastic matrix material Polyether Ether Ketone (PEEK). The polymer has semi-crystalline properties with a glass transition temperature of 145 °C and a melting temperature of 390 °C. All calculations are carried out using temperature-dependent material properties taken from [16]. The elastic and thermal properties at room temperature used are given in Table 1.  A simplified resistance welding process is modeled using the commercial Finite Element code ABAQUS ® . Two coupon specimens according to ASTM D1002 standard are used as adherends. To reduce the computational time for the modeling procedure, both adherends are assumed to behave similarly, so that only one coupon specimen is modeled with adequate boundary conditions. The schematic resistance welding process based on a typical experimental setup and the corresponding simplification for modeling is shown in Figure 1. The part is modeled using quadratic hexahedral, coupled temperature-displacement elements. The geometry illustrated in Figure 1 is similarly implemented in the width direction to create a 3D model. No edge effects in the width direction are taken into account, as the considered specimen is assumed to be cut from a plate so that constant conditions can be applied. The heat generated during welding is applied directly to the joining zone of the coupon specimen and is simplified as a surface heat flux. It results from the joule heating effect of a resistive element (i.e., carbon fibers, steel mesh) and is assumed to be uniformly distributed at the overlap area. The pressure for consolidation is assumed to be The part is modeled using quadratic hexahedral, coupled temperature-displacement elements. The geometry illustrated in Figure 1 is similarly implemented in the width direction to create a 3D model. No edge effects in the width direction are taken into account, as the considered specimen is assumed to be cut from a plate so that constant conditions can be applied. The heat generated during welding is applied directly to the joining zone of the coupon specimen and is simplified as a surface heat flux. It results from the joule heating effect of a resistive element (i.e., carbon fibers, steel mesh) and is assumed to be uniformly distributed at the overlap area. The pressure for consolidation is assumed to be constant in length and width direction and is applied directly to the specimen. The thermal boundary conditions are represented by free convection due to ambient air temperature at the free surfaces, while the area between the specimen and the insulator is modeled by heat transfer from carbon fiber-reinforced polymer (CFRP) to wood, which is a typical insulator material in this context. The joining surface is assumed to be fixed in all degrees of freedom. The specimen consists of a [+45 • /0 • /−45 • /90 • ] 2S quasi-isotropic layup. Every single unidirectional layer is represented as one quadratic hexahedral coupled temperature-displacement element in the thickness direction. This results in an overall thickness of t = 2 mm. The specimen's length is l = 101.6 mm, and the width is w = 25.4 mm. The overlap length is assumed to be l lap = 25.4 mm. Furthermore, the global coordinate system definition ( Figure 2a) as well as the layup used for the simulation (Figure 2b) need to be addressed. While the global system is characterized by the x, y, and z coordinates, the local layer system is represented by the I-coordinate in the fiber direction, the II-coordinate in-plane perpendicular to the fiber direction, and the III-coordinate outof-plane perpendicular to the fiber direction.
J. Compos. Sci. 2021, 5, x FOR PEER REVIEW 4 of 10 constant in length and width direction and is applied directly to the specimen. The thermal boundary conditions are represented by free convection due to ambient air temperature at the free surfaces, while the area between the specimen and the insulator is modeled by heat transfer from carbon fiber-reinforced polymer (CFRP) to wood, which is a typical insulator material in this context. The joining surface is assumed to be fixed in all degrees of freedom. The specimen consists of a [+45°/0°/−45°/90°]2S quasi-isotropic layup. Every single unidirectional layer is represented as one quadratic hexahedral coupled temperature-displacement element in the thickness direction. This results in an overall thickness of t = 2 mm. The specimen's length is l = 101.6 mm, and the width is w = 25.4 mm. The overlap length is assumed to be l lap = 25.4 mm. Furthermore, the global coordinate system definition ( Figure 2a) as well as the layup used for the simulation (Figure 2b) need to be addressed. While the global system is characterized by the x, y, and z coordinates, the local layer system is represented by the I-coordinate in the fiber direction, the II-coordinate in-plane perpendicular to the fiber direction, and the III-coordinate out-of-plane perpendicular to the fiber direction. The investigation of the residual stresses in the components resulting from the cooling after the welding process is based on reduced elastic properties of the matrix at welding temperature. This approach requires the definition of consolidation levels depending on the elastic properties of the matrix which again depend on the processing temperature, thus discretizing the cooling process. Afterward, the matrix-dominated mechanical properties are updated according to the consolidation level. The resulting stress distribution represents the initial stress level in the composite material after welding and right before cooling, respectively. Due to the re-consolidation in the cooling process, heterogeneously distributed residual stresses in the component arise. Linear thermo-mechanical analyses are carried out, and the entire modeling procedure is divided into two steps. Firstly, a heating step based on the resistance welding simulation is carried out. After the melting temperature in the joining zone is surpassed, the second step of the procedure takes place. The input power is turned off or in other words, the component is cooling down due to free convection without external heat flux under pressure. The recalculated mechanical properties are used as the initial material property. The resulting stresses of the laminate are considered to be the residual stress states caused by the welding process. The investigation of the residual stresses in the components resulting from the cooling after the welding process is based on reduced elastic properties of the matrix at welding temperature. This approach requires the definition of consolidation levels depending on the elastic properties of the matrix which again depend on the processing temperature, thus discretizing the cooling process. Afterward, the matrix-dominated mechanical properties are updated according to the consolidation level. The resulting stress distribution represents the initial stress level in the composite material after welding and right before cooling, respectively. Due to the re-consolidation in the cooling process, heterogeneously distributed residual stresses in the component arise. Linear thermo-mechanical analyses are carried out, and the entire modeling procedure is divided into two steps. Firstly, a heating step based on the resistance welding simulation is carried out. After the melting temperature in the joining zone is surpassed, the second step of the procedure takes place. The input power is turned off or in other words, the component is cooling down due to free convection without external heat flux under pressure. The recalculated mechanical properties are used as the initial material property. The resulting stresses of the laminate are considered to be the residual stress states caused by the welding process.

Welding Simulation
As mentioned in the cited publications, the weld's quality is highly affected by the temperature magnitude and distribution in the joining zone. The performed welding process simulation leads to time-temperature data for the considered component so that the temperature of every single element at a certain time during heating and cooling can be assessed. In Figure 3, the temperature data of three selected measuring points (P1, P2, and P3) are presented. All three points are located in the mid-plane in the width direction of the specimen. It can be seen that the melting temperature of the matrix (T melt = 390 • C) is reached after t weld = 25 s in the middle of the joining zone (P1). Furthermore, it can be noticed that with this welding parameter configuration, only a section of the joining surface is melted, as the temperature in P2 does not reach the melting temperature. Afterward, the input power is modeled to be turned off, and the cooling of the specimen due to free convection begins. While there is a clear drop in temperature to observe for P1 and P2 after turning off the input power, a slight increase in temperature can be noticed for P3, which results from the heat flow in the material.
As mentioned in the cited publications, the weld's quality is highly affected by the temperature magnitude and distribution in the joining zone. The performed welding process simulation leads to time-temperature data for the considered component so that the temperature of every single element at a certain time during heating and cooling can be assessed. In Figure 3, the temperature data of three selected measuring points (P1, P2, and P3) are presented. All three points are located in the mid-plane in the width direction of the specimen. It can be seen that the melting temperature of the matrix (T melt = 390 °C) is reached after t weld = 25 s in the middle of the joining zone (P1). Furthermore, it can be noticed that with this welding parameter configuration, only a section of the joining surface is melted, as the temperature in P2 does not reach the melting temperature. Afterward, the input power is modeled to be turned off, and the cooling of the specimen due to free convection begins. While there is a clear drop in temperature to observe for P1 and P2 after turning off the input power, a slight increase in temperature can be noticed for P3, which results from the heat flow in the material. In Figure 4, the temperature distribution in the specimen during welding is shown. It can be noticed that the heat builds up in the joining zone and nearby areas of the specimen due to the local heat flow input and the insulator. Only with a certain delay, as it can be seen in Figure 3 as well, the heat is transported along the specimen's length direction. Although the heating element generates heat uniformly at the overlap area, the temperature distribution is not uniform. Instead, due to different heat transfer conditions between the edge and the middle of the joint, there is a maximum temperature at the middle and lower ones at the edges of the overlap. As a result of the non-uniformity, different consolidation degrees and thus joint quality along the overlap length will be achieved. Similar behavior was investigated in [17]. In Figure 4, the temperature distribution in the specimen during welding is shown. It can be noticed that the heat builds up in the joining zone and nearby areas of the specimen due to the local heat flow input and the insulator. Only with a certain delay, as it can be seen in Figure 3 as well, the heat is transported along the specimen's length direction. Although the heating element generates heat uniformly at the overlap area, the temperature distribution is not uniform. Instead, due to different heat transfer conditions between the edge and the middle of the joint, there is a maximum temperature at the middle and lower ones at the edges of the overlap. As a result of the non-uniformity, different consolidation degrees and thus joint quality along the overlap length will be achieved. Similar behavior was investigated in [17].
The time-temperature curve for two different power input levels for measuring point P1 is shown in Figure 5. While the melting temperature for an input power of 60 kW/m 2 is reached after 25 s, it takes only 18 s for an input power of 80 kW/m 2 . The cooling process follows a similar course. However, comparing these data with the ones shown in Figure 5, it can be assumed that due to the thermal resistance of the material, the heat cannot be transported fast enough, resulting in a smaller melted area, as the material needs to be prevented from degradation. Such data have to be taken into account for investigating an optimal processing window for the resistance welding of thermoplastic composites. When it comes to the influence of crystallization and different heating and cooling rates on mechanical properties, the estimation of transient temperature data for the process in an early design stage plays a significant role. os. Sci. 2021, 5, x FOR PEER REVIEW 6 of 10 The time-temperature curve for two different power input levels for measuring point P1 is shown in Figure 5. While the melting temperature for an input power of 60 kW/m 2 is reached after 25 s, it takes only 18 s for an input power of 80 kW/m 2 . The cooling process follows a similar course. However, comparing these data with the ones shown in Figure  5, it can be assumed that due to the thermal resistance of the material, the heat cannot be transported fast enough, resulting in a smaller melted area, as the material needs to be prevented from degradation. Such data have to be taken into account for investigating an optimal processing window for the resistance welding of thermoplastic composites. When it comes to the influence of crystallization and different heating and cooling rates on mechanical properties, the estimation of transient temperature data for the process in an early design stage plays a significant role.   The time-temperature curve for two different power input levels for measuring point P1 is shown in Figure 5. While the melting temperature for an input power of 60 kW/m 2 is reached after 25 s, it takes only 18 s for an input power of 80 kW/m 2 . The cooling process follows a similar course. However, comparing these data with the ones shown in Figure  5, it can be assumed that due to the thermal resistance of the material, the heat cannot be transported fast enough, resulting in a smaller melted area, as the material needs to be prevented from degradation. Such data have to be taken into account for investigating an optimal processing window for the resistance welding of thermoplastic composites. When it comes to the influence of crystallization and different heating and cooling rates on mechanical properties, the estimation of transient temperature data for the process in an early design stage plays a significant role.

Residual Stress Analysis
The welding process simulation is followed by a layerwise residual stress analysis. As a result of the high processing temperatures needed for the consolidation of thermoplastics, residual stresses arise and thus need to be investigated in a first step. Based on the proposed approach, the melted material sections are stress-free. The calculated residual stresses arise only due to the cooling process. Selected results are shown and discussed in the following. In Figure 6, the residual stress state of the 0 • layers, with respect to the global coordinate system, along the specimen's x-direction are presented. The data were collected in the mid-plane in the y-direction (y = 0.5w) for each layer at the center of each ply. The specimen's free end is located at x/l = 0. The overlap-zone reaches from x/l = 0.75 to x/l = 1.
ual stresses arise only due to the cooling process. Selected results are shown and discussed in the following. In Figure 6, the residual stress state of the 0° layers, with respect to the global coordinate system, along the specimen's x-direction are presented. The data were collected in the mid-plane in the y-direction (y = 0.5w) for each layer at the center of each ply. The specimen's free end is located at x/l = 0. The overlap-zone reaches from x/l = 0.75 to x/l = 1. Looking at the stresses in the fiber direction, there is a strong influence of the applied pressure in the area of the joining zone to notice. Especially in layer 2, there are high compressive stresses at the edge of the joining zone, while the stress state in between shows hyperbolic behavior with predominantly tensile residual stresses. This behavior is similar for all 0° layers. Considering the residual stresses transverse to the fiber direction, a similar stress trend for all layers can be noticed. From the free edge to the beginning of the joining zone, the residual stresses are close to zero. As shown in the welding simulation, the matrix in this area does not reach temperatures above the melting temperature. The more the area was heated, the higher the stresses. In layer 2, a maximum residual transverse stress of σ II = 221 MPa is reached. In Figures 7 and 8, the residual stresses for the 90° layers and +45° layer are depicted. It can be stated that all layers show comparable characteristics. The σ I stresses seem to be highly influenced by the applied pressure, which is needed for the consolidation of the adherends in the joining zone during cooling. Especially the peaks at the edge of the joining zone and the effect of the pressure put on them need to be considered in detail for future design steps. Looking at the stresses in the fiber direction, there is a strong influence of the applied pressure in the area of the joining zone to notice. Especially in layer 2, there are high compressive stresses at the edge of the joining zone, while the stress state in between shows hyperbolic behavior with predominantly tensile residual stresses. This behavior is similar for all 0 • layers. Considering the residual stresses transverse to the fiber direction, a similar stress trend for all layers can be noticed. From the free edge to the beginning of the joining zone, the residual stresses are close to zero. As shown in the welding simulation, the matrix in this area does not reach temperatures above the melting temperature. The more the area was heated, the higher the stresses. In layer 2, a maximum residual transverse stress of σ II = 221 MPa is reached. In Figures 7 and 8, the residual stresses for the 90 • layers and +45 • layer are depicted. It can be stated that all layers show comparable characteristics. The σ I stresses seem to be highly influenced by the applied pressure, which is needed for the consolidation of the adherends in the joining zone during cooling. Especially the peaks at the edge of the joining zone and the effect of the pressure put on them need to be considered in detail for future design steps.
All stresses shown arise due to the thermal history of the component effected by the welding process without external mechanical loads. While the stresses in the fiber direction arise in moderate magnitudes concerning the ultimate strength (R I,ult = 2070 MPa) and structural integrity of such structures, the transverse residual stresses are assumed to be critical. The maximum value is σ II,max = 221 MPa, observed in layers 1 and 2, which are the layers right next to the heating element. This is a multiple of the ultimate transverse strength of the used unidirectional CFRP layer (R II,ult = 90 MPa) and herewith leads to matrix cracks in the loaded area. Based on the determined stresses, matrix cracks will appear in every single layer, especially in the laminate area between the insulator and heating element.  All stresses shown arise due to the thermal history of the component effected by the welding process without external mechanical loads. While the stresses in the fiber direction arise in moderate magnitudes concerning the ultimate strength (R I ,ult = 2070 MPa) and structural integrity of such structures, the transverse residual stresses are assumed to be critical. The maximum value is σ II ,max = 221 MPa, observed in layers 1 and 2, which are the layers right next to the heating element. This is a multiple of the ultimate transverse  All stresses shown arise due to the thermal history of the component effected by the welding process without external mechanical loads. While the stresses in the fiber direction arise in moderate magnitudes concerning the ultimate strength (R I ,ult = 2070 MPa) and structural integrity of such structures, the transverse residual stresses are assumed to be critical. The maximum value is σ II ,max = 221 MPa, observed in layers 1 and 2, which are the layers right next to the heating element. This is a multiple of the ultimate transverse

Discussion and Conclusions
The presented modeling approach of the resistance welding process represents an effective procedure to estimate the transient temperature distribution in welded thermoplastic CFRP components for the application in an early design stage. Moreover, the approach allows being adopted to almost any adherend geometry. However, the simplification of the process in the FE model will inevitably lead to certain deviations from the actual conditions. For instance, the thermal boundary conditions will need to be validated by experimental methods. Furthermore, the power input simplified as surface heat flux is modeled as a thermal load applied immediately to the specimen, while in the case of an experimental investigation, a certain delay regarding the heat flow development by the resistance element is to be expected. For a more profound investigation, a far more resource-intensive Multiphysics Simulation is necessary, which then should model the actual resistance element and the joule heating effect. Based on the presented procedure and the results, a numerical parameterized optimization of the process parameters with regard to welding time, welding temperature, input power, and the melted area should be carried out.
A modeling approach for the investigation of residual stresses due to cooling after welding was proposed, and first results were shown. While the stresses in fiber direction predominantly arise in moderate magnitudes as tensile as well as compression stresses, the stresses in the transverse direction are purely high-tensile stresses. In particular, the stresses near the joining zone, which are a multiple of typical transverse strength, presumably lead to microcracks of the matrix. This pre-damage of the laminate might have had a significant effect on the strength and thus needs to be further investigated in future research. Especially, the fatigue behavior of welded structures possibly is highly affected by the welding-induced microcracks. For further investigation of residual stresses based on numerical methods, the effect of the applied material model needs to be clarified in detail. The interaction of plasticity and residual stresses constitutes a key aspect for thermomechanical simulations and thus has to be considered for ongoing activities by all means.
Finally, future work on this topic should focus on the effects of residual stresses on the strength and stability of welded components and thus investigate the importance of considering residual stresses during the design process of thermoplastic composite structures. Furthermore, the experimental validation of the numerical results needs to be carried out. In particular, the experimental investigation and determination of residual stresses bring along many challenges. This leads to a very small amount of published and accessible data that can be used for validation of FE models calculating residual stresses. Consequently, the design of concepts for the experimental investigation of residual stresses along with the validation of numerical methods represents another significant issue for future research. Funding: The authors thankfully acknowledge the financial and organizational support of the project JoinTHIS by the federal state of Lower Saxony and the European Regional Development Fund (ERDF). Furthermore, we acknowledge support by the German Research Foundation and the Open Access