Triangular Functionally Graded Porous Moderately Thick Plates—Deformations and Free Vibrations

Abstract

Since the finite element analysis of deformations for plates made of functionally graded materials cannot be carried out with the use of commercial FE packages, aconsistent method of analytical analysis is proposed in the paper. The method of the analysis is based on the application of the weighted residuals and the Bubnov–Galerkin method. The 2D formulation of moderately thick plates is adopted herein for classical and transverse shear deformation plate models (first and third order). Plate deformations subjected to uniform normalpressure and free vibrations are considered. The validity of the analytical model was verified by the comparison of results with FE analysis for isotropic plates. Two correction multipliers were proposed in order to take into account the unsymmetric composition of functionally graded porous plate walls.

1. Introduction

Triangular plates are widely used in different applications in many industrial and engineering fields (e.g., mechanical, aerospace, automotive, etc.). In many cases, the rapid and efficient evaluation of the deformations, the natural vibration frequencies, and associated mode shapes is fundamental in their design.

Due to more complicated forms of triangles and/or polygons (various partial restrictions along the edges that should be included in a significant group of practical problems) the classical methods of solutions for boundary value problems of rectangular plates or shells cannot be applied and transferred easily in order to conduct the analysis for triangular plated structures. Approximate analytical methods present adifficulty regardingthe construction of simple and adequate approximation functions that can be applied to the entire domain of the plate (since the mathematical structure of the boundary conditions becomes complex and the generation of approximating functions becomes very difficult).

In addition, it is also important to note that porous functionally graded materials (FGM) are characterized by a smooth variation of the properties from the bottom to the top surface. The classical commercial finite element codes (e.g., Abaqus, Ansys) do not possess the appropriate type of finite elements that can be adopted to the analysis of structures made of functionally graded materials.

The mentioned limitations of the present analysis are demonstrated in the review papers [1,2]. For triangular plates, an excellent source of reference are the works of Leissa [3,4,5,6].

The methods of computations of deformations and natural frequencies are classical and similar to those presented in the pioneering works of Leissa [3,4,5,6]. More information is demonstrated in Section 3.

The fundamental innovations and significanceof the present work areas follows:

• To illustrate that for FGM triangular plates derived with the use of classical Love-Kirchhoff plate theory (CPT), the influence of material properties can be expressed by two multipliers (M and $\widehat{M}$) similarly as it was verified for FGM rectangular plates (see Muc, Flis [2,7]);
• The validity of the method of computations used herein was verified for isotropic structures only (Section 4) since in the open literature I could not find the results for triangular plates made of functionally graded materials.

However, it should be pointed out that the identical remarks cannot be easily transferred for plates made of unidirectional laminates or plates reinforced with nanostructures (Muc, Muc-Wierzgoń [8]) and similarly for spherical FG shells (Flis, Muc [9]).

2. Formulation of the Problem

Let us consider deformations of triangular plates shown in Figure 1.The plate is made of porous FGM, having either a symmetric (Figure 2a,b) or an unsymmetric pores configuration(Figure 2c)

2.1. Kinematic Relations

The 3-D displacements ${\tilde{U}}_i(x,y,z)$ (i = 1,2,3) at any point of the plate in the x, y, and z directions (see also Figure 1), respectively, are expressed in the following form:

$${\tilde{U}}_1(x,y,z)=u(x,y)-z\frac{\partial w}{\partial x}+F\left(z\right){\psi }_1(x,y)$$

$${\tilde{U}}_2(x,y,z)=v(x,y)-z\frac{\partial w}{\partial y}+F\left(z\right){\psi }_2(x,y)$$

(1)

$${\tilde{U}}_3(x,y,z)=w\left(x,y\right)$$

where u, v, w are the displacements of a generic point on the reference mid-surface, ψ₁, ψ₂ are the rotations of normal to the mid-surface about the y- and x-axes, respectively. The explicit form of the function F(z) is presented in

Table 1. Variants of plate theories.

	Classical Plate Theory (CPT)	First Order Shear Deformation Theory (FSDT)	Higher Order Shear Deformation Theory (HSDT)—Reddy [10]	Higher Order Shear Deformation Theory (HSDT)—Hosseini-Hashemi [11]
F(z)	0	z; the second term in Equation(1) is eliminated	$z\left(1-\frac{4}{3}\frac{z^2}{t^2}\right)$	$z\exp \left(-\frac{2z^2}{t^2}\right)$

Where t denotes the total plate thickness.

.

