Variational Solution and Numerical Simulation of Bimodular Functionally Graded Thin Circular Plates under Large Deformation

Abstract

In this study, the variational method and numerical simulation technique were used to solve the problem of bimodular functionally graded thin plates under large deformation. During the application of the variational method, the functional was established on the elastic strain energy of the plate while the variation in the functional was realized by changing undetermined coefficients in the functional. As a result, the classical Ritz method was adopted to obtain the important relationship between load and maximum deflection that is of great concern in engineering design. At the same time, the numerical simulation technique was also utilized by applying the software ABAQUS6.14.4, in which the bimodular effect and functionally graded properties of the materials were simulated by subareas in tension and compression, as well as the layering along the direction of plate thickness, respectively. This study indicates that the numerical simulation results agree with those from the variational solution, by comparing the maximum deflection of the plate, which verifies the validity of the variational solution obtained. The results presented in this study are helpful for the refined analysis and optimization design of flexible structures, which are composed of bimodular functionally graded materials, while the structure is under large deformation.

1. Introduction

In the last ten years, bimodular functionally graded materials have gradually become a new research topic in academic circles. A bimodular material [1] has different elastic moduli in tension or compression, while a single-modulus material has the same modulus in tension or compression. Functionally graded materials [2] (FGMs) are a new type of composite material, generally composed of two materials, and the composition of the two materials presents continuous gradient changes, thus avoiding interface issues effectively. On the basis of functionally graded materials, considering the bimodular characteristics of the material will undoubtedly increase the difficulty of analysis, not to mention the application of this material model to the analysis of flexible structures involving large deformation (for example, flexible thin plates [3,4]). The problem is quite challenging for the combination of nonlinearity of materials and geometrical nonlinearity, especially in terms of the analytical methods. Therefore, in this study, we try to conduct both analytical and numerical research on this problem to enrich and improve existing research works in this field. For this purpose, the review is conducted in the following order to present a complete research background. The first is the bimodular materials, functionally graded materials, and their combination in recent studies; then, we briefly review the basic analytical methods used for plates and shells; finally, the structure of this paper is presented.

Many investigations show that most materials [5,6], such as graphite, rubber, concrete, ceramics, and biomedical materials, present different strains in tension and compression when they are subjected to tensile stress and compressive stress of the same magnitude. These materials have been referred to as bimodular materials by Jones [7]. In the theoretical analysis, there are basically two material models widely adopted: Bert’s model [8] and Ambartsumyan’s model [9]. Bert’s model is mainly used in the analysis of orthotropic materials and laminated composites [10,11,12], and this model is based on the criterion of positive and negative signs in the strain of longitudinal fibers. Ambartsumyan’s model, which is established on the criterion of positive and negatives signs of principal stresses, is mainly applied to isotropic materials. This model is of particular significance in the analysis and design of structures, and our present study is based on this model. In the application of Ambartsumyan’s model, the principal stress is generally obtained as a result of the solution but not as a known quantity, which necessarily brings difficulties for describing the stress state of a point. Moreover, there is also a lack of experimental results to describe the elastic coefficient in complex stress states. Analytical solutions are only available in a few simple cases, most of them dealing with plates and beams [13,14,15]. However, in complex problems, we must turn to the finite element method (FEM) based on an iterative strategy. During each iteration, we need to judge the principal stress state of each element, thus acquiring a new elastic matrix used for the subsequent iteration. In the review of Ye et al. [16], this method is referred to as a direct iterative method of variable stiffness that has widely been used in earlier studies. Thereafter, Ma et al. [17] established a finite element iterative program to obtain buckling critical loads of bimodular rods. Given that the previous iteration of methods struggled because of the convergence difficulty of the constitutive model, Du et al. [18] established a new computational framework. Their works showed that the proposed framework can be successfully applied in solving the problem.

Functionally graded materials are a new type of composite materials, and the characteristic of its composition presents continuous gradient changes along the thickness direction, thus eliminating interface problems and presenting a smooth stress distribution. The material has been successfully used in various engineering fields since its advent, such as micro-electro-mechanical systems [19], aerospace engineering [20], civil engineering [21], and acoustics [22]. There are many works on the analysis of structural elements made of functionally graded materials, most of them dealing with beams and plates (for example, [23]). Among the studies, few consider the bimodular effect from functionally graded materials. As indicated above, most materials will show the bimodular effect (it is just a matter of whether it is obvious or not); thus, functionally graded materials seem to be no exception.

Recently, the bimodular effect of materials was further introduced into the analysis of functionally graded materials, and some works finally emerged, including bimodular FGM beams [24] and bimodular FGM plates [25,26,27,28]. Aiming at the bimodular FGM plates, a simplified theory on the neutral layer under small deflection was established in [25]; thereafter, the governing equations of the large-deflection problem of bimodular FGM thin circular plates was derived in [26]. For large-deformation problems, both the deflection and rotation angle will increase with the increase in external loads. For this purpose, a single-parameter perturbation method was used to solve the Föppl–von Kármán equations without the small-rotation-angle assumption in [27], and the biparametric perturbation method was used to solve the same problem in [28]. From the above review, it can be seen that for the analytical solution of this problem, the method is still mostly limited to the perturbation method, although extending from a single-parameter perturbation to a multiple-parameter perturbation. However, the analytical method for this problem is still relatively single. To solve this problem, various analytical methods must be sought.

In general, for plate and shell problems, there are three analytical methods widely used in theoretical analysis. The first one is the so-called series expansion method, in which the series may take all kinds of functional forms, for example, the power function, exponential function, and trigonometric function; the variational method is the second method, which is established on an energy principle (Galerkin method and Ritz method, for example); and the third is the perturbation technique. Each of the three methods has its advantages and disadvantages, which are not discussed in this article.

In large-deformation problems of plates and shells, the applications of the perturbation method and variational method both show their unique advantages. In the perturbation method, the first step is to establish the governing equation expressed in terms of the unknown displacement and stress. Then, the unknown displacement and stress are expanded in the form of ascending powers with respect to a certain small parameter (perturbation parameter). By substituting the expansions into the governing equations and boundary conditions, a set of equations determining the approximate solution of all levels are obtained. Due to the fact that the perturbation parameter either appears explicitly or is introduced artificially, in 1947, Chien [29] first selected the maximum deflection of thin plates as a perturbation parameter to acquire the perturbation solution successfully. Compared with the experimental data, Chien’s solution is accurate and regarded as a landmark. For a long period of time, Chien’s solution has been cited in subsequent studies. In the variational method, especially in the displacement variational method, that is, the Ritz method, the first step is to prescribe the displacement containing undetermined coefficients, and the prescribed displacement should satisfy the boundary conditions. As the second step, the energy functional is established, in which the total strain energy stored in the elastic body and the work done by external loads are determined in advance. By substituting the prescribed displacement into the energy functional, the so-called variation is realized only by the change in coefficients, thus determining the unknown coefficients and finally obtaining the displacement.

