Optimization of Fiber-Reinforced Polymer Patches for Repairing Fatigue Cracks in Steel Plates Using a Genetic Algorithm

A total of 864 3D FE models were analyzed with the ABAQUS software version 2017 to obtain the SIF values for different combinations of steel plate configurations, crack lengths, materials and the geometric properties of the patch, and the adhesive layer. A 20-node quadratic solid element, C3D20, was used for both the cracked plate and adhesive layer. An 8-node continuum shell element with reduced integration, SC8R, was appointed to the FRP patch. The layer-wise theory of Reddy (1987) was applied for modeling the FRP patch by stacking SC8R elements throughout its thickness using the command *SHELL SECTION, STACK DIRECTION = 3. To capture the square root singularity of the stress and strain fields near the crack tip, the authors assigned a collapsed 3D element to elements in a small region around the crack front.

Surrounding the crack front are strips of wedge-shaped elements. These elements fill a semicylinder, radius, $R_e=[(a/12)-(a/5)]$, with a center at the crack front. The number of strips spanning along the radial direction of the semicylinder depends on the crack front element size, $L_e$. To determine an appropriate value of $L_e$, a sensitivity analysis (Hmidan et al. 2015) was implemented on the $6\times 90\times \mathrm{1,300}\mathrm{mm}$ center-cracked steel plate. In the sensitivity analysis, the calculated SIF values, $K_{FE}$, were obtained from different $a/L_e$ ratios (from 20 to 100) and two normalized crack lengths ($2a/W_s=0.1$ and 0.9). Fig. 2 shows a comparison between the calculated SIFs and referenced solutions, $K_{Ref}$ (Tada et al. 2000). For both normalized crack lengths, the $K_{FE}/K_{Ref}$ ratio converged to one as the $a/L_e$ ratio approached 100, and so $a/L_e$ ratio of 100 was chosen. The element size of $0.5\times 0.5\times 0.5\mathrm{mm}$ was used for the patch, adhesive layer, and steel plate region below the adhesive layer, while the element size of $1.5\times 1.5\times 1.5\mathrm{mm}$ was used in other regions of the steel plate.

Table 1 shows the material and geometric properties of steel plates, FRP patches, and adhesive layers used in FE models. Four different steel plates were used for two different purposes. Plates 1 ($6\times 90\times \mathrm{1,300}\mathrm{mm}$) and 2 ($16\times 180\times \mathrm{2,600}\mathrm{mm}$) were employed to create a database for developing the patching correction factor, $F_2$, while Plates 3 ($10\times 100\times \mathrm{1,500}\mathrm{mm}$) and 4 ($12\times 150\times \mathrm{2,200}\mathrm{mm}$) were used to obtain independent results for verification of the developed SIF solution. The elastic modulus ($E_s$), Poisson’s ratio (${\upsilon }_s$), and fatigue threshold SIF range of steel were assumed to be 200 GPa, 0.3, and 6.6 MPa $\sqrt{\mathrm{m}}$ (taken at a stress ratio of 0.13), based upon Pook (1975), respectively. Three unidirectional FRP materials of different longitudinal ($E_{p1}$) and transverse ($E_{p2}$) elastic moduli [Sika CarboDur M and H (SikaAG 2009) and MBRACE 460/1500 (Wu et al. 2012)] were chosen. Three adhesive materials [FM-73, FM36, and FM400 (Duong and Wang 2010)] of different elastic modulus ($E_a$) and Poisson’s ratio (${\upsilon }_a$) values were chosen. The steel and adhesive materials were assumed to be isotropic and linearly elastic, while the FRP materials were linear orthotropic under a plane stress condition. The FE models assumed a linear elastic behavior, which is reasonable for fatigue loading conditions.

Fig. 3 shows a typical quarter model of the center-cracked steel plate repaired with adhesive-bonded double-sided FRP patches. Displacement symmetric constraints (Xsym and Ysym) were applied to the corresponding symmetric planes. To enforce geometric compatibility conditions along the steel-adhesive and adhesive-patch interfaces, the authors employed the surface tie constraint. The tie constraint is a reasonable choice when the debonding effect is minimal, as observed in some experimental studies (Akbar et al. 2010; Wu et al. 2012; Yu et al. 2012; Wu et al. 2013).

