A Reliable Fracture Angle Determination Algorithm for Extended Puck’s 3D Inter-Fiber Failure Criterion for Unidirectional Composites

Determination of the fracture angle and maximum exposure value of extended Puck’s 3D inter-fiber failure (IFF) criterion is of great importance for predicting the failure mechanism of unidirectional fiber-reinforced composites. In this paper, a reliable semi-analytical algorithm (RSAA) is presented for searching fracture angle and corresponding exposure value for the extended Puck’s failure criterion. One hundred million cases are tested for verifying the accuracy of the present and other algorithms on Python using the strength-value-stress-state combinations more universal than those in previous literatures. The reliability of previous algorithms is discussed and counterexamples are provided for illustration. The statistical results show RSAA is adequate for implementation in extended Puck’s criterion and much more reliable than previous algorithms. RSAA can correctly predict the results with a probability of over 99.999%.


Introduction
Plenty of failure criteria for fiber-reinforced unidirectional (UD) composites have been proposed during the past decades. Among all the failure criteria, the 3D inter-fiber failure (IFF) criterion developed by Puck and Schürmann [1] was ranked high in the "World Wide Failure Exercise I" (WWFE-I) [2] and "World Wide Failure Exercise II" (WWFE-II) [3] due to its good predictions compared with the provided experimental data. Recently, Gu and Chen [4] extended the Puck's failure criterion, which was expected to be applicable for different types of UD composites. Puck's failure criterion is developed based on a physical hypothesis that the fracture plane is determined by the shear stresses and normal stress acting on this plane [5]. However, the shear and normal stresses vary on different potential fracture planes and the plane with the highest exposure value is the actual fracture plane. Hence the criterion can not only predict the onset of failure but also can determine the orientation of fracture plane.
For a UD composites with an arbitrary 3D stress state and uncertain material properties, it is more important to know its potential fracture angle before predicting the onset of failure. Due to this reason, Puck [1] and VDI [6] proposed a stepwise search algorithm (SSA) to locate the fracture angle which relied on the degree-to-degree scanning of the potential fracture plane. Within an accuracy of 1 • , it will lead to 180 points to be calculated. Higher accuracy leads to more computational effort. Apparently, an efficient and reliable algorithm is required to implement the Puck's criterion.
Later, Wiegand et al. [7] proposed an accurate and numerically efficient fracture angle search algorithm called Extended Golden Section Search (EGSS) for Puck's IFF criterion. EGSS first uses the Golden Section Search (GSS) [8] algorithm to quickly bracket the maximum and then apply a curve interpolation technique called Inverse Parabolic Interpolation The Puck's IFF criterion is based on the Mohr-Coulomb theory and assumes transverse isotropy of the UD composites. Failure will occur on the fracture plane where only shear stresses, τ nt and τ nl , and normal stress σ n exist, which is illustrated in Figure 1. These three stresses can be calculated from σ 2 , σ 3 , τ 12 , τ 13 , τ 23 using Equations (1)-(3): σ n (θ) =σ 2 cos 2 θ + σ 3 sin 2 θ + 2τ 23 sin θ cos θ (1) τ nt (θ) = (σ 3 − σ 2 ) sin θ cos θ + τ 23 (cos 2 θ − sin 2 θ) (2) τ nl (θ) = τ 12 cos θ + τ 13 cos sin Puck's IFF criterion for unidirectional composites distinguishes the failure mechanisms between σ n ≥ 0 and σ n < 0. For each condition, Puck introduces the stress exposure value f IFF , given below, to measure the fracture risk and the fracture occurs when f IFF reaches one. For a given strength-value-stress-state combination, the angle θ is a variable. The problem is thus converted to finding the fracture angle which makes exposure value f IFF maximum.
Here, Y c means the transverse compressive strength and S 21 denotes the longitudinal shear strength of UD composite.

Determination of the Parameters
Puck et al. have already provided inclination parameters [12] for typical GFRP/epoxy (glass fiber reinforced plastic) and CFRP/epoxy (carbon fiber reinforced plastic) UD composites. These materials are regarded as intrinsically brittle materials. However, Gu and Chen [4] recently extended the Puck's failure criterion and found that the parameters are not adequate for materials with low transverse compression/tension ratios. UD compos- can be calculated using Equations (7) and (8). Inclination parameters p ⊥ should be obtained from experiments. In this paper, recommended values are utilized for the two parameters.
According to Puck et al., p is suggested due to that these two parameters should be approximately of the same magnitude [5]. Up to now, there are two unknown parameters left, R A+ ⊥ , p ⊥⊥ , which can be determined using the following equations: Semi-brittle materials: Brittle materials: Intrinsically brittle materials: Y T means the transverse tensile strength of UD composite. The shear strength S 23 is difficult to obtain experimentally because it is hard to induce a state of uniform shear and eliminate the geometry effect of specimens [13,14]. Thus, a simple formula provided by Christenson [15] is used to determine S 23 as follows:

