Combined experimental–numerical mode I fracture characterization of the pultruded composite bars

In this paper, pultruded GFRP bars are investigated to determine their fracture properties. The double cantilever beam test (DCB) is used to assess fracture behavior under mode I loading conditions. However, due to the presence of the R-curve effect (variable fracture energy dependent on the length of the crack), it is necessary to introduce a nonstandard approach to determine fracture properties. The mixed experimental–numerical approach is proposed to deal with this issue. Numerical simulations were carried out in Simulia Abaqus, and with Python scripting it was possible to generate models and obtain R-curve for the material. The numerical model built based on the experimental results has very good agreement with it (force–displacement and delamination length–time characteristics) which allows the use of the mentioned model in the analysis of more complex structures. Acoustic emission analysis was introduced as an auxiliary technique. The delamination obtained from both the numerical model and the experiment complies with the registered acoustic emission events. The proposed method can be used in preparing a material model for other composite materials, which display the presence of the R-curve effect.


Introduction
The advantages of using FRP composites are well-known. These materials have many effective applications in aircraft, automobiles, marine structures, and sports equipment due to their good corrosion resistance, electromagnetic transparency, and high strength-weight radio [1]. Different composite materials manufacturing techniques are available, for example, filament winding [2][3][4] or pultrusion [5]. The FRP bars are usually made of glass (GFRP), basalt (BFRP), carbon (CFRP), or aramid (AFRP) fibers. Glass fibers are the most frequently used in polymeric matrix composites, and GFRP bars were the first to be used and continue to be the most frequently used in various applications due to their high strength and lower price compared to the other fibers [1].
To expand the use of FRP composites, a detailed analysis of their mechanical properties is essential. For structural applications, good toughness or resistance to delamination is a major concern [6]. Resistance to delamination can be determined by different load modes, such as tensile opening (mode I), shear (mode II), or twist (mode II) [7]. The double cantilever beam test (DCB) is the method adopted by ASTM D5528 [8] to determine the delamination resistance of composites under loading conditions of mode I. The test defines interlaminar fracture toughness, G IC [9] by tensile opening, and was developed to measure the delamination resistance of unidirectional carbon and glass fiber tape laminates.
In the literature, many articles are available to assess the resistance to delamination of GFRP [10,11] BFRP [12], and CFRP [13][14][15] laminates using the DCB test. However, for FRP bars, the fracture criteria were already addressed only in a couple of papers for the GFRP and BFRP bars [7,9,16], and information on this topic is still very limited. Additionally, there are no standards available to evaluate the delamination resistance of these bars.
The research investigates the energy release rate for mode I using a standardized experimental procedure and numerical analysis. The fracture energy is determined on the basis of the experiment using data reduction schemas using, for example, the concept of compliance. This approach is simple as it requires only conducting experimental tests and data post-processing. In some cases, this approach is enough [17]. However, in many other cases, it is not sufficient. In that case, another solution must be used-for example, the authors successfully used a mixed experimental-numerical approach in the case of specific composites-fiber metal laminates with material and stiffness asymmetry [18,19]. Another issue is the so-called R-curve effect. The R-curve describes the relationship between fracture energy (crack growth resistance) and crack length (Fig. 1).
The R-curve effect is often used to explain the behavior of ceramic material [20][21][22]. However, this effect is also observable in composite materials. For this material group, it can be an effect of fiber bridging and is connected to the so-called large-scale bridging [23]. In the case of some composites, the registered effect is not very noticeable (i.e., [24]); however, in other cases, it can be more significant (i.e., [25]). Catalanotti in [26] points out that the R-curve effect can depend on the kind of resin or the presence of additional effects such as kink bands, etc. In the presence of the R-curve effect, it is also necessary to extend the basic approach of determining fracture parameters. The natural choice is to introduce numerical modeling using the finite element method.
The main approach to modeling the delamination behavior in composite material using the finite element method in recent years is the cohesive zone model [18,19,[23][24][25][26][27]30]. This approach is derived from the work of Dugdale [31] and Barenblatt [32]. The idea of the model is that, after the onset of the fracture, the loss of the ability to transfer stresses is not immediate. Instead, this is a gradual process. The character of these gradual changes can be different for different materials, so the shape of the so-called traction-separation law which describes that behavior can vary, from a simple bilinear law to another like a trapezoid or exponential. Traction represents stresses ( I ) and separation denotes the distance between separated surfaces ( I ) , where subscript I indicates mode I loading. The area under the curve describes the fracture energy (in case of mode I loading-G IC ). In Fig. 2, exemplary traction-separation law is presented.
The cohesive zone model can be applied to finite element analysis software in different ways. The natural choice is to use cohesive elements (cohesive parameters defined as material property). They can be used with finite thickness [19,33], or zero thickness [34]. Another possibility is to use a cohesive surface (cohesive parameter defined as contact property) [35,36]. Finally, one can use eXtended Finite Element Method (xFEM). In this approach, cohesive behavior is defined as the property of the material and the crack can propagate through elements. If the crack path is known a priori, the use of cohesive elements can be a more efficient solution [29]. However, if the solution is transferred to a more complex geometry, the xFEM solution can be better Fig. 1 The increasing R-curve-the fracture energy increases from the initial value with crack growth to reach the so-called plateau and the steady rate of delamination growth. Rising R-curve is an alternative observed in ceramics and some composites to the simplest case with constant fracture energy  [37]. xFEM is used successfully to investigate crack growth in metal and plastic materials to predict failure (see, e.g., [38]). Overall, for DCB-type specimens, all the above-mentioned methods will give very comparable solutions. In this paper, the xFEM approach is presented because of its ability to propagate crack through elements, which will allow one to use one element for the height of the pre-crack cut.
However, including the R-curve effect in the numerical simulation requires a nonstandard approach. There are limited articles dealing with this problem, but the main approach seems to be the use of the modified shape of the traction-separation law, namely trilinear [39][40][41][42]. The trilinear traction-separation law is the superposition of two bilinear laws defined for the initiation energy and the propagation energy. The idea of such a solution is presented in Fig. 3.
This solution shows a satisfactory agreement with the experimental data mentioned in the mentioned articles. However, other methods are also proposed; for example, Latifi et al. [23] used the interface thick level set method. A potential issue with the use of the trilinear traction-separation law is that it can be less reliable in the case of materials with a slow rising R-curve. In this paper, research conducted on such material will be presented, which includes another approach to include the R-curve effect in the numerical model. Moreover, for the main experiment and numerical investigation, an additional technique will also be used, which is acoustic emission analysis. Acoustic emission is a well-known non-destructive method that is also used in the investigation of fractures in composite material [17,19]. Also other supportive techniques as digital image correlation are currently popular (see e.g., [43]); however, they were not necessary in this research.
The paper presented focusses on proposing a new approach to obtain and define the effect of the R-curve in composite materials used as reinforcement for concrete. Introductory investigations without including the mentioned effect show a strong discrepancy with experimental data and data obtained from acoustic emission. The method presented in the paper is an important step in the credible description of the researched materials in the language of numerical models.