A matrix of 864 rows was constructed from the FE models of two steel plates (Plates 1 and 2), three normalized crack lengths ($2a/W_s=0.1$, 0.5, and 0.9), three patch materials, four normalized patch widths ($W_p/W_s=0.2$, 0.6, 0.8, and 1), four normalized patch lengths ($L_p/2a=1$, 6, 11, and 16), and three adhesive materials [Table 1]. The length ($L_a$) and width ($W_a$) of the adhesive layer were equal to those of the patch. The calculated SIF value for each row occupied the last column of the matrix. The $F_2$ database was obtained by normalizing the calculated SIF values with the remote tensile stress, crack length, and finite-width correction factor.

The FE models of FRP-patched cracked plates were validated using data from past analytical and experimental studies (Kumar and Hakeem 2000; Ayatollahi and Hashemi 2007; Lam et al. 2010). Fig. 4 shows that the obtained SIF values are in good agreement with the published results for both single and double patches. In Fig. 4(a), the FE models can capture the singular stress field in cracked aluminum plates repaired with single-sided patches (Ayatollahi and Hashemi 2007), with an average error of 1.6%. For the cracked aluminum plates repaired with double-sided patches (Kumar and Hakeem 2000), the average error was 1.9% [Fig. 4(b)]. As shown in Fig. 4(c), the FE model can reasonably predict axial strains in the vicinity of a crack tip for both patched and unpatched sides of the test specimen R116 in Lam et al. (2010).

Genetic algorithms are iterative numerical solvers for optimization problems that were inspired by natural selection and natural genetics (Mitchell 1998). Each GA operates on a population of binary strings in the computer program. The GA starts by encoding random numeric values of independent variables in solution space to binary strings in the computer program with respect to the 1-to-1 mapping property. Each binary string represents one numeric point and vice versa. After the encoding process, the GA determines the fitness of each string in the current generation. The fitness is a value of the specified objective function at a particular point corresponding to the string under consideration. The algorithm then arranges all strings in a descending order of fitness values. Based on the order, GA sequentially performs the elite transfer, crossover, and mutation operators (Mitchell 1998) to produce a new population for the next generation. These operations are then repeated for each subsequent generation until the algorithm terminates at a specific number of generations.

The SIF of patch-repaired cracked plates depends on several design parameters. Conventional regression methods, which assume a mathematical model with unknown coefficients, are not efficient. A symbolic regression that simultaneously derives a mathematical model and unknown coefficients is preferred. In this study, GP is adopted as a tool for developing the SIF equation because it effectively minimizes residual errors between the source data and prediction. GP is an important application of GA for regression analysis (Koza 1992; Sette and Boullart 2001). A major difference between GP and GA is the representation of candidate solutions in which binary strings are used in GA and tree structures in GP (Fig. 5).

In this study, GP provides the best fit mathematical model for the patching correction factor ($F_2$) in Eq. (3). A mathematical model that correlates the numerical $F_2$ values with the GP approximations is obtained by using GP to maximize the squared Pearson correlation coefficient (Pearson’s $R^2$) in Eq. (5) (Lee Rodgers and Nicewander 1988)where $N$ = number of $F_2$ values normalized from the SIF database (864 in this case); $Y_i$ = $F_2$, obtained from the $i$th FE analysis; $Y_m$ = mean of all $Y_i$ values; $X_i$ = approximate value of $Y_i$ (randomly generated by GP); and $X_m$ = mean of all $X_i$ values.

