Interpretable Machine-Learning-Assisted Stochastic Buckling Assessment of Cylindrical Shells with Random Geometric Imperfections

1  Introduction

Thin-walled cylindrical shells are widely used in civil, aerospace, marine, mechanical, and petrochemical engineering because of their high load-carrying efficiency and low structural weight. Under axial compression, however, these structures are highly sensitive to buckling. Although the classical theoretical buckling loads of cylindrical shells have been established in early analytical studies [1–3], experimental buckling loads are often much lower than the corresponding theoretical predictions. This discrepancy is generally attributed to unavoidable uncertainties in geometry, boundary conditions, thickness, material properties, and loading conditions [4,5]. Among these factors, initial geometric imperfections are widely recognized as one of the dominant sources of buckling-load reduction.

Accurate characterization of geometric imperfections is therefore essential for reliable buckling assessment and knockdown-factor evaluation of cylindrical shells [5]. Existing imperfection descriptions can be broadly classified into measured imperfections, deterministic modal imperfections, and stochastic imperfection fields. Measured imperfections obtained using contact or non-contact techniques provide direct information on real shell geometries [4,6–8], but they usually describe specific tested specimens and are not sufficient by themselves for large-scale stochastic assessment. Deterministic imperfection modes, often based on the critical buckling mode or a combination of eigenmodes [9,10], are convenient for imperfection-sensitivity studies but cannot fully describe the spatial variability of geometric imperfections. Random field based descriptions provide a more general framework for representing distributed geometric uncertainty on shell surfaces [11–13].

Random field simulation has been extensively studied using spectral representation [14], Karhunen–Loève expansion [15], and stochastic harmonic functions [16]. Many of these methods were originally developed for planar or simply parameterized domains. For cylindrical shells, the imperfection field is defined on a curved surface with a periodic circumferential direction. Several studies have generated random imperfections on cylindrical shells or panels using separable axial and circumferential correlation descriptions [17–21]. Other studies have considered geodesic-distance-based correlation modeling for random fields on curved surfaces [22], and manifold-based random field simulation methods have also been developed [23–25]. For closed cylindrical shells, however, the periodicity of the circumferential coordinate must be treated explicitly to avoid artificial discontinuities at the seam of the unwrapped surface.

Once stochastic imperfection fields are generated, their influence on buckling resistance must be evaluated over a large number of realizations. Nonlinear finite-element analysis is commonly used for this purpose [21,26–28], but repeated nonlinear Riks analyses become computationally demanding for large stochastic datasets. Koiter-based reduced-order approaches provide an efficient alternative for imperfection-sensitivity analysis; for example, Barbero et al. [29] used Koiter’s method to characterize modal imperfection sensitivity of cylindrical shells. In contrast, the present study focuses on spatial random field imperfections and their morphology-dependent buckling response. Beyond computational cost, response statistics alone do not directly reveal which imperfection morphology characteristics are most relevant to buckling resistance.

Machine-learning and surrogate-modeling techniques have increasingly been introduced into shell buckling analysis to improve response prediction efficiency. Artificial neural networks and other machine-learning models have been applied to buckling-load prediction of cylindrical shells using experimental, parametric, and finite-element data [30–33], while physics-informed approaches incorporate physical constraints into data-driven buckling models [34,35]. Image-driven models have further predicted buckling behavior directly from imperfection patterns [36], and neural-network surrogates have recently been developed for random field based geometric imperfections [37]. More recent studies have further extended machine-learning-assisted shell analysis toward explainable and uncertainty-aware modeling, including explainable learning for cylindrical-shell buckling and failure [38,39] and machine-learning-assisted uncertainty quantification for imperfect cylindrical shells [40].

To clarify the positioning of the present study, Table 1 summarizes representative research lines and highlights the different aspects emphasized in the present work.

The remainder of this paper is organized as follows. Section 2 introduces the geodesic random imperfection model with circumferential periodicity and the corresponding Gaussian random field generation procedure. Section 3 describes the finite-element cylindrical shell model, the stochastic imperfection database, and the nonlinear Riks analysis used to extract the critical load factors. Section 4 defines the candidate imperfection descriptors and presents the feature-diagnosis and surrogate-modeling methodology. Section 5 reports the stochastic buckling response, feature-diagnosis results, sample-level surrogate performance, and group-level statistical-response trends. Section 6 discusses the role of correlation-length-induced morphology, the interpretation of the RMS edge jump descriptor, and the limitations of the proposed framework. Section 7 concludes the paper.

2  Geodesic Random Imperfection Model with Circumferential Periodicity

2.1 Radial Geometric Imperfection of Cylindrical Shell

As shown in Fig. 1, consider a thin cylindrical shell with radius R, length L, and thickness t. A point on the middle surface of the perfect cylindrical shell is parameterized by the circumferential coordinate θ and the axial coordinate z as

x(θ,z)=[Rcos⁡θRsin⁡θz]T,0≤θ<2π,0≤z≤L(1)

Figure 1: Radial imperfection of the cylindrical shell.

The geometric imperfection is introduced as a radial perturbation of the cylindrical middle surface. Let w0(θ,z) denote the scalar imperfection field. The imperfect middle surface is expressed as

x0(θ,z)=x(θ,z)+w0(θ,z)er(θ),(2)

where er(θ)=[cos⁡θsin⁡θ0]T is the outward radial unit vector. A positive value of w0 denotes an outward radial imperfection, whereas a negative value denotes an inward radial imperfection.

2.2 Geodesic Distance with Circumferential Periodicity

For a cylindrical shell, the circumferential direction is periodic. As shown in Fig. 2, by unfolding the cylindrical surface into a rectangle with periodic boundary treatment in the circumferential direction, the shortest surface distance between two points (θi,zi) and (θj,zj) is written as

Dg,ij=(RΔθij)2+(zi−zj)2,(3)

where Δθij=min(|θi−θj|,2π−|θi−θj|) is the shortest angular separation.