Materials and methods
This section presents the material, equipment, and procedure that are used and developed to provide a useful tool and methodology for investigating the fracture energy in pultruded bars.

Material preparation and experimental setting
The investigated samples were manufactured using the pultrusion method. This technology can be used to provide various unidirectional profiles with constant cross-section. For the purpose of this research, the emphasis was placed on circular bars reinforced with glass fibers in an epoxy matrix (Biresin CR141). Unidirectional circular samples with a 10 mm diameter were applied for investigating the double cantilever beam (DCB) test.
DCB specimen required a specially prepared geometry that allows investigation of the pure energy release rate for mode I, and for this reason, the initial notch was adopted. The initial notches were machined using a diamond wire saw. The provided notch is responsible for the proper direction of delamination. Additionally, to transfer the load, special loading blocks (U-shaped metallic blocks) were manufactured and adhered to the circular composite surface using cyanoacrylate glue. Finally, to observe the delamination growth, markers were placed every 5 mm. Samples with initial delamination length and loading blocks are presented in Fig. 4. The material data of used resin are presented in Table 1.
The double cantilever (DCB) test required a specially prepared geometry that allows investigation of the pure energy release rate for mode I, and for this reason, the initial notch was adopted. The provided notch is responsible for the proper direction of delamination. Additionally, to transfer the load, special loading blocks were manufactured that adhered to the composite surface using cyanoacrylate glue.
Composite materials exhibit a complex nature led by various failure mechanisms. These mechanisms propagate in a rather general fashion, unlike metallic materials, which propagate locally. To analyze the failure events for the particular loading conditions of the structure, acoustic emission is used. The method allows registering acoustic signals that appear in front of the crack tip and are caused by events due to the fiber pullout and breakage. The Vallen AMSYS-6 system was used with two piezoelectric sensors placed as shown in Fig. 4. VS150-M sensors with The DCB test was carried out on the MTS Bionix 793 system with a displacement control mode of 5 mm/min. The experiment followed the methodology proposed by ASTM D5528. It consists of two parts; at the beginning, every 5 mm of the length after an initial cut is marked with black lines (the first thick exhibit the crack tip). The first step of the experiment requires a pre-crack that must be executed. The length of the pre-crack should be in the range of 3-5 mm. The cause of pre-crack is to reduce the radius of the tip of the crack and to provide repeatable conditions for each sample. The experimental setup used in this investigation is shown in Fig. 5.
The experiment was conducted for both cases. First, for the bar only with one sensor and the second case with two sensors and a localization algorithm. This approach allowed confirmation of the influence of the bending moment caused by the sensor on the acoustic events. This analysis showed that the second AE sensor does not influence acoustic events in the investigated samples.

