1 Introduction

Engineering applications may require a combination of various material properties. However, not all of these mechanical properties can be found in one material. If the materials are used together, it can cause large stress levels between the layers. One method to prevent such unfavorable consequences is to use functionally graded materials (FGMs). These material properties are gradually changing. By doing this, the structures' properties can be improved and the benefits of having multiple minerals in a single structure can be realized simultaneously. This kind of material has applications in high-temperature structures found in the nuclear, microelectronic and aerospace industries. In the field of processing techniques, successful examples of FGM were given [1]. Static analysis of FGM cylinders was done using a mesh-free method, and good results were obtained [2]. FGMs consist of two or more materials, and it is considered as an inhomogeneous material [3].

It seems that the material distributions and the strength of the FGM plates were significantly correlated. Research indicates that regulating the material's fracture behavior can lead to an increase in plate strength [4]. By applying a layer of FGM to the wire arc additive manufacturing process of laser metal deposition, the impact of bimetallic structures can be mitigated [5]. The analysis of FGM circular and rectangular plates were done using the FEM method [6]. Applications for composite materials have expanded significantly in recent years[7] because it is among the factors that lessen delamination and stress. The porous material appeared to have a positive effect on the design of sandwich plates, it both reduces the overall mass and increases the flexural strength [8]. It is possible to predict the stresses and deflections of 2D FG porous beams using the FEM approach [9]. On the other hand, the FGP material provides the structures with improved energy absorption, vibration reduction and thermal insulation capabilities [10]. Using the FEM approach, numerical static analysis of FGM plates were completed effectively [11,12,13,14]. Conversely, several studies use the FEM method to address the static response and buckling analysis of structural elements [15,16,17,18]. Additionally, the examination of FGM shells is one of the most prevalent subjects. It holds significance because of its extensive use in mechanical structures. When the numerical simulation results were compared to the literature, they demonstrated remarkable precision [19,20,21]. Several formulas and numerical solutions have been proposed to analyze heated FGM annular plates [22]. Conversely, research is being done on free vibration analysis for FGP plates [23]. Additionally, the researchers concentrated on examining the dynamic reactions and vibrations of FGP beams and plates [24,25,26,27,28,29,30]. The dynamic responses of axisymmetric circular plates have been successfully investigated [31]. The bending response of elliptic plates were examined by [32]. Numerous studies examine the stability analysis of FGM plates when they are loaded electromechanically [33]. Researchers [34] have been interested in the dynamic analysis of annular plates and solid circular and they examined effect of thickness functions and material variation exponents. FGM's nondeterministic structure verifies its responses to static loads in situations where the material's properties are uncertain [35]. A FEM-based ABAQUS tool is used to analyze the castellated beams' static response with different web opening sizes [36]. A parametric study is carried out to demonstrate how the natural frequencies and mode shapes of the structure are affected by various geometric parameters, boundary conditions, weight fraction of graphene platelets, porosity coefficient, distribution of porosity and graphene platelets' dispersion pattern. Furthermore, the natural frequency analyzed of a functionally graded porous joined truncated conical–cylindrical shell reinforced by graphene platelets [37]. Functionally graded (FG) porous rectangular plates' tendency to buckle under various loading scenarios was investigated in an alternate study [38]. Stress wave propagation and free vibration response of functionally graded graphene platelets of the shell are investigated under three different porosity distributions [39]. The performance of an axisymmetric rotating truncated cone subjected to a thermal loading composed of functionally graded (FG) porous materials is examined for (uniform and functionally graded), porosity distribution types [40]. The static response of functionally graded saturated porous rotating truncated cones was examined using three distinct patterns for the porosity distribution [41]. The low-velocity impact study results of the graphene platelet-reinforced circular plate (FG) demonstrate that the porosity distribution affects the structure's instantaneous central deflection and time history of contact force more than other parameters [42]. Static and transient responses of the plate based on Hamilton’s principle for shear deformation plate theory using FEM was analyzed for different porosity coefficients [43]. Apart from the porosity coefficient, the natural frequency of FG is significantly influenced by the porosity distribution [44]. It has been found that a high weight fraction and a high porosity coefficient generate the strongest and lowest effects of FG on the buckling loads of the structure [45]. Three functions for porosity distribution were investigated in order to determine the optimal distribution of porosity for the torsional buckling forces of the shell [46].

The primary purpose of this paper is to suggest an efficient numerical method to examine the axisymmetric bending response of FGP circular plates. The FEM method is used to analyze the model for the numerical solution. The results obtained are compared with those of the available in the literature.

Literature review shows that there are several studies dealing with analyzing the static response of circular plates. Accordingly, the static response of uniform (UMCR), symmetric (SMCR) and monolithic (MMCR) functionally graded porous (FGP) circular plates is studied via the finite element method (FEM) by using eight-node quadratic quadrilateral axisymmetric element in this paper for the first time. The modulus of elasticity of the plate varies continually in the thickness direction. While the Poisson ratio (v) value is assumed to be constant. For verification, the static response of the FGM plate is examined for porosity coefficient (e0 = 0.2, 0.4 and 0.9). Moreover, several examples of (h/r) are analyzed using different values of radius (r = 1–5 m), for various types of porosity coefficients (0.1–0.9). The clamped boundary condition is used for verification, while the clamped and roller-support types of boundary conditions are employed in the analysis of FGP plates. As a consequence, the FGP circular plate's deflection results show excellent agreement when compared to those obtained with the ANSYS APDL software. However, there is a great deal of agreement when comparing the FGM circular plate results with those reported in the literature that is currently in publication, including solutions derived using the complementary functions method (CFM).

2 Materials and Method