The GP analysis consists of the following steps. In the first step, the initial population with a specific number of tree structures is randomly created using atoms from a given function set ($F$) and terminal set ($T$). Atoms of the $F$ set can be arithmetic operations, mathematical functions, Boolean operations, conditional operators, or any user-defined function (Koza 1992). Within a tree structure, $F^{\prime }s$ atoms occupy functional nodes that have two arguments, as shown in Fig. 5. The $T^{\prime }s$ atoms are independent variables and constants (Koza 1992) that are located at terminal nodes with no argument. In the second step, GP determines the $R^2$ values for all existing tree structures in the current generation. Based on the $R^2$ values, the algorithm stores all tree structures in a column vector and arranges them in descending order of their $R^2$ values. To produce a new population of tree structures for the next generation, GP sequentially implements three genetic operators: elite transfer [Fig. 6(a)], crossover [Fig. 6(b)], and mutation [Fig. 6(c)]. As a GA’s family member, GP’s solution is determined as the best-so-far tree structure stored in the cache of the computer program. The solution can be improved with a longer running time of the algorithm. For GP analysis, this study uses the HeuristicLab environment developed by Wagner et al. (2014). Different GP analyses can provide different SIF solutions. However, all solutions are capable of predicting SIF values if they are derived from the same database and GP’s input (Do and Lenwari 2018).

In the HeuristicLab, the $F_2$ database was randomly shuffled into four groups (1, 2, 3, and 4) to avoid any bias of the GP performance on certain groups. Each group had four subgroups (A, B, C, and D), each of which contained 216 data points (25% of the database) (Table 2). Subgroups A, B, C, and D were then sequentially used as a test set of 25% of the database, while in each case, the other three subgroups (75% of the database) formed the training set. Thus, four different $F_2$ functions were obtained from the four data groups, and the function with the highest $R^2$ value was chosen.

The function and terminal sets were $F=\{+,-,*,\mathrm{exponential},\mathrm{square},\mathrm{power}\}$ and $T=\{x_1,x_2,x_3,x_4,[-10,10]\}$, while the independent variable atoms were $x_1$, $x_2$, $x_3$, and $x_4$, as given in Eq. (4). Constant atoms were initially generated within an interval $[-10,10$]. Additional control parameters defined at the beginning of each GP analysis are listed in Table 3. The termination for all GP analyses was at the 2000th generation. All GP analyses were performed using an Intel Core i7-7700HQ CPU @ 2.81 GHz.

Fig. 7 shows that the training Pearson’s $R^2$ for the four GP analyses was either improved or remained constant after each GP generation. At the 2000th generation, the training $R^2$ values corresponding to the four groups were 0.902, 0.901, 0.900, and 0.906, respectively. The value of $R^2$ higher than 0.9 was an indicator of a good correlation for the training set. Fig. 8 shows the scatter diagram of $F_2$ at the 2000th generation in the case of $R^2=0.906$ (Group 4) for training data ($R^2=0.895$ for test data). The mean absolute error for the training data was 0.059. Although the FE and GP results were correlated, a discrepancy between the FE and GP results was observed at high $F_2$ values, when a crack is not entirely covered with the patch. Therefore, the SIF solution is recommended for complete crack coverage.

Fig. 9 shows the chosen GP’s tree structure with training $R^2=0.906$ for the patching correction factor ($F_2$) in Eq. (3). The corresponding mathematical function is given in Eq. (6)where $x_1$, $x_2$, $x_3$, and $x_4$ are the normalized parameters as defined in Eq. (4); and $c_0–c_9$ are the constant coefficients, as shown in Fig. 9. The application range of Eq. (6) is shown in Table 4. When the design parameters are beyond the range of applicability in Table 4, the FE and GP analyses in accordance with the previous sections can be performed to generate the new correction factor function.

Independent FE analyses and fatigue experiments in the literature were used to verify the proposed SIF solution. Plates 3 and 4 with four combinations of crack lengths, material and geometric properties of the FRP patch, and an adhesive layer were used in the FE analyses. Fig. 10 shows a comparison of the $F_2$ factors predicted by Eq. (6) with the independent FE results, where the predicted $F_2$ values were in good agreement with the FE results. For Plate 3, the average differences between the GP and FE solutions for different crack lengths, patch lengths, and patch widths at $2a/W_s=0.5$ or 0.9 were 3%, 5%, 3%, and 5%, respectively. For Plate 4, the differences were 2%, 4%, 2%, and 7%, respectively. However, the proposed SIF solution became less accurate when a crack was not entirely covered with the patch, and a large discrepancy was observed at a crack length of $2a=0.9W_s$, as shown in Fig. 10(d).

