Nonlinear Finite Element Model for Bending Analysis of Functionally-Graded Porous Circular/Annular Micro-Plates under Thermomechanical Loads Using Quasi-3D Reddy Third-Order Plate Theory

Abstract

A nonlinear finite element model for axisymmetric bending of micro circular/annular plates under thermal and mechanical loading was developed using quasi-3D Reddy third-order shear deformation theory. The developed finite element model accounts for a variation of material constituents utilizing a power-law distribution of a two-constituent material, three different porosity distributions through plate thickness, and geometrical nonlinearity. The modified couple stress theory was utilized to account for the strain gradient effects using a single material length scale parameter. Three different types of porosity distributions that have the same overall volume fraction but different enhanced areas were considered as a form of cosine functions. The effects of the material and porosity distribution, microstructure-dependency, the geometric nonlinearity, and various boundary conditions on the bending response of functionally-graded porous axisymmetric microplates under thermomechanical loads were studied using the developed nonlinear finite element model.

1. Introduction

Functionally-graded materials (FGMs) are advanced engineering materials composed of two or more constituents with a continuous variation in their compositions. Unlike FGMs, laminated composites exhibit immediate changes in thermal and mechanical properties of the constituents, resulting in stress concentrations at the interfaces where two discrete materials bond together. This leads to delamination problems and the presence of residual stresses in conventional composites working under severe conditions. FGMs were developed by researchers in Japan in 1984 to overcome these issues encountered in a thermal coating material requirement of a hypersonic space plane project [1]. Since then, FGMs have been used in various fields such as aerospace, automobile, electronic, and medical industries due to their advantages over laminated composites and their flexibility to be designed according to the needs of the application field and working environment. The reader is referred to the following review articles [2,3,4,5] for details of the historical development of these materials, manufacturing techniques, and optimization of their functionality.

The FGMs have great potential for improving the performance of various components in engineering structures, especially circular and annular plates. Over the last few decades, researchers have extensively studied the behavior of functionally-graded (FG) circular and annular plates under thermal, mechanical, and combined thermomechanical loadings. Since FGMs were initially designed to withstand extreme thermal environments, most of the literature focuses on their thermal analysis. Typical FGMs used in these studies are made from a mixture of ceramics for their low thermal conductivity and metals for their ductility and resistance to fracture caused by stresses likely to occur in high-temperature gradients. Additionally, the majority of studies on FG plates employ a power law or exponential distribution of materials through the thickness direction of the plates

In 1998, Reddy and Chin [6] conducted a numerical study to investigate thermomechanical responses of FG cylinders and plates under extreme thermal loading conditions using the first-order shear deformation plate theory (FSDT). In their study, the effects of thermomechanical coupling on the response of FGMs subjected to thermal shock were investigated. For the functionally-graded axisymmetric cylinder subjected to high thermal loading, the temperature distribution obtained from both coupled and uncoupled formulations did not show significant differences. However, it was observed that the radial stresses were more affected than the hoop stresses in the FG cylinder.

Using the FSDT, exact solutions of the static bending analysis of FG circular and annular plates having various boundary conditions were presented by Reddy et al. [7]. They derived the solutions of deflections, forces, and the moments of the FG plates based on FSDT in terms of the associated quantities for the isotropic plates based on the classical plate theory (CPT). Hence, the bending solutions of the FG circular plate became readily available whenever the CPT solution was known. Ma and Wang [8] studied the axisymmetric nonlinear bending and post buckling response of functionally-graded circular plates under thermal, mechanical and combined thermomechanical loading conditions. In this study, governing equations were derived using the von Kármán plate theory and the numerical solutions were obtained with the help of the shooting method. The results of this study showed that temperature distribution, deflection values, critical buckling temperature, and post buckling behavior of the functionally-graded circular plates were significantly affected by the volume fraction index.

Praveen and Reddy [9] introduced the finite element formulations that account for the transverse shear strains, rotary inertia, and von Kármán nonlinear strains to perform static and dynamic thermoelastic analysis of the functionally-graded ceramic-metal plates based on FSDT. In 2000, Reddy [10] presented the formulation and analytical solution of simply-supported rectangular FG plates using third order shear deformation plate theory (TSDT) including thermomechanical coupling, time dependency, and von Kármán geometric nonlinearity. From these two studies, it was concluded that the distribution of material constituents in the functionally-graded plates had a significant influence on the resulting thermoelastic response of FG plates. Najafizadeh and Heydari [11] investigated the thermal buckling analysis of functionally-graded circular plates under both uniform and non-uniform temperature changes by employing the TSDT.

Prakash and Ganapathi [12] employed the finite element method to carry out asymmetric-free vibration and thermoelastic stability analysis of functionally-graded circular plates. Nie and Zhong [13] studied the three-dimensional free and forced vibration analysis of functionally-graded circular plates and found that the lowest nondimensional frequency and circumferential wave number of the plate increased as the thickness-to-width ratio increased. They also observed that the magnitudes of the displacements and stresses became larger as the forcing frequency approached the natural frequency of the FG circular plate.

Efraim and Eisenberger [14] presented the free vibration analysis of variable thickness thick annular plates using the exact element method and the dynamic stiffness method. They used FSDT in their formulations and varied Poisson ratio according to the power law distribution in addition to elastic modulus and mass density. Golmakani and Kadkhodayan [15] presented another study that accounted for the gradation of Poisson ratio. They investigated the nonlinear bending analysis of annular FG plates based on both FSDT and TSDT. The same authors [16] later performed a large deflection analysis of circular and annular FG plates subjected to thermomechanical loading within the framework of FSDT, including von Kármán nonlinearity.

Saidi et al. [17] employed unconstrained third order shear deformation theory to analyze the axisymmetric bending and buckling behavior of thick FG circular plates.

