Buckling Analysis of a Thin-Walled Structure Using Finite Element Method and Design of Experiments

Abstract

In structural engineering, thin-walled structures play an important role in the design of the lightweight structural model. It carries different loading conditions when it exists in any model, and it is designed with thin plates or thin shells. Penetrating thin-walled structures with different kinds of holes can decrease their weight and facilitate repair and maintenance operations, such as those carried out for the wing of an airplane. In such applications, cutouts are often employed as part of the design of composite plates. Therefore, this paper attempted to design and analyse the thin-walled composite structure with a C-cross-section shape. To model and analyse the structures, a finite element method was utilized using the ABAQUS commercial tool, and the results of critical buckling load for different laminate types were obtained. Composite materials and structures have different parameters that can vary the results of analysis; therefore, to optimize the current mode a design of experiments method is used via MINITAB 20 and Design-Expert 13 tools. The selected parameters for this work were the opening ratio, spacing ratio, and shape of the hole for the output response as a critical buckling load was carried out. Based on the current results of simulation and optimization, it was found that the parameters of composite materials and structures will impact the output response, and the current study investigated the optimum parameters for the best possible outcome of the structural analysis.

1. Introduction

Thin-walled structures are applicable for different engineering applications, such as in aerospace, automotive and construction industries, and they have played an important role in reducing the total weight of structures. One of the systems with the most efficient use of materials to prevent buckling is the thin-walled member because it is made up of several thin-walled parts. The thin-walled portion can easily be formed into various shapes with a high shape factor while using less material. However, the thin-walled member has significant disadvantages associated with its formed plates, such as local buckling. The composite plates of the member usually buckle before collapse when a thin-walled column is subjected to compressive loading [1]. Various shapes of the thin-walled structure can be used to suit its application in the mentioned industries. For example, some studies have shown T-shaped thin-walled structures considering various geometries, which are used in aircraft as ribs [2,3,4,5]. Furthermore, Rozylo and Debski [6] experimented on the Z shape of the thin-walled composite structure and showed the enhancement in improving the structural performance of the whole model.

Generally, buckling is a typical phenomenon in thin-walled structures. It is a term that refers to a component’s lack of stability caused by the lateral deflection on a member when an axial force is applied. The weakness of the column causes it to bend. This mode of failure is rapid and, therefore, risky. The length, strength, and other variables of a column determine whether or not it will buckle. Elastic buckling can occur in long columns relative to their thickness [7]. It can also occur when an applied compressive load exceeds the critical allowable load of the thin-walled structure. Pre-buckling, critical buckling, and post-buckling are the three major states that describe a standard thin-walled composite structure. Buckling and carrying capacity determine the proper loading process of the structure. The structure can still function even after buckling because the post-critical equilibrium trend is stable, and the increasing compressive load increases the wall deflection [8,9]. Figure 1 below shows the state of each operation for a thin-walled structure.

Figure 1. Graph of load applied versus deflection [9]. Reprinted under the Creative Commons (CC) License (CC BY 4.0).

Figure 1. Graph of load applied versus deflection [9]. Reprinted under the Creative Commons (CC) License (CC BY 4.0).

Damaged/holed thin-walled structures have become more prominent in recent investigations, and the reason for failure of the whole structure has been found. Composite materials are widely used in the aerospace and automotive industries due to their advantages, such as their lightweight nature and high strength ratio as compared to other types of materials. In recent investigations, composite materials have been utilized in a wide range of applications in response to the growing need for lightweight and efficient structures. Therefore, it is necessary to investigate the influence of the various sizes, shapes, and distances between cutouts of thin-walled composite structures on buckling and other behaviours. This study continued from the study conducted by Khazaal et al. [10], in which the authors used aluminium alloys as the subject of the experiment. Three types of shapes, circular, hexagon, and rectangular, were chosen to optimize the results. They also studied the parameter that contributes most to the buckling behaviour of the aluminium alloy thin-walled structure.

The FE method was used to explore the stability properties of folded structures and analysed the buckling behaviour under a compression load, shear load, and combinations of both loads with varied boundary conditions. As a result, critical buckling loads for open- and closed-box structures with varying loading conditions were investigated [11]. It has been stated that the newly proposed transverse extension modes have a practical advantage in that the increased artificial stiffness that results from constrained assessments of pure global and distortional deformations can be offset by taking into account the appropriate new transverse extension modes in thin-walled structures [12]. For axial–flexural buckling, post-buckling, and geometrically nonlinear studies of thin-walled beams, an effective one-dimensional FE method using a novel structural concept known as equivalent layered composite cross-sectional modelling was proposed, in which nonclassical factors such as transverse shear and normal flexibilities are included in the formulation [13]. A brand-new beam-type FE method has been developed that may be used for web distortion and elastic buckling analysis of thin-walled beams. Under the premise that there are no pre-buckling deformations, the buckling formulation was created based on nonlinear equilibrium equations [14].

