Research on the Buckling Behavior of Functionally Graded Plates with Stiffeners Based on the Third-Order Shear Deformation Theory

This paper presents a finite element formulation to study the mechanical buckling of stiffened functionally graded material (FGM) plates. The approach is based on a third-order shear deformation theory (TSDT) introduced by Guangyu Shi. The material properties of the plate were assumed to be varied in the thickness direction by a power law distribution, but the material of the stiffener was the same as that of the one of the bottom surface where the stiffener was placed. A parametric study was carried out to highlight the effect of material distribution, the thickness-to-width ratio, and stiffener parameters on the buckling characteristics of the stiffened FGM plates. Numerical results showed that the addition of stiffener to the FGM plate could significantly reduce the weight of the FGM plate but that both the FGM plates with and without stiffener had equally high strength in the same boundary condition and compression loading.


Introduction
The common way to achieve higher strength for functionally graded material (FGM) plates and shell without stiffeners is to either increase the thickness of this structure or to add stiffeners. The weight of the unstiffened structure will become higher with increasing thickness, but reinforcement with stiffener will reduce the weight as well as the cost of this structure. For this reason, using stiffeners is the best method in special cases such as ship building, bridge construction, aerospace, marine, and so on.
To give more useful information about the application in practice, the buckling behavior of composite structures has received much attention from scientists. Broekel and Gangadharaprusty [1] used experimental and theoretical solutions to study the mechanical responses of stiffened and unstiffened composite panels subjected to a uniform transverse loading. Liu and Wang [2] explored the elastic buckling of a plate reinforced by stiffener under in-plane loading. ANSYS modeling was applied to find the optimal height, number, and arrangement of stiffeners. Ueda et al. [3] proposed an analytical approach to research the buckling and deflection responses of a stiffened plate, which has a deflection under out-of-plane and in-plane loads. Danielson et al. [4] presented a combination of von Karman and nonlinear beam theories to predict the buckling behavior of a stiffened plate subjected to axial

Finite Element Formulation for Mechanical Buckling of Stiffened FGM Plates
We considered a functionally graded material stiffened plate composed of ceramic and metal phases. The material on the top surface of this plate was full of ceramic and was graded to metal at the bottom surface of the plate by the power law distribution. The stiffener was placed at the bottom surface, and it was full of metal. This meant that the material of the stiffener and the bottom surface was the same.
The thickness, length, and width of the plate have been noted as h, a, and b, respectively, while the depth and width of the stiffener are h s and b s , respectively, as sketched in Figure 1. The material was graded by the power law distribution, and it was used for describing the volume fraction of the ceramic (V c ) and the metal (V m ) as follows [23,26]: and where, h is the thickness of the plate; n is the gradient index (n ≥ 0); z is the thickness coordinate variable (−h/2 ≤ z ≤ h/2); and subscripts c and m represent the ceramic and metal constituents, respectively.  In this study, the Young's modulus E and the Poisson's ratio ν change through the z-direction as follows [23,26]: Using the Shi shear deformation theory [22,23], the FGM plate model has the following displacement field (u, v, w): u(x, y, z) = u 0 (x, y) v(x, y, z) = v 0 (x, y) + 5 4 z − 4 3h 2 z 3 φ y (x, y) + 1 4 z − 5 3h 2 z 3 w 0,y (3b) w(x, y, z) = w 0 (x, y) where u 0 , v 0 , and w 0 represent the displacements at z = 0 (the mid-plane of a plate); φ x and φ y are the transverse normal rotations of the y and x axes; the comma denotes the differentiation with respect to x and y coordinates. Four nodes per element, seven degrees of freedom per node were used for this problem. The displacement vector of node i for plate element is as follows: Because of the degree of freedom, w had the additional first derivative components ∂w 0i ∂x , ∂w 0i ∂y . Therefore, in order to guarantee the continuous condition of displacement w and its first derivative components at each node, we had to approximate the displacement w by Hermite interpolation functions. The other four degrees of freedom were approximated by Lagrange interpolation functions.
The displacements of the plate in this approach may be expressed as follows:  (8) where N i are Lagrange interpolating functions and H i are Hermite interpolating functions. The displacement vector is interpolated through the element's nodal displacement vector as follows: where B H is the interpolation function matrix; u 0 and q e are expressed as follows: The total strain of this plate in the case of the plate subjected to in-plane prebuckling stresses can be written as follows: Materials 2019, 12, 1262 By substituting Equation (4) into Equations (13) and (14), the strain field can be obtained as follows: with The constitutive relations are derived from Hooke's law by the following equation: with In this work, the stiffener was assumed to be parallel to the x-axis (See Figure 2). There was no delamilation phenomenon between the stiffener and the plate during the performance of the structure; the stiffener seemed to be a beam, and it just bent in the zx-plane. The displacement field of the x-stiffener can be expressed as follows: The displacement vector of node i for the x-stiffener element is written as follows: The displacements of the stiffener can now be given in the following form: structure; the stiffener seemed to be a beam, and it just bent in the zx-plane. The displacement field of the x-stiffener can be expressed as follows: The displacement vector of node i for the x-stiffener element is written as follows: The displacements of the stiffener can now be given in the following form:  These functions can be obtained by substituting s = s0 into Ni and Hi of the plate. The displacement vector of x-stiffener is interpolated as follows: with The displacement vector of x-stiffener is interpolated as follows: with q xs e = q xs 1e q xs 2e T The strain components of the stiffener are as follows: with By substituting Equation (24) into Equations (27)- (29), the strain field is as follows: with The relationship between stresses and strains obtained from Hooke's law is as follows: with σ xs = σ xs The x-stiffener is considered to place at the lower surface of the plate. The condition of displacement at the contact line is as follows: Using Equations (3) and (19), Equation (34) becomes the following: or in shorter form: The nodal displacement vector of stiffener element is as follows: The elastic strain energy of the stiffened plate is written as follows: or in matrix form: where The geometric strain energy enforced by in-plane prebuckling stresses is then computed by the following: By substituting the geometric strain of the plate and the stiffeners into Equation (49), we get the following: Equation (50) now becomes the following: with For the buckling problem, we get the following equation: where K p , K xs and K Gp , K Gxs are the global stiffness matrix and global geometric stiffness matrix, respectively. d stands for the vector of unknowns. Equation (57) is solved to obtain the buckling load λ b and the buckling mode shape.

