A Robust Fast Fracture Plane Orientation Angle Search Algorithm for Puck 3D Inter-Fibre Failure Criterion

Nanda Wirawan^{1, 2,}*, Ibrahim H. Abuzayed¹, Mahesa Akbar³ and Jose L. Curiel-Sosa¹

¹Department of Mechanical Engineering, University of Sheffield, Sheffield S1 3JD, United Kingdom
²Research Center for Aeronautics Technology, National Research and Innovation Agency (BRIN), Bogor 16350, Indonesia
³Faculty of Mechanical and Aerospace Engineering, Institut Teknologi Bandung, Bandung 40132, Indonesia
E-mail: nwirawan1@sheffield.ac.uk; nanda.wirawan@brin.go.id
*Corresponding Author

Received 16 October 2023; Accepted 07 May 2024

Abstract

In the present work, a novel fast fracture plane orientation angle (FPOA) search algorithm for the 3D Puck failure criterion is proposed. In the 3D Puck failure criterion, a linear search algorithm is employed to calculate the maximum inter-fibre failure (IFF) value by iterating and comparing the IFF value for each FPOA. This process itself requires a substantial amount of computational resources. The proposed fast FPOA search algorithm is implemented to substitute the linear search algorithm in order to reduce the computational time. A total of $1\times {10}^5$ randomised stress cases are used to analyse the accuracy of the algorithm. The result was then compared with the Puck Stepwise Seach Method (SSM) and other fast FPOA search algorithms. The results show that the proposed fast FPOA search algorithm has better accuracy compared to the other fast FPOA search algorithms and is almost 5 times faster compared to the SSM algorithm by Puck. In addition, a subroutine contains the Puck failure criterion and the proposed fast FPOA search algorithm is embedded into a Finite Element Analysis (FEA) software to simulate the open-hole test (OHT) experiment on the composite material.

Keywords: Composite material, 3D Puck failure criterion, fast fracture plane orientation angle search algorithm.

1 Introduction

In the past few decades, the usage of composite material has increased significantly, particularly in the automotive and aerospace industries, due to its high specific strength characteristic. The composite material is also highly customisable. It could be tailored by using different combinations of fibre/matrix and different ply stacks or orientations to conform to the loading condition requirement. As a consequence, failure modes on composite material are becoming more complex and cannot accurately be represented by simple failure criterion evaluation such as Tresca or Von Misses failure criterion.

Initially, Tsai and Wu [1] developed a failure criterion based on the failure strength of the material in each direction. In the Tsai-Wu failure criterion, the damage does not differentiate between fibre or matrix failure. Hashin [2] then improved the composite material failure evaluation by calculating fibre failure and matrix failure separately. Later, In the World-Wide Failure Exercise (WWFE) [3], a novel failure criterion developed by Puck [4–6], which is based on the Hashin failure criterion, was able to reproduce the failure under several stress conditions accurately compared to other failure criteria. Thus, the Puck failure criterion becomes one of the options for evaluating composite failure.

Puck uses fibre failure (FF) stress exposure ($f_{E,\mathit{FF}}$) and inter-fibre failure (IFF) stress exposure ($f_{E,\mathit{IFF}}$) to determine the damage that occurs across the fibre and parallel to the fibre, respectively. In 3D case stress, a potential IFF fracture plane is required to determine the $f_{E,\mathit{IFF}}$. The potential IFF fracture plane is obtained using the Stepwise Search Method (SSM) algorithm by comparing the highest value of $f_{E,\mathit{IFF}}$ from all of the action planes parallel to the fibre. Puck and VDI [7] propose that the action planes are inclined at 1${}^{\circ }$ intervals. As a result, 180 iterations are required to acquire the IFF fracture plane for each stress case, which leads to high computational costs.

Due to its high computational cost, several fast FPOA search algorithms were developed to reduce the computational cost of the SSM algorithm. Weigand [8] proposes the Extended Golden Section Search (EGSS) algorithm, which utilises the Golden Section Search (GSS) algorithm [9] and successive parabolic interpolation technique [10] to obtain the IFF fracture plane orientation angle (FPOA). However, for several stress conditions, the EGSS algorithm incorrectly identifies the global maximum value of $f_{E,\mathit{IFF}}$ location.

Schirmaier [11] improves the EGSS algorithm accuracy into the Selective Range Golden Section Search (SRGSS) algorithm using the selective range (SR) method, which divides the iteration range into several smaller blocks and then applies the GSS algorithm into blocks that contain the local maximum value of $f_{E,\mathit{IFF}}$. Rezasefat [12] later proposes the Simple Parabolic Interpolation Search (SPIS) algorithm. Instead of using the GSS algorithm, Rezasefat employs the successive parabolic interpolation technique in the SR method.

The fast FPOA search algorithm works by reducing the total amount of iteration required to obtain IFF FPOA. Since the data points from the iteration are less than the SSM algorithm, the accuracy of the fast FPOA search algorithm is also affected. In this paper, a robust, fast FPOA search algorithm is proposed. Several new methods are implemented into the algorithm to alleviate the accuracy problem without increasing the computational time significantly. The result is then compared with the previous fast FPOA search algorithm and then applied as a subroutine in finite element analysis software.

2 Puck Inter-Fibre Failure (IFF) Calculation

The Puck [6] IFF criterion expands the work by Hashin’s idea on matrix mode, which is based on Mohr failure theory [13]. According to Hashin, if the failure plane could be identified, the failure is only caused by the normal stress (${\sigma }_n$), shear stress transverse to the fibre direction (${\tau }_{nt}$) and shear stress parallel to the fibre direction $({\tau }_{n1})$ acting on the oblique plane. Puck uses the normal (${\sigma }_n$) and shear stresses (${\tau }_{nt},{\tau }_{n1}$) as inputs to calculate the failure exposure on IFF. The failure exposure on IFF is obtained by finding the maximum value of the failure exposure for every “action plane” that is rotated parallel to the $X_1$ axis on each inclination angle $\theta$.

