Dynamic Behavior of Two-Layered Beam Subjected to Mechanical Load in Thermal Environment

1. Introduction

In recent decades, the dynamic behavior of thin-walled structures subjected to multiphysics fields has been the topic of increasing investigations. One explanation for this is the fast growth of micro- and nanotechnologies. Many micro-electromechanical and nanoelectromechanical systems are subjected to the influence of temperature, electrical, magnetic, or fluid fields. The interaction between mechanical fields and concomitant temperature, electrical, or other fields could be significant for the design, proper exploitation, and safety of structures. Often, the part of structures like sensors, membranes, and other micro-electronic devices are composite structures consisting of different materials. Such structures could be found in many other technological areas, such as the aerospace industry and chemical and energy technologies.

The thermomechanical issue is one of the most commonly encountered in engineering practice and is extensively researched in scientific literature. The classical books that present the basic problems of thermoelasticity are Boley and Weiner [1] and Nowacki [2].

Jones’ article [3] from 1963 was one of the first to consider the equations of motion for 3D beams, taking into account the temperature propagation along all beam dimensions. It used Euler-Bernoulli and Timoshenko beam theories but lacked numerical examples.

The evolution of dynamic thermoelastic investigations is closely related to the aerospace industry. A series of works studies the phenomena arising due to the combination of airflow, thermal, and mechanical fields acting on aerospace structures. Some of them are systematized in Torton’s book [4].

Other works (referenced as [5,6,7,8]) further explore flutter, buckling, and post-buckling behaviors in panels exposed to thermal environments. In [9], the influence of thermal gradients on flexural waves in bilayer beams, indicating potential for unbounded growth, is discussed.

An important characteristic of many works devoted to dynamic thermoelastic problems is that the subjects of investigation are often composite structures. Most often a homogenization procedure is applied for such structures and the problem is transformed to the one of a structure with generalized properties. For example, in [10,11] beams with functionally graded materials are considered and dynamic and stability problems are solved in the case of thermal loading. The dependence of material gradations on the critical buckling temperature is shown.

11In the previously mentioned references and many other papers treating dynamic thermoelastic issues, it is generally accepted that temperature influences the mechanical response. In contrast, the effect of structural motion on the temperature field is often overlooked. Some researchers have attempted to address this by solving the fully coupled problem, considering the simultaneous evolution of both mechanical and thermal fields [12,13,14,15,16,17]. For instance, the elastic-plastic behavior of beams exposed to various heat impacts is examined in [12], employing finite difference discretization to model both beam motion and temperature change and then solving the resulting system of ordinary differential equations for both fields concurrently. The nonlinear response of a clamped circular plate under thermal loads, including a constant and a harmonically varying temperature component, is analysed in [13]. The study in [14] focuses on the free vibrations and buckling of clamped annular plates subjected to axisymmetric in-plane thermal loads, investigating how different thermal and geometric parameters affect the system’s response. Integral transformations such as Fourier and Laplace are utilized to analyse the coupled thermoelastic vibrations of a Timoshenko beam under a step heat flux. In [15], Henkel and Laplace integral transformations are applied by Trajkovski and Ĉukić to explore the coupled vibrations of circular plates. Considering the differing propagation velocities of mechanical and temperature oscillations, an algorithm for a sequential approach to solving the equations of temperature propagation and beam vibrations is proposed in [16,17].

In [18], the equations of motion of a simply supported rectangular plate subjected to thermal and mechanical loads are transformed to three ordinary differential equations by the Galerkin method. The authors have computed numerically different nonlinear features as Poincare maps, phase plots, power spectra, and bifurcation diagrams. In [19], thermoelastic vibrations of a rotating micro-beam subjected to a laser pulse heat source and sinusoidal heating are considered. Mathematical modelling of the problem is developed using the Euler–Bernoulli beam theory and the generalized theory of thermal conductivity with three relaxation coefficients of time. The closed-form solutions are obtained by Laplace’s transform. Laplace’s transform is again applied in [20] to study the thermoelastic vibration of functionally graded nano-beams interacting with abrupt heat in nonclassical thermoelasticity with phase delays. In [21], the influence of temperature on vibrations of laminated layers made of two different materials is presented. Essentially, the object of the investigation can be considered as a functionally graded beam. The systematic and consistent studies of the nonlinear thermo-elastic dynamic behavior of plates are presented in [22,23,24,25,26,27,28].

