Bending–Extension Coupling at Restrained Ends and the Nonlinear Transient Response of Porous FGM Sandwich Beams in Thermal Environments

1. Introduction

Functionally graded materials (FGMs) vary the volume fractions of ceramic and metal smoothly through one or more directions. This gradation avoids the stress concentration and delamination found in conventional laminated composites, while combining the heat resistance of ceramics with the toughness of metals. FGMs are therefore used in thermal barrier components, aerospace structures and nuclear plants. Sandwich construction adds a high bending-stiffness-to-weight ratio, and FGM sandwich beams and plates have been studied extensively for buckling, vibration and dynamic response [1,2,3,4].

Pores may form inside FGMs during fabrication, for example during sintering. Their effect has been modelled with modified rules of mixture for even and uneven distributions [5]. Porosity has since been included in the linear vibration of FGM beams [6,7], the thermal postbuckling and nonlinear vibration of porous FGM beams with elastically restrained ends [8,9], and the nonlinear dynamic response of porous FGM plates on elastic foundations [10].

The in-plane restraint of beam ends strongly affects thermally loaded structures. When the ends cannot move freely, a uniform temperature rise produces axial compressive forces that reduce bending stiffness. This effect has been modelled through an elastic restraint parameter [11], and a growing body of semi-analytical work has applied it to the buckling, vibration and nonlinear transient response of porous FGM beams, plates, shells and sandwich beams in thermal environments [8,9,12,13,14,15]. These studies have clarified how porosity, gradation, temperature, foundation support, initial imperfection and the degree of end restraint govern the global response; for example, end restraint deepens the dynamic deflection at elevated temperature, and initial imperfection can be beneficial at room temperature but detrimental when heated.

Bending–extension coupling gives asymmetric FGM beams a behaviour with no counterpart in homogeneous beams. Their nonlinear vibration differs for positive and negative amplitudes [16,17]. Under uniform heating, simply supported FGM beams bend instead of buckling by bifurcation, whereas clamped beams buckle [18,19,20,21]. Formulations referred to the physical neutral surface remove the coupling from the constitutive equations [22,23], and for immovable simple supports the mid-plane and neutral-surface formulations give different responses because each fixes the structure at a different plane [24]. None of these effects is new, and the present paper does not claim them as such.

Despite this progress, the attention of transient studies has been directed mainly at the global response, and less at a boundary-level detail with a clear physical effect: the moment that bending–extension coupling and the thermal moment create at restrained ends. A widely used semi-analytical approach, adopted in many studies of this class of structures (see, e.g., [8,9,12,13,14,15]), expands the deflection in a single sine mode and projects it by the Galerkin method onto the strong-form equations of motion. This choice is efficient and accurate for symmetric sections. For asymmetric sections, however, the sine mode has zero curvature at the ends and cannot satisfy the end-moment condition, and the constant end moment vanishes when the strong-form equations are differentiated. The beam then behaves as if it were hinged at its neutral surface and free of thermal bending, whatever support conditions are assumed. To the authors’ knowledge, the consequences of this common simplification for the nonlinear transient response of asymmetric porous FGM sandwich beams have not been quantified, and this is the gap the present paper addresses.

The present work uses a refined higher-order shear deformation theory (HSDT), which gives a parabolic shear strain that vanishes at the surfaces and needs no shear correction factor [25,26]. A preliminary check with a quasi-3D theory showed that thickness stretching changes the peak dynamic deflection of Si₃N₄/SUS304 beams by less than 2% down to L/h = 5, so it is not pursued here.

This paper makes several key contributions: it develops a weak-form, multi-term Galerkin solution based on the HSDT for porous FGM sandwich beams with elastically restrained ends, clearly revealing what the strong-form approximation neglects; it validates this approach using an exact linear solution and an independent nonlinear finite element model; it quantifies the strong-form deviations in frequency and peak deflection by separating the effects of end-reaction coupling stiffness and thermal bending; it properly accounts for the thermal pre-load as an initial equilibrium state to show how the loading sequence affects the response; and it conducts comprehensive parametric studies to determine which previous strong-form trends remain valid for asymmetric sandwich beams.

2. Porous FGM Sandwich Beams

2.1. Geometry

Consider a sandwich beam of length L, unit width b = 1 and total thickness h, resting on a Winkler–Pasternak foundation. The origin of the Cartesian system (x, y, z) lies on the mid-plane at the left end; x runs along the beam and z through the thickness, positive downward. The beam has two face sheets of equal thickness h_f and a core of thickness h − 2h_f, perfectly bonded. The layer boundaries are shown in Figure 1 and given by

$$h_0=-{\displaystyle \frac{h}{2}}, h_1=-{\displaystyle \frac{h}{2}}+h_f, h_2={\displaystyle \frac{h}{2}}-h_f, h_3={\displaystyle \frac{h}{2}}$$

(1)

2.2. Type A: FGM Face Sheets and Metal Core

The ceramic volume fraction decreases from each outer surface to its interface, and the core is pure metal. The beam is symmetric about the mid-plane:

$$V_c\left(z\right)=\left\{\begin{array}{ll}{\left({\displaystyle \frac{z-h_1}{h_0-h_1}}\right)}^N,& h_0\le z\le h_1\\ 0,& h_1\le z\le h_2\\ {\left({\displaystyle \frac{z-h_2}{h_3-h_2}}\right)}^N,& h_2\le z\le h_3\end{array}\right.$$

(2)

2.3. Type B: Metal Top Face, FGM Core, Ceramic Bottom Face

The top face is pure metal, the bottom face is pure ceramic, and the ceramic fraction in the core rises from the upper to the lower interface. The beam is asymmetric about the mid-plane:

$$V_c\left(z\right)=\left\{\begin{array}{ll}0,& h_0\le z\le h_1\\ {\left({\displaystyle \frac{z-h_1}{h_2-h_1}}\right)}^N,& h_1\le z\le h_2\\ 1,& h_2\le z\le h_3\end{array}\right.$$

(3)

In both types the metal fraction is V_m = 1 − V_c, and N ≥ 0 is the volume fraction index.

2.4. Porosity Models

The effective property P_eff (Young’s modulus E, thermal expansion coefficient α, mass density ρ) of the porous FGM follows the modified rule of mixture [5], with total pore volume fraction ξ (0 ≤ ξ ≪ 1):

$$P_{\mathrm{eff}}\left(z,T\right)=\left(P_c\left(T\right)-P_m\left(T\right)\right)V_c\left(z\right)+P_m\left(T\right)-{\displaystyle \frac{\xi }{2}}\left(P_c\left(T\right)+P_m\left(T\right)\right)\hspace{0.17em}\psi \left(z\right)$$

(4)

For even porosity ψ(z) = 1 in every FGM layer. For uneven porosity:

• Type A: ψ(z) = (z − h₀)/(h₁ − h₀) in the top face, ψ(z) = 0 in the metal core, and ψ(z) = (z − h₃)/(h₂ − h₃) in the bottom face. Pores are densest at the interfaces.
• Type B: ψ(z) = 1 − 2|z|/(h₂ − h₁) in the core and ψ(z) = 0 in the homogeneous faces. Pores are densest at the mid-plane.

The porosity distribution is written ψ(z) because f(z) is reserved for the shear shape function of the HSDT in Section 3.

2.5. Temperature-Dependent Properties

The beam is heated uniformly from T₀ = 300 K to T, so ΔT = T − T₀. The properties of each constituent follow the Touloukian form

P(T) = P₀(P₋₁T⁻¹ + 1 + P₁T + P₂T² + P₃T³)

(5)

with the coefficients of Table 1 [27]. The mass densities are ρ_c = 2370 kg/m³ for Si₃N₄ and ρ_m = 8166 kg/m³ for SUS304. Poisson’s ratio is taken as constant, ν = 0.28, a value commonly adopted for this material pair [15], which also allows direct comparison with published FSDT results. In the HSDT used here, ν enters only through the shear modulus.

3. HSDT Formulation

3.1. Displacement Field

The refined higher-order shear deformation theory splits the deflection into bending and shear parts [25,26]:

$$\left.\begin{array}{c}\bar{u}\left(x,z,t\right)=u\left(x,t\right)-z\hspace{0.17em}{\displaystyle \frac{\partial w_b}{\partial x}}-f\left(z\right)\hspace{0.17em}{\displaystyle \frac{\partial w_s}{\partial x}},\\ \bar{w}\left(x,z,t\right)=w_b\left(x,t\right)+w_s\left(x,t\right)==w\left(x,t\right)\end{array}\right\}$$

(6)

$$f\left(z\right)=-{\displaystyle \frac{z}{4}}+{\displaystyle \frac{5z^3}{3h^2}}, g\left(z\right)=1-f\mathrm{\text{'}}\left(z\right)={\displaystyle \frac{5}{4}}-{\displaystyle \frac{5z^2}{h^2}}$$

(7)

The transverse shear strain is parabolic through the thickness and vanishes at z = ±h/2, so no shear correction factor is needed. The initial geometric imperfection is w^*(x).

3.2. Strains3.3. Stresses and Resultants

$${\epsilon }_{\mathit{xx}}={\epsilon }_0+z\hspace{0.17em}{\kappa }_b+f\left(z\right)\hspace{0.17em}{\kappa }_s, {\gamma }_{\mathit{xz}}=g\left(z\right)\hspace{0.17em}{\displaystyle \frac{\partial w_s}{\partial x}}$$

(8)

$${\epsilon }_0={\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial x}}\right)}^2+{\displaystyle \frac{\partial w}{\partial x}}{\displaystyle \frac{\partial w^{*}}{\partial x}}, {\kappa }_b=-{\displaystyle \frac{{\partial }^2{w}_b}{\partial x^2}}, {\kappa }_s=-{\displaystyle \frac{{\partial }^2{w}_s}{\partial x^2}}$$

(9)

The uniform temperature rise ΔT is applied before the load and acts as a pre-existing axial thermal stress:

$${\sigma }_{\mathit{xx}}=E\left(z,T\right)\left({\epsilon }_{\mathit{xx}}-\alpha \left(z,T\right)\hspace{0.17em}\mathit{\Delta T}\right), {\tau }_{\mathit{xz}}={\displaystyle \frac{E\left(z,T\right)}{2\left(1+\nu \right)}}\hspace{0.17em}{\gamma }_{\mathit{xz}}$$

(10)

The resultants N_x, M_x^b, M_x^s (weights 1, z, f) and Q_xz (weight g) are

$$\begin{array}{ll}{N}_x& =A_{11}{\epsilon }_0+B_{11}{\kappa }_b+B_{11}^s{\kappa }_s-N^T\\ M_x^b& =B_{11}{\epsilon }_0+D_{11}{\kappa }_b+D_{11}^s{\kappa }_s-M^T\\ M_x^s& =B_{11}^s{\epsilon }_0+D_{11}^s{\kappa }_b+H_{11}^s{\kappa }_s-P^T\\ Q_{\mathit{xz}}& =A_{55}^s\hspace{0.17em}{\displaystyle \frac{\partial w_s}{\partial x}}\end{array}$$

(11)

With

$$\left.\begin{array}{c}\left(A_{11},B_{11},D_{11},B_{11}^s,D_{11}^s,H_{11}^s\right)={\int }_{-{\displaystyle \frac{h}{2}}}^{{\displaystyle \frac{h}{2}}}E\left(1,z,z^2,f,\mathit{zf},f^2\right)\mathit{dz}\\ A_{55}^s={\int }_{-{\displaystyle \frac{h}{2}}}^{{\displaystyle \frac{h}{2}}}{\displaystyle \frac{E}{2\left(1+\nu \right)}}\hspace{0.17em}{g}^2\hspace{0.17em}\mathit{dz}\\ \left(N^T,M^T,P^T\right)={\int }_{-h/2}^{h/2}E\hspace{0.17em}\alpha \hspace{0.17em}\mathit{\Delta T}\left(1,z,f\right)\mathit{dz}, I_0={\int }_{-h/2}^{h/2}\rho \hspace{0.17em}\mathit{dz}\end{array}\right\}$$

(12)

For the symmetric type A beam B₁₁ = B^s₁₁ = M^T = P^T = 0. For the asymmetric type B beam all coefficients are non-zero, which is the source of the end-moment effect discussed in Section 3.6.

3.4. Hamilton’s Principle

Axial and rotary inertia are neglected, as is common in this class of analyses [9,15]; neglecting axial inertia is the usual static condensation of hinged-beam models [28]. With the foundation on the bottom surface and the load on the top surface, both deflecting by w, the weak statement of motion is

$$\begin{array}{rr}& {\int }_{t_1}^{t_2}{\int }_0^L[I_0\hspace{0.17em}\dot{w}\hspace{0.17em}\delta \dot{w}-N_x\left(\delta u_{,x}+{\left(w+w^{*}\right)}_{,x}\hspace{0.17em}\delta w_{,x}\right)-M_x^b\hspace{0.17em}\delta {\kappa }_b-M_x^s\hspace{0.17em}\delta {\kappa }_s\\ & -Q_{\mathit{xz}}\hspace{0.17em}\delta w_{s,x}-\left(k_{1}w\hspace{0.17em}\mathit{\delta w}+k_2{w}_{,x}\hspace{0.17em}\delta w_{,x}\right)+q\hspace{0.17em}\mathit{\delta w}]\hspace{0.17em}\mathit{dx}\hspace{0.17em}\mathit{dt}=0\end{array}$$

(13)

This weak statement is the starting point of the solution in Section 4.

3.5. Equations of Motion

Integrating by parts gives the strong form:

$${\displaystyle \frac{\partial N_x}{\partial x}}=0$$

(14)

$${\displaystyle \frac{{\partial }^2{M}_x^b}{\partial x^2}}+N_{x0}{\displaystyle \frac{{\partial }^2\left(w+w^{*}\right)}{\partial x^2}}+q-k_{1}w+k_2{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}=I_0\hspace{0.17em}\ddot{w}$$

(15)

$${\displaystyle \frac{{\partial }^2{M}_x^s}{\partial x^2}}+{\displaystyle \frac{\partial Q_{\mathit{xz}}}{\partial x}}+N_{x0}{\displaystyle \frac{{\partial }^2\left(w+w^{*}\right)}{\partial x^2}}+q-k_{1}w+k_2{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}=I_0\hspace{0.17em}\ddot{w}$$

(16)

Equation (14) gives N_x = N_x0, constant along the beam.

3.6. End Conditions and the End Moment

Both ends are hinged and elastically restrained against axial displacement by springs of stiffness c. The hinges and springs act at the height z = z_s: z_s = 0 places them at the mid-plane, as usually assumed in mid-plane formulations, and z_s = z₀ = B₁₁/A₁₁ places them at the physical neutral surface [22,23,24]. Moment equilibrium about the hinge and the spring law give