Figure 1 The stress tensor on an element (left) and the transformation of the stress tensor into the potential fracture plane with inclination angle $\theta$ on the same element with the top half is omitted for clarity (right).

The normal and shear stresses on the action plane are derived from the stress tensor. The stress tensor is transformed according to the current inclination angle. The ${\sigma }_n$ and ${\tau }_{nt}$ are derived from the stress tensor components ${\sigma }_{22},{\sigma }_{33},{\sigma }_{23}$, while ${\tau }_{n1}$ is derived from the stress tensors ${\sigma }_{12},{\sigma }_{13}$, as shown in Equations (1)–(3). Figure 1 shows the normal (${\sigma }_n$) and shear stresses (${\tau }_{nt},{\tau }_{n1}$) that are being transformed on the action plane at the inclination angle $\theta$.

${\sigma }_n(\theta )$	$={\sigma }_{22}\cos {(\theta )}^2+{\sigma }_{33}\sin {(\theta )}^2+\mathrm{\hspace{0.25em}2}{\tau }_{23}\sin (\theta )\cos (\theta )$	(1)
${\tau }_{nt}(\theta )$	$=-{\sigma }_{22}\sin (\theta )\cos (\theta )+{\sigma }_{33}\sin (\theta )\cos (\theta )+{\tau }_{23}(\cos {(\theta )}^2-\sin {(\theta )}^2)$	(2)
${\tau }_{n1}(\theta )$	$={\tau }_{31}\sin (\theta )+{\tau }_{21}\cos (\theta )$	(3)

Puck states that the failure exposure value on the action plane at the current angle can be calculated using Equation (4) for tension load (${\sigma }_n\ge 0)$ and Equation (5) for compression load (${\sigma }_n<0$).

$f_{E,IFF+}(\theta )$	$=\sqrt{{\left[\left({\displaystyle \frac{1}{R_{\perp }^t}}-{\displaystyle \frac{p_{\perp \psi }^t}{R_{\perp \psi }^A}}\right){\sigma }_n\right]}^2+{\left({\displaystyle \frac{{\tau }_{nt}}{R_{\perp \perp }^A}}\right)}^2+{\left({\displaystyle \frac{{\tau }_{n1}}{R_{\perp \parallel }^A}}\right)}^2}$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{p_{\perp \psi }^t}{R_{\perp \psi }^A}}{\sigma }_n\mathit{\hspace{1em}}\mathit{for}{\sigma }_n\ge 0$	(4)
$f_{E,IFF-}(\theta )$	$=\sqrt{{\left({\displaystyle \frac{{\tau }_{nt}}{R_{\perp \perp }^A}}\right)}^2+{\left({\displaystyle \frac{{\tau }_{n1}}{R_{\perp \parallel }^A}}\right)}^2+{\left({\displaystyle \frac{p_{\perp \psi }^c}{R_{\perp \psi }^A}}{\sigma }_n\right)}^2}$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{p_{\perp \psi }^c}{R_{\perp \psi }^A}}{\sigma }_n\mathit{\hspace{1em}}\mathit{for}{\sigma }_n<0$	(5)

with

$$\frac{p_{\perp \psi }^t}{R_{\perp \psi }^A}=\frac{p_{\perp \perp }^t}{R_{\perp \perp }^A}co{s}^2\psi +\frac{p_{\perp \parallel }^t}{R_{\perp \parallel }^A}si{n}^2\psi$$	(6)
$$\frac{p_{\perp \psi }^c}{R_{\perp \psi }^A}=\frac{p_{\perp \perp }^c}{R_{\perp \perp }^A}co{s}^2\psi +\frac{p_{\perp \parallel }^c}{R_{\perp \parallel }^A}si{n}^2\psi$$	(7)
$$co{s}^2\psi =1-si{n}^2\psi =\frac{\tau _{nt}{}^2}{\tau _{nt}{}^2+\tau _{n1}{}^2}$$	(8)
$$R_{\perp \perp }^A=\frac{R_{\perp }^c}{2(1+p_{\perp \perp }^c)}$$	(9)

The $R_{\perp \perp }^A$, $R_{\perp \parallel }^A$, $R_{\perp \psi }^A$ are the fracture strength of the material. The $R_{\perp \perp }^A$ quantifies the fracture strength of the action plane due to shear stress. The $R_{\perp \perp }^A$ is calculated using Equation (9). The $R_{\perp \parallel }^A$ is the shear strength that is obtained from the in-plane shear stress test. The $R_{\perp }^t$ and $R_{\perp }^c$ are the tensile and compressive strength of the material transverse to the fibre direction. The $R_{\perp }^t$ and $R_{\perp }^c$ are obtained using an uniaxial test on the material.

Table 1 Inclination factor for different types of composite material [5]

Material	$p_{\perp \perp }^t$	$p_{\perp \perp }^c$	$p_{\perp \parallel }^t$	$p_{\perp \parallel }^c$
GFRP	0.30	0.25	0.20	0.25
CFRP	0.35	0.30	0.25	0.30

Figure 2 Master fracture body ${\sigma }_n-{\tau }_{n1}$ at $\psi ={90}^{\circ }$ (left) and ${\sigma }_n-{\tau }_{nt}$ at $\psi =0^{\circ }$ (right).

The $p_{\perp \perp }^t$, $p_{\perp \perp }^c$, $p_{\perp \parallel }^t$, $p_{\perp \parallel }^c$ are the inclination factors for composite material. The values of the inclination factor depend on the type of the composite material. The inclination factor values are obtained from the slope of the master fracture body ${\sigma }_n-{\tau }_{n1}$ at $\psi ={90}^{\circ }$ and ${\sigma }_n-{\tau }_{nt}$ at $\psi =0^{\circ }$ when ${\sigma }_n=0$ as illustrated in Figure 2. If the experimental data is unavailable, Puck provides the typical inclination factor values for the GFRP and CFRP in Table 1.