The load-bearing capability and involvement of stress components in the failure analysis of top-hat-shaped composite columns exposed to uniform compression and utilized a digital image correlation (DIC) technique to see the full-field displacements and strains, and failure tests were conducted across the whole load range [15]. For the geometrically nonlinear analysis of thin-walled circular pipes, the formulation of generalized beam theory was used, and the complete geometrically nonlinear analysis was added to the existing generalized beam theory analysis of circular pipes, which is currently confined to buckling analysis [16]. A new technique called the finite strip–Riccati transfer matrix method was established for buckling analysis of thin-walled components with a tree-branched cross-section. The procedure combines the semi-analytical finite strip approach with the Riccati transfer matrix method for a tree multi-body system [17]. To calculate the changes in beam stiffness due to the decrease in thickness in the pre- and post-buckling stages, a nonlinear FE model of a thin-walled beam with variable thickness was introduced. The gauge sensitivity of the beam was then calculated using these stiffness values, and it changed as the beam thickness changed [18]. Residual stresses in thin-walled structures manufactured by the directed energy deposition method were determined; additionally, in situ measurements, fast thermo-mechanical simulation, and buckling were investigated through the experimental and numerical approach [19]. Buckling analysis of thin-walled metal liners of cylindrical composite overwrapped pressure vessels with depressions after autofrettage processing was investigated through the FE approach [20]. Similarly, large deflection and post-buckling of thin-walled structures by finite elements with node-dependent kinematics [21] were investigated.

Sudhirsastry et al. [22] stated that composite materials have been widely used in the automotive and aircraft industries for plates and shells due to their high strength-to-weight ratio and stiffness-to-weight ratio. Due to the widespread usage of composite structures in the form of relatively thin plates, the load-carrying capacity of composite plates against buckling is crucial. Additionally, cutouts are frequently used in composite plates as part of the design of the structures. As a result, it is critical to have a firm grasp of the buckling properties of composite plates with cutouts. Erkliǧ and Yeter [23] found that for composites with a thin-walled, fibre orientation, the angle of the cutouts, the plate size, and lastly the size of holes result in an effect on the buckling behaviour. Nevertheless, the parameter that contributes most towards the behaviour of the buckling of the holed composite plate is not identified. This, together with the work conducted by Khazaal et al. [10], has ignited the motivation to carry out the present study to contribute to the knowledge on thin-walled structures with cutouts specifically for composite material.

The main objective of the current study was to use glass-fibre-reinforced polymer as a composite material. Different laminate types were tested using ABAQUS software upon verification of the work of [24]. The laminates used were quasi-isotropic (0/90/45/−45)s, angle-ply (45/−45/45/−45)s, cross-ply (0/90/0/90/0/90/0/90), and balanced laminate (−45/30/60/30/−30/45/−60/−30). In this study, the implementation of the design of experiment (DOE) seemed necessary to ensure that all parameters affecting the buckling load were optimized. We began with identifying the factors and levels followed by selecting a suitable factorial design. This study used a full factorial L₂₇ orthogonal array where 27 indicates the number of runs (simulation) required. Response surface methodology (RSM) was used to analyse the obtained data through the main effect plot, contour plot, regression, and analysis of variance (ANOVA). Additionally, the present study focused on perforated composite plates, and the supports were fixed at both ends and the orientation of the fibre was as mentioned. The cross-sectional shape considered was a C section. The novel work of this study focused on optimizing the parameters to achieve the highest probable critical buckling load for perforated composites and thin-walled members, through simulation using MINITAB and Design Epxert tools.

2. Finite Element Method

2.1. Designed Model

The model was considered a C-section channel, and the dimensions of the model are included in Figure 2. The material was glass-fibre-reinforced polymer (GFRP), which is shown in the model, and the channel had the following dimensions: a height of 80 mm, a width of 40 mm, and a length of L = 250 mm. The flanges and web had a thickness of 2.08 mm and were made up of eight 0.26 mm laminated layers. The same dimensions were used for subsequent analysis for four types of laminate, named quasi-isotropic, angle-ply, cross-ply, and balanced laminate. In the present work, the model had eight layers of GFRP material [24], and these properties were initialized into ABAQUS software. Additionally, the GFRP material had a mass density of 2200 kg/m³.

Table 1. GFRP properties used for simulation.

E1 (Mpa)	E2 (Mpa)	υ	G12 (Mpa)	G13 (Mpa)	G23 (Mpa)
38,000	8100	0.27	2000	2000	2000

shows the GFRP strut’s strength properties.

2.2. Loading and Boundary Conditions

The current simulation model was developed in ABAQUS software. All mechanical properties can be input into the engineering data before the other step. The size of the mesh was 2 mm and was used based on mesh sensitivity analysis. To have fixed-fixed end supports, all points on the bottom end had restrained rotations (ROTX, ROTY, ROTZ) and translations (UX, UY, UZ). On the upper part of the column, the load was applied at the reference point, which is shown in Figure 3. RP-1 represents the reference point at the bottom, while RP-2 is for the top section.