In this study, circular plates are subjected to the axisymmetric loading. Three-dimensional (3D) problems can be analyzed as simple two-dimensional (2D) problems. The problems are reviewed in the r-z plane. Three different material distributions are assumed in this study as shown in Fig. 1. The materials' mechanical properties are supposed to change along the thickness direction. Uniform porosity type (UMCR) is regularly distributed throughout the plate, whereas the effects of the SMCR type are assumed to be harmonic around the center of the plate. Regarding the monolithic MMCR type, the effects of porosity vary from highest to lowest ratio, along the thickness direction.

Fig. 1
figure 1

Geometry of Circular plate with different porosity distributions

Full size image

2.1 Stress and Strains Components

Considering the unit volume, the potential energy can be obtained as in Eqs. (1–2), which is dealing with the following parameters [47];

dA = area, dz = thickness direction, dr = radial direction, = \({\text{d}}\theta\)angular direction

$$\Pi =\frac{1}{2}\int_{0}^{2\uppi }{{\int_{\text{A}} }}{{\varvec{\upsigma}}}^{{\text{T}}}{\varvec{\upvarepsilon}} r dA {\text{d}}\theta -{\int }_{0}^{2\uppi }{{\int_{\text{A}} }}{\mathbf{m}}^{{\text{T}}}\mathbf{f} r {\text{d}}A d\theta -{\int }_{0}^{2\uppi }{{\int_{\text{L}} }}{\mathbf{m}}^{{\text{T}}}\mathbf{T} r {\text{d}}l d\theta -\sum_{{\text{i}}}{\mathbf{m}}_{{\text{i}}}^{{\text{T}}}{\mathbf{P}}_{\mathrm{i }}$$
(1)
$$\Pi =2\uppi \left(\frac{1}{2}{{\int_{\text{A}} }}{{\varvec{\upsigma}}}^{{\text{T}}}{\varvec{\upvarepsilon}} r dA -{{\int_{\text{A}} }}{\mathbf{m}}^{{\text{T}}}\mathbf{f} r {\text{d}}A-{{\int_{\text{L}} }}{\mathbf{m}}^{{\text{T}}}\mathbf{T} r {\text{d}}l \right)-\sum_{{\text{i}}}{\mathbf{m}}_{{\text{i}}}^{{\text{T}}}{\mathbf{P}}_{\mathrm{i }}$$
(2)

The vertical (m) and horizontal (n) displacements of any point can be written as:

$$m\, = \,m \, \left( {r,z} \right), \, n\, = \,n \, \left( {r,z} \right)$$
(3)

The weight and distributed loads are given as:

$$\{{\text{f}}\} = \{{{f}_{r} , {f}_{z}\}}^{T }, \{t\} = \{{{t}_{r} , {t}_{z}\}}^{T}$$

The relationship between the strain and displacement is given as follows: (4)

$$\left\{\varepsilon \right\}={\left\{{\varepsilon }_{r}, {\varepsilon }_{z}, {\varepsilon }_{\theta }, {\gamma }_{{\text{rz}}}\right\}}^{{\text{T}}}=\left\{\begin{array}{c}\frac{\partial {\text{m}}}{\partial {\text{r}}} , \frac{\partial {\text{n}}}{\partial {\text{z}}} , \frac{m}{{\text{r}}}\\ \end{array} , \frac{\partial {\text{m}}}{\partial {\text{z}}}+ \frac{\partial {\text{n}}}{\partial {\text{r}}}\right\}$$
(5)

The stress vector is as below:

$$\left\{ \sigma \right\}\, = \,\left\{ {\sigma_{r} ,\sigma_{Z} ,\sigma_{\theta } , \tau_{{{\text{rz}}}} } \right\}^{T}$$
(6)

The constitutive relation is shown in Eq. (7).

$$\left\{\sigma \right\}=\left[C\right] \left\{\varepsilon \right\}$$
(7)
$$\begin{aligned}&\left\{\begin{array}{c}\begin{array}{c}{\sigma }_{r}\\ {\sigma }_{z}\end{array}\\ \begin{array}{c}{\sigma }_{\theta }\\ {\tau }_{{\text{rz}}}\end{array}\end{array}\right\}\\ &=\left[\begin{array}{cc}\begin{array}{cc}(1-\nu )p(z)& {\nu p}(z)\\ & (1-\upnu )p(z)\end{array}& \begin{array}{cc}{\nu p}(z) & 0\\ {\nu p}(z) & 0\end{array}\\ \begin{array}{cc} & \\ {\text{Symmetrical}}& \end{array}& \begin{array}{cc}\left(1-\upnu \right){\text{p}}(z)& 0\\ & G(z)\end{array}\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}{\varepsilon }_{r}\\ {\varepsilon }_{z}\end{array}\\ \begin{array}{c}{\varepsilon }_{\theta }\\ {\gamma }_{{\text{rz}}}\end{array}\end{array}\right\}\end{aligned}$$
(8)

The Young’s modulus E(z) and the shear modulus G(z) are defined by Eqs. (9–10).

$$p(z)=\frac{E\left(z\right)}{\left(1+v\right)\left(1-2v\right)}$$
(9)
$$G(z)=\frac{E(z)}{2(1+\nu )}$$
(10)

The porosity for SMCR, UMCR and MMCR models can be described by the relations in Eqs. (11–17) [48].

