Large Deformation Problem of Bimodular Functionally-Graded Thin Circular Plates Subjected to Transversely Uniformly-Distributed Load: Perturbation Solution without Small-Rotation-Angle Assumption

: In this study, the large deformation problem of a functionally-graded thin circular plate subjected to transversely uniformly-distributed load and with different moduli in tension and compression (bimodular property) is theoretically analyzed, in which the small-rotation-angle assumption, commonly used in the classical Föppl–von K á rm á n equations of large deﬂection problems, is abandoned. First, based on the mechanical model on the neutral layer, the bimodular functionally-graded property of materials is modeled as two different exponential functions in the tensile and compressive zones. Thus, the governing equations of the large deformation problem are established and improved, in which the equation of equilibrium is derived without the common small-rotation-angle assumption. Taking the central deﬂection as a perturbation parameter, the perturbation method is used to solve the governing equations, thus the perturbation solutions of deﬂection and stress are obtained under different boundary constraints and the regression of the solution is satisﬁed. Results indicate that the perturbation solutions presented in this study have higher computational accuracy in comparison with the existing perturbation solutions with small-rotation-angle assumption. Specially, the computational accuracies of external load and yield stress are improved by 17.22% and 28.79% at most, respectively, by the numerical examples. In addition, the small-rotation-angle assumption has a great inﬂuence on the yield stress at the center of the bimodular functionally-graded circular plate.


Introduction
The elastic large deformation problem of flexible thin plates has been a focus of attention for scholars all over the world. Because of its strong geometrical nonlinearity, it is generally difficult to establish the effective governing equations used for the solution to the problem. For a long time, the classical Föppl-von Kármán equation has been used to describe the large deflection behavior of thin plates; however, the influence caused by small-rotation-angle assumption commonly used in the derivation of the equation has always been ignored. In some real applications, the deflection magnitude is moderate, while the rotation angle of the plate may achieve a relatively large value, which makes investigation into small-rotation-angle assumption inevitable.
In addition to the geometrical nonlinearity caused by the elastic large deformation problem, the materials that constitute flexible plates also present some nonlinear problems that cannot be solved by the original methods used for elastic, isotropic, and homogeneous materials, for example, functionally graded material (FGM) with different properties in tension and compression. Therefore, the nonlinearity of materials, in combination with the In recent decades, the issues related to FGM have become an important research topic in many engineering and technical fields, such as aerospace [23], civil engineering [24], acoustics [25], and micro-electro-mechanical systems [26]. FGMs are new types of composite material, which are composed of two or more materials, such as alloy and ceramic, and the composition of these two (or more) materials presents continuous gradient changes [27]. There are many works on structural elements made of FGM, most of them dealing with beams and plates, for example, Nguyen et al. [28]. Among the studies, few studies consider the bimodular property of FGM. As indicated above, most materials will present bimodular properties, which can be obvious or not; thus, FGM seems to be no exception.
More recently, the bimodular property of materials has been gradually introduced into the analysis field of FGM and structures. There are also a few works concerning bimodular FGM beam and plate problems [29][30][31]. Focusing on bimodular FGM plates, He et al. [30] established the simplified theory based on a neutral layer under small-deflection; thereafter, He et al. [31] derived the governing equations of the large-deflection problem of bimodular FGM thin circular plate. However, consideration of geometrical nonlinearity seems to be insufficient. As indicated above, for the large deformation problem of thin plates, not only the deflection of the plate but also the rotation angle presents a relatively large value, and the two factors cannot be ignored.
The definition of the small-rotation-angle assumption means the rotation, β, of the deformed thin plate is so small that it satisfies the approximate relationship sinβ ≈ tanβ. In fact, because sinβ is equal to 1/[1 + 1/tan 2 β] 1/2 and is not equal to tanβ, no matter how small the rotation is, the small-rotation-angle assumption would inevitably bring an error to the solution of the large-deflection problem of thin plates. Moreover, with the increase of the transversely uniformly-distributed load applied on the surface of the plate, the rotation of the thin plate will increase, and the error caused by the small-rotation-angle assumption increase accordingly. Therefore, in order to eliminate the influence of the small-rotation-angle assumption on the computational accuracy of the solution of the large deformation problem, we need to re-examine this mechanical problem by using the basic relationship between sinβ and tanβ to replace the original approximate relationship. That is to say, we will solve a large deformation problem based on large-deflection and give up the small-rotation-angle assumption.
As for the solving method, the classical perturbation technique will be followed to solve the established governing equations. Generally speaking, there are two advantages for this technique: one advantage is that we can easily obtain the relationship of load vs. central deflection if the central deflection is selected as a perturbation parameter; another may be from the asymptotic characteristics that the perturbation method itself has. In the past, there have been many successful examples of the application of perturbation in the solving of flexible thin plate. Groundbreaking works may be found in Vincent [32] and Chien [33], in which Vincent adopted the load, whereas Chien adopted the central deflection as a perturbation parameter to successfully solve the large-deflection problem of thin plates. Subsequent researchers focused mainly on the choice of perturbation parameters, for example, a generalized displacement [34], a linear function of Poisson's ratio [35], and an average angular deflection [36]. Chen and Kuang [37] discussed the differences among the possible perturbation parameters.
On the other hand, the magnitude of the perturbation parameter is an important issue that we must face during the perturbation. According to the classical works by Nayfeh [38], perturbation parameters are small parameters; however, this small parameter does not mean the central deflection is also small. Generally speaking, in light of the magnitude of the ratio of the central deflection to plate thickness, if the ratio is less than 1/5, the behavior of the plate may be regarded as rigid, in which the bending effect is dominant and the small-deflection theory is appropriate. If the ratio is greater than 5, the behavior of the plate may be regarded as an absolute flexible plate, or membrane, and the corresponding membrane theory in which the tensile effect is dominant should be adopted in this case. If the ratio falls into the range of 1/5 to 5, the behavior of the plate may be analyzed by the large-deflection theory, which was founded by Föppl-von Kármán. In our study, the ratio is limited in the range of 1/5 to 4, which obviously belongs to the large-deflection problems, although the parameter used for the perturbation attributes to a small parameter. In our next perturbation, the perturbation parameter does not have to be far less than 1, and in some cases it may be far greater than 1. This fact may also be found in many previous studies, including Chien's pioneer work [39] and, later, Shen's review [40], in which the perturbation parameter may be greater than 1 and reaches 4-5 in the post-buckling problems of plates.
Besides the choice and magnitude of perturbation parameters mentioned above, there have been many studies concerning the number and physical meaning of the parameter selected. Given that the application of this method does not belong to our study emphasis, we do not review them in detail.
In this study, the perturbation method is used to solve the large deformation problem of a bimodular FGM thin circular plate, with the emphasis on the influence of a smallrotation-angle assumption on the final results. This paper is arranged as follows. In Section 2, the governing equation of the studied problem is established, in which the equation of equilibrium is derived without small-rotation-angle assumption. In Section 3, taking the central deflection as a perturbation parameter, the perturbation solutions of stress and deflection are obtained, and the yield stress under large deformation is also analyzed. In Section 4, the influences of the small-rotation-angle assumption on the relationship of loads vs. central deflection and on the yield stress are discussed in detail. Section 5 is the concluding remarks.

