Application of Variational Method to Stability Analysis of Cantilever Vertical Plates with Bimodular Effect

In the design of cantilevered balconies of buildings, many stability problems exist concerning vertical plates, in which reaching a critical load plays an important role during the stability analysis of the plate. At the same time, the concrete forming vertical plate, as a typical brittle material, has larger compressive strength but lower tensile strength, which means the tensile and compression properties of concrete are different. However, due to the complexities of such analyses, this difference has not been considered. In this study, the variational method is used to analyze stability problems of cantilever vertical plates with bimodular effect, in which different loading conditions and plate shapes are also taken into account. For the effective implementation of a variational method, the bending strain energy based on bimodular theory is established first, and critical loads of four stability problems are obtained. The results indicate that the bimodular effect, as well as different loading types and plate shapes, have influences on the final critical loads, resulting in varying degrees of buckling. In particular, if the average value of the tensile modulus and compressive modulus remain unchanged, the introduction of the bimodular effect will weaken, to some extent, the bending stiffness of the plate. Among the four stability problems, a rectangular plate with its top and bottom loaded is most likely to buckle; next is a rectangular plate with its top loaded, followed by a triangular plate with its bottom loaded. A rectangular plate with its bottom loaded is least likely to buckle. This work may serve as a theoretical reference for the refined analysis of vertical plates. Plates are made of concrete or similar material whose bimodular effect is relatively obvious and cannot be ignored arbitrarily; otherwise the greater inaccuracies will be encountered in building designs.