Nosier and Fallah [18] reformulated governing equations of the FSDT into interior and edge-zone equations for functionally-graded circular plates. By introducing two sets of equations to define the edge-zone problem, they uncoupled the bending and extension equations, which made it possible to obtain analytical solutions for the asymmetric behavior of functionally-graded circular plates with various boundary conditions under mechanical and thermal loading. Later, they included the von Kármán nonlinear strains into their formulations and investigated the axisymmetric and asymmetric nonlinear bending of functionally-graded circular plates subject to linearly-varying transverse loading [19]. The axisymmetric bending analysis of FG circular plates under arbitrary transverse loads was studied by Yun et al. [20]. They obtained the analytical solutions for the FG circular plates with elastic simple and rigid slipping supports cases when the material property of the FG plate was varying with an exponential distribution. Another analytical study was conducted to solve for in-plane and out-of-plane free vibrations of thick FG circular and annular plates embedded in piezoelectric layers by Talabi and Saidi [21], employing TSDT. The effects of both electrical and mechanical boundary conditions, geometrical parameters of the plate, and in-plane displacements on the middle plane on the natural frequencies of FG circular and annular plates were discussed. Żur [22] applied the Neumann series method to investigate the free vibration behavior of discrete-continuous FG circular plates that may have several ring attachments such as masses, springs and damping elements.

The FG circular and annular plates can be further improved by adding porosity into their composition to decrease the weight of the structure and/or increase the insulation properties. Hence, it is important to examine the mechanical and thermal responses of FG porous plates under different loading and boundary conditions. A general solution of a porous FG circular plate that is supported by a non-uniform Kerr elastic foundation and subjected to non-axisymmetric, non-uniform shear and normal tractions, and a magnetic actuation was developed by Rad and Shariyat [23]. Their results showed that the radial displacement component was more prone to being affected by the induced magnetic actuation. Additionally, because of the presence of incompressible fluid in the pores in this study, as the porosity increased, the plates became stiffer. The buckling behavior of porous circular plate between piezoelectric layers under thermal loading was investigated by Jabbari et al. [24]. They showed that, as the porosity increased, the critical temperature decreased and the plate whose pores were saturated with fluid became unstable. On the other hand, the critical temperature of the plates can be decreased by increasing the thermal expansion coefficient of the fluid filling the pores and the piezoelectric layers. Zhao et al. [25] studied the free and forced vibration analysis of FG porous circular, annular, and sector plates with general elastic restraints using FSDT.

These extensive studies conducted on FG circular and annular plates show that these structures have an intrinsic advantage resulting from the non-homogeneity and smooth variations of the material properties. It is shown that the deflections and tensile stresses of FG circular and annular plates can be lower and critical buckling loads can be higher as compared to the homogeneous ones, depending on the predetermined variation of material properties of FG circular and annular plates. It is also possible to adjust the natural frequencies of these structures by changing the variation of the material distribution. Hence, all these conclusions make it attractive to examine the performance of the FGMs for the micro-scale structures. However, conventional continuum mechanics cannot capture the size dependency that is experimentally observed at the micro-scale [26,27,28,29]. Therefore, a higher order continuum theory is required for the accurate modeling and analysis of these structures. Couple stress theories [30,31,32], Erigen nonlocal elasticity theory [33] and the strain gradient elasticity theories [26,34,35] are some of the higher order continuum theories that take the size dependency into account. The modified couple stress theory is the most commonly employed theory because only a single length scale parameter is needed to include size effect.

Ke et al. [36] investigated the bending, buckling, and free vibration analyses of FG annular microplates with hinged–hinged and clamped–clamped boundary conditions. Their size-dependent annular microplate model was based on the Mindlin plate theory and the modified couple stress theory. This study showed that elastic buckling analysis was more sensitive to size effect than the free vibration analysis. Similar analyses were presented by Ansari et al. [37] for FG circular and annular microplates. They also employed Mindlin plate theory, but different to the previous study, size dependency was included using modified strain gradient elasticity theory. Both studies agreed that the smaller the dimensionless length scale parameter they had, the smaller the deflection but the higher the critical buckling load and natural frequencies that they obtained.

Reddy and Berry [38] presented the classical and the first order plate theories for axisymmetric bending of circular micro-plates including von Kármán nonlinear strains. Size dependency was captured with the modified couple stress theory. Later, Reddy et al. [39] used this theory to develop nonlinear finite element models for FG circular plates.

An analytical solution for the free vibration of FG circular and annular nanoplates was obtained by Hosseini-Hashemi et al. [40] based on Mindlin plate theory and Eringen nonlocal elasticity theory. Beni et al. [41] studied the same problem for FG cylindrical nanoshells using FSDT in conjunction with the modified couple stress theory. They presented the effects of the material length scale, distribution of the FGMs, nanotube thickness, and length on the fundamental frequencies. Eshraghi et al. [42] studied the bending and the free vibration analysis of FG annular and circular microplates subjected to thermal loading using the modified couple stress theory. They unified the displacement fields such that results for Kirchoff plate theory, Mindlin plate theory, and third order shear deformation plate theory can be generated. Additionally, not only the mechanical and thermal properties of the FG plates but also the material length scale parameter were not kept constant but were changed through the thickness direction, obeying a power law distribution. The transverse deflections, normalized circumferential and radial stresses, and the natural frequencies were presented for different thermal loading, material, and geometrical parameters. Ji et al. [43] developed a plate model capturing the size dependency for FG circular micro-plates based on the strain gradient theory of Zhou. They analyzed the bending and free vibration of a simply-supported circular micro-plate and the results were compared with those obtained by employing the strain gradient theory of Lam, the modified couple stress theory, and the CPT.