Formulation Verification
First, we conducted a comparison of the critical buckling loads for a simply supported FGM plate with the analytical results of Meisam et al. [27] and Bodaghi et al. [28], as shown in Table 1. We considered a square plate with a = b = 1 m, the thickness h = a/100, the material properties E c = 380 × 10 9 Pa, E m = 70 × 10 9 Pa, the Poisson ratio 0.3, and the plate as being under an axial compression at two opposite edges. Next, we compared the buckling coefficient of the clamped plate with one central stiffener with a model, as shown in Figure 3. The plate had geometrical parameters a/b = 1 and thickness of h = a/200; the depth of the stiffener was h s = 10.483h, and the width was b s = h s /2.75. The material properties of the plate and the stiffener were E = 68.7 GPa, ν = 0.3. The buckling coefficient was compared with the results of Mukhopadhyay et al. [29] (semianalytical finite difference method), Peng et al. [30] (mesh-free method), and Rikards et al. [31] (finite element method). This comparison is listed in Table 2.   From Tables 1 and 2, we can see clearly that the result of our work compared with other approaches has a very small error, so the calculation program used in this paper is verified.

Long FGM Plate with One Central Stiffener
In this work, two types of FGM plates (Si 3 N 4 /SUS304 plate, ZrO 2 /SUS304 plate) were employed. The material properties are presented in Table 3. We considered a rectangular plate with one central stiffener with b = 0.2 m, a/b = 2.5. Here, the stiffener was made of the same material as the material of the plate. The buckling coefficient is calculated as follows: where E 0 = 5 × 10 7 Pa, a 0 = 10h 0 .