The buckling analysis of the linear perturbation procedure in ABAQUS was used for the elastic buckling analysis. All models utilized multi-layered, four-node elements, S4R, where each node has six degrees of freedom. The S4R shell is a four-node stress–displacement element with reduced integration. With this element, thin-walled structure geometry is expressed by a flat finite element degenerated from 3D finite element formulations of the structure’s mid-surface [24]. The x, y, and z rotation angles of the upper end were all assumed to be constrained (Figure 4). The vertical z axis was set to be unrestricted, while for the bottom end, at the reference point, all degrees of freedom were restrained. Using rigid body-pin (node) constraints, a uniform load was applied to a reference point that represented the upper nodes of the strut. The contact of the strut with the top plate is represented by this type of constraint. The FE model was given a nominal compressive load of 1.0 N. A linear eigenvalue buckling problem was used to determine the critical buckling load. The first buckling mode, which is referred to as the critical buckling load, was used in this study to generate the buckling mode shapes.

2.3. Laminate Model

The GFRP laminate in this study consisted of eight layers of unidirectional fibre-reinforced laminate. It was decided that the laminate would be a quasi-isotropic material that would have symmetrical layups using a (0/45/45/90)s stacking sequence. When subjected to loading in the plane of the laminate, this material results in the laminate behaving in the same manner as an isotropic surface [24]. The boundary condition has a significant impact on the buckling load of any engineering structure subjected to a uniformly compressive in-plane load. There is no lateral displacement or rotation along the edges if all sides of the plate are fixed, so it is more stable than simply supported [25].

In the ply system (x′, y′, z′), or the local axis system, a single laminate can be thought of as an orthotropic layer (Figure 5). The way the fibres are arranged makes an angle with the global coordinate system of the composite plate, and this angle is parallel to the local axis x′. The generalized Hooke’s Law can be used to figure out the relationship between stress and strain, which can then be used to understand the conditions of a quasi-isotropic laminate. For an individual lamina isotropic matrix, the extension and shear are ignored, and its parts do not depend on how the lamina is oriented [24]. This can be expressed as:

$$\left\{\begin{array}{c}{\sigma }_{x^{\prime }}\\ {\sigma }_{y^{\prime }}\\ {\tau }_{x^{\prime }{y}^{\prime }}\end{array}\right\}=\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& 0\\ Q_{21}& Q_{22}& 0\\ 0& 0& Q_{66}\end{array}\right] \left\{\begin{array}{c}{\epsilon }_{x^{\prime }}\\ {\epsilon }_{y^{\prime }}\\ {\gamma }_{x^{\prime }{y}^{\prime }}\end{array}\right\}$$

(1)

where

$$Q_{11}=\frac{E_1}{1-v_{12}{v}_{21}}$$

$$Q_{12}=\frac{v_{12}{E}_1}{1-v_{12}{v}_{21}}$$

$$Q_{22}=\frac{E_2}{1-v_{12}{v}_{21}}$$

$$Q_{66}=G_{12}$$

The web and flange of the strut had the same overall thickness, which was set at 2.08 mm. The ABAQUS composite layup function was used to model the GFRP strut’s laminated composite geometry. To ensure that the web and flange of the strut had identical properties, the strut was made with the same number of plies, ply thickness, lamination stacks, and orientation of fibres, and this has been illustrated in Figure 6 and Figure 7. Both web and flanges had the same orientation of fibres on the model. Three parameters were considered, which were the spacing ratio, the shape, and the ratio of opening. For each parameter, three different configurations were introduced (Figure 8).

The shapes of the hole were circular, square, and hexagonal (Figure 8). The spacing ratio and opening ratio were selected as shown below.

$$1.08<\frac{\mathrm{D}}{{\mathrm{W}}_{\mathrm{o}}}<1.5$$

(2)

$$1.25<\frac{\mathrm{W}}{{\mathrm{W}}_{\mathrm{o}}}<1.75$$

(3)

3. Design of Experiment

In the early stages of work, the design of experiment (DOE) method was used, as a method utilized in various manufacturing industries [26,27] to investigate the factors that can help designers to reach their goals and reduce the human efforts. The parameters and their levels can be changed by the designer to achieve the optimum solution. For the factorial design, at least two parameters or levels must be used to change the elements, and this must be carried out after users pick a factorial design. The tests are conducted with the parametric combination and levels marked as orthogonal arrays for a chosen factor. Following the orthogonal array helped us to determine the effects of each important part. This allowed the analysis to be carried out, which showed which factors were most important for accomplishing objectives and how these could be changed to make them even better [26,28,29].

3.1. Plan of Experiments with FEM Results