$$\begin{aligned}E\left(z\right)&={E}_{1}^{ }\left[1-{\text{cos}}\left(\frac{\pi z}{h}\right)\right]{+{E}_{0}^{ }{\text{cos}}\left(\frac{\pi z}{h}\right)}\\& =E_{1}^{ }\left[1{-e}_{0}^{ }{\text{cos}}\left(\frac{\pi z}{h}\right)\right]\end{aligned}$$
(11)
$$\begin{aligned}G\left(z\right)&={G}_{1}^{ }\left[1-{\text{cos}}\left(\frac{\pi z}{h}\right)\right]{+{G}_{0}^{ }{\text{cos}}\left(\frac{\pi z}{h}\right)}\\&=G_{1}^{ }\left[1{-e}_{0}^{ }{\text{cos}}\left(\frac{\pi z}{h}\right)\right]\end{aligned}$$
(12)
$$\begin{aligned}E\left(z\right)&={E}_{1}^{ }\left[1{-e}_{0}^{ }\left(\frac{1}{{e}_{o}}-\frac{1}{{e}_{o}} {\left(\frac{2}{\pi } \sqrt{1-{e}_{o}}-\frac{2}{\pi }+1\right)}^{2}\right)\right]\\ &={E}_{1}^{ }\left[1{-e}_{0}^{ }\psi \right]\end{aligned}$$
(13)
$$\begin{aligned}{G\left(z\right)}& =G_{1}^{ }\left[1{-e}_{0}^{ }\left(\frac{1}{{e}_{o}}-\frac{1}{{e}_{o}} {\left(\frac{2}{\pi } \sqrt{1-{e}_{o}}-\frac{2}{\pi }+1\right)}^{2}\right)\right]\\ &=G_{1}^{ }\left[1{-e}_{0}^{ }\psi \right]\end{aligned}$$
(14)
$$\psi = \frac{1}{{e}_{o}}-\frac{1}{{e}_{o}} {\left(\frac{2}{\pi } \sqrt{1-{e}_{o}}-\frac{2}{\pi }+1\right)}^{2}$$
(15)
$$E\left(z\right)={E}_{1}\left[1-{e}_{0}{\text{cos}}\left(\frac{\pi z}{2h}+\frac{0.5\pi }{2}\right)\right]$$
(16)
$$G\left(z\right)={G}_{1}^{ }\left[1-{e}_{0}{\text{cos}}\left(\frac{\pi z}{2h}+\frac{0.5\pi }{2}\right)\right]$$
(17)

In these equations, \({e}_{0}\) is the porosity coefficient and  h shows the thickness of the plate.

2.2 FEM Approach

As shown in Fig. 2, an isoparametric eight-node quadratic element is used in this study.

Fig. 2
figure 2

Eight-node quadratic quadrilateral element

Full size image

The coordinates of nodes can be obtained with the aid of Eq. (18), where, ri and zi are the coordinates of the node i,

$$r= \sum_{{\text{i}}=1 }^{8}{N}_{{\text{i}}}{r}_{{\text{i}}} ; z= \sum_{{\text{i}}=1 }^{8}{N}_{{\text{i}}}{z}_{{\text{i}}}$$
(18)

Ni are the quadratic shape functions. Quadratic base polynomials with 8 parameters are chosen for displacements as:

$$m \left(\xi ,\eta \right)={\left\{P\left(\xi ,\eta \right)\right\}}^{T}\left\{a\right\} , n \left(\xi ,\eta \right)={\left\{P\left(\xi ,\eta \right)\right\}}^{T}$$
(19)
$$m= {a}_{1}+\xi {a}_{2}+\eta {a}_{3}+{\xi }^{2}{a}_{4}+\xi \eta {a}_{5}+{\eta }^{2}{a}_{6}+{\xi }^{2}\eta {a}_{7}+\xi {\eta }^{2}{a}_{8}$$
(20)
$$n= {a}_{1}+\xi {a}_{2}+\eta {a}_{3}+{\xi }^{2}{a}_{4}+\xi \eta {a}_{5}+{\eta }^{2}{a}_{6}+{\xi }^{2}\eta {a}_{7}+\xi {\eta }^{2}{a}_{8}$$
(21)
$$m\left({\xi }_{i}, {\eta }_{i}\right)= {m}_{i} (i\hspace{0.17em}=\hspace{0.17em}\mathrm{1,2}, \dots ,8)$$
(22)

The mi displacement for 8 nodes of the reference element is given in matrix form as follows:

$$\underbrace {{\left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} {m_{1} } \\ {m_{2} } \\ \end{array} } \\ {\begin{array}{*{20}c} {m_{3} } \\ {m_{4} } \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {m_{5} } \\ {m_{6} } \\ \end{array} } \\ {\begin{array}{*{20}c} {m_{7} } \\ {m_{8} } \\ \end{array} } \\ \end{array} } \\ \end{array} } \right\}}}_{{M_{d} }} = \underbrace {{\left[ {\begin{array}{*{20}c} 1 & {\xi_{1} } & {\eta_{1} } & {\xi_{1}^{2} } & {\xi_{1} \eta_{1} } & {\eta_{1}^{2} } & {\xi_{1}^{2} \eta_{1} } & {\xi_{1} \eta_{1}^{2} } \\ 1 & {\xi_{2} } & {\eta_{2} } & {\xi_{2}^{2} } & {\xi_{2} \eta_{2} } & {\eta_{2}^{2} } & {\xi_{2}^{2} \eta_{2} } & {\xi_{2} \eta_{2}^{2} } \\ 1 & {\xi_{3} } & {\eta_{3} } & {\xi_{3}^{2} } & {\xi_{3} \eta_{3} } & {\eta_{3}^{2} } & {\xi_{3}^{2} \eta_{3} } & {\xi_{3} \eta_{3}^{2} } \\ 1 & {\xi_{4} } & {\eta_{4} } & {\xi_{4}^{2} } & {\xi_{4} \eta_{4} } & {\eta_{4}^{2} } & {\xi_{4}^{2} \eta_{4} } & {\xi_{4} \eta_{4}^{2} } \\ 1 & {\xi_{5} } & {\eta_{5} } & {\xi_{5}^{2} } & {\xi_{5} \eta_{5} } & {\eta_{5}^{2} } & {\xi_{5}^{2} \eta_{5} } & {\xi_{5} \eta_{5}^{2} } \\ 1 & {\xi_{6} } & {\eta_{6} } & {\xi_{6}^{2} } & {\xi_{6} \eta_{6} } & {\eta_{6}^{2} } & {\xi_{6}^{2} \eta_{6} } & {\xi_{6} \eta_{6}^{2} } \\ 1 & {\xi_{7} } & {\eta_{7} } & {\xi_{7}^{2} } & {\xi_{7} \eta_{7} } & {\eta_{7}^{2} } & {\xi_{7}^{2} \eta_{7} } & {\xi_{7} \eta_{7}^{2} } \\ 1 & {\xi_{8} } & {\eta_{8} } & {\xi_{8}^{2} } & {\xi_{8} \eta_{8} } & {\eta_{8}^{2} } & {\xi_{8}^{2} \eta_{8} } & {\xi_{8} \eta_{8}^{2} } \\ \end{array} } \right]}}_{{P_{{}} }}\underbrace {{\left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} {a_{1} } \\ {a_{2} } \\ \end{array} } \\ {\begin{array}{*{20}c} {a_{3} } \\ {a_{4} } \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {a_{5} } \\ {a_{6} } \\ \end{array} } \\ {\begin{array}{*{20}c} {a_{7} } \\ {a_{8} } \\ \end{array} } \\ \end{array} } \\ \end{array} } \right\}}}_{a}$$
(23)
$$\left\{{M}_{d}\right\}=\left[{P}_{d}\right]\left\{a\right\}$$
(24)
$${\left[{P}_{d}\right]}_{ij}= {{P}_{j}\left({\xi }_{i},{\eta }_{i}\right)} , \left\{{M}_{d}\right\}={\left[{P}_{d}\right]}\left\{a\right\} , \left\{a\right\}={{\left[{P}_{d}\right]}}^{-1}\left\{{M}_{d}\right\}$$
(25)
$$m\left( {\xi ,\eta } \right) = \underbrace {{\left\{ {P\left( {\xi ,\eta } \right)} \right\}^{T} \left[ {P_{d} } \right]_{ }^{ - 1} }}_{{\text{shape functions}}} \left\{ {M_{d} } \right\}$$
(26)

The shape functions are presented as:

$${\left\{{\text{N}}\left(\upxi ,\upeta \right)\right\}}^{{\text{T}}}={\left\{{\text{P}}\left(\upxi ,\upeta \right)\right\}}^{{\text{T}}}{\left[{{\text{P}}}_{{\text{d}}}\right]}^{-1}$$
(27)
$$\begin{aligned}&{\left\{{{\text{N}}}_{1}, {{\text{N}}}_{2}, {{\text{N}}}_{3},{{\text{N}}}_{4},{{\text{N}}}_{5},{{\text{N}}}_{6},{{\text{N}}}_{7},{{\text{N}}}_{8}\right\}}^{{\text{T}}}\\ & ={\left\{1,\upxi ,\upeta ,{\upxi }^{2},\mathrm{\xi \eta },{\upeta }^{2},{\upxi }^{2}\upeta ,\upxi {\upeta }^{2}\right\}}^{{\text{T}}} {\left[{{\text{P}}}_{{\text{d}}}\right]}^{-1}\end{aligned}$$
(28)

The shape functions obtained by multiplying the matrices from Eq. (28) are given in Eqs. (29–36) [49]:

$${{\text{N}}}_{1}=\frac{1}{4}\left(-1+\xi \right) \left(1-\eta \right)\left(1+\xi +\eta \right)$$
(29)
$${N}_{2}=\frac{1}{2}\left(-1+{\xi }^{2}\right) \left(-1+\eta \right)$$
(30)
$${N}_{3}=\frac{1}{4}\left(1+\xi \right) \left(-1+\eta \right)\left(1-\xi +\eta \right)$$
(31)
$${N}_{4}=\frac{1}{2}\left(1+{\xi }^{ }\right) \left(1-{\eta }^{2}\right)$$
(32)
$${N}_{5}=\frac{1}{4}\left(1+\xi \right) \left(1+\eta \right)\left(-1+\xi +\eta \right)$$
(33)
$${N}_{6}=\frac{1}{2}\left(1-{\xi }^{2}\right) \left(1+\eta \right)$$
(34)
$${N}_{7}=\frac{1}{4}\left(-1+\xi \right) \left(1+\eta \right)\left(1+\xi -\eta \right)$$
(35)
$${N}_{8}=\frac{1}{2}\left(-1+{\xi }^{ }\right) \left(-1+{\eta }^{2}\right)$$
(36)

The vertical and horizontal displacements at any point of the element also depend on the shape functions.

$$m= \sum_{{\text{i}}=1 }^{8}{N}_{{\text{i}}}{m}_{{\text{i}}} ;n= \sum_{{\text{i}}=1 }^{8}{N}_{{\text{i}}}{n}_{{\text{i}}}$$
(37)

2.3 Element Stiffness Matrix

Strain and displacements relations are given [47]:

$$\left\{\begin{array}{c}\begin{array}{c}{\varepsilon }_{r}\\ {\varepsilon }_{z}\end{array}\\ \begin{array}{c}{\varepsilon }_{\theta }\\ {\gamma }_{{\text{rz}}}\end{array}\end{array}\right\}=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}\frac{\partial }{\partial {\text{r}}}\\ \begin{array}{c} \\ 0\\ \end{array}\end{array}\\ \begin{array}{c}\frac{1}{{\text{r}}}\\ \frac{\partial }{\partial {\text{z}}}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial }{\partial {\text{z}}}\end{array}\\ \begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial }{\partial {\text{r}}}\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}m\\ n\end{array}\right\}$$
(38)
$$\left\{\varepsilon \right\}=\sum_{{\text{i}}=1 }^{8}\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}\frac{\partial {N}_{{\text{i}}}}{\partial r}\\ \begin{array}{c} \\ 0\\ \end{array}\end{array}\\ \begin{array}{c}\frac{{N}_{{\text{i}}}}{r}\\ \frac{\partial {N}_{{\text{i}}}}{\partial {\text{z}}}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial {N}_{{\text{i}}}}{\partial z}\end{array}\\ \begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial {N}_{{\text{i}}}}{\partial {\text{r}}}\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}{m}_{{\text{i}}}\\ {n}_{{\text{i}}}\end{array}\right\}=\sum_{{\text{i}}=1 }^{8}\left[{B}_{{\text{i}}}\right]\left\{{{\text{d}}}_{{\text{i}}}\right\}$$
(39)
$$\left[{B}_{{\text{i}}}\right]=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}\frac{\partial {N}_{{\text{i}}}}{\partial {\text{r}}}\\ \begin{array}{c} \\ 0\\ \end{array}\end{array}\\ \begin{array}{c}\frac{{N}_{{\text{i}}}}{{\text{r}}}\\ \frac{\partial {N}_{{\text{i}}}}{\partial {\text{z}}}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial {N}_{{\text{i}}}}{\partial z}\end{array}\\ \begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial {N}_{{\text{i}}}}{\partial {\text{r}}}\end{array}\end{array}\end{array}\right] ; \left\{{{\text{d}}}_{{\text{i}}}\right\}= \left\{\begin{array}{c}{m}_{{\text{i}}}\\ {n}_{{\text{i}}}\end{array}\right\}$$
(40)

For all element

$$\left[{B}_{{\text{i}}}\right]=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}\frac{\partial }{\partial {\text{r}}}\\ \begin{array}{c} \\ 0\\ \end{array}\end{array}\\ \begin{array}{c}\frac{1}{{\text{r}}}\\ \frac{\partial }{\partial {\text{z}}}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial }{\partial {\text{z}}}\end{array}\\ \begin{array}{c}\begin{array}{c} \\ 0\\ \end{array}\\ \frac{\partial }{\partial {\text{r}}}\end{array}\end{array}\end{array}\right]\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}{N}_{1}\\ 0\end{array}& \begin{array}{c}0\\ {N}_{1}\end{array}\end{array}& \begin{array}{cc}\begin{array}{ccc}\begin{array}{c}{N}_{2}\\ 0\end{array}& \begin{array}{c}0\\ {N}_{2}\end{array}& \begin{array}{c}\cdots \\ \cdots \end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{N}_{8}\\ 0\end{array}& \begin{array}{c}0\\ {N}_{8}\end{array}\end{array}\end{array}\end{array}\right]$$
(41)
$$\begin{aligned} \left\{ {d_{e} } \right\}^{t} & = \{ m_{1} ,n_{1} ,m_{2} ,n_{2} ,m_{3} ,n_{3} ,m_{4}, \\ & n_{4} ,m_{5} ,n_{5} ,m_{6} ,n_{6} ,m_{7} ,n_{7} ,m_{8} ,n_{8} \}\end{aligned}$$
(42)
$$\frac{{ \partial N_{i} }}{\partial \xi } = \frac{{\partial N_{i} }}{\partial r}\frac{\partial r}{{\partial \xi }} + \frac{{\partial N_{i} }}{\partial z}\frac{\partial z}{{\partial \xi }}$$
(43)
$$\frac{{\partial N_{i} }}{\partial \eta } = \frac{{\partial N_{i} }}{\partial r}\frac{\partial r}{{\partial \eta }} + \frac{{\partial N_{i} }}{\partial z}\frac{\partial z}{{\partial \eta }}$$
(44)
$$\left[ {\begin{array}{*{20}c} {\frac{{\partial N_{i} }}{\partial \xi }} \\ {\frac{{\partial N_{i} }}{\partial \eta }} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\frac{\partial r}{{\partial \xi }}} & {\frac{\partial z}{{\partial \xi }}} \\ {\frac{\partial r}{{\partial \eta }}} & {\frac{\partial z}{{\partial \eta }}} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\frac{{\partial N_{i} }}{\partial r}} \\ {\frac{{\partial N_{i} }}{\partial z}} \\ \end{array} } \right] = \left[ J \right]\left[ {\begin{array}{*{20}c} {\frac{{\partial N_{i} }}{\partial r}} \\ {\frac{{\partial N_{i} }}{\partial z}} \\ \end{array} } \right]$$
(45)

The Jacobian transformation is given as following Eq.:

$$\left[ J \right] = \left[ {\begin{array}{*{20}c} {\frac{\partial r}{{\partial \xi }}} & {\frac{\partial z}{{\partial \xi }}} \\ {\frac{\partial r}{{\partial \eta }}} & {\frac{\partial z}{{\partial \eta }}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\mathop \sum \limits_{i = 1}^{8} \frac{{\partial N_{i} }}{\partial \xi }r_{i} } & {\mathop \sum \limits_{i = 1}^{8} \frac{{\partial N_{i} }}{\partial \xi }z_{i} } \\ {\mathop \sum \limits_{i = 1}^{8} \frac{{\partial N_{i} }}{\partial \eta }r_{i} } & {\mathop \sum \limits_{i = 1}^{8} \frac{{\partial N_{i} }}{\partial \eta }z_{i} } \\ \end{array} } \right]$$
(46)
$$\left[ {\begin{array}{*{20}c} {\frac{{\partial N_{i} }}{\partial r}} \\ {\frac{{\partial N_{i} }}{\partial z}} \\ \end{array} } \right] = \left[ J \right]^{ - 1} \left[ {\begin{array}{*{20}c} {\frac{{\partial N_{i} }}{\partial \xi }} \\ {\frac{{\partial N_{i} }}{\partial \eta }} \\ \end{array} } \right] = \underbrace {{\frac{1}{\det \left[ J \right]}\left[ {\begin{array}{*{20}c} {\frac{\partial z}{{\partial \eta }}} & { - \frac{\partial z}{{\partial \xi }}} \\ { - \frac{\partial r}{{\partial \eta }}} & {\frac{\partial r}{{\partial \xi }}} \\ \end{array} } \right]}}_{{\left[ J \right]^{ - 1} = \left[ {J^{*} } \right]^{ } }}\left[ {\begin{array}{*{20}c} {\frac{{\partial N_{i} }}{\partial \xi }} \\ {\frac{{\partial N_{i} }}{\partial \eta }} \\ \end{array} } \right]$$
(47)
$$\frac{{\partial N_{i} }}{\partial r} = J_{11}^{*} \frac{{\partial N_{i} }}{\partial \xi } + J_{12}^{*} \frac{{\partial N_{i} }}{\partial \eta }$$
(48)
$$\frac{{\partial N_{i} }}{\partial z} = J_{21}^{*} \frac{{\partial N_{i} }}{\partial \xi } + J_{22}^{*} \frac{{\partial N_{i} }}{\partial \eta }$$
(49)