Figure 2: Geodesic distance on the unwrapped cylindrical surface with circumferential periodicity.

This distance is referred to as the geodesic distance, representing the intrinsic shortest-path distance along the cylindrical surface. The periodic treatment accounts only for the circumferential topology of the closed cylinder and prevents an artificial discontinuity at the seam of the unwrapped surface. It does not impose a periodic waveform or prescribed geometric shape on individual imperfection realizations.

For a finite element mesh with n nodes, the pairwise geodesic distances form the matrix

Dg=[Dg,ij]i,j=1n.(4)

The matrix is used to construct the covariance matrix of the random imperfection field.

2.3 Geodesic Random Field Generation with Circumferential Periodicity

The imperfection field is modeled as a zero-mean Gaussian random field, w∼𝒩(0,C), where C is the covariance matrix.

Let wi=w0(θi,zi) denote the radial imperfection amplitude at the i-th finite element node. The nodal imperfection vector is

w=[w1w2⋯wn]T,(5)

In this study, the squared-exponential covariance function is adopted:

Cij=σw2exp⁡[−(Dg,ijc)2],(6)

where c is the correlation length and σw is the standard deviation before RMS normalization. A smaller c produces short-correlation imperfection patterns, whereas a larger c produces smoother long-correlation patterns.

The covariance matrix is then decomposed as

C=ΦΛΦT,(7)

where Λ=diag(λ1,λ2,…,λn) contains the eigenvalues and Φ contains the corresponding eigenvectors. In numerical implementation, small negative eigenvalues may appear because a covariance function constructed from geodesic distances is not guaranteed to be positive semidefinite for all discretizations and correlation lengths. A non-negative spectral projection is therefore applied when necessary to obtain a valid covariance representation for random field generation:

Λ+=diag(λ1+,λ2+,…,λn+),(8)

where λk+=max(λk,0) for k=1,2,…,n. Then a realization of the imperfection field is generated by

w=ΦΛ+1/2ξ,(9)

in which ξ∼𝒩(0,I) is a standard normal random vector.

The generated imperfection field is then scaled to a normalized imperfection amplitude as follows

w∗=wRMS∗wRMSw,(10)

where wRMS∗ is a prescribed RMS amplitude, and wRMS=1/n∑i=1nwi2 is the RMS amplitude of a discrete sample. This normalization keeps the global imperfection amplitude fixed, allowing the subsequent buckling analyses to isolate the effect of correlation-length-induced spatial morphology.

Finally, the scaled imperfection amplitudes are applied to the nodal coordinates in the radial direction:

xi∗=x(θi,zi)+wi∗er(θi),i=1,2,…,n.(11)

The prescribed correlation lengths and RMS amplitude define a controlled numerical morphology space rather than a manufacturing specific imperfection model. Application to manufactured shells therefore requires calibration of the covariance structure and imperfection amplitude from measured surface data.

3  Nonlinear Stochastic Buckling Analysis

3.1 Buckling Analysis Procedure of Finite Element Cylindrical Shell Model

The cylindrical shell model considered in this study is based on an experimental cylindrical shell reported in Ref. [20], as shown in Fig. 3. The shell has a radius of R=101.6 mm, a length of L=203.2 mm, and a thickness of t=0.11597 mm. The material is assumed to be linear elastic and isotropic, with Young’s modulus E=104,410 MPa and Poisson’s ratio ν=0.3. The bottom edge of the cylindrical shell was fully constrained. At the top edge, the axial displacement along the global Z-direction was allowed, while the non-axial displacement components were constrained. An axial compressive load was then applied at the top edge through the prescribed loading condition. The shell was discretized using four-node reduced-integration shell elements (S4R in Abaqus) with about 5 mm element size, resulting in 5376 nodes and 5248 shell elements.

Figure 3: Finite-element model of the cylindrical shell.

The same finite element model settings were used for all imperfect shells, with only the nodal coordinates modified according to the prescribed imperfection field. For each imperfect cylindrical shell model, a geometrically nonlinear buckling analysis was performed in Abaqus 2022 using the Static, Riks procedure with NLGEOM=ON. The maximum number of increments was set to 200. The initial arc-length increment, total arc-length scale factor, minimum arc-length increment, and maximum arc-length increment were set to 0.002, 1.0, 1.0×10−6, and 0.02, respectively. The analyses were performed on a workstation equipped with a 13th Gen Intel Core i9-13900 processor and 32 GB of RAM under 64-bit Windows 11, with four CPU cores assigned to each Riks analysis.

The critical load factor, λcr, was extracted from the Load Proportionality Factor (LPF) history as the first significant local maximum along the nonlinear equilibrium path. To suppress spurious early oscillations, the first three frames were excluded, the peak LPF was required to exceed 0.05, and the detected peak had to be followed by a load drop. If no such peak was found, the maximum LPF over the computed path was recorded and the case was flagged for inspection. The same criterion was applied to all imperfection samples.

3.2 Random Imperfection Database

Using the geodesic random field generation procedure described in Section 2.3, a stochastic imperfection database was constructed for the cylindrical shell. In the present parametric study, the correlation length c is the only prescribed random field parameter varied to control the spatial morphology. To characterize the correlation length relative to the shell size, the normalized correlation length is defined as

αc=cR.(12)

Eight normalized correlation lengths, αc={0.05,0.10,0.15,0.20,0.30,0.50,0.75,1.00}, are considered to investigate the influence of spatial correlation. The selected range spans characteristic correlation scales from approximately the baseline element size at αc=0.05 to the shell-radius scale at αc=1.00, providing a controlled transition from short-range rough fields through intermediate spatial scales to relatively smooth long-range fields. These values define the morphology space investigated in this study rather than a range calibrated from measured manufacturing imperfections.