Establishment of Governing Equations
A bimodular FGM thin circular plate with thickness t and radius a is subjected to a transversely uniformly-distributed load, q, as shown in Figure 1, where the dot-dash line at the peripheral of the plate represents the location of the unknown neutral layer of the plate; this location will be determined later. The polar coordinates plane (r, φ) of the cylindrical coordinates system (r, φ, z) is established on the plane where the neutral layer is located; r, φ, and z represent the radial, circumferential, and transverse coordinates, respectively. Due to the axisymmetric characteristics, φ is not presented in Figure 1; t 1 and t 2 represent the height of the tensile and compressive zones, respectively, and the corresponding modulus in the two zones is the tensile modulus, E + (z), and compressive modulus, E − (z), respectively, as shown in Figure 1. In addition, the special constraint on the plate is not provided initially, as we will solve the problem under multiple constraints. For the obtainment of an explicit solution of the problem, we need to define the mathematical form of the bimodular FGM model first. Considering the convenience of differential and integral operations, we model E + (z) and E − (z) as the following exponent type functions [30,31] E where α 1 and α 2 are the two graded indices in tensile and compressive zones, and E 0 denotes the Young's modulus of elasticity on the neutral layer. From Equation (1), it is easy to see that E + (z) = E − (z) = E 0 when α 1 = α 2 = 0 or z = 0. The bimodular property of FGM is thus simulated mathematically. At the same time, because the influence of Poisson's ratio on the deformation and stress is much less than that of Young's modulus of elasticity, Poisson ratios may be assumed as two constants, υ + and υ − , neglecting the change along the z direction. Let us take a differential element, ABCD, from the deformed circular plate to study the static equilibrium problem of this element, as shown in Figure 2, where surfaces AD and BC are located along the radial direction, while surfaces AB and DC are located along the circumferential direction. In Figure 2, N r and N θ are the radial and circumferential forces, respectively, M r and M θ are the radial and circumferential bending moments, respectively, Q r is the transverse shear force acting on surfaces AB and DC, dθ denotes the angle between the radial surfaces AD and BC, and β denotes the rotation of the radial force. Note that, due to axisymmetric characteristics, there is no increment from surface AD to surface BC, while the increment will take place from surface AB to surface DC due to the change of r. The usually so-called in-plane equilibrium equation is [41] where the body force of the plate is neglected. The equation of equilibrium along the vertical direction perpendicular to the polar coordinates plane can be written as N r + dN r dr dr (r + dr)dθ sin β + d sin β dr dr − N r rdθ sin β where the body force of the plate is also neglected. After ignoring the third order and fourth order differential items of Equation (3), we have r dN r dr drdθ sin β + N r drdθ sin β + rN r dθ d sin β dr dr +r dQ r dr drdθ + Q r drdθ − qrdrdθ = 0.
Dividing Equation (4) where sin β = 1/ 1 + 1/ tan 2 The basic relationship of trigonometric functions is adopted in Equation (6) to replace the original small-rotation-angle assumption that the rotation, β, of the deformed thin plate is so small that sinβ ≈ tanβ holds true.
At the same time, the equation of equilibrium of the moment of the differential element ABCD can be expressed as After ignoring the third order and fourth order differential items and dividing by drdθ, Equation (7) can be simplified as Substituting rQ r in Equation (8) into Equation (5) yields Suppose that the displacements of the neutral layer along the r and z directions are u and w, respectively, and the radii of the radial and circumferential curvatures of the neutral layer are ρ r and ρ θ , respectively. The strain components of any point at the z axis in the plate come from two different aspects, one is the tensile strain on the neutral layer and another from the bending strain; that is [31] where We note that, in the small-deflection theory as well as the common large-deflection theory of plate, the approximate expressions of the radial curvature and circumferential curvatures is 1/ρ r = −d 2 w/dr 2 and 1/ρ θ = −(dw/dr)/r, respectively. They are now replaced with the following exact expressions [33]: Suppose that the radial and circumferential stresses in the tensile and compressive zones is σ +/− r and σ +/− θ , respectively. The strain-stress relations give [30] Substituting Equations (10)- (12) into Equation (14), it is found that The radial and circumferential bending moments, M r and M θ , may be expressed, in terms of the radial stress and circumferential stress [31], as Substituting Equation (15) into Equation (16) yields where The integral of the items containing z in Equation (18) should be determined as zero, i.e., A + 1 + A − 1 = 0, because it is exactly the condition used for determining the position of the unknown neutral layer in [31]. Thus, the expressions of M r and M θ can be further rewritten as Substituting Equation (19) into Equation (9), we may obtain the governing equation of the bimodular FGM thin circular plates without small-rotation-angle assumption, as follows The radial force, N r , and circumferential force, N θ , are the sum of integrals in the tensile zone and the compressive zone, such that [31] After substituting Equation (15) into Equation (21), N r and N θ are rewritten as where It should be mentioned that, in the large-deflection problem of plate, any point of the plate is stretched along the radial and circumferential directions, so all integrals along the thickness direction should be performed only on the tensile components. Besides, the integrals of the items containing z are also equal to zero in Equation (21). From Equations (2) and (22), we have Substituting u in Equation (24) into the first formula of Equation (22), we have Equation (25) is the compatible equation of the large deformation problem of the bimodular FGM thin circular plate. Lastly, the improved equilibrium relation (Equation (20)), in combination with the compatible relation (Equation (25)), constitute the governing equations of the problem.

