Non-Linear Stability of the Step-Variable In-Plane Functionally Graded Plates Subjected to Linear Approaches of the Edges

The analysis of gradations through the thickness in structures are commonly used. It usually refers to the problems of the stability of functionally graded (FG) structures. In this work, rectangular in-plane FG plates built of a material gradation along the transversal direction were assumed. Five-strip FG plates with four cases that were based on the boundary conditions on longitudinal edges and simply supported on transverse loaded edges were considered. The non-linear stability problems of the FG plates that were subjected to linear approaches of the transverse edges for several types of loads were solved. The estimations were executed with two methods: an analytical-numerical way based on Koiter’s theory and finite element method (FEM).


Introduction
Since the 1980s, Functionally Graded Materials (FGMs) have been treated as a new class of composites mainly built of two components. In most cases, there are metallic constituents and ceramic constituents. FGM structures are commonly used in thin-walled construction elements thanks to high thermal properties. A pretty extensive review of theories related FGMs was included in references [1][2][3]. Additionally, in [4], a literature review concerning FG plate structures and analysis of the post-buckling behaviour for a step-variable FG box across the plate thickness is included. At present, different techniques of producing FGMs can be distinguished [5].
Owing to that fact, many analytical, numerical, and experimental studies of FGM structures have been conducted to better understand their behaviour with respect to nonlinear stability and limitations of carrying-load capacity in references [6][7][8]. The analysis of vibrations of FGM structures was performed in references [9,10], among others. In works [11,12], constitutive relations referred to neutral physical surface for functionally graded materials based on nonlinear six-parameter shell theory. A response of FGM shell structures under thermal load within the first-order shear deformation theory (FSDT) was investigated in reference [13]. The behaviour of the FGM plate subjected to in-plane compressive or thermal load and combination of both loads, on the basis of the classic plate theory (CPT), was studied in [14]. In reference [15], an analysis of free vibrations of printed polymer FG plate by employing Higher Order Shear Deformation (HSDT) with consideration of changing in the stiffness and density of the material was presented. Analyses of FGM structures applying isogeometric analysis (IGA) were presented in [16,17].
Non-linear stability of step-variable in-plane functionally graded (FG) square plates is presented in [18]. In this work, the stability of plates that were subjected to a shear and compression were

Description of the Problem
The materials of in-plane FGM plate strips obeying Hooke's law were assumed. For the plates, nonlinear geometrical strain relations were assumed to take the full bending of the plate into account [19,20]: and where: u, v, w-components of the displacement vector of the panel along the x, y, z axis direction, respectively; ε x , ε y , ε xy -in-plane strains (full Green's strain tensor); κ x , κ y , κ xy -plate rotations (i.e., the classical laminate plate theory; CLPT). The issue of the nonlinear stability of thin structures was estimated with Koiter's theory [20][21][22]. Based on all variants of Koiter's approach [23], the most used way is the method that was developed by Byskov and Hutchinson [20][21][22] formulated in a convenient way.
Regarding the solution of the problem, BHT in a range of the first and second order approximation was applied. This procedure of estimations is widely described in references [18,20]. For the thin-walled structures with imperfections (defined as a mode of linear buckling), the equations for uncoupled buckling are written in the form [20][21][22]: where: ε is the applied strain, ε cr the buckling strain of the buckling mode, ζ * the dimensionless amplitude as initial deflection referred to the buckling mode, and finally ζ the dimensionless amplitude of the deflection. The coefficient b 1111 can be defined based on the equations given in references [20][21][22].

Analysis of the Calculations Results
The numerical computations for rectangular plates were carried out. The plate that was composed of five strips with the same width along the transverse direction was taken into account ( Figure 1a). Table 1 sorts out the material properties for each strip. Different stiffness determined each strip. The assumed geometric dimensions of the plate and a sequence of the strip composition are illustrated in Figure 1. Five types of loading were assumed to be analysed. A linear displacement of transverse edges was assumed (as shown in Figure 1). Type 1 refers to uniform displacements of the loaded edges taken into consideration in [18]. Two next types of loading concern a displacement of edges varying triangularly, however type 2 corresponds to the maximal strain on the ceramic side Al 2 O 3 , whereas type 3 relates to the maximal strain on the aluminium side Al. The strains are equal to 1 in other points of the plate corners ( Figure 1b).