The strain matrix [B] is calculated based on the reference element coordinates as follows:

$$\left[ {B_{i} } \right] = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} {J_{11}^{*} } \\ 0 \\ \end{array} } \\ {\begin{array}{*{20}c} 0 \\ {J_{21}^{*} } \\ \end{array} } \\ \end{array} } & {\begin{array}{*{20}c} {\begin{array}{*{20}c} {J_{12}^{*} } \\ 0 \\ \end{array} } \\ {\begin{array}{*{20}c} 0 \\ {J_{22}^{*} } \\ \end{array} } \\ \end{array} } \\ \end{array} } & {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 \\ {J_{21}^{*} } \\ \end{array} } \\ {\begin{array}{*{20}c} 0 \\ {J_{11}^{*} } \\ \end{array} } \\ \end{array} } & {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 \\ {J_{22}^{*} } \\ \end{array} } \\ {\begin{array}{*{20}c} 0 \\ {J_{12}^{*} } \\ \end{array} } \\ \end{array} } & {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 \\ 0 \\ \end{array} } \\ {\begin{array}{*{20}c} {1/r} \\ 0 \\ \end{array} } \\ \end{array} } \\ \end{array} } \\ \end{array} } \\ \end{array} } \right] \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\frac{{\partial N_{i} }}{\partial \xi }} & 0 \\ {\frac{{\partial N_{i} }}{\partial \eta }} & 0 \\ \end{array} } \\ {\begin{array}{*{20}c} 0 & {\frac{{\partial N_{i} }}{\partial \xi }} \\ 0 & {\frac{{\partial N_{i} }}{\partial \eta }} \\ \end{array} } \\ {\begin{array}{*{20}c} {N_{i} } & 0 \\ \end{array} } \\ \end{array} } \right]$$
(50)

Element stiffness matrix is given as following:

$$\left[ {k_{e} } \right] = \iint\limits_{{A_{e} }} {\left\{ {\mathop \smallint \limits_{0}^{2\pi } \left[ B \right]^{T} \left[ C \right]\left[ B \right]rd\theta } \right\}drdz}$$
(51)
$$dA = drdz = \det \left[ J \right]d\xi d\eta$$
(52)

If the element stiffness matrix is integrated over θ (0-2π), it becomes as follow:

$$\left[ {k_{e} } \right] = 2\pi \mathop \smallint \limits_{ - 1}^{ + 1} \mathop \smallint \limits_{ - 1}^{ + 1} \left[ B \right]^{T} \left[ C \right]\left[ B \right]r\det \left[ J \right]d\xi d\eta$$
(53)
$$\left\{ {f_{h} } \right\} = \iint\limits_{{A_{e} }} {\left\{ {\mathop \smallint \limits_{0}^{2\pi } \left\{ N \right\}^{T} \left\{ f \right\}rd\theta } \right\}drdz}$$
(54)
$$\left\{ {f_{e} } \right\} = \rho g 2\pi \mathop \smallint \limits_{ - 1}^{ + 1} \mathop \smallint \limits_{ - 1}^{ + 1} \left\{ N \right\}^{T} \left\{ f \right\}r\det \left[ J \right]d\xi d\eta$$
(55)
$$\left\{ {f_{e} } \right\} = \mathop \smallint \limits_{{L_{e} }}^{ } \left\{ {\mathop \smallint \limits_{0}^{2\pi } \left\{ N \right\}^{T} \left\{ T \right\}rd\theta } \right\}dl$$
(56)
$$\left\{ {f_{e} } \right\} = \mathop \smallint \limits_{{L_{e} }}^{ } \left\{ {\mathop \smallint \limits_{0}^{2\pi } \left\{ N \right\}^{T} \left\{ T \right\}rd\theta } \right\}{\text{d}}l$$
(57)
$$\left\{ {f_{e} } \right\} = 2\pi \mathop \smallint \limits_{ - 1}^{ + 1} \left\{ N \right\}^{T} \left\{ T \right\} r \frac{{L_{e} }}{2} d\xi$$
(58)

The node load vector for the volume forces and the distributed forces acting on the element boundary is also obtained as follows:

$$\begin{aligned} \left\{ {f_{e} } \right\} &= \rho \left( r \right) g 2\pi \mathop \smallint \limits_{ - 1}^{ + 1} \mathop \smallint \limits_{ - 1}^{ + 1} \left\{ N \right\}^{T} \left\{ f \right\}r\det \left[ J \right]d\xi d\eta \\ & + 2\pi \mathop \smallint \limits_{ - 1}^{ + 1} \left\{ N \right\}^{T} \left\{ T \right\} r \frac{{L_{e} }}{2} d\xi\end{aligned}$$
(59)

