Prediction of fatigue crack growth rate, maximum strain energy release rate, and total energy dissipation in adhesive bonds: a machine learning approach and interactive R-Shiny panel development

Computing infrastructure and software environment

All computational analyses and model executions were performed on a workstation equipped with an Intel Core i5-4460 CPU @ 3.20GHz, 32GB of RAM, an AMD Radeon R7 200 Series graphics card (2GB VRAM), and a 240GB SSD. The operating system used for all computations was Windows 10 Education. All data analysis and ML model implementations were carried out using R statistical software (version 4.4.3; R Core Team, 2025) within the RStudio integrated development environment (version 2025.05.0 Build 496; Posit Software, PBC).

Key R packages utilized for the analyses include:

• tidyverse (version 2.0.0) and its components such as ggplot2 (version 3.5.2) for data visualization, dplyr (version 1.1.4) for data manipulation, tidyr (version 1.3.1), readr (version 2.1.5), purrr (version 1.0.4), stringr (version 1.5.1), forcats (version 1.0.0), and tibble (version 3.2.1).

• Boruta (version 8.0.0) for feature selection.

• caret (version 7.0-1) for general machine learning workflows, including model training and evaluation.

• Metrics (version 0.1.4) for performance metrics calculation.

• corrplot (version 0.95) for correlation matrix visualization.

• rpart (version 4.1.24) and rpart.plot (version 3.1.2) for decision tree modeling and plotting.

• DT (version 0.33) for interactive tables.

• shiny (version 1.11.1), shinyWidgets (version 0.9.0), colourpicker (version 1.3.0), and shinycssloaders (version 1.1.0) for potential interactive applications.

• Other utilized packages include RColorBrewer (version 1.1-3), lubridate (version 1.9.4), and lattice (version 0.22-6).

Experimental data and specimens

The foundation of this study is a comprehensive experimental dataset meticulously compiled by John-Alan Pascoe as part of his doctoral research (Pascoe, 2016). This dataset (Pascoe, Alderliesten & Benedictus, 2016) originates from a series of quasi-static and fatigue crack growth experiments conducted on DCB specimens. These specimens are a standard and critical methodology for assessing the crack propagation resistance and long-term durability of adhesive bonds under sustained loads.

The specimens were fabricated by bonding Aluminum 2024-T3 adherends with FM94 epoxy adhesive. Aluminum 2024-T3 is widely utilized in aerospace applications due to its high strength to weight ratio, while FM94 epoxy adhesive is frequently chosen for metal to metal bonding scenarios requiring high strength and fatigue durability.

The dataset encompasses raw measurements including crack length (a), applied force (F), and displacement (d). In addition to these raw measurements, the dataset also provides processed (derived) data, such as strain energy release rates (G) and energy dissipation per cycle. The G is a critical indicator within fracture mechanics for aerospace engineering, and cyclic energy dissipation is used to elucidate the performance of adhesive bonds under repetitive loading conditions. Key target variables extracted from this dataset for prediction were da/dN, Total_Energy, and Gmax.

Descriptions of the 15 attributes present in the dataset are provided in Table 1.

Dependent variables are typically associated with crack propagation (dadN, G_max), energy consumption (Total_Energy, Cyclic_Energy), and deformation (d). Independent variables generally include load parameters (F and R), time and cycle (t and N), and material properties (C and a). Load ratio (R) is one of the fundamental factors influencing crack growth in fatigue tests. Delta_sqrt(G) is utilized as a mathematical transformation of the energy release rate and is typically modeled as an independent variable. In light of this information, other independent variables associated with the dependent variables, along with the purpose of the resulting model, are specified in Table 2.

Data analysis and ML methodology

The analytical framework of this study, implemented within an R-Shiny application, involves several sequential steps, including data preprocessing, feature selection, model training, and performance evaluation. The R-Shiny interface (Fig. S2) further provides an interactive reporting platform for comprehensive visualization and export of the analysis results.

Data preprocessing

Upon loading the raw CSV file, an initial preprocessing step involves converting suitable columns to numeric format, handling potential comma separated decimal values. Missing values (NA) in critical columns, namely N, F, d, a, da/dN, G_max, Delta_sqrt.G., R, Cyclic_Energy, Monotonic_Energy, and Total_Energy, are addressed by removing rows containing incomplete cases to ensure data integrity for subsequent analysis.

Outlier management is an optional but crucial step, controlled by user input, to mitigate the impact of extreme values. Outliers are identified based on a user defined Z-score threshold. The application offers three methods for outlier handling: removing the entire row, capping the outlier values to the Z-score bounds, or replacing them with NA values. This robust preprocessing ensures the quality and reliability of the data fed into the modeling pipeline.

Feature selection

