A Three-Dimensional Layer-Wise Formulation for the Coupled Thermo-Magneto-Elastic Analysis of Multilayered Composite Flat and Curved Panels

1. Introduction

Aerospace structures operate under service conditions where several physical fields interact simultaneously. The mutual interaction among these fields may significantly influence the structural response and generate deformation mechanisms that cannot be captured through uncoupled analyses [1,2]. Among the various multifield phenomena, thermo-magneto-elastic coupling represents a particularly relevant topic for advanced aerospace applications, including spacecraft, satellites, and launch vehicles. Although the scientific and industrial interest in this research area has steadily increased in recent years, only a limited number of contributions have addressed the fully coupled thermo-magneto-elastic problem. Experimental investigations are also of considerable importance, since they provide a deeper understanding of the interactions among the thermal, magnetic, and elastic fields under prescribed loading conditions. Such analyses make it possible to distinguish the contribution of each physical field to the overall coupled response. In addition, the direct piezomagnetic effect offers the possibility of designing actuators capable of producing strains and displacements under externally applied magnetic fields, whereas the inverse piezomagnetic effect provides the basis for sensing devices able to detect deformations and displacements generated by combined thermal and mechanical actions.

Thermo-magneto-elastic coupling has been investigated in the literature for both flat and curved structural configurations. Concerning flat structures, including beams and plates, several contributions have been reported. Guo and Alam [3] investigated the nonlinear bending response and thermal post-buckling behaviour of magneto-electro-elastic (MEE) beams by incorporating surface molecular effects within a nonlocal strain-gradient continuum framework. A node-based smoothed radial point interpolation formulation was proposed by Ren et al. [4] to investigate the free vibration characteristics and thermo-mechanical behaviour of functionally graded MEE structures. The same research group subsequently extended the formulation to the transient thermo-mechanical analysis of functionally graded MEE components containing cut-outs and operating under thermal loading conditions [5]. The influence of size-dependent effects on the buckling and free-vibration responses of functionally graded magneto-electro-thermo-elastic Timoshenko nanobeams was later examined by Ren and Qing [6]. A fully three-dimensional finite element formulation for multilayered MEE beams subjected to thermal loading was introduced by Vinyas and Kattimani [7]. Zhang et al. [8] employed a higher-order shear deformation beam model combined with von Kármán geometric nonlinearity to investigate the nonlinear bending and vibration behaviour of laminated MEE beams under thermal environments. In a subsequent study, the same authors [9] adopted the generalized differential quadrature method together with the Timoshenko beam theory and von Kármán nonlinear kinematics to analyze the thermal post-buckling response of laminated MEE beams. Thuy Anh et al. [10] investigated the nonlocal static behaviour of magneto-electro-elastic sandwich micro- and nano-plates reinforced with functionally graded carbon nanotubes and operating under hygrothermal conditions. Esen and Ozmen [11] analyzed the thermal vibration and buckling behaviour of porous nanoplates made of functionally graded barium titanate and cobalt ferrite materials. The coupled thermo-magneto-electro-elastic response of heterogeneous circular plates subjected to uniformly distributed thermal loads along the boundary was examined by Li et al. [12]. In a subsequent contribution, the same authors [13] developed a coupled cell-based smoothed finite element formulation to study the dynamic behaviour of functionally graded MEE thin-walled structures subjected to combined mechanical loading and thermal environments. Mahesh [14,15] proposed a finite element framework integrated with an artificial neural network algorithm to predict both the static response and the nonlinear pyro-coupled deflection of smart sandwich plates composed of a porous agglomerated carbon nanotube nanocomposite core and piezo-magneto-thermo-electric face sheets. Mohammadimehr et al. [16] carried out bending, buckling, and free-vibration analyses of microcomposite circular-annular sandwich plates subjected to hydro-thermo-magneto-mechanical loading conditions. AlMukahal et al. [17] proposed a numerical procedure for the axisymmetric bending analysis of functionally graded graphene platelet-reinforced annular sandwich nanoplates with porous cores. Talleb and Ren [18] introduced a multiscale analytical magneto-elastic formulation to characterize magnetostrictive structures under coupled magnetic, thermal, and mechanical loadings. Vinyas and Kattimani [19,20] presented a three-dimensional finite element formulation to investigate both the thermo-multiphysics response of multilayered MEE plates and the static behaviour of MEE plates subjected to hygrothermal loading conditions. Vu and Tran [21] investigated the thermo-electro-magneto-mechanical free-vibration response of sandwich plates with a graphene oxide powder-reinforced nanocomposite core and functionally graded MEE outer layers. Wang et al. [22] analyzed the coupled magneto-thermo-elastic behaviour of simply supported soft ferromagnetic plates subjected to magnetic and thermal loading by adopting the linearized theory of magneto-elasticity together with a perturbation approach. Wu and Hsu [23] formulated a unified size-dependent shear deformation theory to investigate the free-vibration characteristics of functionally graded magneto-electro-elastic microplates under simply supported boundary conditions. Xing and Liu [24] employed the differential quadrature method to perform dynamic and quasi-static analyses of rectangular magneto-thermo-elastic plates under arbitrarily varying magnetic fields. Zhou et al. [25] proposed an inhomogeneous magneto-electro-elastic coupling element-free Galerkin formulation for static analyses. Saadatfar et al. [26] investigated the creep response of variable-thickness annular functionally graded MEE plates. Zenkour [27] examined the magneto-thermo-elastic behaviour of a functionally graded annular sandwich disk placed in a magnetic field and subjected to non-uniform steady-state thermal loading. More recently, Brischetto et al. [28,29] introduced a three-dimensional analytical coupled formulation for the static and free-vibration analysis of magneto-electro-elastic plates, where thermal effects and shell curvature were intentionally neglected.

The number of studies devoted to thermo-magneto-elastic coupling in curved structures, particularly shell configurations, is significantly smaller than that available for beams and plates. Among these contributions, Pham et al. [30] analyzed the free-vibration behaviour of a functionally graded porous curved nanobeam resting on an elastic foundation and operating under coupled hygro-thermo-magnetic conditions. An analytical formulation for anisotropic multilayer magneto-piezo-elastic hollow spheres and heterogeneous piezo-magneto-elastic discs subjected to coupled magneto-electro-hygro-thermo-mechanical loading was presented by Barati and Shariyat [31,32]. Bian and Tian [33] examined the thermo-magneto-elastic response of thin conductive spherical segment shells subjected to externally applied magnetic fields. Dai and Wang [34] derived an analytical solution describing the magneto-thermo-electro-elastic behaviour of a piezoelectric hollow cylinder immersed in an axial magnetic field and subjected to arbitrary thermal shocks, transient electrical excitations, and mechanical loading. Gui and Wu [35] analyzed the buckling behaviour of thermo-magneto-electro-elastic nanocylindrical shells embedded in an elastic medium. Higuchi et al. [36,37] investigated the dynamic and quasi-static magneto-thermo-elastic stress fields generated by transient magnetic excitations in both solid and hollow conducting circular cylinders. Hu et al. [38,39] studied the coupled magneto-thermo-elastic vibration response of rotating and stationary ferromagnetic functionally graded cylindrical shells. Kong et al. [40] proposed an analytical framework for evaluating thermo-magneto-elastic stresses together with the perturbation of the magnetic field vector in non-homogeneous hollow cylinders subjected to thermal shock loading. The magneto-thermo-elastic behaviour of a thick double-walled cylinder composed of a functionally graded interlayer and a homogeneous outer layer was investigated by Loghman and Parsa [41]. Ni et al. [42] established a buckling formulation for magneto-electro-elastic composite cylindrical shells subjected to coupled hygro-thermo-magneto-electro-elastic loading conditions. Oveissi et al. [43] investigated hydro-magneto-electro wave propagation in an axially moving circular cylindrical nanoshell resting on an electromagnetic-visco-Pasternak foundation and conveying a magnetic nanofluid under different thermal and hygrothermal loading conditions. Selvamani et al. [44] adopted a multiple-scale perturbation technique to investigate the nonlinear wave propagation characteristics of doubly curved sandwich composite piezoelectric shells with flexible cores operating in hygrothermal environments. Tong et al. [45] introduced a three-phase cylindrical model for the analysis of fiber-reinforced composite structures subjected to in-plane mechanical loading under coupled thermal, electric, magnetic, and elastic fields. Tornabene et al. [46,47,48] proposed refined two-dimensional multifield formulations based on a generalized equivalent-layer-wise approach for the hygro-thermo-magneto-electro-elastic analysis of curved shell structures. Finally, Lang and Xuewu [49] presented a formulation capable of predicting both the buckling and free-vibration responses of functionally graded magneto-electro-thermo-elastic circular cylindrical shells. A series of three-dimensional analytical formulations for the static and free-vibration analysis of spherical shells was proposed by Brischetto et al. [50,51,52,53], addressing thermo-elastic, magneto-elastic, and magneto-electro-elastic coupling phenomena. However, none of these contributions considered the fully coupled thermo-magneto-elastic problem. In addition, Ponnusamy and Selvamani [54] analyzed the three-dimensional propagation of elastic waves in a homogeneous transversely isotropic magneto-thermo-elastic cylindrical panel.

The objective of the present work is the development of a three-dimensional thermo-magneto-elastic shell formulation suitable for both flat and curved multilayered panels. Curved geometries are described by adopting a mixed orthogonal curvilinear reference system ($\alpha$, $\beta$, z), which provides a unified mathematical framework for the analysis of plates, cylinders, cylindrical panels, and spherical shells through the appropriate definition of the curvature radii $R_{\alpha }$ and $R_{\beta }$. One of the main features of the proposed formulation is therefore its ability to analyze different structural configurations within a unique thermo-magneto-elastic framework. The model is governed by five coupled second-order differential equations expressed in terms of the displacement components, magnetic potential, and over-temperature. These equations are obtained by combining the three-dimensional elastic equilibrium equations with the magnetic induction divergence equation and the three-dimensional Fourier heat conduction equation. The resulting system is solved analytically by adopting Navier harmonic expansions in the in-plane directions together with the exponential matrix method along the thickness coordinate. Such a solution procedure enables an exact three-dimensional evaluation of the secondary field variables, including strains, stresses, magnetic induction components, and heat fluxes throughout the thickness. Finally, the multilayered configuration is represented through a layer-wise description in which interlaminar continuity conditions are enforced at every interface between adjacent layers.

The manuscript is organized as follows. Section 2 describes the theoretical framework of the proposed formulation, including the analytical derivation of the three-dimensional multifield governing, constitutive, and geometrical equations. Section 3 first presents a preliminary assessment to validate the proposed 3D thermo-magneto-elastic model for both flat and curved structures. It then introduces a set of new benchmark cases covering different thickness ratios, for which both tabulated results at selected through-the-thickness positions and graphical distributions of the main field variables are reported. These benchmarks illustrate the capability of the formulation to capture the effects associated with thermo-magneto-elastic coupling, material lay-up, thickness ratio, zigzag behaviour, and shell curvature. Finally, the main conclusions and the principal outcomes of the present study are summarized in Section 4.

