Studies on Numerical Buckling Analysis of Cellulose Microfibrils Reinforced Polymer Composites

Scientists are drawn to the new green composites because they may demonstrate qualities that are comparable to those of composites made of synthetic fibers due to concerns about environmental contamination. In this work, the potential for using the produced green composite in different buckling load-bearing structural applications is explored. The work on composite buckling characteristics is vital because one needs to know the composite’s structural stability since buckling leads to structural instability. The buckling properties of composite specimens with epoxy as the matrix and chemically treated cellulose microfibrils as reinforcements are examined numerically in this study when exposed to axial compressive stress. The numerical model is first created based on the finite element method model. Its validity is checked using ANSYS software by contrasting the critical buckling loads determined through research for three samples. The numerical findings acquired using the finite element method are then contrasted with those produced using the regression equation derived from the ANOVA. The utilization of the created green composite in different buckling load-bearing structural applications is investigated in this study. As a result of the green composite’s unaltered buckling properties compared to synthetic composites, it has the potential to replace numerous synthetic composites, improving environmental sustainability.


Introduction
Polymer composites are used in many structural applications. Thus, it is crucial to promote green composites comprised of natural fibers to conserve the environment for future generations. The green composite produces fewer carbon emissions during destruction than a synthetic composite since green ingredients were employed to make it. As a result, there is less environmental contamination and a more stable, green ecosystem. Due to the use of natural fiber composites in load-bearing applications, including automotive, construction and others, researchers have been studying natural fiber composites under static and dynamic loads [1]. Due to the structural instability of composite structures, the structure will attain a failure in various ways, such as fatigue failure and buckling failure. Due to the axial compression pressure or heat load, composite constructions are frequently susceptible to buckling failure [2]. The natural fiber polymer composite constructions may experience an axial compression load throughout their service since they are employed in a variety of load-bearing applications [3]. A natural fiber-reinforced polymer composite Several researchers looked at the effect of the coupling agent sodium hydroxide in alkali treatment. The findings demonstrated that varying concentrations of sodium hydroxide have distinct effects on fiber surface morphology [16][17][18]. Sodium hydroxide was used as an alkali treatment because it converts fibers' smooth surface to a rough surface by making undulation on the surface. It indicated that the matrix in the cell wall was removed after alkali treatment. The author reported that the pore sizes and void ratios vary depending on cellulose origin and as well as treatment history. These pores are important for dissolution because it provides spaces for the diffusion of solvent chemicals into fibers. Na-Cellulose-I lattices were formed when excess sodium hydroxide was removed. Cellulose-II, having a more stable lattice in comparison to Cellulose-I, was formed by rinsing the cellulose with water and removing the sodium ions from the cellulose [19][20][21][22]. The fiber is affected by the alkaline treatment in two ways. As a result, mechanical interlocking was improved, and the surface roughness and cellulose exposure to reaction sites were both increased. It was observed that mechanical interlocking has a long-term influence on fiber mechanical behavior, notably the strength and stiffness of fibers [23][24][25].
Kabir et al. [26] studied the surface treatments of natural fibers to maximize the composites' bonding strength and stress transferability. The overall mechanical properties of natural fiber-reinforced polymer composites were highly dependent on the morphology, aspect ratio, hydrophilic tendency and dimensional stability of the fibers. Harish et al. [27] evaluated the mechanical properties of coir-reinforced composites. The tensile results showed the suitability of coir fiber for low load-bearing applications. Rashed et al. [28] studied the tensile strength of jute fiber-reinforced composites. The effects of parameters, such as alkali treatment (compared to no treatment), fiber size (1, 2 and 4 mm) and fiber loading (5, 10 and 15 wt%) on the tensile strength were considered. They analyzed the tensile behavior and performed fractographic observations. Akil et al. [29] reviewed kenaf fiber-reinforced composites. Kenaf fibers are readily available and are used as reinforcement in various ways. NaOH-treated kenaf fiber increases the tensile and flexural properties of epoxy composites, but the thermal resistance decreases compared to untreated kenaf fiber-reinforced composites. Senthilkumar et al. [30] evaluated the mechanical and vibrational properties of pineapple leaf fiber-reinforced composites. They prepared the PALF polyester composites by the hand lay-up process and then compressed the materials using a compression testing machine. They analyzed the tensile, flexural and vibration results and found that an increase in PALF reinforcement increases the mechanical strength and reduces the damping ratio of the composites. They suggested that 45 wt% PALF composites are better suited for structural applications.
Ozbek [31] investigated the effect of silica nanoparticles (NS) on the buckling characteristics of Kevlar/epoxy fiber-reinforced composite laminates. The results revealed that decreases in the length of samples resulted in significant increases in axial and lateral buckling characteristics. Rozylo et al. [32] studied the failure phenomenon of carbon fiber-reinforced plastic columns with three distinct lay-up patterns of the composite. It was observed that the dominant failure form occurred in the region representing the end sections of the composite structures (short 180 mm columns). Moreover, the delamination phenomenon usually occurs just before structural failure. Debski et al. [33] investigated the effect of eccentric compressive load on the stability, critical states and load-carrying capacity of thin-walled composite Z-profiles. Compressive fiber damage began in Plies 2 and 7 in the squeezed composite structure. The damage process evolved into a complex failure mechanism, causing the structure to lose its load-carrying capability when all mixed damage initiation requirements were met in the damaged area.
Although many researchers studied the mechanical and flexural characteristics of different fiber/polymer combinations, they still have limited knowledge related to the buckling behavior of composites reinforced with natural cellulose microfibrils. However, they have not carried out optimization methods such as the response surface approach to determine the ideal critical buckling load for natural fiber composites. In the present study, the critical buckling stress of a composite reinforced with cellulose microfibrils is analyzed using Box-Behnken response surface design, and the desired parameter result is optimized. The research flowchart is mentioned in Figure 1. The impact of three different parameters, including fiber volume% (w/w), fiber diameter (µm), and sodium hydroxide % (w/w) on buckling behavior (dependent variable), is investigated.
Although many researchers studied the mechanical and flexural characteristics of different fiber/polymer combinations, they still have limited knowledge related to the buckling behavior of composites reinforced with natural cellulose microfibrils. However, they have not carried out optimization methods such as the response surface approach to determine the ideal critical buckling load for natural fiber composites. In the present study, the critical buckling stress of a composite reinforced with cellulose microfibrils is analyzed using Box--Behnken response surface design, and the desired parameter result is optimized. The research flowchart is mentioned in Figure 1. The impact of three different parameters, including fiber volume% (w/w), fiber diameter (µ m), and sodium hydroxide % (w/w) on buckling behavior (dependent variable), is investigated.