Numerical modeling
The investigated material is fiber-reinforced epoxy resin. The DCB test is meant to be a method to determine the energy of the fracture according to ASTM standards. However, the force-displacement curve obtained in such a test has a specific shape, which is a linear part until maximum load and then a hyperbolic drop of the force with the progress of displacement. This is true for materials with fracture energy independent of crack size (delamination). However, in some materials, the socalled R-curve effect is present, which can be described as an increase in crack growth resistance with growing crack. This effect makes it impossible to determine one fracture energy from the DCB test. At the same time, a proper material model Fig. 4 Basic dimensions of the specimen along with the position of the sensors. The force is applied through auxiliary blocks of metal which are glued to the specimen and will be attached to the mounting system. Auxiliary markings were made every 5 mm to allow for growth detection of crack length  5 Experimental setup to investigate the fracture energy in mode I. Special mounting is used to transfer the force from the universal testing machine to the specimen. On the right side of the sample, an acoustic emission sensor can be seen is very useful for modeling the behavior of composite rods in more complex structures. This is the reason why a mixed experimental-numerical approach was introduced to determine a suitable material model for future research. Based on the experiment, the crack growth resistance curve was introduced to the numerical model. The numerical simulations were carried out using Simulia Abaqus. The geometrical model was created in a three-dimensional space using solid elements type C3D8R according to the dimensions shown in Fig. 4. Due to the symmetry of the experiment, only half of the specimen was modeled to reduce computational cost. It is also possible to use only a quarter of the model, but half was chosen due to the stability of the crack growth modeling. The material model was defined using the orthotropic material definition with proper material orientation. The material data were based on experiments conducted during the research. Basic data about a material are presented in Table 2.
The key part of the material model is to define the behavior of the fracture. The cohesive zone model is commonly used to model fracture behavior in composites [18,19,28,29]. In the cohesive zone model, the fracture behavior is introduced via the traction-separation law. The shape of the traction-separation law can vary; however, commonly used is a bilinear shape. The cohesive zone model can be used in various applications, including cohesive elements, surfaces, and xFEM. In this research, cohesive zone model was used along the eXtended Finite Element Method. As underlined in the Introduction, the R-curve effect needs to be included in the model, as the shape of the experimentally determined force-displacement curves indicates that the fracture energy is not constant. In the presented research, the way of including the R-curve effect was inspired by the work highlighted in the Introduction (i.e., [39]); however, instead of using the trilinear cohesive law or similar methods to imitate fracture behavior, in this paper, the bilinear cohesive law was used with gradual parameter change while delamination is growing. For simplicity, the character of the curve was assumed to be hyperbolic based on the available data for similar material. The initial and propagation fracture energies were obtained by assuming that they can be calculated for the initial force loss area and the areas with high displacement values, respectively. Fracture values were determined by performing auxiliary simulations with constant fracture energy, and then the energy values can be associated with the crack length values (a). Examples of force-displacement curves with different energy values are presented in Fig. 6 (left).
Based on the values of the fracture energy used to generate a model that reflects a particular value of the length of the crack, it is possible to define the R-curve as presented in Fig. 6 (right). The initial energy determined as 0.64 N/mm is based on the simulation with this energy (red line in Fig. 7). Further propagation energy can be derived from simulations with higher energy (for example, blue line in Fig. 7.) Determining the R-curve is only the first step in the proposed method. After that, we must transfer this information to the numerical model. As mentioned above, it can be done by altering the traction-separation law used for the cohesive zone model. However, for the purpose of validation, it seems that a good approach is to use multiple sections with the material defined taking into account the location of the section and taking the corresponding fracture energy based on the determined R-curve. In the research presented, the width of the section was assumed to be 1 mm. Due to the large number of sections, it is not reasonable to prepare such a model by hand, so Python scripting was used to automatically generate models and assign the material model or boundary conditions. An example of a model generated using this technique is presented in Fig. 8. Each color reflects a different material section. In the figure, there are also visible boundary conditions-displacement applied via kinematic coupling and restrictions that are the effect of modeling only half of the specimen. For simplicity reasons, without losing accuracy, the specimen is cut in the  plane where the load is acting, and coupling is assigned to these surfaces. The size of the mesh can be an important factor that influences the quality of the mesh. On the other hand, a too fine mesh can lead to significant computational time. The standard technique is to use coarser mesh in regions that contribute little to the results [44]. In this case, the coarser mesh is used in the area where the pre-crack is already introduced. To find the balance between the mesh quality and time consumption, often mesh sensitivity studies are conducted during preliminary studies [45][46][47]. In the case of the presented method, there is also one more factor-size of the section. Both factors were checked in the sensitivity studies (Fig. 9).
The size of the section is found to be an important factor in the proposed approach. Too large sections size causes the final output to be not smooth enough to be a reliable output. In the case given in this study, the optimal section width is set to be 1 mm. Due to the small size of the section, the number of elements in the section is less important to the overall results. However, differences can still be spotted, so four elements (0.25 mm width) will be used in each section in the final models.
The presented method to obtain the proper material model is a good and fast method that can be an alternative to obtaining the same information through the cyclic loading and unloading of the specimen during the DCB test.

