Thermoelastic Stress and Deformation Analyses of Functionally Graded Doubly Curved Shells

Abstract

In this paper, the authors develop Reissner’s mixed variational theorem (RMVT)-based finite layer methods for the three-dimensional (3D) coupled thermoelastic analysis of simply supported, functionally graded, doubly curved (DC) shells with temperature-independent material properties. A two-phase composite material is considered to form the shell, and its material properties are assumed to obey a power–law distribution of the volume fractions of the constituents through the thickness direction of the shell. The effective material properties are estimated using the Mori–Tanaka scheme. The accuracy and convergence rate of these RMVT-based finite layer methods are validated by comparing their solutions with the quasi 3D and accurate two-dimensional solutions available in the literature.

1. Introduction

In 1984, a Japanese scientist, Niino, proposed the concept of functionally graded (FG) materials for use in thermal barrier materials at the National Aerospace Laboratory of Japan [1,2]. Since then, the development of these FG materials has progressed rapidly, and their application in various advanced industries has become increasingly popular. These FG materials are composed of two- or more-phases of dissimilar materials in such a manner that the material properties vary gradually and continuously in one or more directions. One of the most obvious applications of the FG materials is the ceramic–metal two-phase FG composite materials, which can have the low thermal conductivity property of ceramic materials, as well as the high toughness and strength properties of metallic materials. The FG materials also have considerable design flexibility. Engineers can design the most appropriate material properties, according to a variety of structural performance requirements, by adjusting the spatial distributions of the volume fractions of the constituents. Moreover, because of the continuous spatial distribution of the material properties of the FG material, it can be used to form a variety of single-layered, sandwiched, and multi-layered FG structures to prevent delamination failure, which often occurs at the interface between adjacent layers of conventional laminated composite structures, in which the material properties suddenly change.

Since the above-mentioned FG structures are intended for use in severe thermal environments, the thermoelastic analysis of these structures has attracted considerable attention in order to accurately estimate the thermoelastic stress and deformation induced in them. Some theoretical methodologies and numerical models for the assorted analyses of laminated composite structures have been extended to those of FG structures. Based on the classical lamination plate theory (CLPT), Shaw [3] investigated the thermal residual stresses induced in an infinite substrate-coating plate subjected to a thermal load, in which the material properties of the substrate were considered to be homogeneous, while those of the coating were functionally graded in the thickness direction with a power–law material model [4]. Based on the same plate theory (i.e., CLPT), Chi and Chung [5] and Chung and Chang [6] presented the analytical and finite element solutions of stress and deformation induced in the FG plate subjected to thermomechanical loads, in which the through-thickness distributions of material properties were considered as power–law, sigmoid function, and exponential function models. The CLPT was also extended by Dai and Dai [7] to the thermomechanical bending analysis of FG cylindrical shells with fully simply supported and clamped edges. Zenkour and Alghamdi [8,9] developed a unified shear deformable plate theory (USDPT) for the bending analysis of FG sandwich plates under thermomechanical loads, in which the core was a homogeneous ceramic layer, while the top and bottom face sheets were the FG two-phase material (i.e., the metal–ceramic composite material) layers, the material properties of which obeyed a power–law distribution through the thickness direction of the face sheets according to the volume fractions of the constituents. Numerical solutions obtained from various refined and advanced theories reduced from the USDPT, such as the classical lamination theory (CLT), and first-order, third-order, and sinusoidal shear deformation theories (FSDT, TSDT, and SSDT), were presented, and the results obtained using these theories were compared with one another. The above-mentioned theories can be categorized as the two-dimensional (2D) equivalent single layer theories (ESLTs) based on the principle of virtual displacements (PVD), rather than the Reissner mixed variational theorem (RMVT). In addition, Kreja and Sabik [10] evaluated the performance of some selected ESLTs by comparing their solutions with several benchmark solutions for the static problems of multilayered plates.

Grigolyuk and Kulikov [11,12] presented the geometrically nonlinear and the static analyses of multilayered anisotropic shells, in which the corresponding governing equations associated with the possible boundary conditions were derived on the basis of the Hellinger–Reissner variational principle [13] and the through-thickness distributions of the transverse normal and shear stresses induced in the shells were numerically investigated. Subsequently, Grigolyuk and Kulikov [14] developed a generalized model for laminated composite thin-walled structures, in which a variety of multi-layered shell theories were constructed. The transverse displacement components of these multi-layered shell theories were assumed to nonlinearly vary through their thickness direction, which makes them possible to study the actual through-thickness stress and strain distributions induced in these shells.

In 2003, Carrera [15] proposed Carrera’s unified formulation (CUF), which is able to handle numerous PVD- and RMVT-based plate/shell theories in a unified manner. The 2D refined and advanced plate/shell theories (i.e., a variety of ESLTs) and their corresponding layerwise (LW) versions can be deduced by taking a set of appropriate fundamental nuclei in the CUF, and these can thus be regarded as the special cases of the CUF, such as the global CLT, FSDT, TSDT, and SSDT, as well as the layerwise CLT, FSDT, TSDT, and SSDT [16,17,18]. The CUF has been applied to estimate the thermomechanical behavior of single-layered, sandwiched, and multi-layered plates [19] and shells [20,21], as well as FG plates [22] and shells [23,24,25], in which both of the ESLTs and their corresponding LW versions were considered along with variable orders of expansions of the field variables in the thickness direction, from linear to fourth orders. After the implementation of various 2D PVD- and RMVT-based theories, the conclusion was that the LW models are superior to the ESLT ones, and that the RMVT-based models are superior to the PVD-based ones, based on an identical order used for expansions of the field variables through the thickness direction. In order to extend the application of CUF, Cinefra et al. [26,27] developed a shell finite element with variable through-the-thickness kinematic for the coupled thermoelastic and coupled hygro-thermoelastic analyses of laminated structures. Moreover, numerous finite element analyses of various structures based on the CUF were carried out by Carrera et al. [28]. Fazzolari [29] and Fazzolari and Carrera [30] used the CUF combined with a set of variable kinematics and the Ritz method to study the free vibration characteristics of laminated composite and FG doubly curved shells and to examine the coupled thermoelastic effect on the frequency parameters of multi-layered composite and FG plates, respectively. Based on an advanced hierarchical higher-order ESLTs, Fazzolari [31] investigated the free vibration of P- and S-FG plates, where the volume fractions of the constituents of the plates were assumed to be a power–law function and a sigmoid function varying through the thickness direction of the plates, and where the material properties of the plates were assumed to be dependent upon the temperature variable.

Some three-dimensional (3D) analytical and semi-analytical numerical theories were also developed for the thermomechanical analysis of laminated composite plates/shells and multi-layered FG ones. Reddy and Cheng [32] studied the 3D thermomechanical deformations of FG rectangular plates using an asymptotic approach, in which the power–law material model was used, and the effective material properties were estimated using the Mori–Tanaka scheme [33]. Based on the FSDT combined with a 3D heat conduction analysis, Reddy and Chin [34] presented the dynamic thermoelastic responses of FG cylinders and plates. Jabbari et al. [35,36] presented the 3D solutions for the thermomechanical behavior of FG hollow cylinders under the radially symmetric loads and non-axisymmetric steady-state loads using a power series method, which was also extended by Poultangari et al. [37] to the analysis of FG hollow spheres under non-axisymmetric thermomechanical loads. Within the framework of 3D elasticity theory, Kulikov and Plotnikova [38,39,40] proposed a sampling surfaces (SaS) method for the 3D thermal stress analysis of laminated composite plates and shells, as well as FG plates subjected to a steady-state thermal load, in which the field variables were expanded as the Fourier series functions in the in-surface domain, and discretized by the differential quadrature (DQ) method in the thickness direction by selecting a series of nodal surfaces located at the Chebyshev polynomial nodes in order to provide a uniform convergence of the SaS solutions. In conjunction with the state space and DQ methods, Akbari Alashti and Khorsand [41] presented the 3D DQ solutions for the nonlinear thermoelastic analysis of FG cylindrical shells bonded with the piezoelectric layers on the top and bottom surfaces of the shells. This semi-analytical DQ method was also extended to the quasi-3D nonlinear thermoelastic analysis of thick-walled FG piezoelectric cylinders by Arefi [42], and the thermoelastic analysis of multi-directional FG rectangular plates resting on an elastic foundation by Adineh and Kadkhodayan [43]. A comprehensive survey with regard to the exact and approximate 3D methods for assorted structural analyses of laminated composite structures and multi-layered FG ones in the literature was carried out by Wu et al. [44] and Wu and Liu [45].

Based on Reddy’s refined higher-order shear deformation theory (HSDT) [46] coupled with the von Kármán geometrical nonlinearity (VKGN) term, Shen [47,48] investigated the nonlinear thermal bending responses of FG plates subjected to the thermomechanical loads, in which temperature-dependent material properties were considered, and this approach was also extended to the similar analysis of FG cylindrical panels resting on elastic foundations by Shen and Wang [49]. The above issue was also examined by Mahapatra et al. [50] using a TSDT and considering the Green–Lagrange nonlinear strains. Wu and Ding [51] presented the coupled thermoelectromechanical analysis of sandwiched hybrid FG elastic material and piezoelectric plates.