Methodology
Raw banana stem fibers were processed to micron size, acquired from local vendors (ECO Green unit, Coimbatore, TN, India), and used as a filler. Chemicals such as NaOH, HCl and demineralized water were used to remove cellulose microfibrils (SRL Chemicals, Sigma Aldrich, and NICE Chemicals, Chennai, TN, India). Epoxy LY556 and hardener HY951 were purchased from S.M. Composites, Chennai, TN, India.
Raw banana fibers were cut to 4-5 mm lengths and prewashed with demineralized water to remove dirt and impurities. Then, the fibers were air-dried for two days to remove excess moisture. Fully dried banana fibers were powdered by the pulverizing process (Saral Pulverizer, Gujarat, India) and size separation was performed using a sieve shaker. Figure 2 describes the procedure for the chemical treatment of powdered banana fiber to prepare cellulose microfibrils for reinforcement of epoxy composites. The chopped banana fibers were pretreated with different w/w percentages of sodium hydroxide (NaOH) solution for 2 h. Then, the fibers were washed several times with distilled water. The pretreated banana fibers were hydrolyzed using a 1M HCl solution at 80 °C ± 5 °C for 2 h. Then, the fibers were washed several times with demineralized water. The acid-hydrolyzed fibers were treated again with a 2% (w/w) NaOH solution for 2 h at 60 °C ± 5 °C. The acid-alkali-treated fibers were washed several times with demineralized water until the pH reached 7. These acid-alkali-treated fibers had more cellulose microfibrils and less pectin, hemicelluloses and lignin.