Boruta algorithm (Kursa, Jankowski & Rudnicki, 2010) is employed for automatic and robust feature selection. Boruta, a wrapper algorithm built around Random Forest, iteratively compares the importance of original features with that of randomly shuffled “shadow features” to identify statistically significant predictors for the target variable. This step helps to reduce dimensionality, remove irrelevant or redundant features, and improve the interpretability and performance of the predictive models. The specific set of potential predictors fed to Boruta is dynamically determined based on the selected target variable (da/dN, Total_Energy, or Gmax) and the availability of these features in the cleaned dataset.

Data splitting

The preprocessed and feature-selected data is then partitioned into training and testing sets. A user-defined split ratio (typically between 60% and 90%) determines the proportion of data allocated to the training set, with the remainder reserved for testing. A fixed random seed is utilized to ensure reproducibility of the data partition across multiple runs.

ML model training and evaluation

This study employed a comprehensive suite of 16 different ML regression algorithms for predicting the selected target parameters. The selection of these 16 specific models was strategic, designed to provide a comprehensive comparative analysis across a wide spectrum of algorithmic families. The rationale was to evaluate the trade-offs between performance, complexity, and interpretability.

The models include:

• (1)

  simple, interpretable linear baselines (e.g., lm, ridge, lasso) to establish a performance benchmark;

• (2)

  powerful, non-linear tree-based ensembles (e.g., rf, gbm, xgbTree) renowned for their high accuracy;

• (3)

  flexible kernel-based and instance-based methods (e.g., svmRadial, knn); and

• (4)

  advanced approaches like Gaussian Process Regression (gaussprRadial) and Neural Networks (brnn, qrnn), which are adept at capturing complex patterns and quantifying uncertainty. These models range from fundamental linear approaches to complex ensemble and kernel-based methods, allowing for a thorough comparative analysis of their predictive capabilities. Each model was trained using k-fold cross-validation, with the number of folds being a user-defined parameter. It is important to note that hyperparameter tuning was an integral part of this training process. The caret package’s train function inherently performs automated hyperparameter optimization via a grid search approach during the k-fold cross-validation. For each algorithm, the function evaluates a predefined set of key hyperparameters (e.g., the mtry parameter for RF) and selects the combination that yields the best performance (e.g., lowest Root Mean Squared Error) on the cross-validation folds. This ensures that each model is compared in its optimized state, improving the fairness and transparency of the analysis.

Linear Regression (lm) is a fundamental regression algorithm that models the linear relationship between a dependent variable and one or more independent variables, forming the basis of statistical modeling (Gavss, 1809; Su, Yan & Tsai, 2012; Hu et al., 2025). Ridge Regression (ridge) is a linear model that adds L2 regularization to the cost function to address multicollinearity and reduce overfitting (Hoerl & Kennard, 1970; McDonald, 2009). Lasso Regression (lasso) is another linear model that performs feature selection by adding L1 regularization to the cost function, effectively shrinking the coefficients of irrelevant features to zero (Tibshirani, 1996; Abdullah et al., 2024). Elastic Net Regression (glmnet) is a hybrid linear model that combines the advantages of both Ridge and Lasso regression by simultaneously employing L1 and L2 regularization (Kumar & Sinha, 2025). Support Vector Machine Linear (svmLinear) is an extension of Support Vector Machines, developed by Cortes and Vapnik, which performs regression by finding a linear hyperplane to separate data, proving effective particularly in high-dimensional spaces (Cortes & Vapnik, 1995; Schölkopf, Burges & Smola, 1999). Support Vector Machine Radial (svmRadial) is a flexible and powerful variant of SVM that uses kernel functions, specifically the Radial Basis Function (RBF) kernel, to model non-linear relationships in the data (Cortes & Vapnik, 1995; Wu et al., 2024; Rathore et al., 2025). Decision Tree (rpart) is an interpretable and intuitive model that makes predictions by hierarchically partitioning the data based on feature values (Breiman et al., 2017; Konda, Verma & Jayaganthan, 2021). It was popularized as the “CART” (Classification and Regression Trees) algorithm by Breiman, Friedman, Olshen, and Stone (Breiman et al., 2017). Random Forest (rf) is a high-performance, ensemble learning algorithm that operates by constructing a multitude of decision trees at training time and outputting the mean prediction of the individual trees to reduce overfitting (Breiman, 2001; Huang et al., 2025). Gradient Boosting Machine (gbm) is a powerful ensemble algorithm that sequentially combines weak learners, with each new learner correcting the errors of its predecessors (Friedman, 2001; Aydın et al., 2025). Extreme Gradient Boosting Machine (xgbTree) is an optimized implementation of the gradient boosting framework, known for its superior speed and performance, particularly effective with structured and unstructured data (Chen & Guestrin, 2016; Peng et al., 2024; Huang et al., 2025). K-Nearest Neighbors (knn) is a non-parametric model that predicts the value of a new data point based on the average or weighted average of the target values of its knn in the training set (Fix & Hodges, 1951; Wang et al., 2025). Partial Least Squares (pls) is a dimension reduction and regression technique that analyzes multivariate relationships by constructing new components (latent variables) that best explain the variance in both independent and dependent variables (Shan et al., 2025). Generalized Additive Model is a flexible extension of linear models that allows for modeling non-linear but additive relationships between independent variables and the response using smoothing functions (Hastie & Tibshirani, 1986; Salkind, 2013). Bayesian Regularized Neural Network is a neural network model that incorporates a Bayesian approach to conventional artificial neural networks, adding prior information on the distribution of weights and helping to prevent overfitting (Chen et al., 2024). Gaussian Process Regression (Radial) (gaussprRadial) is a powerful non-parametric model that defines a probability distribution over functions by fitting to observed data, particularly strong in estimating uncertainty quantification for small to medium-sized datasets (Rasmussen & Williams, 2005; Han et al., 2025). Quantile Regression Neural Network (QRNN) is a neural network-based model that predicts not only the mean but also specific quantiles of the dependent variable, providing richer information about the conditional distribution of predictions (Koenker & Bassett, 1978; Han et al., 2025).