To find the best possible response, the experiments were run using the full factorial design L₂₇ orthogonal array provided by MINITAB-20 (Coventry, UK) and Design Expert-13 software (Minneapolis, MN, USA). If the orthogonal array degrees of freedom are greater than, or at least equal to the previous ones, then the test parameters and the array are selected. The current experiment’s parameters had a greater impact. An orthogonal array of 27 tests was used, and each row had a parameter assigned to it (each column had a row number). The simulations were carried out following the orthogonal array L₂₇ standard. A planned analysis should be able to examine how different composite laminate types’ critical buckling behaviour was influenced by the opening ratio (W/W_o), spacing ratio (D/W_o), and cutout shape (

Table 2. The parameters and levels.

Parameter	Levels
	1	2	3
W/Wo	1.5	1.6	1.7
D/Wo	1.3	1.4	1.5
Shapes	Circular	Square	Hexagonal

).

Using the L₂₇ orthogonal array, we were able to perform a linear buckling analysis for a variety of parameter combinations. Included here is the response for each laminate type as requested by Minitab and design expert software (

Table 3. FEM results were obtained from simulation for different types of laminate.

Run Order	W/Wo	S/Wo	SHAPE	Critical Buckling Load (N)
				Angle-Ply	Balanced	Quasi-Isotropic	Cross-Ply
1	1.7	1.5	CIRCULAR	9773.7	9077.5	8554	5611.1
2			HEXAGON	9464.4	8765.2	8271.3	5387.8
3			SQUARE	9108.1	8411.4	7963.3	5146.9
4	1.6	1.5	CIRCULAR	9430.4	8742	8330.1	5483
5			HEXAGON	9109	8411	8024.7	5233.9
6			SQUARE	8644.6	8064.3	7727.5	4998.8
7	1.5	1.5	CIRCULAR	9040.9	8351.7	8066.8	5323
8			HEXAGON	8716.2	8011.3	7746.6	5055.8
9			SQUARE	8129	7600.4	7471.9	4844.5
10	1.7	1.4	CIRCULAR	9798.9	9084.6	8495.7	5507.2
11			HEXAGON	9475	8761.3	8201.6	5280
12			SQUARE	9107.2	8397.8	7879.9	5031
13	1.6	1.4	CIRCULAR	9438.7	8738.6	8253.3	5366.6
14			HEXAGON	9100	8393	7937.1	5118.1
15			SQUARE	8639.1	8028.9	7618.4	4866
16	1.5	1.4	CIRCULAR	9029	8332.6	7974	5200.1
17			HEXAGON	8680.1	7970	7637	4927.7
18			SQUARE	8049.5	7559.5	7335.3	4688.2
19	1.7	1.3	CIRCULAR	9844.5	9107.3	8452	5410.4
20			HEXAGON	9505.9	8774.6	8145	5172.5
21			SQUARE	9131.3	8406.4	7822.2	4930.9
22	1.6	1.3	CIRCULAR	9468.9	8752.6	8190	5253.5
23			HEXAGON	9105.6	8387	7854.5	4993.4
24			SQUARE	8653	8011.6	7530.4	4744.7
25	1.5	1.3	CIRCULAR	9035.3	8329.7	7890.6	5076
26			HEXAGON	8648.2	7940.5	7528.6	4789.8
27			SQUARE	7988.4	7521.3	7213.7	4543.7

).

3.2. Response Surface Methodology (RSM)

To obtain the best possible results when dealing with multiple quantitative variables, the response surface methodology (RSM) is one of the significant methods that can be used. In RSM, the dependent variables are referred to as responses, while the independent variables or factors are referred to as predictor variables. In contrast to the statistical significance test for a single point, RSM can be used to optimize variables or factors more practically.

The RSM can be used to figure out which parameters are affected by the many variables by using the smallest possible amount of effort, real proof, and a test design that works for every point of the response variable [26,28,29]. The RSM and the relationship between the factors stated were also used in this study to find the best setting for each one. By using an RSM analysis, we were able to determine the effect of shape, opening ratio, and spacing ratio on critical buckling loads for various laminate orientations. Additionally, Kumari and Gupta [30] stated that predictive models can be developed more accurately due to RSM’s ability to reproduce results and fine-tune processes. Response surfaces in RSM are graphical representations for describing the influence of process variables and their effects on response. Therefore, the present work considered an RSM to optimize the current problem with the selected parameters. Figure 9 below depicts the RSM process in more detail.

3.3. Analysis of Variance (ANOVA)

ANOVA analysis shows the impact of parameters and their relationships on the mean square and buckling load errors at predetermined confidence levels. Using the ANOVA, one can determine the impact of each factor on the overall variance in the findings [29]. Variance analysis is a technique for comparing models based on their resources. If the means of two or more groups differ significantly, this statistical analysis is used. The ANOVA measures the influence of one or more variables by comparing the methods of different samples, and it can be used to determine whether all treatments were equally effective. The ANOVA was used to determine the influence of each parameter on the critical buckling load.

4. Finite Element Results

4.1. Model Validation