Each action plane at the corresponding angle will have its own unique failure exposure value. The angle of the action plane, when the failure exposure value is maximum, will define the FPOA of 3D Puck IFF for the current load. If the failure exposure value reaches 1 ($f_{E,\mathit{IFF}}\ge 1$), the material will have a permanent fracture plane at the corresponding angle.

3 Fracture Plane Orientation Angle Search Algorithm

In this section, the original fracture plane orientation angle (FPOA) search algorithm by Puck and the current fast FPOA search algorithms are explained. An IM7/8552 carbon fibre reinforced polymer (CFRP) material (Table 2) with varying load cases is used as the test case for the FPOA search algorithm. The result from each FPOA search algorithm is plotted into a fracture angle – $f_{E,\mathit{IFF}}$ value graph as a comparison.

Table 2 Material properties for IM7/8552 [14–17]

$E_{11}$	$E_{22}$, $E_{33}$	$G_{12},G_{13}$	${\nu }_{12}$, ${\nu }_{13}$	${\nu }_{23}$
171.42 GPa	9.08 GPa	5.39 GPa	0.32	0.52
$R_{\parallel }^t$	$R_{\parallel }^c$	$R_{\perp }^t$	$R_{\perp }^c$	$R_{\parallel \perp }$
2323.5 MPa	1200.1 MPa	62.3 MPa	199.8 MPa	92.3 MPa
${\mathcal{𝒢}}_{\parallel }^t$	${\mathcal{𝒢}}_{\parallel }^c$	${\mathcal{𝒢}}_{\perp }^t$	${\mathcal{𝒢}}_{\perp }^c$	${\mathcal{𝒢}}_{\parallel \perp }^t$
81.5 N/mm	106.3 N/mm	0.2774 N/mm	1.3092 N/mm	0.7879 N/mm

Table 3 Stress state examples to calculate failure exposure on 3D Puck IFF

Load
Case	${\sigma }_{22}$ (MPa)	${\sigma }_{33}$ (MPa)	${\sigma }_{12}$ (MPa)	${\sigma }_{13}$ (MPa)	${\sigma }_{23}$ (MPa)	Notes
1	0	30	5	10	$-$40	Max value of $f_{E,\mathit{IFF}}$ at ($-$56${}^{\circ }$, 0.929)
2	4	$-$38	20	$-$33	46	Two local max values of $f_{E,\mathit{IFF}}$ at ($-$7${}^{\circ }$, 0.668) and (72${}^{\circ }$, 0.6746)
3	$-$66	$-$55	$-$4	$-$13	70	Two local max values of $f_{E,\mathit{IFF}}$ with the global max value of $f_{E,\mathit{IFF}}$ located near 90${}^{\circ }$/$-$90${}^{\circ }$ at (85${}^{\circ }$, 0.745)

Three different load cases (Table 3) are selected to demonstrate the most common FPOA calculation result. Load case 1 is an example of one maximum value of $f_{E,\mathit{IFF}}$ iteration result. Load case 2 is an example of two local maximum values of $f_{E,\mathit{IFF}}$ results, with the global maximum value of $f_{E,\mathit{IFF}}$ nowhere near 90${}^{\circ }$ or $-$90${}^{\circ }$ angle orientation. Load case 3 is a specific case for two local maximum values of $f_{E,\mathit{IFF}}$ where the global maximum value of $f_{E,\mathit{IFF}}$ is located near 90${}^{\circ }$ or $-$90${}^{\circ }$ angle orientation. The load case 3 is used to demonstrate that other fast FPOAs have difficulties in accurately determining the global maximum value of $f_{E,\mathit{IFF}}$.

3.1 Existing FPOA Search Algorithm

3.1.1 Load case 1 – one maximum value of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$

In order to find the maximum value of the failure exposure and the corresponding action plane for the 3D Puck IFF, Puck utilises the Stepwise Search Method (SSM) algorithm. Puck SSM algorithm calculates all the failure exposure values that are obtained from $-$90${}^{\circ }$ $\le$ $\theta$ $<$ 90${}^{\circ }$ iteration range angle. All the failure exposure values are then compared to obtain the maximum failure exposure value. The SSM algorithm requires 180 supporting points for a 1${}^{\circ }$ inclination angle to predict the corresponding FPOA for the maximum failure exposure value for the given load, which requires a considerable amount of time to compute.

Figure 3 shows the value of $f_{E,\mathit{IFF}}$ and its corresponding inclination angle, totalling 180 supporting points. The inclination angle with the highest $f_{E,\mathit{IFF}}$ value will be used as the FPOA for the current stress state.

Figure 3 SSM algorithm on load case 1 with $f_{E,\mathit{IFF}}$ maximum value at ($-$56${}^{\circ }$, 0.929).

Weigand [8] proposes the EGSS (Extended Golden Section Search) algorithm to reduce the computation time when obtaining the maximum failure exposure and its corresponding FPOA on 3D Puck IFF. The EGSS algorithm utilises the golden section search (GSS) algorithm [9] to reduce the number of supporting points required to find the maximum failure value significantly. The GSS algorithm is used to iterate the supporting points sufficiently close until it reaches the prescribed tolerance value. A successive parabolic interpolation (Equation (10)) [10] is implemented between the last two GSS iteration supporting points to obtain the maximum failure exposure value. The subscripts in Equation (10) refer to the inclination angle and the corresponding $f_{E,\mathit{IFF}}$ value results of the last three supporting points, which are iterated using the GSS algorithm.