Methodology
Raw banana stem fibers were processed to micron size, acquired from local vendors (ECO Green unit, Coimbatore, TN, India), and used as a filler. Chemicals such as NaOH, HCl and demineralized water were used to remove cellulose microfibrils (SRL Chemicals, Sigma Aldrich, and NICE Chemicals, Chennai, TN, India). Epoxy LY556 and hardener HY951 were purchased from S.M. Composites, Chennai, TN, India.
Raw banana fibers were cut to 4-5 mm lengths and prewashed with demineralized water to remove dirt and impurities. Then, the fibers were air-dried for two days to remove excess moisture. Fully dried banana fibers were powdered by the pulverizing process (Saral Pulverizer, Gujarat, India) and size separation was performed using a sieve shaker. Figure 2 describes the procedure for the chemical treatment of powdered banana fiber to prepare cellulose microfibrils for reinforcement of epoxy composites. The chopped banana fibers were pretreated with different w/w percentages of sodium hydroxide (NaOH) solution for 2 h. Then, the fibers were washed several times with distilled water. The pretreated banana fibers were hydrolyzed using a 1M HCl solution at 80 • C ± 5 • C for 2 h. Then, the fibers were washed several times with demineralized water. The acid-hydrolyzed fibers were treated again with a 2% (w/w) NaOH solution for 2 h at 60 • C ± 5 • C. The acid-alkali-treated fibers were washed several times with demineralized water until the pH reached 7. These acid-alkali-treated fibers had more cellulose microfibrils and less pectin, hemicelluloses and lignin.
The Euler critical formula is utilized for predicting the buckling load of the supporting plate with different boundary parameters. Euler's Critical Load Formula for Plate is presented in Equations (1) and (2).
where E represents Young's modulus of Text plates, ν represents Poisson's ratio, k c = a/b Buckling Coefficient, h-thickness, D-Flexural rigidity of the plate per unit length and N xcr -critical buckling load. The plate thickness, width, Young's modulus and Poisson's ratio are a = 100 mm, E = 5249.27 MPa, v = 0.34, h = 0.75 mm and a/b = 1, respectively. One can theoretically calculate the critical buckling load using this Euler critical formula and these property values. The critical value factor in this expression represents the Euler load for a strip of unit width and length 'a'. The second factor, 'b', denotes the proportion of greater stability gained by the continuous plate compared with that of an isolated strip. Table 1 shows that the epoxy reinforced with 20% NaOH-treated filler (run no. 14) exhibits higher tensile stress and modulus compared to the other samples. The Euler critical formula is utilized for predicting the buckling load of the supporting plate with different boundary parameters. Euler's Critical Load Formula for Plate is presented in Equations (1) and (2).
where E represents Young's modulus of Text plates, ν represents Poisson's ratio, kc = a/b Buckling Coefficient, h-thickness, D-Flexural rigidity of the plate per unit length and Nxcrcritical buckling load. The plate thickness, width, Young's modulus and Poisson's ratio are a = 100 mm, E = 5249.27 MPa, v = 0.34, h = 0.75 mm and a/b = 1, respectively. One can theoretically calculate the critical buckling load using this Euler critical formula and these property values. The critical value factor in this expression represents the Euler load for a strip of unit width and length' a'. The second factor, 'b', denotes the proportion of greater stability gained by the continuous plate compared with that of an isolated strip. Table 1 shows that the epoxy reinforced with 20% NaOH-treated filler (run no. 14) exhibits higher