Introduction
In structural engineering, the analysis and design of cantilevered balconies are of huge importance, especially for long-cantilevered balconies. Basically, there are four structural forms for long-cantilevered balconies: (a) O-shaped, (b) L-shaped, (c) C-shaped and (d) U-shaped, as shown in Figure 1, in which h and l are the height and length of the cantilever vertical plate, respectively. In Figure 1, it may be seen that for a cantilever vertical plate, there are two basic shapes: rectangular, as shown in (a) and (c), and triangular, as shown in (b) and (d). In the four structural forms, the rectangular or triangular plate serves as an important structural element in resisting external loads from the top and bottom of the balcony; therefore, it is a key load-bearing component in cantilevered balconies. Among static, dynamic and stability analyses for cantilever vertical plates, the stability problem is particularly prominent because the vertical plate must reach a certain height if it is to meet the design and functional requirements of the building's architecture; compressive stress in the plate becomes inevitable, inducing elastic buckling due to the relatively low plate thickness. On the other hand, most parts constituting the balcony plate are made of concrete, compressive stress in the plate becomes inevitable, inducing elastic buckling due to the relatively low plate thickness. On the other hand, most parts constituting the balcony plate are made of concrete, which has greater compressive strength but lower tensile strength. This disparity between tensile and compressive strength also seems to mean that the tensile and compressive properties of concrete are different. This creates an issue of concern, i.e., what are the effects of different types and degrees of tension and compression? Therefore, this paper provides a stability analysis of cantilever vertical plates from a long-cantilevered balcony in which the applied material has different properties in terms of tension and compression. We will first discuss the elastic buckling problems of plates and the corresponding solutions, before addressing the bimodular problem, as discussed in existing studies. The elastic buckling problem associated with plates has been extensively investigated over the last century. All kinds of plate shapes, different boundary conditions and applied in-plane force distributions have been considered, and buckling critical loads are documented in some stability handbooks [1], standard texts on plates [2][3][4] and in a large number of technical papers (for example, [5,6]). Under different boundary constraints, Wang et al. [7] investigated the elastic buckling problem of vertical plates under body forces (self-weight) in which the extensive buckling body force parameters presented were shown to be useful to engineers for designing structures or machines with vertical plate components. Wang [8] studied the buckling of a standing plate subjected to self-weight and top load, in which an initial value method which was accurate, and more efficient than previous methods, was used. His results indicated that the effects of self-weight and top load both contribute to buckling, although their contributions are not linear. Notably, for narrow plates, if top load dominates, it does not matter whether an end is clamped or simply supported, while if self-weight dominates, the buckling of a narrow plate depends only on the bottom condition. The buckling problems of rectangular plates with various boundary conditions under nonuniform, in-plane loads were investigated in [9][10][11]. It was found that the critical buckling load is associated with loading form and constraint conditions. Jędrysiak and Michalak [12] dealt with the stability problem of thin composite plates with a smooth and slow gradation of macroscopic properties. In their study, The elastic buckling problem associated with plates has been extensively investigated over the last century. All kinds of plate shapes, different boundary conditions and applied in-plane force distributions have been considered, and buckling critical loads are documented in some stability handbooks [1], standard texts on plates [2][3][4] and in a large number of technical papers (for example, [5,6]). Under different boundary constraints, Wang et al. [7] investigated the elastic buckling problem of vertical plates under body forces (self-weight) in which the extensive buckling body force parameters presented were shown to be useful to engineers for designing structures or machines with vertical plate components. Wang [8] studied the buckling of a standing plate subjected to self-weight and top load, in which an initial value method which was accurate, and more efficient than previous methods, was used. His results indicated that the effects of self-weight and top load both contribute to buckling, although their contributions are not linear. Notably, for narrow plates, if top load dominates, it does not matter whether an end is clamped or simply supported, while if self-weight dominates, the buckling of a narrow plate depends only on the bottom condition. The buckling problems of rectangular plates with various boundary conditions under nonuniform, in-plane loads were investigated in [9][10][11]. It was found that the critical buckling load is associated with loading form and constraint conditions. Jędrysiak and Michalak [12] dealt with the stability problem of thin composite plates with a smooth and slow gradation of macroscopic properties. In their study, governing equations of the tolerance models for the stability of thin plates with functionally graded structures were derived. Chróścielewski [13] investigated the torsional buckling phenomenon of thin-walled I-beam columns. In their study, the effect of initial deflection on the torsional buckling load of a thin-walled I-beam column was discussed, and the localization of local buckling modes was studied. Magnucka-Blandzi [14] studied the bending and buckling of a circular plate with symmetrically varying mechanical properties. Magnucki and Magnucka-Blandzi [15] generalized the model of sandwich structure, in which analytical studies of bending and buckling problems of rectangular plates were presented. Through an experiment and numerical simulation, Zhang et al. [16] investigated the compressive buckling behavior of a double-sided, laser-welded Al-Li alloy aircraft fuselage panel. Kolakowski and Jankowski [17,18] investigated the effect of the membrane components of transverse forces on the magnitudes of total transverse forces in the nonlinear stability of plate structures, and noted some inconsistencies in nonlinear buckling in first-order and higher-order shear deformation theory. Cheng [19,20] applied a variational method to analyze the lateral instability problem of cantilever rectangular plates, obtaining the critical loads of plates under the action of a concentrated force, uniformly-distributed loads, triangularly-distributed loads and a concentrated couple. His results indicated that the values of the critical loads are associated with the bending stiffness of the plate and the aspect ratio of rectangular plane, as well as with the choice of the corresponding deflection function. In this study, we focus on determining the critical buckling load. Therefore, the variational method is applied.
Among the various methods of solving the stability problem of plates, there are basically three which apply analytical techniques. The first is the triangle series expansion method, which first requires the buckling differential equation of the plate to be expressed in terms of the deflection and internal forces, and then yields the associated internal forces. The use of this method is limited as, in some cases, it is difficult to establish the buckling differential equation. The second is the differential method, which also makes use of the buckling differential equation. The third is the variational method based on the law of energy, which only requires knowledge of the bending strain energy of the plate and the work done by external forces, thus avoiding the need to include the buckling differential equation. Although the variational method can be used to determine the critical load only, and the displacement and internal forces under the critical load remain unknown, it is still believed that the method may serve as an effective way to determine critical loads, which is particularly important in analyses of plate stability, while the displacement and internal forces are generally obtained as secondary elements, especially in the design of plate.
On the other hand, the structural material also plays an important role in the analysis and design of the structural elements. Generally, the structural material forming the vertical plates in cantilevered balconies is concrete, which is regarded as continuous, elastic, homogeneous and isotropic. However, many studies have indicated that most materials, including ceramics, plastics, concrete, graphite, powder metallurgy materials, polymeric materials and some composites exhibit different tensile and compressive strains under the same tension or compression stress [21,22]. Thus, certain materials exhibit different elastic moduli for tension and compression; these materials are referred to as "bimodular" materials [23]. The concrete used in the present analysis, without exception, is a bimodular material. As a brittle material, the concrete forming vertical plates possesses larger compressive strength but lower tensile strength, which means that its tensile and compression properties are different. Up to now, few reports have been published on concrete with a bimodular effect, even though consideration of the bimodular effect of structural concrete is crucial.
Overall, two bimodular material models are widely used in theoretical analyses within the field of engineering. One is the criterion of positive-negative signs in the longitudinal strain of fibers, as proposed by Bert [24]. This model is mainly applicable to orthotropic materials, and is therefore widely used for research on laminated composites [25][26][27][28]. The other model is the criterion of positive-negative signs of principal stress put forward by Ambartsumyan [29], which is mainly applicable to isotropic materials. In structural engineering, the stress state along certain principal directions is a key issue in stress analyses of components like beams and plates, since it is this factor that determines whether the point is under tension or compression. Due to the fact that this bimodular theory defines the constitutive model based on principal directions, and the principal stress is generally obtained as the final result but not as a known condition before solving, this inevitably gives rise to difficulties in terms of describing the stress state of a given point. This model also lacks the ability to describe the experimental results of elastic coefficients in complex states of stress. Analytical solutions are available in a few cases, although they only concern the static, dynamic and thermal problems of beams and plates [30][31][32][33][34]. In some complex problems, it is necessary to resort to a finite element method (FEM) based on an iterative technique [35][36][37][38].
To date, numerous studies have been published concerning the stability problem for plates under different loads. However, for vertical plates, there is significantly less reported research. From the literature we collected, reports concerning vertical plates mainly included the elastic buckling problem of vertical plates under body forces (self-weight) [7], the buckling of a standing plate subjected to self-weight and top load [8], the lateral instability problem of cantilever rectangular plates under the action of concentrated force, uniformly-distributed loads, triangularly-distributed loads and concentrated couple [19], and the buckling of cantilever rectangular plates under symmetrical edge loading [20]. Among these works, there are obviously two aspects that could be improved upon. One is the fact that none of these studies took into account the bimodular effect of the material, which seems to be problematic. If the material possesses a significant bimodular effect, the error due to the absence of this parameter will be significant. Another may come from the plate shape. For example, triangular plates are seldom discussed (N.B. the originality of the present study is not constituted by our inclusion of triangular plates alone). Broadly speaking, investigations of the stability problems of cantilever vertical plates with different loading types and plate shapes, also considering the bimodular effect of the materials, will be helpful for future analyses of the buckling problem of plates from the perspective of design optimization.
In this study, the variational method is used to solve the stability problem of cantilever vertical plates with bimodular effect, in which different loading types and plate shapes which are widely applied in the construction of real balconies are considered. This paper is arranged as follows. In Section 2, four stability problems of plates with bimodular effect are presented. For the purpose of the effective implementation of the variational method based on bimodular theory, in Section 3, the strain potential energy of a bimodular plate is derived and the work done by external loads is presented. In Section 4, the variational method is applied to determine the critical loads in the four stability problems. The four critical loads without bimodular effect are discussed and the bimodular effect on critical loads is analyzed in Section 5. Some important conclusions are drawn in Section 6.