2. 3D Thermo-Magneto-Elastic Static Analysis of Shells

This section presents the theoretical formulation adopted for the three-dimensional thermo-magneto-elastic static analysis of multilayered shell structures. The presentation is divided into four subsections. The first one introduces the three-dimensional multifield governing equations expressed in the mixed orthogonal curvilinear reference system. The constitutive relations and the multifield congruence equations are subsequently presented in the second and third subsections, respectively. The final subsection describes the analytical solution procedure based on the exponential matrix method.

2.1. 3D Multifield Governing Equations

The thermo-magneto-elastic static problem is governed by a system of five fully coupled three-dimensional equations. These consist of the three equilibrium equations for the elastic field, the divergence equation for the magnetic induction field, and the steady-state Fourier heat conduction equation. The complete system is formulated in the mixed orthogonal curvilinear reference system ($\alpha$, $\beta$, z) illustrated in Figure 1. This coordinate system combines two curvilinear in-plane coordinates ($\alpha$ and $\beta$) with a rectilinear thickness coordinate (z), thus providing a unified framework for the analysis of different structural configurations. By assigning the appropriate values to the curvature radii $R_{\alpha }$ and $R_{\beta }$, the same formulation can be directly applied to plates, cylinders, cylindrical panels, and spherical shells.

The thermo-magneto-elastic governing equations for spherical shells are explicitly expressed as

$$H_{\beta }{\sigma }_{\alpha {\alpha }_{,\alpha }}^k+H_{\alpha }{\sigma }_{\alpha {\beta }_{,\beta }}^k+H_{\alpha }{H}_{\beta }{\sigma }_{\alpha z_{,z}}^k+\left({\displaystyle \frac{2H_{\beta }}{R_{\alpha }}}+{\displaystyle \frac{H_{\alpha }}{R_{\beta }}}\right){\sigma }_{\alpha z}^k=0,$$

(1)

$$H_{\beta }{\sigma }_{\alpha {\beta }_{,\alpha }}^k+H_{\alpha }{\sigma }_{\beta {\beta }_{,\beta }}^k+H_{\alpha }{H}_{\beta }{\sigma }_{\beta z_{,z}}^k+\left({\displaystyle \frac{2H_{\alpha }}{R_{\beta }}}+{\displaystyle \frac{H_{\beta }}{R_{\alpha }}}\right){\sigma }_{\beta z}^k=0,$$

(2)

$$H_{\beta }{\sigma }_{\alpha z_{,\alpha }}^k+H_{\alpha }{\sigma }_{\beta z_{,\beta }}^k+H_{\alpha }{H}_{\beta }{\sigma }_{z{z}_{,z}}^k+{\displaystyle \frac{H_{\beta }}{R_{\alpha }}}{\sigma }_{\alpha \alpha }^k+{\displaystyle \frac{H_{\alpha }}{R_{\beta }}}{\sigma }_{\beta \beta }^k+\left({\displaystyle \frac{H_{\beta }}{R_{\alpha }}}+{\displaystyle \frac{H_{\alpha }}{R_{\beta }}}\right){\sigma }_{zz}^k=0,$$

(3)

$${\displaystyle \frac{{\mathcal{B}}_{{\alpha }_{,\alpha }}^k}{H_{\alpha }}}+{\displaystyle \frac{{\mathcal{B}}_{{\beta }_{,\beta }}^k}{H_{\beta }}}+{\mathcal{B}}_{z_{,z}}^k=0,$$

(4)

$${\displaystyle \frac{q_{{\alpha }_{,\alpha }}^k}{H_{\alpha }}}+{\displaystyle \frac{q_{{\beta }_{,\beta }}^k}{H_{\beta }}}+q_{z_{,z}}^k=0.$$

(5)

where ${\sigma }_{\alpha \alpha }^k$, ${\sigma }_{\alpha \beta }^k$, ${\sigma }_{\alpha z}^k$, ${\sigma }_{\beta \beta }^k$, ${\sigma }_{\beta z}^k$, and ${\sigma }_{zz}^k$ denote the six stress components, whereas ${\mathcal{B}}_{\alpha }^k$, ${\mathcal{B}}_{\beta }^k$, and ${\mathcal{B}}_z^k$ represent the magnetic induction components. Similarly, $q_{\alpha }^k$, $q_{\beta }^k$, and $q_z^k$ are the heat flux components. The subscripts $,\alpha$, $,\beta$, and $,z$ indicate partial differentiation with respect to the corresponding coordinates, while the superscript k identifies the physical layer. The quantities $R_{\alpha }$ and $R_{\beta }$ denote the radii of curvature in the two in-plane directions. The formulation reported in Eqs. (1)–(5) follows the theoretical developments presented in [50,55,56]. The curvature parameters appearing in the governing equations are defined as

$$\left\{\begin{array}{c}{H}_{\alpha }\\ H_{\beta }\\ H_z\end{array}\right\}=\left\{\begin{array}{c}{H}_{\alpha }\left(z\right)\\ H_{\beta }\left(z\right)\\ H_z\left(z\right)\end{array}\right\}=\left\{\begin{array}{c}1+{\displaystyle \frac{z}{R_{\alpha }}}\\ 1+{\displaystyle \frac{z}{R_{\beta }}}\\ 1\end{array}\right\}=\left\{\begin{array}{c}1+{\displaystyle \frac{\tilde{z}-h/2}{R_{\alpha }}}\\ 1+{\displaystyle \frac{\tilde{z}-h/2}{R_{\beta }}}\\ 1\end{array}\right\},-h/2<z<h/2\mathrm{or}0<\tilde{z}<h.$$

(6)

h denotes the total thickness of the structure. The curvature coefficients $H_{\alpha }$ and $H_{\beta }$ differ from unity because the $\alpha$ and $\beta$ coordinates follow curved directions, whereas $H_z$ is equal to one since the thickness coordinate z is rectilinear. For the sake of conciseness, the dependence of the curvature coefficients $H_{\alpha }$ and $H_{\beta }$ on the thickness coordinate is omitted in the remainder of the paper, although it is implicitly retained throughout the formulation.

2.2. 3D Multifield Constitutive Equations

The constitutive relations governing the coupled thermo-magneto-elastic problem are derived from the definition of a multifield free energy density function. This function is expressed as

$$\mathcal{W}({\epsilon }_{ij},{\mathcal{H}}_i,T,\theta )$$

where ${\epsilon }_{ij}$ denotes the strain tensor components, ${\mathcal{H}}_i$ represents the magnetic field components, T is the absolute temperature, and $\theta$ is the over-temperature, defined as $\theta =T-T_0$, with $T_0$ being the reference environmental temperature expressed in Kelvin. The adopted free energy density accounts for the elastic, magnetic, and thermal contributions together with all the associated coupling effects, namely the magneto-elastic, thermo-elastic, and thermo-magnetic interactions. Neglecting the entropy contribution, the free energy density assumes the following quadratic form:

$$\mathcal{W}({\epsilon }_{ij},{\mathcal{H}}_i,\theta )={\displaystyle \frac{1}{2}}{\epsilon }_{ij}{C}_{ijkl}{\epsilon }_{kl}-{\epsilon }_{ij}{r}_{ijk}{\mathcal{H}}_k-{\epsilon }_{ij}{\lambda }_{ij}\theta +{\displaystyle \frac{1}{2}}{H}_i{\mu }_{ij}{H}_j+H_i{g}_i\theta +{\displaystyle \frac{1}{2}}{\theta }^T{k}_{ij}\theta .$$

(7)

Here, $C_{ijkl}$ denotes the elastic stiffness coefficients, $r_{ijk}$ the piezomagnetic coefficients, ${\lambda }_{ij}$ the thermo-elastic coefficients, ${\mu }_{ij}$ the magnetic permeability coefficients, $g_i$ the thermo-magnetic coupling coefficients, and $k_{ij}$ the thermal conductivity coefficients. The constitutive equations are then obtained by differentiating the free energy density with respect to the corresponding field variables:

$$\sigma ={\displaystyle \frac{\partial \mathcal{W}}{\partial \epsilon }},$$

(8)

$$\mathcal{B}={\displaystyle \frac{\partial \mathcal{W}}{\partial \mathbf{H}}},$$

(9)

$$\mathbf{q}=-{\displaystyle \frac{\partial \mathcal{W}}{\partial \theta }}.$$

(10)

By evaluating Eqs. (8)–(10), the explicit form of the 3D thermo-magneto-elastic constitutive equations is obtained:

$${\sigma }_{\alpha \alpha }^k=C_{11}^k{\epsilon }_{\alpha \alpha }^k+C_{12}^k{\epsilon }_{\beta \beta }^k+C_{13}^k{\epsilon }_{zz}^k-r_{31}^k{\mathcal{H}}_z^k-{\lambda }_{11}^k{\theta }^k,$$

(11)

$${\sigma }_{\beta \beta }^k=C_{12}^k{\epsilon }_{\alpha \alpha }^k+C_{22}^k{\epsilon }_{\beta \beta }^k+C_{23}^k{\epsilon }_{zz}^k-r_{32}^k{\mathcal{H}}_z^k-{\lambda }_{22}^k{\theta }^k,$$

(12)

$${\sigma }_{zz}^k=C_{13}^k{\epsilon }_{\alpha \alpha }^k+C_{23}^k{\epsilon }_{\beta \beta }^k+C_{33}^k{\epsilon }_{zz}^k-r_{33}^k{\mathcal{H}}_z^k-{\lambda }_{33}^k{\theta }^k,$$

(13)

$${\sigma }_{\beta z}^k=C_{44}^k{\gamma }_{\beta z}^k-r_{24}^k{\mathcal{H}}_{\beta }^k,$$

(14)

$${\sigma }_{\alpha z}^k=C_{55}^k{\gamma }_{\alpha z}^k-r_{15}^k{\mathcal{H}}_{\alpha }^k,$$

(15)

$${\sigma }_{\alpha \beta }^k=C_{66}^k{\gamma }_{\alpha \beta }^k,$$

(16)

$${\mathcal{B}}_{\alpha }^k=r_{15}^k{\gamma }_{\alpha z}^k+{\mu }_{11}^k{\mathcal{H}}_{\alpha }^k,$$

(17)

$${\mathcal{B}}_{\beta }^k=r_{24}^k{\gamma }_{\beta z}^k+{\mu }_{22}^k{\mathcal{H}}_{\beta }^k,$$

(18)

$${\mathcal{B}}_z^k=r_{31}^k{\epsilon }_{\alpha \alpha }^k+r_{32}^k{\epsilon }_{\beta \beta }^k+r_{33}^k{\epsilon }_{zz}^k+{\mu }_{33}^k{\mathcal{H}}_z^k+g_3^k{\theta }^k,$$