Two types of boundary conditions are used in this study and are presented in Table 1.

Table 1 Boundary condition
Full size table

3 Numerical Results and Discussion

First of all, the numerical results of FGM circular plates for two types of boundary conditions and three different material distribution constants (0, 2 and 6), using the FEM are presented in this paper. For this purpose, computer program is coded in Fortran. To validate the chosen technique, the maximum values of the non-dimensional vertical displacements are compared with the available literature and those of the ANSYS (PLANE183) in Table 2. The modulus of elasticity of the plate varies continuously in thickness directions, while the Poisson ratio (v) is supposed to be constant (v = 0.288). Considering the plate is made of aluminum/zirconia. Er depicts the ratio of the modulus of elasticity of aluminum and zirconia as in [50,51,52], which is equal to (Er = 0.396). The Em describes the modulus of elasticity of aluminum which is equal to (Em = 70 Gpa), while the Ec describes the modulus of elasticity of zirconia. A circular plate subjected to a uniformly distributed load is considered. Radius of this plate, r0 = 1.00 m, thickness, t = 0.10 m and load, pz = 10 N/m2. Equation 62 shows the maximum non-dimensional deflection (\(\overline{w }\)) and Dc refers to the flexural rigidity.

Table 2 Verification of the FGM plate results for clamped type (mm)
Full size table

The result obtained for clamped and roller-supported boundary conditions is applicable to this study.

$${E(z)=\left({E}_{m}-{E}_{c}\right){\left(\frac{h-2z}{2h}\right)}^{n}+{E}_{c}}$$
(60)
$${E}_r= {E}_m/{E}_c$$
(61)
$${D}_{c}=\frac{{E}_{c} {h}^{3}}{12(1-{v}^{2})} ; \overline{w }= \frac{64 Dc}{{P}_{Z} {{r}_{0}}^{4 } {w}_{max}}$$
(62)

Reddy et al. [51] used the analytics method, which depended on the first-order shear deformation (Mindlin plate) theory. Saidi et al. [50] used the closed-form method to employ the solutions based on the unconstrained third-order shear deformation plate theory (UTST), While Noori and Temel [52] used the CFM. It can be seen from Table 2 that the results of the presented approach are in excellent agreement with those of the previous literature. The main purpose of this paper is to investigate the porosity effects and the material properties, which is accepted to be changed in the thickness direction. For the design of the models, the uniform (UMCR) models, symmetric (SMCR)) and monolithic (MMCR) porosity material types are used. The plate is divided into 80 layers in the thickness direction. Taking into consideration, the thickness value of (h) is constant, which is equal to (h = 0.1 m), while the radius value is assumed to vary which is between (1–5 m). Poisson ratio (v) is taking constant with a value equal to (v = 0.3). Moreover, the plate is divided into 30 layers in the radial direction. As a result, the model of the plate consists of 2400 elements. To verify the results, the maximum vertical displacement values obtained from this study for porosity coefficient (e0 = 0.2, 0.4 and 0.9) compared with ANSYS program, which is presented in Table 3, 4, 5, 6 and 7, respectively. Regarding the model in ANSYS, through the thickness of the plate, 80 layers of isotropic materials are defined and PLANE-183 element is used for the model. The plate is divided into 2400 elements. This solid element is defined by a quadratic 8-node element with two degrees of freedom at each node and translations in the nodal x and y directions with a torsion option for the axisymmetric. For more information about this element, check the Mechanical APDL Element Reference [53].

Table 3 Verification of the results with those of ANSYS for (h/r= 1/10) clamped circular plate (mm)
Full size table
Table 4 Verification of the results with those of ANSYS for (h/r = 1/20) clamped circular plate (mm)
Full size table
Table 5  Verification of the results with those of ANSYS for (h/r = 1/30) clamped circular plate (mm)
Full size table
Table 6  Verification of the results with those of ANSYS for (h/r = 1/40) clamped circular plate (mm)
Full size table
Table 7  Verification of the results with those of ANSYS for (h/r = 1/50) clamped circular plate (mm)
Full size table

The results obtained from the Tables 3, 4, 5, 6 and 7 have been examined for the clamped type of boundary condition. In this study, a computer program is coded in FORTRAN. The modulus of elasticity of the plate varies continually in the thickness direction, where (E1 = 210 × 109 N/m2). To verify the results, several examples of (h/r) are analyzed. According to Tables 3, 4, 5, 6 and 7, the deflection of the plate increases as the porosity coefficient and (h/r) increase. For clamped circular plates, excellent agreement is observed between present study and the result of the ANSYS in the UMCR type, compared with the other two types. For MMCR, the largest difference between the ANSYS and the present study is observed in Table 7 for (h/r = 1/50) at porosity coefficients (0.9), while the biggest difference in SMCR type is seen in Table 3 for h/r = 1/10 at porosity coefficients 0.4.

As shown in Tables 8, 9, 10, 11, 12 and 13, the maximum deflection values are examined for two types of boundary conditions (clamped and roller supported), for various types of porosity coefficients (0.1–0.9) taking into consideration the different h/r values (r = 1–5 m).

Table 8 Maximum vertical displacement values of fixed supported UMCR porous circular plate (mm)
Full size table
Table 9 Maximum vertical displacement values of fixed supported SMCR porous circular plate (mm)
Full size table
Table 10 Maximum vertical displacement values of fixed supported MMCR porous circular plate (mm)
Full size table
Table 11 Maximum vertical displacement values of roller-supported UMCR porous circular plate (mm)
Full size table
Table 12 Maximum vertical displacement values of roller-supported SMCR porous circular plate (mm)
Full size table
Table 13 Maximum vertical displacement values of roller-supported MMCR porous circular plate (mm)
Full size table