Response Surface Methodology
In the study, a sophisticated statistically verified prediction model is used to conduct and evaluate trials to determine the optimal combination and the impact of various elements, such as a change in sodium hydroxide %, a change in fiber diameter and a change in volume percentage, on the buckling characteristics. The experimental design in this work employs the Response Surface Methodology's Box-Behnken design. Table 2 lists the process parameters and the three levels at which each parameter is evaluated. Box-Behnken design is used to generate higher-order response surfaces using fewer required runs than a standard factorial technique [34]. The design uses the twelve middle edge nodes and three center nodes to fit a 2nd order equation. NaOH% (w/w) 15 17. 5 20 In MINITAB software, these levels and factors are utilized to frame the different combinations and represent the results, as shown in Table 3. Table 3. Response Surface Methodology.

Buckling Analysis
This examination is performed to find out how well laminated composite samples buckle under axial compressive stresses. Software based on the finite element method, ANSYS V.18, is used to do the numerical analysis. The critical buckling loads and load factors of 15 composite samples are determined. Shell 181, a four-node element with six degrees of freedom at each node, was utilized to simulate the laminated composite in this investigation. The lay-up technique is applied to define the orientation of the sample and the number of layers under sections. An evenly distributed 10 kN compression load is applied. An axial compression load is applied to a mesh-supported boundary condition. The plate is supported with loaded edges, and all the loaded edges are simply supported. Uz = 0, as shown in Figure 3a. Simply supported end conditions are used, as shown in Figure 3b. The critical buckling load is determined by multiplying the load factor by the applied load [35]. It reveals the buckling factor found numerically using ANSYS 18; the results are mentioned in Table 4. The deformation was applied in the longitudinal direction to the upper & side parts of the model, representing the upper grip. All the other degrees of freedom were fixed for this part. The displacement was a linear function up to a chosen maximum deformation limit. The essential output of this model is the critical buckling force, i.e., the force for the onset of buckling. The result of the buckling deformation is shown in Figure 3b. The model allows the whole bottom grip of the specimen to move freely in the X-Y plane. This simplification is expected to give lower buckling force values in the model than in the experiments. In general, it is observed that the calculated values from the model are smaller than the measured values from the experiments.

Buckling Analysis
This examination is performed to find out how well laminated composite samples buckle under axial compressive stresses. Software based on the finite element method, ANSYS V.18, is used to do the numerical analysis. The critical buckling loads and load factors of 15 composite samples are determined. Shell 181, a four-node element with six degrees of freedom at each node, was utilized to simulate the laminated composite in this investigation. The lay-up technique is applied to define the orientation of the sample and the number of layers under sections. An evenly distributed 10 kN compression load is applied. An axial compression load is applied to a mesh-supported boundary condition. The plate is supported with loaded edges, and all the loaded edges are simply supported. Uz = 0, as shown in Figure 3a. Simply supported end conditions are used, as shown in Figure 3b. The critical buckling load is determined by multiplying the load factor by the applied load [35]. It reveals the buckling factor found numerically using ANSYS 18; the results are mentioned in Table 4. The deformation was applied in the longitudinal direction to the upper & side parts of the model, representing the upper grip. All the other degrees of freedom were fixed for this part. The displacement was a linear function up to a chosen maximum deformation limit. The essential output of this model is the critical buckling force, i.e., the force for the onset of buckling. The result of the buckling deformation is shown in Figure 3b. The model allows the whole bottom grip of the specimen to move freely in the X-Y plane. This simplification is expected to give lower buckling force values in the model than in the experiments. In general, it is observed that the calculated values from the model are smaller than the measured values from the experiments.

Numerical Analysis
All 15 samples are numerically simulated and mentioned in Table 5. ANSYS 18 is utilized to calculate the critical buckling loads for each of the 15 samples. An analysis of variance in MINITAB software is carried out using data from the numerical simulation ANSYS 18. Maximum critical buckling loads are discovered in the samples with a sodium hydroxide ratio of 18.03% w/w, a fiber diameter of 333.3 µm, and a fiber volume of 4.14% w/w based on numerical simulation findings.