(19)

$$q_{\alpha }^k=-{\displaystyle \frac{k_{11}^k}{H_{\alpha }}}{\theta }_{,\alpha }^k,$$

(20)

$$q_{\beta }^k=-{\displaystyle \frac{k_{22}^k}{H_{\beta }}}{\theta }_{,\beta }^k,$$

(21)

$$q_z^k=-k_{33}^k{\theta }_{,z}^k,$$

(22)

where ${\epsilon }_{\alpha \alpha }^k$, ${\epsilon }_{\beta \beta }^k$, ${\epsilon }_{zz}^k$, ${\gamma }_{\beta z}^k$, ${\gamma }_{\alpha z}^k$, and ${\gamma }_{\alpha \beta }^k$ denote the six strain components. The elastic behaviour of the cross-ply orthotropic lamina is described by the stiffness coefficients $C_{11}^k$, $C_{12}^k$, $C_{13}^k$, $C_{22}^k$, $C_{23}^k$, $C_{33}^k$, $C_{44}^k$, $C_{55}^k$, and $C_{66}^k$. The piezomagnetic coupling is characterized by the coefficients $r_{31}^k$, $r_{32}^k$, $r_{33}^k$, $r_{24}^k$, and $r_{15}^k$, while ${\mathcal{H}}_{\alpha }^k$, ${\mathcal{H}}_{\beta }^k$, and ${\mathcal{H}}_z^k$ represent the magnetic field components. The thermo-elastic coupling coefficients are denoted by ${\lambda }_{11}^k$, ${\lambda }_{22}^k$, and ${\lambda }_{33}^k$. The over-temperature is indicated by ${\theta }^k$, with $\theta =T-T_0$ and $T_0$ representing the reference environmental temperature expressed in Kelvin. The magnetic permeability coefficients are given by ${\mu }_{11}^k$, ${\mu }_{22}^k$, and ${\mu }_{33}^k$, whereas $g_3^k$ is the thermo-magnetic coupling coefficient relating the over-temperature to the transverse magnetic induction component. Finally, $k_{11}^k$, $k_{22}^k$, and $k_{33}^k$ are the thermal conductivity coefficients.

The constitutive relations reported in Eqs. (11)–(22) provide a fully coupled description of the mechanical, magnetic, and thermal behaviour of the multilayered material.

2.3. 3D Multifield Congruence Equations

The 3D multifield congruence equations establish the kinematic relationships of the thermo-magneto-elastic problem by linking the displacement field with the strain components and the magnetic potential with the magnetic field. Their explicit form is given by

$${\epsilon }_{\alpha \alpha }^k={\displaystyle \frac{u_{,\alpha }^k}{H_{\alpha }}}+{\displaystyle \frac{w^k}{H_{\alpha }{R}_{\alpha }}},$$

(23)

$${\epsilon }_{\beta \beta }^k={\displaystyle \frac{v_{,\beta }^k}{H_{\beta }}}+{\displaystyle \frac{w^k}{H_{\beta }{R}_{\beta }}},$$

(24)

$${\epsilon }_{zz}^k=w_{,z}^k,$$

(25)

$${\gamma }_{\alpha \beta }^k={\displaystyle \frac{v_{,\alpha }^k}{H_{\alpha }}}+{\displaystyle \frac{u_{,\beta }^k}{H_{\beta }}},$$

(26)

$${\gamma }_{\alpha z}^k={\displaystyle \frac{w_{,\alpha }^k}{H_{\alpha }}}+u_{,z}^k-{\displaystyle \frac{u^k}{H_{\alpha }{R}_{\alpha }}},$$

(27)

$${\gamma }_{\beta z}^k={\displaystyle \frac{w_{,\beta }^k}{H_{\beta }}}+v_{,z}^k-{\displaystyle \frac{v^k}{H_{\beta }{R}_{\beta }}}.$$

(28)

$${\mathcal{H}}_{\alpha }^k=-{\displaystyle \frac{{\psi }_{,\alpha }^k}{H_{\alpha }}},$$

(29)

$${\mathcal{H}}_{\beta }^k=-{\displaystyle \frac{{\psi }_{,\beta }^k}{H_{\beta }}},$$

(30)

$${\mathcal{H}}_z^k=-{\psi }_{,z}^k.$$

(31)

Here, ${\psi }^k$ denotes the magnetic potential associated with the k-th physical layer. Together, Eqs. (23)–(31) provide the kinematic description required to couple the mechanical and magnetic fields within the proposed three-dimensional formulation.

2.4. Analytical Solution via the Exponential Matrix Method

This subsection presents the analytical solution procedure via the exponential matrix method. By recursively introducing the 3D multifield congruence equations in Eqs. (23)–(31) into the constitutive equations in Eqs. (11)–(22), and subsequently substituting the resulting expressions into the governing equations in Eqs. (1)–(5), the complete system can be reformulated in terms of the primary variables u, v, w, $\psi$, and $\theta$. These variables respectively represent the three displacement components, the magnetic potential, and the over-temperature, and constitute the primary unknowns of the coupled thermo-magneto-elastic problem.

The analytical solution is obtained by expressing the primary variables through Navier harmonic expansions. Their explicit expressions are given by

$$u^k(\alpha ,\beta ,z)=U^k\left(z\right)cos\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),$$

(32)

$$v^k(\alpha ,\beta ,z)=V^k\left(z\right)sin\left(\overline{\alpha }\alpha \right)cos\left(\overline{\beta }\beta \right),$$

(33)

$$w^k(\alpha ,\beta ,z)=W^k\left(z\right)sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),$$

(34)

$${\psi }^k(\alpha ,\beta ,z)={\mathsf{\Psi }}^k\left(z\right)sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),$$

(35)

$${\theta }^k(\alpha ,\beta ,z)={\mathsf{\Theta }}^k\left(z\right)sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right).$$

(36)

The functions $U^k\left(z\right)$, $V^k\left(z\right)$, $W^k\left(z\right)$, ${\mathsf{\Psi }}^k\left(z\right)$, and ${\mathsf{\Theta }}^k\left(z\right)$ represent the through-the-thickness amplitudes of the corresponding field variables and constitute the unknown functions of the coupled thermo-magneto-elastic problem. The harmonic coefficients are defined in terms of the half-wave numbers m and n and the panel dimensions a and b as

$$\overline{\alpha }={\displaystyle \frac{m\pi }{a}},\overline{\beta }={\displaystyle \frac{n\pi }{b}}.$$

Imposition of Navier harmonic forms permits the satisfaction of simply-supported boundary conditions and the exact computation of the derivatives along $\alpha$ and $\beta$ directions. In this way, the 3D governing equations for the 3D thermo-magneto-elastic problem can be written in explicit form as:

$$\begin{array}{cc}& \left(-{\displaystyle \frac{H_{\beta }{C}_{55}^k}{H_{\alpha }{R}_{\alpha }^2}}-{\displaystyle \frac{C_{55}^k}{R_{\alpha }{R}_{\beta }}}-{\overline{\alpha }}^2{\displaystyle \frac{C_{11}^k{H}_{\beta }}{H_{\alpha }}}-{\overline{\beta }}^2{\displaystyle \frac{C_{66}^k{H}_{\alpha }}{H_{\beta }}}\right)U^k+\left(-\overline{\alpha }\overline{\beta }{C}_{12}^k-\overline{\alpha }\overline{\beta }{C}_{66}^k\right)V^k+\hfill \\ & +\left(\overline{\alpha }{\displaystyle \frac{C_{11}^k{H}_{\beta }}{H_{\alpha }{R}_{\alpha }}}+\overline{\alpha }{\displaystyle \frac{C_{12}^k}{R_{\beta }}}+\overline{\alpha }{\displaystyle \frac{C_{55}^k{H}_{\beta }}{H_{\alpha }{R}_{\alpha }}}+\overline{\alpha }{\displaystyle \frac{C_{55}^k}{R_{\beta }}}\right)W^k+\left({\displaystyle \frac{C_{55}^k{H}_{\beta }}{R_{\alpha }}}+{\displaystyle \frac{C_{55}^k{H}_{\alpha }}{R_{\beta }}}\right)U_{,z}^k+\hfill \\ & +\left(\overline{\alpha }{C}_{13}^k{H}_{\beta }+\overline{\alpha }{C}_{55}^k{H}_{\beta }\right)W_{,z}^k+\left(C_{55}^k{H}_{\alpha }{H}_{\beta }\right)U_{,zz}^k+\left(+2\overline{\alpha }{\displaystyle \frac{r_{15}^k{H}_{\beta }}{H_{\alpha }{R}_{\alpha }}}+\overline{\alpha }{\displaystyle \frac{r_{15}^k}{R_{\beta }}}\right){\mathsf{\Psi }}^k+\hfill \\ & +\left(\overline{\alpha }{r}_{31}^k{H}_{\beta }+\overline{\alpha }{r}_{15}^k{H}_{\beta }\right){\mathsf{\Psi }}_{,z}^k-\left(\overline{\alpha }{\lambda }_{11}^k{H}_{\beta }\right){\mathsf{\Theta }}^k=0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}& \left(-{\displaystyle \frac{H_{\alpha }{C}_{44}^k}{H_{\beta }{R}_{\beta }^2}}-{\displaystyle \frac{C_{44}^k}{R_{\alpha }{R}_{\beta }}}-{\overline{\alpha }}^2{\displaystyle \frac{C_{66}^k{H}_{\beta }}{H_{\alpha }}}-{\overline{\beta }}^2{\displaystyle \frac{C_{22}^k{H}_{\alpha }}{H_{\beta }}}\right)V^k+\left(-\overline{\alpha }\overline{\beta }{C}_{12}^k-\overline{\alpha }\overline{\beta }{C}_{66}^k\right)U^k+\hfill \\ & +\left(\overline{\beta }{\displaystyle \frac{C_{44}^k{H}_{\alpha }}{H_{\beta }{R}_{\beta }}}+\overline{\beta }{\displaystyle \frac{C_{44}^k}{R_{\alpha }}}+\overline{\beta }{\displaystyle \frac{C_{22}^k{H}_{\alpha }}{H_{\beta }{R}_{\beta }}}+\overline{\beta }{\displaystyle \frac{C_{12}^k}{R_{\alpha }}}\right)W^k+\left({\displaystyle \frac{C_{44}^k{H}_{\alpha }}{R_{\beta }}}+{\displaystyle \frac{C_{44}^k{H}_{\beta }}{R_{\alpha }}}\right)V_{,z}^k+\hfill \\ & +\left(\overline{\beta }{C}_{44}^k{H}_{\alpha }+\overline{\beta }{C}_{23}^k{H}_{\alpha }\right)W_{,z}^k+\left(C_{44}^k{H}_{\alpha }{H}_{\beta }\right)V_{,zz}^k+\left(+2\overline{\beta }{\displaystyle \frac{r_{24}^k{H}_{\alpha }}{H_{\beta }{R}_{\beta }}}+\overline{\beta }{\displaystyle \frac{r_{24}^k}{R_{\alpha }}}\right){\mathsf{\Psi }}^k+\hfill \\ & +\left(\overline{\beta }{r}_{32}^k{H}_{\alpha }+\overline{\beta }{r}_{24}^k{H}_{\alpha }\right){\mathsf{\Psi }}_{,z}^k-\left(\overline{\beta }{\lambda }_{22}^k{H}_{\alpha }\right){\mathsf{\Theta }}^k=0,\hfill \end{array}$$