Due to the fact that RMVT-based models have better performance than the PVD-based ones and that the LW models are better than the ESLT ones in capturing the interlaminar stresses and deformation, in this work we thus extend some existing RMVT-based finite layer methods [52,53,54,55] for laminated composite plates/hollow cylinders to the thermoelastic analysis of multilayered FG doubly curved (DC) shells subjected to thermal loads. In the formulation the material properties of the DC shell were assumed to obey the power–law and sigmoid function models, in which the through-thickness distributions of material properties were assigned as certain specific functions according to the volume fractions of the constituents. The accuracy and convergence rate of the RMVT-based finite doubly-curved layer (FDCL) methods with different orders, which were used for expansions of the field variables in the thickness direction, are examined by comparing their solutions with the accurate ones available in the literature. A parametric study with regard to the influence of various factors on the thermoelastic behavior of functionally graded doubly curved (FGDC) shells is conducted, such as the material–property gradient indices, different material models, aspect ratios, and radii of curvature.

2. Effective Material Properties

In this work, we consider a simply supported, FGDC shell subjected to steady-state thermal loads, as shown in Figure 1a, in which the material properties of the shell are thickness-dependent. In the formulation, the shell is artificially divided into N_l layers with a small thickness for each individual layer, as compared with the in-surface dimensions. In order to derive the formulation, a global doubly curved coordinate system (i.e., $\xi ,\eta$ and $\zeta$ coordinates) is located on the mid-surface of the shell, and a set of local thickness coordinates $z_m(m=1,2,3,...,N_l)$ is located at the mid-surface of each individual layer, as shown in Figure 1b, in which N_l is taken to be five. The in-surface dimensions of the shell in the $\xi \mathrm{and}\eta$ directions are defined as $L_{\xi }\mathrm{and}{L}_{\eta }$, respectively, and the curvature radii of the shell are $R_{\xi }\mathrm{and}{R}_{\eta }$. The thicknesses of each individual layer and the shell are $h_m(m=1,2,...,N_l)$ and h, respectively, while $h={\displaystyle \sum _{m=1}^{N_l}{h}_m}$. The relationship between the global and local thickness coordinates in the mth-layer is $\zeta ={\overline{\zeta }}_m+z_m$, in which ${\overline{\zeta }}_m=({\zeta }_m+{\zeta }_{m-1})/2$, as well as ${\zeta }_m\mathrm{and}{\zeta }_{m-1}$ denote the global thickness coordinates measured from the mid-surface of the shell to the top and bottom surfaces of the mth-layer, respectively.

The material properties of the FGDC shell are assumed to obey a specific function distribution according to the volume fractions of the constituents through the thickness coordinate, such as the power–law distribution. The effective bulk modulus, shear modulus, thermal expansion coefficient and thermal conductivity coefficient of the FGDC shell are estimated using the Mori–Tanaka scheme, and are given as follows [23,33]:

$$B(\zeta )=\left[V_c(\zeta )(B_c-B_m)\right]/\left\{1+\left[1-V_c(\zeta )\right](B_c-B_m)/\left[B_m+(4/3)G_m\right]\right\}+B_m,$$

(1)

$$G(\zeta )=V_c(\zeta )(G_c-G_m)/\left[1+\left[1-V_c(\zeta )\right](G_c-G_m)/(G_m+f_m)\right]+G_m,$$

(2)

$$\alpha (\zeta )=B_c({\alpha }_c-{\alpha }_m)/\left\{\left\{1-\left[B(\zeta )-B_c\right]/\left[B(\zeta )-B_m\right]\right\}B(\zeta )\right\}+{\alpha }_m,$$

(3)

$$\lambda (\zeta )=V_c(\zeta )({\lambda }_c-{\lambda }_m)/\left\{1+\left[1-V_c(\zeta )\right]({\lambda }_c-{\lambda }_m)/(3{\lambda }_m)\right\}+{\lambda }_m,$$

(4)

where $f_m=G_m\left[9B_m+8G_m\right]/\left[6(B_m+2G_m)\right]$, and the subscripts c and m are defined as the material properties of certain ceramic and metal materials, respectively. $B_c,G_c,{\alpha }_c,{\lambda }_c$ and $B_m,G_m,{\alpha }_m,{\lambda }_m$ are the bulk modulus, shear modulus, thermal expansion coefficient, and thermal conductivity coefficient of the ceramic and metal materials, respectively. The relationships between the bulk and shear moduli, as well as Young’s modulus (E) and Poisson’s ratio ($\upsilon$), are $B=E/\left[3(1-2\upsilon )\right]$ and $G=E/\left[2(1+\upsilon )\right]$.

As mentioned above, the FGDC shell is considered to be a two-phase composite one, which is composed of the ceramic and metal materials according to a specific function distribution of volume fractions of the constituents through the thickness direction of the shell, such as the power–law. In the power–law model the through-thickness distribution of the volume fraction of the ceramic material, $V_c(\zeta )$ in Equations (1)–(4), is given by

$$V_c(\zeta )={\left[(\zeta +h/2)/h\right]}^{{\kappa }_p},(-h/2)\le \zeta \le (h/2),$$

(5)

where $V_c=1$, when $\zeta =h/2$, such that $B(\zeta =h/2)=B_c$, $G(\zeta =h/2)=G_c$, $\alpha (\zeta =h/2)={\alpha }_c$, and $\lambda (\zeta =h/2)={\lambda }_c$; as well as $V_c=0$, when $\zeta =-h/2$, such that $B(\zeta =-h/2)=B_m$, $G(\zeta =-h/2)=G_m$, $\alpha (\zeta =-h/2)={\alpha }_m$, and $\lambda (\zeta =-h/2)={\lambda }_m$. ${\kappa }_p$ denotes the material–property gradient index for the power–law model. The distributions of the volume fraction $V_c(\zeta )$ along with the thickness direction of the shell are shown in Figure 2, in which the values of κ_p is taken to be κ_p = 0.2, 0.5, 1, 2, and 5.

3. Heat Conduction Analysis

A heat conduction analysis of the FGDC shell subjected to a steady-state thermal load will first be analyzed, and the determined temperature distribution of the shell domain is then used for the thermomechanical analysis of the shell. The material properties of the FGDC shell are considered to be dependent upon the thickness coordinate and temperature.

The equation of steady-state heat conduction of the shell without heat generation is given as

$$(1/{\gamma }_{\xi })H_{\xi }{,}_{\xi }+(1/{\gamma }_{\eta })H_{\eta }{,}_{\eta }+H_{\zeta }{,}_{\zeta }+\left[(1/{\gamma }_{\xi }{R}_{\xi })+(1/{\gamma }_{\eta }{R}_{\eta })\right]H_{\zeta }=0,$$

(6)

where the commas denote partial differentiation with respect to the suffix variables, $H_k(k=\xi ,\eta \mathrm{and}\zeta )$ denote the heat fluxes in the $\xi ,\eta \mathrm{and}\zeta$ directions, and $H_k{,}_k=(\partial H/\partial k)$. ${\gamma }_{\xi }\mathrm{and}{\gamma }_{\eta }$ are the scale factors in the $\xi \mathrm{and}\eta$ directions, respectively, in which ${\gamma }_{\xi }=1+(\zeta /R_{\xi })$ and ${\gamma }_{\eta }=1+(\zeta /R_{\eta })$.

The relations between the heat flux and temperature change in the $\xi ,\eta \mathrm{and}\zeta$ directions are given as

$$H_{\xi }=-({\lambda }_{\xi }/{\gamma }_{\xi })T{,}_{\xi },$$

(7)

$$H_{\eta }=-({\lambda }_{\eta }/{\gamma }_{\eta })T{,}_{\eta },$$

(8)

$$H_{\zeta }=-({\lambda }_{\zeta })T{,}_{\zeta },$$

(9)

where T denotes the temperature change, which is measured from the room temperature (${\widehat{T}}_0$), and ${\widehat{T}}_0$ = 300 K. In addition, the current temperature variable is defined as $\widehat{T}$, such that $T=\widehat{T}-{\widehat{T}}_0$.

The state space equations related to the heat conduction analysis of the FGDC shell is given as

$$\left\{\begin{array}{c}T{,}_{\zeta }\\ H_{\zeta }{,}_{\zeta }\end{array}\right\}=\left[\begin{array}{cc}0& -{\lambda }_{\zeta }^{-1}\\ ({\lambda }_{\xi }/{\gamma }_{\xi }^2){\partial }_{\xi \xi }+({\lambda }_{\eta }/{\gamma }_{\eta }^2){\partial }_{\eta \eta }& -\left[(1/{\gamma }_{\xi }{R}_{\xi })+(1/{\gamma }_{\eta }{R}_{\eta })\right]\end{array}\right]\left\{\begin{array}{c}T\\ H_{\zeta }\end{array}\right\}.$$

(10)

The temperature changes on the top and bottom surfaces of the FGDC shell are given as

$$T={\overline{T}}^{-}\mathrm{at}\zeta =-h/2,$$

(11)

$$T={\overline{T}}^{+}\mathrm{at}\zeta =h/2,$$

(12)

where ${\overline{T}}^{-}$ and ${\overline{T}}^{+}$ are the applied temperature changes at the bottom and top surfaces of the DC shell.

In this work, the temperature change prescribed on the top and bottom surfaces of the DC shell are expanded as the double Fourier series as ${\overline{T}}^{\pm }={\displaystyle \sum _{\widehat{m}=1}^{\infty }{\displaystyle \sum _{\widehat{n}=1}^{\infty }{\overline{T}}_{\widehat{m}\widehat{n}}^{\pm }}}\sin \overline{m}\xi \sin \overline{n}\eta$, in which $\overline{m}=\widehat{m}\pi /L_{\xi }$ and $\overline{n}=\widehat{n}\pi /L_{\eta }$, as well as $\widehat{m}\mathrm{and}\widehat{n}$ are the half-wave numbers, the values of which are positive integers. In addition, the thermal conditions at the edges are $\widehat{T}$ = 300 K (or T = 0 K).