This study used Doan and Thai’s [24] buckling results to verify the proposed FE model, and the specimen’s length and cross-sectional dimension were matched and reproduced by the same model to verify the current work. The quasi-isotropic laminate and angle-ply laminate FE buckling behaviour results were compared to the corresponding numerical analysis results in order to verify the accuracy of the proposed model. The stacking sequences (0/+45/+45/+45/+45)s and (+45/+45/+45/+45)s were used for the quasi-isotropic laminate and angle-ply laminate, respectively. To begin, the model of the beam without holes was compared for both laminates to have better reliability. Table 4 shows the comparison of numerical results for the critical buckling load with those of a previous study [24].

Table 4. Comparison of present results.

Laminate Types	Critical Buckling Load (N)		Percentage Error (%)
	Previous Study [12]	Present Simulation (ABAQUS 2020)
Quasi-isotropic	11,258	11,255	0.026
Angle-ply	12,652	11,382	0.017

Table 4. Comparison of present results.

Laminate Types	Critical Buckling Load (N)		Percentage Error (%)
	Previous Study [12]	Present Simulation (ABAQUS 2020)
Quasi-isotropic	11,258	11,255	0.026
Angle-ply	12,652	11,382	0.017

4.2. Mesh Convergence Test

To investigate the effect of meshing size, this study tested five different element sizes, which ranged from 1 mm to 5 mm. The element sizes that achieved the closest results were chosen for subsequent analysis. The mesh refinement was investigated considering a single case of current work which is similar to existing work [24], and we found the percentage error of 2 mm was close to the reference value. Once satisfactory results were obtained, the same mesh size was used for all models, and slightly modified mesh was used near the holes of the structure for smoothened mesh at the boundaries of holes. Referring to

Table 5. Different element sizes test.

Element Size	Critical Buckling Load (N)	Number of Element	Percentage Error (%)
1	11,232	40,000	0.230947
2	11,255	10,000	0.026648
3	11,313	4399	0.488541
4	11,365	2520	0.950435
5	11,423	1600	1.465624

and Figure 10, the element size of 2 mm had the lowest error, and the result started to converge at this element size. Additionally, the computing time was also quite reasonable for this mesh setup. This is necessary to justify that further refinement is unnecessary and the errors are minimized.

4.3. Displacement of Laminate Structures

The obtained results from the FE analysis prove that the introduction of holes reduces the ability of a beam to withstand axial load before it fails in this case, buckling for all types of laminate. The eigenvalue shown in Figure 11, Figure 12, Figure 13 and Figure 14 represents the critical buckling load that those models can handle before they start to buckle. For all models along the longitudinal axes, the flanges were deformed in a half-wave. This may occur as a result of the compressive stresses induced on the flange plane during the loading process. In the reference model (model with no holes), local buckling can be seen clearly along the web, while in the perforated models, buckling can also be seen near the edges of holes.

Referring to the attached Figure 11, Figure 12, Figure 13 and Figure 14, the angle-ply laminate had the highest value of critical buckling compared to the others despite the presence and shape of holes. Meanwhile, the cross-ply laminate was the weakest configuration among all laminates. All models experienced Euler buckling load with the first mode along the web and flanges. A consistent pattern of a half-wave could be seen despite the existence of perforation and its shape for quasi-isotropic and cross-ply. Hence, the difference in laminate configuration does impact the thin-wall structure’s critical buckling load and also the shape of buckling.

4.3.1. Effect of Shape

Figure 15 shows the effect of cutout shapes with different laminates for a constant spacing and opening ratio. The buckling load decreased as the cutouts were introduced. The circular cutout showed the highest critical buckling load, while the square had the lowest capability of buckling load. The decrement in the load seemed significant as holes were introduced. It is believed that this was due to greater area removal on the web of the thin wall with a circular shape, as this web had the least amount of area removed.

4.3.2. Effect of Opening Ratio

Figure 16 indicates a declining trend as the ratio of opening decreased. From Equations (2) and (3), one can observe that the size of holes increased as the ratio decreased. This finding shows agreement with the previous paragraph, where the buckling load was described to decrease as greater perforation was introduced onto the web of the thin-wall GFRP composite member.

4.3.3. Effect of Spacing Ratio

Lastly, in Figure 17, the graphs show an insignificant decrement in the trend as the spacing ratio decreased. For a constant opening ratio, one can deduce that the distance between the centre of two holes becomes closer as the ratio decreases. From all mentioned figures in this section, it was concluded that for all laminates, the circular, 1.7 opening ratio and 1.5 spacing ratio showed the highest critical buckling load. However, the effect of a combination of these parameters is discussed in the next sections. Furthermore, the cross-ply laminates had the lowest critical buckling load, as shown in Figure 15, Figure 16 and Figure 17.

4.4. Nonlinear Analysis