Each selected model is trained using k-fold cross-validation (CV), with the number of folds being a user-defined parameter. CV is crucial for obtaining robust model performance estimates and mitigating overfitting. For certain models (e.g., gaussprRadial, qrnn), pre-processing steps such as centering and scaling are applied to the data within the caret::train function to optimize model performance.

In a regression model, the term Y_i represents the actual, experimentally measured, or known value for a given data point (or observation), while Ŷ_i denotes the value predicted by the model for that same data point. The variable n signifies the total number of observations in the dataset, and μ is the mean of the actual values of the dependent variable. Given these definitions, the metrics are calculated as follows:

Root Mean Squared Error (RMSE) measures the average magnitude of the errors and gives greater weight to larger errors. It is expressed in the same units as the target variable. A smaller RMSE indicates a better fit of the model to the data, meaning the predictions are closer to the actual values. Conversely, a larger RMSE suggests less accurate predictions. RMSE is calculated using the following Eq. (1).

(1) $$RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n\left(Y_i-{\hat{\mathrm{Y}}}_i\right)}}.$$The coefficient of determination (R²) indicates how much of the variance in the dependent variable is explained by the model. It ranges from 0 to 1, with values closer to 1 suggesting an excellent fit of the model to the data. A higher R² value implies that a larger proportion of the variance in the actual data is accounted for by the model, indicating stronger predictive power. A lower R² value (approaching 0 or negative) suggests that the model does not explain much of the variance, or performs worse than simply predicting the mean. The R² value is determined by the following Eq. (2).

(2) $$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(Y_i-{\hat{\mathrm{Y}}}_i)}^2}{{\sum }_{i=1}^n{(Y_i-\mu )}^2}.}$$Mean Absolute Error (MAE) is the average of the absolute differences between predictions and actual observations. It is less sensitive to outliers than RMSE and directly expresses the average magnitude of errors. A smaller MAE indicates higher accuracy of the predictions, while a larger MAE suggests less accurate predictions on average. MAE is calculated using the following Eq. (3).

(3) $$MAE={\displaystyle \frac{1}{n}\sum _i^n|Y_i-{\hat{\mathrm{Y}}}_i.|.}$$Mean Absolute Percentage Error (MAPE) is the average of the absolute percentage errors of the predictions. This metric is useful for comparing the performance of datasets with different scales but can be problematic when target values are close to zero. A lower MAPE indicates more accurate predictions in percentage terms, making it easier to understand the relative error size. A higher MAPE indicates larger percentage errors. MAPE is determined by the following Eq. (4).

(4) $$MAPE={\displaystyle \frac{1}{n}\sum _i^n\left|{\displaystyle \frac{Y_i-{\hat{\mathrm{Y}}}_i}{Y_i}}\right|.}$$Median Absolute Error (MdAE) is the median of the absolute errors. It is even more robust to outliers than MAE and provides a stable indicator of the typical magnitude of errors. MdAE is expressed by the following Eq. (5).

(5) $$MdAE=Median\left(\left|Y_i-{\hat{\mathrm{Y}}}_i\right|\right).$$Finally, Root Mean Squared Logarithmic Error (RMSLE) is utilized particularly when the target variable spans a wide range and it is desired that smaller errors are penalized less severely than larger errors. It penalizes over predictions more heavily than under-predictions. RMSLE is computed using the following Eq. (6).

(6) $$RMSLE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n(\log \left(1+Y_i\right)-\mathrm{l}\mathrm{o}\mathrm{g}{\left(1+{\hat{\mathrm{Y}}}_i\right)}^2}}.$$These metrics were employed to quantitatively evaluate the training and test performance of the models, as well as the cross-validation results.

Model interpretability and comparison

To enhance the understanding of complex model predictions, interpretable ML techniques are employed, specifically for lm and rpart models. For lm, the coefficients and the resulting equation are displayed, providing insight into the linear relationships learned by the model. For rpart, the tree structure itself is visualized, illustrating the decision rules and feature splits.

