Bending Analysis of Functionally Graded Axisymmetric Circular Plates using the Dual Mesh Finite Domain Method

https://doi.org/10.1590/1679-78256218 Abstract The dual mesh finite domain method (DMFDM) introduced by Reddy employs one mesh for the approximation of the primary variables (primal mesh) and another mesh for the satisfaction of the governing equations (dual mesh). The present study deals with the extension and application of the DMFDM to functionally graded circular plates under axisymmetric conditions. The formulation makes use of the traditional finite element interpolation of the primary variables with a primal mesh and a dual mesh to satisfy the integral form of the governing differential equations, the basic premise of the finite volume method. The method is used to analyze axisymmetric bending of through-thickness functionally graded circular plates using the classical plate theory (CPT) and first-order shear deformation plate theory (FST). The displacement model of the FST and the mixed model of the CPT using the DMFDM are developed along with the displacement and mixed finite element models. Numerical results are presented to illustrate the methodology and a comparison of the generalized displacements and bending moments computed with those of the corresponding finite element models. The influence of the extensional-bending coupling stiffness (due to the through-thickness grading of the material) on the deflections is also brought


Functionally graded structures
Structures whose material properties are continuously varied in one or more coordinate directions are abundant in nature (e.g., sea shells, bones, etc). With increasing demand for novel materials whose properties can be tailored to suit a particular application, such natural structures provided an inspiration for the design and development of an advanced class of composite materials called functionally graded materials (FGMs) (see, e.g., [1,2]). This is usually done by providing in-depth graded compositions, microstructure and properties [3]. With advent of advanced manufacturing and material processing techniques, the applications of FGMs have only increased [4]. For example, FGMs find applications in space structures as thermal barriers and in turbine rotors, flywheels, and gears as thin coatings, among other applications. While any number of materials can be used to achieve functional gradation, two-constituent functionally graded materials are the most widely used in engineering applications.
The variation in material properties of FGM can be expressed by a thickness function describing the spatial variation of volume fractions of FGM constituents. Various such functions are proposed in the literature (see, e.g., [5] for a review). For structural elements like beams and plates, the power-law variation of the modulus of elasticity, while keeping the Poisson's ratio constant, can be used. If the -coordinate is taken along the thickness of the plate, the variation of Young's modulus is given by (see Reddy [6]) where and are the material properties of the top and bottom faces of the plate, respectively, is the power-law index, and is the plate thickness. Note that when = 0, we obtain the single-material plate (with modulus ).

Dual mesh finite domain method (DMFDM)
Recently, Reddy [7] introduced a numerical approach termed the dual mesh finite domain method (DMFDM) for the solution of second-order differential equations in one and two dimensions with a single unknown. In the DMFDM, the domain of the governing equations is represented with a primal mesh of finite elements to introduce an approximation of the dependent variables, and a dual mesh of finite domains is superimposed on the primal mesh such that the finite element nodes are at the center of the finite domains, except for the nodes on the boundary, where half or quarter finite domains cover the boundary (see [7] for the details). Then the governing equations are satisfied in an integral sense over the finite (control) domain, as in the finite volume method (FVM). The second-order terms in the differential equation are integrated by parts (and thereby weakening the differentiability) and expressed as quantities (dual variables) on the interfaces of the dual mesh. Thus, the DMFDM can be viewed as a hybrid method that makes use of two best features of the FEM [8], namely, (a) the interpolation of the variables and (b) imposition of physical boundary conditions, and two salient features of the FVM [9]: (a) satisfaction of the global balance equations over the finite domain and (b) computation of the secondary variables at the boundaries of the finite domains where they are uniquely defined.