For each prescribed value of αc, the covariance matrix was checked for numerical positive semidefiniteness. The magnitude of the negative spectrum before projection was quantified by r−=∑λk<0|λk|/∑λk>0λk. The negative spectral ratio was zero for αc≤0.50. Small negative eigenvalues appeared only for αc=0.75 and αc=1.00, with r−=4.75×10−9 and 1.20×10−5, respectively. These values indicate that the spectral projection has negligible influence on the generated imperfection fields.

For each value of αc, 200 independent random imperfection samples were generated, resulting in a total of 1600 realizations. For reproducibility, the random field generation used a fixed MATLAB random seed of 2026. All generated imperfection fields were scaled to the same RMS amplitude before being applied to the finite element model. The target amplitude was set as wRMS∗/t=0.5. Then for each random imperfection sample, the nodal coordinates of the imperfect shell were generated by imposing the RMS-normalized radial imperfection vector according to Eq. (11). Thus, each random realization produces one geometrically imperfect shell model, while the mesh topology, material properties, loading pattern, and boundary conditions are kept unchanged.

Fig. 4 shows representative RMS-normalized imperfection fields for selected values of αc. Circumferential continuity is evident across the boundary at θ=0 and 2π, without an artificial seam discontinuity. As αc increases, the characteristic spatial scale increases and the imperfection patterns evolve from densely distributed local fluctuations to broader and smoother spatially coherent regions. Since all fields are normalized to the same RMS amplitude, these differences primarily reflect changes in spatial morphology rather than global imperfection amplitude.

Figure 4: Representative random imperfection fields at different correlation lengths with a common RMS amplitude of wRMS∗=0.5t.

3.3 Mesh Sensitivity and Statistical Adequacy Assessment

The shortest prescribed correlation length, c=5.08 mm, is comparable to the characteristic element size of the baseline mesh, making its numerical response potentially sensitive to spatial discretization. A mesh-sensitivity assessment was therefore conducted to evaluate the robustness of the principal buckling response trends under mesh refinement. In addition to the baseline mesh used for the complete 1600-sample database, fine and extra-fine meshes were constructed by successively halving the approximate element dimensions in both surface directions, as summarized in Table 2.

Five realizations were selected at each of the eight correlation lengths, yielding 40 paired baseline–fine comparisons. Six representative realizations were additionally analyzed using the extra-fine mesh, including two cases each at αc=0.05 and 0.10, and one case each at αc=0.30 and 1.00. For each paired comparison, the same underlying imperfection field with wRMS∗/t=0.5 was mapped to the refined mesh by linear interpolation in the circumferential and axial directions. This procedure allows the influence of mesh refinement to be assessed while minimizing changes in the underlying imperfection morphology.

The statistical adequacy of the 200 realizations adopted at each correlation length was further examined through repeated subsampling. Sample sizes of nsample={25,50,75,100,125,150,175} were considered, and 500 random subsets were drawn at each sample size to evaluate the stability of the estimated mean and standard deviation.

4  Candidate Imperfection Descriptors and Surrogate-Modeling Methodology

4.1 Characterization of Candidate Imperfection Descriptors

The random imperfection field generated in Section 3.2 is represented by the nodal imperfection vector w∗, whose dimension is equal to the number of finite element nodes. Directly using w∗ as the surrogate-model input would therefore result in a high-dimensional representation with limited interpretability. To obtain a compact and physically interpretable representation, each imperfection realization is instead characterized by a set of scalar candidate variables, including the prescribed correlation length and morphology descriptors extracted from the imperfection field.

The candidate descriptors considered in this study are summarized in Table 3. For the m-th random imperfection realization, the candidate descriptor vector is defined as

q(m)=[αc(m),w+(m),w−(m),wa(m),wr(m),S(m),K(m),J(m)]T.(13)

These descriptors cover four aspects of the imperfection field: the prescribed spatial scale through αc, extreme amplitudes through w+, w−, wa, and wr, marginal distributional shape through S and K, and local roughness through J. Since all imperfection realizations are normalized to the same RMS amplitude, the field-derived descriptors are introduced to characterize sample-specific morphology.

Because αc is the prescribed parameter controlling the spatial correlation of the random imperfection field, the descriptors extracted from each realization may vary systematically with αc. Therefore, before using these descriptors in the response analysis, their variation across the generated database is first examined. Fig. 5 shows the variation of the candidate descriptors with the prescribed normalized correlation length. This inspection clarifies how the extracted scalar descriptors are distributed across different correlation-length groups.

Figure 5: Variation of candidate imperfection descriptors with the prescribed normalized correlation length αc.

Because J is calculated from differences between adjacent nodal imperfection amplitudes, its magnitude depends on spatial resolution and is therefore treated as a fixed resolution discrete roughness descriptor. All 1600 samples were evaluated on the same baseline mesh, while 40 paired baseline–fine cases were used to assess its resolution dependence, as shown in Fig. 6. Mesh refinement reduced the group-mean J values to 50.04%–56.70% of their baseline values, consistent with the approximately halved nodal spacing of the fine mesh and the corresponding reduction in adjacent-node imperfection differences. Despite this change in magnitude, the ordering of the eight group means was fully preserved. Thus, mesh refinement changes the numerical scale of J but preserves its relative roughness trend.

Figure 6: Resolution dependence of the RMS edge jump.

The candidate descriptors may also contain redundant information because several of them describe related aspects of the same imperfection field. For example, the amplitude-based descriptors characterize different forms of extreme-amplitude information, whereas other descriptors describe distributional shape or local variation. The correlations among the candidate descriptors are therefore examined in Fig. 7. This correlation matrix provides a preliminary view of descriptor redundancy and supports the interpretation of the subsequent feature-diagnosis and surrogate-modeling results.

Figure 7: Correlation matrix among the candidate imperfection descriptors.

4.2 Learning Tasks and Surrogate Models

After the nonlinear Riks analyses described in Section 3.1, each random imperfection sample is associated with a critical load factor λcr. Together with the candidate descriptor vector defined in Section 4.1, the finite element database is written as