-Effect of the Depth of Stiffener (h s )
We evaluated how the depth of the stiffener affects the buckling behavior of the stiffened FGM plate. The volume fraction index n had a value in the range from 0.1 to 10. Considering three cases of plates with b/h = 100, 150, and 200, the geometrical parameters of the stiffener were b s = 2h, h s = h − 4h; the plate was fully simply supported (SSSS). We also determined the mass reduction of the stiffened structure compared with the plate without stiffeners with the same dimensions (the length and the width of this plate). For the stiffened plate with the thickness h = b/100, the buckling coefficients are represented by the horizontal dash-dot lines in the Figure 4. For each of the plates without stiffeners, we let the thickness of the plate change from b/95 to b/70. The buckling coefficients are represented by the solid lines in Figure 4 and are listed in Table 4. The thickness relationship between the Si 3 N 4 /SUS304 plate with and without a central stiffener is listed in Tables 5 and 6. The red point, which is the intersection point between the dash-dot line and the solid line, represents the equivalent buckling coefficient of the plate without stiffeners and the stiffened plate. Following the same procedure for the stiffened plate with the thickness h = b/150 (the thickness of the unstiffened plate changed from b/140 to b/110) and h = b/200 (the thickness of the unstiffened plate changed from b/180 to h = b/155), we obtained the results given in Figure 5 and Tables 7-16. Si3N4/SUS304 plate with and without a central stiffener is listed in Tables 5 and 6. The red point, which is the intersection point between the dash-dot line and the solid line, represents the equivalent buckling coefficient of the plate without stiffeners and the stiffened plate. Following the same procedure for the stiffened plate with the thickness h = b/150 (the thickness of the unstiffened plate changed from b/140 to b/110) and h = b/200 (the thickness of the unstiffened plate changed from b/180 to h = b/155), we obtained the results given in Figure 5 and Tables 7-16.            These numerical results are very interesting as it obviously demonstrates that the addition of stiffener still ensured the stability of the structure under a compressive load and that it could significantly reduce the mass of the structure due to the thickness of the plate decreasing. This point is really significant for many cases in practical engineering, where the mass of the structure must be reduced to make sure the process works under compressive loads, especially in aerospace engineering, shipbuilding and so on.

-Effect of the gradient index (n)
Next, the authors evaluated the influence of the gradient index n on the buckling response of the plate. We considered two cases of plates made of Si 3 N 4 /SUS304 or ZrO 2 /SUS304 with geometrical parameters h s = 2h, b s = b/100 and b/h ratio value of 10, 20, 50, and 100; the plate was fully simply supported (SSSS). By changing the volume fraction index n from 0.1 to 10, we had a variation of buckling coefficients of FGM plates with one central stiffener and without stiffeners, as shown in Figure 6.
Note that in Figure 6, the solid lines and dash-dot lines represent buckling coefficients of stiffened and unstiffened plates, respectively. The solid lines are always higher than the dash-dot lines, meaning that the buckling load of the stiffened plate was larger than that of the unstiffened plate. This proves that the stiffness of the plate will increase by the addition of the stiffener. We can see that when the thickness of the plate increased with a central stiffener, the plate became stiffer, leading to an enhancing of the buckling load of the plate. When the volume fraction index n was increased, the buckling load of the ZrO 2 /SUS304 plate got larger; in contrast, the buckling load of the Si 3 N 4 /SUS304 plate got smaller. Figure 7 presents the first four buckling mode shapes of the FGM plate with and without a central stiffener. The figures show clearly that the stiffener had a strong influence on the buckling mode shape of the plate.
lines, meaning that the buckling load of the stiffened plate was larger than that of the unstiffened plate. This proves that the stiffness of the plate will increase by the addition of the stiffener. We can see that when the thickness of the plate increased with a central stiffener, the plate became stiffer, leading to an enhancing of the buckling load of the plate. When the volume fraction index n was increased, the buckling load of the ZrO2/SUS304 plate got larger; in contrast, the buckling load of the Si3N4/SUS304 plate got smaller.

-Effect of the length-to-width ratio (a/b)
Next, we examined the effects of the length-to-width ratio (a/b) on the buckling load of the plate with one central stiffener. In this case, we constantly maintained the width of the plate as b = 0.2 m and the width of the stiffener as bs = b/100; the plate was under fully clamped (CCCC). By changing the length of the plate a, the a/b ratio reached the value from 2 to 4. We then obtained the diagram of the buckling load, as shown in Figure 8. The figure shows that when the length-to-width ratio (a/b) increased, the plate became "softer", meaning the buckling load decreased. -Effect of the length-to-width ratio (a/b) Next, we examined the effects of the length-to-width ratio (a/b) on the buckling load of the plate with one central stiffener. In this case, we constantly maintained the width of the plate as b = 0.2 m and the width of the stiffener as b s = b/100; the plate was under fully clamped (CCCC). By changing the length of the plate a, the a/b ratio reached the value from 2 to 4. We then obtained the diagram of