Comparison of Numerical and Theoretical Results
In order to validate, the developed numerical model is compared with the theoretical model available in the literature [34]. The finite element method/developed numerical model is validated by comparing numerical and theoretical critical buckling loads for three samples using ANSYS V.18. The validation indicates a positive sign as the error percentage is less than two. The error percentage of a random three from 15 samples is listed in Table 6.

Analysis of Variance Assessment for the Regression Model and Regression Equation
Utilizing Minitab software, the Box-Behnken design of the Response Surface Methodology approach is used to investigate the major impact of three variables: sodium hydroxide percent (w/w), fiber diameter and fiber volume on the critical buckling load of composite plates. The quadratic polynomial regression equation for the critical buckling load is derived from the fitting curve, as shown in Equation (3)

Comparing the Critical Buckling Load of the Finite Element Analysis and the Regression Equation
The resulting regression equation is then used to compute the critical buckling load for each combination. Table 7 details the results of this comparison with those obtained using finite element analysis software.  Table 7 shows that, except for a few circumstances, the average variation is between 3-5 and less than 11 percent. The comparison proves the validity of the created regression model and the derived regression equation. This demonstrates the stability of the finite element method model, and a similar model is utilized for analyses of the same type, as shown in Figure 4. The critical buckling load for each possible combination may also be calculated using this regression equation. The ANOVA is also carried out for the critical buckling load quadratic model. The p-value for the model is <0.02, demonstrating its importance. The modified R 2 value of 0.82 is quite similar to the projected R 2 value of 0.93. The model is noteworthy since the gap between the anticipated R 2 value and the adjusted R 2 value is likewise smaller than 0.11. The model summary values are shown in Table 8. Figure 4 presents the predicted vs. actual values of R 2 . The points in the graph are a bit scattered, which indicates that the predicted and actual values are a bit wider, but the points are still within a range.
3-5 and less than 11 percent. The comparison proves the validity of the created regression model and the derived regression equation. This demonstrates the stability of the finite element method model, and a similar model is utilized for analyses of the same type, as shown in Figure 4. The critical buckling load for each possible combination may also be calculated using this regression equation. The ANOVA is also carried out for the critical buckling load quadratic model. The p-value for the model is <0.02, demonstrating its importance. The modified R 2 value of 0.82 is quite similar to the projected R 2 value of 0.93. The model is noteworthy since the gap between the anticipated R 2 value and the adjusted R 2 value is likewise smaller than 0.11. The model summary values are shown in Table 8. Figure 4 presents the predicted vs. actual values of R 2 . The points in the graph are a bit scattered, which indicates that the predicted and actual values are a bit wider, but the points are still within a range.  Figure 5 depicts the contour plot exhibiting the influence of fiber diameter and volume on buckling load while maintaining a sodium hydroxide concentration of 17.5 percent w/w. According to the graph, the maximum buckling load is attained at a fiber volume of 4% w/w ratio and a fiber diameter of 375 µm. It is found that the fiber volume should be a 4% w/w ratio to have the greatest critical buckling load and that the fiber diameter does not affect the tensile critical buckling load, shown in Figure 5.   Figure 5 depicts the contour plot exhibiting the influence of fiber diameter and volume on buckling load while maintaining a sodium hydroxide concentration of 17.5 percent w/w. According to the graph, the maximum buckling load is attained at a fiber volume of 4% w/w ratio and a fiber diameter of 375 µm. It is found that the fiber volume should be a 4% w/w ratio to have the greatest critical buckling load and that the fiber diameter does not affect the tensile critical buckling load, shown in Figure 5.