转公式直接闪退：
Home-Styles-Manage Styles-Import/Export-全选左边-Delete-全点 OK-退出右下角 MathType Server-注销再进电脑 平时不使用 Mendeley Desktop：Tools 第三个 C:\Users\MDPI\AppData\Roaming\Notepad++\shortcuts body find 录制宏 红色录制 搜索相关标签 停止录制 保存 The two last types refer to loading with differ directions of strains, nevertheless the absolute values of the strains on the opposite edges are the same. In the case of type 4, the strip Al 2 O 3 is stretched, but strip Al is compressed. For type 5, it is reverse (see Figure 1).
In addition, four cases of boundary conditions were assumed: simply supported plate on all of its edges and a mixture of simply supported and clamped edges. With regard to the plates that are built of strips with different stiffness, the following notation for a sequence of the boundary conditions was introduced: line (1) for x = 0, line (2) for y = 0, line (3) for x = L, and line (4) for y = b. The following Four cases of the boundary conditions were taken into account for the study of non-linear stability, namely: • Case A -SSSS: • Case B -SCSC: • Case C -SSSC: The stress and moment resultants (N x , N y , N xy , M x , M y , M xy ) are defined while using CLPT, as in references [20,21], among others.
The buckling strains were determined to compare the post-buckling equilibrium paths for each case. For five cases of load ( Figure 1b), a five-strip rectangular FG plate of the dimensions ( Figure 1a): b = 200mm, plate thickness t = 2mm, and each strip b1 = b2 = b3 = b4 = b5 = 0.2b were assumed. The length of the plate L was different for the considered boundary conditions and determined based on minimal value of bifurcation.
The present results were obtained by two methods: • BHT by using the first-and second-order non-linear approximation-more details can be found in references [18,20]; and, • FEM with the commercial ANSYS ® software [24]-a discrete model is illustrated in Figure 2. elements with high stiffness on the loaded edges were taken into consideration on the loaded edges to satisfy type 4 and type 5 of loads. The assumption of the beam elements on these edges does not influence the total bending stiffness of the plate, because all of edges of the plate were constrained in the perpendicular direction with respect to the in-plane of the plate. The parameter  that is given in Figure 2 denotes a displacement of the point or the edge.

Analytical-Numerical Method Using BHT
This subsection presents the results that were obtained within BHT. Table 2 lists the values of bifurcational strains for the considered cases and types of load. For the cross-section area of the plate for const x  (Figure 1), a displacement of the gravity centre from axis mm b y 100 2 /   due to different density of individual strips (Table 1) occurs. The coordinate of the gravity centre in the crosssection area of the plate is mm y c 41 95.  . The bifurcational load acting in the plate was reduced to this central point. In Table 2, the values of the buckling forces cr P and the buckling bending moments cr M acting in the plate plane were given. A positive value of the moment is assumed when the ceramic strip is compressed, whereas the metal strip stretched. The element was considered to model a strip plate [24]. The element is characterised by eight nodes with six degrees of freedom in each node. The used element is pretty good for linear and non-linear simulation for both isotropic and multilayers structures. This element is based on Mindlin-Reissner's theory (the first-order deformation theory). Figure 2 shows a discrete model with the conditions of simulations. The size of the finite element amounted to 2 mm long. The numerical calculations were carried out for large displacements with Green-Lagrangian equations. The computations were performed in steps by setting a number of steps from 100 to 500. Beam189 elements with high stiffness on the loaded edges were taken into consideration on the loaded edges to satisfy type 4 and type 5 of loads. The assumption of the beam elements on these edges does not influence the total bending stiffness of the plate, because all of edges of the plate were constrained in the perpendicular direction with respect to the in-plane of the plate. The parameter ∆ that is given in Figure 2 denotes a displacement of the point or the edge.