Results and discussion
In this section, the results obtained from the experiment, acoustic emission, and FEM are highlighted and elaborated. Several graphs are presented and a detailed discussion is provided.
The fiber bridging effect was noticeable during the experiment. This mechanism is an important phenomenon that limits the growth of the delamination. This phenomenon is illustrated in Fig. 10. The size of the processing zone has significant dimensions in comparison to the dimensions of the specimen. This confirms that the R-curve effect is observable in that material and is consistent with the data obtained from the machine (force-displacement results).
The numerical model was prepared as described in the previous section based on auxiliary simulations that effect the determination of the R-curve (Fig. 6). The results of this model are presented in Fig. 11 (brown line) along with the experimental results. The model presents very good agreement with the experimental results, so it can be recognized as valid for the future use in simulations. It is visible that models with a bilinear cohesive law with constant energy cannot reliably reflect the experiment (red and blue lines).
In Fig. 12, quantitative results are presented. On the left, stresses (maximum principal stresses, MPa) are shown-the highest tension values are present in the crack area, while the highest compression value is in the corresponding places in the external part of the sample. On the right, the current state of bonding is presented with the internal variable of Abaqus (STATUSXFEM).
The force-displacement characteristic presented in Fig. 11 is a good validation of the model. However, apart from this comparison, it is worth also checking the delamination growth as a function of time or displacement. In Fig. 13, observed delamination from the finite element method (yellow) is shown along with the delamination observed in the experiment (black). Both characteristics are presented with acoustic emission events in the background. There is a good  The colors indicate the sections that refer to different values of fracture energy. The displacement is applied to the reference point and transferred via coupling. Reaction values (F) can be read after simulation for these reference points agreement between the experimental and numerical results, but also with the boundary of acoustic events occurring in a given time. This confirms the correct choice of the R-curve. It is worth underlining that the highest amplitude events occur more frequently in the crack tip area. The presence of acoustic emission events at a greater distance from the delamination front proves that the large-bridging behavior occurs and complies with optical observations. Additional acoustic emission analyses were conducted to better understand the behavior of the material and analyze the possibility of acoustic emission application in the structural health monitoring (SHM) of complex structures constructed with the researched material (pultruded composite rods). The force versus time curve obtained from the universal testing machine has been correlated with the cumulative number of AE events (Fig. 14). It is proposed to highlight three different sections that analyze the force and AE events curve. The following sections were selected; first, region A reflects the increase of progressive damage events, whereas the force changes linear behavior and reaches the maximum value. In this region, many micro-failures occur rapidly, such as matrix micro-cracking. It happens until the delamination is initiated, after which the crack is going to be stable, as well as acoustic events that happen not only in front of the  Additionally, two characteristics derived from models with constant energy, the bilinear law, are shown. These approaches cannot reflect the real behavior of the material crack tip. When region A is exceeded, the crack is stable and propagates gradually; the constant increase of the acoustic events is noticeable in terms of the linear behavior for a cumulative number of events. In this area, the force also is gradually decreasing. The largest number of acoustic events is caused by fiber breakage in the area before the crack. After a certain value of the total length of the delamination is reached, progressive damage accumulation appears with a significant force drop. It might be caused first by a decrease in stiffness leading to a decrease in force. The change in the cumulative number of acoustic emission events may be used as an indicator for the loss of properties of the material used in complex structures. However, the effect is not so observable compared to what is observed in other materials; for example, fiber metal laminates [19]. Further investigation in the area of useful characteristics derived from acoustic emission needs to be performed.
Another acoustic emission parameter that points to important moments in the experimental tests is the socalled root mean squared of the noise. By investigating the noise level, this characteristic can indicate rupture changes in the behavior of the material. In Fig. 15, the agreement between the RMS of noise and the characteristic of force is visible. Each significant fall in the force is accompanied by an increase in the RMS signal. Rupture changes in RMS and force values were also observed during the experiment, as events in the area of the processing zone have such a character.
Steps required to determine fracture energy (in function of delamination/crack size) in materials with present R-curve effect are presented in the form of block schema in Fig. 16