A critical aspect of this study involves comparing the predictive capabilities of the developed ML models against established theoretical or empirical models from the literature, specifically referencing Pascoe’s thesis (Pascoe, 2016). This comparison is performed by generating a dedicated plot and a summary table showing key metrics (MAE, R², Correlation) for both ML predictions and the theoretical model’s predictions on the test set. This direct comparison highlights the potential advantages of data-driven ML approaches over traditional methods.

Reporting and export functionality

The R-Shiny application is designed to generate comprehensive analysis reports in HTML formats. These reports encapsulate all key outputs of the analysis, including data previews, exploratory data analysis plots (histograms, boxplots, correlation matrices), Boruta feature selection results, detailed model performance metrics (training, test, and cross-validation), model training times, actual vs. predicted plots, residual plots, and the comparison with theoretical models. Furthermore, the application allows users to download the cleaned dataset and all generated model metrics in comma-separated values (CSV) format, ensuring data accessibility and reproducibility.

Data preprocessing and outlier management

Upon loading the raw dataset, an initial assessment of missing values was performed. The data showed complete cases for the critical columns after an initial NA cleaning step, indicating that any rows with missing values in the predefined critical columns were removed.

For outlier management, the analysis was performed with the outlier management option enabled, utilizing a Z-score threshold of 3. The chosen method was to “Remove Rows” identified as outliers. The outlier summary indicated that specific numbers of outliers were detected and removed per column. This process ensured that the cleaned dataset (df_clean) used for subsequent modeling contained robust and reliable data points.

Exploratory data analysis for da/dN

Histograms of the feature distributions revealed the general shape and spread of each numeric variable in the cleaned dataset. For da/dN, the distribution of its relevant predictors such as Gmax, R, and Delta_sqrt.G. were visualized. Boxplots provided a visual representation of the data’s central tendency, spread, and the presence of any remaining extreme values after outlier management.

The correlation analysis (Fig. S3) revealed the linear relationships between da/dN and its predictor variables. A strong positive correlation was observed between da/dN and both Gmax (r = 0.56) and Delta_sqrt.G. (r = 0.56). In contrast, da/dN showed a weak negative correlation with R (r = −0.21). Among the predictor variables themselves, Gmax and R exhibited a moderate negative correlation (r = −0.42), as did R and Delta_sqrt.G. (r = −0.47). Gmax and Delta_sqrt.G. demonstrated a perfect positive correlation (r = 1).

Feature selection for da/dN

This Boruta feature importance plot (Fig. S4) illustrates the relative importance of the provided features (G_max, Delta_sqrt.G., and R) for predicting da/dN. All features are labeled as “Confirmed” in green, indicating that the Boruta algorithm identified them all as statistically significant predictors for da/dN. G_max exhibits the highest mean importance, followed by Delta_sqrt.G., and then R. This suggests that G_max is the strongest predictor of da/dN, with Delta_sqrt.G. also being a substantial contributor, and R being significant but less impactful than the other two.

ML model performance for da/dN

A total of 16 different ML regression algorithms were trained and evaluated for the prediction of da/dN. The dataset was split into a training set and a test set (e.g., 80% training, 20% testing), and a 5-fold cross-validation was applied during training. The performance metrics for each model on the training set, test set, and through cross-validation are summarized below.

Training set metrics for da/dN

The performance of each model on the training set, indicating how well the models learned the underlying patterns in the data, is presented in Table S1.

Based on the training set metrics, rf model demonstrated exceptionally strong performance with an outstanding R² value of 0.9999219 and the lowest RMSE of 4.00 × 10⁻⁷. Other high-performing models included knn, gam, brnn, and qrnn, all exhibiting R² values above 0.9997 and similarly low RMSE values around 6.00 × 10⁻⁷. These models show a remarkable ability to capture the complex relationships within the training data. Conversely, rpart, xgbTree, pls, and glmnet generally displayed lower R² values (below 0.92) and higher RMSEs, indicating a less accurate fit to the training data compared to the top-performing models.

Test set metrics for da/dN

The generalization capability of each model was assessed on the unseen test set. These metrics provide a more realistic estimate of model performance on new data. Table S2 summarizes the test set metrics for all models.

Upon evaluation of the test set, rf model consistently demonstrated superior generalization performance, achieving the highest R² value of 0.9999219 and the lowest RMSE of 4.00 × 10⁻⁷. Closely following in performance were knn, gam, brnn, and qrnn, all showing robust R² values exceeding 0.9996 and remarkably low RMSEs around 6.00 × 10⁻⁷. These models exhibited excellent predictive accuracy on unseen data. Conversely, pls, rpart, xgbTree, and glmnet yielded comparatively lower R² values (ranging from 0.827 to 0.915) and notably higher RMSEs, indicating a diminished ability to generalize effectively to new, unobserved data.

Cross-validation metrics for da/dN

Cross-validation metrics provide a robust estimate of model performance by averaging metrics across multiple folds, thus reducing variability. Table S3 displays the cross-validation metrics for the best-tuned parameters of each model.