$${\theta }_{fp}={\theta }_2-\frac{1}{2}\frac{\begin{array}{c}{({\theta }_2-{\theta }_1)}^2(f_{E,IFF(2)}-f_{E,IFF(3)})\hfill \\ -{({\theta }_2-{\theta }_3)}^2(f_{E,IFF(2)}-f_{E,IFF(1)})\hfill \end{array}}{\begin{array}{c}({\theta }_2-{\theta }_1)(f_{E,IFF(2)}-f_{E,IFF(3)})\hfill \\ -({\theta }_2-{\theta }_3)(f_{E,IFF(2)}-f_{E,IFF(1)})\hfill \end{array}}$$

(10)

Figure 4 shows the same load case 1 calculation but using the EGSS algorithm. Compared to the SSM algorithm, the EGSS algorithm requires less supporting point calculation to obtain a similar result. As a consequence, the EGSS algorithm is faster than the SSM algorithm.

Figure 4 EGSS algorithm has similar result as SSM algorithm on load case 1.

3.1.2 Load case 2 – two local maximum values of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$

However, the EGSS algorithm tends to incorrectly identify the global maximum value of $f_{E,\mathit{IFF}}$ when there are two or more local maximum values of $f_{E,\mathit{IFF}}$. As an example, on load case 2, the SSM algorithm identifies two local maximum values of $f_{E,\mathit{IFF}}$ coordinates at ($-$7${}^{\circ }$, 0.668) and (72${}^{\circ }$, 0.6746). The SSM algorithm then determines that the second coordinate point is the global maximum value of $f_{E,\mathit{IFF}}$ (Figure 5) since the second coordinate has a higher $f_{E,\mathit{IFF}}$ value compared to the first coordinate point.

Figure 5 Load case 2 with the global maximum value of $f_{E,\mathit{IFF}}$ at (72${}^{\circ }$, 0.6746) using SSM algorithm.

As shown in Figure 6, the EGSS algorithm incorrectly determines the global maximum value of $f_{E,\mathit{IFF}}$ at ($-$6,4625${}^{\circ }$, 0.668). The result shows that EGSS’s global maximum value of $f_{E,\mathit{IFF}}$ is closer to the SSM algorithm’s first local maximum value of $f_{E,\mathit{IFF}}$ coordinate than the SSM algorithm’s second local maximum value of $f_{E,\mathit{IFF}}$, which is the true global maximum value of $f_{E,\mathit{IFF}}$.

Figure 6 EGSS algorithm fails to identify the correct global maximum value of $f_{E,\mathit{IFF}}$ on two local maximum value of $f_{E,\mathit{IFF}}$ case (load case 2).

Figure 7 The supporting points with its corresponding values calculated using Selective Range (SR) method to localise potential global maximum value of $f_{E,\mathit{IFF}}$.

Schirmaier [11] suggests a new search algorithm called the Selective Range Golden Section Search (SRGSS) algorithm to improve the EGSS algorithm by combining the Selective Range (SR) method and the GSS algorithm. According to Schirmaier, 20% of the $1\times {10}^5$ randomised stress cases evaluated using the EGSS algorithm have a 5${}^{\circ }$ deviation on its FPOA. Schirmaier observes that the $f_{E,\mathit{IFF}}$ – inclination angle result plot has distinct characteristics, such as smooth, may have up to three local maxima and the minimum distance between two local maxima is always larger than 25${}^{\circ }$.

Figure 8 SRGSS successfully identify the correct global maximum value of $f_{E,\mathit{IFF}}$ on two local maximum values of $f_{E,\mathit{IFF}}$ case (load case 2).

Schirmaier proposes the SR method, which reduces the wide 180${}^{\circ }$ search range interval into one or more narrow 20${}^{\circ }$ search range intervals to find the possibility of multiple local maxima on the $f_{E,\mathit{IFF}}$ – inclination angle result plot. The 20${}^{\circ }$ search range interval is determined by calculating the supporting points on every 10${}^{\circ }$ interval inclination angle. If the current supporting point has higher $f_{E,\mathit{IFF}}$ values than the previous and the next supporting point, then the span of inclination angle between the previous and the next supporting point is designated as the SR interval (Figure 7). After the SR interval is obtained, the GSS algorithm is applied to obtain the maximum $f_{E,\mathit{IFF}}$ value in each SR interval. If there are two or more SR intervals, the local maximum value of $f_{E,\mathit{IFF}}$ on each SR interval is compared to obtain the global maximum value of $f_{E,\mathit{IFF}}$ (Figure 8). This method greatly improves the accuracy of the GSS algorithm while maintaining less computational resources since it uses fewer supporting points to obtain the global maximum value of $f_{E,\mathit{IFF}}$ compared to the SSM algorithm.

Rezasefat [12] later improves the SRGSS algorithm by implementing the Simple Parabolic Interpolation Search (SPIS) algorithm (Figure 9). SPIS aims for less computation time compared to the SRGSS. SPIS uses parabolic interpolation (Equation (10)) instead of the GSS algorithm to reduce the required supporting points to obtain the global maximum value of $f_{E,\mathit{IFF}}$. The parabolic interpolation is only implemented into the SR interval range that has the highest average of two smaller $f_{E,\mathit{IFF}}$ values neighbouring supporting points, further reducing the total amount of the calculation time. As a final check, SPIS compares the value obtained from the parabolic interpolation with the $f_{E,\mathit{IFF}}$ value at the inclination angle $-$90${}^{\circ }$ and 90${}^{\circ }$.

Figure 9 SPIS successfully identify the correct global maximum value of $f_{E,\mathit{IFF}}$ on two local maximum value of $f_{E,\mathit{IFF}}$ case (load case 2).

3.1.3 Load case 3 – two local maximum values of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$ with the global maximum value of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$ located near $-$90${}^{\circ }$/90${}^{\circ }$