Additionally, nonlinearity was presented because of the material behaviour. The nonlinear behaviour of the structures was observed to ensure that the FE model is valid for further applications. A geometric imperfection scale of 0.1 was applied to the linear model, acting as a base state of the model before nonlinear analysis. Additionally, the node displacement derived from the linear buckling analysis was then applied to the initial deformation, and the first buckling load was allocated as the initial load for the post-buckling analysis model [24]. As seen in Figure 18. above, the path of post-buckling inhibited relatively stable bifurcation states. Hence, using the straight-line intersection method, a straight line was approximately plotted on both linear and nonlinear equilibrium paths, therefore introducing an intersection point representing critical buckling load [24]. The ${\mathrm{R}}^2$ values on the figure above represent the coefficient of determination for the straight line based on scattered data for both pre- and post-buckling. The approximation equation portrayed was based on the applied load and displacement of the strut.

According to the obtained post-buckling analysis results, the quasi-isotropic model without holes had a critical buckling load of 12,086.53 N with a percentage of error of around 7.39% compared to linear analysis, which is still within the tolerance error of 10%. The other three configurations, which were angle-ply, balanced, and cross-ply, had a percentage of error of 0.3%, 0.03%, and 8.12%, respectively. All of them had a high determination coefficient for approximation of straight lines on both pre-buckling and post-buckling trends. Those R² values indicate the accuracy of the approximation for a curved trend in a post-buckling state. Additionally, the critical buckling load for Figure 19, Figure 20 and Figure 21 was 12,667.85 N, 12,551.1 N, and 7878.47 N, respectively.

5. Optimization Results

5.1. Response Surface Analysis

In this section, the obtained analysis results from MINITAB 20 and Design Expert-13 software will be discussed in detail. The plots that were included in this study were obtained from response surface analysis.

5.1.1. Main Effect Plot

The critical buckling load was studied in relation to several different parameters. MINITAB 20 and Design-Expert 13 software were used to create the plots. The main effects plot shows how controlled process parameters influenced critical buckling load for different laminates. As can be seen in Figure 22a–d, in the almost same plot phenomenon, there were only differences in the numbers which we obtained through the FE results. The opening ratio had a significant effect on the response, except for the cross-ply laminate. A greater ratio of openings resulted in a greater buckling load for all laminates. For all laminates, the plot shows that circular shape led to the highest possible buckling load compared to hexagonal and square shapes. The same margin between circular and hexagon can be seen for all laminates except for angle-ply, where the margin seemed to be closer. The effect of the spacing ratio on the angle-ply did not show any significant result for the buckling load. As for the balanced laminate, the response from the spacing ratio was also insignificant, while for quasi-isotropic and angle-ply laminates, a higher buckling load could be achieved as the ratio was increased.

5.1.2. Contour Plot

On the other hand, the contour plots illustrate the variation in every response with respect to the input parameters with 2D and 3D effects. Therefore, this study examined the critical buckling load variation for all laminate types. In the contour plot, the lowest region is coloured in darker blue, while darker red indicates a stronger response. As seen in Figure 23, Figure 24, Figure 25 and Figure 26, the contour had a similar pattern, which was a circular plot centred somewhere with the quasi-isotropic having a greater radius. The higher opening ratio with a combination of a 1.40 to 1.50 spacing ratio led to a higher response for all shapes of the cutout. However, at a lower value of W/W_o, circular holes showed a better response compared to differently shaped holes. For cross-ply, the plot indicates a more sensitive response as the slope of each contour was steeper compared to the quasi-isotropic laminate. This means a slight adjustment of the input variable will affect the response of the output. Both angle-ply and balanced laminates had an almost linear contour plot. According to this result, one can freely choose the spacing ratio despite the size of the holes to obtain optimum results.

From all given contour plots, one can freely choose the spacing ratio despite the size of holes to obtain optimum results. For example, for a balanced laminate, we can predict the response of different opening ratios and spacing ratios toward critical buckling load. It can be deduced from Figure 19 above that considering the three different shapes observed, the response was highest at the combination of an opening ratio of 1.70 and a spacing ratio of 1.50 and lowest at an opening ratio of 1.50 and spacing ratio 1.30. Additionally, this was applied to the other three different laminates, which showed the same behaviour at the mentioned point. However, the response might change as the input was within the maximum and minimum values of the opening ratio and spacing ratio.

5.1.3. Analysis of Variance (ANOVA)