In [22], vibrations of a laminated von Karman plate with a linear temperature distribution along the plate thickness are considered. The obtained continuous model is subjected to a minimum reduction via the Galerkin procedure, ending with a reduced model of three coupled ODEs. The linear distribution of the temperature along the plate thickness is accepted. In [23,24], the third-order plate theory is used for thermoelastic vibration of the laminated plate, as the variation in plate reference temperature varies according to cubic law. In [24], parametric investigation of the response is accomplished using bifurcation diagrams, phase portraits, and planar cross-sections of the four-dimensional basins of attraction to describe the model’s local and global dynamical behavior. Different scenarios induced by the boundary conditions are highlighted, along with the effect of varying modelling of temperature distribution on the system response. Assessing the possible occurrence, disappearance, or modification of mechanical buckling as a result of thermal influence is discussed in [25]. The effects of coupling are explored, and the crucial role played by the slow thermal transient evolution in modifying the fast, steady mechanical response is studied. The same authors studied nonlinear dynamics of plates subjected to thermal and mechanical loading in many different aspects [26,27]. In [27], the coupled problem of the nonlinear dynamics of composite plates is analytically treated, and the main bifurcation phenomena are studied. The resulting modulation equations and the steady-state mechanical and thermal responses are determined and compared with the numerical solutions.

The nonlinear dynamics of the coupled thermoelastic vibrations of isotropic and orthotropic plates are discussed [28,29], respectively. Various characteristics of nonlinear dynamics, including phase diagrams, Poincaré maps, fractal dimensions, bifurcation diagrams, power spectra, and Lyapunov exponents, are employed to describe the behaviour of these structures.

Laminated composite structures such as plates and beams are mostly used structures concerning dynamic thermoelasticity. One of the approaches in this case is to create a homogenized model of the structure using generalized properties of the material. Then, often reduced models are developed. This is a reasonable approach for layers with equal thickness and symmetrically disposed along the thickness. When homogenization cannot be applied, usually the finite elements method with 3 D elements is used. The cost for the computations, in this case, is very high. Other often-used structures considered objects of thermoelastic vibrations are functionally graded structures.

Reference [30] presents the analysis of the large amplitude thermo-mechanically coupled vibration of the simply supported orthotropic rectangular thin plate. The governing partial differential equation of the large deflection of thermo-mechanical coupling is derived and reduced to a set of three non-linear ordinary differential equations by the Galerkin method.

References [31,32] examine the effects of the coupled field on the response of functionally graded beams. Additionally, reference [33] investigates how moisture and temperature influence wave propagation in porous FG sandwich plates, emphasizing their impact on material expansion.

The authors of the present work modelled thermoelastic vibrations of circular plates in [34,35,36]. The complete model considers transversal and longitudinal displacements, taking into account the geometrical nonlinearity, shear deformation, inertia terms due to cross-section rotation, and thermal and mechanical loadings. Based on that model, temperature influences the first resonance zone and the bifurcation scenario, leading to buckling and chaotic oscillations, as presented in [34,35]. In [36], the reduced model is based on three-mode reduction and considers the effect of temperature, which is regarded as uniformly distributed through the plate. This paper demonstrates that all three modes are involved in the plate dynamics, and the elevated temperature essentially changes the plate’s response. The developed analytical model enables the detection of the buckling bifurcation point, post-buckling oscillations, and then a series of period-doubling bifurcations and zones of multistable solutions. A special case of composite structures is the bi-material beams. Usually, the two layers have quite different properties and thicknesses. Examples of such structures can be found as parts of many micro-electronic devices (MEMS), energy harvesting devices, etc. Bi-material Euler-Bernoulli beams with L and inverted T’s cracks are considered in [37]. An analytical method for determining dimensions of longitudinal cracks based on frequency measurements is developed. When such structures are subjected to mechanical and thermal loads, internal stresses could arise due to the difference in the mechanical and thermal properties of the structure. Most of the studies of the thermomechanical behaviour of the bi-material beams are devoted to their static states. Some consideration of the bending of bimaterial beams can be found as early as 1960 in the book of Boley and Weiner [1]. Later, many studies appear about bi-metallic beams’ deformation, stresses, and temperature distribution. The static thermoelastic deformation of composite beams is studied in [38,39,40,41]. The dynamic instabilities and transient vibrations of a bi-material beam with alternating magnetic fields and thermal loads are investigated in [42]. The authors used the Hamilton principle to deduce the equation of the beam’s vibration on the base of the Euler- Bernoulli theory. The equation of motion and the solution of thermal effect are obtained by superposing specific fundamental linear elastic stress states. The authors of this study investigated coupled vibrations of bi-material beams subjected to mechanical and thermal loads in [43] on the base of a reduced model and successful solution of the beam vibration equation and the temperature propagation equation. The influence of the elevated temperature on the nonlinear vibration, buckling and post-buckling behaviour of the beam has been presented. In the following paper [44], the model of homogenization of the beam has been applied to study in detail the dynamics of the beam in the frequency domain. Bifurcation analysis is performed, and conclusions about the possibility of arising multiple solutions, loss of stability, and bifurcation are presented. The reduced model is based on three modes reduction. The beam was somewhat thick, and getting the nonlinear response required large loads. The attempt for an experimental study approved this.

