Three Dimensional Coupling between Elastic and Thermal Fields in the Static Analysis of Multilayered Composite Shells

1  Introduction

The monitoring of temperature gradients in the aerospace structures is one of the main topics in the stress analysis. The study of temperature contribution on strains and stresses is very important for modern aircraft, spaceships, launchers, high-tech propulsion devices, pressure vessels for nuclear energy applications and other installations for industries. In each proposed application, a dedicated structural thermal analysis must be conducted where several features must be included, e.g., heat transfer conditions, transitory and static state thermal stress components, production and remaining stress components, vibrations of heated plate and shell structures, large deflections, post-buckling behavior and analyses of shell and plate geometries [1–4]. When thermal stress analyses are conducted, both the opportune developing of the mechanical models and the accurate definition of the temperature fields are fundamental. In this work, the coupled structural model is related to the application of the three-dimensional equilibrium relations and the three-dimensional Fourier heat conduction relation for shells. In this way, the thermal field is an unknown variable as the displacements along with the three main directions. For this reason, the proposed model is defined as full coupled. The main features about full coupled thermo-mechanical analyses were discussed in [5–11] using the appropriate developing of divergence and gradient relations, constitutive equations, variational statements for the case of linear coupled thermoelastic problem, boundary conditions, proportional relations between the gradient of the thermal unknown and the heat flux, field equations, balance relations for energy and the starting circumstances.

The most recent works in literature about the thermo-elastic investigation of laminated and one-layered shells and plates are discussed in four different sections: the first one is related to 3D exact solutions, the second one is on 3D numerical models, the third one is related to 2D exact solutions and the last part is on 2D numerical models. Therefore, the innovative points of this new three-dimensional full coupled thermo-elastic shell theory for the thermal stress investigation of sandwich and multilayered structures are properly focused.

Bhaskar et al. [12] proposed a closed-form three-dimensional model for plates having an assumed temperature profile as a linear function in the thickness direction. The model was written as an uncoupled and linear thermo-elastic solution based on the classic pure elastic 3D models by Pagano et al. [13–15]. The closed-form three-dimensional plate models in [16–19] employed a calculated through-the-thickness temperature evaluation accomplished via the study of the typical heat conduction equations in static conditions. Kulikov and Plotnikova [17,18] employed the methodology of sample surfaces and they developed exact 3D solutions for laminated composite shells. Governing relations in displacement form for a plate embedding monoclinic layers were written in [19]. Exact solutions for thermal stress investigations of composite multilayered plates were developed in [20]. Brischetto et al. [21,22] developed a 3D uncoupled thermo-elastic model for composite and functionally graded shells and plates having the temperature evaluation defined via the three-dimensional Fourier heat conduction relation. Monge et al. [23] proposed a 3D semi-analytical thermo-mechanical model where the Differential Quadrature Method (DQM) was employed for the derivatives in the thickness direction analogous to the 3D pure elastic shell theory in [24] and [25]. In [26], a 3D theory of elasticity was proposed for functionally graded and carbon nanotube reinforced composite rectangular plates having simply supported sides and thermo-mechanical loads imposed at the external surfaces. The use of the Fourier series expansion in the in-plane directions and the state space methodology along the thickness direction allowed exact solutions for static analysis of plates.

Numerical 3D models are more general than exact 3D models but even more complicated. They permit the analysis of general boundary, loading and lamination conditions. Moleiro et al. [27] presented a mixed layer-wise theory for the full coupled thermo-elastic bending investigation of multilayered plates embedding functionally graded material layers under thermal and/or mechanical loads. Adineh et al. [28] proposed a three-dimensional solution for plates using the differential quadrature method and a temperature evaluation accomplished from the resolution of the three-dimensional heat conduction problem. Thermal stress investigation of plates was proposed in [29] where finite difference methods for the thermal investigation and 3D shell finite elements for the stress calculation were employed. In [30], a semi-loof finite element model for the thermal stress investigation of laminated shells and plates was obtained. The 3D multifield relations for shells of revolution embedding functionally graded piezoelectric layers and subjected to thermal and mechanical loads were obtained in [31]; the related displacement and heat conduction relations included thermal effects. Governing equations, based on the three-dimensional elasticity and heat conduction equations, were proposed in [32] for the thermoelastic investigation of multilayered composite cylinders. Thermoelastic vibrations of circular panels having free and clamped sides, with a thermal shock applied on the external faces, were discussed in [33]. Yeh [34] proposed a system of coupled thermo-elastic relations for vibrations of plates via the Galerkin method. The finite element static analysis of piezo-magnetic-elastic plates exposed to hygrothermal loads was analyzed in [35]. The numerical modeling is obtained via linear coupled constitutive relations for magneto-electro-elastic materials considering the effects connected with thermal and hygroscopic fields and the principle of total potential energy. A three-dimensional finite element theory based on the zig-zag first order sub-laminate approximations was discussed in [36]. The temperature field was calculated using a thermal theory where the trough-the-thickness temperature profile was imposed to linearly change in each layer.

Important simplified assumptions are performed through the thickness direction of the structures in 2D exact models. The 2D exact model for shells proposed in [37] was based on the assumption of a temperature profile through the thickness and the Sander’s theory for spherical shells. An exact 2D analytical solution was presented in [38,39] where the temperature profile was assumed as linear. In [38], the thermal effects were analysed on plates by using refined theories. The thermal effects on cross-ply shells were presented in [39]. An exact 2D formulation in cylindrical coordinates was discussed in [40]. The exact solution was proposed via the Airy stress function written in a Fourier-series form and the temperature distribution was also written in cylindrical coordinates. A layered shear deformation laminated plate theory was presented in [41]. Both in-plane displacements and temperature were assumed to be linear for each discretized layer. In [42], a new higher-order theory was developed for thick laminated plates subjected to several combinations of thermal and mechanical loads. Other analytical two-dimensional refined equivalent single layer and layer-wise models were presented in [43] for laminated and sandwich shell panels to remark the difference between the linear assumed temperature profile and the three-dimensional one computed in the thickness direction. The temperature profile effects were compared in terms of displacement and stress results. Full coupled 2D thermo-elastic analyses of shell and plate structures (that includes the temperature in the governing relations as a primary unknown of the model) were introduced in [44] and [45].

In the framework of two-dimensional numerical theories, the work [46] solved the one-dimensional heat transfer equation to calculate a temperature profile for cylindrical shells where only the material layer effect was included. The Rayleigh-Ritz method applied to classical shell theories was proposed for the thermal stress analysis of multilayered cylindrical panels. The temperature profile was calculated in [47] by solving the three-dimensional variant of heat conduction relation in order to consider both the effects related to the embedded material and the thickness layer. A third-order theory for shells where the transverse normal displacement is constant through the thickness direction was proposed in [47]. Barut et al. [48] proposed a non-linear shell model based on the Midlin theory for the finite element thermal stress investigation. The proposed finite element was locking free and the temperature distribution was defined as non-uniform. Thermally-induced large-deflection behaviors of laminates were investigated in [49] in the case of doubly-curved panels, cylindrical shells and flat panels. The employed finite element methodology was developed via the first-order shear deformation theory. Sandwiches with honeycomb core and composite faces were studied in [50]. The analysis was developed for sandwich configurations loaded with a through-the-thickness temperature gradient combined with pressure loadings. Raju et al. [51] presented a thermal study for tapered columns with rectangular sections via the finite element analysis; the work was useful to understand the post-buckling behaviour of the columns. Daneshjo et al. [52,53] developed a zigzag third-order shear deformation plate theory for the dynamic coupled thermo-elastic investigation of multilayered structures. The analysis was proposed as a mixed finite element model. Ibrahimbegovic et al. [54] developed a coupled thermo-mechanical first-order shear deformation theory for the finite element investigation of shell panels. Lee [55] investigated the coupled 1D thermo-elastic problem in quasi-static conditions for laminated hollow axisymmetric cylinders having clamped edges and boundary conditions which depend on time. The model was built on the basic relations of thermoelastic analysis using a polar coordinate system. Oh et al. [56] presented a full coupled thermo-electro-mechanical zigzag cubic finite element plate theory. A simple triangular finite element with three nodes was used.