Present study
In [7], the DMFDM was introduced for second-order differential equations involving a single variable, with applications to heat transfer-like problems in one and two dimensions (also see [10]). Recently, Reddy and Nampally [11] and Reddy, Nampally, and Srinivasa [12] extended the method to linear and nonlinear analysis, respectively, of beams with multiple variables (i.e., to systems of coupled second-order differential equations in several dependent unknowns). In this paper, the method is extended to the bending of through-thickness functionally graded circular plates under axisymmetric loads and boundary conditions. Numerical results are presented to illustrate the accuracy of the DMFDM in predicting the deflections and stresses when compared with the finite element solutions. The effect of the extensionalbending stiffness on the deflections is also elucidated.
Following this introduction, a review of the governing equations of functionally graded axisymmetric circulate plates is given in section 2. A reformulation of the governing equations of circular plates (see Reddy [13,14] for the derivation of the governing equations) as applied to functionally graded plates is presented in section 3. For the mixed model of the CPT, the governing equations are expressed in terms of the generalized displacements and bending moment. The dual mesh finite domain discretization of these equations are presented in section 4. Although the formulations are simple, they are not readily available in the literature. Numerical results are presented in section 5 for hinged and clamped boundary conditions, and the numerical results obtained with the FEM and DMFDM models are compared with the exact Latin American Journal of Solids and Structures, 2020, 17(7), e302 3/24 solutions to assess the relative accuracy in predicting the generalized displacements and generalized forces. Finally, some remarks about DMFDM and its extensions are outlined in section 6.

Governing equations
Consider a through-thickness functionally graded circular plate of thickness and radius subjected to axisymmetric distributed load (i.e., independent of the angular coordinate, ) on the top face. If further the boundary conditions are also selected to be axisymmetric, then the plate can be considered as functionally graded axisymmetric circular plate. We select the cylindrical coordinate system ( , , ) to analyze the plate, where is the radial coordinate outward form the center of the plate (0 ≤ ≤ ), denotes the transverse coordinate (− /2 ≤ ≤ /2), and is the angular coordinate (0 ≤ ≤ 2 ) (see figure 1). The displacement field of such a circular plate based on the firstorder shear deformation plate theory (FST) are given by The non-zero linear strains based on the above displacements would then be while the constitutive relations of the functionally graded axisymmetric circular plate are given by We further define the stress and moment resultants on the circular plate (see figure 1) as follows: Latin American Journal of Solids and Structures, 2020, 17(7), e302 4/24 The explicit expressions of the above integrals for the power-law variation of ( ) (see Eq. (1.1)) are given as The equations of equilibrium of the functionally graded (FG) axisymmetric circular plate can be obtained using the principle of virtual work. The virtual work statement for FG axisymmetric circular plate is Now using the stress and moment resultant definitions of Eqs. (2.4a)-(2.4e) and strain definition of Eqs. (2.2a)-(2.2c) we can rewrite the virtual work statement as The governing differential equations of the FG axisymmetric circular plate are obtained by taking the Euler-Lagrange equations of the above variational statement. The resulting governing equations read The governing equations derived above are based on the first-order shear deformation theory, where the transverse shear strain is non-zero and assumes a constant value (and hence require the shear correction factor, in Eq. (2.5)). However, as the radius to thickness ratio of the circular plate increases the shear strains are known to tend to zero. In such cases we may also use classical plate theory (i.e., based on Kirchhoff hypothesis) to predict the plate behavior. To The stress and moment resultants on the FG axisymmetric circular plate based on the classical plate theory would then be In this case the transverse shear strain vanishes; however, the governing equations remain the same as those listed in Eqs. (2.9)-(2.11) with the exception that the shear stress resultant can only be expressed in terms of the moments and as given by Eq. (2.11) rather than directly computing in terms of shear strain (which is zero in this case). The governing differential equations of FG axisymmetric circular plate based on the FST can be expressed in terms of the generalized displacements ( , , ) by substituting the expressions of stress and moment resultants of Eqs. (2.4a)-(2.4e) into Eqs. (2.9)-(2.11). The resulting equations are listed below.
Similarly, the governing differential equations of FG axisymmetric circular plate based on CPT can be expressed in terms of the displacements and by first expressing Eq. (2.10) in terms of moments and using Eq. (2.11) as and then substituting the stress and moment resultants of Eqs. (2.13a)-(2.13d) into the above equation and Eq. (2.9) to give Latin American Journal of Solids and Structures, 2020, 17 (7), e302 6/24

Classical plate theory
The governing equations of FG axisymmetric circular plate based on the classical plate theory when expressed in terms of the displacements, and , would result in fourth order differential equations. However, since the dual mesh finite domain method is applicable only for first-order or second-order differential equations, we will recast the governing equations (2.9)-(2.11) as second-order differential equations such that the primary variables of the resulting equations are { , , }. First, we will consider the equations (2.13a) and (2.13b). From these two equations we can write an expression of in terms of and . Similarly, can be expressed in terms of and using the equations (2.13c) and (2.13d). In the following we list the resulting equations after some algebraic manipulations.
To eliminate from Eq. (3.2) we use Eq. (2.13b) and Eq. (2.13d) to obtain Thus can be written as Finally, Eq. (2.13b) can be expressed as The governing equations of FG axisymmetric circular plate based on classical plate theory in terms of { , , } would then be given by Eq.