Existing Algorithms for Searching the Fracture Angle
The search of fracture angle represents an optimization problem. What we should do is to find the angle which makes the stress exposure in Equation (4) greatest in given strength-value-stress-state case. However, purely analytical solutions are not available due to the complexity of the formulas, and therefore, semi-analytical and numerical approaches are applied. 1.
The curves are smooth; 2.
The curves have at most three local maxima; 3.
The distance between two adjacent local maxima is larger than 25 • .
Based on the above three characteristics, SRGSS is developed, the main idea of which is to isolate all the local maxima and then localize the maximum in the maximum-containing range using GSS algorithm. At first, the whole range (−90 •~9 0 • ) is evenly divided into 18 parts. Owing to that the minimum distance is larger than 25 • , two adjacent maxima cannot be distributed at the neighboring parts. Thus, the maximum-containing range can be found when a supporting point has two smaller neighbors. Then the global maximum can be determined by comparing the previous two or three local maxima.

Semi-Analytical Algorithm Proposed by Thomson et al.
Thomson et al. proposed a semi-analytical solution to the puck's IFF criteria by using an approximation of the potential fracture plane, and in this way, the computational effort of the search is greatly reduced.
First, the expressions of three tractions in Equations (1)-(3) are simplified into a single cosine function using basic trigonometric relations. In this way, the location, namely potential angle listed in Table 1, of its maximum for each traction can be easily determined. The Puck's IFF criterion is a combination of these tractions and the authors believe these angles can be used to approximate the global maximum due to that the fracture angle will tend to the maximum location associated with one of the three tractions when the traction is dominant with respect to the other two tractions. In practice, the simplest approach described in [10] is to consider the potential angles as initial candidates for the resulting global maximum and take the plane with the highest exposure as starting point for a numerical search of EGSS proposed by Wiegand et al. [7]. According to Thomson, this search algorithm can significantly reduce the computational effort and the result can converge on the global maximum instead of the local ones.

IAAGSS Algorithm Proposed by Wang and Zhao
The search algorithm proposed by Wang and Zhao is an extension of the algorithm proposed by Thomson et al. [10]. They take six angles, listed in Table 1, as initial candidates based on the hypothesis that when one of |σ n | max , |τ nt | max , |τ nl | max takes the dominant position with respect to the others, the resulting fracture plane will tend to the angles mentioned above. The difference between the two algorithms is that Wang takes both maximum and minimum values of the three tractions into consideration. The next step is to search the global maximum based on the six angles. In Wang's previous literature [16], an Analytical Approximation Golden Section Search (AAGSS) algorithm is proposed. AAGSS takes the highest two exposures (from the six initial candidates) as starting points for GSS algorithm. While IAAGSS improves the AAGSS algorithm and it takes the two planes neighboring to the plane with the highest exposure as starting points for GSS algorithm. The authors also conclude that IAAGSS can converge on the exact global maximum. In this paper, only IAAGSS is utilized for comparison with the present algorithm.