For the fully simply supported boundary conditions, the thermal field variables can also be expressed as the double Fourier series in the in-surface domain to exactly satisfy the edge conditions, as follows:

$$T(\xi ,\eta ,\zeta )={\displaystyle \sum _{\widehat{m}=1}^{\infty }{\displaystyle \sum _{\widehat{n}=1}^{\infty }{T}_{\widehat{m}\widehat{n}}}}(\zeta )\sin \overline{m}\xi \sin \overline{n}\eta ,$$

(13)

$$H_{\zeta }(\xi ,\eta ,\zeta )={\displaystyle \sum _{\widehat{m}=1}^{\infty }{\displaystyle \sum _{\widehat{n}=1}^{\infty }{H}_{\zeta \widehat{m}\widehat{n}}}}(\zeta )\sin \overline{m}\xi \sin \overline{n}\eta ,$$

(14)

where the symbols of double summations will be omitted in the later work of this paper for brevity.

Substituting Equations (13) and (14) in (10), the authors can obtain

$$\left\{\begin{array}{c}{T}_{\widehat{m}\widehat{n}}{,}_{\zeta }\\ H_{\zeta \widehat{m}\widehat{n}}{,}_{\zeta }\end{array}\right\}=\left[\begin{array}{cc}0& -{\lambda }_{\zeta }^{-1}\\ -\left[({\lambda }_{\xi }{\overline{m}}^2/{\gamma }_{\xi }^2)+({\lambda }_{\eta }{\overline{n}}^2/{\gamma }_{\eta }^2)\right]& -\left[(1/{\gamma }_{\xi }{R}_{\xi })+(1/{\gamma }_{\eta }{R}_{\eta })\right]\end{array}\right]\left\{\begin{array}{c}{T}_{\widehat{m}\widehat{n}}\\ H_{\zeta \widehat{m}\widehat{n}}\end{array}\right\}.$$

(15)

Equation (15) represents two simultaneously first-order differential equations in terms of two variables (i.e., $T_{\widehat{m}\widehat{n}}\mathrm{and}{H}_{\zeta \widehat{m}\widehat{n}}$), while the related coefficients are dependent on the temperature variable and thickness coordinate. A modified Pagano method with an iteration process is thus used to obtain the through-thickness distributions of the temperature changes. The relevant solution process of the modified Pagano method can be found in Wu and Lu [56], and is not repeated here for brevity.

4. Coupled Thermoelastic Analysis

In the former Section, a heat conduction analysis of the FGDC shell under thermal loads is presented. Once the through-thickness distribution of the temperature changes of the DC shell is obtained, an FDCL method can be developed for the thermal deformation and stress analyses of the FGDC shell, and the corresponding formulation is derived as follows.

4.1. Kinematic and Kinetic Assumptions

In an earlier work, Wu and Li [52,54,55] developed an RMVT-based FDCL for the elastic deformation, stress, free vibration, and buckling analyses of simply supported, laminated composite plates and shells. The earlier formulation is then extended to the elastic deformation and stress analyses of simply supported, FGDC shell, in which the primary variables of a typical layer (i.e., the mth-layer) of the shell, of which the domains are in $0\le \xi \le L_{\xi }$ and $0\le \eta \le L_{\eta }$, and $-h_m/2\le z_m\le h_m/2$, are thus given by

$${\overline{f}}^{(m)}(\xi ,\eta ,z_m)={\displaystyle \sum _{i=1}^{n+1}{\psi }_i^{(m)}(z_m)}{\left[f^{(m)}(\xi ,\eta )\right]}_i,$$

(16)

where $f^{(m)}=u_{\xi }^{(m)},u_{\eta }^{(m)},u_{\zeta }^{(m)},{\tau }_{\xi \zeta }^{(m)},{\tau }_{\eta \zeta }^{(m)}\mathrm{and}{\sigma }_{\zeta }^{(m)}$, in which ${\left(u_{\xi }^{(m)}\right)}_i,{\left(u_{\eta }^{(m)}\right)}_i\mathrm{and}{\left(u_{\zeta }^{(m)}\right)}_i$ denote the elastic displacement components at the ith-nodal surface of the mth-layer of the shell, while ${\left({\tau }_{\xi \zeta }^{(m)}\right)}_i$, ${\left({\tau }_{\eta \zeta }^{(m)}\right)}_i$, and ${\left({\sigma }_{\zeta }^{(m)}\right)}_i$ are the transverse shear and normal stress components. ${\psi }_i^{(m)}$ (i = 1, 2,…, and (n + 1)) is the shape function of the corresponding ith-nodal surface. n denotes the related orders used for expansions of each primary variable, and when n = 1, 2, and 3, the FDCL methods are called linear, quadratic, and cubic ones, respectively.

For a typical layer, the linear constitutive equations, which are valid for the orthotropic materials, are given by

$$\left\{\begin{array}{c}{\sigma }_{\xi }^{(m)}\\ {\sigma }_{\eta }^{(m)}\\ {\sigma }_{\zeta }^{(m)}\\ {\tau }_{\eta \zeta }^{(m)}\\ {\tau }_{\xi \zeta }^{(m)}\\ {\tau }_{\xi \eta }^{(m)}\end{array}\right\}=\left[\begin{array}{cccccc}{c}_{11}^{(m)}& c_{12}^{(m)}& c_{13}^{(m)}& 0& 0& 0\\ c_{12}^{(m)}& c_{22}^{(m)}& c_{23}^{(m)}& 0& 0& 0\\ c_{13}^{(m)}& c_{23}^{(m)}& c_{33}^{(m)}& 0& 0& 0\\ 0& 0& 0& c_{44}^{(m)}& 0& 0\\ 0& 0& 0& 0& c_{55}^{(m)}& 0\\ 0& 0& 0& 0& 0& c_{66}^{(m)}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{\xi }^{(m)}\\ {\epsilon }_{\eta }^{(m)}\\ {\epsilon }_{\zeta }^{(m)}\\ {\gamma }_{\eta \zeta }^{(m)}\\ {\gamma }_{\xi \zeta }^{(m)}\\ {\gamma }_{\xi \eta }^{(m)}\end{array}\right\}-\left\{\begin{array}{c}{c}_{\alpha 1}^{(m)}\\ c_{\alpha 2}^{(m)}\\ c_{\alpha 3}^{(m)}\\ 0\\ 0\\ 0\end{array}\right\}{T}^{(m)},$$

(17)

where ${\sigma }_{\xi }^{(m)},{\sigma }_{\eta }^{(m)},...,\mathrm{and}{\tau }_{\xi \eta }^{(m)}$ are the stress components; ${\epsilon }_{\xi }^{(m)},{\epsilon }_{\eta }^{(m)},...,\mathrm{and}{\gamma }_{\xi \eta }^{(m)}$ are the strain components; $c_{ij}^{(m)}$ (i, j = 1–6) are the elastic coefficients, which are variable through the thickness direction of each FGEM individual layer, and are constants in the homogeneous one. $c_{\alpha i}^{(m)}$ (i = 1–3) denote the stress-temperature coefficients, and $c_{\alpha i}^{(m)}=c_{1i}^{(m)}{\alpha }_1^{(m)}+c_{2i}^{(m)}{\alpha }_2^{(m)}+c_{3i}^{(m)}{\alpha }_3^{(m)}$ (i = 1–3), in which ${\alpha }_i^{(m)}$ (i = 1–3) are the thermal expansion coefficients in the $\xi ,\eta \mathrm{and}\zeta$ directions.

The strain–displacement relations for each individual layer are written as follows:

$${\epsilon }_{\xi }^{(m)}=\left[u_{\xi }^{(m)}{,}_{\xi }+\left(u_{\zeta }^{(m)}/R_{\xi }\right)\right]/{\gamma }_{\xi }={\displaystyle \sum _{i=1}^{n+1}\left[{\psi }_i^{(m)}{\left(u_{\xi }^{(m)}\right)}_i{,}_{\xi }+{\psi }_i^{(m)}{\left(u_{\zeta }^{(m)}\right)}_i/R_{\xi }\right]}/{\gamma }_{\xi },$$

(18)

$${\epsilon }_{\eta }^{(m)}=\left[u_{\eta }^{(m)}{,}_{\eta }+\left(u_{\zeta }^{(m)}/R_{\eta }\right)\right]/{\gamma }_{\eta }={\displaystyle \sum _{i=1}^{n+1}\left[{\psi }_i^{(m)}{\left(u_{\eta }^{(m)}\right)}_i{,}_{\eta }+{\psi }_i^{(m)}{\left(u_{\zeta }^{(m)}\right)}_i/R_{\eta }\right]}/{\gamma }_{\eta },$$

(19)

$${\epsilon }_{\zeta }^{(m)}=u_{\zeta }^{(m)}{,}_{\zeta }={\displaystyle \sum _{i=1}^{n+1}\left(D{\psi }_i^{(m)}\right)}{\left(u_{\zeta }^{(m)}\right)}_i,$$

(20)