-Effect of the boundary condition
Next, we investigated the influence of the boundary condition on the buckling load of the stiffened FGM plate. Figures 10 and 11 present the buckling coefficient of the structure. The results in Figure 10 were calculated for plates with various b/h ratios, and the results in Figure 11 were calculated for plates with various a/b ratios. From these figures, we can see that the fully clamped supported plate had the largest buckling load, and the fully simply supported plate had the smallest buckling load.

-Effect of the boundary condition
Next, we investigated the influence of the boundary condition on the buckling load of the stiffened FGM plate. Figures 10 and 11 present the buckling coefficient of the structure. The results in Figure 10 were calculated for plates with various b/h ratios, and the results in Figure 11 were calculated for plates with various a/b ratios. From these figures, we can see that the fully clamped supported plate had the largest buckling load, and the fully simply supported plate had the smallest buckling load.

Long FGM Plate with Two Stiffeners
We considered a rectangular plate with two parallel stiffeners, as shown in Figure 12, with b = 0.2 m, a/b = 2.5. The geometrical parameters of the stiffener were h s = 2h, b s = b/100. The buckling coefficient of the plate can be determined by the following formula: with E 0 = 5 × 10 7 Pa, a 0 = 10h 0 .
-Effect of the width of the stiffener and the distance between two stiffeners The effects of the width of the stiffener and the distance between two stiffeners were studied. We compared the buckling loads between the plate with two stiffeners (b s = b/100, b s = b/200, Figure 13b) and the one with one central stiffener (b s = b/100, Figure 13a).
By varying the volume index fraction n and the distance between the two stiffeners d, we obtained the numerical results of the buckling coefficients of the plate, as plotted in Figure 14.  We considered a rectangular plate with two parallel stiffeners, as shown in Figure 12, with b = 0.2 m, a/b = 2.5. The geometrical parameters of the stiffener were hs = 2h, bs = b/100. The buckling coefficient of the plate can be determined by the following formula: with E0 = 5 × 10 7 Pa, a0 = 10h0.