The cross-validation results largely confirm the superior performance of several models. The rf model demonstrated exceptionally strong generalization capabilities, achieving a very high cross-validated R² of 0.9997595 and a remarkably low RMSE_CV of 6.00 × 10⁻⁷. Similarly, knn, gam, brnn, and qrnn also exhibited excellent and consistent performance, with cross-validated R² values all above 0.99967 and RMSE_CV values at or around 6.00 × 10⁻⁷ to 8.00 × 10⁻⁷. These models show a robust ability to predict da/dN accurately across different data subsets. Conversely, models such as rpart, xgbTree, and glmnet generally yielded lower cross-validated R² values (ranging from 0.882 to 0.919) and higher RMSE_CVs, indicating a less consistent and less accurate predictive performance on unseen data when subjected to cross-validation.

Model training time for da/dN

The time taken for each model to complete its training process is crucial for practical implementation and model selection, particularly for large datasets. Figure S5 illustrates the training times.

As depicted in Fig. S5, there were significant variations in training times across the different models. qrnn (3,463.9 s) and gaussprRadial (459.21 s) were the most computationally intensive, requiring substantially longer training durations. Following them, rf (144.84 s), xgbTree (26.7 s), and also brnn (12.89) exhibited relatively high training times. In contrast, models like lm (0.26 s), pls (0.28 s), lasso (0.35 s), rpart (0.42 s), ridge (0.49 s), glmnet (0.76 s), knn (1.28 s), gam (1.37 s), gbm (1.99 s), svmRadial (3.04 s), and svmLinear (4.94 s) all completed their training within 5 s, suggesting greater efficiency for practical applications where training speed is a critical factor.

Actual vs. predicted values and residual plots for da/dN

Visual inspection of model performance was conducted through scatter plots of actual vs. predicted values and plots of predicted values vs. residuals for the test set.

Actual vs. Predicted Values graph (Fig. 1) displays how well each model’s predictions align with the true experimental values on the test set. A perfect model would show all points lying on the Y = X dashed red line. As observed in Fig. 1, models such as rf, knn, gam, brnn, and qrnn demonstrate excellent predictive accuracy, with their predicted values closely hugging the Y = X dashed red line, indicating a near-perfect alignment with the actual test set values. This suggests minimal bias and high precision for these models. In contrast, models like rpart, xgbTree, glmnet, and pls show a wider dispersion of points around the Y = X line, indicating less accurate predictions and a greater deviation from the true values. The performance of lm, ridge, lasso, svmLinear, svmRadial, and gaussprRadial falls somewhere in between, showing a reasonable, but not always perfect, alignment with the ideal prediction line.

Figure 2 illustrates the distribution of residuals for each model on the test set. For rf, knn, gam, brnn, and qrnn, the residuals are tightly clustered and randomly scattered around the zero line, with no discernible patterns. This indicates that these models produce unbiased predictions and their errors are random, suggesting a robust fit. Conversely, models such as rpart, xgbTree, glmnet, and pls exhibit more pronounced patterns or wider spreads in their residuals, particularly at certain predicted values, which suggests systematic errors or heteroscedasticity in their predictions. The remaining models, including lm, ridge, lasso, svmLinear, svmRadial, and gaussprRadial, show varying degrees of random residual distribution, but generally with wider spreads compared to the top-performing models, indicating less consistent error patterns.

Linear regression model for da/dN

The equation of the lm model, derived from its coefficients, provides a direct mathematical relationship between the target variable (da/dN) and the selected features. The Linear Regression model’s estimated equation for predicting da/dN, designated as Eq. (7), is presented as follows:

(7) $$da/dN=-9.1039e-04+1.3750e-02Delta\mathrm{\_}sqrt.G.+3.0945e-03R-6.8287e-03G\mathrm{\_}max.$$This equation, representing the lm model, reveals the linear relationships between da/dN and its predictors. The positive coefficient for Delta_sqrt.G. (1.3750e−02) indicates that an increase in Delta_sqrt.G. is associated with an increase in da/dN, assuming other variables are held constant. Similarly, R exhibits a small positive association (3.0945e−03) with da/dN. Conversely, Gmax shows a negative coefficient (−6.8287e−03), suggesting that as Gmax increases, da/dN tends to decrease when other predictors remain unchanged. The intercept (−9.1039e−04) represents the predicted da/dN when all predictor variables are zero.

Decision tree model for da/dN

Figure 3 illustrates the rpart model used to predict da/dN based on hierarchical decision rules derived from the variable Gmax. The root node represents the entire dataset (n = 5,401), which is first split based on the threshold Gmax < 0.42. The majority of samples (88%) fall into this category, resulting in a lower predicted da/dN value (31e−6). The remaining 12% of the data, with Gmax ≥ 0.42, is further divided at Gmax < 0.48, leading to intermediate (114e−6) and higher (204e−6) da/dN predictions. This structure indicates that higher Gmax values are associated with increased crack growth rates, demonstrating the model’s capacity to separate observations based on the progressive influence of Gmax.