Proposed Changes for Searching the Fracture Angle
As mentioned in the introduction, the main purpose of this paper is to improve reliability of the algorithm for determining the fracture angle and exposure value of extended Puck's 3D IFF criterion. The algorithm proposed in this paper is based on that proposed by Thomson et al. [10]. The basic idea of this algorithm is using a semi-analytical approximation to localize all the local maxima near the global maximum and then using the GSS to determine the exact global maximum.
As discussed in Section 3.2, the three tractions in Equations (1)-(3) have to be simplified into single cosines using the trigonometric relations listed below: In this way, the expressions of the three tractions are given below: It should be noted that the expressions of τ nt , τ nl are not the same with those provided in Thomson et al. [10] and Wang et al. [11,16]. Most probably, the trigonometric formulas are used by the authors and there are misprints in their papers. After obtaining the simplified expressions of tractions, the locations of maximum and minimum can easily be defined, listed as θ 1 ∼ θ 5 in Table 2. Some explanations should be stated here that the frequency of τ nl is 180 • , which is half of that of σ n and τ nt . While the fracture plane is searched within the interval (−90 • , 90 • ). Hence, only one potential angle, θ 5 , for τ nl is obtained. For σ n and τ nt , both the maximum and minimum locations are calculated. In this way, the exposure curves are split in order to prevent the simultaneous occurrence of local maximum and minimum values in the same interval, which may affect the following GSS algorithm for search local maximum. It should be noticed that all these potential angles should be in the range of (−90 • , 90 • ), so the angles are calculated using different formulas according to different situations.
However, in actual cases, the final results are related to the interaction of the three tractions and the global maximum lies in between the angles for almost all the state cases. Hence, several initial candidates, θ 6 ∼ θ 10 , should be added in consideration of the stress weights. An operator, {}, is defined here to calculate the new added angles. An illustrative example for {θ 1 , θ 3 } is shown in Equation (18). After the operator calculation, the newly added angle should be transferred to range (−90 • , 90 • ) according to the frequency of 180 • if it is out of range. Here: In order to test the reliability of the fracture angle search algorithm, Python scripts are used to implement all the algorithms due to its outstanding performance in numerical calculation. The process of RSAA for extended Puck's IFF criterion is shown in Figure 2. Inputs should be generated randomly in stage-1 for the material properties and stress states, ranges of which are shown in Table 3. The types of composites include the semibrittle materials with Y C /Y T ∈ (1, 2.5), brittle materials with Y C /Y T ∈ (2.5, 3.45) and intrinsically brittle materials with Y C /Y T > 3.45 according to Gu and Chen [4]. Each of stress components are set to be in the range of [−100, 100] MPa which can represent all the possible stress state combinations. All the strength and stress parameters are randomly selected from the given ranges. It should be mentioned that the value of stress component may be greater than the corresponding strength and thus may lead to the stress exposure f IFF greater than one. However, f IFF and the stress component combination are linearly positively correlated, thus the location with the highest exposure, namely the fracture plane, remains the same. In step 3, the related parameters required in extended Puck's IFF criterion are determined based on the material types using equations in Section 2.2. In stage 2, RSAA is implemented to find the angle that makes the exposure value f IFF highest in the range between −90 and 90 degrees. In the first two steps, the exposure values of the ten initial supporting points, θ 1 ∼ θ 10 , are calculated, and then select the highest exposure, the value of which is assumed to be f 0 . In step 3, to localize the potential angle ranges that contain the global maximum, the supporting points whose exposures greater than 0.96 f 0 are recorded. In step 4, each supporting point, recorded in step 3, together with the two neighboring points form two search ranges. Then, all local maxima can be found using GSS algorithm. At last, global maximum and corresponding fracture angle are determined by comparing the aforementioned local maxima. An illustration example for present algorithm is shown in Figure 3. Exposure value greater than 0.96 f 0 is depicted in the gray area. Ranges selected for GSS algorithm are depicted in light blue area. Monte Carlo simulations were implemented in Python to obtain the characteristics due to the complexity of the Puck's criterion and stress state. To verify the accuracy of the present algorithm, outcomes obtained from SRGSS, SAA, and IAAGSS are compared with the outcomes obtained from Puck's original algorithm (SSA) with an accuracy of 0.01 • regardless of the search efficiency under the same strength-stress combinations. A total of 100,000,000 combinations were tested. Monte Carlo simulations were implemented in Python to obtain the characteristics due to the complexity of the Puck's criterion and stress state. To verify the accuracy of the present algorithm, outcomes obtained from SRGSS, SAA, and IAAGSS are compared with the outcomes obtained from Puck's original algorithm (SSA) with an accuracy of 0.01° regardless of the search efficiency under the same strength-stress combinations. A total of 100,000,000 combinations were tested.

SRGSS Algorithm
For SRGSS algorithm, its main prerequisites are those that the stress exposure curves have at most three local maxima and the minimum distance between two adjacent maxima is no less than 25 • . Besides, the method of localizing the local maxima is to compare the supporting points and obtain the points with two smaller neighbors. However, after thirty million random tests, the results show that the curves may have up to four local maxima and the minimum distance can be lower than 25 • (minimum distance is about 3 • based on our results). Details of the numbers of local maxima are listed in Table 4.
Hence, no local maximum exists in range (C, E) according to SRGSS. Therefore, the search range is (A, C) and it results in the local maximum using GSS algorithm. Actually, exposure values obtained from SSA and SRGSS algorithms are 5.4027 and 5.3941 respectively. Besides, there is another aspect which may also not be able to determine the global maximum. Since the minimum distance between two adjacent supporting points could be smaller than the conclusion from Schirmaier et al. [9], two local maxima may exist in one search range and the same problem may arise using GSS as is discussed by Schirmaier et al. range (C, E) according to SRGSS. Therefore, the search range is (A, C) and it results in the local maximum using GSS algorithm. Actually, exposure values obtained from SSA and SRGSS algorithms are 5.4027 and 5.3941 respectively. Besides, there is another aspect which may also not be able to determine the global maximum. Since the minimum distance between two adjacent supporting points could be smaller than the conclusion from Schirmaier et al. [9], two local maxima may exist in one search range and the same problem may arise using GSS as is discussed by Schirmaier et al.