(38)

$$\begin{array}{cc}& \left({\displaystyle \frac{C_{13}^k}{R_{\alpha }{R}_{\beta }}}+{\displaystyle \frac{C_{23}^k}{R_{\alpha }{R}_{\beta }}}-{\displaystyle \frac{C_{11}^k{H}_{\beta }}{H_{\alpha }{R}_{\alpha }^2}}-{\displaystyle \frac{2C_{12}^k}{R_{\alpha }{R}_{\beta }}}-{\displaystyle \frac{C_{22}^k{H}_{\alpha }}{H_{\beta }{R}_{\beta }^2}}-{\overline{\alpha }}^2{\displaystyle \frac{C_{55}^k{H}_{\beta }}{H_{\alpha }}}-{\overline{\beta }}^2{\displaystyle \frac{C_{44}^k{H}_{\alpha }}{H_{\beta }}}\right)W^k+\hfill \\ & +\left(\overline{\alpha }{\displaystyle \frac{C_{55}^k{H}_{\beta }}{H_{\alpha }{R}_{\alpha }}}-\overline{\alpha }{\displaystyle \frac{C_{13}^k}{R_{\beta }}}+\overline{\alpha }{\displaystyle \frac{C_{11}^k{H}_{\beta }}{H_{\alpha }{R}_{\alpha }}}+\overline{\alpha }{\displaystyle \frac{C_{12}^k}{R_{\beta }}}\right)U^k+\left(\overline{\beta }{\displaystyle \frac{C_{44}^k{H}_{\alpha }}{H_{\beta }{R}_{\beta }}}-\overline{\beta }{\displaystyle \frac{C_{23}^k}{R_{\alpha }}}+\overline{\beta }{\displaystyle \frac{C_{22}^k{H}_{\alpha }}{H_{\beta }{R}_{\beta }}}+\overline{\beta }{\displaystyle \frac{C_{12}^k}{R_{\alpha }}}\right)V^k+\hfill \\ & +\left(-\overline{\alpha }{C}_{55}^k{H}_{\beta }-\overline{\alpha }{C}_{13}^k{H}_{\beta }\right)U_{,z}^k+\left(-\overline{\beta }{C}_{44}^k{H}_{\alpha }-\overline{\beta }{C}_{23}^k{H}_{\alpha }\right)V_{,z}^k+\left({\displaystyle \frac{C_{33}^k{H}_{\beta }}{R_{\alpha }}}+{\displaystyle \frac{C_{33}^k{H}_{\alpha }}{R_{\beta }}}\right)W_{,z}^k+\hfill \\ & +\left(-{\overline{\alpha }}^2{\displaystyle \frac{r_{15}^k{H}_{\beta }}{H_{\alpha }}}-{\overline{\beta }}^2{\displaystyle \frac{r_{24}^k{H}_{\alpha }}{H_{\beta }}}\right){\mathsf{\Psi }}^k+\left(-{\displaystyle \frac{r_{31}^k{H}_{\beta }}{R_{\alpha }}}-{\displaystyle \frac{r_{32}^k{H}_{\alpha }}{R_{\beta }}}+{\displaystyle \frac{r_{33}^k{H}_{\beta }}{R_{\alpha }}}+{\displaystyle \frac{r_{33}^k{H}_{\alpha }}{R_{\beta }}}\right){\mathsf{\Psi }}_{,z}^k+\hfill \\ & +\left(C_{33}^k{H}_{\alpha }{H}_{\beta }\right)W_{,zz}^k+\left(r_{33}^k{H}_{\alpha }{H}_{\beta }\right){\mathsf{\Psi }}_{,zz}^k+\left({\displaystyle \frac{{\lambda }_{11}^k{H}_{\beta }}{R_{\alpha }}}+{\displaystyle \frac{{\lambda }_{22}^k{H}_{\alpha }}{R_{\beta }}}-{\displaystyle \frac{{\lambda }_{33}^k{H}_{\beta }}{R_{\alpha }}}-{\displaystyle \frac{{\lambda }_{33}^k{H}_{\alpha }}{R_{\beta }}}\right){\mathsf{\Theta }}^k+\hfill \\ & -\left({\lambda }_{33}{H}_{\alpha }{H}_{\beta }\right){\mathsf{\Theta }}_{,z}^k=0,\hfill \end{array}$$

(39)

$$\begin{array}{cc}& \left(\overline{\alpha }{\displaystyle \frac{r_{15}^k}{H_{\alpha }^2{R}_{\alpha }}}\right)U^k+\left(\overline{\beta }{\displaystyle \frac{r_{24}^k}{H_{\beta }^2{R}_{\beta }}}\right)V^k+\left(-{\overline{\alpha }}^2{\displaystyle \frac{r_{15}^k}{H_{\alpha }^2}}-{\overline{\beta }}^2{\displaystyle \frac{r_{24}^k}{H_{\beta }^2}}\right)W^k+\left(-\overline{\alpha }{\displaystyle \frac{r_{15}^k}{H_{\alpha }}}-\overline{\alpha }{\displaystyle \frac{r_{31}^k}{H_{\alpha }}}\right)U_{,z}^k+\hfill \\ & +\left(-\overline{\beta }{\displaystyle \frac{r_{24}^k}{H_{\beta }}}-\overline{\beta }{\displaystyle \frac{r_{32}^k}{H_{\beta }}}\right)V_{,z}^k+\left({\displaystyle \frac{r_{31}^k}{H_{\alpha }{R}_{\alpha }}}+{\displaystyle \frac{r_{32}^k}{H_{\beta }{R}_{\beta }}}\right)W_{,z}^k+\left(r_{33}^k\right)W_{,zz}^k+\left({\overline{\alpha }}^2{\displaystyle \frac{{\mu }_{11}^k}{H_{\alpha }^2}}+{\overline{\beta }}^2{\displaystyle \frac{{\mu }_{22}^k}{H_{\beta }^2}}\right){\mathsf{\Psi }}^k+\hfill \\ & -\left({\mu }_{33}^k\right){\mathsf{\Psi }}_{,zz}^k+\hfill \\ & +\left(g_3^k\right){\mathsf{\Theta }}_{,z}^k=0,\hfill \end{array}$$

(40)

$$\left(-{\overline{\alpha }}^2{\displaystyle \frac{k_{11}^k}{H_{\alpha }^2}}-{\overline{\beta }}^2{\displaystyle \frac{k_{22}^k}{H_{\beta }^2}}\right){\mathsf{\Theta }}^k+\left(k_{33}^k\right){\mathsf{\Theta }}_{,zz}^k=0.$$

(41)

Subscripts $,z$ and $,zz$ indicate first and second derivatives in z, respectively. For the sake of compactness, all the coefficient groups enclosed within parentheses in Eqs. (37)–(41) are replaced by the coefficients $A_i^k$ and $J_i^k$. Accordingly, the governing system can be rewritten in the following compact form:

$$\begin{array}{cc}\hfill & A_1^k{U}^k+A_2^k{V}^k+A_3^k{W}^k+A_4^k{\mathsf{\Psi }}^k+J_1^k{\mathsf{\Theta }}^k+A_5^k{U}_{,z}^k+A_6^k{W}_{,z}^k+A_7^k{\mathsf{\Psi }}_{,z}^k+A_8^k{U}_{,zz}^k=0,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & A_9^k{U}^k+A_{10}^k{V}^k+A_{11}^k{W}^k+A_{12}^k{\mathsf{\Psi }}^k+J_2^k{\mathsf{\Theta }}^k+A_{13}^k{V}_{,z}^k+A_{14}^k{W}_{,z}^k+A_{15}^k{\mathsf{\Psi }}_{,z}^k+A_{16}^k{V}_{,zz}^k=0,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill & A_{17}^k{U}^k+A_{18}^k{V}^k+A_{19}^k{W}^k+A_{20}^k{\mathsf{\Psi }}^k+J_3^k{\mathsf{\Theta }}^k+A_{21}^k{U}_{,z}^k+A_{22}^k{V}_{,z}^k+A_{23}^k{W}_{,z}^k+A_{24}^k{\mathsf{\Psi }}_{,z}^k+J_4^k{\mathsf{\Theta }}_{,z}^k+\hfill \\ \hfill & +A_{25}^k{W}_{,zz}^k+A_{26}^k{\mathsf{\Psi }}_{,zz}^k=0,\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & A_{27}^k{U}^k+A_{28}^k{V}^k+A_{29}^k{W}^k+A_{30}^k{\mathsf{\Psi }}^k+A_{31}^k{U}_{,z}^k+A_{32}^k{V}_{,z}^k+A_{33}^k{W}_{,z}^k+J_5^k{\mathsf{\Theta }}_{,z}^k+\hfill \\ \hfill & +A_{34}^k{W}_{,zz}^k+A_{35}^k{\mathsf{\Psi }}_{,zz}^k=0,\hfill \end{array}$$

(45)

$$J_6^k{\mathsf{\Theta }}^k+J_7^k{\mathsf{\Theta }}_{,zz}^k=0.$$

(46)

Eqs. (42)–(46) form a coupled system of second-order ordinary differential equations along the thickness coordinate z. The exponential matrix method requires both constant coefficients and a first-order differential formulation. However, the coefficients $A_i^k$ depend on the curvature functions $H_{\alpha }$ and $H_{\beta }$, which vary through the thickness. To overcome this issue, each physical layer is subdivided into a set of sufficiently thin mathematical layers in which the thickness coordinate can be assumed constant. Consequently, the curvature terms, together with the associated coefficients $A_i^k$ and $J_i^k$, become constant within each mathematical layer.

Let M denote the total number of mathematical layers used to discretize the thickness, and let j identify a generic mathematical layer. Since the governing equations have the same form in every mathematical layer, the following derivation is presented for a generic layer j. The second-order system is subsequently transformed into an equivalent first-order one by introducing the first derivatives of the primary variables through a standard doubling procedure. The resulting equations are

$$A_8^j{U}_{,z}^j=A_8^j{U}_{,z}^j,$$

(47)

$$A_{16}^j{V}_{,z}^j=A_{16}^j{V}_{,z}^j,$$

(48)

$$P_1^j{W}_{,z}^j=P_1^j{W}_{,z}^j,$$

(49)