Verification of Regression and Simplification of Equations
Next, we will verify that Equation (20) can be regressed to the classical Föppl-von Kármán equation. If the mechanical properties of different modulus and FGM effect are not taken into account, the bending stiffness of the bimodular FGM thin plate D * is regressed to the bending stiffness of the isotropic plate D. The limit of D * when α 1 , Because t 1 = t 2 = t/2 when α 1 = α 2 = 0 and υ + = υ − = υ, the bending stiffness of isotropic and homogeneous plate can be obtained Besides, when the term 1 + (−dw/dr) 2 in Equation (20) is approximately equal to 1, Equation (20) may be further simplified as which is the classical Föppl-von Kármán equation in the case of large-deflection. The satisfaction of regression indicates that the improved governing equation in large deformation may regress to the classical Föppl-von Kármán equation in large-deflection. Note that there are many nonlinear items concerning 1 + (−dw/dr) 2 in Equation (20), which further makes its solving complicated; therefore, moderate simplification is necessary. For this purpose, the denominator terms in Equation (20) are expanded in the form of the power series with respect to (−dw/dr) 2 , such that and Substituting the first and second terms of Equations (29) and (30) into Equation (20), we have Integrating Equation (31) yields It can be seen that C = 0 according to the boundary conditions dw/dr = 0 and N r = 0 at r = 0. The items containing υ +/− are then all offset by calculating the first derivation in Equation (32), and the simplified equilibrium equation is now While, at the same time, the counterpart in the classical Föppl-von Kármán equation is [33] D r