Semi-Analytical Algorithm Proposed by Thomson et al.
The SAA may locate the incorrect global maximum, even the stress exposure curve has only one local maximum and an illustration example is shown in Figure 5. It should be stated that in the definition the θ 1 may locate the minimum of σ n in certain cases and we have used the form of θ 1 proposed in this paper, hence it can exactly locate the maximum of σ n in all stress states. Other initial candidates such as θ 4 are also modified based on our results. Another point that should also be mentioned here is that according to Thomson et al., the four initial candidates are calculated and take the plane with the highest exposure as the starting point for a numerical search of EGSS described in [7]. However, an ending point is also required to implement the EGSS algorithm. Hence, we take the planes with the highest two exposures as starting and ending intervals for EGSS, the method of which is also conducted in literature [17], which has used the algorithm proposed by Thomson et al. [10].
As is observed from the figure, the real global maximum is not in the range (A, B), the highest two exposures. The EGSS algorithm is then implemented and after several times GSS, the three points are determined for IPI, which results in the incorrect global maximum. θ are also modified based on our results. Another point that should also be mentioned here is that according to Thomson et al., the four initial candidates are calculated and take the plane with the highest exposure as the starting point for a numerical search of EGSS described in [7]. However, an ending point is also required to implement the EGSS algorithm. Hence, we take the planes with the highest two exposures as starting and ending intervals for EGSS, the method of which is also conducted in literature [17], which has used the algorithm proposed by Thomson et al. [10]. As is observed from the figure, the real global maximum is not in the range (A, B), the highest two exposures. The EGSS algorithm is then implemented and after several times GSS, the three points are determined for IPI, which results in the incorrect global maximum.

IAAGSS algorithm
The IAAGSS algorithm is an extension of SAA and used the initial candidates related to the three tractions. Similarly, the initial angles may not be capable to locate the maximum or minimum of the tractions or may be outside the range of (−90°, 90°). So, the modifications have been made so as to achieve what the authors expect. An illustration example is depicted in Figure 6. The stress exposure values at point A and B are 3.3958 and 3.3756 respectively. According to IAAGSS, candidate at point A is with the highest exposure and the range used for GSS is depicted in blue region in Figure 6. Obviously, a local maximum is obtained rather than a global one.

IAAGSS Algorithm
The IAAGSS algorithm is an extension of SAA and used the initial candidates related to the three tractions. Similarly, the initial angles may not be capable to locate the maximum or minimum of the tractions or may be outside the range of (−90 • , 90 • ). So, the modifications have been made so as to achieve what the authors expect. An illustration example is depicted in Figure 6. The stress exposure values at point A and B are 3.3958 and 3.3756 respectively. According to IAAGSS, candidate at point A is with the highest exposure and the range used for GSS is depicted in blue region in Figure 6. Obviously, a local maximum is obtained rather than a global one.

IAAGSS algorithm
The IAAGSS algorithm is an extension of SAA and used the initial candidates related to the three tractions. Similarly, the initial angles may not be capable to locate the maximum or minimum of the tractions or may be outside the range of (−90°, 90°). So, the modifications have been made so as to achieve what the authors expect. An illustration example is depicted in Figure 6. The stress exposure values at point A and B are 3.3958 and 3.3756 respectively. According to IAAGSS, candidate at point A is with the highest exposure and the range used for GSS is depicted in blue region in Figure 6. Obviously, a local maximum is obtained rather than a global one.