For large-deflection problems of thin circular plates, both Chien’s perturbation solution [29] and the corresponding variational solution [30] have given satisfactory results in the literature. Compared with the perturbation method, the variational method has a distinct advantage. Due to the fact that the displacement variational equation itself represents the equilibrium equation and stress boundary conditions, it naturally avoids the consideration for the equilibrium condition of thin plates, while in the perturbation method, the establishment of an equation of equilibrium is necessary and somewhat complicated. Recent studies [31,32] also indicated that the variational method can be successfully used in the analysis of plates and shells. First, Xue et al. [31] adopted the variational method to obtain the critical loads of cantilever vertical plates with different moduli. Thereafter, He et al. [32] also used the variational method to solve bimodular thin shells under large deformation. The studies indicated that the variational method can be used for the analytical investigation of flexible plate and shell structures, but the introduction of nonlinearity of materials will increase the complexity in the analysis. From the currently collected literature, it seems that there is still no application of the variational method to bimodular functionally graded thin plates under large deformation.

In this study, the variational method of displacement is applied to solve the large deformation problem of bimodular functionally graded thin circular plates. The purpose and scope of this work are to seek a feasible analytical method for this problem and, at the same time, this analytical method is verified by the appropriate numerical simulation technique. From point of view of the nonlinearity of problems, the analytical solution for bimodular functionally graded thin plates under large deformation is challenging because the nonlinearities of materials and geometry that are intertwined further makes the obtainment of analytical solution more complicated. To this end, the whole article is organized as follows. In Section 2, the variational method and the bimodular FGM thin circular plate problem are briefly described. In Section 3, the physical equations of bimodular functionally graded materials and the geometrical relations under large deformation are presented. In Section 4, the total strain energy of the plate is derived first and then the Ritz method is adopted to solve the large deformation problem of bimodular functionally graded thin circular plates. Section 5 shows the numerical simulation and the comparisons with the variational solution and the previous perturbation solution. Section 6 shows the corresponding results and discussion, and the concluding remarks are given in Section 7.

2. Method and Problem

2.1. Displacement Variational Method

In a spatial axisymmetric problem of elasticity, the cylindrical coordinate system is established as O-rθz, in which O denotes the origin; r and θ denote the radial and circumferential direction, respectively; and z denotes the direction normal to the rOθ plane, as shown in Figure 1. Let σ_r, σ_θ, and σ_z be the normal stress along the radial, circumferential, and z directions, respectively; and τ_rθ = τ_θr, τ_rz = τ_zr, and τ_zθ = τ_θz be the three shearing stress components. Due to the axisymmetric characteristic, τ_rθ = τ_θr = 0 and τ_zθ = τ_θz = 0, and there are four stress components in total remaining, σ_r, σ_θ, σ_z, and τ_rz = τ_zr, which are the functions of r and z (please refer to Figure 1). Let ε_r, ε_θ, and ε_z be the normal strain along the radial, circumferential, and z directions, respectively; and let γ_zr be the shearing strain of r and z directions. In addition, let u_r, u_θ, and w be the radial, circumferential, and z direction displacements, respectively. Note that, due to the axisymmetry, u_θ = 0.

The geometrical equation of the spatial axisymmetric problem gives [33]

$${\epsilon }_r=\frac{\partial u_r}{\partial r}, {\epsilon }_{\theta }=\frac{u_r}{r}, {\epsilon }_z=\frac{\partial w}{\partial z}, {\gamma }_{zr}=\frac{\partial u_r}{\partial z}+\frac{\partial w}{\partial r}.$$

(1)

At the same time, the physical equation of the spatial axisymmetric problem gives [33]