A free vibration and thermal buckling analysis of an FG porous circular micro-plate was conducted by Shojaeefard et al. [44] based on CPT and FSDT with modified couple stress theory. The effects of the temperature change, distribution of the material properties, size-dependency, and porosity on the fundamental frequencies and critical temperature were investigated. Kim et al. [45] presented the analytical solutions of bending, free vibration, and the buckling problem for FG porous micro-plates using CPT and FSDT in conjunction with the modified couple stress theory. Recently, Wang and Zhang [46] studied the thermal buckling and postbuckling responses of GPL-reinforced nanocomposite beams using the higher order shear deformation theory with temperature-dependent properties. Zhang et al. [47] carried out analytical studies on thermo-mechanical responses of porous functionally-graded, graphene-reinforced cylindrical panels based on a third order shear deformation theory. The acoustic characteristics of functionally-graded porous graphene reinforced nano composite plates (FG-PGRC) were studied by Xu et al. [48]. In their study, a higher order shear deformation theory was utilized to study the vibration and noise reduction of an FG-PGRC plate.

This study aimed to investigate the behavior of FG porous circular microplates under thermal and mechanical loadings, which has not been studied in the literature. To this end, a nonlinear finite element model was developed based on quasi-3D Reddy third-order shear deformation theory and the modified couple stress theory, taking into account von Kármán nonlinear strains to consider geometrical nonlinearity. The FGM was composed of two constituents based on a power law distribution through the thickness direction, and three different porosity profiles were considered. Parametric analyses were conducted to investigate the effects of the distribution of material properties and porosity, size-dependency, geometric nonlinearity, and different boundary conditions on the static bending analysis of FG porous circular microplates.

2. Constitutive Models

2.1. Functionally-Graded Porous Materials

The model considers isotropic axisymmetric plates composed of two constituents with varying material properties and internal porosity through the thickness, modeled using a power-law distribution and cosine variation, respectively. The typical material properties of functionally-graded porous materials (FGPM) are thus captured in the model, as shown in Equation (1).

$$P\left(z\right)=\left[\left(P_t-P_b\right)f\left(z\right)+P_b\right]\left(1-\psi \left(z\right)\right),f\left(z\right)={\left(\frac{z}{h}+\frac{1}{2}\right)}^n,$$

(1)

where $P_t$ and $P_b$ are material properties on the top and bottom surfaces of plates, n is power-index, $f\left(z\right)$ is a volume fraction function, and $\psi \left(z\right)$ is a porosity distribution function. Three different types of porosity distributions were considered in this study.

$$\begin{array}{cc}\hfill \mathrm{Type}1:\psi \left(z\right)& =\varphi \cos \left[\pi \left(\frac{z}{h}\right)\right]\hfill \\ \hfill \mathrm{Type}2:\psi \left(z\right)& =\varphi \cos \left[\frac{\pi }{2}\left(\frac{z}{h}+0.5\right)\right]\hfill \\ \hfill \mathrm{Type}3:\psi \left(z\right)& =\varphi \cos \left[\frac{\pi }{2}\left(\frac{z}{h}-0.5\right)\right],\hfill \end{array}$$

(2)

where $\varphi$ is the maximum porosity value along thickness direction. The distribution of porosity through the thickness of the plates was normalized to have the same porous volume, and it is important to investigate the effect of different porosity distributions [45]. Figure 1a displays the normalized porosity distribution throughout the thickness of the plate. Figure 1b–d show the effects of porosity distributions on the variation of typical materials properties. As an example, a porosity value was set to $\varphi =0.5$, three different power-law index values $n=0,0.5$, and 5.0 were set. The ratio of material properties on the top and bottom surfaces was assumed to be $\frac{E_t}{E_b}=10$. The Type 1, Type 2, and Type 3 porosity distributions are a symmetric and center-enhanced, a bottom area enhanced, and a top area enhanced porosity distributions, respectively.

2.2. Modified Couple Stress Theory

The motion of the material particles in classical couple stress theory [30,49] is described to rotate the material particles in addition to forces in the classical continuum mechanics. The size-dependent effect was captured using two additional material constants in the classical couple stress theory. These two material constants are difficult to determine because of their indeterminacy. Eringen [33] proposed a micropolar theory and defined the motion of a particle using the location vector and inner product of a rigid vector. A modified couple stress theory using the concept of the representative volume elements and a higher order equilibrium condition was proposed by Yang et al. [32]. According to the modified couple stress theory, the deviatoric part of a couple stress tensor is only associated with the symmetric part of rotation gradient and it contributes to the total strain energy along with the classical strain energy. The strain energy potential of an axisymmetric plate based on the modified couple stress theory can be expressed as