The distributions of the transverse shear corrections functions F(z) used in the literature [10,11,12,13] are demonstrated for first and third order transverse shear deformation theory. Let us note that for HSDT there is no need to use shear correction factors since the form of the functions leads to the parabolic distributions of transverse shear stresses. The above relations are valid for moderately thick FG plates (i.e., as the thickness, t, to the length of the plate ratio is equal to 0.05–0.1). For thick FG plates it is necessary to take into account deformations of the normal to the plate mid-surface—see [2].

Based on the above expressions, the total, linear 3-D strain tensor can be written as follows:

$${\epsilon }_{ij}\left(x,y,z\right)=\frac{1}{2}\left(\frac{\partial {\tilde{U}}_i}{\partial x_j}+\frac{\partial {\tilde{U}}_j}{\partial x_i}\right), x_1=x, x_2=y,x_3=z, i,j=1,2,3$$

(2)

2.2. Constitutive Equations and Variational Formulation

Eliminating from the considerations the normal strains ε₃₃(x,y,z) and stresses σ₃₃(x,y,z), the constitutive stress-strain relations can be expressed in the following way (the Kelvin-Voigt notation) [14]:

σ_i(x,y,z) = C_ijε_j(x,y,z), i,j = 1,2,…5

(3)

where,

$$[C]=\frac{E(z)}{1-{\nu }^2}\left[\begin{array}{ccccc}1& \nu & 0& 0& 0\\ \nu & 1& 0& 0& 0\\ 0& 0& 0.5(1-\nu )& 0& 0\\ 0& 0& 0& 0.5(1-\nu )& 0\\ 0& 0& 0& 0& 0.5(1-\nu )\end{array}\right],$$

(4)

Young’s modulus, E, varies through thickness direction z and Poission’s ratio, is assumed as constant (see Leissa [15]).

The control of global properties (for instance, density ρ and elastic modulus E variations) are described by the identical function P(z) characterizing the pore distributions (Figure 2):

$$P(z)/P_b=[(P_t/P_b-1)f(z)+1](1-\phi (z))$$

(5)

$$f(z)={\left(\frac{z}{t}+\frac{1}{2}\right)}^n,$$

Type 1: $\phi \left(z\right)=c_p\cos \left(\frac{\pi z}{t}\right)$, Type 2: $\phi \left(z\right)=c_p\cos \left[\frac{\pi }{2}\left(\frac{z}{t}+0.5\right)\right]$, Type 3: $\phi \left(z\right)=c_p\cos \left[\frac{\pi }{2}\left(\frac{z}{t}-0.5\right)\right]$ where the symbols t and b refer to the material properties on the top and bottom surfaces, n is the power index, φ(z) is the distribution of porosity along the thickness direction z, and c_p denotes the chosen variation of porosity.

2.3. Governing Relations

Governing relations are derived with the method discussed by Muc [16]. For 3D structures, let us consider the following functional:

$$H={\displaystyle \underset{{\tau }_1}{\overset{{\tau }_2}{\int }}\left(K-\mathsf{\Pi }\right)}d\mathsf{\Omega }d\tau =0$$

(6)

where: K denotes the total kinetic energy, Π is the total potential energy and τ means the physical time. The total kinetic energy can be written in the following form:

$$K=\frac{1}{2}{\displaystyle {\int }_{\mathsf{\Omega }}\rho (z)\frac{d{\tilde{U}}_i}{d\tau }}\frac{d{\tilde{U}}_i}{d\tau }d\mathsf{\Omega }$$

(7)

and the total potential energy is defined as follows:

$$\mathsf{\Pi }({\tilde{U}}_i)={\displaystyle {\int }_{\mathsf{\Omega }}Wd\mathsf{\Omega }}-{\displaystyle {\int }_{S_T}p}{\tilde{U}}_{3}dS={\mathsf{\Pi }}_{\mathrm{int}}-{\displaystyle {\int }_{S_T}p}{\tilde{U}}_{3}dS$$

(8)

where W is the potential of internal forces (the strain energy density function):

$$W=\frac{1}{2}{\sigma }_{\alpha \beta }{\epsilon }_{\alpha \beta }=\frac{1}{2}\left({\sigma }_{11}{\epsilon }_{11}+{\sigma }_{22}{\epsilon }_{22}+{\sigma }_{33}{\epsilon }_{33}\right)+{\sigma }_{12}{\epsilon }_{12}+{\sigma }_{13}{\epsilon }_{13}+{\sigma }_{23}{\epsilon }_{23}$$

(9)

Ω defines the volume occupied by the deformable body. p is the uniform normal pressure in the global reference system. S_T is the portion of the surface on which the pressure is specified. ρ(z) is the density distribution defined with the use of the relation (5).