Stability Problem of Bimodular Cantilever Vertical Plates
Given the structural forms of the balconies presented in Figure 1, two rectangular and two triangular types of cantilever vertical plates are analyzed, as shown in Figure 2. The cantilever length of the plate and the height of the plate are denoted by l and h, respectively. In Figure 2, the xoy plane coordinate system is established, and o is the origin of the coordinate system. According to the force transmission characteristic of the components, the loads are transferred from top to bottom in turn, so that all the loads act downward. Among the stability problems shown in Figure 2, Case (a) is typical, since the vertical plate is subjected to loads from the top and bottom simultaneously. Cases (b) and (c) are supplementary problems, intended to serve as comparisons with Case (a). Comparing Figures 1 and 2, it is easy to see that Cases (a), (b) and (c) in Figure 2 correspond to the O-shaped and C-shaped balconies in Figure 1, while Case (d) corresponds to the L-shaped and U-shaped balconies in Figure 1, in which the cantilever vertical plate is triangular, with its bottom loaded only. For simplification and ease of comparison, we regard all loads q as being uniformly distributed, and for Cases (a), the magnitudes of the two loads are the same when the vertical plate is under the combined action of the top and bottom loads. with its bottom loaded only. For simplification and ease of comparison, we regard all loads q as being uniformly distributed, and for Cases (a), the magnitudes of the two loads are the same when the vertical plate is under the combined action of the top and bottom loads. Additionally, the vertical plate is made of concrete which, in the analysis, is treated as continuous, elastic, isotropic and homogenous, but with the bimodular characteristic indicated above. Generally, in reality, the concrete is reinforced by steel bars in order to make it more resistant to tensile force. In the theoretical analysis presented here, however, we only consider plain concrete without reinforcements. Therefore, our results may serve as a theoretical reference for subsequent reinforcements. Alternatively, from the point of view of structural design, the present theoretical analysis is based on plain concrete and any reinforcements are considered as a safety reserve. Additionally, the vertical plate is made of concrete which, in the analysis, is treated as continuous, elastic, isotropic and homogenous, but with the bimodular characteristic indicated above. Generally, in reality, the concrete is reinforced by steel bars in order to make it more resistant to tensile force. In the theoretical analysis presented here, however, we only consider plain concrete without reinforcements. Therefore, our results may serve as a theoretical reference for subsequent reinforcements. Alternatively, from the point of view of structural design, the present theoretical analysis is based on plain concrete and any reinforcements are considered as a safety reserve.

Variational Method Based on Bimodular Materials Theory
The total potential energy of the system is where U and W are the strain potential energy and the potential of external loads, respectively. From the point of view of energy, the critical load can be obtained on condition that the work done by the variation of external loads, δW, is equal to the variation of strain potential energy, δU, when the plate deforms from the in-plane state to the adjacent bending state. According to Rayleigh-Ritz method, when the system reaches the limit of stable equilibrium, the total potential energy is the minimum, that is which may be used to determine the critical load. It should be noted that in the application of this variational method, U and W must be expressed in terms of the z-axis displacement, w, and that w should satisfy all boundary conditions of displacement. As for the stress boundary conditions, w need not be satisfied, but if a part or all of it can be satisfied, the solution accuracy is greatly improved. In addition, we also note that the calculation of U is based on the theory of elasticity from uniform modulus. If the bimodular theory of elasticity is introduced here, some important improvements may be incorporated. In determining W, however, the bimodular effect has no any influence.