The present new 3D coupled thermo-elastic shell model has the following innovative points. It is general for cylindrical and spherical shells, cylinders and plates including different composite, orthotropic and isotropic layers. The closed-form solution of equations is written using simply-supported hypotheses for all the sides as boundary conditions and harmonic structures for thermal and elastic primary unknowns. The system of differential relations is solved using the exponential matrix methodology, along the thickness direction. In this way, a layer-wise theory is possible where equilibrium and compatibility conditions can be easily imposed. Moreover, the included temperature evaluation permits taking into account the thermal thickness layer effect and the thermal material layer effect without separately solving the 3D Fourier heat conduction equation. In the open literature, there are not 3D exact full coupled thermo-elastic models based on the exponential matrix methodology able to analyze several benchmarks using the same tool. The proposed three-dimensional exact coupled thermo-elastic shell theory can be seen as a generalization of the pure elastic model already proposed by Brischetto in [57–59] for free vibration and static analyses of plates and shells with constant radii of curvature. The inclusion of the three-dimensional Fourier heat conduction relation in perpendicular mixed curved coordinates (see [60–63]) to 3D equilibrium relations gives a homogeneous differential equation system that can be solved via the procedure shown in [64] and [65]. The new results are proposed as displacements, stresses and temperature profiles through the thickness direction. They can be employed as benchmarking results for those researchers interested in the working up and testing phase of new three-dimensional and two-dimensional numerical shell theories for the thermal stress investigation of sandwich and laminated composite structures.

2  3D Exact and Coupled Thermoelastic Shell Model

The proposed three-dimensional methodology is able to consider plate and shell geometries by means of a generic and unique theory where several structures can be analysed by means of the same tool. The innovative point with respect to the similar 3D model proposed by Brischetto et al. in [21] is the full coupling considered for the elastic and thermal fields. Therefore, the primary variables of the model are the three displacement components and the scalar temperature profile.

2.1 3D Equilibrium Relations for Spherical Panels

Fig. 1 shows the different proposed geometries in mixed curved orthogonal reference coordinates (α, β, z).

Figure 1: Geometry and applied temperatures for the eight benchmarks

The value h for the thickness is constant in each structure. The origin of these curved coordinates is positioned in an angle, where α and β are parallel to the sideways and located on the centre reference surface Ω0. The middle surface Ω0 is the reference surface for the evaluation of all the geometrical parameters. z is the coordinate in the thickness direction and it is perpendicular to Ω0 and directed towards the top surface. Rα and Rβ are the radii of curvature in the in-plane directions α and β, respectively, and they are considered as constant. a and b are the shell dimensions in α and β directions, respectively. A specific parametric coefficient is calculated for α, β and z directions:

Hα=(1+zRα)=(1+z~−h/2Rα),Hβ=(1+zRβ)=(1+z~−h/2Rβ),Hz=1.(1)

Hα and Hβ are functions of the thickness coordinate z (it varies in the range between −h/2 and +h/2 where the origin is located in correspondence of the Ω0 surface) or z~ (it oscillates in the range between 0 and h with the zero located at the bottom surface).

For shells with constant curvature radii, the parameters Hα and Hβ showed in Eq. (1) are linear functions in the thickness direction z. For a spherical panel having a given lamination sequence and composed by a number NL of physical layers, the governing equations are the three equilibrium equations and the 3D Fourier heat conduction relation. They are valid for each k physical layer:

Hβ∂σααk∂α+Hα∂σαβk∂β+HαHβ∂σαzk∂z+(2HβRα+HαRβ)σαzk=0,(2)

Hβ∂σαβk∂α+Hα∂σββk∂β+HαHβ∂σβzk∂z+(2HαRβ+HβRα)σβzk=0,(3)

Hβ∂σαzk∂α+Hα∂σβzk∂β+HαHβ∂σzzk∂z−HβRασααk−HαRβσββk+(HβRα+HαRβ)σzzk=0,(4)

κ1∗k∂2θ∂α2+κ2∗k∂2θ∂β2+κ3∗k∂2θ∂z2=0.(5)

Eqs. (2)–(5) are valid for all the structures proposed in Fig. 1. The relations are proposed for spherical panels having constant radii of curvature because they are the most general. They allow the formulation for plates, cylinders and cylindrical panels by assuming appropriate values for the two radii of curvatures. Eq. (5) is the 3D Fourier heat conduction equation in steady-state conditions specifically written for spherical shells and already presented by Brischetto et al. in [21]. In a full coupled thermoelastic problem, the three-dimensional Fourier heat conduction relation is directly coupled with the three-dimensional equilibrium relations. The three conduction coefficients in the α, β and z directions and employed for the three-dimensional Fourier heat conduction relation are:

κ1∗k=κ1kHα2,κ2∗k=κ2kHβ2,κ3∗k=κ3k.(6)

where the dependence by z is due to Hα and Hβ.

2.2 Geometrical and Constitutive Equations

The geometries proposed are loaded by means of a sovratemperature field θ(α,β,z)=T−T0, where T0 is a reference temperature (commonly, the environment temperature). The geometrical equations defined in orthogonal mixed curved reference coordinates (α, β, z) are:

[εααkεββkεzzkγβzkγαzkγαβk]=[1Hα∂∂α01HαRα01Hβ∂∂β1HβRβ00∂∂z0∂∂z−1HβRβ1Hβ∂∂β∂∂z−1HαRα01Hα∂∂α1Hβ∂∂β1Hα∂∂α0][ukvkwk]−[μαkμβkμzk000]θk,(7)

where [ɛααk ɛββk ɛzzk γβzk γαzk γαβk]T are the six components of the strain vector for each k layer and they are defined as the algebraic sum of displacement strain and thermal strain components. The strain components are defined via the vector of displacements uk, vk and wk in α, β, z directions and via the sovra-temperature θk using the thermal expansion coefficient μαk, μβk and μzk components described in the structural reference coordinates (α, β, z). The coefficients for thermal expansion definition in the structural reference coordinates are calculated from thermal expansion coefficients μ1k, μ2k and μ3k in the material reference coordinates (1, 2, 3) via the appropriate rotations.

The constitutive equation allows the six strain components to be linked with the six stress components via the elastic coefficient matrix C. The constitutive equation is valid for each layer k:

σk=Ckϵk=Ck(ϵmk−ϵθk),(8)

the stress vector σk=[σααk σββk σzzk σβzk σαzk σαβk]T has 6×1 dimension, the matrix Ck including elastic coefficients has 6×6 dimension and the strain vector (algebraic summation of displacement part ϵmk and thermal part ϵθk) has dimension 6×1. The matrix of elastic coefficients in the global reference coordinates α, β and z for the orthotropic case, with related 0∘ or 90∘ angles, has the following structure:

Ck=[C11kC12kC13k000C12kC22kC23k000C13kC23kC33k000000C44k000000C55k000000C66k].(9)

The analytical form of Eqs. (2)–(5) is only developed for angles of orthotropy 0∘ or 90∘; in fact, these features permit simplifying the mathematical formulation thanks to the feature C16k=C26k=C36k=C45k=0. The substitution of Eq. (7) into Eq. (8) allows to link the stresses with the primary variables of the problem u and θ. The explicit form of Eq. (8) defines the thermo-mechanical coupling coefficients λαk, λβk and λzk in the global coordinates α, β and z. They are defined as follows:

λαk=C11kμαk+C12kμβk+C13kμzk,(10)

λβk=C12kμαk+C22kμβk+C23kμzk,(11)

λzk=C13kμαk+C23kμβk+C33kμzk.(12)

2.3 Solution Procedure

The closed-form solution of the coupled relations (2)–(5) is defined via the imposition of the harmonic form for displacement components and sovra-temperature. This feature means simply supported boundary conditions for all the analyses that will be carried out. The harmonic relations for displacement components and the temperature field are:

uk(α,β,z)=Uk(z)cos(α¯α)sin(β¯β),(13)

vk(α,β,z)=Vk(z)sin(α¯α)cos(β¯β),(14)

<wk(α,β,z),θk(α,β,z)>=<Wk(z),Θk(z)>sin(α¯α)sin(β¯β),(15)

where the two parameters α¯ and β¯ are written as α¯=mπa and β¯=nπb. a and b indicate the shell dimensions in in-plane directions. m and n indicate the half-wave numbers in α and β directions, respectively. Uk(z), Vk(z) and Wk(z) stand for the displacement maximum values. Θk(z) is the sovra-temperature maximum value.