On load case 3, two local maximum values of $f_{E,\mathit{IFF}}$ coordinate at (9${}^{\circ }$, 0.735) and (85${}^{\circ }$, 0.745) are identified using the SSM algorithm. The SSM algorithm determines the second local maximum value of $f_{E,\mathit{IFF}}$ as the global maximum value of $f_{E,\mathit{IFF}}$ since it has a higher $f_{E,\mathit{IFF}}$ value compared to the other point. The location of the global maximum value of $f_{E,\mathit{IFF}}$ is near the edge of the Puck search range interval ($-$90${}^{\circ }$ $\le$ $\theta$ $\le$ 90${}^{\circ }$). The result of the fracture angle – $f_{E,\mathit{IFF}}$ value plot using the SSM algorithm is shown in Figure 10.

Figure 10 SSM algorithm on load case 3. The global maximum value of $f_{E,\mathit{IFF}}$ is located at the second local maximum value of $f_{E,\mathit{IFF}}$ at (85${}^{\circ }$, 0.745).

Even though load case 3 is similar to load case 2 (two local maximum values of $f_{E,\mathit{IFF}}$), the current fast FPOA search algorithm (EGSS, SRGSS and SPIS) incorrectly determines the global maximum value of $f_{E,\mathit{IFF}}$, as shown in Figure 11.

Figure 11 Current fast algorithm, EGSS, SRGSS and SPIS are fail to correctly identify the global maximum value of $f_{E,\mathit{IFF}}$.

3.2 Proposed Fast FPOA Search Algorithm

A novel fast FPOA search algorithm is proposed to improve the reliability of the current fast FPOA search algorithm in finding the global maximum value of $f_{E,\mathit{IFF}}$ and its corresponding fracture angle. The proposed fast FPOA search algorithm is named the Improved Selective Range Brent Method (ISRBM).

In the ISRBM search algorithm, an improved version of the SR method is implemented to rectify the unintentional bug in the SRGSS or SPIS algorithm. The SR method itself works by comparing the current supporting point with two points on both preceding and following supporting points. If the current supporting point has a higher $f_{E,\mathit{IFF}}$ value, the preceding and the following supporting points are regarded as the starting and the stopping points for the SR interval.

Figure 12 The supporting points are calculated within $-$90${}^{\circ }$ $\le$ $\theta$ $\le$ 90${}^{\circ }$ inclination angle on SRGSS or SPIS algorithm.

On the SRGSS and SPIS algorithm, the calculation of the supporting points for the SR method is limited to $-$90${}^{\circ }$ $\le$ $\theta$ $\le$ 90${}^{\circ }$ inclination angle [11, 12], as shown in Figure 3.10. Since the SR interval evaluation requires both preceding and following points, the SR interval evaluation will never initialise at the edge points ($-$90${}^{\circ }$/90${}^{\circ }$ inclination angle) since there is only one point to compare (either preceding or following point, not both).

Figure 13 The supporting points are calculated within $-$90${}^{\circ }$ $\le$ $\theta$ $\le$ 100${}^{\circ }$ inclination angle on ISRBM search algorithm. Note the additional calculation at 100${}^{\circ }$ inclination angle.

Figure 14 Brent Method is applied on the SR interval on ISRBM search algorithm. Brent Method calculated points are omitted for clarity.

The ISRBM search algorithm improves the SR method by appending an additional supporting point to ensure that the SR interval evaluation can be initialised even at the edge point ($-$90${}^{\circ }$/90${}^{\circ }$ inclination angle). The additional supporting point is calculated at a 100${}^{\circ }$ inclination angle. Since $f_{E,\mathit{IFF}}$ is periodic for every 180${}^{\circ }$ (the $f_{E,\mathit{IFF}}$ value at $\theta$ and $\theta$ $+$ 180${}^{\circ }$ will be identical), the $f_{E,\mathit{IFF}}$ calculation result that is obtained at $-$90${}^{\circ }$ and $-$90${}^{\circ }$ inclination angle can be reused as the $f_{E,\mathit{IFF}}$ result at 90${}^{\circ }$ and 100${}^{\circ }$ inclination angle respectively, thus further reducing unnecessary computation time. As a result, the SR method is implemented into a $-$90${}^{\circ }$ $\le$ $\theta$ $\le$ 100${}^{\circ }$ inclination angle on the ISRBM search algorithm (Figure 13) instead of the usual $-$90${}^{\circ }$ $\le$ $\theta$ $\le$ 90${}^{\circ }$ inclination angle.

The ISRBM search algorithm obtains the maximum value of $f_{E,\mathit{IFF}}$ on each SR interval using the Brent Method (Figure 14). Brent method [18] combines the more reliable GSS algorithm and the fast successive parabolic interpolation technique. The successive parabolic interpolation technique has a superlinear convergence rate but is prone to fail when the evaluated points are collinear. On the contrary, the GSS algorithm has a linear convergence rate but is more robust. In the Brent method, the default calculation uses successive parabolic interpolation, but when the evaluated points are collinear, the calculation will fall back to the GSS algorithm.

In the final step, when the corresponding fracture angle of the global maximum value of $f_{E,\mathit{IFF}}$ is between 90${}^{\circ }$ $\le$ $\theta$ $\le$ 100${}^{\circ }$ inclination angle, the fracture angle needs to be subtracted by 180${}^{\circ }$ due to the $f_{E,\mathit{IFF}}$ periodicity characteristic. The pseudocode of the ISRBM search algorithm is summarised in Algorithm 1.

---

Algorithm 1 Pseudocode of the ISRBM search algorithm

---

1: Calculate the $f_{E,\mathit{IFF}}$ at supporting points from $-$90${}^{\circ }$ $\le$ $\theta$ $\le$ 100${}^{\circ }$ inclination angle with 10${}^{\circ }$ interval. Reuse the $f_{E,\mathit{IFF}}$ calculation result of the supporting points to reduce the computation time since the $f_{E,\mathit{IFF}}$ result is periodic for every 180${}^{\circ }$ inclination angle