This work applies a similar approach but for a much thinner beam. In addition, the case of the temperature linearly distributed along the beam thickness is studied. This makes the reduced model more complicated and allows the discovery of new phenomena that arise due to temperature influence.

2. Mathematical Model

Figure 1 shows a composite beam consisting of two layers of different materials (Material 1 and Material 2) with thicknesses h₁ and h₂. The length of the beam is denoted with l, the width with b, and the thickness is h where h= h₁+ h₂ . The coordinate along $z$ axis of the interface between the layers is denoted with $z_m$.

Assuming large displacements and applying the Timoshenko beam theory the strain and curvature-displacement relationships associated with the mid-axes can be expressed as:

$${\epsilon }_x^0={\displaystyle \frac{\mathsf{\partial }u}{\mathsf{\partial }x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\mathsf{\partial }w}{\mathsf{\partial }x}}\right)}^2{\epsilon }_{xz}^0={\displaystyle \frac{\mathsf{\partial }w}{\mathsf{\partial }x}}-\psi {\kappa }^0=-{\displaystyle \frac{\mathsf{\partial }\psi }{\mathsf{\partial }x}}$$

(1)

and the strain vector is expressed as follows:

$$\epsilon ={\left\{{\epsilon }_x^0+z{\kappa }_x^0, f(z){\epsilon }_{xz}^0\right\}}^T$$

(2)

where $w\left(x,t\right)$is the transverse displacement, $\psi \left(x,t\right)$ is the rotation angle,$f(z)$ , according to the Timoshenko beam theory is a function describing the distribution of the shear strain along the beam thickness. The upper index ⁰ associates the strains ${\epsilon }_x, {\epsilon }_{xz}$and the curvature $\kappa$ with the mid-axes.

Constitutive Equations

The relationships between the stress $S=\{{\sigma }_{\mathrm{x}},{\sigma }_{\mathrm{xz}}\}$  and strains $\epsilon =\{{\epsilon }_{\mathrm{x}},{\epsilon }_{\mathrm{xz}}\}$, considering the temperature changes, are presented by the well-known equations:

$${\sigma }_x^{(i)}=E^{(i)}\left[{\epsilon }_x-{\alpha }_T^{(i)}(T-T_i)\right], {\sigma }_{xz}^{(i)}=G^{(i)}{\epsilon }_{xz}\hspace{1em}i=1,2$$

(3)

where E⁽ⁱ⁾ is the Young modulus of i^th layer, G⁽ⁱ⁾ is the corresponding shear modulus and ${\alpha }_T^{(i)}$ is the coefficient of thermal expansion, T is current temperature, $T_i$ is the initial (environmental) constant temperature.

Equations of Motion

Starting from the equilibrium equations written for each layer and then following the procedure for homogenization described in [44] the following equations are obtained:

$${\displaystyle \frac{{\mathsf{\partial }}^{2}u}{\mathsf{\partial }{x}^2}}=G_u+G_u^T$$

(4)