$$w_b=w_s=0, M_x^b-z_s{N}_x=0, M_x^s-f\left(z_s\right)N_x=0, N_x=N_{x0} \mathrm{at}x=0,\hspace{0.17em}L$$

(17)

$$\begin{array}{ll}{N}_{x0}& =-{\displaystyle \frac{c}{L}}{\int }_0^L{\displaystyle \frac{\partial \bar{u}\left(x,z_s,t\right)}{\partial x}}\hspace{0.17em}\mathit{dx}, \mathrm{which\; gives}\\ N_{x0}& =-\lambda N^T+{\displaystyle \frac{\lambda }{L}}{\int }_0^L\left[\left(B_{11}-A_{11}{z}_s\right){\kappa }_b+\left(B_{11}^s-A_{11}f\left(z_s\right)\right){\kappa }_s\right]\mathit{dx}\\ & +{\displaystyle \frac{\lambda A_{11}}{L}}{\int }_0^L\left({\displaystyle \frac{1}{2}}{w}_{,x}^2+w_{,x}\hspace{0.17em}{w}_{,x}^{*}\right)\mathit{dx}\end{array}$$

(18)

with λ = c/(c + A₁₁): λ = 0 for movable, 0 < λ < 1 for partially movable and λ = 1 for immovable ends. Eliminating ε₀ with N_x = N_x0 gives the moment about the hinge as

$$M_x^b-z_s{N}_{x0}=\underset{{\widehat{M}}_0^b}{\underbrace{{\displaystyle \frac{B_{11}}{A_{11}}}\left(N_{x0}+N^T\right)-M^T-z_s{N}_{x0}}}+{\tilde{D}}_{11}{\kappa }_b+{\tilde{D}}_{11}^s{\kappa }_s$$

(19)

The term M̂^b₀ is constant along the beam. In an asymmetric beam hinged at the mid-plane it has two parts: the end reaction acts at the mid-plane while the neutral axis lies at the offset z₀, and the thermal moment M^T is non-zero. With hinges at the neutral surface the first part vanishes, but M̂^b₀ = −(M^T − z₀N^T), the thermal moment about the neutral axis, remains. In either case the end condition (17) requires non-zero curvature at the ends, which a sine mode cannot provide.

If Galerkin’s method is applied to the strong form with a one-term sine expansion, as in the approach discussed in Section 1, the constant M̂^b₀ disappears when differentiated twice and its effect is lost. If it is applied to the weak form (13) with the spring acting at z_s, M̂^b₀ does work on the end rotations of each mode and is retained. Section 5 confirms by an exact solution and by finite elements that the weak form is correct.

4. Solution Procedure

The deflections are expanded in n symmetric sine modes, which satisfy the essential end conditions:

$$w_b=\sum _i{W}_{\mathit{bi}}\left(t\right)\sin {\beta }_{i}x, w_s=\sum _i{W}_{\mathit{si}}\left(t\right)\sin {\beta }_{i}x, w^{*}=\mathit{\mu h}\sin {\beta }_{1}x, {\beta }_i={\displaystyle \frac{\mathit{i\pi }}{L}},$$

i = 1,3,…,2n − 1

(20)

Even modes are not excited by a uniform load and a symmetric imperfection, so they are omitted. The modal deflection is W_i = W_bi + W_si, and the mid-span deflection is w(L/2) = Σ_iW_isin(iπ/2). The one-term case (n = 1) corresponds to the solution form commonly used in the literature [9,15]. Unless stated otherwise, z_s = 0.

With N_x = N_x0, the membrane strain is ε₀ = (N_x0 + N^T − B₁₁κ_b − B^s₁₁κ_s)/A₁₁, so the moments about the hinge are

$$M_x^b-z_s{N}_{x0}={\widehat{M}}_0^b+{\tilde{D}}_{11}{\kappa }_b+{\tilde{D}}_{11}^s{\kappa }_s, M_x^s-f\left(z_s\right)N_{x0}={\widehat{M}}_0^s+{\tilde{D}}_{11}^s{\kappa }_b+{\tilde{H}}_{11}^s{\kappa }_s$$

(21)

$${\widehat{M}}_0^b={\displaystyle \frac{B_{11}}{A_{11}}}\left(N_{x0}+N^T\right)-M^T-z_s{N}_{x0}, {\widehat{M}}_0^s={\displaystyle \frac{B_{11}^s}{A_{11}}}\left(N_{x0}+N^T\right)-P^T-f\left(z_s\right)N_{x0}$$

(22)

$${\tilde{D}}_{11}=D_{11}-{\displaystyle \frac{B_{11}^2}{A_{11}}}, {\tilde{D}}_{11}^s=D_{11}^s-{\displaystyle \frac{B_{11}{B}_{11}^s}{A_{11}}}, {\tilde{H}}_{11}^s=H_{11}^s-{\displaystyle \frac{{\left(B_{11}^s\right)}^2}{A_{11}}}$$

(23)

Substituting the trial functions into Eq. (18) gives