Strain Potential Energy of Bimodular Plates
First, let us derive the expression of the increase of strain potential energy U, which may originate from two different aspects. The first is from the bending deformation in the vertical plane; however, this is negligible, since the bending stiffness in the vertical plane is very large. The second is the torsion and transverse bending under the critical load, thus making the plate deviate from its original vertical plane. The potential energy generated by this deformation represents the vast majority of the total potential energy, so we only consider the deformation of torsion and bending when calculating the increase of potential energy. This practice is similar to the buckling problem of the straight bar in compression, in which only the bending strain energy is considered, while the strain energy from axial compression deformation is ignored. Figure 3 shows a rectangular thin plate with different moduli under tension and compression which is being subjected to the vertical downward uniformly-distributed loads, p, and which is thus bending downward and forming the tensile and compressive areas which are bounded by an unknown neutral layer. In Figure 3, the xoy plane is established on the unknown neutral layer and the z axis is vertical. The thickness of the plate is t, while t 1 and t 2 stand for the tensile and compressive thickness, respectively, both of which will be determined later. E + and µ + are the tensile Young's modulus of elasticity and Poisson ratio, respectively; E − and µ − are the compressive modulus and Poisson ratio, as shown in Figure 3. The constraints of the four sides of the plate are variable; for example, four sides may be fixed or simply-supported, two opposite sides may be fixed or simply-supported, or other, mixed constraint modes may be applied, so long as they cause downward bending under uniformly-distributed load, p. Based on the Kirchhoff-Love hypothesis, the small-deflection bending problem of a thin plate may be analyzed by two-dimensional thin plate bending theory, in which the physical equation will adopt the equation in the state of plane stress, neglecting the effects of the strain components εz, γyz and γzx. Meanwhile, the introduction of the bimodular effect of the material has no influence on the geometrical relation based on deformation characteristics, and therefore, the geometrical equation can still be written as Based on the Kirchhoff-Love hypothesis, the small-deflection bending problem of a thin plate may be analyzed by two-dimensional thin plate bending theory, in which the physical equation will adopt the equation in the state of plane stress, neglecting the effects of the strain components ε z , γ yz and γ zx . Meanwhile, the introduction of the bimodular effect of the material has no influence on the geometrical relation based on deformation characteristics, and therefore, the geometrical equation can still be written as where w is the deflection along the z axis, and ε x , ε y and γ xy are the strain components in two-dimensional thin plate bending theory. Note that due to subarea under tension and compression, the stress components in the tensile and compressive areas are different, which may be expressed as σ + x , σ + y , τ + xy for the tensile area and σ − x , σ − y , τ − xy for the compressive area. Thus the physical equation in the state of plane stress in the tensile area (0 ≤ z ≤ t 1 ,) may be given as and in the compressive area (−t 2 ≤ z ≤ 0) as Substituting Equation (3) into Equations (4) and (5), for 0 ≤ z ≤ t 1 , we have and The bending moment per unit length on the section normal to x and y axes is M x and M y , respectively. These values may be computed by integrating the segments of the tensile and compressive thicknesses, as follows: and Additionally, the torsion moment per unit length on the section normal to x axis, M xy , is In the small-deflection bending problem of thin plates under the action of the uniformlydistributed load, p, the equation of equilibrium expressed in terms of the bending moments, M x and M y , and the torsion moment, M xy , gives Substituting Equations (8) to (10) into Equation (11), we have where 4 is the biharmonic operator. If we let D be the bending stiffness of a bimodular plate, then A familiar form may be obtained, such as D 4 w = p.
The strain potential energy U in a two-dimensional thin plate with bimodular effect may be written, by neglecting the effects of strain components ε z , γ yz and γ zx , as Substituting Equation (3) as well as Equations (6) and (7) into Equation (14), and also integrating z, in which 0 ≤ z ≤ t 1 for the first tensile term and −t 2 ≤ z ≤ 0 for the second compressive term, we obtain Lastly, we have where In fact, the total bending stiffness of the bimodular plate is D = D + +D − . Upon returning to the uniform modulus problem, we have E + = E − = E, µ + =µ − =µ, and t 1 = t 2 = t/2; , which satisfies the regression of the solution.
Note that up to now, the tensile section thickness and the compressive one have not been determined. To address this, we use the condition that the normal forces per unit length on the section normal to x and y axes, N x and N y , are zero in order to determine the location of the unknown neutral layer, that is, and Substituting σ + x , σ + y and σ − x , σ − y from Equations (6) and (7) into Equations (18) and (19), respectively, we have and By adding the left end and right end of Equations (20) and (21), we obtain the following relation after simplification Combining the relation t 1 + t 2 = t, we have , which determines the location of the unknown neutral layer. Obviously, when E + = E − = E, µ + = µ − = µ, the relation t 1 = t 2 = t/2, may be easily obtained. Finally, we obtain the expressions of the strain potential energy and the neutral layer of a bimodular plate.