$${\gamma }_{\eta \zeta }^{(m)}=\left(u_{\zeta }^{(m)}{,}_{\eta }/{\gamma }_{\eta }\right)+u_{\eta }^{(m)}{,}_{\zeta }-\left[u_{\eta }^{(m)}/\left({\gamma }_{\eta }{R}_{\eta }\right)\right]={\displaystyle \sum _{i=1}^{n+1}\left[{\psi }_i^{(m)}{\left(u_{\zeta }^{(m)}\right)}_i{,}_{\eta }/{\gamma }_{\eta }+\left(D{\psi }_i^{(m)}\right){\left(u_{\eta }^{(m)}\right)}_i-{\psi }_i^{(m)}{\left(u_{\eta }^{(m)}\right)}_i/\left({\gamma }_{\eta }{R}_{\eta }\right)\right]},$$

(21)

$${\gamma }_{\xi \zeta }^{(m)}=\left(u_{\zeta }^{(m)}{,}_{\xi }/{\gamma }_{\xi }\right)+u_{\xi }^{(m)}{,}_{\zeta }-\left[u_{\xi }^{(m)}/\left({\gamma }_{\xi }{R}_{\xi }\right)\right]={\displaystyle \sum _{i=1}^{n+1}\left[{\psi }_i^{(m)}{\left(u_{\zeta }^{(m)}\right)}_i{,}_{\xi }/{\gamma }_{\xi }+\left(D{\psi }_i^{(m)}\right){\left(u_{\xi }^{(m)}\right)}_i-{\psi }_i^{(m)}{\left(u_{\xi }^{(m)}\right)}_i/\left({\gamma }_{\xi }{R}_{\xi }\right)\right]},$$

(22)

$${\gamma }_{\xi \eta }^{(m)}=\left(u_{\eta }^{(m)}{,}_{\xi }/{\gamma }_{\xi }\right)+\left(u_{\xi }^{(m)}{,}_{\eta }/{\gamma }_{\eta }\right)={\displaystyle \sum _{i=1}^{n+1}\left[\left({\psi }_i^{(m)}\right){\left(u_{\eta }^{(m)}\right)}_i{,}_{\xi }/{\gamma }_{\xi }+\left({\psi }_i^{(m)}\right){\left(u_{\xi }^{(m)}\right)}_i{,}_{\eta }/{\gamma }_{\eta }\right]},$$

(23)

where the commas denote partial differentiation with respect to the suffix variables, and $\left(D{\psi }_i^{(m)}\right)=d\left({\psi }_i^{(m)}\right)/d{z}_m$.

4.2. An RMVT-Based Weak-form Formulation

The Reissner mixed variational theorem is used to derive the system equations of the FGDC shell for the FDCL methods, and its corresponding energy functional is written in the form of

$$\begin{array}{l}{\Pi }_R={\displaystyle \sum _{m=1}^{N_l}}{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\displaystyle {\iint }_{\Omega }\left[{\sigma }_{\xi }^{(m)}{\epsilon }_{\xi }^{(m)}+{\sigma }_{\eta }^{(m)}{\epsilon }_{\eta }^{(m)}+{\sigma }_{\zeta }^{(m)}{\epsilon }_{\zeta }^{(m)}+{\tau }_{\xi \zeta }^{(m)}{\gamma }_{\xi \zeta }^{(m)}+{\tau }_{\eta \zeta }^{(m)}{\gamma }_{\eta \zeta }^{(m)}+{\tau }_{\xi \eta }^{(m)}{\gamma }_{\xi \eta }^{(m)}-B({\sigma }_{ij})\right]}}{\gamma }_{\xi }{\gamma }_{\eta }d\xi d\eta d{z}_m\\ -{\displaystyle \sum _{m=1}^{N_l}}{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\displaystyle {\int }_{{\Gamma }_{\sigma }}\left({\overline{t}}_{\xi }^{(m)}{u}_{\xi }^{(m)}+{\overline{t}}_{\eta }^{(m)}{u}_{\eta }^{(m)}+{\overline{t}}_{\zeta }^{(m)}{u}_{\zeta }^{(m)}\right)}}d\Gamma d{z}_m\\ -{\displaystyle \sum _{m=1}^{N_l}}{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\displaystyle {\int }_{{\Gamma }_u}\left[(u_{\xi }^{(m)}-{\overline{u}}_{\xi }^{(m)})t_{\xi }^{(m)}+(u_{\eta }^{(m)}-{\overline{u}}_{\eta }^{(m)})t_{\eta }^{(m)}+(u_{\zeta }^{(m)}-{\overline{u}}_{\zeta }^{(m)})t_{\zeta }^{(m)}\right]}}d\Gamma d{z}_m,\end{array}$$

(24)

where $\Omega$ denotes the DC shell domain on the $\xi -\eta$ surface, and ${\Omega }^{\pm }$ denotes the top surface ($\zeta =h/2$) and the bottom one ($\zeta =-h/2$) of the shell, in which ${\Gamma }_{\sigma }\mathrm{and}{\Gamma }_u$ denote the portions of the edge boundary, in which the surface traction and elastic displacement components are prescribed, respectively (i.e., $t_i^{(m)}={\overline{t}}_i^{(m)}$ and $u_i^{(m)}={\overline{u}}_i^{(m)}$, in which i = $\xi$, $\eta$, and $\zeta$); $B({\sigma }_{ij})$ is the complementary energy density function.

In this formulation, we artificially divide the shell into a series of N_l finite layers and take the elastic displacement components and transverse shear and normal stress components as the primary variables subject to variation. A weak-form formulation can be derived by expressing the energy functional in terms of the primary variables, and then letting the first-order variation of the Reissner energy functional be zero, which yields

$$\begin{array}{l}\delta {\Pi }_R={\displaystyle \sum _{m=1}^{N_l}{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\displaystyle {\iint }_{\Omega }}\{{\left(\delta {\mathit{\epsilon }}_p^{(m)}\right)}^T{\mathit{\sigma }}_p^{(m)}+{\left(\delta {\mathit{\epsilon }}_s^{(m)}\right)}^T{\mathit{\sigma }}_s^{(m)}+\delta {\epsilon }_{\zeta }^{(m)}{\sigma }_{\zeta }^{(m)}+{\left(\delta {\mathit{\sigma }}_s^{(m)}\right)}^T\left({\mathit{\epsilon }}_s^{(m)}-Q_{\tau }^{(m)}{\mathit{\sigma }}_s^{(m)}\right)}}\\ +\delta {\sigma }_{\zeta }^{(m)}\left[{\epsilon }_{\zeta }^{(m)}-\left(1/c_{33}^{(m)}\right){\sigma }_{\zeta }^{(m)}+{\left(Q_{\zeta }^{(m)}\right)}^{{\rm T}}{\mathit{\epsilon }}_p^{(m)}-\left(c_{\alpha 3}^{(m)}/c_{33}^{(m)}\right)T^{(m)}\right]\}{\gamma }_{\alpha }{\gamma }_{\beta }d\alpha d\beta d{z}_m\\ -{\displaystyle \sum _{m=1}^{N_l}}{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\displaystyle {\int }_{{\Gamma }_{\sigma }}\left({\overline{t}}_{\xi }^{(m)}\delta u_{\xi }^{(m)}+{\overline{t}}_{\eta }^{(m)}\delta u_{\eta }^{(m)}+{\overline{t}}_{\zeta }^{(m)}\delta u_{\zeta }^{(m)}\right)}}d\Gamma d{z}_m\\ -{\displaystyle \sum _{m=1}^{N_l}}{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\displaystyle {\int }_{{\Gamma }_u}\left[(u_{\xi }^{(m)}-{\overline{u}}_{\xi }^{(m)})\delta t_{\xi }^{(m)}+(u_{\eta }^{(m)}-{\overline{u}}_{\eta }^{(m)})\delta t_{\eta }^{(m)}+(u_{\zeta }^{(m)}-{\overline{u}}_{\zeta }^{(m)})\delta t_{\zeta }^{(m)}\right]}}d\Gamma d{z}_m\\ =0,\end{array}$$

(25)

where the superscript of T denotes the transposition of the matrices or vectors; and ${\Gamma }_u$ and ${\Gamma }_{\sigma }$ stand for the boundary edges, in which the essential and natural conditions are prescribed. In addition,

$${\mathit{\epsilon }}_p^{(m)}={\left[\begin{array}{ccc}{\epsilon }_{\xi }^{(m)}& {\epsilon }_{\eta }^{(m)}& {\gamma }_{\xi \eta }^{(m)}\end{array}\right]}^T=B_1^{(m)}{u}^{(m)}+B_2^{(m)}{w}^{(m)},$$

(26)

$${\mathit{\epsilon }}_s^{(m)}={\left[\begin{array}{cc}{\gamma }_{\xi \zeta }^{(m)}& {\gamma }_{\eta \zeta }^{(m)}\end{array}\right]}^T=B_3^{(m)}{u}^{(m)}+B_4^{(m)}{w}^{(m)},{\epsilon }_{\zeta }^{(m)}=B_5^{(m)}{w}^{(m)},$$

(27)

$${\mathit{\sigma }}_s^{(m)}={\left[\begin{array}{cc}{\tau }_{\alpha \zeta }^{(m)}& {\tau }_{\beta \zeta }^{(m)}\end{array}\right]}^T=B_6^{(m)}{\mathit{\tau }}^{(m)},{\sigma }_{\zeta }^{(m)}=B_7^{(m)}{\mathit{\sigma }}^{(m)},$$

(28)

$${\mathit{\sigma }}_p^{(m)}={\left[\begin{array}{ccc}{\sigma }_{\xi }^{(m)}& {\sigma }_{\eta }^{(m)}& {\tau }_{\xi \eta }^{(m)}\end{array}\right]}^T=Q_p^{(m)}{B}_1^{(m)}{u}^{(m)}+Q_p^{(m)}{B}_2^{(m)}{w}^{(m)}+Q_{\zeta }^{(m)}{B}_7^{(m)}{\mathit{\sigma }}^{(m)}-Q_{\alpha }^{(m)}{T}^{(m)},$$