A key objective of this study was to compare the predictive capabilities of the developed ML models with a theoretical model derived from existing literature, specifically referencing Pascoe’s thesis. This comparison provides valuable context on the added value of ML approaches.

Table S4 provides a quantitative summary of the comparison metrics (MAE, R², and Correlation) for the ML models and the theoretical model. This table allows for a clear assessment of which models best capture the underlying relationships in the data. The table clearly shows that models such as qrnn, rf, and brnn outperform others with nearly perfect scores across all metrics (MAE ≈ 0, R² ≈ 0.9998, Correlation ≈ 0.9999), indicating exceptional predictive accuracy. In contrast, models like rpart and xgbTree exhibit the lowest performance, with significantly lower R² and higher MAE, suggesting they are less effective at capturing the true relationship in the data.

Prediction of Total_Energy

This subsection details the results obtained from applying the ML methodologies for predicting the Total_Energy in adhesive bonds.

Exploratory data analysis for Total_Energy

Histograms (Fig. S6) revealed the distributions of Total_Energy and its potential predictors, including N, a, F, d, C, and R. Boxplots provided insights into the spread and presence of outliers for these features.

The correlation matrix for Total_Energy and its relevant predictors (N, a, F, d, C, R) provided insights into their linear relationships. The correlation matrix (Fig. S7) reveals strong negative correlations between Total_Energy and both a (−0.97) and C (−0.96), suggesting that increases in crack length and compliance are associated with a significant reduction in total energy. Additionally, N also shows a strong negative correlation (−0.80) with Total_Energy. On the other hand, F has a weak positive correlation (0.28), while d shows virtually no linear relationship. These patterns highlight which variables most influence energy dissipation during crack growth.

Feature selection for Total_Energy

The Boruta algorithm was applied to identify important features for predicting Total_Energy. The specified predictors fed into Boruta were N, a, F, d, C, and R. As illustrated in Fig. S8, the Boruta algorithm identified N, a, F, d, C, and R as “Confirmed” important features for predicting Total_Energy, indicating their statistical significance. Among these, N exhibited the highest mean importance, followed by a, C, F, R, and d in descending order.

ML model performance for total energy

The same suite of 16 ML regression algorithms was employed to predict Total_Energy. The models were trained and evaluated using the previously defined training/test split and 5-fold cross-validation.

Training set metrics for Total_Energy

Table S5 summarizes the performance of each model on the training set for Total_Energy prediction. Among all models, the rf model achieved the best results with the lowest RMSE (0.2723), MAE (0.1979), and MAPE (0.0008), as well as the highest R² score (0.9999), indicating exceptional predictive accuracy. In contrast, models like gam and xgbTree showed relatively higher error values, though all models maintained R² values above 0.999, suggesting strong overall fit.

Test set metrics for Total_Energy

Table S6 presents the generalization performance of the models on the unseen test set for Total_Energy prediction. The qrnn model outperformed others with the lowest RMSE (0.454) and MAE (0.359), maintaining a high R² score of 0.9998, indicating strong generalization ability. While the rf model showed top performance during training, its test errors were slightly higher, suggesting potential overfitting. Overall, all models retained excellent R² values above 0.999, confirming their robustness on unseen data.

Cross-validation metrics for Total_Energy

Table S7 displays the robust cross-validation metrics for Total_Energy prediction models. Among them, the qrnn model achieved the best overall performance with the lowest RMSE_CV (0.4451) and highest R²_CV (0.9998), indicating strong consistency and generalization across folds. While brnn and gam also performed well, models like rf and xgbTree showed relatively higher cross-validated errors, suggesting more variability in their predictions across different data subsets.

Model training time for Total_Energy

Figure S9 illustrates the computational time required for each model to train for Total_Energy prediction. As depicted in Fig. S9, there were significant variations in training times across the different models. qrnn (5,132.98 s) and gaussprRadial (481.23 s) were the most computationally intensive, requiring substantially longer training durations. Following them, rf (331.84 s), xgbTree (38.08 s), and also brnn (27.45) exhibited relatively high training times. In contrast, models like lm (0.3 s), pls (0.35 s), lasso (0.37 s), rpart (0.52 s), ridge (0.63 s), svmLinear (0.68 s) glmnet (0.81 s), knn (1.72 s), svmRadial (2.05 s), gbm (2.72 s), and gam (7.83 s) all completed their training within 10 s, suggesting greater efficiency for practical applications where training speed is a critical factor.

Actual vs. predicted values and residual plots for Total_Energy