Analytical-Numerical Method Using BHT
This subsection presents the results that were obtained within BHT. Table 2 lists the values of bifurcational strains for the considered cases and types of load. For the cross-section area of the plate for x = const (Figure 1), a displacement of the gravity centre from axis y = b/2 = 100mm due to different density of individual strips (Table 1) occurs. The coordinate of the gravity centre in the cross-section area of the plate is y c = 95.41mm. The bifurcational load acting in the plate was reduced to this central point. In Table 2, the values of the buckling forces P cr and the buckling bending moments M cr acting in the plate plane were given. A positive value of the moment is assumed when the ceramic strip is compressed, whereas the metal strip stretched.
A negative value P cr denotes the tension force, however a negative value of the moment M cr means an opposite value of the moment. An influence of the boundary conditions on the stability and post-buckling equilibrium paths was compared for each of five types of load. As is known from the literature, the lowest values of bifurcational loads are for A (SSSS), however they are the greatest ones for B (SCSC). For the case C (SSSC) and D (SCSS), the critical loads are intermediate. The strip that is built of the FG plate can introduce some kind of modification, because the ceramic strip is significantly stiffer than the aluminium one. For type 1, the above formulated rules for the buckling strains are concerned. The longest plate L is for the case A but the shortest one is for case B.
In the case C, the bifurcational values are lower than for the case D, but it is opposite in the case of the plate length L. Conclusions are very similar for types 2 and 3. The exception is type 3 for the case C and D. Here, a contrary sequence of bifurcational values occurs. In type 3 for the case D, the stiffest ceramic strip is less loaded and the condition of its fixing does not matter meaningfully as the fixing of the aluminium strip at minor stiffness. Thus, a change in bifurcational values follows due to this reason. The critical values of moments are low and negative, which points out to a change in direction. For type 4, the ceramic strip is elongated, but the aluminium strip is compressed, which decides about the stability loss of the plate for the assumed boundary conditions. The bifurcational values for the cases A and D are the lowest and close to each other, similarly as the length L; however, for cases B and C, the values are the highest and the same refers to the length. The assumed distribution of strains for type 4 explains the negative values of forces and moments that are caused by considerable changes in the stiffness of utmost strips. The negative force means that the plate is stretched, but a negative moment results from a change of its direction of operation. For type 5, the ceramic strip is compressed, which decides about the loss stability, but the aluminium strip is stretched. For cases A and C, the lowest critical values are achieved contrary to the values of the cases B and D, which are the highest and differ slightly from each other. The lengths L for the cases A and C are greater than for the cases B and D. In Figures 3-7, post-buckling paths (Equation (3)) for the FG plates subjected to five types of loads 1-5 are shown. On each of the charts, post-buckling curves for the cases A-D are depicted.       The post-buckling paths for the perfect plates (i.e., for ζ * = 0.0 in Equation (3)) are presented in an untraditional system of coordinates, namely ε/ε cr as a function of the non-dimensional square amplitude of the deflection ζ 2 , which causes the post-buckling paths to be straight lines, according to Equation (3). For types 2 (Figure 4), 3 ( Figure 5), and 5 (Figure 7), the lowest post-buckling equilibrium paths respond to the case C (SSSC), next D (SCSS) and the path for B (SCSC) and the highest path for A (SSSS). For type 1 (Figure 3), the lowest one is for C, next for the tapes A and D, and the highest one for B. In Figure 6, for type 4, the order of the paths is as follows: D, C, B, and A. The type of load and the boundary conditions on the stiffest edge determine the order of the post-buckling paths (i.e., strip 1).  The strip that is built of the FG plate can introduce some kind of modification, because the ceramic strip is significantly stiffer than the aluminium one. For type 1, the above formulated rules for the buckling strains are concerned. The longest plate L is for the case A but the shortest one is for case B. In the case C, the bifurcational values are lower than for the case D, but it is opposite in the case of the plate length L . Conclusions are very similar for types 2 and 3. The exception is type 3 for the case C and D. Here, a contrary sequence of bifurcational values occurs. In type 3 for the case D, the stiffest ceramic strip is less loaded and the condition of its fixing does not matter meaningfully as the fixing of the aluminium strip at minor stiffness. Thus, a change in bifurcational values follows due to this reason. The critical values of moments are low and negative, which points out to a change in direction. For type 4, the ceramic strip is elongated, but the aluminium strip is compressed, which decides about the stability loss of the plate for the assumed boundary conditions. The bifurcational values for the cases A and D are the lowest and close to each other, similarly as the length L ; however, for cases B and C, the values are the highest and the same refers to the length. The assumed distribution of strains for type 4 explains the negative values of forces and moments that are caused  The strip that is built of the FG plate can introduce some kind of modification, because the ceramic strip is significantly stiffer than the aluminium one. For type 1, the above formulated rules for the buckling strains are concerned. The longest plate L is for the case A but the shortest one is for case B. In the case C, the bifurcational values are lower than for the case D, but it is opposite in the case of the plate length L . Conclusions are very similar for types 2 and 3. The exception is type 3 for the case C and D. Here, a contrary sequence of bifurcational values occurs. In type 3 for the case D, the stiffest ceramic strip is less loaded and the condition of its fixing does not matter meaningfully as the fixing of the aluminium strip at minor stiffness. Thus, a change in bifurcational values follows due to this reason. The critical values of moments are low and negative, which points out to a change in direction. For type 4, the ceramic strip is elongated, but the aluminium strip is compressed, which decides about the stability loss of the plate for the assumed boundary conditions. The bifurcational values for the cases A and D are the lowest and close to each other, similarly as the length L ; however, for cases B and C, the values are the highest and the same refers to the length. The assumed distribution of strains for type 4 explains the negative values of forces and moments that are caused