$$\begin{array}{cc}\hfill U& =\frac{1}{2}{\int }_V\left(\mathit{\sigma }:\mathit{\epsilon }+\mathit{m}:\mathit{\chi }\right)dV\hfill \\ \hfill & =\frac{1}{2}{\int }_{r_i}^{r_o}\left[{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\left(\mathit{\sigma }:\mathit{\epsilon }+\mathit{m}:\mathit{\chi }\right)dz\right]rdr,\hfill \end{array}$$

(3)

where $r_i$ and $r_o$ are the inner and outer radii of the plate, $\mathit{\sigma }$ and $\mathit{\epsilon }$ are the Cauchy stress tensor and Von Kámán nonlinear strain tensor, $\mathit{m}$ and $\mathit{\chi }$ are the deviatoric part of the symmetric couple stress tensor and the symmetric curvature tensor. Note that the differential volume element $dV$ can be written as $dV=rdrd\theta dz$ and $2\pi$ from the integration with respect to $\theta$ being omitted in Equation (3). The curvature tensor and the deviatoric part of the symmetric couple stress tensor are defined as [32]

$$\begin{array}{cc}\hfill \mathit{\chi }& =\frac{1}{2}\left[\nabla \mathit{\omega }+{\left(\nabla \mathit{\omega }\right)}^{\mathrm{T}}\right]\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \mathit{m}& =2\mu {\ell }^2\mathit{\chi },\hfill \end{array}$$

(5)

where $\mathit{\omega }$ is the rotation vector, $\mathit{\omega }=\frac{1}{2}\nabla \times \mathit{u}$, $\mu$ is the shear modulus, and ℓ is a length scale parameter.

In this study, an isotropic linear elastic material was assumed and the stress and strain relation [50] for an axisymmetric plate is

$$\left\{\begin{array}{c}{\sigma }_{rr}\\ {\sigma }_{\theta \theta }\\ {\sigma }_{zz}\\ {\sigma }_{rz}\end{array}\right\}=\mathsf{\Lambda }\left[\begin{array}{cccc}1-\nu & \nu & \nu & 0\\ \nu & 1-\nu & \nu & 0\\ \nu & \nu & 1-\nu & 0\\ 0& 0& 0& \frac{1}{2}\left(1-2\nu \right)\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{rr}-\alpha \left(T-T_0\right)\\ {\epsilon }_{\theta \theta }-\alpha \left(T-T_0\right)\\ {\epsilon }_{zz}-\alpha \left(T-T_0\right)\\ {\epsilon }_{rz}\end{array}\right\},$$

(6)

where $\mathsf{\Lambda }=E/\left[\right(1+\nu \left)\right(1-2\nu \left)\right]$, E is Young’s modulus, which varies along the plate’s thickness, $\nu$ is a constant Poisson’s ratio in the elastic stiffness matrix. $\alpha$ is the thermal expansion coefficient, and T and $T_0$ are the temperature at a material point and the reference temperature of the undeformed body.

The nonzero curvatures and modified couple stresses are

$$\left\{\begin{array}{c}{m}_{r\theta }\\ m_{\theta z}\end{array}\right\}=2{\ell }^2\mu \left\{\begin{array}{c}{\chi }_{r\theta }\\ {\chi }_{\theta z}\end{array}\right\}.$$

(7)

3. Quasi-3D Reddy Third-Order Plate Theory

3.1. Displacement and Strains

The displacement field of quasi-3D Reddy third-order plate theory can be derived from an assumption of a cubic variation of in-plane displacements and a quadratic variation of deflection (i.e., out-of-plane displacement) with zero tangential traction on top and bottom surfaces. The displacement field of cubic variation of in-plane displacement and a quadratic variation of deflection through thickness direction for axisymmetric plates takes the form of

$$\begin{array}{cc}\hfill u_r\left(r,z,t\right)& =u_0\left(r,t\right)+z{\theta }_r\left(r,t\right)+z^2{\varphi }_r\left(r,t\right)+z^3{\psi }_r\left(r,t\right)\hfill \\ \hfill u_z\left(r,z,t\right)& =w_0\left(r,t\right)+z{\theta }_z\left(r,t\right)+z^2{\varphi }_z\left(r,t\right).\hfill \end{array}$$

(8)

With the assumption of zero tangential traction on top and bottom surfaces, the displacement (8) can be written in the form of

$${\epsilon }_{rz}\left(r,\frac{h}{2},t\right)={\epsilon }_{rz}\left(r,-\frac{h}{2},t\right)=0.$$

(9)

The form of quasi-3D Reddy third-order plate theory for axisymmetric plates takes

$$\begin{array}{cc}\hfill u_r\left(r,z,t\right)& =u_0\left(r,t\right)+z{\theta }_r\left(r,t\right)-\frac{z^2}{2}\frac{\partial {\theta }_z}{\partial r}-z^3{c}_1\left[{\theta }_r\left(r,t\right)+\frac{\partial \lambda \left(r,t\right)}{\partial r}\right]\hfill \\ \hfill u_z\left(r,z,t\right)& =w_0\left(r,t\right)+z{\theta }_z\left(r,t\right)+z^2{\varphi }_z\left(r,t\right),\hfill \end{array}$$

(10)

where $u_0$ is the membrane displacement, ${\theta }_r$ is the rotation of a transverse normal about $\theta$ direction, $w_0$ is the deflection, ${\theta }_z$ and ${\varphi }_z$ are the thickness stretch, $\lambda =w_0+\frac{h^2}{4}{\varphi }_z$ and $c_1=\frac{4}{3h^2}$.

Based on the assumption of small strains and moderate rotations, nonzero von Kámán nonlinear strain for the axisymmetric plate is given by [39].

$$\begin{array}{cc}\hfill {\epsilon }_{rr}& =\frac{\partial u_r}{\partial r}+\frac{1}{2}{\left(\frac{\partial u_z}{\partial r}\right)}^2\hfill \\ \hfill {\epsilon }_{\theta \theta }& =\frac{u_r}{r}\hfill \\ \hfill {\epsilon }_{zz}& =\frac{\partial u_z}{\partial z}\hfill \\ \hfill {\epsilon }_{rz}& =\frac{1}{2}\left(\frac{\partial u_r}{\partial z}+\frac{\partial u_z}{\partial r}\right).\hfill \end{array}$$

(11)

The non-zero strains with the displacement field (10) of quasi-3D Reddy third-order plate theory are defined as

$$\begin{array}{cc}\hfill \left\{\begin{array}{c}{\epsilon }_{rr}\\ {\epsilon }_{\theta \theta }\\ {\epsilon }_{zz}\\ 2{\epsilon }_{rz}\end{array}\right\}& =\left\{\begin{array}{c}\frac{\partial u_0}{\partial r}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial r}\right)}^2\\ \frac{u_0}{r}\\ {\theta }_z\\ {\theta }_r+\frac{\partial w_0}{\partial r}\end{array}\right\}+z\left\{\begin{array}{c}\frac{\partial {\theta }_r}{\partial r}\\ \frac{{\theta }_r}{r}\\ 2{\varphi }_z\\ 0\end{array}\right\}-z^2\left\{\begin{array}{c}\frac{1}{2}\frac{{\partial }^2{\theta }_z}{\partial r^2}\\ \frac{1}{2r}\frac{\partial {\theta }_z}{\partial r}\\ 0\\ -c_2\left({\theta }_r+\frac{\partial \lambda }{\partial r}\right)+\frac{\partial {\varphi }_z}{\partial r}\end{array}\right\}\hfill \\ \hfill & -z^3\left\{\begin{array}{c}{c}_1\left(\frac{\partial {\theta }_r}{\partial r}+\frac{{\partial }^2\lambda }{\partial r^2}\right)\\ \frac{c_1}{r}\left({\theta }_r+\frac{\partial \lambda }{\partial r}\right)\\ 0\\ 0\end{array}\right\},\hfill \end{array}$$