(29)

$$u^{(m)}=\left[\begin{array}{c}{\left(u_{\xi }^{(m)}\right)}_i\\ {\left(u_{\eta }^{(m)}\right)}_i\end{array}\right],w^{(m)}=\left[{\left(u_{\zeta }^{(m)}\right)}_i\right],{\mathit{\tau }}^{(m)}=\left[\begin{array}{c}{\left({\tau }_{\xi \zeta }^{(m)}\right)}_i\\ {\left({\tau }_{\eta \zeta }^{(m)}\right)}_i\end{array}\right],{\mathit{\sigma }}^{(m)}=\left[{\left({\sigma }_{\zeta }^{(m)}\right)}_i\right],$$

(30)

$$Q_{\tau }^{(m)}=\left[\begin{array}{cc}{\left(c_{55}^{(m)}\right)}^{-1}& 0\\ 0& {\left(c_{44}^{(m)}\right)}^{-1}\end{array}\right],Q_p^{(m)}=\left[\begin{array}{ccc}{Q}_{11}^{(m)}& Q_{12}^{(m)}& 0\\ Q_{12}^{(m)}& Q_{22}^{(m)}& 0\\ 0& 0& Q_{66}^{(m)}\end{array}\right],Q_{\zeta }^{(m)}=\left[\begin{array}{c}\left(c_{13}^{(m)}/c_{33}^{(m)}\right)\\ \left(c_{23}^{(m)}/c_{33}^{(m)}\right)\\ 0\end{array}\right],$$

(31)

$$Q_{\alpha }^{(m)}=\left[\begin{array}{c}{Q}_{\alpha 1}^{(m)}\\ Q_{\alpha 2}^{(m)}\\ 0\end{array}\right],B_1^{(m)}=\left[\begin{array}{cc}\left({\psi }_i^{(m)}/{\gamma }_{\xi }\right){\partial }_{\xi }& 0\\ 0& \left({\psi }_i^{(m)}/{\gamma }_{\eta }\right){\partial }_{\eta }\\ \left({\psi }_i^{(m)}/{\gamma }_{\eta }\right){\partial }_{\eta }& \left({\psi }_i^{(m)}/{\gamma }_{\xi }\right){\partial }_{\xi }\end{array}\right],B_2^{(m)}=\left[\begin{array}{c}\left({\psi }_i^{(m)}/R_{\xi }{\gamma }_{\xi }\right)\\ \left({\psi }_i^{(m)}/R_{\eta }{\gamma }_{\eta }\right)\\ 0\end{array}\right],$$

(32)

$$B_3^{(m)}=\left[\begin{array}{cc}-\left({\psi }_i^{(m)}/R_{\xi }{\gamma }_{\xi }\right)+D{\psi }_i^{(m)}& 0\\ 0& -\left({\psi }_i^{(m)}/R_{\eta }{\gamma }_{\eta }\right)+D{\psi }_i^{(m)}\end{array}\right],B_4^{(m)}=\left[\begin{array}{c}\left({\psi }_i^{(m)}/{\gamma }_{\xi }\right){\partial }_{\xi }\\ \left({\psi }_i^{(m)}/{\gamma }_{\eta }\right){\partial }_{\eta }\end{array}\right],$$

(33)

$$B_5^{(m)}=\left[D{\psi }_i^{(m)}\right],B_6^{(m)}=\left[\begin{array}{cc}{\psi }_i^{(m)}& 0\\ 0& {\psi }_i^{(m)}\end{array}\right],B_7^{(m)}=\left[{\psi }_i^{(m)}\right],$$

(34)

where i = 1, 2, …, (n + 1), $Q_{kl}^{(m)}=c_{kl}^{(m)}-\left(c_{k3}^{(m)}{c}_{l3}^{(m)}/c_{33}^{(m)}\right)$ (k, l = 1 and 2), $Q_{66}^{(m)}=c_{66}^{(m)}$, $Q_{\alpha l}^{(m)}=c_{\alpha l}^{(m)}-\left(c_{\alpha 3}^{(m)}{c}_{l3}^{(m)}/c_{33}^{(m)}\right)$ (l = 1 and 2).

4.3. System Equations and Boundary Conditions

The thermoelastic static behavior of a simply supported, FGDC shell under thermal loads is studied in the following illustrative examples, in which the temperature-independent material properties are used. The applied thermal conditions on the top and bottom surfaces of the shell are mentioned in Equations (13) and (14), and the traction free conditions are given as follows:

$$\left[{\tau }_{\xi \zeta }^{(N_l)}\left(\xi ,\eta ,h/2\right){\tau }_{\eta \zeta }^{(N_l)}\left(\xi ,\eta ,h/2\right){\sigma }_{\zeta }^{(N_l)}\left(\xi ,\eta ,h/2\right)\right]=\left[000\right]\mathrm{on}\mathrm{the}\mathrm{top}\mathrm{surface},$$

(35)

$$\left[{\tau }_{\xi \zeta }^{(1)}\left(\xi ,\eta ,-h/2\right){\tau }_{\eta \zeta }^{(1)}\left(\xi ,\eta ,-h/2\right){\sigma }_{\zeta }^{(1)}\left(\xi ,\eta ,-h/2\right)\right]=\left[000\right]\mathrm{on}\mathrm{the}\mathrm{bottom}\mathrm{surface}.$$

(36)

As mentioned before, the thermal conditions at the edges are $\widehat{T}$ = 300 K (or T = 0 K). In addition, the edge boundary conditions of each individual layer are considered as fully simple supports, which requires that the following quantities are satisfied.

$$u_{\eta }^{(m)}=u_{\zeta }^{(m)}={\sigma }_{\xi }^{(m)}=0\mathrm{at}\xi =0\mathrm{and}\xi =L_{\xi },$$

(37)

$$u_{\xi }^{(m)}=u_{\zeta }^{(m)}={\sigma }_{\eta }^{(m)}=0\mathrm{at}\eta =0\mathrm{and}\eta =L_{\eta },$$

(38)

where $m=1,2,...,N_l$.

The primary variables of each individual layer are expanded as the double Fourier series as follows:

$$\left(u_{\xi }^{(m)},{\tau }_{\xi \zeta }^{(m)}\right)={\displaystyle \sum _{\widehat{m}=1}^{\infty }{\displaystyle \sum _{\widehat{n}=1}^{\infty }\left(u_{\widehat{m}\widehat{n}}^{(m)},{\tau }_{13\widehat{m}\widehat{n}}^{(m)}\right)}}\cos \overline{m}\xi \sin \overline{n}\eta ,$$

(39)

$$\left(u_{\eta }^{(m)},{\tau }_{\eta \zeta }^{(m)}\right)={\displaystyle \sum _{\widehat{m}=1}^{\infty }{\displaystyle \sum _{\widehat{n}=1}^{\infty }\left(v_{\widehat{m}\widehat{n}}^{(m)},{\tau }_{23\widehat{m}\widehat{n}}^{(m)}\right)}}\sin \overline{m}\xi \cos \overline{n}\eta ,$$

(40)

$$\left(u_{\zeta }^{(m)},{\sigma }_{\zeta }^{(m)}\right)={\displaystyle \sum _{\widehat{m}=1}^{\infty }{\displaystyle \sum _{\widehat{n}=1}^{\infty }\left(w_{\widehat{m}\widehat{n}}^{(m)},{\sigma }_{3\widehat{m}\widehat{n}}^{(m)}\right)}}\sin \overline{m}\xi \sin \overline{n}\eta ,$$

(41)

where the edge boundary conditions of the FGDC shell are exactly satisfied.

Substituting Equations (39)–(41) in (25) and then imposing the stationary principle of the Reissner energy functional (i.e., $\delta {\Pi }_R=0$), we finally obtain the system equations of the FGDC shell, as follows:

$${\displaystyle \sum _{m=1}^{N_l}\left[\begin{array}{cccc}{K}_{\mathrm{I}\mathrm{I}}^{(m)}& K_{\mathrm{I}\mathrm{II}}^{(m)}& K_{\mathrm{I}\mathrm{III}}^{(m)}& K_{\mathrm{I}\mathrm{IV}}^{(m)}\\ K_{\mathrm{II}\mathrm{I}}^{(m)}& K_{\mathrm{II}\mathrm{II}}^{(m)}& K_{\mathrm{II}\mathrm{III}}^{(m)}& K_{\mathrm{II}\mathrm{IV}}^{(m)}\\ K_{\mathrm{III}\mathrm{I}}^{(m)}& K_{\mathrm{III}\mathrm{II}}^{(m)}& K_{\mathrm{III}\mathrm{III}}^{(m)}& 0\\ K_{\mathrm{IV}\mathrm{I}}^{(m)}& K_{\mathrm{IV}\mathrm{II}}^{(m)}& 0& K_{\mathrm{IV}\mathrm{IV}}^{(m)}\end{array}\right]}\left[\begin{array}{c}{\overline{u}}^{(m)}\\ {\overline{w}}^{(m)}\\ {\overline{\mathit{\tau }}}^{(m)}\\ {\overline{\mathit{\sigma }}}^{(m)}\end{array}\right]={\displaystyle \sum _{m=1}^{N_l}\left[\begin{array}{c}{F}_{\mathrm{I}}^{(m)}\\ F_{\mathrm{II}}^{(m)}\\ 0\\ F_{\mathrm{IV}}^{(m)}\end{array}\right]},$$