The introduction of the harmonic expression of the displacements and temperature profile (Eqs. (13)–(15)), the geometrical relations and the constitutive equations (Eqs. (7)–(12)) into the coupled thermo-mechanical Eqs. (2)–(5) gives a set of four differential relations written as functions of maximum values for the displacement components and for the sovra-temperature with associated derivatives defined in z. The derivatives performed in α and β directions are exactly defined via the harmonic expressions (13)–(15): therefore, they are known algebraic values. The final system is composed of four differential relations of second order in z. They do not have constant terms because of Hα and Hβ which depend on z. Therefore, each k physical layer must be divided into an appropriate number of artificial layers. As a consequence, a further index j is introduced for the global artificial layers and it varies in the range between 1 and the total number of artificial or mathematical layers M. The number of artificial layers can be calculated from the number of subdivisions of each physical layer. Assuming that each mathematical layer is sufficiently thin, the coefficients Hα and Hβ can be exactly calculated in its middle point. In this way, coefficients Asj (having s from 1 to 19) and coefficients Jrj (with r from 1 to 7) become constant terms. The compact form of this system of partial differential relations of second order in z is:

A1jUj+A2jVj+A3jWj+J1jΘj+A4jU,zj+A5jW,zj+A6jU,zzj=0,(16)

A7jUj+A8jVj+A9jWj+J2jΘj+A10jV,zj+A11jW,zj+A12jV,zzj=0,(17)

A13jUj+A14jVj+A15jWj+J3jΘj+A16jU,zj+A17jV,zj+A18jW,zj+J4jΘ,zj+A19jW,zzj=0,(18)

(J5j+J6j)Θj+J7jΘ,zzj=0.(19)

This system of equations can be decreased to a first order one via the methodology discussed in [66] and [67]. The order of the derivatives performed in z is decreased when the number of unknowns in all j layers is redoubled from 4 (Uj, Vj, Wj, Θj) to 8 (Uj, Vj, Wj, Θj, Uj′, Vj′, Wj′, Θj′). The superscript ′ indicates the partial derivatives in z. Terms Θj and Θj′ are defined in analogy with the displacement amplitudes (and the related derivatives) via the solution of the system of Eqs. (16)–(19). In this case, the temperature evaluation through the thickness direction does not must be “a priori” defined by separately solving the Fourier heat conduction relation or by assuming a linear behaviour. The matrix form of the problem is:

[A6j00000000A12j00000000A19j00000000J7j00000000A6j00000000A12j00000000A19j00000000J7j][UjVjWjΘjUj′Vj′Wj′Θj′]′=[0000A6j00000000A12j00000000A19j00000000J7j−A1j−A2j−A3j−J1j−A4j0−A5j0−A7j−A8j−A9j−J2j0−A10j−A11j0−A13j−A14j−A15j−J3j−A16j−A17j−A18j−J4j000−(J5j+J6j)0000][UjVjWjΘjUj′Vj′Wj′Θj′],(20)

Eq. (20) can be compacted as:

DjXj' =AjXj,(21)

where Xj=[UjVjWjΘjUj′Vj′Wj′Θj′]T and Xj=∂Xj∂z. Superscript T indicates a transposed vector. Eq. (21) can be developed as:

Xj'=Dj−1AjXj,(22)

Xj'=A∗jXj,(23)

with A∗j=Dj−1Aj. The inclusion of the proposed formulation into a Matlab tool is easy because it is sufficient to properly define the different matrices presented in the above formulation. The implemented code can compute both the displacement components and the temperature profile along the thickness direction.

A general homogeneous system of first order differential relations can be proposed as:

dxdt=Ax,(24)

x indicates a vector of dimension M×1 and A indicates a matrix of dimension M×M (consisting of constant terms). Eq. (24) has the following solution:

x(t)=eAtx0.(25)

In analogy, the resolution of Eq. (23) becomes:

Xj(z~j)=e(A∗jz~j)Xj(0).(26)

Eq. (26) is employed to define the unknown vector Xj at the highest point of the layer j when the term A∗∗j=e(A∗jhj) has been calculated in each j artificial layer having thickness hj. The exponential matrix is expanded and calculated for each j layer having thickness hj as:

A∗∗j=e(A∗jhj)=I+A∗jhj+A∗j22!hj2+A∗j33!hj3+…+A∗jNN!hjN,(27)

where I has 8×8 dimension and it is the typical identity matrix. Therefore, Eq. (27) in Eq. (26) gives:

Xj(hj)=A∗∗jXj(0),(28)

an alternative form of Eq. (28) is:

Xtj=A∗∗jXbj,(29)

where Xj(hj) is Xtj and Xj(0) is Xbj. t and b indicate the highest point (top) and the lowest point (bottom) of the j layer, respectively.

Eq. (29) connects the top and bottom unknown vector (primary variables and the related derivatives in z) within the corresponding artificial j layer. To move from the j mathematical layer to the next one, a set of interlaminar requirements must be fulfilled. These interlaminar continuity conditions consider the displacement components, the transverse shear and transverse normal stresses, the temperature and the transverse normal heat flux at each interface. The continuity of the displacement components and sovra-temperature at each interface can be imposed as:

ubj=utj−1,vbj=vtj−1,wbj=wtj−1,θbj=θtj−1.(30)

The requirements specified in Eq. (30) are obviously imposed for the displacement maximum values Uj, Vj, Wj and for the sovra-temperature amplitude Θj. In the same way, the continuity of transverse shear stresses, transverse normal stress and transverse normal heat flux in the z direction are imposed as:

σzzbj=σzztj−1,σαzbj=σαztj−1,σβzbj=σβztj−1,qzbj=qztj−1.(31)

In Eqs. (30) and (31), for each considered unknown, the continuity is established between the evaluation made at the lowest point (b) of the general j-th layer and the evaluation made at the highest point (t) of the (j−1)-th layer. The use of constitutive Eq. (8) and the harmonic forms (13)–(15) allows for obtaining an amplitude form of Eqs. (30) and (31). The methodology has already been discussed in [57–59], the only difference lies in the dimension of the involved matrices and vectors. The amplitude form of Eqs. (30) and (31) is re-proposed in compact matrix composition by using a new transfer matrix:

[UVWΘU′V′W′]bj=[10000000010000000010000000010000T10T20T30000T4T500T600T7T8T9τ100T1000000000τ2]j,j−1[UVWΘU′V′W′]tj−1(32)

The diagonal part composed of 1 indicates the congruence conditions expressed as displacement components and the continuity of temperature written in Eq. (30). The other coefficients Ti and τi indicate the stress and heat flux continuity of Eq. (31) expressed as displacement components and temperature (and related derivatives in z). A matrix expression of Eq. (32) is:

Xbj=Tj,j−1Xtj−1.(33)

Eq. (33) allows linking unknows and related derivatives in z defined at the bottom of the j layer with unknows (and related derivatives in z) defined at the top of the (j−1) layer.

All structures taken into account have simply supported sides, this boundary condition is automatically imposed via the harmonic forms for all the primary unknowns:

θ=w=v=0,σαα=0forα=0,a,(34)

θ=w=u=0,σββ=0forβ=0,b.(35)

In addition, load boundary requirements must be imposed at the top and bottom of the structures. These conditions can be written as:

σzz=0,σαz=0,σβz=0,Θ=T−T0forz=±h/2,(36)

the last equations allow to impose the sovra-temperature values directly at the top and bottom of the structure for the pure thermal stress analysis. Eq. (36) can be rewritten, in matrix form as:

[−C13MHαtMα¯−C23MHβtMβ¯C13MHαtMRα+C23MHβtMRβ−λzM00C33M00−C44MHβtMRβC44MHβtMβ¯00C44M00−C55MHαtMRα0C55MHαtMα¯0C55M00000010000][UVWΘU′V′W′Θ′]tM=[000Θt],(37)

[−C131Hαb1α¯−C231Hβb1β¯C131Hαb1Rα+C231Hβb1Rβ−λz100C33100−C441Hβb1RβC441Hβb1β¯00C44100−C551Hαb1Rα0C551Hαb1α¯0C55100000010000][UVWΘU′V′W′Θ′]b1=[000Θb].(38)

Eqs. (37) and (38) can be compacted as:

BtMXtM=𝒫t,(39)

Bb1Xb1=𝒫b,(40)

subscripts t and b indicate the highest and lowest point, respectively. Superscript M shows the last artificial layer and superscript 1 shows the first artificial layer. Vector 𝒫 includes the impositions related to the mechanical load conditions in the three directions α, β and z and the sovra-temperature. Assuming a classical thermal stress analysis, the mechanical load conditions are set equal to zero.

In order to include Eqs. (39) and (40) into a matrix expression of an algebraic system, it is useful to indicate XtM=XM(hM) in terms of Xb1=X1(0) (displacements, temperature and related derivatives in z~ at the highest point of the last mathematical layer connected with the displacement components, temperature and related derivatives in z~ at the lowest point of the first mathematical layer). These conditions are obtained by introducing recursively Eq. (33) into Eq. (29):