-Effect of the width of the stiffener and the distance between two stiffeners
The effects of the width of the stiffener and the distance between two stiffeners were studied. We compared the buckling loads between the plate with two stiffeners (bs = b/100, bs = b/200, Figure  13b) and the one with one central stiffener (bs = b/100, Figure 13a). By varying the volume index fraction n and the distance between the two stiffeners d, we obtained the numerical results of the buckling coefficients of the plate, as plotted in Figure 14. From Figure 14, we can reach the following conclusions. First, Figure 14 clearly shows that the distance between the two stiffeners strongly influences the buckling coefficient of the structure. In both cases, i.e., Si 3 N 4 /SUS304 and ZrO 2 /SUS304 plates, when d increased, the buckling coefficient of the plate decreased. This is explained by the fact that the strain energy focuses on the central area of the structure. Therefore, when the stiffeners are set there, the stiffness of this plate is larger than other places, thus leading to enhanced buckling coefficient of the plate.
Second, similar to the previous sections, when the volume index fraction n increased, the buckling coefficient of the Si 3 N 4 /SUS304 plate decreased, and the buckling coefficient of the ZrO 2 /SUS304 plate increased. By varying the volume index fraction n and the distance between the two stiffeners d, we obtained the numerical results of the buckling coefficients of the plate, as plotted in Figure 14. From Figure 14, we can reach the following conclusions. First, Figure 14 clearly shows that the distance between the two stiffeners strongly influences the buckling coefficient of the structure. In both cases, i.e., Si3N4/SUS304 and ZrO2/SUS304 plates, when d increased, the buckling coefficient of the plate decreased. This is explained by the fact that the strain energy focuses on the central area of the structure. Therefore, when the stiffeners are set there, the stiffness of this plate is larger than other places, thus leading to enhanced buckling coefficient of the plate.
Second, similar to the previous sections, when the volume index fraction n increased, the buckling coefficient of the Si3N4/SUS304 plate decreased, and the buckling coefficient of the ZrO2/SUS304 plate increased. Third, we explored the different buckling coefficients between a plate with a central stiffener and a plate with two stiffeners. We considered a plate with two stiffeners and a plate with a central stiffener with the same geometrical parameters except the width of the stiffener bs. In case 1, we considered a plate with two stiffeners, which had the total bs of two stiffeners (each stiffener had bs/2) equal to the bs of the plate with a central stiffener (Figure 13a,b). For both computed results of Si3N4/SUS304 and ZrO2/SUS304 plates (Figure 14a,b), the buckling coefficient of the plate with one central stiffener was always larger than the buckling coefficient of the plate with two stiffeners. When the distance d was small, the plate with two stiffeners had the same buckling coefficient as a plate with one central stiffener. In case 2, we considered a plate with two stiffeners, with the width of each stiffener bs equal to bs of the plate with a central stiffener (Figure 13a,c). As can be seen in Figure 14b,d, for both Si3N4/SUS304 and ZrO2/SUS304 plates, when d = 0.5b, the buckling coefficient of the structure was the same as the plate with one stiffener. When the distance of stiffeners was d < 0.5b, the buckling coefficient of the plate with two stiffeners was larger than that of the plate with a central stiffener. When d > 0.5b, the buckling coefficient of the plate with two stiffeners was smaller than that of the plate with a central stiffener. These interesting results have a good significance in practical Third, we explored the different buckling coefficients between a plate with a central stiffener and a plate with two stiffeners. We considered a plate with two stiffeners and a plate with a central stiffener with the same geometrical parameters except the width of the stiffener b s . In case 1, we considered a plate with two stiffeners, which had the total b s of two stiffeners (each stiffener had b s /2) equal to the b s of the plate with a central stiffener (Figure 13a,b). For both computed results of Si 3 N 4 /SUS304 and ZrO 2 /SUS304 plates (Figure 14a,b), the buckling coefficient of the plate with one central stiffener was always larger than the buckling coefficient of the plate with two stiffeners. When the distance d was small, the plate with two stiffeners had the same buckling coefficient as a plate with one central stiffener. In case 2, we considered a plate with two stiffeners, with the width of each stiffener b s equal to b s of the plate with a central stiffener (Figure 13a,c). As can be seen in Figure 14b,d, for both Si 3 N 4 /SUS304 and ZrO 2 /SUS304 plates, when d = 0.5b, the buckling coefficient of the structure was the same as the plate with one stiffener. When the distance of stiffeners was d < 0.5b, the buckling coefficient of the plate with two stiffeners was larger than that of the plate with a central stiffener. When d > 0.5b, the buckling coefficient of the plate with two stiffeners was smaller than that of the plate with a central stiffener. These interesting results have a good significance in practical engineering. For case 1, all plates had the same volume of the stiffener, but using the plate with one central stiffener was much better due to its buckling load always being higher than that of the plate with two stiffeners. In case 2, the plate with two stiffeners was only better than that with one stiffener in the buckling problem when the distance between the two stiffeners was suitable.

-Effect of the plate thickness
In this study, the thickness of the plate (the plate model is given in Figure 12) had values from b/50 to b/5, and the numerical results are shown in Figure 15. We can again see the same phenomenon as in the case of the plate with a central stiffener. When the volume fraction index n increased, the buckling coefficient of Si 3 N 4 /SUS304 plate increased, and the buckling coefficient of ZrO 2 /SUS304 plate decreased. The reason for this is that the plate became softer.
buckling coefficient of Si3N4/SUS304 plate increased, and the buckling coefficient of ZrO2/SUS304 plate decreased. The reason for this is that the plate became softer.