$$N_{x0}=-\lambda N^T+\lambda \sum _i{\displaystyle \frac{2}{\mathit{i\pi }}}\hspace{0.17em}{\beta }_i^2\left[\left(B_{11}-A_{11}{z}_s\right)W_{\mathit{bi}}+\left(B_{11}^s-A_{11}f\left(z_s\right)\right)W_{\mathit{si}}\right]+{\displaystyle \frac{\lambda A_{11}}{4}}\sum _i{\beta }_i^2\hspace{0.17em}{W}_i\left(W_i+2{\mu }_{i}h\right)$$

(24)

where μ₁ = μ and μ_i = 0 for i > 1. The modes are coupled only through N_x0.

The weak statement (13) is applied with each mode as test function, and every equation is divided by L/2. With Π_i = 4/(iπ):

$$I_0{\ddot{W}}_i+{\tilde{D}}_{11}{\beta }_i^4{W}_{\mathit{bi}}+{\tilde{D}}_{11}^s{\beta }_i^4{W}_{\mathit{si}}+{\Pi }_i{\beta }_i^2\hspace{0.17em}{\widehat{M}}_0^b+N_{x0}{\beta }_i^2\left(W_i+{\mu }_{i}h\right)+\left(k_1+k_2{\beta }_i^2\right)W_i={\Pi }_i\hspace{0.17em}q$$

(25)

$$I_0{\ddot{W}}_i+{\tilde{D}}_{11}^s{\beta }_i^4{W}_{\mathit{bi}}+\left({\tilde{H}}_{11}^s{\beta }_i^4+A_{55}^s{\beta }_i^2\right)W_{\mathit{si}}+{\Pi }_i{\beta }_i^2\hspace{0.17em}{\widehat{M}}_0^s+N_{x0}{\beta }_i^2\left(W_i+{\mu }_{i}h\right)+\left(k_1+k_2{\beta }_i^2\right)W_i={\Pi }_i\hspace{0.17em}q$$

(26)

The terms Π_iβ_i²M̂₀ are the work of the constant end moments on the end rotations. A strong-form projection drops them. Through N_x0, M̂^b₀ contains a term linear in the deflection, so it changes the linear stiffness of an immovable asymmetric beam, not only its nonlinear response. It also contains −M^T, which acts as an equivalent thermal load.

Because rotary inertia is neglected, Eqs. (25) and (26) have the same inertia term I₀Ẅ_i. Subtracting them gives n algebraic equations:

$$\left({\tilde{D}}_{11}^s-{\tilde{D}}_{11}\right){\beta }_i^4{W}_{\mathit{bi}}+\left[\left({\tilde{H}}_{11}^s-{\tilde{D}}_{11}^s\right){\beta }_i^4+A_{55}^s{\beta }_i^2\right]W_{\mathit{si}}+{\Pi }_i{\beta }_i^2\left({\widehat{M}}_0^s-{\widehat{M}}_0^b\right)=0$$

(27)

With W_bi = W_i − W_si, these equations are linear in the vector of W_si for given W_i, since N_x0 is linear in W_si. They are solved at every evaluation, which leaves n second-order equations I₀Ẅ_i = −F_bi(W), where F_bi is the non-inertial part of Eq. (25).

The temperature rise is applied before the load. When heating alone deflects the beam (type B, or μ ≠ 0 with λ > 0), the thermal equilibrium state is found by Newton–Raphson iteration with q = 0 and used as the initial state. The load is then applied suddenly at t = 0. Some studies apply heating and loading together from a flat beam at rest, W(0) = Ẇ(0) = 0; for comparison with such results, zero initial conditions are also used where stated.

4.1. Time Integration and Frequencies

The 2n first-order equations are integrated with the classical fourth-order Runge–Kutta scheme (RK4). The step is the smaller of 0.5/ω_max and T₁/400, where ω_max is the highest linear frequency in the expansion and T₁ the fundamental period. The results are cross-checked with an implicit Radau solver. Natural frequencies follow from the eigenvalues of the linearised system about the thermal equilibrium state. The number of terms n is fixed by the convergence study of Section 5.

5. Verification

The weak-form solution is verified against published results (Section 5.1 and Section 5.3), an exact solution (Section 5.2) and an independent nonlinear finite element model (FE)(Section 5.5). Single-layer FGM beams are obtained from type B with h_f = 0. Frequencies are non-dimensionalised as ω̄ = ω(L²/h)√(ρ_m/E_m), with metal properties at the current temperature.

5.1. Linear Frequencies Against Published Results

Table 2 compares the fundamental frequency of Al/Al₂O₃ beams (E_m = 70 GPa, ρ_m = 2702 kg/m³, E_c = 380 GPa, ρ_c = 3960 kg/m³, ν = 0.3; L/h = 20, N = 2, ΔT = 0). In the published studies these values are reported for immovable ends; as Section 5.2 shows, they coincide with the movable-end value. They are therefore compared with the present HSDT at λ = 0. The agreement between the results is remarkably high.

5.2. Exact Solution for Immovable Ends