(42)

$$\mathrm{where}{K}_{\mathrm{I}\mathrm{I}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left({\overline{B}}_1^{(m)}\right)}^T{Q}_p^{(m)}{\overline{B}}_1^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m},K_{\mathrm{I}\mathrm{II}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left({\overline{B}}_1^{(m)}\right)}^T{Q}_p^{(m)}{B}_2^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m},$$

$$K_{\mathrm{I}\mathrm{III}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left(B_3^{(m)}\right)}^T{B}_6^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m},K_{\mathrm{I}\mathrm{IV}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left({\overline{B}}_1^{(m)}\right)}^T{Q}_{\zeta }^{(m)}{B}_7^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m},$$

$$K_{\mathrm{II}\mathrm{II}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left(B_2^{(m)}\right)}^T{Q}_p^{(m)}{B}_2^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m},K_{\mathrm{II}\mathrm{III}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left({\overline{B}}_4^{(m)}\right)}^T{B}_6^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m}$$

$$K_{\mathrm{II}\mathrm{IV}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}\left[{\left(B_5^{(m)}\right)}^T{B}_7^{(m)}+{\left(B_2^{(m)}\right)}^T{Q}_{\zeta }^{(m)}{B}_7^{(m)}\right]{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m}$$

$$K_{\mathrm{III}\mathrm{III}}^{(m)}=-{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left(B_6^{(m)}\right)}^T{Q}_{\tau }^{(m)}{B}_6^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m}$$

$$K_{\mathrm{IV}\mathrm{IV}}^{(m)}=-{\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left(B_7^{(m)}\right)}^T{\left(c_{33}^{(m)}\right)}^{-1}{B}_7^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m}$$

$$K_{ij}^{(m)}={\left(K_{ji}^{(m)}\right)}^T(i,j=\mathrm{I},\mathrm{II},\mathrm{III},\mathrm{IV}),$$

$${\overline{B}}_1^{(m)}=\left[\begin{array}{cc}-\overline{m}\left({\psi }_i^{(m)}/{\gamma }_{\xi }\right)& 0\\ 0& -\overline{n}\left({\psi }_i^{(m)}/{\gamma }_{\eta }\right)\\ \overline{n}\left({\psi }_i^{(m)}/{\gamma }_{\eta }\right)& \overline{m}\left({\psi }_i^{(m)}/{\gamma }_{\xi }\right)\end{array}\right],{\overline{B}}_4^{(m)}=\left[\begin{array}{c}\overline{m}\left({\psi }_i^{(m)}/{\gamma }_{\xi }\right)\\ \overline{n}\left({\psi }_i^{(m)}/{\gamma }_{\eta }\right)\end{array}\right],$$

$${\overline{u}}^{(m)}=\left[\begin{array}{c}{\left(u_{\widehat{m}\widehat{n}}^{(m)}\right)}_i\\ {\left(v_{\widehat{m}\widehat{n}}^{(m)}\right)}_i\end{array}\right],{\overline{w}}^{(m)}=\left[{\left(w_{\widehat{m}\widehat{n}}^{(m)}\right)}_i\right],{\overline{\mathit{\tau }}}^{(m)}=\left[\begin{array}{c}{\left({\tau }_{13\widehat{m}\widehat{n}}^{(m)}\right)}_i\\ {\left({\tau }_{23\widehat{m}\widehat{n}}^{(m)}\right)}_i\end{array}\right],\hspace{1em}{\overline{\mathit{\sigma }}}^{(m)}=\left[{\left({\sigma }_{3\widehat{m}\widehat{n}}^{(m)}\right)}_i\right],$$

$$F_{\mathrm{I}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left({\overline{B}}_1^{(m)}\right)}^T{Q}_{\alpha }^{(m)}{T}_{\widehat{m}\widehat{n}}^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m},F_{\mathrm{II}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left(B_2^{(m)}\right)}^T{Q}_{\alpha }^{(m)}{T}_{\widehat{m}\widehat{n}}^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m}$$

$$F_{\mathrm{IV}}^{(m)}={\displaystyle {\int }_{-h_m/2}^{h_m/2}{\left(B_7^{(m)}\right)}^T\left(c_{\alpha 3}^{(m)}/c_{33}^{(m)}\right)T_{\widehat{m}\widehat{n}}^{(m)}{\gamma }_{\xi }{\gamma }_{\eta }d{z}_m}$$

According to the above weak-form formulation, the authors can analyze the thermoelastic static behavior of simply supported, FGDC shells when the thermal loads are applied, which means the primary variables at each nodal surface can be determined, and subsequently the in-surface stress components at the nodal surfaces can be obtained using the resulting primary variables.

5. Illustrative Examples

In the following numerical examples the orders for expansions of the field variables in the RMVT-based FDCL methods are defined to be identical to one another, which means the in- and out-of-surface elastic displacements, and the transverse shear and normal stresses, are expanded as the nth-order Lagrange polynomials in the thickness coordinate of each layer, and n is taken to be 1, 2, and 3 for the linear, quadrature, and cubic FLDC methods, respectively.

5.1. Sandwiched Composite Spherical Shells

The thermoelastic behavior of simply supported, sandwiched composite spherical shells under the thermal loads has been investigated by Brischetto and Carrera [20] using various advanced and refined shell theories deduced from the CUF, while Khare et al. [57] used the FSDT and HSDT, in which the face sheets are made of the graphite-epoxy composite material, and the core is made of an orthotropic elastic material. The through-thickness distribution of the temperature change was directly given as $T=\left(\zeta T_{11}/h\right)\sin \overline{m}\xi \sin \overline{n}\eta$, and $\widehat{m}=\widehat{n}=1$. These solutions are thus used to validate the accuracy and convergence rate of assorted RMVT-based FDCL methods. The lay-up of the sandwiched spherical shell is [0°/core/0°]. The material properties of face sheets are given as $E_1=172.37\mathrm{GPa}$, $E_2=E_3=6.89\mathrm{GPa}$, $G_{12}=G_{13}=3.45\mathrm{GPa}$, $G_{23}=1.38\mathrm{GPa}$, ${\upsilon }_{12}={\upsilon }_{13}={\upsilon }_{23}=0.25$, ${\alpha }_{11}={\alpha }_{33}=0.1\times {10}^{-5}(1{/}^{\circ }\mathrm{C})$, and ${\alpha }_{22}=2\times {10}^{-5}(1{/}^{\circ }\mathrm{C})$, and those of the core are $E_1=E_2=0.28\mathrm{GPa}$, $E_3=3.45\mathrm{GPa}$, $G_{12}=0.11\mathrm{GPa}$, $G_{13}=G_{23}=0.41\mathrm{GPa}$, ${\upsilon }_{12}={\upsilon }_{13}={\upsilon }_{23}=0.02$, ${\alpha }_{11}={\alpha }_{33}=0.1\times {10}^{-6}(1{/}^{\circ }\mathrm{C})$, and ${\alpha }_{22}=2\times {10}^{-5}(1{/}^{\circ }\mathrm{C})$. The material properties in this example are assumed to be temperature-independent.

The geometric parameters of the shell are $L_{\xi }=L_{\eta }=1\mathrm{m}$, $h=0.25\mathrm{m}$, $h_1:h_2:h_3=0.1h:0.8h:0.1h$ and $R_{\xi }=R_{\eta }=5\mathrm{m},10\mathrm{m},20\mathrm{m}\mathrm{and}\infty (\mathrm{plate}).$ A set of dimensionless variables is defined as follows:

$$\left[{\widehat{u}}_{\xi },{\widehat{u}}_{\eta },{\widehat{u}}_{\zeta }\right]=\left[u_{\xi }\left(0,L_{\eta }/2,\zeta \right),u_{\eta }\left(L_{\xi }/2,0,\zeta \right),u_{\zeta }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right)\right]\left(10h/{\widehat{\alpha }}_1{T}_{11}{L}_{\xi }^2\right),$$

(43)

$$\left[{\widehat{\sigma }}_{\xi },{\widehat{\sigma }}_{\eta },{\widehat{\tau }}_{\xi \eta }\right]=\left[{\sigma }_{\xi }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right),{\sigma }_{\eta }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right),{\tau }_{\xi \eta }\left(0,0,\zeta \right)\right]\left(10h/{\widehat{\alpha }}_1{T}_{11}{\widehat{E}}_1{L}_{\xi }\right),$$

(44)

$$\left[{\widehat{\tau }}_{\xi \zeta },{\widehat{\tau }}_{\eta \zeta },{\widehat{\sigma }}_{\zeta }\right]=\left[{\tau }_{\xi \zeta }\left(0,L_{\eta }/2,\zeta \right),{\tau }_{\eta \zeta }\left(L_{\xi }/2,0,\zeta \right),{\sigma }_{\zeta }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right)\right]\left(10h/{\widehat{\alpha }}_1{T}_{11}{\widehat{E}}_1{L}_{\xi }\right),$$

(45)

where ${\widehat{\alpha }}_1=1\times {10}^{-6}(1{/}^{\circ }\mathrm{C})$ and ${\widehat{E}}_1=1\mathrm{GPa}.$

Table 1. Results of convergence and accuracy studies for the dimensionless deflections induced in thick sandwiched spherical shells ([0°/core/0°]) under a thermal load.