Boundary Conditions
The following four edge constraints are considered, as shown in Figure 3, that is, for Edge (1), rigidly clamped, we have for Edge (2), movably clamped for Edge (3), simply hinged and for Edge (4), simply supported From Equation (24), the radial displacement, u, of the middle plane can be expressed in terms of the radial force N r , i.e., For the radial bending moment, M r , in Equation (19), we introduce 1 + (−dw/dr) 2 ≈ 1 into the expression of M r , thus expressing the M r = 0 only in terms of the deflection w, i.e.,
In addition, the axisymmetric conditions used for solving gives The dimensionless deflection at the center is chosen as the perturbation parameter [33] where w 0 is the central deflection of the circular plate, i.e., the maximum deflection. The dimensionless quantities of P, S, and W are represented by the perturbation parameter Thus, P, S, and W are expressed in the form of the power series of W m where p i (i = 1,2,3 . . . ) are undetermined constants, and ω i (η) and s i (η) (i = 1,2,3 . . . ) are undetermined functions with respect to η. The introduction of 1/16 into the expansion of P can make the numerical calculation easier.
Step 1: p 1 , ω 1 (η) and s 1 (η) The differential equation used for the solution of s 1 (η) can be obtained from the coefficient of W m in Equation (40) which should satisfy the boundary conditions as The solution of Equation (50) under the boundary conditions of Equation (51) is s 1 (η) = 0. In the same way, the differential equation used for the solution of ω 1 (η) and p 1 can be obtained from the coefficient of W m in Equation (41) The boundary conditions are Thus, p 1 = 4K/(2λ 1 + 1) and ω 1 (η) = (η 2 + 2λ 1 η)/(2λ 1 + 1) may be obtained. It is found that ω 1 (η) and p 1 is exactly the solution in the case of small-deflection.
Step 2: p 2 , ω 2 (η) and s 2 (η) The differential equation used for the solution of s 2 (η) can be obtained from the coefficient of W 2 m in Equation (40), that is which should satisfy the boundary conditions as Using ω 1 (η), obtained in the previous step, the solution of Equation (54) under the boundary conditions of Equation (55) is The differential equation used for the solutions of ω 2 (η) and p 2 from the coefficient of The boundary conditions are The results ω 2 (η) = 0 and p 2 = 0 are easily obtained.
Step 3: p 3 , ω 3 (η) and s 3 (η) The differential equation used for the solution of s 3 (η) can be obtained from the coefficient of W 3 m in Equation (40), that is which should satisfy the boundary conditions as The result s 3 (η) = 0 is easily obtained. The differential equation used for the solution of ω 3 (η) and p 3 from the coefficient of W 3 m is The boundary conditions are Using the ω 1 (η) and s 2 (η) obtained in the previous steps, the solutions of Equation (61) under the boundary conditions of Equation (62) are Step 4: p 3 , ω 3 (η) and s 3 (η) The differential equation used for the solution of s 4 (η) from the coefficient of W 4 m is which should satisfy the boundary conditions as Using the ω 1 (η) and ω 3 (η) obtained in the previous steps, the solution of Equation (65) under the boundary conditions of Equation (66) is (67) The differential equation used for the solutions of ω 4 (η) and p 4 from the coefficient of The boundary conditions are The results ω 4 (η) = 0 and p 4 = 0 are easily obtained. Thus, the unknown constants p 1 and p 3 and functions ω 1 (η), ω 3 (η), s 2 (η), and s 4 (η) have been determined. The remaining solutions of p i , ω i (η), and s i (η) (i = 5,6,7 . . . ) can be obtained through similar calculations. The further calculations will not be described in detail.