According to the classical macro-mechanical approach for composite structures (the homogenization of ply properties), each component of the functional in Equation (8) is evaluated as the sum of contributions from FG material.

Through integrating along the plate thickness z the stiffness matrix components are defined as follows:

$$\begin{gathered} \left[A\right]={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}}\left[C\right]dz, \left[B\right]={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}}\left[C\right]zdz, \left[D\right]={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}}\left[C\right]z^{2}dz, \left[F\right]={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}}\left[C\right]z^{3}dz, \left[H\right]={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}}\left[C\right]z^{4}d,z, \\ \left[I\right]={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}}\left[C\right]z^{5}dz, \left[J\right]={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}}\left[C\right] z^{6}dz \end{gathered}$$

(10)

The terms [F], [H], [I], [J] characterize the influence of the higher order transverse shear deformation theories. The material data in the matrix [C] refer to the principal material coordinates.

Using Hamilton’s Principle

$${\displaystyle \underset{{\tau }_1}{\overset{{\tau }_2}{\int }}\delta \left(K-\mathsf{\Pi }\right)}=0$$

(11)

one can derive the fundamental system of equations describing the deformations of the plate.

Let us notice that the above operations can be carried out in a symbolic way searching for the possible variations of unknown functions u, v, w, ψ₁, ψ₂ with the aid of the operation “Variational Calculus” that can be found in the package Mathematica (the single command Euler Equations). Finally, the variations of the Hamilton functional (11) leads to five differential equations. The first two characterize the in-plane effects and are functions of in-plane displacements u and v. The next three describe bending effects and are directly connected with the description of the bending and free vibrations effects. The final system of equations is expressed in the following way:

$$L_{ij}\left[\begin{array}{c}u\\ v\\ w\\ {\psi }_1\\ {\psi }_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ pw-\widehat{\rho }{\omega }^{2}w\\ 0\\ 0\end{array}\right],i,j=1,2,\dots ,5\widehat{\rho }={\displaystyle \underset{-t/2}{\overset{t/2}{\int }}\rho \left(z\right)}dz$$

(12)

The symbols L_ij denote linear differential operators obtained with the use of the Mathematica package. For FSDT, the explicit form of the operators is presented in Ref. [16].

Two of the frequencies being the solution of free vibrations relations are much higher than the third one. They are referred toas the motion of thickness-shear and thickness-twist, respectively. The primarily lateral mode corresponds to the lowest mode. Therefore, simplified frequency equations which neglect the effects of the rotatory inertia are written in Equation (12).

3. Method of the Solution

In the literature various numerical methods were combined for predicting the behaviors of FGM plates. A broad review of them is presented (e.g., in Refs [17,18]). The authors of Refs [17,18] proposed novel numerical methods for the analysis of FGM and piezoelectric FGM structure (square plates, skew plates, and toruses).

In the present paper the system of differential Equation (12) is considered as the set of linear equations with five unknowns (u, v, w, ψ₁, ψ₂) and constants L_ij. Using the classical Kramer method forsolving linear equations, the system of Equation (12) is reduced to one linear differential variable with respect to normal deflection w. The broader discussion of those problems is presented by Muc, Flis [19]. Since for HSDT and FSDT the final form of the above relation is too complex, we demonstrate the result for CPT:

$${\nabla }^{4}w=-pM-{\omega }^2\widehat{M},{\nabla }^4={\nabla }^2{\nabla }^2,{\nabla }^2=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}M=\frac{A}{-B^2+AD}, \widehat{M}=\widehat{\rho }M$$

(13)

For thin-walled plates, the influence of material properties is described by two parameters, with M and $\widehat{M}$ being the function of the material properties on top and bottom surfaces, power index n, the distribution of porosity along the thickness direction φ(z) and the chosen variation of porosity c_p (see Equation (5)). The variations of the multipliers with the parameters characterizing the porosity are illustrated in Figure 3.

As it may be seen in Figure 3a–d the influence of porosity is similar for different parameters c_p and the P_t/P_b ratio. Thisdemonstrates that the description of the coupling effects via the correction factors M and $\widehat{M}$ is general.

For uniform or symmetric distribution of pores (Figure 2a,b B = 0) the relation (13)is reduced to the case of isotropic structures. Let us note that equating M to 1/D₀ Equation (13) is equivalent to the formulation of isotropic plates. Thus, computing numerically with the help of finite elements deformations or eigenvalues of polygons with the bending stiffness D₀ one can evaluate the appropriate values for structures with non-zero values of the B matrix components.