2: IF inclination angle between $-$90${}^{\circ }$ $\le$ $\theta$ $<$ 90${}^{\circ }$ THEN

3: |Calculate failure exposure at corresponding inclination angle $f_{E,\mathit{IFF}}(\theta )$

4: ELSE IF inclination angle is equal or more than 90${}^{\circ }$ THEN

5: |Assign $f_{E,\mathit{IFF}}(\theta -{180}^{\circ })$ result into $f_{E,\mathit{IFF}}(\theta )$

6: Find the appropriate selective range interval that contains local maximum value of $f_{E,\mathit{IFF}}$

7: Applying Brent Method into the selective range interval

8: Compare the local $f_{E,\mathit{IFF}}$ maximum values obtained from the previous step

9: Assign the highest local maximum value of $f_{E,\mathit{IFF}}$ into global maximum value of $f_{E,\mathit{IFF}}$

10: Assign the corresponding inclination angle of global maximum value of $f_{E,\mathit{IFF}}$ as the FPOA, ${\theta }_{fp}$

11: Check the ${\theta }_{fp}$ result

12: IF ${\theta }_{fp}$ is more than 90${}^{\circ }$ THEN

13: |Substract the fracture angle by 180${}^{\circ }$

---

4 Results

In order to evaluate the performance of the FPOA search algorithm, $1\times {10}^5$ randomised stress cases are tested. IM7/8552 CFRP material (Table 2) is used as the material input. Each FPOA search algorithm is evaluated in terms of the calculation time and reliability to obtain the correct calculation result. Table 4 shows the average calculation time to obtain the global maximum value of $f_{E,\mathit{IFF}}$ and its corresponding FPOA ${\theta }_{fp}$ for each stress case. In this table, the performance increase of each fast FPOA is also compared with the standard Puck SSM algorithm.

Table 4 Average calculation time for each stress case and the performance comparison of FPOA search algorithm to the Puck SSM 1${}^{\circ }$ increment search

Search	Average Calculation	Average Calculation Time
Algorithm	Time (Millisecond)	Decrease Compared to SSM (%)
SSM	0.59	–
EGSS	0.029	1934.48
SRGSS	0.079	646.84
SPIS	0.071	730.99
ISRBM	0.099	495.96

Next, the reliability of each fast FPOA search algorithm is also evaluated by comparing the global maximum value of $f_{E,\mathit{IFF}}$ and its corresponding FPOA ${\theta }_{fp}$ using the Puck SSM algorithm. However, instead of 1${}^{\circ }$, a more precise 0.01${}^{\circ }$ increment search is used to iterate the interval angle. The average difference result obtained from each fast FPOA search algorithm with the Puck SSM algorithm is shown in Table 5. In addition, the distribution of the FPOA result difference on each fast FPOA search algorithm for all $1\times {10}^5$ randomised stress cases is visualised in Figure 15. The FPOA result difference is obtained by taking the absolute value of the difference between the fast FPOA search algorithm result and the FPOA result from SSM algorithm with 0.01${}^{\circ }$ increment. Since the $f_{E,\mathit{IFF}}$ value is periodic every 180${}^{\circ }$, the maximum FPOA difference between the fast FPOA search algorithm and the Puck SSM algorithm is always less or equal to 90${}^{\circ }$.

Table 5 The comparison of the FPOA and global maximum value of $f_{E,\mathit{IFF}}$ difference for each fast FPOA search algorithm to the Puck SSM 0.01${}^{\circ }$ increment search

Search	Average Fracture Plane	Average Failure Exposure
Algorithm	Orientation Angle Difference (${}^{\circ }$)	Value Difference
EGSS	9.485	0.0156
SRGSS	2.991	0.0027
SPIS	0.5801	1.3322e-04
ISRBM	0.0048	1.0875e-07

Figure 15 The distribution of $1\times {10}^5$ randomised stress cases on the FPOA difference for each fast FPOA search algorithm compared to the Puck SSM 0.01${}^{\circ }$ increment search.

According to Table 4, the ISRBM search algorithm has the slowest computation time among the other fast FPOA search algorithms. However, Table 4 shows that the ISRBM search algorithm has the most reliable and consistent FPOA and failure exposure result since it has the average lowest difference result compared to the other FPOA search algorithms. Figure 15 captures the consistency of the ISRBM search algorithm to obtain FPOA. The figure shows that most of the ISRBM search algorithm calculation result has less than a 3${}^{\circ }$ FPOA difference when compared to the Puck SSM algorithm with a 0.01${}^{\circ }$ increment. In contrast, the other FPOA search algorithms have more scattered results. The slower computation time in the ISRBM search algorithm is the result of a more reliable and rigorous error-checking algorithm.

Before applying the ISRBM search algorithm to the finite element analysis (FEA) software, another validation is conducted by comparing the stress tensor at the material failure point between the Puck failure criterion and the experimental test. For this validation, the biaxial failure envelope ${\sigma }_{22}-{\tau }_{12}$ test result on e-glass/LY556 [19] is used as a benchmark. For comparison, both standard Puck SSM and ISRBM search algorithms are used alternatively to calculate the IFF failure exposure. At the IFF failure exposure $f_{E,\mathit{IFF}}=1$ (the irreversible damage on the material is started), the ${\sigma }_{22}-{\tau }_{12}$ result values are overlaid on top of the experimental test result. Figure 16 shows that the IFF Puck failure criterion is able to predict the result from the experimental result accurately.

Figure 16 The failure envelope ${\sigma }_{22}-{\tau }_{12}$ of 0${}^{\circ }$ uniderectional e-glass/LY566 epoxy biaxial test result.