The results show that when the value of the property coefficient (e0) increases, so does the value of the maximum deflection for all types of porosity distribution increases. In addition, that when the value of (h/r) increases, the value of maximum deflection increases for all types of porosity distribution. Also, in Tables 8, 9, 10, 11, 12 and 13, for all types of boundary conditions. It is observed, that the highest value of maximum vertical displacement is recorded by roller-supported type in Table 11 for UMCR of porosity coefficient 0.9, h/r (1/50) and the lower value is recorded by clamped type in Table 9 for SMCR of porosity 0.1, h/r (1/10).

As for the clamped circular plate, the highest value of the maximum vertical displacement is observed in Table 8 for UMCR with porosity coefficient (0.9), while the lowest value is seen in Table 9 for (SMCR) type with the porosity coefficient (0.1). Regarding the roller-supported type, the highest value of maximum vertical displacement is recorded by the UMCR porous type in Table 11, while the lowest value is observed for SMCR porous type in Table 12.

As shown in Figs. 3, 4, 5, 6 and 7, a comparison of the maximum vertical displacement of the FGP circular plate for two types of boundary conditions is presented: clamped and roller-supported, respectively. The model is tested for different values of (h/r).

Fig. 3
figure 3

Comparison of maximum vertical displacements for FGP clamped circular plate (h/r = 1/10)

Full size image
Fig. 4
figure 4

Comparison of maximum vertical displacements for FGP clamped circular plate (h/r = 2/10)

Full size image
Fig. 5
figure 5

Comparison of maximum vertical displacements for FGP clamped circular plate (h/r = 3/10)

Full size image
Fig. 6
figure 6

Comparison of maximum vertical displacements for FGP clamped circular plate (h/r = 4/10)

Full size image
Fig. 7
figure 7

Comparison of maximum vertical displacements for FGP clamped circular plate (h/r = 5/10)

Full size image

For the clamped FGP circular plate type, the deflection values obtained for UMCR, SMCR and MMCR material types are presented in Figs. 8, 9, 10, 11 and 12. From the results, an increase in the maximum vertical displacement is observed, in the value of UMCR compared to MMCR and SMCR types, after crossing the porosity coefficient of 0.4. However, it is consistent with the values before this value.

Fig. 8
figure 8

Comparison of the maximum vertical displacements for FGP roller-support circular plate (h/r = 1/10)

Full size image
Fig. 9
figure 9

Comparison of the maximum vertical displacements for FGP roller-support circular plate (h/r = 2/10)

Full size image
Fig. 10
figure 10

Comparison of the maximum vertical displacements for FGP roller-support circular plate (h/r = 3/10)

Full size image
Fig. 11
figure 11

Comparison of the maximum vertical displacements for FGP roller-support circular plate (h/r = 4/10)

Full size image
Fig. 12
figure 12

Comparison of the maximum vertical displacements for FGP roller-support circular plate (h/r = 5/10)

Full size image

For the roller-supported FGP circular plate, as shown in Figs. 13, 14 and 15, it is observed that the maximum vertical displacement of the MMCR is higher than those of the UMCR before the porosity coefficient of 0.6, while it showed the opposite of its performance, with a decrease in its value and an increase in the value of UMCR type after crossing this the point.

Fig. 13
figure 13

Comparison of modulus of elasticity between SMCR, UMCR and MMCR types with porosity coefficients of (0.6)

Full size image
Fig. 14
figure 14

Comparison of modulus of elasticity between SMCR, UMCR and MMCR types with porosity coefficients of (0.4)

Full size image
Fig. 15
figure 15

Comparison between modulus of elasticity of (0.2), (0.4), (0.6) and (0.8) porosity coefficients for uniform porosity material type

Full size image

The comparison of modulus of elasticity between SMCR, UMCR, and MMCR porosity material types with porosity coefficients of 0.6 and 0.4 are presented in Figs. 14 and 15. As it can be seen in Figs. 14 and 15, the lowest elasticity is observed in the uniform porosity type, so the highest expected displacement occurs in this type of material. This harmony is achieved in the results.

The comparison between modulus of elasticity of 0.2, 0.4, 0.6 and 0.8 porosity coefficients for uniform porosity material type show that the lowest elasticity value appeared in the 0.8 porosity coefficient. Consequently, the maximum anticipated displacement occurs in this type of material. Therefore, the deflection of the plate increases as the porosity coefficient (e0) increases.

4 Conclusion

This paper aims to investigate the axisymmetric bending response of functionally graded porous plates with the aid of the finite element method based on the Galerkin’s approach. An eight-node quadrilateral element with two degrees of freedom on each node is implemented in the analysis. The Gauss quadrature is employed for the process of numerical integration. Results of the presented approach are verified with a comparison between the numerical results of (FGM) circular plates and the previous literature, and good agreement is observed. After the validation of the suggested scheme, it is used to investigate the static response of functionally graded porous circular plates for various thickness/radius ratios, porosity coefficients, three different material distributions (UMCR, MMCR, SMCR) and two boundary conditions. Based on the results obtained from the static response of FGP plates, the following important points can be concluded.

  • The deflection of the plate increases as the porosity coefficient (e0) increases.

  • The maximum displacement observed in uniform porosity material type is greater than that of other porosity material types for the values of e0 greater than 0.4 (clamped) and 0.6 (roller supported). For porosity coefficients less than 0.4 (clamped) and 0.6 (roller supported), the MMCR plates have higher vertical displacements.

  • It can be concluded from the results that an increase in the ratio (h/r) leads to an increase in maximum displacement.

  • The effect (h/r) on the maximum displacement of the clamped plate is less than on roller-supported FGP plates.

  • It is believed that the results of the present study can be used as benchmark solutions for future studies on functionally graded porous circular plates.