Different approximation techniques exist that areused to solve deformations of plates (see Leissa [3]) due to the complex form of governing equations in the form of Equation (13) (valid for FG CPT) or much more complicated for FSDT. For HSDT FG plates, the method of weighted residuals is adopted. The basic idea of the method of the weighted residuals is to use a trial function with a set of unknown parameters to approximate the solution:

$$w(x,y)={\displaystyle \sum _{s=1}^S{r}_s{\mathsf{\Phi }}_s}(x,y)$$

(14)

The trial function is evaluated using the polynomial expansion:

$$\mathsf{\Phi }(x,y)={\displaystyle \sum _{m=0}^M{\displaystyle \sum _{n=0}^N{a}_{mn}{x}^n{y}^m}}$$

(15)

Having solved all the boundary conditions, we obtainthe values of the coefficients a_ij. Substituting the values of the coefficients in the general polynomial, we end up with trial function (14). The residual function R is determined as:

$$R(x,y)={\nabla }^{4}w(x,y)+pM+{\omega }^2\widehat{M} \mathrm{for} \mathrm{CPT} \mathrm{or} \mathrm{more} \mathrm{complicated} \mathrm{for} \mathrm{FSDT} \mathrm{and} \mathrm{TSDT}$$

(16)

For the residual function (16) the coefficients r_i in Equation (14) are derived with the use of the Bubnov-Galerkin method.

4. Numerical Example

Let us consider a triangular plate clamped on all sides. For such boundary conditions, the expansion (15) is symmetric with respect to the line y = x and m and n should start from 2.

For plates made of FGM the pressure multiplier M (Figure 3a) reduces the value of the maximal displacement as the thickness t/L ratio increases. The reduction factor is higher for HSDT and it increases with more accurate approximations of deformations.

Let us note that the first mode of natural vibrations is identical to these plotted in Figure 4 and Figure 5 (see Equation (13)). The comparison of natural modes is presented in Figure 6. The results are almost identical for those derived for isotropic plates.

The natural frequency reduction factor is a function of the coefficient n and natural frequencies can be higher or lower than for isotropic plates (Figure 3b) similarly as for rectangular plates [2,7].

The validity of the analytical model was compared with the FE analysis for isotropic plates. The triangular plate deflections and natural vibrations can be evaluated with the use of finite element package NISA II in

Table 2. Comparison of analytical and FE results for different 2D models of triangular clamped isotropic plates (t/L = 0.1, D^t/p = 1, ${\rho }^t/D^t$ = 1).

2D Model	Normal Dimensionless Maximal Displacement– ${\mathit{w}}_{\mathbf{max}}{\mathit{E}}^{\mathit{t}}{\mathit{t}}^3/\left[12\mathit{p}\left(1-{\mathit{\nu }}^2\right)\right]$		Dimensionless Eigenfrequency ${\mathit{\omega }}^2\sqrt{{\mathit{\rho }}^{\mathit{t}}/{\mathit{D}}^{\mathit{t}}}$–First Mode
	Analytical	FE (Nisa II)	Analytical	FE (Nisa II)
CPT	9.516 × 10−5	9.744 × 10−5 NKTP 20	89.4	95.4 NKTP 20
FSDT	9.464 × 10−5	9.681 × 10−5 NKTP32	86.3	93.2 NKTP 32
TSDT	9.412 × 10−5	9.594 × 10−5 NKTP32	85.1	90.4 NKTP 32

. Two types of finite elements were used, namely NKTP 20 (three kinematic variables) for thin-walledstructures (CPT) and NKTP 32 (five kinematic variables (FSDT and TSDT)). The agreement of results between analytical and FE analysis is very good. For moderately thick plates (t/L ≤ 0.1), the influence of transverse shear effects is less significant (see also Figure 6).

Both the results shown in Table 2 and in Figure 7 demonstrate that transverse shear effects reduce the values of maximal deformations and natural frequencies. It is in a good agreement (qualitative not quantitative) with the results obtained in refs [17,20].

5. Concluding Remarks

In the present paper, triangular moderately thick plates made of functionally graded porous materials are considered. The proposed model is based on three variants of 2D approximations of plate deformations (i.e., classical theory, first order shear theory, and third order shear theory). The problem deals with deformations and free vibration analysis. It is solved with the use of the weighting residuals and the Bubnov-Galerkin method. The validity and correctness of the results are verified by the comparison with the FE results derived with the use of NISA II package for isotropic structures.

The asymmetry of material characteristics of FG materials can be taken into account by the proposed two multiplication factors evaluated for classical plate theories.