Theories $({\mathit{N}}_{\mathit{l}}=2{\mathit{N}}_{\mathit{f}}+{\mathit{N}}_{\mathit{c}})$	${\mathit{R}}_{\mathit{\xi }}/{\mathit{L}}_{\mathit{\xi }}=5$	${\mathit{R}}_{\mathit{\xi }}/{\mathit{L}}_{\mathit{\xi }}=10$	${\mathit{R}}_{\mathit{\xi }}/{\mathit{L}}_{\mathit{\xi }}=20$	${\mathit{R}}_{\mathit{\xi }}/{\mathit{L}}_{\mathit{\xi }}=\mathit{\infty }$ (Plates)
Linear FDCL method
($N_f=1$, $N_c=1$)	4.3444	4.3683	4.3744	4.3764
($N_f=1$, $N_c=2$)	4.3518	4.375	4.3809	4.3828
($N_f=1$, $N_c=4$)	4.3422	4.3656	4.3715	4.3734
($N_f=2$, $N_c=8$)	4.3497	4.3732	4.3791	4.3811
($N_f=4$, $N_c=8$)	4.3498	4.3733	4.3792	4.3811
($N_f=4$, $N_c=16$)	4.3496	4.373	4.3789	4.3809
Quadratic FDCL method
($N_f=1$, $N_c=1$)	4.3493	4.3728	4.3787	4.3806
($N_f=1$, $N_c=2$)	4.3494	4.3728	4.3787	4.3807
($N_f=1$, $N_c=4$)	4.3495	4.373	4.3789	4.3808
($N_f=2$, $N_c=8$)	4.3496	4.373	4.3789	4.3809
Cubic FDCL method
($N_f=1$, $N_c=1$)	4.3497	4.3731	4.379	4.381
($N_f=1$, $N_c=2$)	4.3496	4.373	4.3789	4.3809
($N_f=1$, $N_c=4$)	4.3496	4.373	4.3789	4.3809
($N_f=2$, $N_c=8$)	4.3496	4.373	4.3789	4.3809
CLT [20]	1.8043	1.8025	1.8021	1.8019
FSDT [20]	3.1472	3.1632	3.1672	3.1685
FSDT [57]	3.2618	3.2745	3.2775	3.2784
HSDT [57]	4.2032	4.2343	3.2422	4.2448
ED1 [20]	3.1466	3.1631	3.1672	3.1685
ED2 [20]	3.0306	3.0471	3.0512	3.0525
ED3 [20]	4.1867	4.2308	4.2419	4.2456
ED4 [20]	4.1928	4.236	4.2469	4.2505
EDZ1 [20]	4.3705	4.419	4.4312	4.4352
EDZ2 [20]	4.3228	4.372	4.3843	4.3885
EDZ3 [20]	4.3261	4.3754	4.3878	4.3919
LD1 [20]	4.3417	4.3653	4.3712	4.3732
LD2 [20]	4.342	4.3651	4.3709	4.3729
LD3 [20]	4.3427	4.3658	4.3716	4.3736
LD4 [20]	4.3426	4.3657	4.3715	4.3735

shows the linear-, quadratic-, and cubic-order solutions of the RMVT-based FDCL methods for the out-of-surface displacement components induced at certain positions of the thick sandwiched spherical shell, in which $N_f$ and $N_c$ denote the number of divided layers of face sheets and core, respectively, such that $N_l=2N_f+N_c$. It can be seen in Table 1 that various RMVT-based FDCL methods converge rapidly, and the convergent solutions are obtained when $(N_f,N_c)$ are taken as $(N_f,N_c)=(4,16)$, (2, 8), and (1, 2) for linear, quadratic, and cubic FDCL methods, respectively, on the basis of the same four digits after the decimal point. The performance of these FDCL methods is the cubic FDCL method > quadrature FDCL method > linear FDCL method, in which the symbol “>” means a more accurate result and a faster convergence rate. These convergent solutions are also compared with those obtained using the CLT [20], FSDT [20,57], HSDT [57], ESLTs with first-, second-, third-, and fourth-orders (i.e., ED1, ED2, ED3, and ED4) [20], various ESLTs combined with Murakami’s zig-zag function [58] (i.e., EDZ1, EDZ2, and EDZ3) [20] and layer-wise linear, quadratic, cubic, and quartic theories (i.e., LD1, LD2, LD3, and LD4) [20]. The results show the convergent solutions of FDCL methods are in excellent agreement with the closed-form solutions of layer-wise cubic and quartic theories, and the accuracy of the advanced and refined theories is LD theories > EDZ theories > ED theories when the orders used to expand the field variables through the thickness direction remain the same, in which the symbol “>” means more accurate. In addition, the out-of-surface displacement decreases when the curvature radius of the sandwiched spherical shell becomes smaller, which also means when the shell becomes a deeper shell.

Figure 3 shows the through-thickness distributions of displacement and stress components induced in the sandwiched spherical shells, in which $L_{\xi }=L_{\eta }$, $R_{\xi }=R_{\eta }$, $R_{\xi }/L_{\eta }=5$, and $L_{\xi }/h=5$ (thick shell), 10 (moderately thick shell), and 20 (thin shell). In the case of homogeneous sandwiched spherical shells, the through-thickness distributions of out-of-surface and in-surface displacements appear to be global linear and layer-wise linear functions, respectively, while those of in-surface and out-of-surface stress components are layer-wise linear and layer-wise higher-order polynomial functions. The in-surface stresses are discontinuous when crossing over the interfaces between adjacent layers, because the material properties are dissimilar to each other for the bottom surface of the upper layer and the top surface of the lower layer, which might result in a matrix cracking failure at the face sheet-core interfaces. In addition, the transverse shear and normal stresses are shown to be concentrated at the face sheet-core interfaces, which might cause delamination damage.

5.2. FG Cylindrical Shells

The thermoelastic behavior of a simply supported, FG cylindrical shell under a thermal load has been investigated by Cinefra et al. [25] using various layer-wise theories and a quasi-3D theory deduced from the CUF, in which the through-thickness distribution of the temperature change of the shell was obtained by performing a heat conduction analysis of it. The specified temperature changes on the top and bottom surfaces are given as ${\overline{T}}^{+}=T_{11}\sin \overline{m}\xi \sin \overline{n}\eta$ and ${\overline{T}}^{-}=0$, respectively, in which $\widehat{m}=\widehat{n}=1$. The FG cylindrical shell was made of a two-phase composite material, with the particulate-phase material being ceramic zirconia (ZrO₂) and the matrix-phase material being nickel-based alloy (Monel, 70Ni-30Cu). The volume fractions of the ceramic and metal materials (i.e., $V_c(\zeta )$ and $V_m(\zeta )$) are taken to be $V_c(\zeta )={\left[(\zeta +0.5h)/h\right]}^{{\kappa }_p}$ and $V_m(\zeta )=1-V_c(\zeta )$, in which the subscripts m and c denote the metal and ceramic materials, respectively, while the effective material properties are estimated using the Mori–Tanaka micromechanics scheme [33]. The material properties of these two phase materials [23,32] are given as $B_m=227.24\mathrm{GPa}$, $G_m=65.55\mathrm{GPa}$, ${\alpha }_m=15\times {10}^{-6}(1/\mathrm{K})$, ${\lambda }_m=25\mathrm{W}/\mathrm{mK}$, and $B_c=125.83\mathrm{GPa}$, $G_c=58.08\mathrm{GPa}$, ${\alpha }_c=10\times {10}^{-6}(1/\mathrm{K})$, ${\lambda }_c=2.09\mathrm{W}/\mathrm{mK}$. The material properties in this example are assumed to be temperature-independent. The geometric parameters of the shell are $L_{\xi }=1\mathrm{m}$, $L_{\eta }=(10\pi /3)\mathrm{m}$, $R_{\xi }\to \infty \mathrm{m}$, $R_{\eta }=10\mathrm{m}$, and $R_{\eta }/h=50\mathrm{and}1000$.

The dimensionless variables are defined as follows:

$$\left[{\overline{u}}_{\xi },{\overline{u}}_{\eta },{\overline{u}}_{\zeta }\right]=\left[u_{\xi }\left(0,L_{\eta }/2,\zeta \right),u_{\eta }\left(L_{\xi }/2,0,\zeta \right),u_{\zeta }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right)\right]\left(1/{\overline{\alpha }}_1{T}_{11}{L}_{\xi }\right),$$

(46)

$$\left[{\overline{\sigma }}_{\xi },{\overline{\sigma }}_{\eta },{\overline{\tau }}_{\xi \eta }\right]=\left[{\sigma }_{\xi }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right),{\sigma }_{\eta }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right),{\tau }_{\xi \eta }\left(0,0,\zeta \right)\right]\left(1/{\overline{\alpha }}_1{T}_{11}\overline{B}\right),$$

(47)

$$\left[{\overline{\tau }}_{\xi \zeta },{\overline{\tau }}_{\eta \zeta },{\overline{\sigma }}_{\zeta }\right]=\left[{\tau }_{\xi \zeta }\left(0,L_{\eta }/2,\zeta \right),{\tau }_{\eta \zeta }\left(L_{\xi }/2,0,\zeta \right),{\sigma }_{\zeta }\left(L_{\xi }/2,L_{\eta }/2,\zeta \right)\right]\left(1/{\overline{\alpha }}_1{T}_{11}\overline{B}\right),$$

(48)

where ${\overline{\alpha }}_1=1\times {10}^{-6}(1/\mathrm{K})$ and $\overline{B}=1\mathrm{GPa}$.

Table 2. Results of cubic FDCL methods for various field variables induced in a single-layered FGDC shell under a thermal load.