$$P_1^j{\mathsf{\Psi }}_{,z}^j=P_1^j{\mathsf{\Psi }}_{,z}^j,$$

(50)

$$J_7^j{\mathsf{\Theta }}_{,z}^j=J_7^j{\mathsf{\Theta }}_{,z}^j,$$

(51)

$$\begin{array}{cc}\hfill A_8^j{U}_{,zz}^j=& -A_1^j{U}^j-A_2^j{V}^j-A_3^j{W}^j-A_4^j{\mathsf{\Psi }}^j-J_1^k{\mathsf{\Theta }}^k-A_5^j{U}_{,z}^j-A_6^j{W}_{,z}^j-A_7^j{\mathsf{\Psi }}_{,z}^j,\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill A_{16}^j{V}_{,zz}^j=& -A_9^j{U}^j-A_{10}^j{V}^j-A_{11}^j{W}^j-A_{12}^j{\mathsf{\Psi }}^j-J_2^k{\mathsf{\Theta }}^k-A_{13}^j{V}_{,z}^j-A_{14}^j{W}_{,z}^j-A_{15}^j{\mathsf{\Psi }}_{,z}^j,\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill P_1^j{W}_{,zz}^j=& -P_2^j{U}^j-P_3^j{V}^j-P_4^j{W}^j-P_5^j{\mathsf{\Theta }}^j-P_6^j{\mathsf{\Psi }}^j-P_7^j{U}_{,z}^j-P_8^j{V}_{,z}^j-P_9^j{W}_{,z}^j-P_{10}^j{\mathsf{\Psi }}_{,z}^j-P_{11}^j{\mathsf{\Theta }}_{,z}^j,\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill P_1^j{\mathsf{\Psi }}_{,zz}^j=& -P_{12}^j{U}^j-P_{13}^j{V}^j-P_{14}^j{W}^j-P_{15}^j{\mathsf{\Psi }}^j-P_{16}^j{U}_{,z}^j-P_{17}^j{V}_{,z}^j-P_{18}^j{W}_{,z}^j-P_{19}^j{\mathsf{\Theta }}_{,z}^j,\hfill \end{array}$$

(55)

$$J_7^j{\mathsf{\Theta }}_{,zz}^j=-J_6^j{\mathsf{\Theta }}^j,$$

(56)

The above first-order representation is directly suitable for the application of the exponential matrix method. Accordingly, Eqs. (47)–(56) can be assembled into the following state-space form:

$${\mathbf{D}}^j{\mathbf{X}}_{,z}^j={\mathbf{A}}^j{\mathbf{X}}^j\Rightarrow {\mathbf{X}}_{,z}^j={\mathbf{A}}^{*j}{\mathbf{X}}^j,$$

(57)

where ${\mathbf{D}}^j$ and ${\mathbf{A}}^j$ are two $10\times 10$ coefficient matrices.

${\mathbf{X}}^j={\left\{U^j{V}^j{W}^j{\mathsf{\Psi }}^j{\mathsf{\Theta }}^j{U}_{,z}^j{V}_{,z}^j{W}_{,z}^j{\mathsf{\Psi }}_{,z}^j{\mathsf{\Theta }}_{,z}^j\right\}}^T$ is the $10\times 1$ state vector, and ${\mathbf{X}}_{,z}^j$ denotes its derivative with respect to the thickness coordinate. Owing to the doubling procedure, the first derivatives of the original unknowns are incorporated into the state vector, leading to a first-order system with ten independent variables. The first-order state-space system in Eq. (57) admits an analytical solution based on the exponential matrix method. The solution relates the state vector at the bottom of the generic mathematical layer ($z=0$) to that at its top ($z=h_j$), and can be written as

$${\mathbf{X}}^j\left(h_j\right)={\mathbf{A}}^{**j}{\mathbf{X}}^j\left(0\right),\mathit{where}{\mathbf{A}}^{**j}=exp\left({\mathbf{A}}^{*j}{h}_j\right)={\displaystyle \sum _{i=0}^N}{\displaystyle \frac{{\left({\mathbf{A}}^{*j}\right)}^i{h}_j^i}{i!}}.$$

(58)

The matrix ${\mathbf{A}}^{**j}$ is the $10\times 10$ exponential matrix associated with the j-th mathematical layer, whereas $h_j$ denotes its thickness. In the present formulation, the exponential matrix is evaluated through the truncated Taylor series expansion reported in Eq. (58), where N is the order of the approximation. The first term of the series corresponds to the identity matrix, namely ${\left({\mathbf{A}}^{*j}\right)}^0=\mathbf{I}$, with $\mathbf{I}$ denoting the $10\times 10$ identity matrix.

The layer-wise description of the multilayered structure is achieved by enforcing interlaminar continuity conditions at every interface between two adjacent mathematical layers. These continuity conditions involve both primary variables (displacements, magnetic potential, and over-temperature) and secondary variables (transverse shear and normal stresses, transverse magnetic induction, and transverse heat flux). Their enforcement guarantees the continuity of the thermo-magneto-elastic fields throughout the thickness direction. The interlaminar continuity conditions can be expressed as

$$u_b^j=u_t^{j-1},v_b^j=v_t^{j-1},w_b^j=w_t^{j-1},{\psi }_b^j={\psi }_t^{j-1},{\theta }_b^j={\theta }_t^{j-1},$$

(59)

$${\sigma }_{\alpha z_b}^j={\sigma }_{\alpha z_t}^{j-1},{\sigma }_{\beta z_b}^j={\sigma }_{\beta z_t}^{j-1},{\sigma }_{z{z}_b}^j={\sigma }_{z{z}_t}^{j-1},{\mathcal{B}}_{z_b}^j={\mathcal{B}}_{z_t}^{j-1},q_{z_b}^j=q_{z_t}^{j-1},$$

(60)

where the subscripts b and t denote the bottom surface of the j-th mathematical layer and the top surface of the adjacent $(j-1)$-th mathematical layer, respectively. By introducing the Navier harmonic expansions in Eqs. (32)–(36), together with the curvature relations in Eq. (6) and the constitutive equations (11)–(22), the interlaminar continuity conditions can be rewritten in terms of the harmonic amplitudes as follows:

$$U_b^j=U_t^{j-1},$$

(61)

$$V_b^j=V_t^{j-1},$$

(62)

$$W_b^j=W_t^{j-1},$$

(63)

$${\mathsf{\Psi }}_b^j={\mathsf{\Psi }}_t^{j-1},$$

(64)

$${\mathsf{\Theta }}_b^j={\mathsf{\Theta }}_t^{j-1},$$

(65)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{C_{55}^j}{H_{{\alpha }_b}^j}}\overline{\alpha }{W}_b^j+C_{55}^j{U}_{,z_b}^j-{\displaystyle \frac{C_{55}^j}{H_{{\alpha }_b}^j{R}_{\alpha }}}{U}_b^j+{\displaystyle \frac{r_{15}^j}{H_{{\alpha }_b}^j}}\overline{\alpha }{\mathsf{\Psi }}_b^j=\hfill \\ & +{\displaystyle \frac{C_{55}^{j-1}}{H_{{\alpha }_t}^{j-1}}}\overline{\alpha }{W}_t^{j-1}+C_{55}^{j-1}{U}_{,z_t}^{j-1}-{\displaystyle \frac{C_{55}^{j-1}}{H_{{\alpha }_t}^{j-1}{R}_{\alpha }}}{U}_t^{j-1}+{\displaystyle \frac{r_{15}^{j-1}}{H_{{\alpha }_t}^{j-1}}}\overline{\alpha }{\mathsf{\Psi }}_t^{j-1},\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{C_{44}^j}{H_{{\beta }_b}^j}}\overline{\beta }{W}_b^j+C_{44}^j{V}_{,z_b}^j-{\displaystyle \frac{C_{44}^j}{H_{{\beta }_b}^j{R}_{\beta }}}{V}_b^j+{\displaystyle \frac{r_{24}^j}{H_{{\beta }_b}^j}}\overline{\beta }{\mathsf{\Psi }}_b^j=\hfill \\ & +{\displaystyle \frac{C_{44}^{j-1}}{H_{{\beta }_t}^{j-1}}}\overline{\beta }{W}_t^{j-1}+C_{44}^{j-1}{V}_{,z_t}^{j-1}-{\displaystyle \frac{C_{44}^{j-1}}{H_{{\beta }_t}^{j-1}{R}_{\beta }}}{V}_t^{j-1}+{\displaystyle \frac{r_{24}^{j-1}}{H_{{\beta }_t}^{j-1}}}\overline{\beta }{\mathsf{\Psi }}_t^{j-1},\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill & -{\displaystyle \frac{C_{13}^j}{H_{{\alpha }_b}^j}}\overline{\alpha }{U}_b^j+{\displaystyle \frac{C_{13}^j}{H_{{\alpha }_b}^j{R}_{\alpha }}}{W}_b^j-{\displaystyle \frac{C_{23}^j}{H_{{\beta }_b}^j}}\overline{\beta }{V}_b^j+{\displaystyle \frac{C_{23}^j}{H_{{\beta }_b}^j{R}_{\beta }}}{W}_b^j-{\lambda }_{33}^j{\mathsf{\Theta }}_b^j+C_{33}^j{W}_{,z_b}^j+r_{33}^j{\mathsf{\Psi }}_{,z_b}^j=\hfill \\ & -{\displaystyle \frac{C_{13}^{j-1}}{H_{{\alpha }_t}^{j-1}}}\overline{\alpha }{U}_t^{j-1}+{\displaystyle \frac{C_{13}^{j-1}}{H_{{\alpha }_t}^{j-1}{R}_{\alpha }}}{W}_t^{j-1}-{\displaystyle \frac{C_{23}^{j-1}}{H_{{\beta }_t}^{j-1}}}\overline{\beta }{V}_t^{j-1}+{\displaystyle \frac{C_{23}^{j-1}}{H_{{\beta }_t}^{j-1}{R}_{\beta }}}{W}_t^{j-1}-{\lambda }_{33}^{j-1}{\mathsf{\Theta }}_t^{j-1}+\hfill \\ & +C_{33}^{j-1}{W}_{,z_t}^{j-1}+r_{33}^{j-1}{\mathsf{\Psi }}_{,z_t}^{j-1},\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill & -{\displaystyle \frac{r_{31}^j}{H_{{\alpha }_b}^j}}\overline{\alpha }{U}_b^j+{\displaystyle \frac{r_{31}^j}{H_{{\alpha }_b}^j{R}_{\alpha }}}{W}_b^j-{\displaystyle \frac{r_{32}^j}{H_{{\beta }_b}^j}}\overline{\beta }{V}_b^j+{\displaystyle \frac{r_{32}^j}{H_{{\beta }_b}^j{R}_{\beta }}}{W}_b^j+g_3^j{\mathsf{\Theta }}_b^j+r_{33}^j{W}_{,z_b}^j-{\mu }_{33}^j{\mathsf{\Psi }}_{,z_b}^j=\hfill \\ & -{\displaystyle \frac{r_{31}^{j-1}}{H_{{\alpha }_t}^{j-1}}}\overline{\alpha }{U}_t^{j-1}+{\displaystyle \frac{r_{31}^{j-1}}{H_{{\alpha }_t}^{j-1}{R}_{\alpha }}}{W}_t^{j-1}-{\displaystyle \frac{r_{32}^{j-1}}{H_{{\beta }_t}^{j-1}}}\overline{\beta }{V}_t^{j-1}+{\displaystyle \frac{r_{32}^{j-1}}{H_{{\beta }_t}^{j-1}{R}_{\beta }}}{W}_t^{j-1}+g_3^{j-1}{\mathsf{\Theta }}_t^{j-1}+\hfill \\ & +r_{33}^{j-1}{W}_{,z_t}^{j-1}-{\mu }_{33}^{j-1}{\mathsf{\Psi }}_{,z_t}^{j-1},\hfill \end{array}$$