Work Done by External Loads
Next, let us derive the work done by external loads, W. From Figure 2, it is easy to see that the external loads are uniformly-distributed, q. Supposing the uniformly-distributed loads are acting on the top of the plate, that is y = h, the work done by external loads, W, may be given as where ∆ y (x) is the displacement of the action point of load q when the plate undergoes lateral bending. To calculate ∆ y (x), we may refer to Figure 4. In Figure 4, any cross section of the plate represents a derivation along the z axis, that is, the plate thickness direction from AO to MN, where the tangent at point M intersects the y axis at point C, and θ is the angle between lines AC and MC. Note that θ is a small quantity, and thus, we have the relations AC = MC, θ = (∂w/∂y) M and MB = w M , in which subscript M denotes the quantity which will take the coordinate value of point M. Based on the above relations, ∆ y may be computed as uted loads are acting on the top of the plate, that is y = h, the work done by external loads, W, may be given as where Δy(x) is the displacement of the action point of load q when the plate undergoes lateral bending. To calculate Δy(x), we may refer to Figure 4. In Figure 4, any cross section of the plate represents a derivation along the z axis, that is, the plate thickness direction from AO to MN, where the tangent at point M intersects the y axis at point C, and θ is the angle between lines AC and MC. Note that θ is a small quantity, and thus, we have the relations AC = MC, θ = (əw/əy)M and MB = wM, in which Also noting that MC = w M /sinθ= w M /(∂w/∂y) M , by substituting it into Equation (25), we have Lastly, W may be computed as

Application of Variational Method
In this section, we will derive the critical loads for the four stability problems with bimodular effect shown in Figure 2.
The whole procedure may be described as follows: 1. The first step is to select a deflection function with two unknown parameters. The deflection function should satisfy all boundary conditions of displacement.

2.
Substituting the deflection function into Equations (16) and (27), the potential energy of strain, U, and the potential energy of external forces, W, may be computed. 3.
Using Equations (1) and (2), the total potential energy of the system and its variation may be determined, resulting in two homogeneous linear equations with respect to the previously selected two parameters.

4.
Letting the two linear equations have a nonzero solution, the coefficient determinant of the two equations must be zero; thus, another quadratic equation with respect to q may be obtained; 5.
Finally, by solving the quadratic equation of q, its minimum real root will give the formulas of the corresponding critical load.

Bimodular Rectangular Plate with Top and Bottom Loaded
In Figure 5, the length, height and thickness of the plate are l, h and t, respectively; points A and B stand for the bottom and top corners of the right free side of the plate. Note that, as indicated in Section 3.1, due to the introduction of the bimodular effect, the neutral layer of the plate under bending is not located at the geometrical middle plane, and thus, the y axis deviates from the middle plane, as shown in Figure 5.
The following deflection function is selected [19] where f 1 and f 2 are the z-direction deflection of points A and B of the plate, respectively, and are unknown parameters. Note that a difference remains between Equation (28) in this study and the deflection function in [19] due to the different establishment modes of the coordinate system. From Equation (28), it is easy to see that when x = 0, we have w = 0 and ∂w/∂x = 0, which satisfies the boundary conditions of the entire fixed side; at the same time, when x = l, we have w = f 1 + (f 2 − f 1 ) y/h, which also satisfies the displacement conditions of corner points A and B, that is, for point A, when y = 0, we have w = f 1 . For point B, when y = h, w = f 2 may be obtained. Therefore, Equation (28) is rational and may be used to obtain the solution. The following deflection function is selected [19] where f1 and f2 are the z-direction deflection of points A and B of the plate, respectively, and are unknown parameters. Note that a difference remains between Equation (28) in this study and the deflection function in [19] due to the different establishment modes of the coordinate system. From Equation (28), it is easy to see that when x = 0, we have w = 0 and ∂w/∂x = 0, which satisfies the boundary conditions of the entire fixed side; at the same time, when x = l, we have w = f1 + (f2 − f1) y/h, which also satisfies the displacement conditions of corner points A and B, that is, for point A, when y = 0, we have w = f1 For point B, when y = h, w = f2 may be obtained. Therefore, Equation (28) is rational and may be used to obtain the solution.
Substituting Equation (28) into Equation (16) and integrating with respect to y, 0 to h, and with respect to x, 0 to l, we have The work done by the external loads is divided into two parts: one is the uniformlydistributed load acting on the top of the plate, and the other is that acting on the bottom of the plate. For the work done by the former load, Wt, by substituting Equation (28) into Equation (27), we obtain For the work done by the bottom load, Wb, The total work done by the external loads is Substituting Equation (28) into Equation (16) and integrating with respect to y, 0 to h, and with respect to x, 0 to l, we have The work done by the external loads is divided into two parts: one is the uniformlydistributed load acting on the top of the plate, and the other is that acting on the bottom of the plate. For the work done by the former load, W t , by substituting Equation (28) into Equation (27), we obtain For the work done by the bottom load, W b , The total work done by the external loads is Substituting Equations (29) and (32) into Equation (1), we determine the total potential energy of the system Variations of the two parameters f 1 and f 2 give the following two homogeneous linear equations for f 1 and f 2 and To allow the system of Equations (34) and (35) to yield a nonzero solution, the coefficient determinant must be zero. Thus, we have which gives the minimum real root, that is, the critical load for a bimodular rectangular plate with its top and bottom loaded