In addition, the fatigue life of four specimens in the previous experimental studies was estimated using the proposed SIF solution. Fig. 11 shows a comparison between the predicted and measured fatigue crack growth curves. Specimens CS-1 and CS-2 (stress $\mathrm{range}=120\mathrm{MPa}$) were tested by Wang et al. (2016b), while specimens D1L and D2L (stress range $\mathrm{I}=135\mathrm{MPa}$) were tested by Emdad and Al-Mahaidi (2015). The Paris law constants used in the fatigue life calculation were $C=2.427\times {10}^{-12}$ (for MPa and m units), and $m=3.3$ (Maiti 2015). The fatigue life of a specimen with a crack growing from an initial crack length, $a_0$, to a critical one, $a_c$, was predicted by Eq. (7)where $\mathrm{\Delta }{K}_{eff}$ = effective range of SIF that accounts for crack closure effects (MPa $\sqrt{\mathrm{m}}$); and $\mathrm{\Delta }{K}_{eff}$ was assumed to be equal to $\mathrm{\Delta }K$.

As shown in Fig. 11, a reasonable agreement was observed when the debonding effect was minimal. Therefore, the proposed SIF solution is applicable for predicting the fatigue life of FRP-patched cracked steel plates. The failure mode of specimen D1L was patch rupture, while that for specimen D2L was debonding. Debonding of the FRP sheet surrounding the fatigue crack was reported for both specimens CS-1 and CS-2, where the debonding area increased with an increasing number of cycles, but the CFRP sheets did not separate completely from the steel plates when the specimens fractured.

When the crack length, FRP and adhesive material properties, and adhesive thickness are specified, the correction factor $F_2$ in Eq. (6) will be a function of the FRP patch geometry. The SIF in Eq. (3) is rewritten as Eq. (8)where the terms $X_1$ to $X_3$ are derived as shown in Eq. (9)

The optimization statement minimizes the patch volume, $V_p(\mathbf{X})=X_1{X}_2{X}_3$, with respect to the constraint $\mathrm{\Delta }K(\mathbf{X})\le \mathrm{\Delta }{K}_{th}$, where $\mathrm{\Delta }{K}_{th}$ is the fatigue threshold SIF range of steel. The vector of design parameters is within the lower (${\mathbf{X}}_L$) and upper (${\mathbf{X}}_U$) bounds shown in Eq. (10)

In Eq. (10), the lower bound of the patch length is $2a$, which is the lower bound used in the GP analysis. It is recommended that the effective bond length be used as the lower bound in the optimization problem. The effective bond length is a function of the thickness and elastic modulus of the steel plate, FRP patch, and adhesive layer (Nozaka et al. 2005; Bocciarelli et al. 2007; Duong and Wang 2010).

Fig. 12(a) shows three $F_2$ surfaces corresponding to $t_p=2\mathrm{mm}$ (Surf 1), $t_p=1.2\mathrm{mm}$ (Surf 2), and fatigue threshold SIF range (Surf 3), respectively. The surfaces can help visualize the inequality constraint for the SIF range. When the patch thickness increases, SIF decreases. Therefore, the solution space, i.e., the underthreshold area, expands from the smallest space [Fig. 12(b)] to the largest one [Fig. 12(c)] as the patch thickness increases from 1.2 to 2 mm. If the crack length increases, Surf 3 may lie under Surf 1 and Surf 2, and no feasible solution can be obtained.

The performance of GA (command ga) was compared with nonlinear programming (command fmincon) in MATLAB version 2018a (MathWorks 2018). Both optimization solvers do not warrant the global minima. For the purpose of direct comparison, both solvers used an identical input of an initial point of the upper bound, function tolerance of ${10}^{-20}$, constraint tolerance of ${10}^{-15}$, and a maximum number of generations of 100.

A center-cracked steel plate ($t_s=10\mathrm{mm}$, $W_s=90\mathrm{mm}$, and $L_s=500\mathrm{mm}$) was considered in the design example. The elastic modulus ($E_s$) and Poisson’s ratio (${\upsilon }_s$) of steel were 200 GPa and 0.3, respectively. The fatigue threshold SIF range ($\mathrm{\Delta }{K}_{th}$) was 6.6 MPa $\sqrt{\mathrm{m}}$ or 209 MPa $\sqrt{\mathrm{mm}}$, (Pook 1975). The steel plate was subjected to constant amplitude cyclic loading with a tensile stress range of $\mathrm{\Delta }\sigma ={\sigma }_{\max }-{\sigma }_{\min }=55\mathrm{MPa}$ (${\sigma }_{\max }=55\mathrm{MPa}$ and ${\sigma }_{\min }=0\mathrm{MPa}$), and two normalized crack lengths of $2a=0.2W_s$ and $0.3W_s$ were considered.