Effect of Parameter on Critical Buckling Load
The effects of sodium hydroxide percentage, fiber size and fiber volume on the buckling load of chemically treated cellulose microfibril-reinforced epoxy composites are investigated using a three levels/factors Box-Behnken design model. Contour plots are created using the model to display the primary and interactive impacts between the response and independent variables. To better comprehend the interactions between the variables and to identify the ideal conditions, these graphs are created between two independent variables while keeping one variable constant. Figure 3 depicts the contour plot influence of fiber diameter and volume on buckling load while maintaining a sodium hydroxide percentage of 17.5 percent w/w. According to the graph, the maximum buckling load is attained at a fiber volume of 4% w/w ratio and a fiber diameter of 375 µm. Figure 6 depicts the contour response graph indicating sodium hydroxide percentage and fiber diameter influence buckling load when fiber volume is maintained at 4% w/w. Maximum buckling load is attained at 18.5% sodium hydroxide percentage by weight. Figure 7 depicts the contour response plot of the sodium hydroxide percentage and fiber volume influence on the buckling load with a fiber diameter of 375 µm. Figure 7 shows that the maximum buckling load is reached at a sodium hydroxide percentage of 18.5% and a fiber volume of approximately 4%. Surface plot effects demonstrate that the natural surface area of fiber grows as the sodium hydroxide percentage increases. This increase in the surface area helps the matrix to wet the fiber effectively and creates an excellent connection with the natural fiber [36], as shown in Figure 6. Additionally, chemical processing enhances the possibility of exposing more cellulose fibrils to the matrix material, increasing the mechanical strength of composites reinforced with cellulose fibrils, as shown in Figure 7. NaOH treatment is a process of alkalization on natural fibers, which may improve the bonding. It is indicated that the best chemical treatment is modifying the surface of natural fibers offering better adhesion with resin polymer composites. In addition, the mechanical properties were observed to enhance with the increasing NaOH percentage [37].

Verification and Optimization of the Model
Derringer's desirability function optimization process is used to obtain the optimal parameters to produce the largest critical buckling load. It depicts the desirability ramp for maximizing the input variables to provide the optimum results, as shown in Figure 7. In order to obtain the best results for the maximum critical buckling load of 1361.61 N, it is recommended to modify the input variables of fiber diameter to 333.3 m, sodium hydroxide ratio to 18.03% w/w, and fiber volume to 4.14% w/w ratio, as shown in Tables 9 and 10.

Verification and Optimization of the Model
Derringer's desirability function optimization process is used to obtain the optimal parameters to produce the largest critical buckling load. It depicts the desirability ramp for maximizing the input variables to provide the optimum results, as shown in Figure 7. In order to obtain the best results for the maximum critical buckling load of 1361.61 N, it is recommended to modify the input variables of fiber diameter to 333.3 m, sodium hydroxide ratio to 18.03% w/w, and fiber volume to 4.14% w/w ratio, as shown in Tables 9 and 10.  Tables 11 and 12 show the critical buckling load calculated with a numerical model and the percentage departure from the optimal condition to validate the optimization result. The deviation between the model and the experiments is 1% between the average of experiments and the model. The model shows a significantly better correlation for the analysis. In general, it is observed that the calculated values from the model are smaller than the measured values from the experiments. The verification of optimization solutions is conducted via the finite element method. The buckling factor is 135.11, and the buckling load is 1351.1 N. The critical buckling load was significant depending on the thickness of the plate. As is evident, the critical buckling load was significantly larger with an increase in the thickness of the plate [38]. Table 13 shows the variation of the buckling factor depending on the sodium hydroxide percentage and fiber volume percentage of the reinforced composite. The buckling factor is high at higher concentrations of the sodium hydroxide and fiber volume percentage.

Conclusions
The Box-Behnken design approach is used in the present work to achieve optimum process parameters for the critical buckling load of polymer composite. Using a threelevel/parameter model, the variables sodium hydroxide percent (w/w), fiber diameter (µm) and fiber volume percent (w/w) are taken into account. Buckling load results are tabulated and analyzed as well as optimized with Minitab software. The analytical results and the expected results show excellent agreement. The contour diagram provides a clear explanation of how the factors affect the critical buckling load. Its optimized value is 1361.64 N at the input variables of 18.30% w/w ratio of sodium hydroxide, 333.3 µm of fiber diameter and 4.14% w/w of fiber volume. This work demonstrates how well Box-Behnken design can be used to model cellulose microfibril-reinforced polymer composites and achieve optimal outcomes quickly and with a small number of experimental runs.