XtM=(A∗∗MTM,M−1A∗∗M−1TM−1,M−2…A∗∗2T2,1A∗∗1)Xb1.(41)

Eq. (41) defines the 8×8 matrix Hm for multilayered structures. This matrix is different from that in [57,59] for the pure elastic investigation; the difference is due to the coupling of the 3D Fourier heat conduction relation with the 3D equilibrium equations. In a compact form, it can be rewritten as:

XtM=HmXb1.(42)

Using Eqs. (42), (39) can be rewritten in terms of Xb1:

BtMHmXb1=𝒫t.(43)

Eqs. (43) and (40) can be now compacted as:

EXb1=𝒫,(44)

where

E=[BtMHmBb1],(45)

and

𝒫=[𝒫t𝒫b].(46)

Matrix E is always 8×8 even if the number M of artificial/mathematical layers changes and a layer wise approach is employed. This matrix changes only in size with respect to the pure elastic investigation in [57]. For this reason, this new model can be defined as the generalization of the pure mechanical model proposed in [57] by Brischetto. The vector 𝒫 now contains the load impositions. The matrix form in Eq. (44) is substantially the same seen in [57–59] for the pure elastic investigation, but the addition of the three-dimensional Fourier heat conduction relation in the 3D equilibrium equations is now considered. The vector of unknows Xb1 has 8×1 dimension.

After the calculation of the displacements at the bottom, Eqs. (33) and (29) can be subsequently employed to define the displacement components and temperature (and connected derivatives made with respect to z) through the global thickness of the shell or plate. Therefore, the unknown trends along the z direction are obtained.

3  Results

The first part of this section shows three preliminary validation results to verify the developed general three-dimensional exact coupled thermo-elastic shell theory. Comparisons with [21] are performed to define the appropriate order of expansion N for the exponential matrix in Eq. (27) and the correct number of artificial layers M for the exact definition of constant curvature terms related to shell geometry. After the authentication of the model, new benchmarking analyses are proposed in the second part of this section. N and M parameters are opportunely chosen and the new results discuss the effects of thickness ratios, structure geometries, lamination sequences, embedded materials and sovra-temperature impositions when a 3D thermal stress investigation of shells and plates is performed. All the geometrical data for the validation results and benchmarks proposed in this section are available in Tables 1 and 2, material data are shown in Table 3. In Table 3, Ei, νij and Gi means, respectively, the Young modulus, Poisson ratio and shear modulus in the three directions of the material reference system (1,2,3). μi and κi are the thermal expansion coefficients and the thermal conduction coefficients in the three directions of the material reference system, respectively. For isotropic materials, Gi is computed thanks to the equation written in the specific row.

Tables 4–14 present four different acronyms used to indicate the different models compared in this section: 3D-u- θ means the 3D full coupled thermo-elastic model where the displacements (indicated as u ) and the sovra-temperature profile (indicated as θ ) are fully coupled and they must be considered as primary unknowns of the problem. 3D() is the general indication of the 3D uncoupled thermo-elastic models presented by Brischetto et al. in [21], where the arguments inside the parentheses denote the variable separately solved and the used Fourier heat conduction relation. Specifically, 3D( θ c , 3D) indicates the 3D uncoupled model where the three-dimensional Fourier heat conduction relation is separately solved to calculate θ c , 3D( θ c , 1D) denotes the 3D uncoupled model where the one-dimensional Fourier heat conduction relation is separately solved to calculate θ c and 3D( θ a ) denotes the “a priori” assumed linear temperature evaluation. In the following subsection, the models will be denoted via the proper acronyms.

3.1 Preliminary Validation Results

A number of three preliminary validation results are employed for the validation of this new three-dimensional general exact coupled thermo-elastic shell model, indicated as 3D-u- θ . A spherical panel, a cylindrical panel and a square plate are analysed for several thickness ratio values, temperature profiles, materials and lamination stacking sequences. These validation results allow the verifiction of the 3D-u- θ model when N = 3 and M = 300 . Therefore, a valid exact three-dimensional thermo-elastic analysis is possible for the eight further benchmark cases discussed and analysed in Section 3.2. In general, as can be seen in the cases presented in this section, convergence occurs for fewer mathematical layers than those employed in the 3D( θ c , 3D) model presented in [21]. Table 1 summarises all the configurations analysed in this section.

The first preliminary validation result considers a square plate having simply-supported sides. The geometrical data can be seen in the second column of Table 1 and the material data can be seen in Table 3, in correspondence to the Composite A column. The results considered as a reference are those obtained via the 3D( θ c , 3D) model proposed by Brischetto et al. [21]. The new 3D-u- θ model includes directly the three-dimensional Fourier heat conduction relation in the system given by three-dimensional equilibrium relations. The proposed new solution employs N = 3 and an appropriate number of artificial layers M for the definition of constant curvature parameters. The temperature profile is considered a primary unknown as the displacements. Table 4 shows the no-dimensional transverse normal displacement and in-plane normal and in-plane shear stress components for several thickness ratio values a / h . The new 3D-u- θ model is perfectly consistent with the reference solution [21] for each thickness ratio a / h ; in fact, the material and thickness layer effects are evaluated in the same way. Only the in-plane normal stress in a very thick plate ( a / h  = 2) shows a small difference between the results. This feature is probably due to the mathematical formulation presented in [21] that used hyperbolic functions to solve the 3D Fourier heat conduction relation. This methodology could fail for very thick structures. However, the difference is always lower than 2%. In addition, for all the cases proposed, M = 102 is always sufficient to obtain exact results, except for the in-plane normal stress for thickness ratio a / h = 2 . The acquired temperature profile considers both the effects connected with the material embedded in the layer and the related thickness.

The second preliminary validation result shows a two-layered isotropic cylindrical shell having simply-supported sides. The geometrical data can be seen in the third row of Table 1 and the material data can be seen in Table 3, in correspondence to the Aluminium column (for the bottom physical layer) and the Titanium column (for the top physical layer). The two physical layers have completely different thermal and elastic properties. The employed reference results come from the 3D( θ c , 3D) model presented in [21]. The 3D-u- θ model adds the 3D Fourier heat conduction equation to the 3D equilibrium equations and it possesses all the peculiarities already discussed in the first validation result. Table 5 indicates transverse and in-plane displacement components for several thickness ratio values R α / h . The 3D-u- θ theory is always identical to the reference theory [21] because both employ the 3D Fourier heat conduction relation: coupled with the 3D equilibrium equations or separately solved in an external tool. For each thickness ratio R α / h , the results are perfectly coincident. For all the proposed cases, a number of artificial layers M = 100 is always sufficient to obtain exact results.

The last and third preliminary validation result proposes a simply-supported spherical shell. The configuration and the geometrical data can be visualized in Table 1. The material used for this validation result is the Composite B material shown in Table 3. The non-dimensional transverse displacement is written as w ¯ = 10 h b 2 μ 1 T 1 w . In it, T 1 = 1 K denotes the gradient of the temperature. The solutions calculated via the 3D-u- θ model are compared with the 3D( θ c , 3D) model, always proposed by Brischetto et al. [21] and already described in the previous two validation results. Table 6 shows no-dimensional transverse normal displacement for several thickness ratio values. The coupled 3D-u- θ model is always coincident with the uncoupled reference solution [21] for each thickness ratio because both models use the three-dimensional Fourier heat conduction relation by including both the effects due to the material and the thickness layer. M = 100 fictitious layers are always enough to calculate the exact results for each investigated thickness ratio.

The next section proposes new benchmarking analyses where different geometrical data, thickness ratios, lamination stacking sequences, materials, temperature profiles and displacement-stress studies are proposed. Even if the preliminary results suggest M = 100 for all the validation results, the conservative choice of M = 300 and N = 3 will be used for these new cases.

3.2 New Benchmarking Analyses

Eight benchmarking analyses are here discussed for cylindrical and spherical panels, cylinders and plates. Temperature is enforced at the outer faces with different amplitudes and half-wave numbers (see Fig. 1 for further information about geometrical data and temperature applications). In these new results, several lamination schemes and materials are considered. Therefore, each configuration can be one-layered or/and multi-layered. The 3D-u- θ model will be used for these new analyses. N = 3 terms in the expansion of the exponential matrix and M = 300 artificial layers will be employed for each proposed 3D coupled result (see the preliminary validation result section where these parameters have been discussed and chosen). The 3D coupled thermal stress analysis will be discussed and compared with past uncoupled 3D results where 3D calculated temperature profiles ( θ c , 3D), 1D calculated temperature profiles ( θ c , 1D) and assumed temperature profiles ( θ a ) are externally defined. These new results will be test cases for those researchers interested in the creation of two-dimensional and three-dimensional analytical and/or numerical shell models for thermal stress analysis. Table 1 summarises all the geometrical configurations and temperature impositions proposed in this section.