A total of nine combinations of patch and adhesive materials were used to investigate the effects of the material properties on the optimum solutions at two different crack lengths. For each combination, GA was performed twice to examine its stochastic property.

Based on the lower and upper bounds of the design parameters, the GA solver randomly generated a total of 500 binary strings in which each string had a particular fitness value (patch volume value). Therefore, there were 500 different fitness values at each GA generation.

Fig. 13 shows the convergence histories of the 1st and 2nd GA, and the fmincon analyses for the case when $2a=0.2W_s$, $E_{p1}=460\mathrm{GPa}$, and $E_a=\mathrm{2,944}\mathrm{MPa}$. The best and mean fitness refers to the lowest and mean values of 500 fitnesses, respectively. Although the mean fitness values at the beginning of both analyses were much higher than the optimum ones, the GA could quickly achieve the solutions by the 8th and 7th generation for the 1st and 2nd GA analyses, respectively. The difference between the results of both GA analyses was 0.5%. Meanwhile, fmincon required 23 iterations to achieve a better solution (0.6% lower than the GA solutions).

In the case of a longer crack ($2a=0.3W_s$, $E_{p1}=460\mathrm{GPa}$, and $E_a=\mathrm{2,944}\mathrm{MPa}$), the same conclusions could be drawn. Fig. 14 shows that the way GA searched for the optimum solution was smoother than fmincon. Again, fmincon required more iterations to achieve the optimum solution, which was 2.2% lower than the GA solutions. Note that although Figs. 13 and 14 show that the patch volumes obtained from fmincon were the lowest after the 7th and 2nd iteration, respectively, these are not the minima because the first Karush-Kuhn-Tucker condition (Rao and Rao 2009) was violated.

Table 5 summarizes the optimum patch solutions for different patch and adhesive materials at two normalized crack lengths. It was found that the GA solution provided an equal or slightly higher patch volume than the fmincon solution. The maximum differences between the two solvers were 0.6% and 2.2% for $2a$ values of 0.2 $W_s$ and 0.3 $W_s$, respectively. Also, the GA’s stochastic property was not pronounced, and the difference between the two GA solutions was trivial.

In the case of a shorter crack ($2a=0.2W_s$), the optimum patch width varied from 0.7- to 0.9-fold the width of the steel plate, while the optimum patch length varied from 4.4- to 6.1-fold the crack length. The optimum patch thickness was 1.2 mm. The optimum patch volume was $\mathrm{10,566}{\mathrm{mm}}^3$ when $E_{p1}$ was 210 GPa. It decreased by 18.4% and 40.1% when $E_{p1}$ were 300 and 460 GPa, respectively. The lowest patch volume ($\mathrm{6,347}{\mathrm{mm}}^3$) was found when $E_{p1}$ was 460 GPa and $E_a$ was 2,944 MPa. In addition, the optimum patch volume was found to be relatively insensitive to the adhesive modulus. The optimum patch volume slightly decreased from 0.1% to 0.2% when $E_a$ increased from 89% (from 959 to 1,815 MPa) to 207% (from 959 to 2,944 MPa).

In the case of a longer crack ($2a=0.3W_s$), no feasible solution existed when $E_{p1}$ was 210 GPa because the threshold surface Surf 3 lay below Surf 1 and Surf 2 (Fig. 12). However, an optimum patch volume could be achieved if the patch modulus $E_{p1}$ increased from 210 to 300 or 460 GPa. As $E_{p1}$ increased from 300 to 460 GPa, the patch volume decreased about two-fold. Similarly, the effect of the adhesive modulus on the optimum patch volume was negligible.