(69)

$$-k_{33}^j{\mathsf{\Theta }}_{,z_b}^j=-k_{33}^{j-1}{\mathsf{\Theta }}_{,z_t}^{j-1}.$$

(70)

The explicit interlaminar continuity conditions reported in Eqs. (61)–(70) can be assembled into the following matrix equation:

$${\mathbf{X}}^j\left(0\right)={\mathbf{T}}^{j,j-1}{\mathbf{X}}^{j-1}\left(h_{j-1}\right),$$

(71)

where ${\mathbf{T}}^{j,j-1}$ is the constant $10\times 10$ transfer matrix relating the state vector at the bottom (${\tilde{z}}_j=0$) of the j-th mathematical layer to the state vector at the top (${\tilde{z}}_{j-1}=h_{j-1}$) of the adjacent $(j-1)$-th layer.

By recursively combining the exponential matrix solution in Eq. (58) with the transfer relation given in Eq. (71), the analytical solution can be propagated through the entire thickness of the multilayered structure. The resulting global relation is

$${\mathbf{X}}^M\left(h_M\right)={\mathbf{A}}^{**M}{\mathbf{T}}^{M,M-1}\dots {\mathbf{T}}^{2,1}{\mathbf{A}}^{**1}{\mathbf{X}}^1\left(0\right)={\mathbf{H}}_m{\mathbf{X}}^1\left(0\right),$$

(72)

where ${\mathbf{H}}_m$ is the global transfer matrix obtained by successively assembling the exponential matrices ${\mathbf{A}}^{**j}$ and the transfer matrices ${\mathbf{T}}^{j,j-1}$ corresponding to all the mathematical layers. Regardless of the number of mathematical layers M adopted to discretize the thickness or the Taylor expansion order N employed to evaluate the exponential matrix, ${\mathbf{H}}_m$ always remains a $10\times 10$ matrix. Consequently, the computational effort is significantly reduced while preserving the full three-dimensional description of the multilayered structure. The matrix ${\mathbf{H}}_m$ directly relates the unknown state vector at the bottom surface of the first mathematical layer (${\tilde{z}}_1=0$) to that at the top surface of the last mathematical layer (${\tilde{z}}_M=h_M$).

Simply-supported boundary conditions are adopted along the panel edges since they are the only boundary conditions compatible with the Navier harmonic expansions and therefore allow the derivation of a closed-form analytical solution. For the generic mathematical layer j, the edge boundary conditions are expressed as

$$\begin{array}{cc}\hfill & v^j=0,w^j=0,{\psi }^j=0,{\theta }^j=0,{\sigma }_{\alpha \alpha }^j=0\mathrm{for}\alpha =0,a,\hfill \\ \hfill & u^j=0,w^j=0,{\psi }^j=0,{\theta }^j=0,{\sigma }_{\beta \beta }^j=0\mathrm{for}\beta =0,b.\hfill \end{array}$$

(73)

Once the edge boundary conditions have been satisfied, loading boundary conditions must be prescribed at the lower and upper surfaces of the multilayered structure. Depending on the physical problem under investigation, mechanical, magnetic and thermal quantities can be assigned in terms of transverse normal stress, magnetic potential, transverse magnetic induction and over-temperature. Four possible combinations are considered:

$${\sigma }_{zz}={\overline{P}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),\theta =\overline{\mathsf{\Theta }}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),\psi =\overline{\mathsf{\Psi }}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right)\mathrm{for}\tilde{z}=0,h,$$

(74)

$${\sigma }_{zz}={\overline{P}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),\theta =\overline{\mathsf{\Theta }}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),{\mathcal{B}}_z={\overline{\mathcal{B}}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right)\mathrm{for}\tilde{z}=0,h,$$

(75)

$${\sigma }_{zz}={\overline{P}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),q_z={\overline{Q}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),\psi =\overline{\mathsf{\Psi }}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right)\mathrm{for}\tilde{z}=0,h,$$

(76)

$${\sigma }_{zz}={\overline{P}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),q_z={\overline{Q}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),{\mathcal{B}}_z={\overline{\mathcal{B}}}_{z}sin\left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right)\mathrm{for}\tilde{z}=0,h.$$

(77)

In the previous expressions, the overbar denotes a prescribed boundary value. By substituting the Navier harmonic expansions together with the thermo-magneto-elastic constitutive equations into Eqs. (74)–(77), the loading boundary conditions can be assembled into the following matrix form:

$${\mathbf{B}}_b^1{\mathbf{X}}^1\left(0\right)={\mathbf{P}}_b,$$

(78)

$${\mathbf{B}}_t^M{\mathbf{X}}^M\left(h_M\right)={\mathbf{P}}_t,$$

(79)

where ${\mathbf{B}}_b^1$ and ${\mathbf{B}}_t^M$ are the $5\times 10$ loading matrices associated with the bottom surface of the first mathematical layer and the top surface of the last mathematical layer, respectively. Likewise, ${\mathbf{P}}_b$ and ${\mathbf{P}}_t$ are the corresponding $5\times 1$ loading vectors applied at the lower ($\tilde{z}=0$) and upper ($\tilde{z}=h$) surfaces of the multilayered structure.

By combining Eqs. (78) and (79) with the global transfer relation given in Eq. (72), the final algebraic system becomes

$$\mathbf{E}{\mathbf{X}}^1\left(0\right)=\mathbf{P},\mathit{where}\mathbf{E}=\left[\begin{array}{c}{\mathbf{B}}_t^M{\mathbf{H}}_m\\ {\mathbf{B}}_b^1\end{array}\right],\mathbf{P}=\left\{\begin{array}{c}{\mathbf{P}}_t\\ {\mathbf{P}}_b\end{array}\right\}.$$

(80)

The solution of Eq. (80) provides the unknown state vector at the bottom surface of the first mathematical layer. Subsequently, the state vector throughout the thickness can be reconstructed by recursively applying Eqs. (58) and (71). Once the primary variables have been determined, the secondary variables, namely strains, stresses, magnetic induction components and heat flux components, are readily evaluated from the constitutive and congruence relations reported in Eqs. (11)–(22) and Eqs. (23)–(31), respectively.

3. Results

This section presents the numerical results obtained using the proposed 3D layer-wise coupled thermo-magneto-elastic shell model for both flat and curved multilayered structures. Throughout the paper, the formulation is referred to as the 3D-u-$\theta$-$\psi$ model. This designation highlights the main characteristics of the proposed analytical approach, namely a fully three-dimensional formulation (3D) in which the displacement components (u), the over-temperature ($\theta$), and the magnetic potential ($\psi$) are selected as the primary unknowns.

The Results section is organized into two subsections. The first one is devoted to preliminary assessments, where the proposed 3D-u-$\theta$-$\psi$ model is validated through comparisons with thermo-magneto-elastic solutions available in the literature. These comparisons also allow the identification of suitable values for both the exponential matrix order and the number of mathematical layers required to achieve accurate predictions. The second subsection presents a series of new benchmark results for multilayered plates, cylinders, cylindrical panels, and spherical shells with sandwich configurations and different thickness ratios. The geometries considered throughout this section, together with the corresponding mixed curvilinear orthogonal reference systems, are illustrated in Figure 1.

3.1. Preliminary Results

This subsection presents four preliminary assessments aimed at validating the proposed 3D layer-wise coupled thermo-magneto-elastic shell model (3D-u-$\theta$-$\psi$). Each assessment has been designed to verify specific features of the formulation by comparison with reference solutions available in the literature.

The first assessment considers a multilayered rectangular plate and is intended to validate the thermo-magneto-elastic coupling together with the material layer and thickness effects for a flat geometry. The second assessment investigates a multilayered spherical shell, extending the validation to curved structures and demonstrating the capability of the proposed model to accurately account for curvature effects in addition to thermo-magneto-elastic coupling and layer-dependent phenomena. The third assessment is devoted to a single-layered plate in which only the magneto-elastic coupling is considered, while thermal effects are neglected. This case provides an independent validation of the magnetic formulation as well as the material and thickness effects. Finally, the fourth assessment examines a two-layered square plate in which only the thermo-elastic coupling is retained, whereas magnetic effects are neglected. This last case is employed to validate the thermal formulation together with the corresponding layer-dependent effects.

The last two assessments have been introduced because, to the best of the authors’ knowledge, no three-dimensional analytical solutions fully coupling the elastic, magnetic, and thermal fields are currently available in the literature. Consequently, the thermo-elastic and magneto-elastic couplings are validated separately through these benchmark problems.

For the third and fourth assessments, the numerical parameters are fixed to $M=300$ mathematical layers and a Taylor expansion order $N=3$ for the exponential matrix, consistently with the convergence analysis previously proposed for three-dimensional static elastic problems by Brischetto [57]. For the first and second assessments, an additional convergence study is reported to justify the adopted values of M and N for the coupled thermo-magneto-elastic problem.

Both the first and the second assessments consider a three-layered Adaptive Wood/Soft Adaptive Wood/Adaptive Wood sandwich configuration. The two outer Adaptive Wood layers have a thickness equal to $0.3h$, whereas the Soft Adaptive Wood core has a thickness of $0.4h$, with h denoting the total thickness of the structure. Both Adaptive Wood and Soft Adaptive Wood are composite materials obtained by combining wood with barium titanate ($BaTi{O}_3$) and cobalt ferrite ($CoF{e}_2{O}_4$). Their elastic, magnetic, and thermal material properties are summarized in Table 1. Following the convention adopted in the reference solutions of [58,59], the magnetic permeability coefficients ${\mu }_{11}$ and ${\mu }_{22}$ are assumed to be negative in the preliminary assessments. In contrast, positive values will be adopted for the new benchmark problems presented in Section 3.2 in order to better represent the physical behavior of the material.