$${\displaystyle \frac{{\mathsf{\partial }}^2\psi }{\mathsf{\partial }{x}^2}}-{\displaystyle \frac{\overline{Gh}}{\overline{EIZ}}}{k}^2\left({\displaystyle \frac{\mathsf{\partial }w}{\mathsf{\partial }x}}-\psi \right)-c_2{\displaystyle \frac{\mathsf{\partial }\psi }{\mathsf{\partial }t}}-{\displaystyle \frac{\overline{\rho I}}{\overline{EIZ}}}{\displaystyle \frac{{\mathsf{\partial }}^2\psi }{\mathsf{\partial }{t}^2}}=G_1^T$$

$$k^2{\displaystyle \frac{\overline{Gh}}{\overline{Eh}}}\left({\displaystyle \frac{{\mathsf{\partial }}^{2}w}{\mathsf{\partial }{x}^2}}-{\displaystyle \frac{\mathsf{\partial }\psi }{\mathsf{\partial }x}}\right)-c_1{\displaystyle \frac{\mathsf{\partial }w}{\mathsf{\partial }x}}-{\displaystyle \frac{\overline{\rho h}}{\overline{Eh}}}{\displaystyle \frac{{\mathsf{\partial }}^{2}w}{\mathsf{\partial }{t}^2}}=-p_1(x,t)+G_2^L+G_2^T$$

where c₁ and c₂ denote damping coefficients which are assumed to be proportional to the mass terms in front of the inertia terms. However, the inertia term in the longitudinal direction is neglected.

In Eq. (4) the following notations are introduced:

$$\overline{E}\overline{I}={\displaystyle \frac{1}{3}}\left\{E^{(1)}\left(z_1^3+{\displaystyle \frac{h^3}{8}}\right)+E^{(2)}\left({\displaystyle \frac{h^3}{8}}-z_1^3\right)\right\}$$

(5)

$$\overline{E}\overline{I}\overline{Z}=\left[\overline{E}\overline{I}-{\displaystyle \frac{{\left(E^{(1)}-E^{(2)}\right)}^2}{\overline{Eh}}}{\displaystyle \frac{{(h^{(1)}{h}^{(2)})}^2}{4}}\right]$$

$$\begin{array}{l}{G}_u=-\left({\displaystyle \frac{\mathsf{\partial }w}{\mathsf{\partial }x}}{\displaystyle \frac{{\mathsf{\partial }}^{2}w}{\mathsf{\partial }{x}^2}}\right)-{\displaystyle \frac{\left(E^{(1)}-E^{(2)}\right)}{\overline{Eh}}}{\displaystyle \frac{h^{(1)}{h}^{(2)}}{2}}{\displaystyle \frac{{\mathsf{\partial }}^2\psi }{\mathsf{\partial }{x}^2}}\\ G_u^T={\displaystyle \frac{1}{\overline{Eh}}}(E^{(1)}{\alpha }_T^{(1)}{\displaystyle \frac{\mathsf{\partial }{\gamma }_T^{(1)}}{\mathsf{\partial }x}}+E^{(2)}{\alpha }_T^{(2)}{\displaystyle \frac{\mathsf{\partial }{\gamma }_T^{(2)}}{\mathsf{\partial }x}})\end{array}{p}_1(x,t)={\displaystyle \frac{p(x,t)}{\overline{Eh}}}$$

(6)

$$G_1^T==-\left(E^{(1)}{\alpha }_T^{(1)}{\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }x}}{\chi }_T^{(1)}+E^{(2)}{\alpha }_T^{(2)}{\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }x}}{\chi }_T^{(2)}\right){\displaystyle \frac{1}{\overline{EIZ}}}$$

(7)

$$\begin{array}{l}{G}_2^L=-\left[\left({\displaystyle \frac{\mathsf{\partial }u}{\mathsf{\partial }x}}+0.5{\left({\displaystyle \frac{\mathsf{\partial }w}{\mathsf{\partial }x}}\right)}^2\right)+{\displaystyle \frac{\left(E^{(1)}-E^{(2)}\right)h^{(1)}{h}^{(2)}}{2\overline{Eh}}}{\displaystyle \frac{\mathsf{\partial }\psi }{\mathsf{\partial }x}}\right]{\displaystyle \frac{{\mathsf{\partial }}^{2}w}{\mathsf{\partial }{x}^2}}\\ G_2^T=\left[E^{(1)}{\alpha }_T^{(1)}{\gamma }_T^{(1)}+E^{(2)}{\alpha }_T^{(2)}{\gamma }_T^{(2)}\right]{\displaystyle \frac{{\mathsf{\partial }}^{2}w}{\mathsf{\partial }{x}^2}}{\displaystyle \frac{1}{\overline{Eh}}}\end{array}$$