𝒟={(q(m),λcr(m))}m=1Ns,(14)

where Ns is the total number of random imperfection samples, q(m) is the candidate descriptor vector of the m-th imperfection field, and λcr(m) is the corresponding critical load factor obtained from the nonlinear finite element analysis.

Two learning tasks are considered based on this database:

{λ^cr(m)=ℳs(fs(m)),sample-level response prediction,r^k=ℳg(fg,k),group-level statistical trend approximation.(15)

here, fs(m) denotes a low-dimensional input for the m-th imperfection realization, and ℳs predicts the corresponding critical load factor λ^cr(m). Similarly, fg,k denotes a group-level input for the k-th correlation-length group, and ℳg predicts its statistical response r^k. The inputs are constructed from the prescribed generation parameter and/or candidate imperfection descriptors, with different input configurations examined in the subsequent analyses.

For sample-level prediction, five commonly used nonlinear regression algorithms were considered: bagged regression trees (BT), random forest (RF), least-squares boosting (LSBoost), support vector regression (SVR), and Gaussian process regression (GPR). The BT and RF models each used 200 regression trees trained by bootstrap resampling; BT considered all available predictors at each split, whereas RF used random predictor subsampling. LSBoost used 200 regression trees with a learning rate of 0.05. SVR used a Gaussian kernel with a kernel scale of 1, a box constraint of 1, and an epsilon-insensitive margin of 0.01. GPR used a constant basis function and an automatic-relevance-determination squared-exponential kernel with exact fitting and prediction. Predictor standardization was applied to SVR and GPR using training-set statistics only.

For group-level trend approximation, GPR was adopted, with a separate regression model established for each response statistic. The group-level input was constructed from either the prescribed generation parameter or a group summary of a candidate imperfection descriptor.

4.3 Feature Diagnosis

The feature-diagnosis step examines the response relevance of the candidate variables before surrogate validation. A one-way variance decomposition with respect to the αc groups is first used to distinguish between-group and within-group contributions to the total variation of λcr.

The diagnosis then consists of three complementary analyses. First, single-descriptor screening is performed using each component of q(m) individually as the input of a bagged-tree surrogate to assess its standalone predictive information. Second, permutation feature importance is evaluated using the field-derived morphology-descriptor set ℱm={w+,w−,wa,wr,S,K,J}, excluding αc because it is the prescribed random field parameter rather than a descriptor extracted from an individual realization. Third, within-αc Pearson correlations between each field-derived descriptor and λcr are evaluated to examine sample-to-sample response associations within each fixed correlation length group.

Permutation feature importance was evaluated using a bagged-tree surrogate trained with the morphology-descriptor set ℱm. For each descriptor, its values were randomly permuted among the test samples while the other descriptors were kept unchanged, and the resulting increase in RMSE relative to the original prediction was used as the importance measure. A larger RMSE increase indicates greater predictive relevance of the corresponding descriptor.

4.4 Validation Strategy

Two validation protocols are used for sample-level surrogate modeling. The first applies a stratified 80%/20% holdout to the complete 1600-sample database using a fixed MATLAB random seed of 1, with 160 training and 40 test samples from each αc group, giving 1280 training and 320 test samples in total. The same data partition was used for all compared models and feature sets to ensure a consistent comparison. Since both the training and test sets contain samples from the same correlation length groups, this validation evaluates the prediction accuracy for unseen random realizations within the constructed database. The second protocol is leave-one-αc-out validation. In each validation round, all samples corresponding to one value of αc are excluded from training and used as the test set. Thus, 1400 samples from seven groups are used for training and 200 samples from the held-out group for testing. This validation is used to assess the prediction performance for an unseen correlation length group. The sample-level prediction performance is evaluated using the root mean squared error (RMSE), mean absolute error (MAE), and coefficient of determination (R2), respectively. For leave-one-αc-out validation, RMSE and MAE are averaged over the eight held-out groups.

Leave-one-group-out validation is used for the group-level statistical response trend approximation, in which one correlation length group is excluded from training and then predicted by the group-level trend approximation model. In each round, seven groups are used for fitting and the remaining group for validation, with all eight groups serving once as the held-out group. The predicted and Riks-derived values from the eight held-out groups are then compared using R2 for each response statistic.

The overall analysis workflow, including feature diagnosis, sample-level prediction, group-level trend approximation, and the corresponding validation procedures, is summarized in Fig. 8.

Figure 8: Workflow of feature diagnosis and surrogate-assisted stochastic buckling assessment.

5  Results

5.1 Correlation-Length-Dependent Buckling Response and Robustness

Fig. 9 shows the distribution of the critical load factors for all random imperfection samples in each αc group, with the corresponding group mean and standard deviation superimposed. Since all samples have the same prescribed RMS imperfection amplitude, the observed differences in λcr mainly reflect the influence of the spatial morphology induced by the correlation length, rather than variations in the global imperfection amplitude.

Figure 9: Distribution of λcr vs. αc.

As shown in Fig. 9a, a clear but non-monotonic dependence of λcr on αc is observed. The mean critical load factor decreases from αc=0.05 to a minimum at αc=0.10, and then increases markedly with increasing correlation length. For αc≥0.50, the mean response reaches a relatively high and nearly stable level. Thus, the lowest buckling resistance does not occur at the shortest correlation length, indicating that the response is governed by more than the magnitude of short-scale local fluctuations. The scatter of λcr also varies with αc. The dispersion is relatively small at the shortest and longest correlation lengths and becomes more pronounced in the intermediate range, reaching its maximum around αc=0.30.

The response quantiles in Fig. 9b show a consistent overall trend, remaining low at short correlation lengths, increasing rapidly over the intermediate range, and approaching a high and relatively stable level for αc≥0.50. The reduced separation between the lower and upper quantiles at large αc further indicates decreased response scatter for smoother long-correlation imperfections. Taken together, these results indicate that the buckling response depends on the spatial organization and characteristic scale of the imperfection field rather than on local fluctuation magnitude alone.

The mesh-sensitivity results are presented in Fig. 10. Mesh refinement changed the absolute values of λcr, particularly in the short- and intermediate-correlation ranges, but retained the principal non-monotonic dependence and the low-response regime around αc=0.10. For the selected three-level cases, the baseline–fine differences were approximately 5%–25%, whereas the fine–extra-fine differences remained below approximately 4.5%, indicating substantially reduced sensitivity under further refinement. These results support the robustness of the principal response trend, while the complete baseline-mesh database is used as mesh-consistent comparative data rather than as a mesh-independent reference.

Figure 10: Mesh-sensitivity assessment of λcr.

The statistical adequacy of the database is assessed in Fig. 11 using the repeated subsampling procedure described in Section 3.3. At n=175, the 95th-percentile relative errors are 0.548% for the mean and 8.918% for the standard deviation. The substantially greater stability of the mean supports the use of 200 realizations per group for comparing the principal response levels and morphology dependent trends, while the larger uncertainty in the standard deviation suggests that small differences in response dispersion should be interpreted cautiously. The zero error at n=200 results from using the full 200 sample estimates as the reference and does not represent independent evidence of Monte Carlo convergence.

Figure 11: Subsample-based stability of the mean and standard deviation of λcr relative to the corresponding 200-sample estimates.

5.2 Machine-learning Feature Diagnosis

The results in Section 5.1 show a strong dependence of the buckling response on the prescribed normalized correlation length. A one-way variance decomposition further shows that 98.67% of the total variance of λcr is associated with between-αc variation, whereas only 1.33% arises from within-group sample-to-sample variation. Thus, the dominant response variation in the present database is associated with the morphology change induced by the prescribed correlation length, while the residual within-group variation is comparatively small.

The response relevance of the candidate descriptors is evaluated using the procedures defined in Section 4.3. Fig. 12a shows that αc and J yield the lowest RMSE values in the single-descriptor screening, while Fig. 12b identifies J as the dominant field-derived descriptor in permutation importance. Together, these results indicate that J is the most informative field-derived descriptor among the candidates considered.

Figure 12: Predictive relevance of the candidate imperfection descriptors.

While Fig. 12 identifies J as the most informative field-derived descriptor for the overall response trend, Fig. 13 shows that the within-αc correlations between the morphology descriptors and λcr are generally weak or inconsistent. Although J exhibits relatively stronger negative correlations in several groups, particularly in the intermediate αc range, these associations remain moderate. This indicates that the scalar descriptors primarily represent the dominant cross-group morphology trend, whereas the residual within-group response variation depends on spatial information not fully captured by the present descriptors.

Figure 13: Within-αc correlation coefficients between morphology descriptors and λcr.

5.3 Sample-Level Surrogate Modeling of Critical Load Factors

Based on the feature-diagnosis results in Section 5.2, three compact input configurations are compared for sample-level prediction of λcr: the prescribed correlation length αc, the field-derived RMS edge jump J, and their combination (αc,J). Their predictive performance is benchmarked across the five regression algorithms described in Section 4.2 under the validation protocols defined in Section 4.4.

As summarized in Table 4, all three input configurations achieve similarly high accuracy under the stratified holdout, whereas clearer differences emerge under leave-one-αc-out validation. The J-based models show relatively consistent cross-group performance across the five algorithms, and replacing αc with J reduces the mean cross-group error for BT, RF, LSBoost, and GPR, while SVR performs better with αc. Combining αc and J does not consistently improve cross-group prediction, indicating that the additional prescribed parameter does not necessarily provide complementary information for extrapolation to an unseen correlation-length group.

Among the J-only models, GPR-J gives the lowest random-holdout error and cross-group performance comparable to the other benchmark algorithms, and is therefore used as a representative smooth morphology-based surrogate. Fig. 14 shows close agreement between its predictions and the Riks results under the stratified holdout, while the remaining scatter reflects spatial information in the full imperfection field that is not represented by J alone. Fig. 15 further shows that the leave-one-αc-out errors vary among the held-out groups, but the overall cross-group information carried by J is consistently observed across the benchmark algorithms.

Figure 14: Predicted vs. Riks λcr for the GPR-J surrogate.

Figure 15: Leave-one-αc-out RMSE of the J-based sample-level surrogate models.

Table 5 further examines the variation among held-out groups for the three GPR input configurations. The J-only model has a lower mean RMSE than the αc-only model, although the improvement is not uniform across the correlation-length range. The αc-only model performs better for αc=0.05, 0.10, and 0.15, whereas J generally gives lower errors for αc≥0.20. In particular, the large J-only error at αc=0.05 indicates that the shortest-correlation case is difficult to extrapolate using a single scalar morphology descriptor.

These group-wise results also show that neither J nor any single input configuration is uniformly superior for every held-out group. Accordingly, GPR-J is used here as a representative morphology-based surrogate rather than as a universally optimal predictor.

Table 6 compares the computational cost of the finite element analyses and the representative GPR-J surrogate. Once the finite element database is available, surrogate training and prediction require negligible time relative to a nonlinear Riks analysis, while evaluation of J adds little additional cost. The surrogate therefore enables efficient repeated prediction within the calibrated structural configuration and morphology range, although its construction still relies on the underlying finite element database.

5.4 Group-Level Statistical-Response Trend Representation

The sample-level analysis is further complemented by a group-level representation of the mean, standard deviation, and selected quantiles of λcr. The prescribed normalized correlation length αc and the group-averaged RMS edge jump J¯ are considered separately as low-dimensional coordinates for these statistical trends. Because only eight correlation-length groups are available, the analysis is intended to summarize trends within the investigated database rather than to establish a general extrapolative statistical surrogate.

Fig. 16a,b compare the Riks-derived statistics with the corresponding group-level trend representations. Both αc and J¯ reproduce the principal variation of the mean and quantile responses across the investigated groups. The J¯-based representation therefore provides a field-derived coordinate for the dominant morphology-dependent trend, although deviations remain at the largest values of J¯. This indicates that J¯ captures the principal group-level morphology information without replacing the full random field description.

Figure 16: Group-level statistics of λcr.

Table 7 shows that the mean and quantile trends are represented more accurately than the standard deviation. The J¯-based representation performs better for the mean and median, whereas αc performs better for the lower and upper quantiles. The relatively low R2 values for the standard deviation indicate that its variation is less readily represented by either low-dimensional input. Thus, the group-level analysis is best viewed as a compact representation of the central and quantile response trends, whereas the sample-level surrogates are used for prediction of individual imperfection realizations.

6  Discussion

6.1 Mechanical Interpretation of the Unfavorable Morphology Scale

The non-monotonic response provides insight into the role of imperfection morphology beyond local roughness alone. The lowest mean critical load factor occurs at αc=0.10, whereas the maximum RMS edge jump occurs at αc=0.05. Under the common RMS normalization, this difference suggests that severe buckling is governed not simply by the magnitude of short-scale fluctuations but also by their spatial organization. At αc=0.05, the imperfection field is dominated by highly fragmented short-scale fluctuations, whereas αc=0.10 retains appreciable local variation over more spatially coherent regions that may interact more effectively with localized nonlinear buckling deformation. As the correlation length increases further, the field becomes progressively smoother and the local variation decreases. The minimum near αc=0.10 is therefore interpreted as a configuration-dependent unfavorable spatial scale between fragmented short-wave and smooth long-wave morphologies, rather than as a universal critical correlation length.

6.2 Interpretation and Practical Role of the Morphology Descriptors

The RMS edge jump J and the prescribed correlation length αc play fundamentally different roles. The former is extracted directly from an imperfection realization and quantifies local roughness, whereas the latter controls the correlation structure used to generate the random field. The results indicate that J provides a useful field-derived representation of the dominant cross-group morphology variation, but it cannot encode the complete spatial structure of an individual imperfection field. Differences in imperfection location, orientation, concentration, and spatial arrangement can therefore contribute to the residual within-group response variation.

The limited predictive relevance of the amplitude-based descriptors w+, w−, wa, and wr and the distributional descriptors S and K is closely related to the RMS-normalized database. These descriptors characterize extreme amplitudes, amplitude ranges, or the marginal distribution of nodal values, but they do not directly quantify the spatial rate of variation of the imperfection field. Their predictive relevance may increase when the RMS amplitude is allowed to vary or when non-Gaussian, localized, or systematically biased imperfection patterns are considered.

The practical role of J is therefore as a compact morphology coordinate for comparative screening and surrogate modeling within a calibrated structural configuration and spatial resolution. For measured imperfection fields, J can be evaluated at a specified resolution and compared with a configuration-specific calibrated descriptor–response relation for rapid morphology screening and buckling-response assessment. It should not be interpreted as a universal manufacturing acceptance threshold, because both its numerical scale and its relation to buckling resistance depend on the structural configuration and the resolution at which the imperfection field is represented.

6.3 Applicability and Limitations

The present results are obtained for a cylindrical shell with L/R=2.0 and R/t≈876.3 under axial compression using RMS-normalized Gaussian random imperfection fields. The reported response trends and trained surrogates are therefore specific to the adopted geometry, boundary conditions, loading condition, material model, imperfection amplitude, and stochastic imperfection representation. The proposed workflow can be extended to other configurations by constructing the corresponding finite element database and recalibrating the descriptor–response relation. Extensions to post-buckling or collapse behavior would similarly require response labels appropriate to those stages.

Because both λcr and the discrete descriptor J depend on spatial discretization, the reported surrogate performance is conditional on the baseline resolution. The refined-mesh analyses support the robustness of the principal response trend and relative roughness ordering, but the numerical scale of J and the associated surrogate relation should be recalibrated when the mesh or measurement resolution changes.

The present database represents a controlled stochastic morphology study rather than a manufacturing specific probabilistic model, because the random field parameters were not calibrated from measured imperfection surveys. Application to manufactured shells would require measured imperfection fields to establish the relevant amplitude and spatial correlation characteristics and to evaluate the morphology descriptors at a defined resolution. Such data would also be needed to assess the transferability of J and recalibrate the morphology buckling response relation for the structural configuration of interest.

7  Conclusions

This study developed a machine-learning-assisted framework for stochastic buckling assessment of cylindrical shells with random geometric imperfections. Geometric random fields with different normalized correlation lengths were generated, scaled to the same RMS imperfection amplitude, and introduced into nonlinear Riks analyses. Based on feature diagnosis, sample-level surrogate modeling, and statistical trend assessment, the following conclusions can be drawn.

•   Under the common RMS imperfection amplitude, the dominant variation of the buckling response in the present database is associated with correlation-length-induced changes in imperfection morphology, with between-αc variation accounting for most of the total response variance.

•   The influence of correlation length on buckling resistance is non-monotonic. The lowest mean critical load factor does not occur at the shortest correlation length, and the largest response scatter appears in the intermediate range. Mesh refinement preserves the principal non-monotonic trend and the low-resistance regime, although the absolute critical load factors remain mesh sensitive.

•   The RMS edge jump J is the most informative scalar descriptor among the considered imperfection features and effectively represents the dominant cross-group morphology variation. Its numerical scale is resolution dependent, although the relative roughness ordering is preserved under mesh refinement.

•   Compact inputs based on αc and/or J accurately reproduce the dominant sample-level response variation under random holdout validation. The J-based models also retain useful cross-group predictive information across the benchmark algorithms, while the group-level models provide compact representations of the main statistical trends.

The present conclusions are restricted to the investigated shell configuration, loading condition, RMS imperfection amplitude, stochastic imperfection model, and spatial resolution, and the resulting surrogates should therefore be regarded as database-dependent response approximations rather than universal buckling predictors. Future work should assess the framework against measured imperfection fields, calibrate the morphology–response relation for manufactured shells, and extend the analysis to broader structural configurations, imperfection amplitudes, and more complex non-Gaussian or localized imperfection patterns.

Acknowledgement: None.

Funding Statement: This research was funded by National Natural Science Foundation of China (grant numbers 52508233, 52578606), the Zhejiang Provincial Natural Science Foundation (grant number LQN26E080052), the visiting scholar funding of the State Key Laboratory of Disaster Reduction in Civil Engineering (grant number SLDRCE25-F-10), and the Alexander von Humboldt Foundation of Germany.

Author Contributions: The authors confirm contribution to the paper as follows: conceptualization, methodology, software, validation, formal analysis, investigation, writing—original draft preparation, funding acquisition, Yan-Ping Liang; writing—review and editing, supervision, funding acquisition, Zhiqiang Wan. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Lorenz R. Achsensymmetrische verzerrungen in dünnwandigen hohlzylindern. Z Ver Dtsch Ing. 1908;52(431706–13. (In German). [Google Scholar]

2. Timoshenko SP. Einige stabilitätsprobleme der elastizitätstheorie. Zeitschrift Für Mathematik Und Physik. 1910;58(4):337–85. (In German). doi:10.1002/zamm.19220020408. [Google Scholar] [CrossRef]

3. Wilson WM, Newmark NM. The strength of thin cylindrical shells as columns. In: Selected papers by Nathan M. Newmark: civil engineering classics. Reston, VA, USA: ASCE; 1933. p. 1–42. [Google Scholar]

4. Arbocz J, Babcock CD Jr. Experimental investigation of the effect of general imperfections on the buckling of cylindrical shells. Washington, DC, USA: NASA; 1968. [Google Scholar]

5. NASA/SP-8007-2020/REV 2. Buckling of thin-walled circular cylinders. Hampton, VA, USA: NASA Langley Research Center; 2020. [Google Scholar]

6. Sebek R. Imperfection surveys and data reduction of ARIANE interstages I/II and II/III [Ir thesis]. Delft, The Netherlands: Delft University of Technology; 1981. [Google Scholar]

7. Vamsi Krishna G, Narayanamurthy V, Viswanath C. Buckling behaviour of FRP strengthened cylindrical metallic shells with cut-outs. Compos Struct. 2022;300(4):116176. doi:10.1016/j.compstruct.2022.116176. [Google Scholar] [CrossRef]

8. Castro SG, Almeida JHS Jr, St-Pierre L, Wang Z. Measuring geometric imperfections of variable–angle filament–wound cylinders with a simple digital image correlation setup. Compos Struct. 2021;276(8):114497. doi:10.1016/j.compstruct.2021.114497. [Google Scholar] [CrossRef]

9. Koiter W. The effect of axisymmetric imperfections on the buckling of cylindrical shells under axial compression. Mechanics. 1963:265–79. doi:10.1090/qam/99838. [Google Scholar] [CrossRef]

10. Catellani G, Pellicano F, Dall’Asta D, Amabili M. Parametric instability of a circular cylindrical shell with geometric imperfections. Comput Struct. 2004;82(31–32):2635–45. doi:10.1016/j.compstruc.2004.07.006. [Google Scholar] [CrossRef]

11. Lauterbach S, Fina M, Wagner W. Influence of stochastic geometric imperfections on the load-carrying behaviour of thin-walled structures using constrained random fields. Comput Mech. 2018;62(5):1107–25. doi:10.1007/s00466-018-1554-0. [Google Scholar] [CrossRef]

12. Chen G, Zhang H, Rasmussen KJ, Fan F. Modeling geometric imperfections for reticulated shell structures using random field theory. Eng Struct. 2016;126(46):481–9. doi:10.1016/j.engstruct.2016.08.008. [Google Scholar] [CrossRef]

13. Fina M, Weber P, Wagner W. Polymorphic uncertainty modeling for the simulation of geometric imperfections in probabilistic design of cylindrical shells. Struct Saf. 2020;82(6):101894. doi:10.1016/j.strusafe.2019.101894. [Google Scholar] [CrossRef]

14. Shinozuka M, Deodatis G. Simulation of multi-dimensional Gaussian stochastic fields by spectral representation. Appl Mech Rev. 1996;49(1):29–53. doi:10.1115/1.3101883. [Google Scholar] [CrossRef]

15. Phoon K, Huang S, Quek S. Simulation of second-order processes using Karhunen–Loeve expansion. Comput Struct. 2002;80(12):1049–60. doi:10.1016/s0045-7949(02)00064-0. [Google Scholar] [CrossRef]

16. Chen J, He J, Ren X, Li J. Stochastic harmonic function representation of random fields for material properties of structures. J Eng Mech. 2018;144(7):04018049. doi:10.1061/(asce)em.1943-7889.0001469. [Google Scholar] [CrossRef]

17. Schenk C, Schuëller G. Buckling analysis of cylindrical shells with random geometric imperfections. Int J Non Linear Mech. 2003;38(7):1119–32. doi:10.1016/s0020-7462(02)00057-4. [Google Scholar] [CrossRef]

18. Papadopoulos V, Papadrakakis M. The effect of material and thickness variability on the buckling load of shells with random initial imperfections. Comput Methods Appl Mech Eng. 2005;194(12–16):1405–26. doi:10.1016/j.cma.2004.01.043. [Google Scholar] [CrossRef]

19. Schenk C, Schuëller G. Buckling analysis of cylindrical shells with cutouts including random boundary and geometric imperfections. Comput Methods Appl Mech Eng. 2007;196(35–36):3424–34. doi:10.1016/j.cma.2007.03.014. [Google Scholar] [CrossRef]

20. Papadopoulos V, Stefanou G, Papadrakakis M. Buckling analysis of imperfect shells with stochastic non-Gaussian material and thickness properties. Int J Solids Struct. 2009;46(14–15):2800–8. doi:10.1016/j.ijsolstr.2009.03.006. [Google Scholar] [CrossRef]

21. Majumder R, Chakraborty S, Mishra SK. Reliability analysis and design of randomly imperfect thin cylindrical shells against post-critical drops. Thin Walled Struct. 2023;185(11):110576. doi:10.1016/j.tws.2023.110576. [Google Scholar] [CrossRef]

22. Scarth C, Adhikari S, Cabral PH, Silva GHC, Prado AP. Random field simulation over curved surfaces: applications to computational structural mechanics. Comput Methods Appl Mech Eng. 2019;345:283–301. [Google Scholar]

23. Feng DC, Liang YP, Ren X, Li J. Random fields representation over manifolds via isometric feature mapping-based dimension reduction. Comput Aided Civ Infrastruct Eng. 2021;37(5):593–611. doi:10.1111/mice.12752. [Google Scholar] [CrossRef]

24. Liang YP, Ren X, Feng DC. Efficient stochastic finite element analysis of irregular wall structures with inelastic random field properties over manifold. Comput Mech. 2022;69(1):95–111. doi:10.1007/s00466-021-02084-4. [Google Scholar] [CrossRef]

25. Liang YP, Feng DC, Ren X, Li J. Three-stage non-Gaussian homogeneous random field representation over manifolds. Comput Aided Civ Infrastruct Eng. 2023;38(11):1462–82. doi:10.1111/mice.12959. [Google Scholar] [CrossRef]

26. Zhang D, Chen Z, Li Y, Jiao P, Ma H, Ge P, et al. Lower-bound axial buckling load prediction for isotropic cylindrical shells using probabilistic random perturbation load approach. Thin Walled Struct. 2020;155(2):106925. doi:10.1016/j.tws.2020.106925. [Google Scholar] [CrossRef]

27. Mahidan F, Ifayefunmi O. The imperfection sensitivity of axially compressed steel conical shells-lower bound curve. Thin Walled Struct. 2021;159(2):107323. doi:10.1016/j.tws.2020.107323. [Google Scholar] [CrossRef]

28. Wagner HNR, Hühne C, Elishakoff I. Probabilistic and deterministic lower-bound design benchmarks for cylindrical shells under axial compression. Thin Walled Struct. 2020;146(79):106451. doi:10.1016/j.tws.2019.106451. [Google Scholar] [CrossRef]

29. Barbero EJ, Madeo A, Zagari G, Zinno R, Zucco G. Imperfection sensitivity analysis of composite cylindrical shells using Koiter’s method. Int J Comput Methods Eng Sci Mech. 2017;18(1):105–11. doi:10.1080/15502287.2016.1276359. [Google Scholar] [CrossRef]

30. Tahir ZR, Mandal P. Artificial neural network prediction of buckling load of thin cylindrical shells under axial compression. Eng Struct. 2017;152(4):843–55. doi:10.1016/j.engstruct.2017.09.016. [Google Scholar] [CrossRef]

31. Tahir ZuR, Mandal P, Adil MT, Naz F. Application of artificial neural network to predict buckling load of thin cylindrical shells under axial compression. Eng Struct. 2021;248(11):113221. doi:10.1016/j.engstruct.2021.113221. [Google Scholar] [CrossRef]

32. Lin X, Jiao P, Xu H, Li X, Chen Z. A machine learning-driven prediction of lower-bound buckling design load for cylindrical shells under localized axial compression. Thin Walled Struct. 2025;209:112960. doi:10.1016/j.tws.2025.112960. [Google Scholar] [CrossRef]

33. Lee HG, Sohn JM. A comparative analysis of buckling pressure prediction in composite cylindrical shells under external loads using machine learning. J Mar Sci Eng. 2024;12(12):2301. doi:10.3390/jmse12122301. [Google Scholar] [CrossRef]

34. Tao F, Liu X, Du H, Yu W. Physics-informed artificial neural network approach for axial compression buckling analysis of thin-walled cylinder. AIAA J. 2020;58(6):2737–47. doi:10.2514/1.j058765. [Google Scholar] [CrossRef]

35. Liu F, Chen H, Yang J, Wang X. Application of physics-informed machine learning methods in buckling design of axially compressed cylindrical shells. Thin Walled Struct. 2024;200(5):111963. doi:10.1016/j.tws.2024.111963. [Google Scholar] [CrossRef]

36. Hao P, Duan Y, Liu D, Yang H, Liu D, Wang B. Image-driven intelligent prediction of buckling behavior for geometrically imperfect cylindrical shells. AIAA J. 2023;61(5):2266–80. doi:10.2514/1.j062470. [Google Scholar] [CrossRef]

37. Schweizer M, Fina M, Wagner W, Kasic S, Freitag S. Artificial neural networks for random fields to predict the buckling load of geometrically imperfect structures. Comput Mech. 2025;76(1):181–204. doi:10.1007/s00466-024-02595-w. [Google Scholar] [CrossRef]

38. Wafi MG, Palar PS, Robani MD, Jusuf A, Zuhal LR, Morlier J. Revisiting cylindrical buckling under axial compression using explainable machine learning. In: Proceedings of the AIAA SCITECH, 2024 Forum; 2024 Jan 8–12; Orlando, FL, USA. [Google Scholar]

39. Zhan M, Wang H, Li Y, Wang L. Explainable machine-learning-assisted hydrostatic failure analysis of moderately thick composite cylindrical shells with ovality. Ocean Eng. 2026;357:125597. doi:10.1016/j.oceaneng.2026.125597. [Google Scholar] [CrossRef]

40. Chen H, Chen G, Yang D, Fu ZJ. Neural network-based DPIM for uncertainty quantification of imperfect cylindrical stiffened shells with multiple random parameters. Eng Anal Bound Elem. 2024;166(12):105795. doi:10.1016/j.enganabound.2024.105795. [Google Scholar] [CrossRef]