Dual mesh finite domain formulations 4.1 Mixed DMFDM for CPT
As noted earlier, the dual mesh finite domain method is only applicable to equations of second order or less. Thus we develop dual mesh finite domain formulation of the FG axisymmetric circular plate based on classical plate theory using the governing equations (3.6)-(3.8). The resulting formulation is called mixed dual mesh finite domain method (see [11]) since the formulation has both generalized displacements ( , ) and the moment ( ) as the primary variables. This is similar to the mixed formulations of finite element method. To obtain the discretized equations, we consider the computational domain = (0, ) and divide it into finite elements (primal mesh) such that each node has its associated finite domain (dual mesh). Except for the finite domains corresponding to the boundary nodes, all the interior finite domains encompass two finite elements (primal mesh) such that one half of each of these two finite elements lie within the finite domain (see figure 2). Thus, for a uniform primal mesh the interior nodes lie at the center of their corresponding finite domains.
The dependent variables are approximated on each finite element using Lagrange interpolation functions (see figure 3). Although different sets of interpolation functions can be used for different primary variables, in the present paper we use same set of linear 1-D Lagrange interpolation functions for all the primary variables. For a typical finite element ( ) = ( , +1 ), the primary variables can thus be approximated using the set of linear 1-D Lagrange Here , , represent the values of , and respectively at node , while +1 , +1 , +1 represent the values of , and respectively at node + 1. linear Lagrange interpolation is needed on [8]. However, for the mixed formulations, akin to the one considered here, a minimum of linear Lagrange interpolations on all the primary variables would suffice [11]. The linear Lagrange interpolation functions 1 ( ) and 2 ( ) on a typical finite element ( ) are given by  Now we proceed to develop the discretized equations corresponding to Eqs. (3.6)-(3.8) by writing their integral statements on a typical interior finite domain ( ( ) , ( ) ), associated with node (see figure 4).
The first step in developing the discretized equations is to carry out the indicated integration of Eqs.
In the above equations the secondary variables are defined as follows  .7) can be either evaluated numerically (see, e.g., Gauss quadrature rule) or exactly. In the present paper we use Gauss quadrature to evaluate the integral expressions. Thus for an interior node associated with interior finite domain numbered the discretized equations would read Latin American Journal of Solids and Structures, 2020, 17 (7), e302 10/24 In the above equations −1 represents the value of in the finite element ( −1) and represents the value of in the finite element ( ) . Similar meaning applies to other coefficients , , and . In a similar fashion, the discretized equations for the finite domain 1 (see figure 5) would be Finally, the discretized equations on ( + 1) th finite domain (see figure 5) would be This completes the mixed dual mesh finite domain formulation of FG axisymmetric circular plate based on CPT.

Displacement DMFDM for FST
Since the governing differential equations of axisymmetric circular plate based on FST in terms of generalized displacements are second-order differential equations, we can developed dual mesh finite domain method directly from the equations (2.14)-(2.16). The integral statements of the governing equations on a typical interior finite domain ( ( ) , ( ) ), associated with node , would be The indicated integration in the above equations is carried out such that the resulting boundary terms constitute the secondary variable, dual to the primary variable of the equation considered.

Numerical Results
To illustrate the workings of the dual mesh finite domain models presented in the previous sections, we consider two cases: (a) hinged (or simply supported) and (b) clamped functionally graded circular plates, subjected to uniformly distributed load of intensity . A comparison between the numerical results obtained with FEM (see Appendix B for various models) and DMFDM models for each case considered is presented. There are three models of FEM and two models of DMFDM, as summarized below: • We investigate the effect of mesh and power-law index on the deflections and stresses. We consider a functionally graded circular plate of radius = 10 in, thickness = 0.1 in and subjected to axisymmetric uniformly distributed load intensity = 0.5 lb/in 2 . The functionally graded material properties are taken to be = 3 × 10 7 psi, = 3 × 10 6 psi, and = 0.3 (units are stated here only to give a physical meaning; one may replace them with appropriate metric units).
The stresses and stress resultants can be post-computed in the numerical models described. For the displacement models, once the generalized displacements are obtained, their derivatives can be easily determined. Then the stresses can be directly computed from the constitutive relations of Eq. (2.3), while the stress resultants and moments can be computed using Eqs. (2.4a)-(2.4) or Eqs. (2.13a)-(2.13d) depending on whether FST is used or CPT is used. For the mixed models, the second derivative of the transverse deflection needs to be calculated from the moment , which is determined as a nodal variable in mixed models. The second derivative of transverse deflection can be obtained from the following equation: Then the stresses can be computed using Eqs. (2.3) while the stress resultants can be computed using Eqs. (3.1) and (3.4).