(8)

$$p_1(x,t)={\displaystyle \frac{p(x,t)}{\overline{Eh}}}$$

(9)

$$\begin{array}{l}{\gamma }_T^{(1)}={\displaystyle \underset{-h/2}{\overset{z_1}{\int }}Tdz},\hspace{1em}{\gamma }_T^{(2)}={\displaystyle \underset{z_1}{\overset{h/2}{\int }}Tdz};{\chi }_T^{(1)}={\displaystyle \underset{-h/2}{\overset{z_1}{\int }}Tzdz},\hspace{1em}{\chi }_T^{(2)}={\displaystyle \underset{z_1}{\overset{h/2}{\int }}Tzdz},\\ {\gamma }_T^{(1)}={\displaystyle \underset{-h/2}{\overset{z_1}{\int }}Tdz},\hspace{1em}{\gamma }_T^{(2)}={\displaystyle \underset{z_1}{\overset{h/2}{\int }}Tdz};{\chi }_T^{(1)}={\displaystyle \underset{-h/2}{\overset{z_1}{\int }}Tzdz},\hspace{1em}{\chi }_T^{(2)}={\displaystyle \underset{z_1}{\overset{h/2}{\int }}Tzdz},\end{array}$$

(10)

The equation describing the temperature propagation in the beam , without considering full coupling with the mechanical field is:

$${\displaystyle \frac{c_p^{(i)}}{{\lambda }_T^{(i)}}}{\displaystyle \frac{\mathsf{\partial }T}{\mathsf{\partial }t}}={\displaystyle \frac{{\mathsf{\partial }}^{2}T}{\mathsf{\partial }{x}^2}}+{\displaystyle \frac{{\mathsf{\partial }}^{2}T}{\mathsf{\partial }{z}^2}}i=1,2$$

(11)

where $c_p^{(i)} \mathrm{and} {\lambda }_T^{(i)}$ are the heat capacity per unit volume and the thermal conductivity of i^th material .

In this work, the temperature field is established and the temperature is uniformly distributed along the beam length. This means that the temperature does not change in time and the derivatives along x axis are zero. Therefore,

$$G_u^T=0\hspace{1em}\mathrm{and} G_1^T=0.$$

Concerning the distribution of the temperature along beam thickness two cases are considered: when $\Delta T=T-T_i$—the temperature is elevated/lowered at the whole beam and the second case when the temperature is distributed along the beam thickness following the linear low:

$$T(z)=T_0+T_{1}z$$

(12)

The constants T₀ and T₁ are defined by the temperature boundary conditions.

The following boundary conditions about the temperature are defined:

$$\begin{array}{l}T(-h/2)=T_d\\ T(h/2)=T_{up}\end{array}$$

(13)

where T_d, T_up denote temperature on bottom (down) and upper beam surfaces.

In this case the coefficients in equation (12) are:

$$T_0={\displaystyle \frac{T_{up}+T_d}{2}},\hspace{1em}{T}_1={\displaystyle \frac{T_{up}-T_d}{h}}$$

(14)

Considering that (see Figure 1):

$$\begin{array}{l}{h}_2=h/2-z_m\\ h_1=z_m+h/2\end{array}$$

(15)

it is easy to obtain that in the case of linear distribution of the temperature along the beam thickness:

$$\begin{array}{l}{\gamma }_T^{(1)}={\displaystyle \underset{-h/2}{\overset{z_1}{\int }}(T_0+T_{1}z)dz}=T_0{h}_1-0.5T_1{h}_1{h}_2\\ {\gamma }_T^{(2)}={\displaystyle \underset{z_1}{\overset{h/2}{\int }}(T_0+T_{1}z)dz}=T_0{h}_2+0.5T_1{h}_1{h}_2\end{array}$$

(16)

while in the case of uniform temperature:

$${\gamma }_T^{(i)}=h_i\Delta Ti=1,2$$

(16)

$$\begin{array}{l}{\gamma }_T^{(1)}=\Delta T{h}_1=(T-T_i)h_1\\ {\gamma }_T^{(2)}=\Delta T{h}_2=(T-T_i)h_2\end{array}$$

In the further text the notation ΔT means a difference between the beam and ambient temperatures.