(12)

where $c_2=\frac{4}{h^2}$. The symmetric part of the curvature tensor for axisymmetric plates is defined as

$$\begin{array}{cc}\hfill {\chi }_{r\theta }& =\frac{1}{2}\left(\frac{\partial {\omega }_{\theta }}{\partial r}-\frac{{\omega }_{\theta }}{r}\right)\hfill \\ \hfill {\chi }_{\theta z}& =\frac{1}{2}\frac{\partial {\omega }_{\theta }}{\partial z},\hfill \end{array}$$

(13)

where ${\omega }_{\theta }=\frac{1}{2}\left(\frac{\partial u_r}{\partial z}-\frac{\partial u_z}{\partial r}\right)$. Thus, the ${\chi }_{r\theta }$ and ${\chi }_{\theta z}$ in terms of the displacements in Equation (10) take the form of

$$\begin{array}{cc}\hfill \left\{\begin{array}{c}{\chi }_{r\theta }\\ {\chi }_{\theta z}\end{array}\right\}& =\frac{1}{4}\left\{\begin{array}{c}\frac{\partial {\theta }_r}{\partial r}-\frac{{\partial }^2{w}_0}{\partial r^2}-\frac{1}{r}\left({\theta }_r-\frac{\partial w_0}{\partial r}\right)\\ -2\frac{\partial {\theta }_z}{\partial r}\end{array}\right\}-\frac{z}{2}\left\{\begin{array}{c}\frac{{\partial }^2{\theta }_z}{\partial r^2}-\frac{1}{r}\frac{\partial {\theta }_z}{\partial r}\\ c_2\left({\theta }_r+\frac{\partial \lambda }{\partial r}\right)+\frac{\partial {\varphi }_z}{\partial r}\end{array}\right\}\hfill \\ \hfill & -\frac{z^2}{4}\left\{\begin{array}{c}{c}_2\left[\frac{\partial {\theta }_r}{\partial r}+\frac{{\partial }^2\lambda }{\partial r^2}-\frac{1}{r}\left({\theta }_r+\frac{\partial \lambda }{\partial r}\right)\right]+\frac{{\partial }^2{\varphi }_z}{\partial r^2}-\frac{1}{r}\frac{\partial {\varphi }_z}{\partial r}\\ 0.\end{array}\right\}.\hfill \end{array}$$

(14)

3.2. Governing Equations

In this study, the soft-coupled thermoelastical behavior of functionally-graded porous materials was analyzed using the finite element method. The equations of equilibrium and the weak form finite element model for static bending problems of axisymmetric plate were obtained using the principle of virtual displacement.

$$\begin{array}{cc}\hfill 0=-& {\int }_V\left({\sigma }_{rr}\delta {\epsilon }_{rr}+{\sigma }_{\theta \theta }\delta {\epsilon }_{\theta \theta }+{\sigma }_{zz}\delta {\epsilon }_{zz}+2{\sigma }_{rz}\delta {\epsilon }_{rz}+2m_{r\theta }\delta {\chi }_{r\theta }+2m_{\theta z}\delta {\chi }_{\theta z}\right)dV\hfill \\ \hfill & +{\int }_V\left({\overline{f}}_i\delta u_i+{\overline{c}}_i\delta {\omega }_i\right)dV+{\int }_S\left({\overline{t}}_i\delta u_i+{\overline{s}}_i\delta {\omega }_i\right)dS\hfill \\ \hfill & +{\int }_{{\mathsf{\Omega }}^t}\left(q_i^t\delta u_i+p_i^t\delta {\omega }_i\right)d{\mathsf{\Omega }}^t+{\int }_{{\mathsf{\Omega }}^b}\left(q_i^b\delta u_i+p_i^b\delta {\omega }_i\right)d{\mathsf{\Omega }}^b,\hfill \end{array}$$

(15)

where ${\sigma }_{ij}$ and $m_{ij}$ are the symmetric part of the stress tensor and the deviatoric part of the couple stress tensor. ${\overline{f}}_i$ and ${\overline{c}}_i$ are the body forces and couples. ${\overline{t}}_i$, ${\overline{s}}_i$, and $\overline{d}$ are the surface forces and couples on the side surfaces. $q_i^{\alpha }$ and $p_i^{\alpha }$ are the surface forces and couples on top $\left(\alpha =t\right)$ and bottom $\left(\alpha =b\right)$ surfaces.

The governing equations of quasi-3D Reddy third order theory are