Conclusion
In this paper, the results of fracture behavior under the mode I loading conditions of pultruded composite rods are presented. The method to describe this behavior was proposed by determining the fracture parameters that include the R-curve (Fig. 6. The block diagram of this method is shown in Fig. 16. As an auxiliary research method, acoustic emission analysis was used. The main reason was to provide reliable information regarding the damage and failure modes that occur in the material during the experiment. Additionally, the acoustic emission technique proved to be capable of tracking the delamination growth because it shows satisfactory agreement with both the optical method and numerical simulation. In summary, the following conclusions may be drawn:  The observed delamination front based on optical observations during the experimental DCB test (black line) is presented along with the location of the delamination front derived from finite element analysis (yellow line). There is a good agreement between them and also the presence of acoustic emission events. The wide band of high-amplitude acoustic emission events at any given time complies with the assumed large processing zone • The fracture energy in delamination (mode I fracture) in the investigated material depends on the size of the crack (R-curve effect). • The proposed algorithm (Fig. 16) can be used to determine the appropriate material model for composites in the presence of an R-curve effect using the DCB test and numerical modelling. To facilitate the FEA work, the automatically generated section using Python scripting may be used to assess the R-curve effect in the numerical model. • Numerical modeling based on the eXtended Finite Element Method (xFEM) and cohesive zone model (CZM) can be used to model the delamination behavior in a pultruded composite rod with satisfactory agreement with experimental data. • Acoustic emission (AE) measurements were found to be useful for composite materials due to the variety of  16 Block diagram of the proposed method, which can be used to determine fracture parameters of composite material that are affected by the effect of the R-curve damage mechanisms in these structures. The analysis of the intensity of measurements, duration, and frequency of the signal provides valuable information on the failure modes. AE provided information about the location of the acoustic event. It is a very practical tool, which enables researchers to specify the place of destruction. It is also possible to follow the development of the crack.
The main limitation of the presented methods is the necessity of using another approach to implement the obtained R-curve model to more complex geometries. This also indicates the further work direction, which can include other parts of the framework based, for example, on own written Abaqus subroutines. Data availability Data used in this study is part of an ongoing project and cannot be shared at the moment of the publication.

Declarations
Conflict of interest All authors declare that they have no conflicts of interest.

Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.