Visual assessment of the models’ performance for Total_Energy was carried out. As depicted in Fig. 4, which illustrates the correlation between predicted and actual Total_Energy values, all models except rpart (specifically, the qrnn, brnn, gam, and rf models) demonstrate nearly perfect agreement. The data points for these models remarkably cluster around the ideal diagonal line, indicating that their predictions are exceptionally accurate and align almost precisely with the true Total_Energy values. This visual representation unequivocally confirms the outstanding performance of these models, particularly highlighting their near-perfect fit in predicting Total_Energy. This visual representation further confirms the excellent performance of these models which appear to have the tightest clusters, reinforcing the numerical findings from Tables S6 and S7 regarding their superior predictive capabilities.

Furthermore, a detailed examination of the residual plots in Fig. 5, which visualizes the distribution of errors for each model, provides additional evidence of their predictive accuracy for Total_Energy. For all models, the residuals are tightly clustered around zero, displaying a symmetrical distribution with minimal scatter. This indicates that these models produce unbiased predictions with small and randomly distributed errors, which is characteristic of models with excellent predictive power and a very good overall fit. Conversely, the rpart model’s residual plot, if available for Total_Energy, typically exhibits a wider spread of errors and a less centralized distribution, suggesting a comparatively lower accuracy and higher variability in its predictions when compared to the other advanced models.

Linear regression model for total energy

The derived equation for the lm model predicting Total_Energy is presented below, illustrating the linear relationships between features and the target variable in Eq. (8).

(8) $$\begin{array}{rl}Total\mathrm{\_}Energy=& 1.0964e+03+2.6452e+02R-1.1894e+02C-2.3716e-01d\\ & +2.8517e-03F-1.0671e+01a+8.6403e-05N.\end{array}$$

Decision tree model for Total_Energy

Figure 6 illustrates the hierarchical decision rules of the rpart model for Total_Energy prediction.

Table S8 provides a quantitative comparison of the performance metrics (MAE, R², and Correlation) between the ML models and the theoretical model for Total_Energy prediction. While the correlation coefficients for all models, except rpart, are approximately 99.9%, the rpart model yielded a correlation of 89.89%. A similar pattern is observed for the R² values.

Prediction of Gmax

This subsection presents the results obtained from the application of ML methodologies for predicting the Gmax in adhesive bonds.

Exploratory data analysis for Gmax

Histograms (Fig. S10) provided insights into the distributions of Gmax and its potential predictors, specifically a and F. Boxplots visually represented the spread and presence of outliers for these features.

The correlation matrix (Fig. S11) for Gmax and its relevant predictors (a, F) provided insights into their linear relationships. A strong negative correlation was observed between Gmax and a (r = −0.96), and a low positive correlation was observed between Gmax and F (r = 0.28). These patterns highlight which variables most influence Gmax during its determination.

Feature selection for Gmax

The Boruta algorithm was applied to identify important features for predicting Gmax. The specified predictors fed into Boruta were a and F. The Boruta algorithm’s feature importance analysis (Fig. S12) reveals that feature a is a much more significant predictor for G_max compared to feature F. Both features have been “Confirmed” by Boruta, indicating that they are more important than random noise, but a holds a dominant importance.

ML model performance for Gmax

The same suite of 16 ML regression algorithms was employed to predict Gmax. The models were trained and evaluated using the previously defined training/test split and 5-fold cross-validation.

Training set metrics for Gmax

Table S9 summarizes the performance of each model on the training set for Gmax prediction.

Test set metrics for Gmax

Table S10 presents the generalization performance of the models on the unseen test set for Gmax prediction.

Cross-validation metrics for Gmax

Table S11 displays the robust cross-validation metrics for Gmax prediction models.

Model training time for Gmax

Figure S13 illustrates the computational time required for each model to train for Gmax prediction. As depicted in Fig. S13, there were significant variations in training times across the different models. qrnn (4,342.64 s) and gaussprRadial (483.2 s) were the most computationally intensive, requiring substantially longer training durations. Following them, rf (86.07 s), and also xgbTree (34.83 s) exhibited relatively high training times. In contrast, models like pls (0.25 s), lm (0.27 s), lasso (0.31 s), rpart (0.38 s), ridge (0.53 s), glmnet (0.87 s), gam (1.04 s), knn (1.1 s), svmRadial (2.05 s), gbm (1.88 s), svmLinear (6.3 s), and brnn (7.31) all completed their training within 10 s, suggesting greater efficiency for practical applications where training speed is a critical factor.

Actual vs. predicted values and residual plots for Gmax

Visual assessment of the models’ performance for Gmax was carried out. Figure 7 shows the scatter plots of actual vs. predicted Gmax values on the test set, indicating the models’ accuracy. Figure 8 displays the residuals plotted against predicted Gmax values, which should ideally show a random distribution around zero.

Linear regression model for Gmax

The derived equation for the lm model predicting Gmax is presented below, illustrating the linear relationships between features and the target variable in Eq. (9).

(9) $$Gmax=1.8074e+00+2.0913e-06F-1.6816e-02a.$$

Decision tree model for Gmax

Figure 9 illustrates the hierarchical decision rules of the rpart model for Gmax prediction.

General performance of ML models