$$\begin{array}{cc}\hfill 0=& \frac{1}{r}\left[\frac{\partial }{\partial r}\left(r{N}_{rr}^{\left(0\right)}\right)-N_{\theta \theta }^{\left(0\right)}\right]+F_r^{\left(0\right)}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill 0=& \frac{1}{r}\left\{\frac{\partial }{\partial r}\left(r{N}_{rr}^{\left(1\right)}\right)-N_{\theta \theta }^{\left(1\right)}-c_1\left(\frac{\partial }{\partial r}\left(r{N}_{rr}^{\left(3\right)}\right)-N_{\theta \theta }^{\left(3\right)}\right)-r\left(N_{rz}^{\left(0\right)}-c_2{N}_{rz}^{\left(2\right)}\right)\right.\hfill \\ \hfill & \left.+\frac{1}{2}\left[\frac{\partial }{\partial r}\left(r{M}_{r\theta }^{\left(0\right)}\right)+M_{r\theta }^{\left(0\right)}-c_2\left(\frac{\partial }{\partial r}\left(r{M}_{r\theta }^{\left(2\right)}\right)+M_{r\theta }^{\left(2\right)}-2r{M}_{\theta z}^{\left(1\right)}\right)\right]\right\}\hfill \\ \hfill & +F_r^{\left(1\right)}-c_1{F}_r^{\left(3\right)}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill 0=& \frac{1}{r}\left\{\frac{\partial }{\partial r}\left[r{N}_{rr}^{\left(0\right)}\left(\frac{\partial w}{\partial r}\right)\right]+c_1\frac{\partial }{\partial r}\left[\frac{\partial }{\partial r}\left(r{N}_{rr}^{\left(3\right)}\right)-r{N}_{\theta \theta }^{\left(3\right)}\right]+\frac{\partial }{\partial r}\left(r{N}_{rz}^{\left(0\right)}\right)-c_2\frac{\partial }{\partial r}\left(r{N}_{rz}^{\left(2\right)}\right)\right.\hfill \\ \hfill & \left.+\frac{1}{2}\left[\frac{\partial }{\partial r}\left(\frac{\partial }{\partial r}\left(r{M}_{r\theta }^{\left(0\right)}\right)+M_{r\theta }^{\left(0\right)}\right)+c_2\frac{\partial }{\partial r}\left(\frac{\partial }{\partial r}\left(r{M}_{r\theta }^{\left(2\right)}\right)+M_{r\theta }^{\left(2\right)}\right)-2c_2\frac{\partial }{\partial r}\left(r{M}_{\theta z}^{\left(1\right)}\right)\right]\right\}\hfill \\ \hfill & +F_z^{\left(0\right)}+\frac{c_1}{r}\frac{\partial }{\partial r}\left(r{F}_r^{\left(3\right)}\right)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill 0=& \frac{1}{r}\left\{\frac{1}{2}\frac{\partial }{\partial r}\left[\frac{\partial }{\partial r}\left(r{N}_{rr}^{\left(2\right)}\right)-N_{\theta \theta }^{\left(2\right)}\right]-N_{zz}^{\left(0\right)}\right.\hfill \\ \hfill & \left.+\left[\frac{\partial }{\partial r}\left(\frac{\partial }{\partial r}\left(r{M}_{r\theta }^{\left(1\right)}\right)+M_{r\theta }^{\left(1\right)}-r{M}_{\theta z}^{\left(0\right)}\right)\right]\right\}+F_z^{\left(1\right)}+\frac{1}{2}\frac{\partial }{\partial r}\left(r{F}_r^{\left(2\right)}\right)\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill 0=& \frac{1}{r}\left\{\frac{1}{3}\frac{\partial }{\partial r}\left[\frac{\partial }{\partial r}\left(r{N}_{rr}^{\left(3\right)}\right)-N_{\theta \theta }^{\left(3\right)}\right]-2r{N}_{zz}^{\left(1\right)}\right.\hfill \\ \hfill & \left.+\frac{\partial }{\partial r}\left[\frac{\partial }{\partial r}\left(r{M}_{r\theta }^{\left(2\right)}\right)+M_{r\theta }^{\left(2\right)}-2r{M}_{\theta z}^{\left(1\right)}\right]\right\}+F_z^{\left(2\right)}+\frac{1}{3}\frac{\partial }{\partial r}\left(r{F}_r^{\left(3\right)}\right),\hfill \end{array}$$

(20)

where $\left\{N_{ij}^{\left(k\right)},M_{ij}^{\left(k\right)}\right\}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{z}^k\left\{{\sigma }_{ij},m_{ij}\right\}dz$ and $F_i^{\left(k\right)}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{z}^k\left\{f_i+\left[q_i^t-1^k{q}_i^b\right]\right\}dz$. Note that the body couple ${\overline{c}}_{\theta }$ is omitted in the governing equation.

The temperature distribution through thickness direction can be determined by solving the steady state energy equation,

$$-\frac{d}{dz}\left(k\left(z\right)\frac{dT}{dz}\right)=0,$$

(21)

where $k\left(z\right)$ is heat conductivity and T is the temperature. The effective thermal conductivity is defined using the Maxwell–Eucken model described by Deng et al. [51]:

$$k\left(z\right)=k_s\left(z\right)\left[{\displaystyle \frac{k_f+2k_s\left(z\right)+2\mathsf{\Phi }\left(k_f-k_s\left(z\right)\right)}{k_f+2k_s\left(z\right)-\mathsf{\Phi }\left(k_f-k_s\left(z\right)\right)}}\right],$$

(22)

where $k_s\left(z\right)$ and $k_f$ are the thermal conductivity of the solid and fluid phases, respectively, and $\mathsf{\Phi }$ is the porosity. In this study, the thermal conductivity of the solid is obtained using a power-law distribution described in previous section.

4. Finite Element Model

A weak from Galerkin finite element model for the circular plate bending is developed using the principle of virtual displacement (15) and a weak form is directly developed from the energy Equation (21) for steady state heat conduction problem. The details of weak form Galerkin finite element model can be found in Reddy [52].

The temperature T and generalized displacements $\left(u_0,{\theta }_r,w_0,{\theta }_z,{\varphi }_z\right)$ are approximated in following form:

$$\begin{array}{cc}\hfill T\left(z\right)=& \sum _{j=1}^n{T}_j{\widehat{\psi }}_j\left(z\right)\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill u_0\left(r\right)=& \sum _{j=1}^n{u}_j{\psi }_j\left(r\right)\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\theta }_r\left(r\right)=& \sum _{j=1}^n{\theta }_j{\psi }_j\left(r\right)\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill w_0\left(r\right)=& \sum _{J=1}^{2n}{\Delta }_J^{\left(1\right)}{\varphi }_J\left(r\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\theta }_z\left(r\right)=& \sum _{J=1}^{2n}{\Delta }_J^{\left(2\right)}{\varphi }_J\left(r\right)\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\varphi }_z\left(r\right)=& \sum _{J=1}^{2n}{\Delta }_J^{\left(3\right)}{\varphi }_J\left(r\right),\hfill \end{array}$$

(28)

where $T_j$ are nodal temperatures through thickness direction; $u_j,{\theta }_j,$ and $w_j$ are nodal displacements in the radial direction; ${\widehat{\psi }}_j$ and ${\psi }_j$ are the Lagrange interpolation functions; ${\varphi }_J$ are the Hermite interpolation functions; ${\Delta }_J^{\left(i\right)}$ are generalized deflections and $i=1,2,3$ correspond to $w_0,{\theta }_z,{\varphi }_z$, respectively; n is the number of nodes in an element.