To settle which treatment of the restrained end is correct, the linear free vibration of a hinged Euler–Bernoulli FGM beam is solved exactly. With u = 0 at both ends, the axial force is constant, N = −(B₁₁/L)[w′(L) − w′(0)], and the end condition M = 0 becomes

$${\tilde{D}}_{11}\hspace{0.17em}w⁗=\rho \hspace{0.17em}{\omega }^{2}w, w=0, w″=-{\displaystyle \frac{B_{11}^2}{A_{11}{\tilde{D}}_{11}\hspace{0.17em}L}}\left[w\mathrm{\text{'}}\left(L\right)-w\mathrm{\text{'}}\left(0\right)\right] \mathrm{at}x=0,L$$

(28)

The characteristic equation is solved for the lowest root. Setting B₁₁ = 0 in the end condition gives the movable-end value. The results are listed in Table 3.

Three observations follow:

The strong-form result is the same with and without end restraint. It lies within 0.3% of the exact movable value, the gap being the shear deformation that the exact Euler–Bernoulli solution ignores.

The weak-form result converges from above to the exact immovable value, as expected of a Ritz-type approximation. With five terms it is within 0.1% of the exact solution.

The strong-form deviation grows with porosity, from 10% to 19% for even pores, because pores lower A₁₁ more than they lower the coupling ratio B₁₁²/(A₁₁D̃₁₁).

With the hinges at the neutral surface (z_s = z₀), the weak form gives ω̄ = 3.8407 for ξ = 0, identical to the movable-end value, because the coupling term of Eq. (18) vanishes. Immovable-end frequencies obtained with one-term strong-form solutions therefore correspond to hinges at the neutral surface, or equivalently to movable ends, rather than to mid-plane hinges. This is a property of the approximation rather than of any particular calculation; it affects only asymmetric beams, and it also explains why independently obtained benchmark values agree closely with one another (Table 2).

5.3. Thermal Environment

Table 4 repeats the Si₃N₄/SUS304 case of Ebrahimi and Jafari [6] (L/h = 20, ΔT = 80 K, λ = 1) with temperature-dependent properties. The strong-form results agree with both references, which confirms the material data and the thermal model. The weak form diverges only marginally for graded beams (N > 0). The effect is smaller than for Al/Al₂O₃ because the stiffness contrast of Si₃N₄/SUS304 is lower, which reduces the coupling ratio.

5.4. Nonlinear Transient Response and Convergence

The nonlinear transient response is checked against published strong-form results for FGM sandwich beams [15] (L/h = 20, N = 2, h_f/h = 0.2, even porosity, λ = 1, T = 300 K, q = 10⁶ N/m², zero initial conditions) in Table 5. The present strong-form solution reproduces the published curves, which confirms that the two implementations are equivalent. For the symmetric type A beam both forms agree, as expected. For type B, the weak form slightly reduces the peak at room temperature.

Table 6 shows convergence in the number of terms for a type B beam (N = 2, h_f/h = 0.2, uneven, ξ = 0.2, λ = 1, zero initial conditions). Three terms give converged peaks, and one term remains extremely close to the converged value. The strong form overestimates the peak by a minor margin at 300 K and by 23% at 400 K, where the thermal end moment adds to the coupling effect. The RK4 solution with n = 3 agrees with an implicit Radau solution (with a very tight relative tolerance) within an imperceptible margin over 5 ms. All results in Section 6 use n = 3.

5.5. Independent Finite Element Verification

The exact solution of Section 5.2 is linear and isothermal. To verify the nonlinear transient response in a thermal environment, an independent model was built that does not use sine modes. The beam is discretised with 40 Euler–Bernoulli elements with Hermite cubic deflection. Because axial inertia is neglected, N_x is constant and the axial displacement is condensed exactly through Eq. (18). The end-moment conditions (17) are natural conditions of the finite element weak form and are satisfied by the solution itself. The thermal equilibrium state is found by Newton–Raphson iteration, and the response to the suddenly applied load is integrated with the Newmark average-acceleration method (Δt = 2 μs) with Newton iterations at each step. The peak deflection is unchanged to five digits between 10 and 80 elements.

For a like-for-like comparison, the Ritz solution is evaluated in the same Euler–Bernoulli limit by making the shear stiffness A^s₅₅ very large. Table 7 and Figure 2 show the results. The finite element and weak-form solutions exhibit an exceptionally close agreement in peak deflection and thermal pre-deflection, whereas the strong-form approximation departs by a notable margin for type B. For type A the three solutions coincide. The last column shows that shear deformation adds an imperceptible effect (less than a minor fraction) a L/h = 20.

6. Results and Discussion