Figure 17 Flowchart of the Puck failure criterion with proposed fast FPOA search algorithm subroutine.

In order to implement the Puck failure criterion in the FEA simulation, a separate subroutine is developed using Fortran with the proposed fast FPOA search algorithm embedded to reduce the computation time without sacrificing the calculation result. The subroutine flowchart for the Puck failure criterion with a fast FPOA search algorithm is shown in Figure 17. The subroutine is called for every Gauss/integration point.

An open-hole test (OHT) using AS4/PEEK CFRP composite material (Table 6) with [0/45/90/-45]${}_{2s}$ stacking sequence by Maa and Cheng [20] is used to validate the Puck failure criterion with fast FPOA search algorithm subroutine in the FEA simulation. The OHT specimen dimension has length l $=$ 100 mm, width w $=$ 20 mm, thickness t $=$ 2 mm, hole diameter d $=$ 5 mm and loading rate $\delta$ $=$ 2 mm/min. During the OHT experiment, Maa and Cheng put a clip gauge at the centre of the specimen with E $=$ 30 mm spacing to monitor the deformation of the specimen. For the FEA simulation, instead of using the full OHT model, a quarter OHT model with symmetric boundary condition applied at the xy-plane and xz-plane of the specimen is used to reduce the computation time [21], as shown in Figure 18.

The FEA simulation result demonstrates a good agreement with the experimental result, as shown in Figure 19. In addition, the ultimate strength difference result between the experimental (15.31 kN) and the FEA simulation (14.94 kN) is about 2.42%.

An example of the FEA simulation result on the OHT can be seen in Figure 20. This figure shows the plot of matrix damage in tension. The damage value in the plot was calculated using the Puck failure criterion and updated on each time increment on each Gauss/integration point. If the damage value reaches 1, the corresponding element will be deleted and will not be used for the next iteration.

Table 6 Material properties for AS4/PEEK [22, 23]

$E_{11}$	$E_{22}$, $E_{33}$	$G_{12},G_{13}$	${\nu }_{12}$, ${\nu }_{13}$	${\nu }_{23}$
127.6 GPa	10.3 GPa	6.0 GPa	0.32	0.49
$R_{\parallel }^t$	$R_{\parallel }^c$	$R_{\perp }^t$	$R_{\perp }^t$	$R_{\parallel \perp }$
2023.0 MPa	1234.0 MPa	92.7 MPa	176.0 MPa	82.6 MPa
	${\mathcal{𝒢}}_{\parallel }^t$, ${\mathcal{𝒢}}_{\parallel }^c$	${\mathcal{𝒢}}_{\perp }^t,{\mathcal{𝒢}}_{\perp }^c$	${\mathcal{𝒢}}_{\parallel \perp }^t$
	128.0 N/mm	5.6 N/mm	4.93N/mm

Figure 18 The dimension and the boundary conditions of a quarter OHT model.

Figure 19 The comparison between the experimental and FEA simulation results for the OHT specimen.

Figure 20 The plot of matrix damage in tension using Puck failure criterion on AS4/PEEK for 90${}^{\circ }$ laminate.

5 Conclusion

This paper proposes a reliable and fast FPOA search algorithm, the ISRBM search algorithm, to substitute the standard Puck SSM algorithm. Several improvements have been implemented into the ISRBM search algorithm, such as additional supporting points for more robust selective range interval testing, particularly at the FPOA near $-$90${}^{\circ }$/90${}^{\circ }$, and a more robust Brent method for finding the failure exposure maximum value within selective range interval. Albeit the computational performance of the ISRBM search algorithm is slower compared to other fast FPOA search algorithms, the ISRBM search algorithm calculation result is exceptionally reliable.

Data Availability

The subroutine is available upon request by emailing the first author.

Acknowledgment

This work was supported by the Indonesia Endowment Fund for Education (LPDP).

References

[1] S. W. Tsai and E. M. Wu, “A general theory of strength for anisotropic materials,” Journal of Composite Materials, vol. 5, pp. 58–80, 1971.

[2] Z. Hashin, “Failure criteria for unidirectional fiber composites,” Journal of Applied Mechanics, vol. 47, no. 2, pp. 329–334, 1980.

[3] M. J. Hinton, A. S. Kaddour and P. D. Soden, Failure criteria in fibre-reinforced-polymer composites, Elsevier, 2004.

[4] A. Puck, “Festigkeitsanalyse von Faser-Matrix-Laminaten, Modelle für die Praxis,” Carl Hanser Verlag, Munich, 1996.

[5] A. Puck and H. Schürmann, “Failure analysis of FRP laminates by means of physically based phenomenological models,” Composites Science and Technology, vol. 58, no. 7, pp. 1045–1067, 1998.

[6] A. Puck, J. Kopp and M. Knops, “Guidelines for the determination of the parameters in Puck’s action plane strength criterion,” Composites Science and Technology, vol. 62, no. 3, pp. 371–378, 2002.

[7] Verein Deutscher Ingeniere., “Development of fibre-reinforced plastics components – analysis,” Beuth Verlag GmbH, Berlin, 2006.

[8] J. Weigand, N. Petrinic and B. Elliott, “An algorithm for determination of the fracture angle for the three-dimensional Puck matrix failure criterion for UD composites,” Composites Science and Technology, vol. 68, no. 12, pp. 2511–2517, 2008.

[9] J. Kiefer, “Sequential minimax search for a maximum,” Proceedings of the American Mathematical Society, vol. 4, no. 3, pp. 502–506, 1953.

[10] P. Jarratt, “An iterative method for locating turning points,” The Computer Journal, vol. 10, no. 1, pp. 82–84, 1967.

[11] F. J. Schirmaier, J. Weiland, L. Kärger and F. Henning, “A new efficient and reliable algorithm to determine the fracture angle for Puck’s 3D matrix failure criterion for UD composites,” Composites Science and Technology, vol. 100, pp. 19–25, 2014.