Comparison between the Present and Existing Algorithms
In this part, the final results, including the fracture angle and stress exposure value, are compared using present and the existing algorithms. The referenced results, namely the real fracture angle and exposure value, are determined using SSA. Thus, the distance (θ d ) between real fracture angle (θ r ) and the angles obtained from different algorithms (θ a ) and the relative error (r f ) between the real exposure value ( f r ) and exposures obtained using different algorithms ( f a ) are defined as: The angle distance and the relative error of exposure value obtained from the SSA and other algorithms were compared. Considering that the results from the SSA may not be the exact global maximum due to the incremental step of 0.01 • , angle distances larger than 1 • and relative errors of exposure value larger than 0.01% are shown in Tables 5 and 6. The angle distance and relative error are classified into four separate groups. It should be highlighted that the RSAA has obvious advantages over other algorithms in predicting both the fracture angle and related exposure value. The incorrect rates of the RSAA for angle distance and exposure relative error are only 5.7% and 1.74%, respectively, of those of SRGSS. Both the RSAA and IAAGSS are extended from the SAA, but the RSAA can correctly predict the results with a probability of over 99.999%.  It can be seen from Table 5 that there is still a small probability that the RSAA will miss the accurate fracture angle with angle distances greater than 1 • . Nonetheless, the relative error of stress exposure value shows no case greater than 1% for the RSAA (actually, the maximum is 0.27%). This phenomenon is attributed to the fact that the stress exposure curve may have two or more very close local maxima in certain cases even though the corresponding angles may vary obviously (e.g., in the uniaxial compression stress state, the exposure curve has two identical local maxima, but the corresponding fracture angles are 53 • and −53 • ).
The efficiency of all the algorithms including the Puck's original algorithm (SSA) under same strength-value-stress-state combinations with accuracy of 0.1 • is compared. The cost time for ten thousand random tests is listed in Table 7. It can be seen the cost time of present algorithm is only about 11% of cost time of SSA. However, the calculation effort is larger than those for SRGSS, SAA, and IAAGSS, but present algorithm has the essential advantage of giving reliable outcomes. In general, RSAA is adequate for predicting the fracture angles and exposures for the extended Puck's 3D IFF criterion and it is a balance of reliability and efficiency.

Application of Present Algorithm in FE Analysis
The present algorithm is implemented in the FE analysis using the material subroutine UMAT in ABAQUS standard. Considering the predicted strengths of multidirectional laminates depend much on the degradation model, UD composites are utilized in this part. A cubic FE model is established to represent the representative volume model of UD composites and x-dir denotes the longitudinal direction, which is shown in Figure 7d. Different types of UD composite materials including semi-brittle [18], brittle [19], and intrinsically brittle [20] materials are considered. Biaxial strength properties in σ y − τ xy stress space are obtained by applying periodical boundary conditions to FE model. Details for RVE model, periodical boundary conditions, and loadings can refer to our previous literature [21]. Figure 7a-c indicates that the FE results using present algorithm correlate well with the theoretical results using SSA algorithm with accuracy of 0.1 • .
It should be mentioned that the FE results using present algorithm and SRGSS algorithm are the same for the above simulation. The reason of this phenomenon can be attributed to that the stress states of the FE model are very simple, but the stress states can be significantly complicated and differences may appear between different algorithms for laminate and fabric composite material. If the damage location and the related criteria value cannot be accurately determined, it will greatly affect the prediction of strength property and the subsequent damage evolution process. However, as is discussed in Section 5.2, these stress states exist and the present algorithm shows superiority than other algorithms. It can be concluded that using the present algorithm can obtain the accurate results with higher probability compared with others.

Conclusions
In the present work, the extended Puck's 3D IFF criterion is utilized for verifying the reliability of the present algorithms for determining the fracture angle and maximum stress exposure value. The algorithm in this paper is based on the SAA proposed by Thomson et al. and the basic idea is to use semi-analytical approximation to localize the local maxima near global maximum and then obtain the global maximum using GSS. Thus, ten initial candidates are required for the present algorithm.

Conclusions
In the present work, the extended Puck's 3D IFF criterion is utilized for verifying the reliability of the present algorithms for determining the fracture angle and maximum stress exposure value. The algorithm in this paper is based on the SAA proposed by Thomson et al. and the basic idea is to use semi-analytical approximation to localize the local maxima near global maximum and then obtain the global maximum using GSS. Thus, ten initial candidates are required for the present algorithm.
The reliability and shortcomings of the existing algorithms are discussed and counterexamples of their algorithms are plotted for illustration. Present algorithm avoids these shortcomings. One hundred million strength-stress combination cases were tested to verify the reliability of the existing algorithms. Statistical results indicate that the RSAA is much more reliable than SRGSS, the SAA, and IAAGSS, which make it a practical alternative for the application of Puck's 3D failure criterion in composite analysis. The computational effort is larger than the other algorithms, however, it is only about one-tenth the time cost of the SSA. Present algorithm is implemented in ABAQUS and biaxial strength properties in σ y − τ xy stress space are calculated, which are in good agreement with theoretical results using SSA. In general, present algorithm is a balance of reliability and efficiency and adequate for evaluating the Puck's 3D IFF criterion.