Unless stated otherwise: Si₃N₄/SUS304, h = 0.02 m, L/h = 20, q = 10⁶ N/m², N = 2, h_f/h = 0.2, uneven porosity with ξ = 0.2, λ = 1, (K₁, K₂) = (100, 10), μ = 0 and T = 400 K, where K₁ = 12k₁L⁴/(E_m⁰h³) and K₂ = 12k₂L²/(E_m⁰h³) as in [15]. The foundation keeps the beams below thermal buckling up to 700 K; without it, type A buckles thermally near 450 K. The weak form uses n = 3 terms and starts from the thermal equilibrium state (Section 4.6). The strong form denotes the one-term strong-form projection [9,15], whose equilibrium state is flat when μ = 0. The hinges are at the mid-plane except in Section 6.3. The peak is the maximum mid-span deflection W/h over 6 ms, measured from the flat beam.

6.1. Temperature

For the symmetric type A beam the two forms coincide at every temperature (Figure 3a), because M̂^b₀ = M̂^s₀ = 0. For type B they separate as the temperature rises (Figure 3b). The strong form overestimates the peak by by a minor margin at 300 K, and by a progressively larger and notable margin at 400 K, 500 K and 600 K. The gap narrows to a lower margin at 700 K. Near the thermal stability limit the load-induced amplitude of the pre-bent beam grows faster than that of the strong form (0.49 against 0.36), which partly offsets the pre-deflection. In the weak form, heating alone bends the type B beam upward, by 0.009h at 400 K, 0.024h at 500 K and 0.18h at 700 K. The strong form has no such pre-deflection.

6.2. End Restraint and the Source of the Difference

At 300 K the restraint hardly matters, because no thermal force develops (Figure 4a). At 500 K the peak grows with λ in both forms, though the strong-form results remain higher across all values of λ, including λ = 0 (0.076 against 0.036). With movable ends N_x0 = 0, so the remaining end moment is the thermal moment alone: a free asymmetric beam bends under uniform heating, as a bimetal strip does.

Figure 4b separates the two parts of the weak-form response. The load-induced amplitude (peak minus W₀) exhibits an exceptionally close alignment with the strong-form predictions across all values of λ. Therefore, the total discrepancy observed in Figure 4a is entirely attributed to the thermal pre-deflection W₀, which decreases in magnitude from 0.040h at λ = 0 to 0.024h at λ = 1 as the restraint resists bending.

For Si₃N₄/SUS304 beams the omitted end moment thus acts mainly through the thermal moment. Its linear stiffening effect is small here, consistent with the minimal differences shown in in Table 4, whereas it dominates for the higher-contrast Al/Al₂O₃ beams of Table 3. Which mechanism matters depends on the stiffness contrast and on the temperature.

6.3. Support Location

Figure 5 repeats the temperature study of type B with the hinges at the neutral surface (z_s = z₀). The thermal moment about the neutral axis is larger than the end moment of the mid-plane case, and the upward pre-deflection is about three times larger: −0.024h at 400 K, −0.069h at 500 K and −0.35h at 700 K. The total peak therefore falls with temperature, to 0.069 at 400 K and 0.059 at 500 K, while the strong form predicts 0.092 and 0.125. The load-induced amplitude remains within an exceptionally close margin of the strong-form value up to 550 K and deviates near the stability limit (Figure 5b).

Therefore, the strong-form deviation does not disappear when the supports are placed at the neutral surface; it becomes larger. Because the response depends on where the hinges act, a study of asymmetric beams with restrained ends must state the support plane, as also observed for plates in [24].

6.4. Initial Conditions

How the thermal state is introduced changes the response substantially (Figure 6). If heating and loading start together from a flat beam at rest, the thermal bending is itself excited dynamically and the weak-form peak is 0.077. Starting from the thermal equilibrium state it is 0.101, oscillating between W₀ = −0.024 and 0.101. The amplitude of 0.125 is the same as that of the strong form, but the strong form starts from zero. A study must therefore state which loading sequence it models. Since the problem statement applies the temperature before the load, the equilibrium start is used throughout.

6.5. Porosity

At 400 K with the base foundation, pores increase the peak of type A, from 0.099 at ξ = 0 to 0.105 at ξ = 0.3 for even pores, and less for uneven pores (0.101). They slightly decrease it for type B, from 0.085 to 0.081 (Figure 7a). Separating the response shows why. For uneven pores, the thermal pre-deflection of type B hardly changes (−0.0092h to −0.0086h as ξ goes from 0 to 0.3), whereas the load-induced amplitude falls from 0.094 to 0.090. At 300 K the same pores raise the peak slightly (0.0726 to 0.0729 for uneven and to 0.0745 for even pores). Pores lower E and with it both the bending stiffness and the axial thermal force N^T; at 400 K the reduction of the thermal softening outweighs the loss of stiffness. This is consistent with earlier findings that porosity is detrimental at room temperature without a foundation [15], and shows that the sign of the porosity effect depends on the thermal state.

The strong-form deviation in Figure 7b remains within a minor margin at 300 K and expands to a notable and significant margin at 400 K, remaining almost entirely independent of ξ. This contrasts with the linear Al/Al₂O₃ results of Table 3, where porosity amplified the deviation, because here the deviation comes mostly from the thermal pre-deflection rather than from coupling stiffness.

6.6. Volume Fraction Index and Face-Sheet Thickness