[12] M. Rezasefat, D. B. Torres, A. Gonzales-Jimenez, M. Giglio and A. Manes, “A fast fracture plane orientation search algorithm for Puck’s 3D IFF criterion for UD composites,” Materials Today Communications, vol. 28, p. 102700, 2021.

[13] B. Paul, “A modification of the Coulomb-Mohr theory of fracture,” Journal of Applied Mechanics, vol. 28, no. 2, pp. 259–268, 1961.

[14] P. P. Camanho, P. Maimí and C. G. Dávila, “Prediction of size effects in notched laminates using continuum damage mechanics,” Composites Science and Technology, vol. 67, no. 13, pp. 2715–2727, 2007.

[15] P. P. Camanho, M. A. Bessa, G. Catalanotti, M. Vogler and R. Rofles, “Modeling the inelastic deformation and fracture of polymer composites – Part II: Smeared crack model,” Mechanics of Materials, vol. 59, pp. 36–49, 2013.

[16] P. Maimí, P. P. Camanho, J. A. Mayugo and C. G. Dávila, “A Continuum Damage Model for Composite Laminates: Part I – Constitutive Model,” Mechanics of Material, vol. 29, no. 10, pp. 897–908, 2007.

[17] P. Maimí, P. P. Camanho, J. A. Mayugo and C. G. Dávila, “A continuum damage model for composite laminates: Part II – Computational implementation and validation,” Mechanics of Material, vol. 39, no. 10, pp. 909–919, 2007.

[18] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, FORTRAN numerical recipes, New York: Cambridge University Press, 1996.

[19] P. D. Soden, M. J. Hinton and A. S. Kaddour, “Biaxial test results for strength and deformation of a range of E-glass and carbon fibre reinforced composite laminates: failure exercise benchmark data,” Composites Science and Technology, vol. 62, no. 12, pp. 1489–1514, 2002.

[20] R.-H. Maa and J.-H. Cheng, “A CDM-based failure model for predicting strength of notched composite laminates,” Composites Part B: Engineering, vol. 33, no. 6, pp. 479–489, 2002.

[21] A. M. Girão Coelho, J. T. Mottram and K. A. Harries, “Finite element guidelines for simulation of fibre-tension dominated failures in composite materials validated by case studies,” Composite Structures, vol. 126, pp. 299–313, 2015.

[22] J. F. Chen, E. V. Morozov and K. Shankar, “A combined elastoplastic damage model for progressive failure analysis of composite materials and structures,” Composite Structures, vol. 94, no. 12, pp. 3478–3489, 2012.

[23] O. Völkerink, E. Petersen, J. Koord and C. Hühne, “A pragmatic approach for a 3D material model considering elasto-plastic behaviour, damage initiation by Puck or Cuntze and progressive failure of fibre-reinforced plastics,” Computers & Structures, vol. 236, p. 106280, 2020.

Biographies

Nanda Wirawan is currently pursuing his PhD in mechanical engineering at the University of Sheffield. He is also working as a researcher at the Aeronautics and Space Research Organisation at the National Research and Innovation Agency of Indonesia. His research focuses on computational solid mechanics, primarily on the composite material.

Ibrahim H. Abuzayed is currently pursuing a Ph.D. in Mechanical Engineering at the University of Sheffield, having earned his Master’s in Advanced Mechanical Engineering from the same institution in 2020. His research focuses on computational mechanics, specifically the fracture of hybrid composite materials. Ibrahim has actively contributed to the academic community through teaching activities at the University of Sheffield since 2021.

Mahesa Akbar received the philosophy of doctorate degree in Mechanical Engineering from The University of Sheffield in 2019. He is currently working as a Researcher at the Centre of Defence and Security Technology and Lightweight Structure Research Group at Institut Teknologi Bandung. His expertise is in the field of computer aided engineering (CAE) with interest in structural dynamics, fluid dynamics, aeroelasticity, fluid-structure interaction, energy harvesting and smart composites.

Jose L. Curiel-Sosa is a lecturer at the University of Sheffield. He studied at the Escuela Superior de Ingenieros Industriales, University of Seville. He has a PhD in engineering modelling and simulation from Swansea University. He has previously worked for both industry (composites engineer in manufacturing and design) and academia (aerostructures postdoc – University of Oxford). His research is polymathic spanning from computational analysis to fundamental research; being applied in a range of industrial sectors, from aerospace (interaction aerostructures-aerodynamics, aeroelasticity, aviation safety) to manufacturing and machining (composite structures).

European Journal of Computational Mechanics, Vol. 33_3, 227–254.
doi: 10.13052/ejcm2642-2085.3332
© 2024 River Publishers

Abstract

1 Introduction

2 Puck Inter-Fibre Failure (IFF) Calculation

3 Fracture Plane Orientation Angle Search Algorithm

3.1 Existing FPOA Search Algorithm

3.1.1 Load case 1 – one maximum value of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$

3.1.2 Load case 2 – two local maximum values of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$

3.1.3 Load case 3 – two local maximum values of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$ with the global maximum value of ${\mathit{f}}_{\mathit{E},\mathit{IFF}}$ located near $-$90${}^{\circ }$/90${}^{\circ }$

3.2 Proposed Fast FPOA Search Algorithm

\01Kumaravel\07_pdf-html\Set-141\01out\ejcm\EJCM_33-3-Article-2\art2.html:

---

\01Kumaravel\07_pdf-html\Set-141\01out\ejcm\EJCM_33-3-Article-2\art2.html:

---

\01Kumaravel\07_pdf-html\Set-141\01out\ejcm\EJCM_33-3-Article-2\art2.html:

---

4 Results

5 Conclusion

Data Availability

Acknowledgment

References

Biographies