Comparison of the BHT Results with FEM Results
On the basis of the conclusions resulting from [18], which confirms close correlation of the results from the BHT and the FEM for a uniform displacement of the plate edges, the BHT results were verified by FEM in 25 percent, as shown in the last column of Table 2 in the present work. A special

Comparison of the BHT Results with FEM Results
On the basis of the conclusions resulting from [18], which confirms close correlation of the results from the BHT and the FEM for a uniform displacement of the plate edges, the BHT results were verified by FEM in 25 percent, as shown in the last column of Table 2 in the present work. A special attention was paid to types 4 and 5. As it can be easily seen, very close results were received with both the methods. Figures 8-12 present the post-buckling equilibrium paths for the verified cases in a traditional system of coordinates ε/ε cr as a function of ∆w/t for the imperfection ζ * = w 0 /t = 0.1. In general, the paths go through similarly, though some greater discrepancies are visible for the boundary condition A, type 4 ( Figure 9). For the remaining variants, differences in the plate deflections turn out in the vicinity of the buckling load ( Figure 8) or at circa 1.5 · ε/ε cr (Figures 10-12). It should be noted that the BHT results were obtained for the uncoupled buckling mode (the modulation of buckling modes within a load increase was not considered). It can be reason that some small differences between FEM and BHT appear.

Comparison of the BHT Results with FEM Results
On the basis of the conclusions resulting from [18], which confirms close correlation of the results from the BHT and the FEM for a uniform displacement of the plate edges, the BHT results were verified by FEM in 25 percent, as shown in the last column of Table 2 in the present work. A special attention was paid to types 4 and 5. As it can be easily seen, very close results were received with both the methods. Figures 8-12 present the post-buckling equilibrium paths for the verified cases in a traditional system of coordinates for the imperfection . In general, the paths go through similarly, though some greater discrepancies are visible for the boundary condition A, type 4 ( Figure 9). For the remaining variants, differences in the plate deflections turn out in the vicinity of the buckling load ( Figure 8) or at circa ( Figures. 10-12). It should be noted that the BHT results were obtained for the uncoupled buckling mode (the modulation of buckling modes within a load increase was not considered). It can be reason that some small differences between FEM and BHT appear.

Summary
The analysis of non-linear stability of graded plate that was subjected to linear approaches of the transverse loaded edges was carried out. Four cases of the supports and five different loads were assumed. The estimations were conducted mainly based on Koiter's approach with the use of BHT.

Summary
The analysis of non-linear stability of graded plate that was subjected to linear approaches of the transverse loaded edges was carried out. Four cases of the supports and five different loads were assumed. The estimations were conducted mainly based on Koiter's approach with the use of BHT.

Summary
The analysis of non-linear stability of graded plate that was subjected to linear approaches of the transverse loaded edges was carried out. Four cases of the supports and five different loads were assumed. The estimations were conducted mainly based on Koiter's approach with the use of BHT. The second used method was relied on FEM by applying commercial ANSYS ® software. The obtained results with these two methods showed good correlation.
Author Contributions: L.C., Z.K. elaborated this analysis. L.C. performed the numerical simulations. Z.K. conducted the analytical-numerical calculations. Z.K. prepared the first version of this paper. Both authors took a part in revision of the manuscript. All authors have read and agreed to the published version of the manuscript.