A larger N lowers the ceramic content and raises the peak of both types (Figure 8a): from 0.089 to 0.106 for type A and from 0.080 to 0.085 for type B as N goes from 0 to 10. The strong form overestimates the type B peak by a notable and steady margin over this range.

Thicker face sheets lower the peak of type A, from 0.106 at h_f/h = 0.05 to 0.095 at 0.35, in line with earlier results [15] (Figure 8b). For type B the two forms disagree qualitatively. The strong form predicts almost no sensitivity to h_f/h (0.093 to 0.092), consistent with previously reported strong-form results [15]. The weak form gives a a moderate decrease, from 0.086 to 0.081. Thicker metal and ceramic faces increase the asymmetry and hence the thermal moment and the upward pre-deflection. The weak sensitivity of type B predicted by the strong-form approximation therefore reflects the omitted thermal moment rather than a property of the beam.

6.7. Initial Imperfection

The imperfection lowers the peak at room temperature and raises it at 400 K for both types (Figure 9). For type A it goes from 0.077 to 0.063 at 300 K and from 0.100 to 0.136 at 400 K as μ increases from 0 to 0.3. For type B the corresponding values are 0.073 to 0.056 and 0.083 to 0.107. The strong form shows the same trends remaining higher by a notable and consistent margin at 400 K. The reversal reported in earlier strong-form studies [15] therefore also holds in the weak-form formulation: at room temperature the imperfection stiffens the beam through the membrane coupling, while at elevated temperature the axial thermal force acting on the imperfection behaves as an additional load.

6.8. Foundation Stiffness

Table 8 lists the peak deflection for several foundations at 400 K. Stiffer foundations lower the peak of both types, by a factor of about 6.5 for type A and 5.3 for type B between no foundation and (K₁, K₂) = (150, 20). The shear layer is particularly effective: adding K₂ = 10 to K₁ = 100 reduces the peak by by a substantial and marked margin. The strong-form overestimate for type B stays within a notable and consistent margin across all foundations.

7. Conclusions

The nonlinear transient response of porous FGM sandwich beams with elastically restrained ends in thermal environments has been analysed with a refined HSDT and a weak-form, multi-term Galerkin solution that retains the moments at the beam ends. The main conclusions are:

• In asymmetric beams, the moment at a restrained hinged end depends on the axial reaction, through bending–extension coupling, and on the thermal moment. These effects are well documented in static and free-vibration studies, but a one-term Galerkin projection of the strong-form equations, a common simplification in transient analyses, does not capture them; the beam then behaves as if hinged at the neutral surface and free of thermal bending. For immovable Al/Al₂O₃ FGM beams hinged at the mid-plane this underestimates the fundamental frequency by a notable and progressive margin, with the deviation growing with porosity. The weak form converges to the exact value within an exceptionally tight margin using five term, and agrees with an independent nonlinear finite element model within a similarly close margin in the transient response.
• For symmetric sandwich beams (type A) the end moment vanishes, and the strong and weak forms give identical results. Published results for such beams are unaffected.
• For asymmetric Si₃N₄/SUS304 sandwich beams (type B), the strong form overestimates the peak dynamic deflection by a minor margin at 300 K, notable margin at 400 K and and a substantial margin between 500 and 650 K. The load-induced amplitude is almost the same in both forms; the difference is the upward thermal pre-deflection caused by the thermal moment, which the strong form omits. In higher-contrast material pairs the coupling stiffness contributes as well.
• The support plane matters. With hinges at the neutral surface the coupling term vanishes but the thermal moment about the neutral axis remains, the thermal pre-deflection triples, and the strong-form deviation grows (0.059 against 0.125 in peak deflection at 500 K). Studies of asymmetric beams with restrained ends should state where the supports act.
• The order in which heat and load are applied matters. Starting from the thermal equilibrium state, rather than applying heating and loading together, changes the peak deflection of type B by a marked and significant margin at 500 K. Studies should state which sequence they model.
• For asymmetric beams, some trends obtained with strong-form approximations change while others remain valid. The weak sensitivity of type B beams to the face-sheet thickness does not persist: the weak-form peak falls by a moderate margin as h_f/h increases from 0.05 to 0.35. Other trends remain valid: initial imperfection still reduces the dynamic deflection at room temperature and increases it at elevated temperature, and a larger volume fraction index still increases it.
• The effect of porosity depends on the thermal state. At 400 K with a supporting foundation, pores increase the peak deflection of type A beams but slightly reduce that of type B beams, because they reduce the axial thermal force more than the bending stiffness.
• Three sine terms yield converged nonlinear transient responses. Thickness stretching, checked using a quasi-3D theory, changes the peak deflection of the Si₃N₄/SUS304 beams by an imperceptible and minimal margin down to L/h = 5, indicating that a refined HSDT is fully adequate for this class of problems.

The current findings assume an isothermal temperature distribution through the thickness, the absence of damping, and neglected axial or rotary inertia. Incorporating through-thickness temperature gradients—which introduce thermal moments to symmetric beams as well—alongside damping effects and experimental validation represent natural directions for future extensions.

8. Patents