$$\{\begin{array}{l}{\epsilon }_r=\frac{1}{E}[{\sigma }_r-\mu ({\sigma }_{\theta }+{\sigma }_z)]\\ {\epsilon }_{\theta }=\frac{1}{E}[{\sigma }_{\theta }-\mu ({\sigma }_z+{\sigma }_r)]\\ {\epsilon }_z=\frac{1}{E}[{\sigma }_z-\mu ({\sigma }_r+{\sigma }_{\theta })]\\ {\gamma }_{zr}=\frac{2(1+\mu )}{E}{\tau }_{zr}\end{array},$$

(2a)

where μ and E are the Poisson’s ratio and modulus of elasticity, respectively. Alternatively, we may give another form of the physical equation as

$$\{\begin{array}{l}{\sigma }_r=\frac{E}{1+\mu }\left[\frac{\mu }{1-2\mu }({\epsilon }_r+{\epsilon }_{\theta }+{\epsilon }_z)+{\epsilon }_r\right]\\ {\sigma }_{\theta }=\frac{E}{1+\mu }\left[\frac{\mu }{1-2\mu }({\epsilon }_r+{\epsilon }_{\theta }+{\epsilon }_z)+{\epsilon }_{\theta }\right]\\ {\sigma }_z=\frac{E}{1+\mu }\left[\frac{\mu }{1-2\mu }({\epsilon }_r+{\epsilon }_{\theta }+{\epsilon }_z)+{\epsilon }_z\right]\\ {\tau }_{zr}=\frac{E}{2(1+\mu )}{\gamma }_{zr}\end{array}.$$

(2b)

The total strain energy stored in the whole elastic body, U, is expressed in stress and strain as

$$U=\frac{1}{2}{\displaystyle \int {\displaystyle \int {\displaystyle \int ({\sigma }_r{\epsilon }_r+{\sigma }_{\theta }{\epsilon }_{\theta }+{\sigma }_z{\epsilon }_z+{\tau }_{zr}{\gamma }_{zr})dV}}}.$$

(3)

Obviously, via the geometrical equation and physical equation, the strain potential energy may be expressed in terms of the displacement components, u_r and w, which opens possibilities for the application of the displacement variational method.

Under a cylindrical coordinate system, we suppose that a spatial axisymmetric elastic body is subjected to external forces including the body force and surface force along the r, θ, and z directions; that is, F_r, F_θ, and F_z, as well as $\overline{F_r}$, $\overline{F_{\theta }}$, and $\overline{F_z}$, and the elastic body are now in equilibrium. The resulting displacements, u_r, u_θ, and w, should satisfy the equation of equilibrium, displacement boundary conditions, as well as stress boundary conditions. If the displacements cause minor changes allowed by the boundary conditions, the new displacements will become (note that, due to the axisymmetry, u_θ = 0)

$$u_r^{\ast }=u_r+\delta u_r, w^{\ast }=w+\delta w,$$

(4)

where δu_r and δw are the virtual displacements that occur. During the virtual displacement, if there are no changes in thermal and kinetic energies, according to the principle of energy conservation, the increment in strain potential energy, δU, is equal to the work done by the external forces; therefore, the displacement variational equation may be obtained as follows [33]

$$\delta U={\displaystyle \iiint (F_r\delta u_r+F_z\delta w)}dV+{\displaystyle \iint (\overline{F_r}\delta u_r+\overline{F_z}\delta w})dS,$$

(5)

which is referred to as the Lagrangian variational equation. This variational equation provides an approximate solution to elastic problems. More specifically, if a group of displacements containing a series of unknown coefficients satisfy the displacement boundary conditions, Equation (5) may be used for the determination of these coefficients, thus obtaining the displacement.

The displacement expression is taken as

$$u_r=u_r{}_0+{\displaystyle \sum _m{A}_m{u}_{rm}}, w=w_0+{\displaystyle \sum _m{C}_m{w}_m},$$

(6)

where A_m and C_m are the independent coefficients; u_r0 and w₀ are the specified functions whose boundary value is equal to the known quantity at the boundary; and u_rm and w_m are given functions that are equal to zero at the boundary. Therefore, regardless of how A_m and C_m are taken, the displacements u_r and w always satisfy the boundary conditions. Note that because the displacement variation is obtained only by changing A_m and C_m, according to Equation (6), the variation of displacement is

$$\delta u_r={\displaystyle \sum _m{u}_r{}_m\delta A_m}, \delta w={\displaystyle \sum _m{w}_m\delta C_m}.$$

(7)

The variation of strain energy gives

$$\delta U={\displaystyle \sum _m\left(\frac{\partial U}{\partial A_m}\delta A_m+\frac{\partial U}{\partial C_m}\delta C_m\right)}.$$

(8)

Substituting Equations (7) and (8) into Equation (5) will yield

$$\begin{array}{l}{\displaystyle \sum _m\left(\frac{\partial U}{\partial A_m}-{\displaystyle \iiint F_r{u}_{rm}}dV-{\displaystyle \iint \overline{F_r}{u}_{rm}}dS\right)\delta A_m}\\ +{\displaystyle \sum _m\left(\frac{\partial U}{\partial C_m}-{\displaystyle \iiint F_z{w}_m}dV-{\displaystyle \iint \overline{F_z}{w}_m}dS\right)\delta C_m=0}\end{array}.$$

(9)

Because the variations δA_m and δC_m are arbitrary and independent from one another, the coefficients of these variations in Equation (9) must be zero, thus obtaining the following two relations:

$$\{\begin{array}{l}\frac{\partial U}{\partial A_m}={\displaystyle \iiint F_r{u}_{rm}}dV+{\displaystyle \iint \overline{F_r}{u}_{rm}}dS\\ \frac{\partial U}{\partial C_m}={\displaystyle \iiint F_z{w}_m}dV+{\displaystyle \iint \overline{F_z}{w}_m}dS\end{array},$$

(10)

which are used for solving the undermined coefficients; thus, the displacement may be obtained via Equation (6). In many references [30,33,34], the variational method based on displacement is also called the Ritz method.

2.2. Description of Problem

As shown in Figure 2, a bimodular FGM thin circular plate is subjected to a transversely uniformly distributed load q, in which t is the plate thickness and a denotes the radius of the plate. The origin O of cylindrical coordinates system (O-rθz) is set at the plate center on the neutral layer; r, θ, and z denote the radial, circumferential, and transverse coordinates, respectively. For the reason of axisymmetry, θ is not depicted in Figure 2.

Note that in Figure 2, the dot-dashed line stands for the location of the unknown neutral layer of the plate, which is determined next. In general, due to the introduction of bimodular functionally graded materials, the neutral layer does not coincide with the geometrical middle plane of the plate. In Figure 2, t₁ and t₂ are the tensile thickness and compressive thickness, respectively, and the corresponding modulus of the material on the two thicknesses is the tensile modulus E⁺(z) and compressive modulus E⁻(z), which is a function of z since the functionally graded property is varied along the thickness direction. To facilitate the application of the displacement variational method, the prescribed displacement should satisfy all displacement boundary conditions. Thus, the constraints for the thin circular plate are considered as fully fixed on its peripheral.

For the convenience of differential and integral operations, E⁺(z) and E⁻(z) are defined as exponent-type functions [25], such that

$$E^{+}(z)=E_0{e}^{{\alpha }_{1}z/t}, E^{-}(z)=E_0{e}^{{\alpha }_{2}z/t},$$

(11)

where α₁ and α₂ are two graded indices of the tensile zone and compressive one, respectively, and E₀ is the elastic modulus on the neutral layer. From Equation (11), it is found that when α₁ = α₂ = 0 or z = 0, E⁺(z) = E⁻(z) = E₀. Meanwhile, according to common practice, the Poisson’s ratio is assumed as two constants, μ⁺ and μ⁻, ignoring the gradient change along the direction of z.

In addition, the determination of the unknown neutral layer (t₁ and t₂) may be via two different conditions, according to our previous study [25]. One is the equilibrium condition, that is, the radial and circumferential normal forces acting on the differential element are zero; the other is the continuity condition, that is, the stresses acting on the neutral layer are continuous. In [25], the equilibrium condition gives

$$\frac{A_1^{+}}{1-{\mu }^{+}}+\frac{A_1^{-}}{1-{\mu }^{-}}=0,$$

(12a)

while also in [25], the continuity condition gives

$$\frac{A_1^{+}}{1-{({\mu }^{+})}^2}+\frac{A_1^{-}}{1-{({\mu }^{-})}^2}=0,$$

(12b)

in which

$$\{\begin{array}{l}{A}_1^{+}={\displaystyle {\int }_0^{t_1}z{e}^{{\alpha }_{1}z/t}dz}=\left[\left(\frac{t{t}_1^{}}{{\alpha }_1}-\frac{t^2}{{\alpha }_1^2}\right)e^{{\alpha }_1{t}_1/t}+\frac{t^2}{{\alpha }_1^2}\right]\\ A_1^{-}={\displaystyle {\int }_{-t_2}^{0}z{e}^{{\alpha }_{2}z/t}dz}=\left[\left(\frac{t{t}_2^{}}{{\alpha }_2}+\frac{t^2}{{\alpha }_2^2}\right)e^{-{\alpha }_2{t}_2/t}-\frac{t^2}{{\alpha }_2^2}\right]\end{array}.$$

(13)

The obtainment process in detail for Equation (12a,b) may refer to our previous study [25], so there is no need to repeat the derivation process again. During the obtainment of Equation (12a,b), we first need to give the functional forms of E⁺(z) and E⁻(z) for the subsequent integral operation, which is also the reason why Equation (11) is given first.

In addition, we note that for the two solutions concerning the neutral layer, it is obvious that the difference between them is slight due to the fact that the values for the Poisson’s ratio generally fall into the range of 0 to 0.3; also in Equation (12a,b), there is always a larger number compared to it (here, it is unit 1) before the Poisson’s ratio, thus making the influence of Poisson’s ratio on the solution small. Moreover, according to our previous study [35], if the influence of Poisson’s ratio is completely neglected, both Equation (12a,b) are reduced to $A_1^{+}$ + $A_1^{-}$ = 0, which is exactly the solution used for the determination of the unknown neutral axis of bimodular FGM beams.

3. Geometrical and Physical Equations of Thin Circular Plates

Note that in a spatial axisymmetric problem, there exist four stress components in total, σ_r, σ_θ, σ_z, and τ_rz, and their corresponding strain components, ε_r, ε_θ, ε_z, and γ_rz; thus, the geometrical and physical relations will involve these stresses and strains. In the bending problem of thin plates, according to the classical Kirchhoff hypothesis, ε_z is negligibly small and γ_rz may be regarded as zero; thus, the geometrical and physical equations will finally involve the two main stresses, σ_r and σ_θ, as well as the corresponding strains, ε_r and ε_θ. Note that σ_z and τ_rz are not zero; they will participate in the equilibrium conditions; however, only the geometric and physical relations are discussed here.

3.1. Geometrical Equations under Large Deformation

The geometrical relation under small-deflection bending may be expressed in terms of the curvature [33] as follows, according to the classical Kirchhoff hypothesis:

$${\epsilon }_{rb}=\frac{z}{{\rho }_r}, {\epsilon }_{\theta b}=\frac{z}{{\rho }_{\theta }},$$

(14)

where ε_rb and ε_θb are the radial and circumferential strain under small-deflection bending, respectively; ρ_r and ρ_θ are the curvature radius along the radial and circumferential directions, respectively; and in the case of small rotation angle, they are the following familiar forms [33]:

$$\frac{1}{{\rho }_r}=-\frac{d^{2}w}{d{r}^2}, \frac{1}{{\rho }_{\theta }}=-\frac{1}{r}\frac{dw}{dr},$$

(15)

where w is the transverse displacement or the deflection. Thus, the geometrical relation under small-deflection bending may be expressed in terms of w as follows:

$${\epsilon }_{rb}=-z\frac{d^{2}w}{d{r}^2}, {\epsilon }_{\theta b}=-z\frac{1}{r}\frac{dw}{dr}.$$

(16)

At the same time, the geometrical relation between in-plane displacements and in-plane strain will give [33]

$${\epsilon }_{rm}=\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2, {\epsilon }_{\theta m}=\frac{u_r}{r},$$

(17)

where ε_rm and ε_θm are the in-plane strains along the radial and circumferential directions, respectively; u_r and w are the radial displacement and deflection, as indicated above. Finally, the total geometrical relation will give

$$\{\begin{array}{l}{\epsilon }_r={\epsilon }_{rm}+{\epsilon }_{rb}=\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2-z\frac{d^{2}w}{d{r}^2}\\ {\epsilon }_{\theta }={\epsilon }_{\theta m}+{\epsilon }_{\theta b}=\frac{u_r}{r}-z\frac{1}{r}\frac{dw}{dr}\end{array}.$$

(18)

If the plate is under small deflection, the above relation will change into the form of Equation (16). In the case of small deflection, only the bending effect is considered while the membrane effect is neglected.

3.2. Physical Equations

We suppose the radial and circumferential bending stresses in tensile and compressive zones are ${\sigma }_{rb}^{+/-}$ and ${\sigma }_{\theta b}^{+/-}$, respectively, in which the subscript b stands for the bending, and the stress–strain relations, in the tensile zone, will give

$$\{\begin{array}{l}{\sigma }_{rb}^{+}=\frac{E^{+}(z)}{1-{({\mu }^{+})}^2}({\epsilon }_{rb}+{\mu }^{+}{\epsilon }_{\theta b})\\ {\sigma }_{\theta b}^{+}=\frac{E^{+}(z)}{1-{({\mu }^{+})}^2}({\epsilon }_{\theta b}+{\mu }^{+}{\epsilon }_{rb})\end{array} , at 0\le z\le t_1,$$

(19a)

and in the compressive zone,

$$\{\begin{array}{l}{\sigma }_{rb}^{-}=\frac{E^{-}(z)}{1-{({\mu }^{-})}^2}({\epsilon }_{rb}+{\mu }^{-}{\epsilon }_{\theta b})\\ {\sigma }_{\theta b}^{-}=\frac{E^{-}(z)}{1-{({\mu }^{-})}^2}({\epsilon }_{\theta b}+{\mu }^{-}{\epsilon }_{rb})\end{array} , at -t_2\le z\le 0.$$

(19b)

At the same time, we note that under large deformation, the in-plane stresses acting on the whole thickness of the cross-section are always tensile; thus. the membrane stress may be changed to ${\sigma }_{rm}^{+}$ and ${\sigma }_{\theta m}^{+}$, in which the subscript m stands for the membrane stress, and the elastic modulus and Poisson’s ratio may also be changed to E⁺(z) and μ⁺, respectively. The physical equation of in-plane deformation may be expressed, along the whole thickness direction, as

$$\{\begin{array}{l}{\sigma }_{rm}^{+}=\frac{E^{+}(z)}{1-{({\mu }^{+})}^2}({\epsilon }_{rm}+{\mu }^{+}{\epsilon }_{\theta m})\\ {\sigma }_{\theta m}^{+}=\frac{E^{+}(z)}{1-{({\mu }^{+})}^2}({\epsilon }_{\theta m}+{\mu }^{+}{\epsilon }_{rm})\end{array} , at -t_2\le z\le t_1.$$

(20)

Next, the geometrical and physical equations obtained above are used for the establishment of the functional of energy.

4. Displacement Variational Method

4.1. Total Strain Energy

The total strain energy, U, consists of the energy produced by the bending deformation, U_b, and the energy produced by the deformation of middle surface, U_m, that is [30],

$$U=U_b+U_m,$$

(21)

where the subscript b denotes the bending deformation and the subscript m denote the in-plane membrane deformation, which is consistent with the above notational conventions in geometrical and physical equations.

First, the strain energy produced by the middle surface deformation, U_m, is computed as follows [30]:

$$U_m=\frac{1}{2}{\displaystyle \underset{V}{\iiint }\left({\sigma }_{rm}^{+}{\epsilon }_{rm}+{\sigma }_{\theta m}^{+}{\epsilon }_{\theta m}\right)dV}.$$

(22)

Substituting Equation (20) into Equation (22) and also noticing that the lower limit and upper limit of the integral along the z direction are −t₂ and t₁, respectively, and dV = dzdS, Equation (22) may be written as

$$U_m=\frac{1}{2}\frac{1}{1-{({\mu }^{+})}^2}{\displaystyle {\int }_{-t_2}^{t_1}{E}^{+}(z)dz}{\displaystyle \underset{S}{\iint }\left[{({\epsilon }_{rm})}^2+{({\epsilon }_{\theta m})}^2+2{\mu }^{+}{\epsilon }_{rm}{\epsilon }_{\theta m}\right]dS}.$$

(23)

Also, substituting Equation (17) into Equation (23) will yield

$$U_m=\frac{1}{2}\frac{1}{1-{({\mu }^{+})}^2}{\displaystyle {\int }_{-t_2}^{t_1}{E}^{+}(z)dz}{\displaystyle \underset{S}{\iint }\left\{{\left[\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2\right]}^2+{\left(\frac{u_r}{r}\right)}^2+2{\mu }^{+}\frac{u_r}{r}\left[\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2\right]\right\}dS}.$$

(24)

Thus, U_m is expressed in terms of u_r and w. Considering Equation (11), we have

$$A_0={\displaystyle {\int }_{-t_2}^{t_1}{E}^{+}(z)}dz={\displaystyle {\int }_{-t_2}^{t_1}{E}_0{e}^{{\alpha }_{1}z/t}dz}=\frac{E_{0}t}{{\alpha }_1}\left(\frac{e^{{\alpha }_1}-1}{e^{{\alpha }_1{t}_2/t}}\right),$$

(25)

and also noting dS = rdθdr, U_m may be further computed as

$$U_m=\frac{\pi A_0}{1-{({\mu }^{+})}^2}{\displaystyle \int \left\{r{\left[\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2\right]}^2+\frac{u_r{}^2}{r}+2{\mu }^{+}{u}_r\left[\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2\right]\right\}dr}.$$

(26)

The energy produced by the bending deformation, U_b, can be derived by the above subareas in tension and compression, that is,

$$U_b=\frac{1}{2}{\displaystyle \underset{V^{+}}{\iiint }\left({\sigma }_{rb}^{+}{\epsilon }_{rb}+{\sigma }_{\theta b}^{+}{\epsilon }_{\theta b}\right)dV}+\frac{1}{2}{\displaystyle \underset{V^{-}}{\iiint }\left(\left({\sigma }_{rb}^{-}{\epsilon }_{rb}+{\sigma }_{\theta b}^{-}{\epsilon }_{\theta b}\right)\right)dV}.$$

(27)

Substituting Equation (19a,b) into Equation (27), we have

$$\begin{array}{r}{U}_b=\frac{1}{2}\frac{1}{1-{({\mu }^{+})}^2}{\displaystyle \underset{V^{+}}{\iiint }{E}^{+}(z)\left[{\left({\epsilon }_{rb}\right)}^2+{\left({\epsilon }_{\theta b}\right)}^2+2{\mu }^{+}{\epsilon }_{rb}{\epsilon }_{\theta b}\right]dV}\\ +\frac{1}{2}\frac{1}{1-{({\mu }^{-})}^2}{\displaystyle \underset{V^{-}}{\iiint }{E}^{-}(z)\left[{\left({\epsilon }_{rb}\right)}^2+{\left({\epsilon }_{\theta b}\right)}^2+2{\mu }^{-}{\epsilon }_{rb}{\epsilon }_{\theta b}\right]dV}\end{array}.$$

(28)

Substituting Equation (16) into Equation (28), and also noticing that the range of integrals in the tensile term is from 0 to t₁ while the range in the compressive term is from −t₂ to 0, we have

$$\begin{array}{r}{U}_b=\frac{1}{2}\frac{1}{1-{({\mu }^{+})}^2}{\displaystyle {\int }_0^{t_1}{z}^2{E}^{+}(z)dz}{\displaystyle \iint \left[{\left(\frac{d^{2}w}{d{r}^2}\right)}^2+{\left(\frac{1}{r}\frac{dw}{dr}\right)}^2+2{\mu }^{+}\frac{1}{r}\frac{dw}{dr}\frac{d^{2}w}{d{r}^2}\right]}dS\\ +\frac{1}{2}\frac{1}{1-{({\mu }^{-})}^2}{\displaystyle {\int }_{-t_2}^0{z}^2{E}^{-}(z)dz}{\displaystyle \iint \left[{\left(\frac{d^{2}w}{d{r}^2}\right)}^2+{\left(\frac{1}{r}\frac{dw}{dr}\right)}^2+2{\mu }^{-}\frac{1}{r}\frac{dw}{dr}\frac{d^{2}w}{d{r}^2}\right]}dS\end{array}.$$

(29)

If we let

$$\{\begin{array}{ll}{A}_2^{+}& =\frac{1}{1-{({\mu }^{+})}^2}{\displaystyle {\int }_0^{t_1}{z}^2{E}^{+}(z)dz}=\frac{E_0}{1-{({\mu }^{+})}^2}{\displaystyle {\int }_0^{t_1}{z}^2{e}^{{\alpha }_{1}z/t}dz}\\ & =\frac{E_0}{1-{({\mu }^{+})}^2}\left[\left(\frac{2t^3}{{\alpha }_1^3}+\frac{t_1^{2}t}{{\alpha }_1}-\frac{2t^2{t}_1}{{\alpha }_1^2}\right)e^{{\alpha }_1{t}_1/t}-\frac{2t^3}{{\alpha }_1^3}\right]\\ A_2^{-}& =\frac{1}{1-{({\mu }^{-})}^2}{\displaystyle {\int }_{-t_2}^0{z}^2{E}^{-}(z)dz}=\frac{E_0}{1-{({\mu }^{-})}^2}{\displaystyle {\int }_{-t_2}^0{z}^2{e}^{{\alpha }_{2}z/t}dz}\\ & =\frac{E_0}{1-{({\mu }^{-})}^2}\left[-\left(\frac{2t^3}{{\alpha }_2^3}+\frac{t_2^{2}t}{{\alpha }_2}+\frac{2t^2{t}_2}{{\alpha }_2^2}\right)e^{-{\alpha }_2{t}_2/t}+\frac{2t^3}{{\alpha }_2^3}\right]\end{array},$$

(30)

and also noting dS = rdθdr, U_m may be finally computed as

$$\begin{array}{r}{U}_b=\pi A_2^{+}{\displaystyle \int \left[r{\left(\frac{d^{2}w}{d{r}^2}\right)}^2+\frac{1}{r}{\left(\frac{dw}{dr}\right)}^2+2{\mu }^{+}\frac{dw}{dr}\frac{d^{2}w}{d{r}^2}\right]dr}\\ +\pi A_2^{-}{\displaystyle \int \left[r{\left(\frac{d^{2}w}{d{r}^2}\right)}^2+\frac{1}{r}{\left(\frac{dw}{dr}\right)}^2+2{\mu }^{-}\frac{dw}{dr}\frac{d^{2}w}{d{r}^2}\right]dr}\end{array}.$$

(31)

Further, Equation (31) may be simplified if the peripheral of the circular plate is fully fixed. Note that the last term of the integrand in Equation (31) may be written as

$${\displaystyle {\int }_0^a\frac{d^{2}w}{d{r}^2}}\frac{dw}{dr}dr={\displaystyle {\int }_0^a\frac{dw}{dr}d\left(\frac{dw}{dr}\right)=\frac{1}{2}}{{\left(\frac{dw}{dr}\right)}^2|}_0^a,$$

(32)

in which a is the radius of the circular plate, as shown in Figure 2. If the peripheral of the circular plate is fully fixed, we may have dw/dr = 0 at r = a; at the same time, the axisymmetric condition also gives dw/dr = 0 at r = 0; it is obvious that, lastly, we have

$${\displaystyle {\int }_0^a\frac{d^{2}w}{d{r}^2}}\frac{dw}{dr}dr=0.$$

(33)

Thus, Equation (32) is simplified as

$$U_b=\pi \left(A_2^{+}+A_2^{-}\right){\displaystyle \int \left[r{\left(\frac{d^{2}w}{d{r}^2}\right)}^2+\frac{1}{r}{\left(\frac{dw}{dr}\right)}^2\right]dr}=\pi D^{\ast }{\displaystyle \int \left[r{\left(\frac{d^{2}w}{d{r}^2}\right)}^2+\frac{1}{r}{\left(\frac{dw}{dr}\right)}^2\right]dr},$$

(34)

in which D^* is exactly the bending stiffness of the bimodular FGM plate,

$$D^{\ast }=A_2^{+}+A_2^{-},$$

(35)

which indicates that the bending stiffness is still obtained via the derivation of bending strain energy, not via the equilibrium relation in our previous study [25].

Finally, we obtain the total strain potential energy, U,

$$\begin{array}{ll}U=U_b+U_m& =\pi D^{\ast }{\displaystyle \int \left[r{\left(\frac{d^{2}w}{d{r}^2}\right)}^2+\frac{1}{r}{\left(\frac{dw}{dr}\right)}^2\right]dr}\\ & +\frac{\pi A_0}{1-{({\mu }^{+})}^2}{\displaystyle \int \left\{r{\left[\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2\right]}^2+\frac{u_r{}^2}{r}+2{\mu }^{+}{u}_r\left[\frac{d{u}_r}{dr}+\frac{1}{2}{\left(\frac{dw}{dr}\right)}^2\right]\right\}dr}\end{array},$$

(36)

which is expressed in terms of the displacement components, u_r and w. Note that for the case of small deflection, in Equation (36), only the bending term U_b is retained, while the membrane force term U_m is omitted.

4.2. Ritz Method

For the large-deformation problem of thin circular plates, we take the following displacement components (note that due to the axisymmetry, u_θ = 0):

$$u_r={\displaystyle \sum _m{A}_m{u}_{rm}}, w={\displaystyle \sum _m{C}_m{w}_m},$$

(37)

where A_m and C_m are the independent coefficients, and u_rm and w_m are the specified functions that are equal to zero on the boundaries. Thus, the displacements, u_r and w, always satisfy displacement boundary conditions. According to the variational method in Section 2.1 and also considering Equation (10), the following two variational equations may be obtained:

$$\frac{\partial U}{\partial A_m}={\displaystyle \iiint F_r{u}_{rm}}dV+{\displaystyle \iint \overline{F_r}{u}_{rm}}dS=0$$

(38)

and

$$\frac{\partial U}{\partial C_m}={\displaystyle \iiint F_z{w}_m}dV+{\displaystyle \iint \overline{F_z}{w}_m}dS={\displaystyle \iint q{w}_m}dS=2\pi {\displaystyle \int q{w}_{m}rdr},$$

(39)

where $\overline{F_z}=q$. Equations (38) and (39) are used for solving A_m and C_m. We adopt the following forms of the radial displacement and deflection:

$$u_r=\left(1-\frac{r}{a}\right)\frac{r}{a}\left[A_0+A_1\frac{r}{a}+A_2{\left(\frac{r}{a}\right)}^2+\cdot \cdot \cdot \right]$$

(40)

and

$$w={\left(1-\frac{r^2}{a^2}\right)}^2\left[C_0+C_1\left(1-\frac{r^2}{a^2}\right)+C_2{\left(1-\frac{r^2}{a^2}\right)}^2+\cdot \cdot \cdot \right].$$

(41)

Obviously, regardless of how the coefficients are chosen, the displacements satisfy the boundary conditions of displacement at the peripheral:

$$u_r=0, w=0, \frac{dw}{dr}=0 at r=a,$$

(42)

and the axisymmetric conditions at the center:

$$u_r=0, \frac{dw}{dr}=0 at r=0.$$

(43)

According to the conclusion from [30,33,34], by taking the first few terms, the variational method can give satisfactory results. Thus, in the next computation, for convenience, we take A₀ and A₁ in Equation (40) and C₀ in Equation (41) and then substitute these two displacement formulas into Equation (34):

$$U_b=\frac{32\pi D^{\ast }}{3a^2}{C}_0^2.$$

(44)

And substituting these displacement formulas into Equation (26), we have

$$\begin{array}{ll}{U}_m& =\frac{E_0\pi t}{1260[1-{({\mu }^{+})}^2]{\alpha }_1{a}^2{e}^{{\alpha }_1{t}_2/t}}\left(e^{{\alpha }_1}-1\right)\\ & \times \left(\begin{array}{l}328{\mu }^{+}a{A}_0{C}_0^2+176{\mu }^{+}a{A}_1{C}_0^2+315a^2{A}_0^2+378a^2{A}_0{A}_1\\ -184a{A}_0{C}_0^2+147a^2{A}_1^2+16a{A}_1{C}_0^2+384C_0^4\end{array}\right)\end{array}.$$

(45)

In addition,

$$2\pi {\displaystyle \int q{w}_{m}rdr}=2\pi q{\displaystyle {\int }_0^a{\left(1-\frac{r^2}{a^2}\right)}^{2}rdr=\frac{\pi }{3}}q{a}^2.$$

(46)

According to Equations (21), (38), and (39), we have

$$\frac{\partial }{\partial A_0}(U_b+U_m)=0,\frac{\partial }{\partial A_1}(U_b+U_m)=0,$$

(47)

and

$$\frac{\partial }{\partial C_0}(U_b+U_m)=\frac{\pi }{3}q{a}^2.$$

(48)

Substituting Equations (44) and (45) into Equation (47), we have

$$\{\begin{array}{l}\frac{\partial U}{\partial A_0}=\frac{E_0\pi t}{1260[1-{({\mu }^{+})}^2]{\alpha }_1{a}^2{e}^{{\alpha }_1{t}_2/t}}\left(e^{{\alpha }_1}-1\right)\left(\begin{array}{l}328{\mu }^{+}a{C}_0{}^2+630a^2{A}_0\\ +378a^2{A}_1-184a{C}_0{}^2\end{array}\right)=0\\ \frac{\partial U}{\partial A_1}=\frac{E_0\pi t}{1260[1-{({\mu }^{+})}^2]{\alpha }_1{a}^2{e}^{{\alpha }_1{t}_2/t}}\left(e^{{\alpha }_1}-1\right)\left(\begin{array}{l}176{\mu }^{+}a{C}_0{}^2+378a^2{A}_0\\ +294a^2{A}_1+16a{C}_0{}^2\end{array}\right)=0\end{array},$$

(49)

and then we express A₀ and A₁ with C₀:

$$A_0=-\frac{C_0^2}{126a}\left(89{\mu }^{+}-179\right), A_1=\frac{C_0^2}{42a}\left(13{\mu }^{+}-79\right).$$

(50)

At the same time, substituting Equations (44)–(46) into Equation (48), we have

$$\begin{array}{l}\frac{E_0\pi t}{1260[1-{({\mu }^{+})}^2]{\alpha }_1{a}^2{e}^{{\alpha }_1{t}_2/t}}\left(e^{{\alpha }_1}-1\right)\left(\begin{array}{l}656{\mu }^{+}a{A}_0{C}_0+352{\mu }^{+}a{A}_1{C}_0\\ -368a{A}_0{C}_0+32a{A}_1{C}_0+1536C_0^3\end{array}\right)\\ +\frac{64\pi D^{\ast }}{3a^2}{C}_0-\frac{\pi }{3}q{a}^2=0\end{array}.$$

(51)

Substituting Equation (50) into Equation (51), we finally have

$$H{C}_0^3+\frac{64D^{\ast }}{3}{C}_0=\frac{1}{3}q{a}^4,$$

(52)

where

$$H=-\frac{2E_{0}t(e^{{\alpha }_1}-1)}{19845[1-{({\mu }^{+})}^2]{\alpha }_1{e}^{{\alpha }_1{t}_2/t}}\left[2791{({\mu }^{+})}^2-4250{\mu }^{+}-7505\right].$$

(53)

Thus, C₀ may be solved and, according to Equation (50), A₀ and A₁ may also be obtained accordingly. Once the displacements become known, the corresponding stresses and strains may be obtained. Note that C₀ also stands for the central deflection of the thin circular plate; therefore, Equation (52) presents the important relationship of load vs. central deflection. In order to better clarify the solution process of the variational method, we add a flow block diagram for reference (please see Figure 3).

In addition, we note that the above solution is derived on the simple form of the displacement functions u_r and w; that is, A₀ and A₁ are taken in Equation (40), and only C₀ is taken in Equation (41). Although the next comparison with the numerical simulation will show its reliability in the case of small-number terms, for higher precision, more terms in displacement functions u_r and w are necessary; thus, the computation for more terms was also conducted. But, for the sake of coherence of this study, we put this part into Appendix A for interested readers, in which A₀, A₁, and A₂ are taken in Equation (40) and C₀ and C₁ are taken in Equation (41).

5. Numerical Simulation and Comparison with Variational Solution

We use the software ABAQUS6.14.4 to conduct the numerical simulation. When constructing the computational model of a thin circular plate, the first step is to create the three-dimensional solid diagram based on the real shape and size. In our study, the radius of the thin plate, a, is taken as 10 m, and the thickness of the plate, t, is taken as 0.2 m. The peripheral of the circular plate is fully fixed and the load magnitudes take different values, ranging from 10 kPa to 200 kPa, with an interval of 10 kPa. The given values in the numerical simulation are listed in

Table 1. Given values in numerical simulation.

Physical Quantities	Taken Values
plate radius a	10 m
plate thickness t	0.2 m
neutral layer modulus E0	2 × 1010 Pa
load magnitudes q	10 kPa to 200 kPa
tensile grade index α1	0.5
compressive grade index α2	0.1
tensile Poisson’s ratio μ+	0.35
compressive Poisson’s ratio μ−	0.25

. A three-dimensional solid element with eight nodes, C3D8, is adopted to conduct the numerical computation. Figure 4 shows the grid division, the loading, and the boundary conditions.

In our numerical computation, the simulation for bimodular functionally graded materials presents a slight degree of complexity. In ABAQUS6.14.4 software, the realization for functionally graded properties of materials is by layering, but before layering, in consideration of the bimodular effect, we must determine the position of the neutral layer first, that is, only after the subareas are in tension and compression can we effectively layer. To this end, we need to use Equation (12a) or (12b) and the data from Table 1 to determine the tensile thickness and compressive one first, which give t₁ = 0.4917t and t₂ = 0.5083t, respectively, where t is the thickness of the plate. Then, the thin circular plate is divided into eight layers along the thickness direction (see Figure 5), taking the middle modulus as the average modulus of this layer. The moduli of elasticity and Poisson’s ratios for these eight layers are computed and listed in

Table 2. Modulus of elasticity of 8-layer plate (t is the plate thickness, m).

Distance from Plate Top (m)	Modulus of Elasticity (×1010 Pa)	Poisson’s Ratio
0.0625t	1.913	0.25
0.1875t	1.937	0.25
0.3125t	1.961	0.25
0.4375t	1.986	0.25
0.5625t	2.055	0.35
0.6875t	2.186	0.35
0.8125t	2.326	0.35
0.9375t	2.475	0.35

, which are used for the property module during the material editing.

Figure 6 shows some representative displacement nephograms under different magnitudes of load, including 20 kPa, 40 kPa, 60 kPa, 80 kPa, 100 kPa, and 120 kPa. From Figure 6, it is easy to see that the maximum deflection occurs at the center of the thin circular plate, as we predicted; the central deflection gradually increases with the increase in load magnitude; and from the center of the plate to the edge of the plate, the deflection gradually decreases, as predicted in our theoretical solution.

Table 3. Numerical results of central deflection of three solutions.

q (kPa)	Central Deflection w0 (m)
	Result from Analytical Calculations		Result from [26]	Result from FEM
10	0.0898 1	0.0897 2	0.0895	0.0885
20	0.1538	0.1514	0.1516	0.1491
30	0.2003	0.1953	0.1963	0.1925
40	0.2368	0.2293	0.2311	0.2263
50	0.2670	0.2571	0.2602	0.2542
60	0.2931	0.2808	0.2851	0.2781
70	0.3160	0.3015	0.3070	0.2991
80	0.3366	0.3200	0.3268	0.3179
90	0.3554	0.3367	0.3447	0.3350
100	0.3727	0.3519	0.3613	0.3507
110	0.3887	0.3660	0.3766	0.3653
120	0.4037	0.3790	0.3910	0.3789
130	0.4178	0.3913	0.4045	0.3917
140	0.4312	0.4028	0.4173	0.4038
150	0.4438	0.4137	0.4294	0.4152
160	0.4558	0.4240	0.4409	0.4262
170	0.4674	0.4338	0.4519	0.4366
180	0.4784	0.4432	0.4625	0.4466
190	0.4890	0.4522	0.4726	0.4562
200	0.4992	0.4608	0.4824	0.4654

¹ Results from the variational solution, in which A₀ and A₁ are taken in Equation (40) and C₀ is taken in Equation (41) (please refer to Section 4.2). ² Results from the variational solution, in which A₀, A₁, and A₂ are taken in Equation (40), and C₀ and C₁ are taken in Equation (41) (please refer to Appendix A).

lists the central deflection values under different load magnitudes (from q = 10 kPa to q = 200 kPa, with an interval of 10 kPa). For an effective comparison, in Table 3, we also list two other groups of value from different theoretical solutions, the variational solution in this study and the perturbation solution in our previous study [26], in which the results from the variational solution are obtained via Equation (52) in this study; the results from the perturbation solution are based on Equation (104) in [26]. By comparing the values of central deflection in Table 3, it is easily found that the values from the three solutions are basically consistent, but there still exist small differences between them. However, the differences are negligibly small and generally acceptable, which verifies the validity of the variational method.

6. Results and Discussion

6.1. Numerical Comparision of Three Solutions

From the current results, it is found that if the numerical solution can be regarded as a standard to test the validity of two theoretical solutions, the perturbation solution is closest, followed by the variational solution. The values from the variational solution are slightly larger than the values from the perturbation solution, which agrees with the results of classical large-deflection solutions from [29,30]. As indicated in the Introduction, in the literature, both Chien’s perturbation solution [29] and the variational solution [30] give satisfactory results. In Chien’s perturbation solution [29], the relationship between load and maximum deflection gives

$$\frac{q{a}^4}{64D}=w_0\left[1+0.544{\left(\frac{w_0}{t}\right)}^2\right],$$

(54)

while, in the variational solution [30], the same relationship is

$$\frac{q{a}^4}{64D}=w_0\left[1+0.486{\left(\frac{w_0}{t}\right)}^2\right],$$

(55)

where a is the radius of a thin circular plate, t is the plate thickness, w₀ is the central deflection, q is the uniformly distributed loads, and D is the bending stiffness of the plate. It is readily found that they are quite close. In addition, we note that if other quantities in this relationship take the same values, the value of the maximum deflection from the variational solution is a little greater than the value from the perturbation solution.

The above conclusion is drawn on the basis of the variation solution with fewer terms. If the variation solution with more terms is taken, another phenomenon should be noticed. From the data of the second column marked with footer 2, it is easy to see that if more terms in the displacement functions are taken, the precision of the variational solution will be significantly improved, even exceeding that of the perturbation solution.

It should be noted here that the observed discrepancies resulting from the variational method and FEM method may come down to the fact that, apart from differences in the calculation methods themselves, it comes from different simulated models of the material properties. In the Ritz method, the bimodular functionally graded properties are considered as, for the bimodular effect, the subarea in tension and compression is used, while for the functionally graded property, in each tensile or compressive area, two smooth and continuous functions ((Equation (11)) are adopted. At the same time, in the FEM method, the subarea in tension and compression is still used for the bimodular effect, but the functionally graded property is realized by the layering along the direction of plate thickness; thus, the difference is inevitable. We can speculate that if more layers are adopted, the simulation of materials is likely to be closer to the continuous function change, like Equation (11). However, considering the computational efforts and time, only eight layers were adopted in our present study. In future work, more layers can be adopted to obtain a more precise result.

6.2. Stress Variation along Plate Thickness

In order to investigate the influence of plate thickness on the radial and circumferential stresses, we take three different thickness values, 0.1 m, 0.2 m, and 0.3 m, to carry out the numerical computation, and other taken values may refer to Table 1, in which the load intensity is taken as 10 kPa. At the same time, we take two different survey locations, that is, r = a/4 and r = 3a/4, where a is the radius of the circular plate, to investigate the influence of different radial locations on stresses. The numerical results are plotted in Figure 7, in which Figure 7a–c correspond to the thickness cases of t = 0.1 m, t = 0.2 m, and t = 0.3 m.

From Figure 7, the following two trends may be found. (i) The stress distribution trend of radial and circumferential stresses under different plate thicknesses is basically the same, that is, when r = a/4 (near the plate center), the differences between the radial and circumferential stresses tend to be smaller; when r = 3a/4 (near the plate edge), the differences between the radial and circumferential stresses tend to be larger, which may be due to the influence of boundary constraints. (ii) When the plate becomes thinner (t = 0.1 m), the radial and circumferential stresses both tend to be tensile, indicating that the membrane stress is dominant in thinner plates; when the plate becomes thicker (t = 0.3 m), the tensile area and compressive area appear distinct, showing that the bending stress is dominant in thicker plates. This phenomenon is also consistent with our expectations.

7. Conclusions

In this paper, the displacement variational method is used to solve the large-deformation problem of bimodular functionally graded thin plates. In order to facilitate the application of the variational method, the physical equations of the bimodular functionally graded material and the geometric equation under large deformation are first given. The total strain potential energy is expressed as the displacement component, which opens up the possibility for the realization of the Ritz method. Finally, the analytical result is verified by numerical simulation. The following three conclusions can be drawn.

(i)

The numerical simulation results verify the validity of the perturbation solution obtained in our previous study and the variational solution presented in this study.

(ii)

The perturbation method and variational method are both, in terms of nature, theoretical, being able to give useful analytical expressions that are convenient for use in the analysis and design. However, the variational method based on the energy principle avoids the establishment of an equation of equilibrium, which is necessary in the perturbation method yet.

(iii)

Compared with the traditional variational method, the improvement on this method in this study lies mainly in such a fact that the derivation of total strain energy is somewhat complicated due to the introduction of bimodular functionally graded materials and structural large deformation. In addition, the bending stiffness of the bimodular FGM plate may also be obtained from the derivation of total strain energy, but not necessarily from the conditions of equilibrium.

The results presented in this study are helpful for the refined analysis and optimized design of flexible thin plate structures, which are composed of functionally graded materials, while at the same time, the bimodular effect of materials is relatively obvious and cannot be ignored.

In the end, it should be pointed out again that Ambartsumyan’s bimodular model is established on the criterion of positive–negative signs of principal stresses. This fact makes it very difficult to use this model in structural analysis, because, except for a very few cases, the state of principal stress at any point in the structure is different each other under the action of external load. Fortunately, the proposal of a simplified mechanical model on a subarea in tension and compression makes it possible to use Ambartsumyan’s bimodular model in structural analysis. While, at the same time, the material model also deviates from the original definition, this may be seen as an imperfection of the method. In the future, we will try to establish a simplified mechanical model that satisfies the requirements of structural analysis and is closer to the original material model. This work is in progress.