Hinged plates
The boundary conditions for the primary variables of hinged axisymmetric circular plate in various models are as follows: Latin American Journal of Solids and Structures, 2020, 17 (7) A comparison between the transverse deflections at the center of a hinged homogeneous circular plate obtained from various models is presented in Table 1 When 32 elements are used, all models give the same results as the exact solution (see [13][14][15]) up to the fourth decimal point. The mixed finite element model has a slow mesh convergence rate at the origin for transverse deflection, , and hence a non-uniform finite element mesh with a finer mesh at the origin is used in reporting the values of (0). The non-uniform mesh is chosen such that for two element mesh, the lengths of the elements are in ratio 1:9, with the shortest element being at origin. All the subsequent refinements are made by breaking the finite elements into half. Thus for four element mesh the element lengths would be {0.5 , 0.5 , 4.5 , 4.5 } and further refinement is made by breaking these elements into their halves and so on. Since the radius to thickness ratio / = 100, both classical plate theory and first order shear deformation theory predict the same results. Figure 7(a) contains a comparison of the moment, , obtained from displacement and mixed finite element models and mixed dual mesh finite domain model of the CPT when a uniform mesh of 32 elements is used, while Fig. 7(b) includes a comparison of the moment, , obtained from displacement finite element model and displacement dual mesh finite domain model of the FST when a uniform mesh of 32 elements is used. While the moments obtained from all the models are very close, the mixed finite element model and mixed dual mesh finite domain model deviate from the displacement finite element model at the origin ( = 0). Both mixed models (i.e., mixed FEM and mixed DMFDM models) are prone to give erroneous moments at the origin. It should be noted that this anomalous behavior at the origin is not a characteristic of the dual mesh finite domain method; both the mixed finite element model and mixed dual mesh finite domain model exhibit such behavior. However, the displacements models do not exhibit such a behavior. Table 1: Center transverse deflection, (0), of hinged homogeneous axisymmetric circular plate as predicted by various models.
The results are given in inches.

Clamped plates
The boundary conditions on the primary variables of clamped axisymmetric circular plate in various models are as follows: A comparison of the stress values, , on the top face of the clamped homogeneous circular plate as predicted by various numerical models based on CPT with the analytical solution [14] are presented in Fig. 8. Although all the models are in good agreement with the analytical solution for most part of the computational domain, the mixed models (i.e., mixed FEM model and mixed DMFDM model) deviate from the analytical solution at and very close to the origin. This is due to the erroneous prediction of moment near the origin by the mixed models. Further, Fig. 9 contains a comparison of stress on the top face of clamped homogeneous circular plate as predicted by various numerical models based on the FST with that of the analytical solution [13,14]. Since both the finite element model and dual mesh finite domain model of the FST are displacement based, the stress values are in good agreement with the analytical solution through out the computational domain.
Finally, to illustrate the working of the dual mesh finite domain method in analyzing the functionally graded axisymmetric circular plates, we compare the transverse deflection at the center, (0), from the analytical solution of the CPT [13][14][15] with the transverse deflections obtained from various models described earlier. Various values of the power law index, , are considered to bring out the effect of on the bending deflections of the plate. Table 2 contains a comparison between the various models with that of the analytical solution (of the FST). For this thin plate ( / = 100), both the CPT and the FST predict solutions close to each other [13][14][15]; 32 elements are used for the results reported in Table 2. Uniform mesh is used in all models except for the mixed FEM model of the CPT. A non-uniform mesh described earlier is used for FE-CP(M) to achieve faster mesh convergence. It can be seen that as the power-law index, , increases the plate stiffness decreases (since plate material tends toward < ).