The finite element model of the steady state heat conduction problem is given by

$${\left[\begin{array}{c}\widehat{K}\end{array}\right]}^e\left\{\begin{array}{c}{T}^e\end{array}\right\}={\left\{\begin{array}{c}\widehat{F}\end{array}\right\}}^e,$$

(29)

where the stiffness matrix and external heat flux are

$$\begin{array}{cc}\hfill {\widehat{K}}_{ij}& ={\int }_{-\frac{h}{2}}^{\frac{h}{2}}k\left(z\right)\left(\frac{d{\widehat{\psi }}_i}{dz}\frac{d{\widehat{\psi }}_i}{dz}\right)dz\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\widehat{F}}_i& =q_i.\hfill \end{array}$$

(31)

The finite element model of an axisymmetric plate static bending is given by

$${\left[\begin{array}{ccccc}{K}^{11}& K^{12}& K^{13}& K^{14}& K^{15}\\ K^{21}& K^{22}& K^{23}& K^{24}& K^{25}\\ K^{31}& K^{32}& K^{33}& K^{34}& K^{35}\\ K^{41}& K^{42}& K^{43}& K^{44}& K^{45}\\ K^{51}& K^{52}& K^{53}& K^{54}& K^{55}\end{array}\right]}^e{\left\{\begin{array}{c}\left\{u_0\right\}\\ \left\{{\theta }_r\right\}\\ \left\{w_0\right\}\\ \left\{{\theta }_z\right\}\\ \left\{{\varphi }_z\right\}\end{array}\right\}}^e={\left\{\begin{array}{c}\left\{F^1\right\}\\ \left\{F^2\right\}\\ \left\{F^3\right\}\\ \left\{F^4\right\}\\ \left\{F^5\right\}\end{array}\right\}}^e.$$

(32)

The elements of the stiffness matrix, $K^{lm}$, and the elements of force vector, $F^l$, are defined in Appendix A.

The solution of the nonlinear finite element model (32) is obtained using Newton’s iteration procedure. The linearized element equations take the form of

$${\mathbf{T}}^e\left({\Delta }^{(i-1)}\right)\delta {\Delta }^{\left(i\right)}=-{\mathbf{R}}^e\left({\Delta }^{(i-1)}\right),$$

(33)

where ${\mathbf{T}}^e$ is the tangent stiffness matrix, $\delta {\Delta }^{\left(i\right)}$ is incremental displacements at the ith iteration, and ${\mathbf{R}}^e$ is the residual vector. The tangent matrix and residual are defined as [52]

$${\mathbf{T}}^e=\frac{\partial {\mathbf{R}}^e}{\partial {\Delta }^e},{\mathbf{R}}^e={\left({\mathbf{K}}^e{\Delta }^e-{\mathbf{F}}^e\right)}^{(i-1)}.$$

(34)

By solving the assembled global system equation, the global incremental displacement vector at ith iteration, $\delta {\mathbf{U}}^{\left(i\right)}$ is obtained.

$$\delta {\mathbf{U}}^{\left(i\right)}=-{\left[\mathbf{T}\left({\mathbf{U}}^{(i-1)}\right)\right]}^{-1}{\mathbf{R}}^{(i-1)}.$$

(35)

The total displacement at the ith iteration is obtained by adding the incremental solution at the ith iteration to the previous solution at the $(i-1)$th iteration [39].

$${\mathbf{U}}^{\left(i\right)}={\mathbf{U}}^{(i-1)}+\delta {\mathbf{U}}^{\left(i\right)}.$$

(36)

In this study, we considered geometrical nonlinearity with elastic material behavior. For this purpose, the Newton’s iteration is sufficient to obtain the converged solutions. However, when limit load, softening branches, or snap-through behavior are considered, another solution procedure, such as the arc length method, should be considered. These solution procedures can be used in conjunction with various numerical methods such as isogeometric techniques [53,54] or the Rayleigh Ritz method [55] in addition to the finite element method.

5. Numerical Results

In the numerical examples, we considered several examples of annular circular plates with various boundary conditions such as simply-supported and clamped boundary conditions. To validate the developed finite element model, we compared our results with available studies in the literature. We also conducted convergence studies to obtain optimal mesh size and different quadrature rules to make sure we avoided any locking phenomena. In this study, we used 16 elements and full quadrature rules for linear parts of the stiffness matrix and reduced quadrature rules for shear, nonlinear, and couple stress parts of the stiffness matrix.

Figure 2 shows the annular plate we studied. The numerical parameters for the validation study were adapted from the study of Reddy et al. [39]: h = 0.1, $r_o=$ 10h, $r_i=0.25r_0$, $E_1={10}^6$, and $E_2={10}^5$.

Figure 3, Figure 4 and Figure 5 show that maximum deflection versus the load parameter $P=\frac{q_0{h}^4}{E_c{r}_o^4}$ at the free edges, where $q_0$ is a distributed load on the top surface, h is the plate thickness, $E_c$ is the Young’s modulus of ceramic materials on the top surface, and $r_o$ is the outer radius of the annular plate. The developed finite element model shows a good agreement with the study of Reddy et al. [39]. In this figure, the effects of the material variations based on the change of the power-law index, and the effect of length scale parameter are presented with various boundary conditions.

With the validated finite element model, we evaluated the effects of various parameters such as the length scale parameters, the shape of porosity distribution, power law index, and boundary conditions. In this study, we considered a porous functionally-graded material with Monel and zirconia and the material properties of them follow [56]:

$$\begin{array}{cc}\hfill K_m& =227.24\mathrm{GPa},{\mu }_m=65.55\mathrm{GPa},{\alpha }_m=15\times {10}^{-6}/\mathrm{K},k_m=25\mathrm{W}/\mathrm{mK}\mathrm{for}\mathrm{Monel}\hfill \\ \hfill K_c& =125.83\mathrm{GPa},{\mu }_c=58.077\mathrm{GPa},{\alpha }_c=10\times {10}^{-6}/\mathrm{K},k_c=2.09\mathrm{W}/\mathrm{mK}\mathrm{for}\mathrm{zirconia},\hfill \end{array}$$

where $K_i$ is the bulk modulus, ${\mu }_i$ is the shear modulus, ${\alpha }_i$ is the thermal expansion coefficient, $k_i$ is the thermal conductivity, and the subscription m and c indicate metal and ceramic, respectively. We assumed that the porous is filled with the air and the thermal conductivity of the air is assumed to be $k_a=0.025572\mathrm{W}/\mathrm{mK}$.

To induce thermal load, two different temperatures were applied on the top and bottom surfaces; 500 K was applied on the top surface and 300 K was applied on the bottom surface. Figure 6 shows the temperature distribution through plate thickness depending on the variation of material constituents and porosity distribution. The temperature distribution was obtained by solving the energy Equation (21). Three different types of porosity distributions and the variation of material constituents were considered. In the area where the volume fraction for porosity is larger, the thermal resistance becomes larger and the temperature change through the thickness is less than the area where the volume fraction of porosity is lesser. With a larger power-law index, the effective thermal conductivity is increased and thermal resistance becomes smaller because the volume fraction of metal is increased.

For illustration purposes, the same plate geometry as Reddy et al. [39] was used, and the effects of various parameters with clamped and simply-supported boundary conditions were considered.

Figure 7 and Figure 8 show effects of the length scale parameter on the maximum deflection of FGM plates. When the power-law index is larger, the volume fraction of Monel is larger and the FGM plate becomes stiffer. The length scale parameter can capture the stiffening behavior in micro scale structures. The deflections of FGM plates with various length scale parameters are shown in Figure 9 and Figure 10.

The overall volume fraction of porosity in three porosity distributions is the same, but the enhanced porous areas are mid, bottom, and top surfaces with Type 1, Type 2, and Type 2 porosity distributions. The porous FGM is softer than non-porous FGM, and Type 1 results in the stiffest plates because the materials on the top and bottom surfaces remain. There are no differences in the plate bending stiffness between Type 1 and Type 2 distribution when a homogeneous material is assumed. When the power-law index is larger than zero, the volume fraction of stiffer material becomes larger in the FGMs. In the case of Type 2 distribution, the volume fraction of the stiffer material is decreased, and in the case of Type 3, the volume fraction of the softer material is decreased. Therefore, Type 2 will be softer than Type 3 in the case of FGMs. Figure 11 and Figure 12 show the effects of three different porosity distributions and material variations on the maximum deflections with clamped and simply-supported outer edges, respectively. The deflections along the radial direction are shown in Figure 13 and Figure 14.

Figure 15, Figure 16, Figure 17 and Figure 18 show a normal stress distribution through the plate thickness. In the case of the clamped outer edge, the normal stress in the area where the volume fraction of porosity is larger is smaller than the area where the volume fraction of porosity is smaller because the area with larger porosity is softer than the other areas. It is clearly shown that the normal stress at the bottom surface $(z=-h/2)$ with porosity distribution Type 3 is larger than porosity distribution Type 2, which enhances the porosity distribution in the lower area of the plates. In the case of the simply-supported outer edge, the normal stress distribution is a parabolic shape unlike the case of the clamped outer edge. This is because the thickness stretch is not constrained in the case of simply-supported boundary conditions. Only the mid plane deflection, $w_0$, is constrained. The nonzero length scale parameters make the FGM plate stiffer, but there are no material property changes, which results in smaller stresses.

Figure 19, Figure 20, Figure 21 and Figure 22 show transverse shear stress distributions through the plate thickness. Similar effects of porosity distribution and length scale parameter on the transverse shear stresses are observed. The length scale parameter makes the FGM plate stiffer and smaller stress values are obtained. However, the porosity distribution affects the material properties and larger stress values are present in stiffer areas. The proposed quasi-3D Reddy third-order plate theory can capture a parabolic variation of the transverse shear stresses, and it does not require a shear correction factor that is present in low order shear deformation theories.

Figure 23, Figure 24, Figure 25 and Figure 26 show the effect of thermal load. The thermal load is induced by temperature boundary conditions; 500 K is applied on top surface and 300 K is applied on bottom surfaces. In the case of the clamped outer edge, the deflection is due to thermal load in the negative direction because the plate bends down due to the thermal load. This is clearly shown in Figure 25. In the case of the simply-supported outer edge, the plate bends down due to the thermal load at the same place; the plate rotates about the outer edge which results in the positive deflection due to the thermal load.

6. Conclusions

In this study, a nonlinear finite element (FE) model for axisymmetric circular/annular plates was presented. The developed finite element model accounts for geometric nonlinearity, variation of material constituents, microstructure size effects, and effects of porosity distributions. Using the developed FE model, the bending behavior of functionally-graded axisymmetric annular plates under thermomechanical loads was analyzed.

Numerical analysis results for an axisymmetric annular plate with various boundary conditions were presented. A parametric study was conducted to understand the effects of porosity distributions, the variation of material properties, and microstructure size on the bending behavior of axisymmetric annular plates. In summary, the following results were observed:

• The presence of pores results in higher thermal resistance and reduces the temperature variation;
• With a larger power law index, the plate becomes stiffer because the stiffer material is placed on the bottom surface;
• The length scale parameter can capture stiffening effects in microstructures. The stiffening effect does not change the material properties, so stress values are decreased with nonzero length scale parameters;
• The thermal and mechanical behavior of FGM plates highly depends on the porosity distribution type. The presence of pores makes the plate softer by reducing the moduli, resulting in smaller stress values;
• Depending on the boundary conditions, thermal loads can result in opposite deflections due to constrained rotational degrees of freedom.

The presented finite element model can be extended to an asymmetric circular/annular plate.