The first benchmark represents a square one-layered composite plate with simply-supported sides (see Fig. 1). The geometrical and material data used for this case are written in Table 2 (column B1) and in Table 3 (the material used is Composite C). The considered thickness ratios are highlighted in Table 2. Fig. 2 gives the temperature evaluation in the thickness direction of a thick and a thin plate. In the thin one-layered composite plate, all the temperature profiles (calculated via ( θ c , 3D), calculated via ( θ c , 1D), linearly assumed ( θ a ) and obtained from u- θ as a primary variable) are in accordance. For the thick structure case, the u- θ profile is coincident with the ( θ c , 3D) profile; these two models have different behaviours through the thickness with respect to the ( θ c , 1D) profile and the ( θ a ) profile because they also include the effect connected with the thickness layer and not only with the material properties. Therefore, for a thick plate, the temperature profile is not linear. Table 7 shows the in-plane and transverse displacement components and the in-plane normal, transverse shear and transverse normal stress components at different thickness coordinate values when thickness ratios vary. For the thin structure cases, the four proposed three-dimensional models give coincident results because both effects connected with the thickness layer and the embedded material are zero (in addition, it can be seen the accordance between the 3D full coupled model and the 3D uncoupled model based on the 3D Fourier heat conduction relation). For very thick and moderately thick structures, the 3D( θ a ) model and the 3D( θ c , 1D) are completely in accordance: for both the thickness layer effect is negligible. The 3D( θ c , 3D) and the 3D-u- θ models give identical results because they are able to also include the effect of the thickness layer. Fig. 3 gives the in-plane and transverse displacement components and two stress components in the thickness direction of a thick (a/h = 10) structure. In-plane displacement is different from the linear one and transverse displacement is different from constant because of the considered big value for the thickness. The transverse shear stress is continuous and symmetric. In fact, only one layer is employed in the plate. The transverse normal stress fulfills the boundary load requirements at the external surfaces ( σ z z t = σ z z b = P z = 0 ).

Figure 2: Benchmark one, temperature profiles for two different a / h ratios (thick and thin) of the composite square plate embedding one 0 ∘ lamina. Maximum value of θ ( α , β , z ) is investigated in the centre of the square plate ( a /2, b /2)