Stress Analysis
For the purpose of studying the stress of each point on the plate, it is necessary to derive the stretching stress from membrane force and the bending stress from bending moment. Let ∑ r (η) be the dimensionless stretching stress at the middle plane of the plate σ r , and let ∑ r (η) be the dimensionless bending stress on the convex surface σ , i.e., Allowing ∑ r (η) ≈ s 2 (η)W 2 m + s 4 (η)W 4 m , the dimensionless stretching stress at the edge of the plate is and the stretching stress at the center of the plate is ∑ r (1) = VW 2 m 6(2λ 1 +1) 2 The dimensionless bending stress on the convex surface is where t 1 is the height of tensile section in the plate and, for convenience in the next computation, we define its dimensionless form as T 1 = t 1 /t; at the same time, T 2 = t 2 /t. Substituting Equations (19) and (77) into Equation (74) gives Substituting the expressions of W(η) ≈ ω 1 (η)W m + ω 3 (η)W 3 m into Equation (78), the dimensionless bending stress at the edge of the plate is obtained, i.e., and the dimensionless bending stress at the center of the plate is Via the above stress expressions, we may analyze the yield problem of the bimodular FGM plate in large deformation. Let three principal stresses be σ 1 , σ 2 , and σ 3 . The yield condition gives where σ y is the yield stress. At the edge of the plate, due to u = r(N θ − ν + N r )/A 3 = 0, it is found that where σ re and σ te are the radial and circumferential stresses on the convex surface of the plate edge, respectively. Supposing σ 3 = 0, we can obtain Therefore, the dimensionless yield stress calculated with the radial stress is Substituting Equations (75) and (79) into Equation (84), the dimensionless yield stress at the edge is (85) Next, we will derive the yield stress at the center of the plate. The principal stresses should be equal at the center, i.e., where σ re is the radial stress on the convex surface of the plate center. Similarly, letting σ 3 = 0, we have σ rc = σ y .
The dimensionless yield stress calculated with the radial stress is Substituting Equations (76) and (80) into Equation (88), the dimensionless yield stress at the center is (89) In addition, under the Edges (1), (2), and (4), the dimensionless yield stresses at the center and at the peripheral edge of the plate can be easily obtained by letting λ 1 and λ 2 equal the corresponding values.

Results and Discussions
The regression of the perturbation solutions of S, P, and W obtained in Section 3.1 need to be discussed. Comparing Equations (70)-(72) in this paper with Equations (104)-(106) in [31], it is found that, if T = 0 in Equations (70)-(72), the perturbation solutions of S, P, and W will be the same as those presented in [31], which verifies, to some extent, the reliability of the perturbation solutions obtained in this paper.
and substituting the first two terms into Equation (91) yields Additionally combining T 1 + T 2 = 1, the expressions of the tensile and compressive heights are, respectively It should be noted here that, from the above equation, it seems that the tensile and compressive heights depend only on the Poisson ratio and are independent of the modulus of elasticity. In fact, both two physical quantities participate in the determination of the heights, in which the influences caused by E + (z) and E − (z) is embodied by α 1 and α 2 , which existed in the original Equation (84). However, after taking the first two items of the expansion T 1 and T 2 , the influences of the modulus of elasticity happens to disappear, which is the reason why Equations (94) and (95) contains only the Poisson ratio. If we take more items from Equation (92), a more complicated expression containing α 1 and α 2 will be obtained, and the solution in this case cannot be explicitly expressed.