${\mathit{R}}_{\mathit{\eta }}/\mathit{h}$	Cubic FDCL Methods	${\widehat{\mathit{u}}}_{\mathit{\xi }}(\mathit{\zeta }=0.5\mathit{h})$	${\widehat{\mathit{u}}}_{\mathit{\xi }}(\mathit{\zeta }=-0.5\mathit{h})$	${\widehat{\mathit{u}}}_{\mathit{\zeta }}(\mathit{\zeta }=0.5\mathit{h})$	${\widehat{\mathit{\sigma }}}_{\mathit{\eta }}(\mathit{\zeta }=0.5\mathit{h})$	${\widehat{\mathit{\tau }}}_{\mathit{\xi }\mathit{\zeta }}(\mathit{\zeta }=0)$	${\widehat{\mathit{\sigma }}}_{\mathit{\zeta }}(\mathit{\zeta }=0)$
50	$N_l=10$	−3.6818	0.5008	7.4042	−1444.4	27.552	5.2791
	$N_l=20$	−3.5948	0.4855	7.2241	−1461.3	26.87	5.1512
	$N_l=40$	−3.5685	0.4807	7.1693	−1466.5	26.643	5.0905
	$N_l=80$	−3.5614	0.4794	7.1546	−1467.8	26.5804	5.0744
	$N_l=400$	−3.5591	0.4789	7.1498	−1468.3	26.56	5.0692
	LD2 [25]	−4.162	0.9074	8.8684	−1409.8	−7.3846	319.18
	LD8 [25]	−3.5545	0.488	7.1548	−1470.8	26.664	7.5271
	LD14 [25]	−3.5477	0.4837	7.1361	−1470.4	26.459	5.1982
	Quasi-3D [25]	−3.5466	0.4833	7.1337	−1481.4	26.448	5.0753
1000	$N_l=10$	−1.8435	−0.4301	45.001	−1127.2	−5.4085	0.2023
	$N_l=20$	−1.8033	−0.4214	44	−1150	−5.2768	0.245
	$N_l=40$	−1.7913	−0.4188	43.7	−1156.9	−5.2373	0.24
	$N_l=80$	−1.7881	−0.4181	43.62	−1158.7	−5.2268	0.2393
	$N_l=400$	−1.787	−0.4179	43.593	−1159.3	−5.2233	0.2391
	LD2 [25]	−1.8872	−0.3785	48.034	−1098.6	−6.6837	259.6
	LD8 [25]	−1.7886	−0.4176	43.653	−1159.3	−5.2415	−1.7681
	LD14 [25]	−1.7871	−0.4178	43.6	−1159.2	−5.2262	0.3165
	Quasi-3D [25]	−1.7868	−0.4178	43.59	−1170.2	−5.2242	0.2428

shows the results of cubic FDCL methods for various field variables induced in the shell, in which ${\kappa }_p=2$. The convergent solutions are obtained when the number of divided layers is taken to be 40, and the relative errors of various field variables are below 0.9% as compared with the 400-layer solutions. The convergence rate of various field variables in the case of the shell with $L_{\xi }/h=5\mathrm{and}{R}_{\xi }/h=50$ is slightly slower than that in the case of case of the shell with $L_{\xi }/h=5\mathrm{and}{R}_{\xi }/h=1000$. These convergent solutions are also shown to be in excellent agreement with the solutions of the LD14 method [25] and the quasi-3D theory [25].

Figure 4 shows the through-thickness distributions of assorted field variables induced in the single-layered FG cylindrical shell, in which $L_{\xi }/h=5,R_{\xi }/h=50$, and ${\kappa }_p=2,1\mathrm{and}0.5$. It can be seen in Figure 4a that the distributions of the temperature change of the FG cylindrical shell appear to be higher-order polynomial functions through the thickness direction, rather than the linear function in the case of homogeneous sandwiched shell. The results in Figure 4b–f also show the through-thickness distributions of thermal deformations and thermal stresses are higher-order polynomial functions, and the trend of the distributions of transverse shear and normal stresses is much steeper than that of the distributions of thermal deformations and in-surface stresses.

5.3. FGDC Shells

In the previous two Sections, the FDCL methods with linear-, quadratic-, and cubic-orders are validated by comparing their solutions with the quasi-3D and accurate 2D solutions available in the literature, in which the material properties of the shells are considered to remain at room temperature ($\widehat{T}=300\mathrm{K}$). In this section, the authors investigate the thermoelastic behavior of simply supported, FGDC shells. The DC shell is made of the two-phase composite material, in which the particulate-phase material is ceramic zirconia and the matrix-phase material is metal Ti-6Al-4V, the TD material properties of which are given in the form of

$$P=P_0\left[\left(P_{-1}/\widehat{T}\right)+1+\left(P_1\widehat{T}\right)+\left(P_2{\widehat{T}}^2\right)+\left(P_3{\widehat{T}}^3\right)\right],$$

(49)

where P denotes various material properties, such as E, $\upsilon$, $\alpha$, and $\lambda$, as well as the values of coefficients $P_i(i=-1,0,1,2\mathrm{and}3)$ are listed in

Table 3. Temperature-dependent material properties of metal and ceramic materials (Ti-6Al-4V and zirconia), in which $P=P_0\lfloor \left(P_{-1}/\widehat{T}\right)+1+\left(P_1\widehat{T}\right)+\left(P_2{\widehat{T}}^2\right)+\left(P_3{\widehat{T}}^3\right)\rfloor$, and $\widehat{T}$ denotes the current temperature and its unit is K.

Materials	Properties P(T)	P0	P−1	P1	P2	P3
Zirconia [34]	E (GPa)	244.27	0	−1.371 × 10−3	1.214 × 10−6	−3.681 × 10−10
	ν	0.2882	0	1.133 × 10−4	0	0
	α (1/K)	12.766 × 10−6	0	−1.491 × 10−3	1.00 × 10−5	−6.778 × 10−11
	λ (W/m K)	1.7000	0	1.276 × 10−4	6.648 × 10−8	0
Ti-6Al-4V [34]	E (GPa)	122.56	0	−4.586 × 10−4	0	0
	ν	0.2884	0	1.121 × 10−4	0	0
	α (1/K)	7.5788 × 10−6	0	6.638 × 10−4	−3.147 × 10−6	0
	λ (W/m K)	1.0000	0	1.704 × 10−2	0	0

[34].

The surface thermal conditions of the shell and dimensionless variables are the same as those used in the previous Section, i.e., Equations (39)–(41). The thermal conditions of the top and bottom surfaces of the shell are ${\overline{T}}^{+}=T_{11}\sin \overline{m}\xi \sin \overline{n}\eta$ and ${\overline{T}}^{-}=0$, respectively, in which $\widehat{m}=\widehat{n}=1$ and $T_{11}=300\mathrm{K}$, and the cubic FDCL method is used.

Figure 5 shows the through-thickness distributions of assorted field variables induced in the FGDC shell, in which the geometric parameters of the shell are $L_{\xi }=L_{\eta }=1\mathrm{m},R_{\xi }=R_{\eta }=5\mathrm{m}$, and $h=0.2\mathrm{m}$, in which the power–law model is used and the values of ${\kappa }_p=0.5,1,\mathrm{and}2$. It can be seen in Figure 5 that various displacement and stress components induced in the FGDC shell appear to be higher-order polynomial functions, which obviously differs from those distributions assumed in CST and FSDT.

Figure 6 shows the through-thickness distributions of thermal and elastic field variables induced in the FGDC shells, in which ${\kappa }_p=1.5$, $L_{\xi }=L_{\eta }=1\mathrm{m},R_{\xi }=R_{\eta }=5\mathrm{m}$, and $L_{\xi }/h=5,10\mathrm{and}20$, which are the thick, moderately thick, and thin shells, respectively. It can be seen in Figure 6 that the effects of length-to-thickness ratio on the thermal field variables are minor, while those on the elastic field variables are significant, especially on the elastic deformation, transverse shear stress, and transverse normal stress components. The magnitude order of the various elastic field variables is in-surface stress components (${\overline{\sigma }}_{\xi }\mathrm{and}{\overline{\tau }}_{\xi \eta }$) > transverse shear stress component (${\overline{\tau }}_{\xi \zeta }$) > transverse normal stress component (${\overline{\sigma }}_{\zeta }$). It is clear that the transverse shear and normal stress components increase when the shell becomes thicker. The effects of the transverse shear and normal stresses must thus be taken into account in the thermoelastic analysis of the FGDC shell when it becomes thicker.

6. Concluding Remarks

In this work, the authors develop a unified formulation of RMVT-based FDCL methods with various orders to investigate the thermoelastic behavior of simply supported, laminated composite and FGDC shells, and subjected to thermal loads. Implementation of assorted RMVT-based FDCL methods shows the current FDCL methods converge rapidly, and their convergent solutions are in excellent agreement with the quasi-3D and accurate 2D solutions available in the literature. The through-thickness distributions of thermal and elastic field variables induced in the laminated composite shells appear to be layer-wise linear functions for the displacement and in-surface stress components, while they are higher-order polynomial functions for the transverse shear and stress components. In contrast, the through-thickness distributions of thermal and elastic field variables induced in the FGDC shells appear to be higher-order polynomial functions for all field variables. The results reveal the above-mentioned through-thickness distributions of field variables induced in FGDC shells are inconsistent with the kinematic assumptions of the existing 2D equivalent single-layered shell theories, and development of a specific advanced 2D shell theory for the FG shells is thus needed. Because the current cubic FDCL methods are demonstrated to be superior to the quadrature and linear ones, the layerwise third-order deformation theory may be a potential candidate for the future study of FGDC shells.