Bimodular Rectangular Plate with Top Loaded
In this case, Equation (28) is still selected as the deflection function, in which the meanings of f 1 and f 2 remain unchanged, as shown in Figure 6. Note that here, the work done by the external load is W t in Equation (30). The total potential energy of the system is Similarly, the variations in f1 and f2 give the following two homogeneous linear equations for f1 and f2  The total potential energy of the system is Similarly, the variations in f 1 and f 2 give the following two homogeneous linear equations for f 1 and f 2 and The coefficient determinant of Equations (39) and (40) must be zero; thus which will give the critical load for a bimodular rectangular plate with its top loaded, as follows

Bimodular Rectangular Plate with Bottom Loaded
Equation (28) may still be used, and all other parameters remain the same as those shown in Figure 7. Note that the work done by the external load is W b in Equation (31); thus, we have

Bimodular Rectangular Plate with Bottom Loaded
Equation (28) may still be used, and all other parameters remain the same as those shown in Figure 7. Note that the work done by the external load is Wb in Equation (31)  Similarly, the variations of f1 and f2 give the following two homogeneous linear equations for f1 and f2: and Similarly, the variations of f 1 and f 2 give the following two homogeneous linear equations for f 1 and f 2 : The coefficient determinant of Equations (44) and (45) must be zero; thus which will give the critical load for a bimodular rectangular plate with its bottom loaded as follows:

Bimodular Triangular Plate with Bottom Loaded
In this case, the following deflection function is selected where f 1 and f 2 are the deflection of middle point A of the right angular edge and middle point B at oblique edge of the plate, respectively, as shown in Figure 8. Other quantities are the same as those for the rectangular plate. Note that there is a slight difference between Equation (48) and Equation (28)

Bimodular Triangular Plate with Bottom Loaded
In this case, the following deflection function is selected where f1 and f2 are the deflection of middle point A of the right angular edge and middle point B at oblique edge of the plate, respectively, as shown in Figure 8. Other quantities are the same as those for the rectangular plate. Note that there is a slight difference between Equation (48) and Equation (28)   Substituting Equation (48) into Equation (16) (noting that the integration area is triangular), and thus integrating with respect to y, 0 to h, and with respect to x, 0 to l(1 − y/h), we have Substituting Equation (48) into Equation (16) (noting that the integration area is triangular), and thus integrating with respect to y, 0 to h, and with respect to x, 0 to l(1 − y/h), we have Substituting Equation (48) into Equation (27), the work done by the external load is Substituting Equations (49) and (50) into Equation (1), the total potential energy of the system is The variations of f 1 and f 2 give the following two homogeneous linear equations for f 1 and f 2 The coefficient determinant of Equations (52) and (53) must be zero; thus giving the critical load for a bimodular triangular plate with its bottom loaded, as follows

Comparision of Four Critical Loads without Bimodular Effect
First, let us compare the critical loads of four stability problems without the bimodular effect, which are denoted by q cr1 *, q cr2 *, q cr3 * and q cr4 *. For this purpose, let D + = D − = D/2 and µ + = µ − = µ in Equations (37), (42), (47) and (55). Thus, for a rectangular plate with its top and bottom loaded, we have  Due to the fact that the magnitude of the critical load is closely related to tendency of the plate to buckle, the following conclusions may be drawn from Figure 9: (1) Among the four cases, the order is, by the degree of ease of buckling, qcr1* < qcr2* < qcr4* < qcr3*; (2) For the rectangular plates, case (a) with top and bottom loaded is most likely to buckle, followed by case (b) with top loaded, and then case (c) with bottom loaded; (3) By comparing Due to the fact that the magnitude of the critical load is closely related to tendency of the plate to buckle, the following conclusions may be drawn from Figure 9: (1) Among the four cases, the order is, by the degree of ease of buckling, q cr1 * < q cr2 * < q cr4 * < q cr3 *; (2) For the rectangular plates, case (a) with top and bottom loaded is most likely to buckle, followed by case (b) with top loaded, and then case (c) with bottom loaded; (3) By comparing (c) with (d), it is easy to see that, under the same load, the triangular plate is more prone to buckling than the rectangular plate.
Cheng [19] also determined the critical load of a cantilever vertical plate with its top loaded without considering the bimodular effect. His solution was as follows where (q 0 ) cr is the critical load, a and b are the height and length of the cantilever vertical plate, respectively, and D is the bending stiffness of the plate. Equation (57) in this study and Equation (60) in [19] are the same. This verifies, to some extent, the correctness of the analytical solutions derived in this study; however, in our study, the bimodular effect is introduced and other loading types and plate shapes are also considered. Of course, a demonstration of validity should include a numerical simulation based on FEM, in which a real material model may be incorporated. Nonetheless, as indicated in the Introduction, for a bimodular problem, it is difficult and time-consuming to use the finite element method (FEM) based on an iterative technique. In addition, there is no bimodular materials model in existing commercial FEM software, e.g., ABAQUS. New material types may be introduced through the use of external modules, for example, by applying the user materials subroutine, UMAT, to evaluate the mechanical behavior of a given type of material. Given that the main aim of this study is to use the variation method to determine the critical loads of various components, more accurate numerical simulations will be considered in future.