To assess the possibility of patch rupture and debonding failures, the authors analyzed the optimum patch-repaired cracked plates in Table 5 at the maximum fatigue load. The rupture failure of each FRP patch was assessed using the Tsai-Hill failure criterion (Staab 2015), as shown in Eq. (11)where $T_s$ = Tsai-Hill failure index; ${\sigma }_{11}$, ${\sigma }_{22}$, and ${\sigma }_{12}$ = longitudinal, transverse, and shear stresses in the patch, respectively; and $S_{11}$, $S_{22}$, and $S_{12}$ = ultimate longitudinal, transverse, and shear strength of the patch, respectively. The CFRP material strength in each direction is given in Table 6. The FRP rupture was deemed to have occurred when Eq. (11) was satisfied.

In recent years, the cohesive zone modeling approach has been employed to analyze adhesive-bonded joints, especially when the debonding behavior is simulated (Chen and Qiao 2009; Campilho et al. 2011, 2013; Chaves et al. 2014; Colombi et al. 2015; Zheng and Dawood 2016; Budhe et al. 2017). Typically, test data are required for cohesive parameters used in the traction-separation law. Due to its simplicity, the strength-based criterion was chosen to assess the debonding failure in this study. This criterion has been used in the assessment under static and fatigue loading conditions (Papanikos et al. 2005; Bocciarelli et al. 2009). However, the prediction can be sensitive to the element size due to the singularity condition. Small element sizes are required when this criterion is used (Da Silva et al. 2011). The debonding failure was assumed when at least one of the following three conditions was satisfied (Papanikos et al. 2005). For the first condition, the failure in the adhesive layer based on the maximum shearing stress at the steel-adhesive interface occurred when the Tresca reached a maximum, determined from Eq. (12)

For the second condition, the debonding at the steel-adhesive interface, based on the maximum normal stress, occurred when Eq. (13) was satisfied

For the third condition, the debonding at the adhesive-patch interface, based on the maximum normal stress, occurred when Eq. (14) was satisfied

In Eqs. (12)–(14), the terms $T_{ay}$, $T_{as}$, and $T_{ap}$ = adhesive failure indexes; ${\sigma }_1$ and ${\sigma }_3$ = maximum and minimum principal stresses, respectively; ${\sigma }_{33}$ = normal stress (peeling stress); and $p_{ay}$, $p_{as}$, and $p_{ap}$ = shear strength of adhesive material, peeling strength of steel-adhesive interface, and peeling strength of adhesive-patch interface, respectively. The values for $p_{ay}$ and $p_{as}$ (Table 6) were taken from Deng et al. (2018). When $p_{as}=p_{ap}$ was assumed, $T_{as}$ was higher than $T_{ap}$ because the maximum normal stresses at the steel-adhesive interface obtained from the FE analyses were higher than ones at the adhesive-patch interface.

Table 7 summarizes the patch and debonding failure indexes for the optimum patch solutions. Because all failure indexes ($T_s$, $T_{ay}$, and $T_{as}$) were less than one, the FRP rupture and debonding failures were not possible for optimum patch solutions at the maximum fatigue load. The Tsai-Hill failure indexes ($T_s$) for patch rupture were low in most cases, while the Tresca index for the adhesive layer [Eq. (12)] became pronounced as the adhesive modulus increased [$T_{ay}=0.58$ for $2a=0.2W_s$ and 0.64 for $2a=0.3W_s$, respectively]. In this example, the $T_{ay}$ values were higher than $T_{as}$ values. If the loading magnitude increased, debonding failure at the adhesive layer could be possible. For the optimum patch solutions, a good correlation between the SIF ranges obtained from ABAQUS and Eq. (3) was observed.

Figs. 15 and 16 show the effects of the adhesive modulus on Tresca and interfacial stresses in the adhesive layer for the two normalized crack lengths of $2a=0.2W_s$ and $0.3W_s$, respectively. The maximum Tresca, shear, and normal stresses in the adhesive layer increased as the adhesive modulus increased.

For this design example, the use of a high modulus patch with a low modulus adhesive was recommended. A high modulus patch caused more load to be transmitted to the patch, while it did not significantly influence the maximum Tresca in the adhesive layer (Table 7). Meanwhile, a low modulus adhesive reduced the possibility of debonding.