Overall, the suite of 16 ML regression algorithms demonstrated strong capabilities in predicting the critical fatigue parameters of adhesive bonds. Across all three target variables (da/dN, Total_Energy, and Gmax), many ML models achieved high R² values and low RMSE, MAE, MAPE, MdAE, and RMSLE scores on the unseen test data. This indicates their robust generalization ability and potential for accurate predictions in real-world engineering applications. The consistency in performance metrics across a diverse range of models (from linear methods like Ridge and Lasso to complex ensemble techniques like XGBoost and non-linear approaches like gpr and qrnn) suggests that the underlying relationships within the experimental dataset are well-captured by data-driven approaches.

The variability in performance observed among different ML models highlights the importance of comparative analysis. While some models, particularly complex ensemble methods (e.g., rf, xgbTree, gbm), generally exhibited superior predictive accuracy (lower RMSE, higher R²) due to their ability to capture complex non-linear relationships and interactions, simpler models (e.g., lm) still provided valuable insights, especially regarding feature interpretability. The gaussprRadial and qrnn models, known for their capabilities in handling non-linear data and quantifying uncertainty, showed competitive performance, which is particularly relevant in material science where data scarcity and inherent variability are common challenges. Their inclusion in the model portfolio underscores the comprehensive nature of this study in exploring diverse modeling paradigms. Overall, the suite of 16 ML regression algorithms demonstrated strong capabilities in predicting the critical fatigue parameters of adhesive bonds.

Interpretation of model performance and physical links

The superior performance of non-linear and tree-based ensemble models (e.g., rf, gbm) over simpler linear methods is a key finding of this study. This can be attributed to their inherent ability to capture the complex, non-linear relationships that govern fatigue and fracture mechanics. Physical phenomena such as crack tip plasticity, viscoelastic behavior, and complex energy dissipation mechanisms create interactions that are not easily modeled by linear equations. Furthermore, the feature importance analysis performed by the Boruta algorithm provides a direct and powerful link between the data-driven findings and established physical principles. For instance, the algorithm’s identification of Gmax and R as the most significant predictors for da/dN serves as a data-driven validation of foundational concepts in Linear Elastic Fracture Mechanics (LEFM), where the stress intensity factor range (related to G) and load ratio are known to be primary drivers of fatigue crack growth. The fact that the algorithm, without any prior physical knowledge, independently identified these variables as the most critical predictors validates both the quality of the experimental data and the efficacy of the ML approach. A more detailed physical interpretation can also be derived from the coefficients of the simpler, though less accurate, Linear Regression model for da/dN (Eq. (7)).

• The positive coefficient for Delta_sqrt.G. aligns perfectly with the physical principle that a larger cyclic energy release rate range acts as a stronger driving force, thus accelerating crack propagation.

• Similarly, the positive coefficient for R is consistent with experimental observations where a higher load ratio (often leading to less crack closure) typically results in faster crack growth.

• Interestingly, the model assigned a negative coefficient to Gmax when the other variables are held constant. While seemingly counterintuitive, this is likely an artifact of the multicollinearity within the linear model, especially given the perfect correlation between Gmax and Delta_sqrt.G. (r=1) in the dataset. With the primary driver (the range of G) already in the model, the negative coefficient for the maximum value of G may be capturing complex secondary effects or statistical interdependencies. This highlights both the interpretability and the potential limitations of using linear models for inherently non-linear physical phenomena.

Generalizability, limitations, and synergy with physics-based models

It is important to acknowledge that the specific models trained in this study are tailored to the experimental data from the FM94 epoxy/Aluminum 2024-T3 system. Therefore, their direct application to other adhesive and adherend material systems would require re-training or transfer learning. However, the overarching conclusion that advanced, non-linear ML models are highly effective for this prediction task is likely generalizable, as the fundamental physics of fatigue are complex across most engineering materials.

Furthermore, a significant opportunity for future work lies in creating a powerful synergy between data-driven and physics-based approaches to overcome the common limitation of data scarcity in experimental studies. While the ML models developed in this article rely on existing experimental data, they are inherently “data-hungry”. Physics-based simulations, such as Finite Element Method (FEM) models incorporating constitutive laws from CDM, can be used to generate large, high-fidelity synthetic datasets under a wide range of conditions. This computationally generated data can then be used to augment sparse experimental datasets, allowing for the training of more robust, accurate, and widely generalizable ML models. This hybrid approach would leverage the theoretical strength of mechanics-based models and the pattern-recognition power of machine learning, representing a significant step forward in predictive materials science.

Future perspective: physics-informed machine learning (PIML)

An even more integrated approach for future research is the use of Physics-Informed Machine Learning (PIML). PIML frameworks embed known physical laws (e.g., differential equations governing energy balance or crack growth) directly into the loss function of a neural network. This would compel the model to not only fit the experimental data but also adhere to the fundamental principles of fracture mechanics. Such an approach could lead to models with superior generalization capabilities, especially when extrapolating to conditions outside the training data range, and represents a powerful direction for creating truly hybrid predictive tools.