The ANOVA made it possible to examine the influence of each variable on the overall variation in the outcomes. The ANOVA tables (

Table 6. Quasi-isotropic ANOVA.

Source	DF	Adj SS	Adj MS	F-Value	p-Value
Model	11	3,257,334	296,121	11,702.66	0.000
Linear	4	3,244,740	811,185	32,057.89	0.000
W/Wo	1	134,5073	1,345,073	53,157.07	0.000
S/Wo	1	129,914	129,914	5134.18	0.000
SHAPE	2	176,9753	884,876	34,970.16	0.000
Square	2	1738	869	34.35	0.000
W/Wo*W/Wo	1	1478	1478	58.41	0.000
S/Wo*S/Wo	1	260	260	10.29	0.006
2-Way Interaction	5	10,856	2171	85.80	0.000
W/Wo*S/Wo	1	6674	6674	263.76	0.000
W/Wo*SHAPE	2	1532	766	30.28	0.000
S/Wo*SHAPE	2	2649	1324	52.34	0.000
Error	15	380	25
Total	26

,

Table 7. Angle-ply ANOVA.

Source	DF	Adj SS	Adj MS	F-Value	p-Value
Model	11	6,640,780	603,707	2127.51	0.000
Linear	4	6,542,370	1,635,592	5763.96	0.000
W/Wo	1	3,460,554	3,460,554	12,195.28	0.000
S/Wo	1	69	69	0.24	0.629
SHAPE	2	3,081,747	1,540,873	5430.17	0.000
Square	2	8385	4192	14.77	0.000
W/Wo*W/Wo	1	7896	7896	27.83	0.000
S/Wo*S/Wo	1	488	488	1.72	0.209
2-Way Interaction	5	90,025	18,005	63.45	0.000
W/Wo*S/Wo	1	10,191	10,191	35.91	0.000
W/Wo*SHAPE	2	75,981	37,991	133.88	0.000
S/Wo*SHAPE	2	3853	1926	6.79	0.008
Error	15	4256	284
Total	26	6,645,036

,

Table 8. Cross-ply laminate ANOVA.

Source	D	Adj SS	Adj MS	F-Value	p-Value
Model	11	1,870,705	170,064	6065.49	0.000
Linear	4	1,864,835	466,209	16,627.76	0.000
W/Wo	1	509,713	509,713	18,179.39	0.000
S/Wo	1	261,581	261,581	9329.54	0.000
SHAPE	2	1,093,541	546,770	19,501.06	0.000
Square	2	681	340	12.14	0.001
W/Wo*W/Wo	1	664	664	23.69	0.000
S/Wo*S/Wo	1	17	17	0.59	0.454
2-Way Interaction	5	5189	1038	37.01	0.000
W/Wo*S/Wo	1	2754	2754	98.23	0.000
W/Wo*SHAPE	2	1702	851	30.36	0.000
S/Wo*SHAPE	2	732	366	13.06	0.001
Error	15	421	28
Total	26

and

Table 9. Balanced laminate ANOVA.

Source	DF	Adj SS	Adj MS	F-Value	p-Value
Model	11	5,237,892	476,172	7286.76	0.000
Linear	4	5,217,484	1,304,371	19,960.51	0.000
W/Wo	1	2,855,333	2,855,333	43,694.55	0.000
S/Wo	1	2307	2307	35.31	0.000
SHAPE	2	2,359,843	1,179,922	18,056.09	0.000
Square	2	8271	4136	63.28	0.000
W/Wo*W/Wo	1	7942	7942	121.54	0.000
S/Wo*S/Wo	1	329	329	5.03	0.040
2-Way Interaction	5	12,137	2427	37.15	0.000
W/Wo*S/Wo	1	3540	3540	54.17	0.000
W/Wo*SHAPE	2	6514	3257	49.84	0.000
S/Wo*SHAPE	2	2084	1042	15.94	0.000
Error	15	980	65
Total	26	5,238,872

) illustrate the interactive effect of each parameter for all possible factor combinations. The obtained ANOVA analysis was taken directly from MINITAB software. The W/W_o, which had the greatest F-value for all laminate types, indicates that it is the parameter that most influenced one-way interaction. The F-value indicates how much the parameter is associated with the response. This once again supports the previously stated hypothesis that the size of the hole affects the eigenvalue load significantly.

The combination of W/W_o and S/W_o had the highest F-value for quasi-isotropic, cross-ply, and balanced laminates. In other words, when it comes to creating a composite thin-wall structure, combining these two elements had the greatest impact. Angle-ply, on the other hand, can utilize a variety of different combinations. There was a strong correlation between the W/W_o ratio and cutout shape, which affects buckling loads.

5.1.4. Regression Equation

The L₂₇ orthogonal array through MINITAB software was used to identify an equation-based statistical model for the current study. A linear polynomial model (regression equations) was used to represent each component. The equation was obtained separately for each cutout shape as a categorical factor and different laminate types. The linear regression equation included below can be used to calculate the expected results and percentage variance for all test cases.

When determining the extent to which particular independent variables are influencing dependent variables, regression analysis is a useful statistical method that can be used to obtain optimized data. For each laminate type, we can see that there was a consistency in the multiplier for S/W_o*S/W_o, W/W_o*S/W_o, and W/W_o*W/W_o. It is believed that the interaction between those parameters affects the critical buckling load of the strut.