Figure 3: Benchmark one, displacement components and stress components for a moderately thick composite square plate embedding one 0 ∘ lamina ( a / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w and σ z z in ( a /2, b /2); u and σ α z in (0, b /2)

The second benchmark shows a square composite plate, equally divided in three layers, having simply-supported sides (in Fig. 1 the structure is shown). Geometrical data are available in Table 2 and material data are available in the column Composite C of Table 3. The sovra-temperature profiles in the thickness direction shown in Fig. 4 do not change with respect to the previous benchmark because the two sides have the same value ( a = b ). Therefore, when the fibre orientation moves from 0 ∘ to 90 ∘ , the κ 3 coefficient does not change and, even if the κ 1 and κ 2 coefficients are exchanged, they multiply the same in-plane dimension value and the same half-wave number, because of the squared geometry. Table 8 gives the results stated as displacement and stress components. For thin structures, the four three-dimensional models are comparable. For thicker plates, the 3D( θ c , 3D) and 3D-u- θ models are in accordance and they show bigger differences with respect to the other two models because they are able to consider the thickness layer effects (which means non linear temperature profile for thick layers). Fig. 5 gives the displacement and stress components for the moderately thick ( a / h = 10 ) plate. In-plane displacement has the conventional zigzag behaviour seen in laminated orthotropic plates with three changes of slope; the transverse displacement is not constant along the thickness direction because of the thick geometry and the different material properties at the interface. The continuity of displacements at each layer interface was possible thanks the correct compatibility requirements. In-plane stress σ α α is not continuous at each interface between two adjacent layers. The transverse shear stress σ α z is continuous; in fact, the equilibrium conditions were correctly imposed. Moreover, it satisfies the free mechanical load requirements at the outer surfaces.

Figure 4: Benchmark two, temperature profiles for two different a / h ratios (thick and thin) of the composite square plate embedding three 0 ∘ / 90 ∘ / 0 ∘ laminae. Maximum value of θ ( α , β , z ) is investigated in the centre of the square plate ( a /2, b /2)

Figure 5: Benchmark two, displacement components and stress components for a moderately thick composite square plate embedding three 0 ∘ / 90 ∘ / 0 ∘ laminae ( a / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w and σ α α in ( a /2, b /2); u and σ α z in (0, b /2)

The third benchmark shows an isotropic homogeneous cylinder having simply supported sides (see Fig. 1). All the characteristics are written in Table 2, for what concerns the geometrical data, and in Table 3, for what concerns the material employed (Aluminium Al2024). The temperature profiles are shown in Fig. 6 for a thick and a thin structure. In this benchmark, the effect connected with the thickness layer for thick cylinders is negligible due to the symmetry and the rigidity of the geometry. These conclusions are still valid for the displacement and stress components discussed in Table 9; here, 3D( θ a ), 3D( θ c , 1D), 3D( θ c , 3D) and 3D-u- θ models give the same results for thin cylinders. Small discrepancies are remarked for very thick structures because the thickness layer effect evaluated via the 3D( θ c , 3D) and the 3D-u- θ model is not so important. The main reason of this behaviour is the big stiffness of the closed and symmetric cylinder. Fig. 7 gives displacements and stresses through the thickness direction of the one-layered isotropic cylinder. The considered cylinder is thick and this feature gives not linear in-plane displacement and not constant transverse normal displacement. Transverse normal stress fulfills the imposed outer mechanical load requirements.

Figure 6: Benchmark three, temperature profiles for two different R α / h ratios (thick and thin) of the cylinder embedding a single isotropic Aluminium (Al2024) layer. Maximum value of θ ( α , β , z ) is investigated in the centre of the cylinder at ( a /2, b /2)

Figure 7: Benchmark three, displacement components and stress components for a moderately thick cylinder embedding a single isotropic Aluminium (Al2024) layer ( R α / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w, σ β β and σ z z in ( a /2, b /2); u in (0, b /2)

The benchmark number four discusses an isotropic two-layered cylinder with simply supported sides (see Fig. 1 to visualize the case). Geometrical data are available in Table 2 in the column B4 and the mechanical properties of the material are written in Table 3 (for the bottom physical layer the reference column is Aluminium, for the top the reference column is Titanium). The results in Fig. 8 demonstrate that the temperature evaluation is never linear even when the cylinder is thin (except for the ( θ a ) case). This behavior is connected with the different thermal properties of the two embedded physical layers, which means different slopes of the temperature profile. The thick cylinder gives a small effect connected with the thickness of the layer because of the same reasons already proposed in the benchmark three about the simmetry and rigidity of this type of structure. As in the previous case, the 3D-u- θ model provides the same results calculated with the 3D( θ c , 3D) model because they consider both the thickness layer and the material layer effects. Displacement and stress components for several thickness ratios R α / h are shown in Table 10: the 3D( θ a ) results are always not correct due to the non-linearity of the sovra-temperature evaluation in the thickness direction of the structure. Displacement and stress components evaluation through the thickness direction shown in Fig. 9 remarks the effects connected with the presence of two layers embedding different materials. In-plane displacement, transverse displacement and transverse shear stress components are continuous; in fact, congruence and equilibrium requirements have been appropriately included in the models. In-plane stress shows an interruption at the boundary between the Aluminium alloy layer and the Titanium alloy layer. Typical zigzag effects due to the transverse anisotropy are shown in Fig. 9.

Figure 8: Benchmark four, temperature profiles for two different R α / h ratios (thick and thin) of the cylinder embedding two isotropic Aluminium/Titanium (Al2024/Ti22) layers. Maximum value of θ ( α , β , z ) is investigated in the centre of the cylinder at ( a /2, b /2)

Figure 9: Benchmark four, displacement components and stress components for a moderately thick cylinder embedding two isotropic Aluminium/Titanium (Al2024/Ti22) layers ( R α / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w, σ β β in ( a /2, b /2); v and σ β z in ( a / 2 , 0)

The fifth benchmark analyses an isotropic one-layered cylindrical panel having simply supported sides (see Fig. 1). The mechanical and thermal properties E i , ν i j , G i , μ i and k i are shown in the column Titanium of Table 3. The geometrical data can be visualized in Table 2, under the column B5. The four temperature profiles in Fig. 10 are formally the same for thin shells because only one isotropic layer is embedded. For the thick structure cases, the calculated temperature profiles by means of the ( θ c , 3D) model and the u- θ model are coincident. They are not like the linear behaviour due to the effect connected with the thickness layer. This behaviour is also proved by the results shown in Table 11 for several thickness ratio values: the results obtained via the four three-dimensional theories are closely associated for thin structures. For this specific case, the imposed temperature has a cylindrical behaviour. Therefore, displacement v and stresses σ α β and σ β z are zero at each coordinate through the thickness direction. Fig. 11 discusses the displacement and stress components that are different from zero for the thick ( R α / h = 10 ) structure. The typical behaviors related to the presence of a single layer in isotropic material and a very thick structure are distinctly demonstrated (continuity for displacements and in-plane normal stress σ α α ). Stress σ α z also fulfils the outside boundary load requirements.

Figure 10: Benchmark five, temperature profiles for two different R α / h ratios (thick and thin) of the cylindrical shell embedding one isotropic Titanium (Ti22) layer. Maximum value of θ ( α , β , z ) is investigated in the centre of the cylindrical shell at ( a /2, b /2)

Figure 11: Benchmark five, displacement components and stress components for a moderately thick cylindrical shell embedding one isotropic Titanium (Ti22) layer ( R α / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w, σ α α in ( a /2, b /2); u and σ α z in (0, b /2)

The benchmark number six analyses the thermal stress investigation of a sandwich cylindrical panel (in Fig. 1 is represented the specific geometry). The geometrical data of the case are presented in Table 2. The skins material is made of Aluminium and the core material is PVC (the mechanical properties are written in Table 3, both for skins and core). Results of Fig. 12 show how the calculated temperature behaviour is mandatory for both thick and thin sandwich shells. The sovra-temperature evaluation is always different from the linear behaviour and it has a typical zigzag form because the thermal properties of the faces are different from those of the core. When a thick sandwich panel is considered, the sovra-temperature evaluation obtained from the 3D( θ c , 1D) model is not adequate to capture the effects connected with the thickness layer. In this benchmark, the use of the 3D( θ c , 3D) model or the 3D-u- θ model is mandatory. The differences between the four proposed temperature profiles are evident displacement and stress results at different thickness positions for several thickness ratios are presented in Table 12. The two models, based on the 3D Fourier heat conduction relation, provide the same results. The 3D( θ c , 1D) model gives the same results of 3D( θ c , 3D) and 3D-u- θ models for thin shells. However, 3D( θ c , 1D) results are not correct for thick shell structures because the effect connected with the thickness layer has been discarded. Results proposed via the 3D( θ a ) model are always wrong because the sovra-temperature evaluation is not linear in any case. Displacement u and stresses σ α β and σ β z are zero because the half-wave number m in the applied sovra-temperature evaluation is zero. Fig. 13 for thick sandwich panel ( R α /h = 10) gives the classical zigzag behaviour for displacement and stress components in the case of sandwich structures. In-plane normal stress shows discontinuities between each skin/core interface. Transverse normal stress components continuously vary at each interface because the equilibrium requirements have been correctly implemented. They also satisfy the conditions connected with the external boundary loads.

Figure 12: Benchmark six, temperature profiles for two different R α / h ratios (thick and thin) of the sandwich cylindrical shell with Aluminium faces and PVC core. Maximum value of θ ( α , β , z ) is investigated in the centre of the cylindrical shell at ( a /2, b /2)

Figure 13: Benchmark six, displacement components and stress components for a moderately thick sandwich cylindrical shell with Aluminium faces and PVC core ( R α / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w, σ α α and σ z z in ( a /2, b /2); u in (0, b /2)

The benchmark seven proposes an isotropic one-layered spherical panel having simply-supported sides (see Fig. 1). Table 2 contains the geometrical data and the column Steel of Table 3 contains the mechanical and thermal properties. Sovra-temperature profiles in Fig. 14 demonstrate how their behaviours are linear in the thickness direction for a thin structure. The temperature has to be defined via the 3D-u- θ model or via the 3D( θ c , 3D) solution in the case of a thick structure because taking into account the effect connected with the thickness layer is mandatory to obtain exact results. These conclusions are attested by the displacement and stress components discussed in Table 13; in it, 3D( θ c , 1D) and 3D( θ a ) theories are always very similar for each R α / h and explored quantity. 3D-u- θ and 3D( θ c , 3D) are always in accordance for each thickness ratio and they always give exact results. Fig. 15 gives the classical displacement and stress behaviors in the thickness direction of a thick spherical panel made of a single isotropic layer.

Figure 14: Benchmark seven, temperature profiles for two different R α / h ratios (thick and thin) of the spherical shell embedding one isotropic Steel layer. Maximum value of θ ( α , β , z ) is investigated in the centre of the spherical shell at ( a /2, b /2)

Figure 15: Benchmark seven, displacement components and stress components for a moderately thick spherical shell embedding one isotropic Steel layer ( R α / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w and σ α α in ( a /2, b /2); u and σ α z in (0, b /2)

The last and eighth benchmark shows an isotropic three-layered spherical panel having simply-supported sides (consult Fig. 1 to visualize the structure). The geometrical characteristics are contained in Table 2 and the properties for the three used materials are written in Table 3: Aluminium, Titanium and Steel (from bottom to top of the structure). The global thickness h is equally divided in three portions. Fig. 16 confirms that the temperature profile has to be calculated with 3D-u- θ or 3D( θ c , 3D) model for each possible thickness ratio because of the changing thermal properties of the three physical layers. For the thick structures cases, both the effects connected with the thickness layer and the embedded material are shown and the applications of the three-dimensional Fourier heat conduction relation is fundamental for exact results. These features are further verified by displacement and stress results in Table 14. The ( θ a ) profile is always inappropriate and the ( θ c , 1 D ) profile has several limitations for very thick structures. Only the 3D-u- θ and the 3D( θ c , 3D) models always guarantee satisfactory results for each R α / h , lamination stacking sequence, material layer and analyzed quantity. The results obtained from the thick structure ( R α / h = 10 ) in Fig. 17 highlight the classical zigzag behaviour of obtained variables in the thickness direction, connected with high transverse anisotropy. In addition, congruence conditions in terms of displacement components and equilibrium conditions in terms of transverse shear/normal stress components have been appropriately introduced in the proposed theory.

Figure 16: Benchmark eight, temperature profiles for two different R α / h ratios (thick and thin) of the spherical shell embedding three isotropic layers (Al2024/Ti22/Steel). Maximum value of θ ( α , β , z ) is investigated in the centre of the spherical shell at ( a /2, b /2)

Figure 17: Benchmark eight, displacement components and stress components for a moderately thick spherical shell embedding three isotropic layers (Al2024/Ti22/Steel) ( R α / h = 10 ) investigated via the 3D full coupled model (3D-u- θ ). Maximum values: w, σ α α and σ z z in ( a /2, b /2); v in ( a /2, 0)

4  Conclusions

An overall full coupled thermo-elastic 3D shell model in closed-form solution for the thermal stress investigation of single-and multi-layered isotropic, sandwich and composite plate and shell geometries has been discussed. The sovra-temperature maximum values have been directly applied at the highest and lowest point of outer faces in static conditions. The temperature profile is then evaluated through the thickness direction. The temperature profile is a primary variable as the displacements thanks to the coupling between the 3D Fourier heat conduction equation and the 3D equilibrium equations for shells. This sovra-temperature profile considers both the effects connected with the thickness layer and the embedded material for each possible case, without the necessity to separately solve the 3D or 1D form of the Fourier heat conduction equation. The global coupled system given by the three-dimensional equilibrium equations for shells and the three-dimensional Fourier heat conduction relation for shells is exactly solved via the exponential matrix methodology. Different results, given as displacement components, in-plane and out-of-plane stress components and sovra-temperature evaluations have been discussed for several thickness ratio values, geometrical properties, lamination sequences, temperature impositions and materials. These analyses showed a complete match between the 3D uncoupled model that uses the 3D Fourier heat conduction relation and the 3D full coupled one. This last new method has the advantage that takes into account both the thickness and the material layer effects using a mathematical formulation that is simpler and more elegant because the three-dimensional Fourier heat conduction relation is not separately solved. Moreover, a reduced number of artificial layers M is requested in comparison with the uncoupled 3D model.

Funding Statement: The authors received no specific funding for this study.

Conflict of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Librescu, L., Marzocca, P. (2003). Thermal stresses ’03, vol. 1. Blacksburg, VA, USA: Virginia Polytechnic Institute and State University. [Google Scholar]

2. Librescu, L., Marzocca, P. (2003). Thermal stresses ’03, vol. 2. Blacksburg, VA, USA: Virginia Polytechnic Institute and State University. [Google Scholar]

3. Noor, A. K., Burton, W. S. (1992). Computational models for high-temperature multilayered composite plates and shells. Applied Mechanics Reviews, 45, 419–446. https://doi.org/10.1115/1.3119742 [Google Scholar] [CrossRef]

4. Nowinski, J. L. (1978). Theory of thermoelasticity with applications, vol. 3. The Netherlands: Sijthoff & Noordhoff. [Google Scholar]

5. Altay, G. A., Dökmeci, M. C. (1996). Fundamental variational equations of discontinuous thermopiezoelectric fields. International Journal of Engineering Science, 34(7), 769–782. https://doi.org/10.1016/0020-7225(95)00133-6 [Google Scholar] [CrossRef]

6. Altay, G. A., Dökmeci, M. C. (1996). Some variational principles for linear coupled thermoelasticity. International Journal of Solids and Structures, 33(26), 3937–3948. https://doi.org/10.1016/0020-7683(95)00215-4 [Google Scholar] [CrossRef]

7. Altay, G. A., Dökmeci, M. C. (2001). Coupled thermoelastic shell equations with second sound for high-frequency vibrations of temperature-dependent materials. International Journal of Solids and Structures, 38(16), 2737–2768. https://doi.org/10.1016/S0020-7683(00)00179-7 [Google Scholar] [CrossRef]

8. Cannarozzi, A. A., Ubertini, F. (2001). A mixed variational method for linear coupled thermoelastic analysis. International Journal of Solids and Structures, 38(4), 717–739. https://doi.org/10.1016/S0020-7683(00)00061-5 [Google Scholar] [CrossRef]

9. Das, N. C., Das, S. N., Das, B. (1983). Eigenvalue approach to thermoelasticity. Journal of Thermal Stresses, 6(1), 35–43. https://doi.org/10.1080/01495738308942164 [Google Scholar] [CrossRef]

10. Kosínski, W., Frischmuth, K. (2001). Thermomechanical coupled waves in a nonlinear medium. Wave Motion, 34(1), 131–141. https://doi.org/10.1016/S0165-2125(01)00064-6 [Google Scholar] [CrossRef]

11. Wauer, J. (1996). Free and forced magneto-thermo-elastic vibrations in a conducting plate layer. Journal of Thermal Stresses, 19(7), 671–691. https://doi.org/10.1080/01495739608946201 [Google Scholar] [CrossRef]

12. Bhaskar, K., Varadan, T. K., Ali, J. S. M. (1996). Thermoelastic solution for orthotropic and anisotropic composite laminates. Composites. Part B: Engineering, 27(5), 415–420. https://doi.org/10.1016/1359-8368(96)00005-4 [Google Scholar] [CrossRef]

13. Pagano, N. J. (1969). Exact solutions for composite laminates in cylindrical bending. Journal of Composite Materials, 3(3), 398–411. https://doi.org/10.1177/002199836900300304 [Google Scholar] [CrossRef]

14. Pagano, N. J. (1970). Exact solutions for rectangular bidirectional composites and sandwich plates. Journal of Composite Materials, 4(1), 20–34. https://doi.org/10.1177/002199837000400102 [Google Scholar] [CrossRef]

15. Pagano, N. J., Wang, A. S. D. (1971). Further study of composite laminates under cylindrical bending. Journal of Composite Materials, 5(4), 521–528. https://doi.org/10.1177/002199837100500410 [Google Scholar] [CrossRef]

16. Kulikov, G. M., Plotnikova, S. V. (2015). A sampling surfaces method and its implementation for 3D thermal stress analysis of functionally graded plates. Composite Structures, 120, 315–325. https://doi.org/10.1016/j.compstruct.2014.10.012 [Google Scholar] [CrossRef]

17. Kulikov, G. M., Plotnikova, S. V. (2015). Three-dimensional thermal stress analysis of laminated composite plates with general layups by a sampling surfaces method. European Journal of Mechanics A/Solids, 49, 214–226. https://doi.org/10.1016/j.euromechsol.2014.07.011 [Google Scholar] [CrossRef]

18. Kulikov, G. M., Plotnikova, S. V. (2014). 3D exact thermoelastic analysis of laminated composite shells via sampling surfaces method. Composite Structures, 115, 120–130. https://doi.org/10.1016/j.compstruct.2014.04.019 [Google Scholar] [CrossRef]

19. Savoia, M., Reddy, J. N. (1995). Three-dimensional thermal analysis of laminated composite plates. International Journal of Solids and Structures, 32(5), 593–608. https://doi.org/10.1016/0020-7683(94)00146-N [Google Scholar] [CrossRef]

20. Tungikar, V. B., Rao, K. M. (1995). Three dimensional exact solution of thermal stresses in rectangular composite laminate. Composite Structures, 27(4), 419–430. https://doi.org/10.1016/0263-8223(94)90268-2 [Google Scholar] [CrossRef]

21. Brischetto, S., Torre, R. (2018). Thermo-elastic analysis of multilayered plates and shells based on 1D and 3D heat conduction problems. Composite Structures, 206, 326–353. https://doi.org/10.1016/j.compstruct.2018.08.042 [Google Scholar] [CrossRef]

22. Brischetto, S., Torre, R. (2019). 3D shell model for the thermo-mechanical analysis of FGM structures via imposed and calculated temperature profiles. Aerospace Science and Technology, 85, 125–149. https://doi.org/10.1016/j.ast.2018.12.011 [Google Scholar] [CrossRef]

23. Monge, J. C., Mantari, J. L. (2020). Exact solution of thermo-mechanical analysis of laminated composite and sandwich doubly-curved shell. Composite Structures, 245, 1–14. https://doi.org/10.1016/j.compstruct.2020.112323 [Google Scholar] [CrossRef]

24. Tornabene, F., Brischetto, S. (2018). 3D capability of refined GDQ models for the bending analysis of composite and sandwich plates, spherical and doubly-curved shells. Thin-Walled Structures, 129, 94–124. https://doi.org/10.1016/j.tws.2018.03.021 [Google Scholar] [CrossRef]

25. Brischetto, S., Tornabene, F. (2018). Advanced GDQ models and 3D stress recovery in multilayered plates, spherical and double-curved panels subjected to transverse shear loads. Composites Part B: Engineering, 146, 244–269. https://doi.org/10.1016/j.compositesb.2018.04.019 [Google Scholar] [CrossRef]

26. Alibeigloo, A., Liew, K. M. (2013). Thermoelastic analysis of functionally graded carbon nanotube-reinforced composite plate using theory of elasticity. Composite Structures, 106, 873–881. https://doi.org/10.1016/j.compstruct.2013.07.002 [Google Scholar] [CrossRef]

27. Moleiro, F., Correia, V. M. F., Ferreira, A. J. M., Reddy, J. N. (2019). Fully coupled thermo-mechanical analysis of multilayered plates with embedded FGM skins or core layers using a layerwise mixed model. Composite Structures, 210, 971–996. https://doi.org/10.1016/j.compstruct.2018.11.073 [Google Scholar] [CrossRef]

28. Adineh, M., Kadkhodayan, M. (2017). Three-dimensional thermo-elastic analysis and dynamic response of a multi-directional functionally graded skew plate on elastic foundation. Composites Part B: Engineering, 125, 227–240. https://doi.org/10.1016/j.compositesb.2017.05.070 [Google Scholar] [CrossRef]

29. Rolfes, R., Noack, J., Taeschner, M. (1999). High performance 3D-analysis of thermo-mechanically loaded composite structures. Composite Structures, 46(4), 367–379. https://doi.org/10.1016/S0263-8223(99)00101-4 [Google Scholar] [CrossRef]

30. Thangaratnam, R. K., Palaninathan, R., Ramachandran, J. (1988). Thermal stress analysis of laminated composite plates and shells. Computers & Structures, 30(6), 1403–1411. https://doi.org/10.1016/0045-7949(88)90204-0 [Google Scholar] [CrossRef]

31. Dehghan, M., Nejad, M. Z., Moosaie, A. (2016). Thermo-electro-elastic analysis of functionally graded piezoelectric shells of revolution: Governing equations and solutions for some simple cases. International Journal of Engineering Science, 104, 34–61. [Google Scholar]

32. Padovan, J. (1976). Thermoelasticity of cylindrically anisotropic generally laminated cylinders. Journal of Applied Mechanics, 43, 124–130. [Google Scholar]

33. Trajkovski, D., Cukic, R. (1999). A coupled problem of thermoelastic vibrations of a circular plate with exact boundary conditions. Mechanics Research Communications, 26(2), 217–224. [Google Scholar]

34. Yeh, Y. L. (2005). The effect of thermo-mechanical coupling for a simply supported orthotropic rectangular plate on non-linear dynamics. Thin-Walled Structures, 43(8), 1277–1295. [Google Scholar]

35. Vinyas, M., Kattimani, S. C. (2017). Hygrothermal analysis of magneto-electro-elastic plate using 3D finite element analysis. Composite Structures, 180, 617–637. [Google Scholar]

36. Pantano, A., Averill, R. C. (2000). A 3D zig-zag sublaminate model for analysis of thermal stresses in laminated composite and sandwich plates. Journal of Sandwich Structures & Materials, 2(3), 288–312. [Google Scholar]

37. Khare, R. K., Kant, T., Garg, A. K. (2003). Closed-form thermo-mechanical solutions of higher-order theories of cross-ply laminated shallow shells. Composite Structures, 59(3), 313–340. [Google Scholar]

38. Khdeir, A. A., Reddy, J. N. (1991). Thermal stresses and deflections of cross-ply laminated plates using refined plate theories. Journal of Thermal Stresses, 14(4), 419–438. [Google Scholar]

39. Khdeir, A. A., Reddy, J. N. (1991). Thermal effects on the response of cross-ply laminated shallow shells. International Journal of Solids and Structures, 29(5), 653–667. [Google Scholar]

40. Kalam, M. A., Tauchert, T. R. (1978). Stresses in an orthotropic elastic cylinder due to a plane temperature distribution T(r,θ). Journal of Thermal Stresses, 1(1), 13–24. [Google Scholar]

41. He, J. F. (1995). Thermoelastic analysis of laminated plates including transverse shear deformation effects. Composite Structures, 30(1), 51–59. https://doi.org/10.1016/0263-8223(94)00026-3 [Google Scholar] [CrossRef]

42. Ali, J. S. M., Bhaskar, K., Varadan, T. K. (1999). A new theory for accurate thermal/mechanical flexural analysis of symmetric laminated plates. Composite Structures, 45(3), 227–232. https://doi.org/10.1016/S0263-8223(99)00028-8 [Google Scholar] [CrossRef]

43. Brischetto, S. (2009). Effect of the through-the-thickness temperature distribution on the response of layered and composite shells. International Journal of Applied Mechanics, 1(4), 581–605. https://doi.org/10.1142/S1758825109000393 [Google Scholar] [CrossRef]

44. Brischetto, S., Carrera, E. (2010). Coupled thermo-mechanical analysis of one-layered and multilayered isotropic and composite shells. Computer Modeling in Engineering & Sciences, 56(3), 249–302. https://doi.org/10.3970/cmes.2010.056.249 [Google Scholar] [CrossRef]

45. Brischetto, S., Carrera, E. (2010). Coupled thermo-mechanical analysis of one-layered and multilayered plates. Composite Structures, 92(8), 1793–1812. https://doi.org/10.1016/j.compstruct.2010.01.020 [Google Scholar] [CrossRef]

46. Miller, C. J., Millavec, W. A., Kicher, T. P. (1981). Thermal stress analysis of layered cylindrical shells. AIAA Journal, 19(4), 523–530. https://doi.org/10.2514/3.7790 [Google Scholar] [CrossRef]

47. Cheng, Z. Q., Batra, R. C. (2001). Thermal effects on laminated composite shells containing interfacial imperfections. Composite Structures, 52(1), 3–11. https://doi.org/10.1016/S0263-8223(00)00197-5 [Google Scholar] [CrossRef]

48. Barut, A., Madenci, E., Tessler, A. (2000). Nonlinear thermoelastic analysis of composite panels under non-uniform temperature distribution. International Journal of Solids and Structures, 37(27), 3681–3713. https://doi.org/10.1016/S0020-7683(99)00119-5 [Google Scholar] [CrossRef]

49. Huang, N. N., Tauchert, T. R. (1991). Large deflections of laminated cylindrical and doubly-curved panels under thermal loading. Computers & Structures, 41(2), 303–312. https://doi.org/10.1016/0045-7949(91)90433-M [Google Scholar] [CrossRef]

50. Noor, A. K., Peters, J. M. (1999). Analysis of curved sandwich panels with cutouts subjected to combined temperature gradient and mechanical loads. Journal of Sandwich Structures and Materials, 1(1), 42–59. https://doi.org/10.1177/109963629900100103 [Google Scholar] [CrossRef]

51. Raju, K. K., Rao, G. V. (1984). Finite element analysis of thermal postbuckling of tapered columns. Computer & Structures, 19(4), 617–620. https://doi.org/10.1016/0045-7949(84)90108-1 [Google Scholar] [CrossRef]

52. Daneshjoo, K., Ramezani, M. (2002). Coupled thermoelasticity in laminated composite plates based on green-lindsay model. Composite Structures, 55(4), 387–392. https://doi.org/10.1016/S0263-8223(01)00164-7 [Google Scholar] [CrossRef]

53. Daneshjoo, K., Ramezani, M. (2004). Classical coupled thermoelasticity in laminated composite plates based on third-order shear deformation theory. Composite Structures, 64(3–4), 369–375. https://doi.org/10.1016/j.compstruct.2003.09.039 [Google Scholar] [CrossRef]

54. Ibrahimbegovic, A., Colliat, J. B., Davenne, L. (2005). Thermomechanical coupling in folded plates and non-smooth shells. Computer Methods in Applied Mechanics Engineering, 194(21–24), 2686–2707. https://doi.org/10.1016/j.cma.2004.07.052 [Google Scholar] [CrossRef]

55. Lee, Z. Y. (2006). Generalized coupled transient thermoelastic problem of multilayered hollow cylinder with hybrid boundary conditions. International Communications in Heat and Mass Transfer, 33(4), 518–528. https://doi.org/10.1016/j.icheatmasstransfer.2006.01.001 [Google Scholar] [CrossRef]

56. Oh, J., Cho, M. (2004). A finite element based on cubic zig-zag plate theory for the prediction of thermo-electric-mechanical behaviors. International Journal of Solids and Structures, 41(5–6), 1357–1375. https://doi.org/10.1016/j.ijsolstr.2003.10.019 [Google Scholar] [CrossRef]

57. Brischetto, S. (2014). An exact 3D solution for free vibrations of multilayered cross-ply composite and sandwich plates and shells. International Journal of Applied Mechanics, 6(6), 1–42. https://doi.org/10.1142/S1758825114500768 [Google Scholar] [CrossRef]

58. Brischetto, S. (2017). Exact three-dimensional static analysis of single-and multi-layered plates and shells. Composites Part B: Engineering, 119, 230–252. https://doi.org/10.1016/j.compositesb.2017.03.010 [Google Scholar] [CrossRef]

59. Brischetto, S. (2017). A closed-form 3D shell solution for multilayered structures subjected to different load combinations. Aerospace Science and Technology, 70, 29–46. https://doi.org/10.1016/j.ast.2017.07.040 [Google Scholar] [CrossRef]

60. Özişik, M. N. (1993). Heat conduction. New York: John Wiley & Sons, Inc. [Google Scholar]

61. Povstenko, Y. (2015). Fractional thermoelasticity. Switzerland: Springer International Publishing. [Google Scholar]

62. Moon, P., Spencer, D. E. (1988). Field theory handbook: Including coordinate systems, differential equations and their solutions. Berlin: Springer-Verlag. [Google Scholar]

63. Mikhailov, M. D., Özişik, M. N. (1984). Unified analysis and solutions of heat and mass diffusion. New York: Dover Publications Inc. [Google Scholar]

64. Arfken, G. B., Weber, H. J. (2005). Mathematical methods for physicists, Sixth edition. San Diego, USA: Elsevier Academic Press. [Google Scholar]

65. Morse, P. M., Feshbach, H. (1953). Methods of theoretical physics. USA: McGraw-Hill. [Google Scholar]

66. Boyce, W. E., Prima, R. C. D. (2001). Elementary differential equations and boundary value problems. New York: John Wiley & Sons, Ltd. [Google Scholar]

67. Systems of differential equations (2013). http://www.math.utah.edu/gustafso/ [Google Scholar]