In the first preliminary assessment (PA1), a multilayered simply-supported moderately thick rectangular plate is considered. The corresponding geometrical parameters, half-wave numbers, and amplitudes of the prescribed loading boundary conditions are reported in the first column of Table 2. As reference solution, the Equivalent Layer-Wise (ELW) theory developed by Tornabene et al. [58] is adopted.

Table 3 reports the convergence analysis performed on the transverse displacement w, whereas Table 4 compares the results obtained with the proposed 3D-u-$\theta$-$\psi$ model against the reference solution. The convergence study demonstrates that an exponential matrix order equal to $N=3$ together with $M=300$ mathematical layers provides fully converged results. Therefore, these values can be regarded as conservative choices, ensuring an accurate evaluation of the thermo-magneto-elastic response.

The reference values reported in Table 4 were extracted from the through-the-thickness distributions presented by Tornabene et al. [58] at the in-plane coordinates $\alpha =0.25a$ and $\beta =0.25b$. An excellent agreement is observed between the proposed formulation and the reference model for all the investigated quantities, including displacement components, magnetic potential, over-temperature, in-plane normal stress, transverse shear stress, transverse normal stress, in-plane normal magnetic induction, and transverse normal heat flux. The comparison confirms the accuracy of the proposed 3D-u-$\theta$-$\psi$ model and its capability to correctly capture the thermo-magneto-elastic coupling together with the material layer and thickness effects in multilayered plate structures.

For the second preliminary assessment (PA2), a multilayered simply-supported moderately thick spherical shell is investigated. The corresponding geometrical parameters, half-wave numbers, and prescribed loading boundary conditions are reported in the second column of Table 2. As in PA1, the reference solution is provided by the Equivalent Layer-Wise (ELW) theory developed by Tornabene et al. [58].

Table 5 presents the convergence analysis performed on the transverse displacement w for the spherical shell configuration. Also for this curved geometry, an exponential matrix order equal to $N=3$ combined with $M=300$ mathematical layers provides converged results. A slightly slower convergence is observed with respect to the plate case because of the additional curvature terms appearing in the governing equations. Nevertheless, the selected values of N and M guarantee an accurate evaluation of the thermo-magneto-elastic response.

The reference values listed in Table 6 were extracted from the through-the-thickness distributions reported by Tornabene et al. [58] at the in-plane coordinates $\alpha =0.25a$ and $\beta =0.25b$. An excellent agreement is achieved for the in-plane and transverse displacement components, magnetic potential, over-temperature, transverse shear and transverse normal stresses, in-plane normal stress, and in-plane normal magnetic induction. These comparisons validate the capability of the proposed 3D-u-$\theta$-$\psi$ model to accurately reproduce the thermo-magneto-elastic coupling together with the material layer, thickness, and curvature effects in multilayered spherical shells.

The small discrepancies observed in Table 4 and Table 6 mainly arise from the fact that the ELW reference values were extracted from the graphical distributions reported by Tornabene et al. [58]. Consequently, the tabulated reference data are affected by the unavoidable approximations associated with the graphical reading procedure.

For the third preliminary assessment (PA3), a single-layered simply-supported thick square plate made of $CoF{e}_2{O}_4$ is considered under open-circuit magnetic boundary conditions. In this configuration, the loading conditions are prescribed by imposing the transverse magnetic induction components ${\mathcal{B}}_{z_t}$ and ${\mathcal{B}}_{z_b}$ at the top and bottom surfaces, respectively. The corresponding geometrical parameters, half-wave numbers, and prescribed loading boundary conditions are reported in the third column of Table 2, whereas the material properties of $CoF{e}_2{O}_4$ are summarized in Table 1. As reference solution, the three-dimensional magneto-elastic model developed by Heyliger and Pan [59] is adopted.

The comparison is carried out at different through-the-thickness locations $\tilde{z}/h$. As shown in Table 7, the proposed 3D-u-$\theta$-$\psi$ model exhibits an excellent agreement with the reference solution in terms of in-plane displacement, transverse displacement, and magnetic potential. Since thermal effects are neglected in this assessment, the obtained results provide an independent validation of the magneto-elastic coupling together with the capability of the proposed formulation to accurately reproduce the material layer and thickness effects.

The fourth preliminary assessment (PA4) considers a two-layered simply-supported thick plate with a Cadmium Selenide/PZT-5A configuration. The corresponding geometrical parameters, half-wave numbers, and prescribed loading boundary conditions are reported in the fourth column of Table 2, whereas the material properties of Cadmium Selenide and PZT-5A are summarized in Table 1. As reference solution, the three-dimensional Sampling Surface (SaS) method proposed by Kulikov and Plotnikova [60] is employed. The comparison is carried out at different through-the-thickness locations $\tilde{z}/h$ for the following dimensionless variables:

$$\begin{array}{cc}\hfill & \overline{u}={\displaystyle \frac{10u}{a{\alpha }_0{\mathsf{\Theta }}_t}};\overline{w}={\displaystyle \frac{100w}{\frac{a}{h}a{\alpha }_0{\mathsf{\Theta }}_t}};\overline{\theta }={\displaystyle \frac{\theta }{{\mathsf{\Theta }}_t}};{\overline{\sigma }}_{\alpha \alpha }={\displaystyle \frac{10{\sigma }_{\alpha \alpha }}{E_1{\alpha }_0{\mathsf{\Theta }}_t}};{\overline{\sigma }}_{\alpha \beta }={\displaystyle \frac{10{\sigma }_{\alpha \beta }}{E_1{\alpha }_0{\mathsf{\Theta }}_t}};\hfill \end{array}$$

(81)

$$\begin{array}{cc}& {\overline{\sigma }}_{\alpha z}={\displaystyle \frac{10\frac{a}{h}{\sigma }_{\alpha z}}{E_1{\alpha }_0{\mathsf{\Theta }}_t}};{\overline{\sigma }}_{zz}={\displaystyle \frac{10{\left(\frac{a}{h}\right)}^2{\sigma }_{zz}}{E_1{\alpha }_0{\mathsf{\Theta }}_t}};{\overline{q}}_z={\displaystyle \frac{10a{q}_z}{\frac{a}{h}{k}_0{\mathsf{\Theta }}_t}}.\hfill \end{array}$$

(82)

${\mathsf{\Theta }}_t$ denotes the prescribed over-temperature at the top surface of the structure, $k_0=9W/mK$, ${\alpha }_0=4.396·{10}^{-6}1/K$, and $E_1$ is the Young’s modulus of Cadmium Selenide expressed in Pascal. As shown in Table 8, the proposed 3D-u-$\theta$-$\psi$ model exhibits an excellent agreement with the three-dimensional reference solution for all the normalized variables considered. These results provide a further validation of the proposed formulation, confirming its capability to accurately capture the thermo-elastic coupling together with the material layer and thickness effects.

Thanks to these four preliminary assessments, all the key features involved in the thermo-magneto-elastic response of both flat and curved multilayered structures have been successfully validated. In particular, the proposed 3D-u-$\theta$-$\psi$ model accurately captures the material layer effect, the thickness effect, the thermo-magneto-elastic coupling, and the influence of structural curvature. The validation procedure has been carried out by adopting $M=300$ mathematical layers and an exponential matrix order equal to $N=3$. The convergence studies reported in Table 3, Table 4 and Table 5 demonstrate that these values provide fully converged solutions for both flat and curved geometries. Consequently, the same numerical parameters ($M=300$ and $N=3$) are employed in the four new benchmark cases presented in the following subsection.

3.2. New Results

This subsection presents four new benchmark problems involving a sandwich plate, a cylinder, a cylindrical panel, and a spherical shell. According to the convergence analyses discussed in the previous subsection, all the results are obtained by adopting an exponential matrix order equal to $N=3$ and $M=300$ mathematical layers.

All the investigated structures consider an Adaptive Wood/Soft Adaptive Wood/Adaptive Wood sandwich configuration. Each outer Adaptive Wood layer has a thickness equal to $0.1h$, whereas the Soft Adaptive Wood core has a thickness of $0.8h$, with h denoting the total thickness of the multilayered structure. The corresponding elastic, magnetic, and thermal material properties are listed in Table 1. Unlike the preliminary assessments, the magnetic permeability coefficients ${\mu }_{11}$ and ${\mu }_{22}$ are assumed to be positive, consistently with the physical interpretation proposed by Pan [61].

Two different loading boundary configurations are considered. Both configurations prescribe the over-temperature at the outer surfaces, whereas the magnetic loading is assigned either in terms of magnetic potential or transverse magnetic induction. The two loading boundary conditions (LB) are defined as follows:

$$LB1:{\overline{\mathsf{\Theta }}}_t=1.5K,{\overline{\mathsf{\Theta }}}_b=0K;{\overline{\mathsf{\Psi }}}_t=10A,{\overline{\mathsf{\Psi }}}_b=0A,$$

(83)

$$LB2:{\overline{\mathsf{\Theta }}}_t=1.5K,{\overline{\mathsf{\Theta }}}_b=0K;{\overline{\mathcal{B}}}_{z_t}=1T,{\overline{\mathcal{B}}}_{z_b}=0T.$$

(84)

Subscripts b and t denote the bottom and top surfaces of the multilayered structure, respectively.

For each benchmark problem, eight representative variables are reported in tabular form at selected through-the-thickness locations $\tilde{z}/h$ and for different thickness ratios. In addition, through-the-thickness distributions are presented for a representative thickness ratio under both loading boundary configurations in order to illustrate the main characteristics of the thermo-magneto-elastic response.

In the first new case (NC1), a sandwich square plate is investigated. The corresponding geometrical parameters and half-wave numbers are reported in the first column of Table 9. Table 10 and Table 11 collect representative results for the two loading boundary configurations by considering the variables u (or v), w, $\psi$, $\theta$, ${\mathcal{B}}_z$, $q_z$, ${\sigma }_{\alpha \beta }$, ${\sigma }_{\beta \beta }$, and ${\sigma }_{zz}$ at selected $\tilde{z}/h$ locations and for different thickness ratios.

As the slenderness ratio $a/h$ increases, both the transverse displacement w and the over-temperature $\theta$ increase because of the progressive reduction in the structural stiffness. The magnetic potential $\psi$ exhibits a sign change for the LB1 configuration at $\tilde{z}/h=0.25$, whereas it remains negative for the LB2 configuration at $\tilde{z}/h=0.75$, owing to the prescribed transverse magnetic induction that maintains the magnetic field oriented along the negative direction. A similar behavior is observed for the transverse magnetic induction ${\mathcal{B}}_z$, which changes from positive to negative values under the LB1 configuration ($\tilde{z}/h=1.00$), while remaining positive for the LB2 configuration ($\tilde{z}/h=0.75$). Conversely, the transverse heat flux $q_z$ assumes identical values for both loading configurations, indicating that the imposed magnetic boundary conditions do not influence the heat transfer through the thickness.