Effect of Small-Rotation-Angle Assumption on Loads vs. Central Deflection
There are four cases of the positive or negative signs of graded parameters α 1 and α 2 in Equation (1): case (a) α 1 > 0, α 2 > 0; case (b) α 1 < 0, α 2 < 0; case (c) α 1 > 0, α 2 < 0; and case (d) α 1 < 0, α 2 > 0. If the positive or negative signs of α 1 and α 2 are determined, the relative magnitude among the modulus of the neutral layer, E 0 , the tensile modulus, E + (z), and compressive modulus, E − (z), will be accordingly defined, i.e., for case (a), we have Some numerical examples were carried out to discuss the effect of the small-rotationangle assumption on the response of the bimodular FGM plate. The FGM we are familiar with may be manufactured using powder metallurgy techniques from a mixture of ceramic and metal. Generally, metals exhibit good tensile properties and can be considered as materials in the tensile zone, whereas ceramics exhibit good compressive properties and can be considered as materials in the compressive zone. The Poisson ratios of the tensile and compressive zones is taken as ν + = 0.35 (manganese bronze) and ν − = 0.25 (silica ceramic), respectively. Moreover, in order to be able to describe the four cases of the relative magnitude among E + (z), E 0 , and E − (z), we chose four groups of data: case (a), α 1 = 0.5, α 2 = 0.1; case (b), α 1 = −0.5, α 2 = −0.1; case (c), α 1 = 0.5, α 2 = −0.1; and case (d), α 1 = −0.5, α 2 = 0.1. From Equations (94) and (95), the dimensionless heights of the tensile and compressive zones are T 1 = 0.4821 and T 2 = 0.5179, respectively. Thus, using Equations (18), (39), and (43) and also based on the known ν + , ν − , α 1 , α 2 , T 1 , and T 2 , the values of λ 1 , λ 2 , K, and V may be determined, and they are listed in Table 1. At the same time, by substituting λ 1 , λ 2 , K, and V into Equations (71) and (72), the expressions of P and W with respect to W m are obtained, and they are listed in Table 2.
Generally speaking, the bending stiffness of the plate depends partly on the magnitude of the modulus of elasticity of the materials; the bigger the modulus of elasticity, the stronger the bending-resistant capacity. Obviously, among the four cases of the tensile modulus, compressive modulus, and the neutral layer modulus, case (c) has the strongest bending-resistant capacity because there exists the relations E + (z) > E 0 and E − (z) > E 0 . On the contrary, case (d) has the weakest capacity resisting deformation due to E + (z) < E 0 and E − (z) < E 0 . In addition, the deformation of the thin plate is greatly affected by the boundary constraint. Next, a numerical example is given to illustrate the effect of the four boundary constraints on the deformation of thin plates. The deflection curves of the bimodular FGM thin plate under four edge constraint conditions when the central deflection of the thin plate, W m , reaches 2, are then plotted, as shown in Figure 4, where Case (a) represents E + (z) > E 0 > E − (z); Case (b) represents E + (z) < E 0 < E − (z); Case (c) represents E + (z) > E 0 , E − (z) > E 0 ; and Case (d) represents E + (z) < E 0 , E − (z) < E 0 . Note that a coordinate transform is carried out, i.e., the coordinate variable (η = 1 − r 2 /a 2 ) in the deflection expressions in Table 2 should be transformed into the ratio of the radial coordinate to the radius of the thin plate (r/a) to show the deflection curves. Table 1. Values of K, V, λ 1 , and λ 2 under Edges (1)-(4) and cases (a)-(d).