-Effect of the length-to-width ratio (a/b)
We next explored the effects of the length-to-width ratio (a/b). The geometrical parameters were b = 0.2 m, h = b/10. We changed the length a and the volume fraction index n, and the results are plotted in Figure 17. As can be seen from the figure, the length-to-width ratio (a/b) had a robust influence on the buckling coefficient of the structure. This point can be explained by the fact that -Effect of the length-to-width ratio (a/b) We next explored the effects of the length-to-width ratio (a/b). The geometrical parameters were b = 0.2 m, h = b/10. We changed the length a and the volume fraction index n, and the results are plotted in Figure 17. As can be seen from the figure, the length-to-width ratio (a/b) had a robust influence on the buckling coefficient of the structure. This point can be explained by the fact that when the a/b ratio decreased, the length a correspondingly increased. This led to the plate becoming "softer", and the buckling coefficient was therefore reduced.   Figure 18 presents the first four buckling mode shapes of the FGM plate for two cases, i.e., a/b = 2 and a/b = 5; the plate is fully simply supported (SSSS). We can see that both the boundary condition and the a/b ratio strongly affected not only the buckling loads but also the buckling mode shapes of the structure. when the a/b ratio decreased, the length a correspondingly increased. This led to the plate becoming "softer", and the buckling coefficient was therefore reduced.  Figure 18 presents the first four buckling mode shapes of the FGM plate for two cases, i.e., a/b = 2 and a/b = 5; the plate is fully simply supported (SSSS). We can see that both the boundary condition and the a/b ratio strongly affected not only the buckling loads but also the buckling mode shapes of the structure.

Conclusions
In this study, a linear finite element formulation based on the Shi shear deformation theory was employed to study the linear mechanical buckling of stiffened functionally graded material plates. Different parameters were examined, including the geometrical and mechanical properties of the plate. Based on the numerical results presented in the above sections, some major conclusions are listed as follows: - The volume index fraction n has a strong influence on the buckling load of the structure. Due to the different type of material, the buckling coefficient of the Si3N4/SUS304 plate decreased as the volume index fraction n increased. In contrast, the buckling coefficient of the ZrO2/SUS304 plate increased as the volume index fraction n increased. - The geometrical parameters of both the plate and the stiffener also strongly affect the buckling coefficient of the structure. For the Si3N4/SUS304 and ZrO2/SUS304 plates with one central stiffener or two stiffeners, the buckling load increased as the thickness of the stiffener increased. - The addition of stiffener to the plate can significantly reduce the total mass of the structure, but it still maintains the stability of the plate with the same buckling load as the plate without stiffeners. The maximum mass reduction could be reached at 20% for the Si3N4/SUS304 plate and 23.26% for the ZrO2/SUS304 plate. This point is very significant in practical engineering when we need to have lighter structures. Moreover, using a plate with two stiffeners was better than using a plate with one central stiffener when we had the right distance (d < 0.5b) between the two stiffeners. -Boundary conditions also deeply affect the buckling load of the structure. The plate with CCCC boundary condition had much more buckling load than the others.

Conclusions
In this study, a linear finite element formulation based on the Shi shear deformation theory was employed to study the linear mechanical buckling of stiffened functionally graded material plates. Different parameters were examined, including the geometrical and mechanical properties of the plate. Based on the numerical results presented in the above sections, some major conclusions are listed as follows: - The volume index fraction n has a strong influence on the buckling load of the structure. Due to the different type of material, the buckling coefficient of the Si 3 N 4 /SUS304 plate decreased as the volume index fraction n increased. In contrast, the buckling coefficient of the ZrO 2 /SUS304 plate increased as the volume index fraction n increased. - The geometrical parameters of both the plate and the stiffener also strongly affect the buckling coefficient of the structure. For the Si 3 N 4 /SUS304 and ZrO 2 /SUS304 plates with one central stiffener or two stiffeners, the buckling load increased as the thickness of the stiffener increased. - The addition of stiffener to the plate can significantly reduce the total mass of the structure, but it still maintains the stability of the plate with the same buckling load as the plate without stiffeners. The maximum mass reduction could be reached at 20% for the Si 3 N 4 /SUS304 plate and 23.26% for the ZrO 2 /SUS304 plate. This point is very significant in practical engineering when we need to have lighter structures. Moreover, using a plate with two stiffeners was better than using a plate with one central stiffener when we had the right distance (d < 0.5b) between the two stiffeners. -Boundary conditions also deeply affect the buckling load of the structure. The plate with CCCC boundary condition had much more buckling load than the others.