Figure 2 and Figure 3 report the through-the-thickness distributions of the selected variables for representative thickness ratios ($a/h=4$ for LB1 and $a/h=10$ for LB2). The zigzag behavior associated with the multilayered configuration is clearly visible through the slope changes occurring at each physical interface. Moreover, all the proposed variables satisfy the interlaminar continuity conditions, as no discontinuities are observed across the interfaces.

The prescribed loading boundary conditions are correctly recovered at the outer surfaces. In the LB1 configuration, the over-temperature varies from $1.5K$ at the top surface to $0K$ at the bottom surface, while the magnetic potential decreases from $10A$ to $0A$. In the LB2 configuration, the same thermal distribution is obtained, whereas the transverse magnetic induction varies from $1T$ at the top surface to $0T$ at the bottom surface. Since no external mechanical pressure is applied, the transverse normal stress ${\sigma }_{zz}$ is equal to zero at both outer surfaces for the two loading configurations.

Overall, the proposed benchmark clearly highlights the capability of the present formulation to accurately capture the thermo-magneto-elastic coupling together with the material layer effect, the thickness effect, the zigzag behavior, and the correct fulfillment of both boundary and interlaminar continuity conditions.

In the second new case (NC2), a sandwich cylinder is considered. The corresponding geometrical parameters and half-wave numbers are reported in the second column of Table 9. Table 12 presents the variables v, w, $\psi$, $\theta$, ${\mathcal{B}}_z$, $q_z$, ${\sigma }_{\beta \beta }$, and ${\sigma }_{zz}$ at selected $\tilde{z}/h$ locations for thickness ratios ranging from $R_{\alpha }/h=2$ to $R_{\alpha }/h=100$ under the LB1 loading configuration. Similarly, Table 13 reports u, w, $\psi$, $\theta$, ${\mathcal{B}}_z$, $q_z$, ${\sigma }_{\alpha \alpha }$, and ${\sigma }_{zz}$ for the LB2 configuration.

For both loading configurations, the variation of the transverse displacement w with the thickness ratio is not monotonic because of the curvature effect associated with the cylindrical geometry. The magnetic loading strongly influences the magnetic variables. Under the LB1 configuration, the magnetic potential $\psi$ changes sign from negative values for thick cylinders to positive values for thin cylinders ($\tilde{z}/h=0.25$), whereas it remains negative for all thickness ratios under the LB2 configuration ($\tilde{z}/h=0.75$). This behavior is associated with the different magnetic boundary conditions: under LB1 the magnetic potential reverses its distribution through the thickness as the cylinder becomes thinner, whereas under LB2 its direction remains unchanged. Conversely, the transverse magnetic induction ${\mathcal{B}}_z$ remains positive for both loading configurations. It decreases with increasing $R_{\alpha }/h$ under LB1 ($\tilde{z}/h=1.00$), whereas it increases under LB2 ($\tilde{z}/h=0.25$). The transverse heat flux $q_z$ assumes identical values for the two loading configurations, confirming that the prescribed magnetic boundary conditions do not affect the heat transfer through the thickness.

Figure 4 and Figure 5 show the through-the-thickness distributions of the selected variables for representative thickness ratios. In both cases, the zigzag behavior is clearly visible through the slope discontinuities at the material interfaces, while all variables remain continuous across the interfaces, confirming the correct enforcement of the interlaminar continuity conditions.

The prescribed loading boundary conditions are correctly satisfied at the outer surfaces. In Figure 4, this is demonstrated by the magnetic potential $\psi$ and over-temperature $\theta$ distributions, whereas in Figure 5${\mathcal{B}}_z$ and $\theta$ recover the imposed boundary values. Since no external mechanical pressure is applied, the transverse normal stress ${\sigma }_{zz}$ vanishes at both outer surfaces of the cylinder. Overall, the proposed benchmark confirms the capability of the present formulation to accurately reproduce the thermo-magneto-elastic coupling together with the curvature effect, the material layer effect, the thickness effect, the zigzag behavior, and the interlaminar continuity conditions.

The third new case (NC3) is devoted to a sandwich cylindrical panel. The corresponding geometrical parameters and half-wave numbers are reported in the third column of Table 9. Table 14 and Table 15 collect representative results for the two loading boundary configurations by considering the variables u, w, $\psi$, $\theta$, ${\mathcal{B}}_z$, $q_z$, ${\sigma }_{\alpha \alpha }$, and ${\sigma }_{zz}$ at selected $\tilde{z}/h$ locations and for different thickness ratios.

The influence of the shell curvature is evident, as several variables exhibit a non-monotonic variation with increasing $R_{\alpha }/h$. Moreover, the two loading boundary configurations produce significantly different responses for the transverse displacement w and the magnetic potential $\psi$, particularly for slender cylindrical panels. These differences arise from the combined effects of structural curvature, the progressive reduction in stiffness as the panel becomes thinner, and the different magnetic boundary conditions prescribed in the LB1 and LB2 configurations.

Figure 6 and Figure 7 present the through-the-thickness distributions of the selected variables for representative thickness ratios. In both loading configurations, the zigzag behavior is clearly identified through the slope changes occurring at each material interface, while the continuity of all variables confirms the correct enforcement of the interlaminar continuity conditions.

The prescribed loading boundary conditions are correctly recovered at the outer surfaces. In Figure 6, the imposed values are satisfied by the over-temperature $\theta$ and the magnetic potential $\psi$, whereas Figure 7 shows the correct fulfillment of the boundary conditions for $\theta$ and the transverse magnetic induction ${\mathcal{B}}_z$. As no external mechanical pressure is applied, the transverse normal stress ${\sigma }_{zz}$ vanishes at both outer surfaces. Overall, the proposed benchmark confirms the capability of the present formulation to accurately reproduce the thermo-magneto-elastic coupling together with the curvature effect, the material layer effect, the thickness effect, the zigzag behavior, and the interlaminar continuity conditions.

In the fourth new case (NC4), a sandwich spherical shell is investigated. The corresponding geometrical parameters and half-wave numbers are reported in the fourth column of Table 9. Table 16 and Table 17 collect representative results for the two loading boundary configurations by considering the in-plane and transverse displacements, magnetic potential $\psi$, over-temperature $\theta$, transverse magnetic induction ${\mathcal{B}}_z$, transverse heat flux $q_z$, in-plane stress components, and transverse normal stress ${\sigma }_{zz}$ at selected $\tilde{z}/h$ locations and for different thickness ratios.

Similarly to the cylindrical panel discussed in NC3, the magnetic variables are strongly influenced by the prescribed loading boundary conditions. The magnetic potential $\psi$ and the transverse magnetic induction ${\mathcal{B}}_z$ exhibit different behaviors as the structure changes from thick to thin spherical shells because the magnetic field distribution varies according to the adopted magnetic boundary conditions. In particular, under the LB1 configuration, the magnetic potential changes sign from negative values for thick shells to positive values for thin shells, whereas it remains negative for all thickness ratios under the LB2 configuration. Conversely, the transverse magnetic induction ${\mathcal{B}}_z$ remains positive for both loading configurations. The in-plane displacement v, transverse displacement w, over-temperature $\theta$, and transverse heat flux $q_z$ exhibit the same overall trend, with increasing values as the shell becomes thinner owing to the progressive reduction in the structural stiffness.

Figure 8 and Figure 9 present the through-the-thickness distributions of w, $\psi$, $\theta$, $q_z$, ${\sigma }_{zz}$, and ${\mathcal{B}}_z$ for representative thickness ratios under the LB1 and LB2 loading configurations, respectively. In both cases, the zigzag behavior is clearly identified through the slope changes occurring at the material interfaces, while the continuity of all variables confirms the correct enforcement of the interlaminar continuity conditions.

The prescribed loading boundary conditions are correctly recovered at the outer surfaces for both loading configurations. Since no external mechanical pressure is applied, the transverse normal stress ${\sigma }_{zz}$ vanishes at the top and bottom surfaces of the spherical shell. Overall, the proposed benchmark confirms the capability of the present formulation to accurately reproduce the thermo-magneto-elastic coupling together with the curvature effect, the material layer effect, the thickness effect, the zigzag behavior, and the interlaminar continuity conditions.

4. Conclusions

The present study has introduced a three-dimensional layer-wise coupled thermo-magneto-elastic shell model for the static analysis of multilayered flat and curved smart structures. The governing formulation consists of the three-dimensional elastic equilibrium equations, the magnetic induction divergence equation, and the Fourier heat conduction equation, all fully coupled within a unified multifield framework. The analytical solution combines Navier harmonic expansions in the in-plane directions with the exponential matrix method along the thickness coordinate. By adopting a mixed orthogonal curvilinear reference system, the proposed formulation provides a unified framework for the analysis of simply supported plates, cylinders, cylindrical panels, and spherical shells through the appropriate definition of the curvature radii. The layer-wise strategy is implemented by enforcing interlaminar continuity conditions on displacements, magnetic potential, over-temperature, transverse shear and normal stresses, transverse magnetic induction, and transverse heat flux. Different thermal, magnetic, and mechanical loading configurations can be prescribed at the outer surfaces through over-temperature, magnetic potential, and transverse magnetic induction boundary conditions.

The numerical investigation has been divided into two complementary stages. First, the proposed formulation has been validated through comparisons with thermo-magneto-elastic solutions available in the literature. The obtained results demonstrate the capability of the model to accurately reproduce thermo-magneto-elastic coupling, curvature effects, and the influence of both the material sequence and the structural thickness. The validation campaign has been carried out by employing an exponential matrix order equal to $N=3$ and $M=300$ mathematical layers, whose adequacy has been confirmed through dedicated convergence analyses. Subsequently, four new benchmark problems involving simply supported sandwich plates, cylinders, cylindrical panels, and spherical shells have been presented. Thick and thin configurations have been investigated under two different loading boundary configurations. The structural response has been evaluated in terms of displacements, magnetic potential, over-temperature, magnetic induction, heat flux, and stress components by means of both tabulated values at selected through-the-thickness locations and graphical distributions along the thickness direction.

The proposed benchmark results clearly demonstrate the capability of the formulation to reproduce all the main physical features characterizing thermo-magneto-elastic multilayered structures. In particular, thermo-magneto-elastic coupling, curvature effects, zigzag behavior, material layer effects, thickness effects, and interlaminar continuity conditions are accurately captured for both flat and curved geometries. Furthermore, the distributions of the elastic, thermal, and magnetic fields are shown to be strongly affected by the prescribed loading boundary conditions, emphasizing the importance of fully coupled three-dimensional analyses for the design and optimization of smart aerospace structures.

The benchmark solutions generated in this work provide a reliable reference database for researchers developing advanced multifield formulations for smart multilayered structures. Owing to its fully three-dimensional analytical nature, the proposed model can also be employed as a rigorous validation tool for future two-dimensional and three-dimensional analytical or numerical formulations involving coupled thermo-magneto-elastic phenomena.