Bimodular Effect on Four Critical Loads
For the investigation of the bimodular effect on the critical loads, let us introduce the following relations [39] where E is the average value of the tensile and compressive moduli, and µ is the average value of the tensile and compressive Poisson ratio. β is an important parameter which may be positive or negative. In addition, we introduce the following dimensionless quantities Using the relations in Equations (61) and (62), we may obtain the dimensionless forms of Equations (37), (42), (47) and (55). For the rectangular plate with its top and bottom loaded, we have For the rectangular plate with its top loaded, we have For the rectangular plate with its bottom loaded, we have For the triangular plate with its bottom loaded, we have Up to now, we have determined the dimensionless critical loads which may be numerical determined via known parameters such as λ h , λ t , µ + , µ − , K + and K − , in which the T 1 and T 2 contained in K + and K − are, based on Equation (23) Table 1 lists the value of β taken from Equation (61), ranging from 0.3 to −0.3, with an interval of 0.05. Thus, there are, in total, thirteen data groups. Table 1 also lists the corresponding computational value of related physical quantities including µ + and µ − from Equation (61), in which the average modulus, µ, is first set to 0.25, and also includes T 1 and T 2 from Equation (67), as well as K + and K − from Equation (62). β = 0 indicates that the tensile modulus is equal to the compressive one, and also that, in this case, we have µ + = µ − , T 1 = T 2 and K + = K − . From the data of Table 1, it is easy to observe a regular pattern, indicating that when the abstract values of β remain unchanged, only a change of positive or negative sign occurs for µ + and µ − , T 1 and T 2 , and K + and K − , which exchange their values. Specifically, if β changes from 0.05 to −0.05, µ + = 0.2625 when β = 0.05 is exactly µ − = 0.2625 when β = −0.05, while at the same time, µ − = 0.2375 when β = 0.05 is exactly µ + = 0.2375 when β = −0.05. A similar pattern may be found for T 1 and T 2 and K + and K − . This fact reminds us that if these quantities are symmetrically distributed in the expressions for the critical load (the so-called symmetry may be found in Equations (61) to (66)), as long as the abstract value of β remains unchanged, whether this value is positive or negative, the final calculated value of the critical loads remains unchanged, that is, the relation of Q cri (i = 1,2,3,4) and β satisfies Q cri (−β) = Q cri (β), (i = 1, 2, 3, 4).
(68) It is easy to see from Table 1 that the dimensionless bending stiffness, K= K + + K − , also forms a regular pattern, that is These results may be helpful to the following analysis and discussions. Tables 2-8 list four critical loads under different λ h and λ t when β = 0, ±0.05, ±0.1, ±0.15, ±0.2, ±0.25 and ±0.3, and λ h ranges from 0.9 to 1.1 and λ t from 0.02 to 0.04, both according to real conditions in the cantilever vertical plates used for the design of balconies. Specifically, one designer may refer to the dimensionless value of critical loads to determine the real critical loads without much effort via the expression q cr = EtQ cr , where E is the average modulus and t is the thickness of the plate.
From Tables 2-8, it is easy to find the influences of parameters λ h and λ t on the critical loads. As indicated in Equation (62), λ h = h/l, ranging from 0.9 to 1.1. When λ h is 0.9, 1.0 and 1.1, this corresponds to the cases of h < l, h = l (square plate) and h > l. It may be observed that, according to the magnitudes of the critical loads, h < l is most likely to buckle, followed by the square plate, while h > l is least likely to buckle. For the thickness variation of λ t = t/l, it is obvious that the increase in thickness improves instability; thus, λ t = 0.02 is most likely to buckle, followed by λ t = 0.03, while λ t = 0.04 is the least likely to buckle.
From Tables 2-8, we can find the basic tendency is still, according to the degree of buckling, Q cr1 > Q cr2 > Q cr4 > Q cr3 , indicating that regardless of the values of β, λ h and λ t , a rectangular plate with its top and bottom loaded is most likely to buckle, followed by the rectangular plate only with its top loaded, the triangular plate with its bottom loaded, and finally, the rectangular plate only with its bottom loaded. The difference magnitudes among the critical loads of four groups are also different. From Tables 2-8, it is easy to see that there are smaller differences between Q cr1 and Q cr2 , as well as slightly bigger differences between Q cr3 and Q cr4 , while larger differences may be observed between Q cr1 (or Q cr2 ) and Q cr3 (or Q cr4 ).  Figure 10 shows the variation curves of four critical loads with β values. Although it is plotted for a given case, for example, λ h = 1.0 and λ t = 0.03, it is applicable to other cases of values of λ h and λ t . It is easy to see from Figure 10 that Q cr1 > Q cr2 > Q cr4 > Q cr3 , as well as the relative differences among the four critical loads, that is, Q cr1 is very close to Q cr2 , Q cr3 is relatively close to Q cr4 , but there is larger difference between Q cr1 (or Q cr2 ) and Q cr3 (or Q cr4 ). In addition, the two curves of Q cr3 and Q cr4 have obvious arch shapes, indicating that β has a greater influence on Q cr3 and Q cr4 , while the two curves of Q cr1 and Q cr2 approach a flat line, indicating a smaller influence of β on Q cr1 and Q cr2 . It may concluded therefore that the introduction of bimodular effect has a greater influence on rectangular and triangular cantilever vertical plates with their bottoms loaded, and a smaller influence on rectangular plates with either their top and bottom loaded or only the top loaded. variation of λt = t/l, it is obvious that the increase in thickness improves instability; thus, λt = 0.02 is most likely to buckle, followed by λt = 0.03, while λt = 0.04 is the least likely to buckle.
From Tables 2-8, we can find the basic tendency is still, according to the degree of buckling, Qcr1 > Qcr2 > Qcr4 > Qcr3, indicating that regardless of the values of β, λh and λt, a rectangular plate with its top and bottom loaded is most likely to buckle, followed by the rectangular plate only with its top loaded, the triangular plate with its bottom loaded, and finally, the rectangular plate only with its bottom loaded. The difference magnitudes among the critical loads of four groups are also different. From Tables 2-8, it is easy to see that there are smaller differences between Qcr1 and Qcr2, as well as slightly bigger differences between Qcr3 and Qcr4, while larger differences may be observed between Qcr1 (or Qcr2) and Qcr3 (or Qcr4). Figure 10 shows the variation curves of four critical loads with β values. Although it is plotted for a given case, for example, λh = 1.0 and λt = 0.03, it is applicable to other cases of values of λh and λt. It is easy to see from Figure 10 that Qcr1 > Qcr2 > Qcr4 > Qcr3, as well as the relative differences among the four critical loads, that is, Qcr1 is very close to Qcr2, Qcr3 is relatively close to Qcr4, but there is larger difference between Qcr1 (or Qcr2) and Qcr3 (or Qcr4). In addition, the two curves of Qcr3 and Qcr4 have obvious arch shapes, indicating that β has a greater influence on Qcr3 and Qcr4, while the two curves of Qcr1 and Qcr2 approach a flat line, indicating a smaller influence of β on Qcr1 and Qcr2. It may concluded therefore that the introduction of bimodular effect has a greater influence on rectangular and triangular cantilever vertical plates with their bottoms loaded, and a smaller influence on rectangular plates with either their top and bottom loaded or only the top loaded.  Figure 11 shows the variation of dimensionless bending stiffness K with β, suggesting the same arch shape as the critical loads Qcr, which indicates that the variation of Qcr with β is similar to that of K with β. From Equations (56-59), it is easy to see that qcr* is directly proportional to the bending stiffness D in a uniform modulus case, but what happens if the bimodular effect is introduced? From Equations (63)-(66), it is found that Qcr is not proportional to the bending stiffness K in a bimodular case; in contrast to Equations (56)-(59), the relation between Qcr and K seems to be nonlinear. However, by considering Equations (63)-(66), it may be seen that if we neglect the difference of μ + and μ − , and further take μ + = μ − = μ, then the last item in the square root will become μ. Thus, the resulting Qcr is also proportional to the bending stiffness K in a bimodular problem. This conclusion is rational, since the influence of Poisson ratio μ is generally regarded as secondary, compared with the influence of the modulus of elasticity. It is for this reason that, at present, most studies do not take into account the influence of Poisson ratio on stress and deformation; of course, this also includes the critical loads discussed here.  Figure 11 shows the variation of dimensionless bending stiffness K with β, suggesting the same arch shape as the critical loads Q cr , which indicates that the variation of Q cr with β is similar to that of K with β. From Equations (56)-(59), it is easy to see that q cr * is directly proportional to the bending stiffness D in a uniform modulus case, but what happens if the bimodular effect is introduced? From Equations (63)-(66), it is found that Q cr is not proportional to the bending stiffness K in a bimodular case; in contrast to Equations (56)-(59), the relation between Q cr and K seems to be nonlinear. However, by considering Equations (63)-(66), it may be seen that if we neglect the difference of µ + and µ − , and further take µ + = µ − = µ, then the last item in the square root will become µ. Thus, the resulting Q cr is also proportional to the bending stiffness K in a bimodular problem. This conclusion is rational, since the influence of Poisson ratio µ is generally regarded as secondary, compared with the influence of the modulus of elasticity. It is for this reason that, at present, most studies do not take into account the influence of Poisson ratio on stress and deformation; of course, this also includes the critical loads discussed here.