Conditions
Case (

Cases Formulas
(1) Rigidly clamped By comparing the four deflection diagrams in Figure 4, it can be seen that the two graded indices in the tensile and compressive zones, α 1 and α 2 , do not change the real shape of the deflection curves of the thin plate. In the Cases a, b, c, and d, the rotation angle of the deflection curve of Edge (2) at r/a = 0 is the largest among the four constraint conditions, while the rotation angle of the deflection curve of Edge (3) at r/a = 0 is the smallest. Moreover, for the same central deflection, when 0 < r/a < 1, the deformation degree of the thin plate under the constraint of the Edge (2) is the smallest, and the deformation of the thin plate under the constraint of the Edge (3) is the largest. In addition, when 0 < r/a < 0.4, the deflection curves of Edge (1) and Edge (4) are relatively close in all four Cases (a, b, c, d).
In the following content, we will focus our discussion on just one case, i.e., Case (a) E + (z) > E 0 > E − (z), which corresponds to the parameters α 1 = 0.5, α 2 = 0.1, and consider the bimodular FGM circular plate with the different dimensionless thicknesses, T = 0, 0.05 and 0.2, respectively, where T = 0 represents the perturbation solution from [31] with the small-rotation-angle assumption. Figure 5 shows the variations of the dimensionless external loads, P, with central deflection W m under the constraints of Edges (1)-(4). The concrete values of external loads and the errors between the conditions T = 0 and T = 0 (i.e., the errors caused by the small-rotation-angle assumption) are listed in Table 3, in which the dimensionless central deflection take different values, i.e., W m = 1, 2, 3, 4.  From Figure 5, it is easy to see that, under any of the four boundary constraints, for the curves of external loads vs. central deflection, the curve T = 0.05 is very close to the curve T = 0; however, there is an obvious distance between curve T = 0.2 and curve T = 0, especially under the Edge (1) and Edge (2) constraints. Moreover, for the same central deflection, the external load acting on the thin plate is always the smallest under the Edge (4) constraint, and it is always the largest under the Edge (1) constraint. Besides, from Table 3, it may be readily seen that, for the same magnitude value of the central deflection, the external loads calculated by the perturbation solutions presented here (for T = 0) are uniformly smaller than those calculated by the perturbation solutions presented in [31] (for T = 0). For example, under the Edge (4) constraint, when W m = 4, the external load calculated by the solution presented in [31] is 32.6609 (T = 0), while the external loads calculated by the solution presented here are 32.3610 for T = 0.05 and 27.8628 for T = 0.2. In this case, the relative error between the external load for T = 0 and that for T = 0.05 is only 0.93%, which is less than the allowable error of 1% for precision instruments. The relative error between the present solution and the existing solution is very small, which can, to some extent, be taken as a criterion to judge the reliability of the perturbation solutions presented here. Furthermore, under the Edge (4) constraint, when W m = 4 and T = 0.2, the relative error between the present solution and the existing solution has reached 17.22%, which is caused by the small-rotation-angle assumption and is even greater than the allowable error in engineering (<15%). In other words, in this case, the error of external load calculated by the perturbation solution with the small-rotation-angle assumption in [31] is 17.22%; accordingly, the computational accuracy of the external load without the small-rotation-angle assumption in this paper is improved by 17.22% in comparison with the external load for [31].
On the other hand, no matter how small rotation angle β is, because sinβ = tanβ, adopting the small-rotation-angle assumption (assuming that sinβ ≈ tanβ) will always create the obtained perturbation solution with a certain error. It can also be found from Table 3 that, as the central deflection, W m , increases, the relative error between the external loads calculated by the perturbation solution in this paper and in [31] also increases. With the increase of the central deflection, W m , the rotation angle, β, will also increase; therefore, the error brought by the small-rotation-angle assumption will also increase accordingly.
At the same time, we also find that, for different boundary constraints, the effects of the small-rotation-angle assumption on the response of plates are different. From Table 3, it can be seen that, for the same central deflection, the errors of external load under the Edges (1) and (3) constraints are less than those under the Edge (2) and (4) constraints; for example, when W m = 3 and T = 0.2, the errors of Edges (1), (3), (2), and (4) are 5.22%, 2.50%, 12.76%, and 14.68%, respectively. Because the small-rotation-angle assumption is abandoned in this paper, the established governing equation exists with five terms related to the rotation angle that are not available in the classic Föppl-von Kármán equation. Besides, when the perturbation method is adopted to solve the governing equations, the external load is expanded to the power series of the perturbation parameter (i.e., W m ). It can be considered that the computational accuracy of the perturbation solution is not only related to the perturbation parameter, but also to the rotation angle at the location of the expansion point. From the four cases in Figure 4, it can be seen that, at the center part of the plate, the rotation angles of the thin plate under the Edge (1) and (3) constraints are smaller than those under the Edge (2) and (4) constraints; thus, the computational accuracies of the perturbation solutions of the Edge (1) and (3) constraints are higher than those of Edges (2) and (4). Therefore, it is necessary to abandon the small-rotation-angle assumption in the large deformation problem of bimodular FGM thin circular plates when the central deflection or the dimensionless thickness is large, especially for Edge (2), movably clamped, and Edge (4), simply supported.

Effect of Small-Rotation-Angle Assumption on Yield Stress
Next, under a logarithmic coordinates system, we plot the curves of yield stress at the edge of the plate under three edges constraints when T = 0 and T = 0.2, as shown in Figure 6. Note that σ y = 0 at the edge of the plate under Edge (4), simply supported, constraint. The concrete values of the yield stress at the edge of the plate under Edges (1), (2), and (3) and the errors between the conditions T = 0 and T = 0 (i.e., the errors induced by the small-rotation-angle assumption) are listed in Table 4. We also plot the curves of yield stress at the center of the plate under the four edge constraints under a logarithmic coordinates system, as shown in Figure 7, and the local magnification of center yield stresses is clearly shown in Figure 8, in which the center deflection is equal to 1 to 10.   It is easy to see from Figure 6 that the yield stress at the edge of the plate increases with the increase of the central deflection W m . When 0.1 < W m < 1, the yield stress at the edge of the plate under the Edge (3) constraint is much smaller than that under the Edge (1) and (2) constraints. The difference of yield stress curves at the edge of the plate between T = 0 and T = 0.2 under the logarithmic coordinate systems is not obvious in Figure 6. From Table 4, it can be seen that the errors caused by the small-rotation-angle assumption under the Edge (1) and (3) constraints are less than 8%, but under the Edge (2) constraint, the errors caused by the small-rotation-angle assumption are quite large, and the maximum error has reached 28.8% (The allowable error for civil engineering is within 15%). For the Edge (2) constraint, substituting λ 1 = λ 2 = 0 into Equation (85), there is only one item related to the parameter T, which is an important parameter that affects the computational accuracy of the perturbation solution, but for the Edge (1) and (3) constraints, after substituting their corresponding λ 1 and λ 2 into Equation (85), the positive and negative signs of the terms related to parameter T are opposite and may be offset during the calculation. Therefore, the yield stress at the edge of the plate under the Edge (2) constraint has a huge difference between T = 0 and T = 0. This indicates that the effect of the small-rotation-angle assumption on the yield stress at the edge of the plate under the Edge (2) constraint is far greater than those under the Edge (1) and (3) constraints.
From Figures 7 and 8, it can be easily seen that, for the maximum deflection with the same magnitude, the yield stresses at the center of plate under the Edge (1), (2), and (4) constraints of T = 0.2 are obviously smaller than those of T = 0. Moreover, for the same central deflection, the yield stresses at the center of the plate under the Edge (1) and (2) constraints are greater than those under the Edge (3) and (4) constraints. Comparing Figure 7 with Figure 6, it can be seen that the effect of the small-rotation-angle assumption on the yield stress at the center is greater than that on the yield stress at the edge. The above numerical examples are sufficient to prove that the computational accuracy of the perturbation solution abandoning the small-rotation-angle assumption is improved.

Conclusions
In this study, the large deformation problem of a bimodular FGM thin circular plate subjected to a transversely uniformly-distributed load is investigated, under the condition of giving up the small-rotation-angle assumption. The effects of the small-rotation-angle assumption on the responses of the bimodular FGM thin circular plate are discussed using specific numerical examples. The following three conclusions can be drawn: (i) The improved governing equations for the large deformation problem of the bimodular functionally-graded thin circular plate without small-rotation-angle assumption can be regressed to the classical Föppl-von Kármán equations, and the corresponding perturbation solutions can also be regressed to the perturbation solutions with small-rotation-angle assumption. This verifies, to some extent, the reliability of the perturbation solutions obtained in this paper.
(ii) The perturbation solutions presented in this paper have higher computational accuracy in comparison with the existing perturbation solutions. For example, in the numerical examples carried out in Section 4, the computational accuracies of external load and yield stress are improved by 17.22% and 28.79% at most, respectively. For the maximum deflection with the same magnitude, the external load calculated by the perturbation solution presented in this paper would be less than the external load from the previous solution, and the errors of the perturbation solutions caused by the small-rotationangle assumption under the movably clamped and simply supported constraints are larger than those under the rigidly clamped and the simply hinged constraints.
(iii) The small-rotation-angle assumption has little influence on the yield stress at the edge of plate, while it has great influence on the yield stress at the center of the plate, especially for the constraint cases with edges rigidly clamped, movably clamped, and simply supported.
The problem and method presented in this study can be used for further research work. In particular, the problem concerning small-rotation-angle assumption may be extended into the establishment of mathematical models with buckling and contact phenomena for elastic plates [42], because the buckling problem of an elastic plate is usually accompanied by its large deformation. In addition, the single-parameter perturbation method based on the central deflection in this study may also be extended to the so-called multi-parameters perturbation method [43] that has been successfully used for investigating functionallygraded, thin, circular piezoelectric plates. The relative work is in progress.