For angle-ply,
CIRCULAR	P_cr(N)	=	−10522 + 19542 W/Wo + 1964 S/Wo—3628 W/Wo*W/Wo + 902 S/Wo*S/Wo—2914 W/Wo*S/Wo
HEXAGON	P_cr(N)	=	−11411 + 19690 W/Wo + 2186 S/Wo—3628 W/Wo*W/Wo + 902 S/Wo*S/Wo—2914 W/Wo*S/Wo
SQUARE	P_cr(N)	=	−14156 + 20988 W/Wo + 2318 S/Wo—3628 W/Wo*W/Wo + 902 S/Wo*S/Wo—2914 W/Wo*S/Wo
For balanced laminate,
CIRCULAR	P_cr(N)	=	−8939 + 17806 W/Wo + 645 S/Wo—3638 W/Wo*W/Wo + 740 S/Wo*S/Wo—1717 W/Wo*S/Wo
HEXAGON	P_cr(N)	=	−9856 + 18013 W/Wo + 818 S/Wo—3638 W/Wo*W/Wo + 740 S/Wo*S/Wo—1717 W/Wo*S/Wo
SQUARE	P_cr(N)	=	−10769 + 18271 W/Wo + 904 S/Wo—3638 W/Wo*W/Wo + 740 S/Wo*S/Wo—1717 W/Wo*S/Wo
For cross-ply,
CIRCULAR	P_cr(N)	=	−4456 + 7037 W/Wo + 3088 S/Wo—1052 W/Wo*W/Wo + 166 S/Wo*S/Wo—1515 W/Wo*S/Wo
HEXAGON	P_cr(N)	=	−5179 + 7266 W/Wo + 3162 S/Wo—1052 W/Wo*W/Wo + 166 S/Wo*S/Wo—1515 W/Wo*S/Wo
SQUARE	P_cr(N)	=	−5442 + 7209 W/Wo + 3244 S/Wo—1052 W/Wo*W/Wo + 166 S/Wo*S/Wo—1515 W/Wo*S/Wo
For quasi-isotropic,
CIRCULAR	P_cr(N)	=	−4921 + 10941 W/Wo + 2626 S/Wo—1569 W/Wo*W/Wo + 659 S/Wo*S/Wo—2358 W/Wo*S/Wo
HEXAGON	P_cr(N)	=	−5825 + 11167 W/Wo + 2786 S/Wo—1569 W/Wo*W/Wo + 659 S/Wo*S/Wo—2358 W/Wo*S/Wo
SQUARE	P_cr(N)	=	−6162 + 11065 W/Wo + 2922 S/Wo—1569 W/Wo*W/Wo + 659 S/Wo*S/Wo—2358 W/Wo*S/Wo

Regarding the contour plot, we needed both parameters to obtain the optimal response. This obtained equation can be used to predict the response for any point of input parameters. To illustrate this, all three equations for different shapes of a cutout from the angle-ply laminate were tested to determine their reliability and compared with the results of previously obtained data. Figure 27 clearly shows that there was a relatively small error between the result obtained from the equation and the FE model. Therefore, it was concluded that these equations are useful to predict the response for any random point of input parameters.

6. Conclusions

In short, ABAQUS software was used to numerically study thin-walled composite structure buckling. The critical buckling load was calculated using eigenvalue linear buckling. Four laminates were studied to see how parameters affect response and critical buckling load. Then, the results were optimized. This study used a full factorial orthogonal array first, then obtained eigenvalue buckling data for all cases to better understand the data. The obtained data were then analysed using MINITAB and Design Expert software using response surface methodologies to verify model estimation, graphs, and variable contours and levels. The W/W_o parameter had the highest F-value in all laminate types, indicating it is the most influential for one-way interaction. In quasi-isotropic, cross-ply, and balanced laminates, W/W_o and S/W_o had the highest F-value. The design of the composite thin-wall structure was influenced by these two combinations. Angle-ply can be combined in many ways. W/W_o and the cutout shape had the greatest impact on the buckling load. This study also found the best parameter combinations for different laminates in terms of critical buckling. An optimization analysis based on parameters showed that range changes affect buckling load values.

All laminates had an abrupt trend except the cross-ply laminate, which had a significant opening ratio. The opening ratio increased buckling loads, affecting response. Round laminates buckled more than hexagonal or square laminates, except for angle-ply, which had a closer margin between the circular and hexagonal shapes. The aspect ratio did not affect buckling load. Increased spacing ratio increased buckling load for quasi-isotropic and cross-ply laminates. The trend was almost flat for both angle-ply and balanced laminate, so the spacing ratio had little effect. Because both factors are steeper, the cross-ply laminate is adjustable. The QIS had a moderate effect on buckling load capacity as it increased gradually. These variables combine to determine a strut’s buckling load. A contour plot requires both parameters to work properly. These equations can be used to predict the response at any given point in the input parameters. To demonstrate this, all three equations for different shapes of angle-ply laminate cutouts were tested. The results were then compared to existing data. Finally, the response surface methodologies showed optimized results from the study’s analysis. As shown in the graph, the regression equation yielded an optimal response. It shows how the factors affect the critical buckling load and also shows the best combination for a higher critical buckling load.