In many applications, coupling between thermal and mechanical domains can significantly influence structural dynamics. Analytical approaches to study such problems have previously used assumptions such as a proscribed temperature distribution or one-way coupling to enable assessments. In contrast, time-stepping numerical simulations have captured more detailed aspects of multiphysics interactions at the expense of high computational demands and lack of insight of the underlying physics. To provide a new tool that closes the knowledge gap and broadens potential for analytical techniques, this research formulates and analytically solves a thermomechanical beam model considering a combination of thermal and mechanical excitations that result in extreme nonlinear behaviors. Validated by experimental evidence, the analytical framework facilitates the prediction of the nonlinear dynamics of multi-degree-of-freedom structures exhibiting two-way thermomechanical coupling. The analysis enables the investigation of mechanical and thermomechanical impedance metrics as a means to forecast future nonlinear dynamic behaviors such as extreme bifurcations. For the first time, characteristics of mechanical impedance previously reported to predict the onset of dynamic bifurcations in discrete systems are translated to illuminate the nearness of distributed parameter structures to bifurcations. In addition, fundamental connections are discovered in the thermomechanical evaluations between nonlinear low amplitude dynamics of the postbuckled beam and the energetic snap-through vibration that are otherwise hidden by studying displacement amplitudes alone.

Introduction

Interactions between thermal and mechanical domains have considerable relevance in applications ranging from microcantilever sensors [1] to shape memory material systems [2] to skin panels of hypersonic aircraft [3], among others. Temperature variations may greatly affect the static configuration and subsequent displacement amplitudes of structures, leading to phenomena such as changing resonant frequencies [4], buckling [5], and temperature-dependent material properties [6]. Likewise, mechanical response may subsequently influence the temperature distribution in the structure, such as through changes to the convective heat transfer coefficient [7,8]. Because of this, studying one- and two-way coupling between mechanical and thermal domains is important to assess how the interactions alter the structural dynamics.

To this end, extensive numerical investigations have been conducted by Culler and McNamara [9] [10], Miller and McNamara [11], and Miller et al. [12] to assess the interaction between coupled structural, thermal, and fluid physics. Finite element methods have been formulated by Daneshjo and Ramezani [13] and Carrera et al. [14] to study the linear dynamics of laminate plates exhibiting rich coupling between thermal and mechanical domains. Reduced-order models have been shown by Matney et al. [15], Perez et al. [16], and Settimi et al. [17] to characterize the intricate thermal and mechanical coupling while requiring less computational expense than numerical integration of a finite element model. Yet, the ability to obtain fundamental insight into thermal–structural interactions via parametric studies may be limited by the case study-dependent nature and computational costs of numerical methods. To surmount such shortcomings, analytical techniques may be employed to study fundamental aspects of thermomechanical interactions to obtain insights.

For example, the nonlinear structural vibrations of plates with prescribed surface temperature distributions have been studied by Pal [18] and Lee [19]. A single mode approximation for two-way interaction between structural and thermal responses of plates undergoing large deflections was studied by Chang and Wan [20], Chang and Jen [21], and Shu et al. [22]. The authors of this work recently proposed an analytical framework by which the near- and far-from-equilibrium nonlinear dynamics of systems with multiple degrees-of-freedom may be evaluated using an equivalent linearization scheme [23]. Additionally, mechanical impedance metrics were revealed to uncover energy transfer mechanisms that help to predict the onset of dynamic bifurcations [23]. Yet, the analytical undertaking and impedance studies lack ability to provide insight into thermomechanical interactions. This research rectifies the limitations by building up a new analytical framework that accounts for thermomechanical coupling between the structural dynamics of a flat beam and the thermal environment. In the process, impedance metrics are newly revealed to elucidate the coupled dynamic behaviors and potentially forecast future dynamic response. Furthermore, Kovacic et al. [24], Amabili [25], and Yamaki et al. [26] [27] have shown that initial imperfections in geometry lead to a connection between intrawell and interwell dynamic regimes at nonzero frequencies, in contrast with the nonlinear response of symmetric structures [28]. The influence of such static asymmetries on the nonlinear structural dynamics is observed in this research experimentally, and a comparison of analytical and experimental impedances examines how imperfections may also affect impedance metrics.

The remainder of this paper is organized as follows: In Sec. 2, the thermomechanical beam model and solution procedure are developed and the experimental system is summarized. Then, experimental and analytical results are presented to verify the accuracy and utility of the analytical model. New characteristics of impedance are examined to enlighten the connections between impedance changes and the onset of bifurcations. The effects of slowly changing mean thermal load on the beam are discussed in detail, providing insight into underlying ties among nonlinear dynamic regimes. Section 6 summarizes the advancements made in the analytical efforts and the insights uncovered.

Thermomechanical Model Formulation and Solution

Governing Partial Differential Equations of Motion.

In this research, a thermomechanical model of a flat, thin beam is employed to elucidate fundamental characteristics of the interactions between the nonlinear structural dynamics and the thermal domain. The equation of motion for the transverse deflection of a beam w(x,t) of length L is given by Eq. (1). The equation accounts for large beam curvature [29] and an applied axial force PM. Note that the axial force PM can be either positive (tensile force) or negative (compressive force). Equation (1) is coupled into the thermal domain by a thermally induced bending moment MT(x,t) and a thermally induced axial load PT(x,t) 
EI2x2[2wx2(x,t)(1+12(wx(x,t))2)]x[(PMPT(x,t))wx(x,t)]+ρA2wt2(x,t)=f(x,t)2MTx2(x,t)
(1)

The parameters in Eq. (1) are Young's modulus E, moment of inertia I, volumetric density ρ, and cross-sectional area A. A transverse force per unit length f(x,t) is applied to the beam. As discussed in Sec. 3, a spring steel beam is considered here. Non-negligible changes to the material properties of steel occur for temperature changes on the order of hundreds of degrees Celsius [30]. Since the range of temperature changes considered in this research is around 3 °C, temperature-dependent material properties are neglected in the present model formulation.

Equation (1) is coupled with the equation of motion (2) for the thermal moment MT(x,t) along the beam length 
(2x21κt)MT(x,t)+ηt[2wx2(x,t)]+μ[q0(x,t)p0(x,t)]=0
(2)
In Eq. (2), κ, η, and μ are parameters that depend upon the material properties and beam geometry [31], and q0(x,t) and p0(x,t) are the heat fluxes over the top and bottom surfaces of the beam, respectively. The parameters κ, η, and μ are defined in Eq. (3) 
κ=Eαbh2λ012cε
(3a)
 
η=Eα2bh2Tref12λ0[3(B23G)+2G]
(3b)
 
μ=Eαbh212λ0
(3c)
In Eq. (3), λ0 is the thermal conductivity, cε is the specific heat at constant strain, Tref is the reference temperature, B is the bulk modulus, G is the shear modulus of the beam, α is the coefficient of linear thermal expansion, b is the width of the beam, and h is the thickness of the beam [32]. The thermal moment MT(x,t) in Eqs. (1) and (2) and thermally induced axial load PT(x,t) are given by Eqs. (4) and (5), respectively [33] 
MT(x,t)=Eαbh/2h/2zT(x,z,t)dz
(4)
 
PT(x,t)=Eαbh/2h/2T(x,z,t)dz
(5)
In Eqs. (4) and (5), T(x,z,t) is the temperature distribution in the beam. Assuming a linear temperature distribution through the thickness of the thin beam in the form of Eq. (6), which is accurate for one-dimensional steady-state heat conduction and temperature-independent thermal properties [34], Eqs. (4) and (5) are expressed using Eqs. (7) and (8), respectively 
T(x,z,t)=12(Tf(x,t)+T0(x,t))+zh(Tf(x,t)T0(x,t))+Tm
(6)
 
MT(x,t)=Eαbh212θ(x,t)
(7)
 
PT(x,t)=Eαbh2θ(x,t)+EαbhTm
(8)

In Eqs. (7) and (8), θ(x,t)=Tf(x,t)T0(x,t) is the temperature difference between the top and bottom surfaces of the beam through the thickness. In Eqs. (6) and (8), Tm is the mean temperature of the beam relative to a reference temperature Tref, about which the temperature difference θ(x,t) varies.

The Ritz method is used to derive governing equations that are solved to obtain approximate solutions to the system described by Eqs. (1) and (2). The Ritz method expansions used here for the beam displacement and temperature difference are given by Eqs. (9) and (10), respectively 
w(x,t)=j=1Nψj(t)βj(x)
(9)
 
θ(x,t)=j=1Mτj(t)ϕj(x)
(10)

The trial functions βj(x) and ϕj(x) are chosen to satisfy the essential mechanical and thermal boundary conditions, respectively. As discussed in Sec. 3, the experimental system used to verify this model has clamped and isothermal boundary conditions. The hierarchical function set proposed by Beslin and Nicolas [35] is used to meet the clamped boundary conditions for the beam. The isothermal boundary conditions correspond to a known temperature at both ends equal to Tm. Sinusoids are often used in analysis to meet isothermal boundary conditions [19,20], so the temperature trial functions are presumed to be of the form sin2(jπx/L) where j is the index.

Derivation of Ordinary Differential Equations for Generalized Displacements.

This Sec. 2.2 derives the ordinary differential equations (ODEs) that govern the response of the generalized displacements ψj(t). Substituting the thermal moment (7) and thermal axial load (8) into Eq. (1) yields the following equation: 
EI2x2[2wx2(x,t)(1+12(wx(x,t))2)](PMEαbhTm)2wx2(x,t)+Eαbh2x[θ(x,t)wx(x,t)]+ρA2wt2(x,t)=f(x,t)Q2θx2(x,t)
(11)
In Eq. (11), Q=Eαbh2/12 is a constant of proportionality between the thermal moment MT and the temperature difference θ. Substituting the expansions (9) and (10) into Eq. (11), multiplying by the ith trial function βi(x), and integrating over the length L of the beam results in Eq. (12). 
Mijψ̈j(t)+γKijψ̇j(t)+Kijψj(t)+Kij(1)ψj(t)+Kijk(th)τj(t)ψk(t)+Kijkl(3)ψj(t)ψk(t)ψl(t)=Fi(t)χijτj(t)
(12)
In Eq. (12), the repeated indices entail summation over the respective index, using standard index notation. The coefficients in Eq. (12) are defined by the integrals (13)(19) 
Mij=ρA0Lβi(x)βj(x)dx
(13)
 
Kij=EI0Lβi(x)βj(x)dx
(14)
 
Kij(1)=(PMEαbhTm)0Lβi(x)βj(x)dx
(15)
 
Kijk(th)=Eαbh2[0Lβi(x)ϕj(x)βk(x)dx+0Lβi(x)ϕj(x)βk(x)dx]
(16)
 
Kijkl(3)=EI2[0L2βj(x)(βk(x)βl(x)+βk(x)βl(x))βi(x)dx+0Lβj(x)(βl(x)βk(x)+2βk(x)βl(x))βi(x)dx+0Lβk(x)(βl(x)βk(x)+βj(x)βl(x))βi(x)dx]
(17)
 
Fi(t)=0Lf(x,t)βi(x)dx
(18)
 
χij=Q0Lβi(x)ϕj(x)dx
(19)

Damping proportional to the linear stiffness matrix, using the relationship [C]=γ[K], is introduced in Eq. (12). Note that throughout the results presented in Secs. 4 and 5, the proportionality constant γ is taken to be 15 ms. The selection of γ may greatly affect the analytically predicted nonlinear dynamics. The proportionality constant γ = 15 ms in the analysis yields predictions in agreement with the experimental data, thus supporting its selection.

Equation (12) is the ith equation of the N ODEs governing the motion of the generalized displacements ψn(t). Nonlinearities are manifest via terms with matrices Kijk(th) and Kijkl(3) that signify the roles of nonlinear thermomechanical coupling and nonlinear beam bending.

Derivation of Ordinary Differential Equations for Generalized Temperatures.

This Sec. 2.3 derives the ODEs that govern the response of the thermal generalized temperatures τn(t) and couple to Eq. (12). Substituting the thermal moment (7) into Eq. (2) yields the following equation: 
Q(2x21κt)θ(x,t)+ηt[2wx2(x,t)]+μ[q0(x,t)p0(x,t)]=0
(20)
Substituting the expansions (9) and (10) into Eq. (20), multiplying by the ith trial function ϕi(x), and integrating over the length L of the beam yields the following equation: 
Bijτ̇j(t)+K̃ijτj(t)+Λijψ̇j(t)=Pi(t)
(21)
The coefficients in Eq. (21) are defined by the integrals (22)(25). 
Bij=Qκ0Lϕi(x)ϕj(x)dx
(22)
 
K̃ij=Q0Lϕi(x)ϕj(x)dx
(23)
 
Λij=η0Lϕi(x)βj(x)dx
(24)
 
Pi(t)=μ0L(q0(x,t)p0(x,t))ϕi(x)dx
(25)

Equation (21) is the ith equation of the M ODEs governing the motion of the generalized temperature τm(t). Equation (21) is solved with Eq. (12) to obtain an approximate solution to Eqs. (1) and (2), and thus to yield a thermomechanical characterization of the beam deflection and temperature distribution when subjected to external forces and heat fluxes.

Solution Procedure and Impedance Derivation.

Equations (12) and (21) are used to form a system of N+M equations that are expressed in Eq. (26) and condensed into Eq. (27) 
[[M][0][0][0]][ψ¯̈τ¯̈]+[γ[K][0][Λ][B]][ψ¯̇τ¯̇]+[[K][χ][0][K̃]][ψ¯τ¯]+[[K(1)][0][0][0]][ψ¯τ¯]+[K¯jk(th)τjψk0¯]+[K¯jkl(3)ψjψkψl0¯]=[F¯P¯]
(26)
 
[M̂]ijÿj+[Ĉ]ijẏj+[K̂]ijyj+[K̂]ij(1)yj+[K̂]ijk(th)yjyk+[K̂]ijkl(3)yjykyl=F̂i
(27)

In Eq. (27), y¯=[ψ¯τ¯]T is a vector of the mechanical and thermal generalized coordinates. The matrices or vectors in Eq. (27) with over-hat markers are identified from the respective terms in Eq. (26).

Principles of harmonic or stochastic linearization [36] are utilized to transform the nonlinear equations (27) into equivalent linear equations (28). The equivalent stiffness matrices of Eq. (28) are defined by Eq. (29). 
[M̂]y¯̈+[Ĉ]y¯̇+{[K̂]+[K̂e](1)+[K̂e](th)+[K̂e](3)}y¯=F¯̂
(28)
 
[K̂e]ih(1)=([K̂]ij(1)yj)yh
(29a)
 
[K̂e]ih(th)=([K̂]ijk(th)yjyk)yh
(29b)
 
[K̂e]ih(3)=([K̂]ijkl(3)yjykyl)yh
(29c)
In this research, the nonlinear dynamics of a beam subjected to harmonic excitation of the form F¯̂=f̂¯s+f¯̂ejωt is studied. The heat flux applied to the beam is constant and is therefore contained in the static component of force f̂¯s. To calculate the entries of Eq. (29), solutions of the form y¯=q¯+r¯ejωt are assumed, where q¯ are the offsets and r¯ are the complex-valued amplitudes of the generalized coordinates. This assumption provides single-harmonic approximate solutions to the steady-state forced vibrations problem and is accurate for most experimentally observed dynamics [28]. Substituting the assumed solution into Eq. (27) and integrating over one period of the harmonic excitation force 2π/ω results in Eq. (30). Substituting the assumed solution into Eq. (28) and isolating the coefficient of parameters dependent upon ejωt yields Eq. (31) 
{[K̂]+[K̂](1)}q¯+[K̂]ijk(th)qjqk+12[K̂]ijk(th)rjrk+[K̂]ijkl(3)qjqkql+12([K̂]ijkl(3)+[K̂]ikjl(3)+[K̂]iklj(3))qjrkrl=0¯
(30)
 
{ω2[M̂]+jω[Ĉ]+[K̂]}r¯+([K̂]ik(1)δjk+{[K̂]ijk(th)+[K̂]ikj(th)}qk+{[K̂]ijkl(3)+[K̂]ijlk(3)+[K̂]iljk(3)}qkql+14{[K̂]ijkl(3)+[K̂]ijlk(3)+[K̂]iljk(3)}rkrl)rj=F¯
(31)

In Eq. (31), δjk is the Kronecker delta function. The coupled algebraic equations (30) and (31) are solved for the offsets q¯ and complex amplitudes r¯. Then, the physical displacement w(x,t) and thermal difference θ(x,t) are reconstructed by Eqs. (9) and (10), respectively. Thus, the solution procedure articulated in this Sec. 2.4 newly enables the analytical prediction of the nonlinear structural and thermal dynamics of a flat beam wherein two-way thermomechanical coupling is prevalent. The analysis provides a bridge between the previous nonlinear analyses [20,21] and the fidelity of numerical simulations [9,10] and reduced-order models formulations [17] that exhibit two-way thermal–structural interactions.

Finally, the mechanical and thermomechanical impedances, denoted by their respective MM and TM subscripts, are determined through Eq. (32) 
ZMM=i=1NFiβijωri
(32a)
 
ZTM=(Tm+i=1Mqi+Nϕi)/(hi=1Njωriβi)
(32b)
Mechanical and thermomechanical impedances are reconstructed at a point along the beam through the trial functions βi and ϕi. Mechanical impedance ZMM, defined in Eq. (32a), is the complex ratio of the input force to the transverse velocity of the beam at a point location. Thermomechanical impedance ZTM, defined in Eq. (32b), is the complex ratio of the static thermal gradient at a point on the beam to the transverse velocity of the beam at another location.

Experimental Setup

The experimental system used to verify the analytical modeling effort of Secs. 2.1 to 2.4 consists of an initially flat spring steel beam constrained at either end by stainless steel clamps. A photograph of the experimental setup is shown in Fig. 1(a), while Fig. 1(b) provides a schematic. The beam is 0.386 mm thick and 12.7 mm wide. The thickness of the beam is defined as the smallest dimension and is parallel to the direction of the shaker table acceleration, indicated by the white arrow in Fig. 1(a). The length of the beam between the clamped ends is 291 mm. The beam is subjected to a uniform base acceleration with slowly varying frequency via the shaker table (APS 400) and amplifier (Crown XLS2500). The shaker acceleration is measured by an accelerometer (PCB 333B40) and controlled by a vibration controller (Vibration Research VR9500). A laser (Micro-Epsilon optoNCDT 2300) measures the transverse displacement of the beam at a location 41% along the beam length. The laser data sampling rate is 4096 Hz. A heat lamp applies a thermal load, focused on the center of the beam, which causes the beam to buckle. Since the beam is heated primarily in the center, the static temperature loading can be reasonably approximated as ΔTsin2(πx/L). The static temperature distribution leads to the mean temperature of the beam to be estimated as Tm=ΔT/2, where ΔT is the maximum transverse temperature rise that occurs at the beam center. Four type K thermocouples are positioned in a horizontal plane with the beam as shown in Fig. 1(b). The sampling rate of the thermocouples is 2 Hz. Thermocouples labeled 1, 2, and 4 in Fig. 1 are offset 30 mm from the unbuckled position of the beam. The transverse temperature difference θ is estimated by thermocouples 1 and 2. A handheld infrared temperature probe is used to validate the temperature gradient along the beam length and to confirm the accuracy of the approximations made by thermocouple measurement. The positioning of the thermocouples from the heat lamp source provides an approximation of the heat flux applied to the beam.

Excitation frequency sweeps of base acceleration are conducted at a rate of 0.05 Hz/s. Such experiments are conducted with increasing and decreasing rates from 10 to 30 Hz and with constant base acceleration amplitude. In Sec. 4, the experimental results at two amplitudes of base acceleration are compared with corresponding analytical results, and new insight into aspects of impedance is obtained.

Experimental Validation of Analytical Formulation

Two cases are considered to validate the analytical formulation of Sec. 2. The heat lamp applies a heat flux of approximately 0.98 W/m2 to the surface of the beam, resulting in a static ΔT increase of 2.2 K with respect to the reference temperature 298 K. The thermal load leads the beam to buckle and induces a fundamental linear natural frequency near 22 Hz. Three mechanical generalized coordinates (N=3) are used in the analysis, since this selection provides good quantitative agreement between the experimental and analytical fundamental natural frequency, which is the focus in this research. Two thermal generalized coordinates (M=2) capture the thermomechanical coupling effects manifest in the beam dynamics and demonstrate the ability of the analysis to account for a multimode expansion in the thermal domain. Table 1 summarizes the parameters utilized in the analytical model. In the first validation case, a base acceleration amplitude of 0.1 m/s2 is applied to the beam.

Figure 2 compares the analytical and experimental results for the first validation case. The time-domain displacement measurement of the beam is converted into the frequency domain via fast Fourier transform. The experimental displacement amplitudes presented in Fig. 2(a) are the frequency response of the relative beam displacement with respect to the motion of the shaker table at the excitation frequency. The 0.1 m/s2 amplitude of base acceleration induces softening in the first resonance frequency of the beam, seen in the experimental and analytical displacement amplitudes of Figs. 2(a) and 2(d), respectively. The displacement amplitudes presented in Figs. 2(a) and 2(d) and in all subsequent figures display only the amplitude of the steady-state oscillations, and do not refer to the static equilibria of the beam. Only intrawell dynamics, oscillations about a single stable equilibrium of a postbuckled structure, are manifest for this amplitude of base acceleration. The softening leads to coexisting low- and high-amplitude intrawell dynamics around 20 Hz that are separated by bifurcations. In Fig. 2, stars denote the frequency at which a vertical tangency appears in the low amplitude intrawell regime and corresponds to a saddle-node bifurcation. Although dynamic stability is not assessed in the results of Fig. 2, it is known that analytical solutions in the low amplitude intrawell regime of Fig. 2(d) with amplitudes larger than that of the vertical tangency bifurcation are unstable [28]. This observation is supported by the lack of such results in the data, Fig. 2(a). The experimental and analytical results exhibit qualitative and quantitative agreement, although the highest intrawell amplitude exhibited in the experiment is twice as large as that predicted analytically. The discrepancy is likely due to imperfections in the experimental beam, which is explained later in Sec. 4.

Figures 2(b) and 2(e) present the experimental and analytical mechanical impedances ZMM, while Figs. 2(c) and 2(f) show the thermomechanical impedances ZTM. An important trend exhibited in the experimental impedances is predicted by the analysis. The bifurcations that occur during the transition between low and high amplitude intrawell dynamics correspond to jumps across zero reactance, the imaginary component of impedance. This behavior is clearly borne out in both analysis and experiment in Figs. 2(b), 2(c), 2(e), and 2(f). Specifically, the high-amplitude intrawell dynamic approaches zero reactance prior to bifurcation, shown by the light gray arrows in Fig. 2. Meanwhile, the bifurcation in the low amplitude dynamic denoted by stars in Fig. 2 occurs at a negative reactance. The difference of this bifurcation from zero reactance is marked in the impedance plots of Fig. 2 by black, double-sided arrows. The experimental and analytical evidence underscores that the bifurcation from low to high amplitude intrawell oscillations occurs as a jump across the zero reactance line from absolute values of reactance greater than the opposite bifurcation from high to low oscillation amplitudes. Thus, an overall comparison of the results in Fig. 2 between experiment (a)–(c) and analysis (d)-(f) indicates that the analysis accurately characterizes the nonlinear behavior of the softening intrawell dynamic. The findings also confirm that dynamic bifurcations of distributed parameter systems between intrawell regimes correspond to vanishing reactance of both mechanical and thermomechanical impedance, previously reported for mechanical impedance of discrete systems [23,37].

Results from the second experimental case utilized in this research for analytical validation are presented in Fig. 3. The same thermal conditions from the results of Fig. 2 are used, while the amplitude of base acceleration by the shaker table is 1.0 m/s2. The increased excitation level leads to greater softening of the 22 Hz resonance. Additionally, snap-through dynamics, or high amplitude oscillations between multiple static equilibria of a postbuckled structure, are manifest and are denoted in the displacement amplitude plots of Figs. 3(a) and 3(d) by solid black curves. The greatest qualitative difference between the experimental and analytical displacement amplitudes in Figs. 3(a) and 3(d) is that the higher amplitude intrawell dynamic in the experiments bends toward the snap-through regime. The snap-through dynamic appears around 16 Hz. These nuances are in contrast to the analytically predicted snap-through dynamic, denoted in Fig. 3(d) by solid black curves, which appears theoretically at zero frequency [28], although the frequency span of Fig. 3(d) is bound between 10 and 30 Hz to correspond to the experimental results.

The analytical mechanical and thermomechanical impedance results of Figs. 3(e) and 3(f) reveal that for decreasing frequency, intrawell and snap-through regimes approach each other. The nearness is observed by noting the proximity of the square markers in Figs. 3(e) and 3(f), which denote dynamics at 1 Hz. The convergence of intrawell and snap-through dynamics in the impedance plane suggests an underlying connection between these otherwise distinct dynamic regimes at zero frequency, or the condition of rigid body motion. The experimental results, meanwhile, are consistent with results reported for asymmetric Duffing oscillators [24] and plates with initial imperfections [25]. In other words, an imperfection or initial curvature in the beam produces asymmetry and results in a veering of intrawell responses toward the snap-through displacement amplitudes in Fig. 3(a). Initial beam curvature is not considered in the analysis at this stage of development, which explains the less confined analytical impedance predictions in Figs. 3(e) and 3(f) compared to the measurements in Figs. 3(b) and 3(c). Yet, the comparison between Figs. 3(b) and 3(c) and Figs. 3(e) and 3(f) indicates that impedances are influenced by initial imperfections in a manner analogous to the displacement amplitudes of asymmetric Duffing oscillators [24] by virtue of the non-zero frequency connection between intrawell and snap-through dynamics. Despite the deviation associated with asymmetry, comparisons between the experimental and analytical displacement amplitudes and impedance results of Figs. 2 and 3 confirm that the analytical model suitably characterizes the rich, nonlinear dynamic responses of the harmonically forced beam in a thermomechanical environment that is analyzed here for the first time. Future efforts will account for the nuances associated with asymmetry in the mechanical configuration of the beam.

The qualitative discrepancy observed in Fig. 3 between the experimental (a)–(c) and analytical (d)–(f) results is attributed to initial imperfections in the experimental beam. Other sources of error, such as the damping model and proportionality constant discussed in Sec. 2.2, may result in deviations. Yet, it is difficult to decompose the relative errors associated with damping and asymmetry. The latter is anticipated to be the greatest contributor to the qualitative difference since the experimental results of Fig. 3(a) more closely resemble those of an imperfect duffing oscillator [24].

The analytical impedance results of Figs. 2 and 3 display deviations from the underlying linear impedances that correlate with the deviation of the displacement amplitude from the linear dynamics. In particular, the softening of the resonant peak corresponds to a “stretching out” of the analytical results along the line of zero reactance, as particularly observed in Figs. 2(e) and 2(f). The deviation from the underlying linear impedances becomes more significant as the deviation from the underlying linear displacement amplitude increases due to increasing base acceleration. Eventually, the low- and high-amplitude intrawell regimes separate and diverge in the impedance plane, denoted in Fig. 3 by the gray double-sided arrows. The onset of separation of intrawell regimes coincides with the emergence of snap-through dynamics that connect to the two intrawell dynamics in the complex impedance plane of Figs. 3(e) and 3(f). The observed changes from the linear impedance results indicate that when interwell dynamics are present, the saddle-node bifurcation that terminates snap-through corresponds to vanishing reactance. In addition, the bifurcations in the intrawell regime [38] may be forecasted by observing the deviations from the underlying linear impedance as zero reactance is approached in the complex impedance plane.

The mechanical and thermomechanical impedances of Figs. 2 and 3 exhibit similar trends, and therefore may both be used to forecast dynamic bifurcations. Yet, both characterize distinct aspects of the dynamics. Mechanical impedance relates the input force to the velocity amplitude and implicates energy management mechanisms [23], while thermomechanical impedance is a multiphysics transfer function. The similarity stems from their shared dependence on beam velocity and the constant phase of the base acceleration. The unique influence of thermal effects on the steady-state dynamics and energy of the beam is explored in Sec. 5.

Influence of Thermal Loads

In this Sec. 5, the maximum transverse temperature elevation ΔT is varied for a constant base acceleration amplitude and frequency. The effects of this changing temperature on the steady-state displacement amplitude, mechanical impedance, and thermomechanical impedance are examined.

Influence of Temperature Variation on Displacement Amplitude.

Figure 4 shows the analytical displacement amplitude predicted for a change in thermal load. The base acceleration is held constant at 1 m/s2 and 16 Hz. A dashed vertical line is shown for ΔT of 2.2 K, which corresponds to the dashed vertical line depicted in Fig. 3(d) at 16 Hz. Figure 4 indicates the presence of two linear resonances at 0.87 K and 2.03 K, which correspond to pre- and postbuckled conditions, respectively. A cusp is observed in the linear dynamic and lower amplitude nonlinear behaviors at 1.76 K, corresponding to the buckling temperature for the beam and denoted by a triangle in Fig. 4. Yet, this cusp does not manifest in the coexisting large amplitude nonlinear dynamic regime. Instead, the higher amplitude of nonlinear displacement exhibits a more continuous change in amplitude as the temperature difference changes around the buckling temperature, suggesting a close underlying connection between pre- and post-buckling large amplitude dynamics of the beam.

The analytical displacement amplitude shown in Fig. 4 indicates that the pre- and postbuckled “thermal resonances” exhibit hardening characteristics. The hardening of the prebuckled thermal resonance, originally around 0.87 K, extends into the postbuckled regime. For the thermomechanically buckled beam, such large amplitude deflections become the snap-through dynamic. Such phenomenon is observed in Fig. 4 by noting that the intersection of the hardening of the prebuckled resonance and the dashed line corresponding to 2.2 K manifests as snap-through at 16 Hz in the frequency sweep of Fig. 3(d). The snap-through dynamic is marked in both Fig. 3(d) and Fig. 4 by circles. The emergence of snap-through from the large amplitude vibrations of the prebuckled beam is consistent with the notion that snap-through is a nonresonant phenomenon and requires increased amplitudes of input force to sustain the behavior at higher frequencies. Furthermore, these findings suggest a fundamental connection between the high amplitude dynamic of the prebuckled structure and the interwell regime of the subsequent multistable structure. Such an underlying connection may also help elucidate the way by which the snap-through dynamics emerge from and connect the separate intrawell regimes in the impedance results of Fig. 3.

The potential energy of the unforced beam is utilized to uncover the connection between the high amplitude dynamics of the pre- and postbuckled structure and to further illuminate the cusp denoted by the triangle in Fig. 4. The expression for this potential energy is given by Eq. (33), which considers the energy resulting from mechanical strain, the potential energy of the induced thermal axial load PT, and from the induced thermal moment MT. The total potential energy of the unforced beam is shown in Fig. 5(a). 
U=UM+UT=120LEI[2wx2{1+12(wx)2}]2dx+120L[Eαbh212ΔT2wx2{1+12(wx)2}Eαbh2ΔT(wx)2]dx
(33)

The cusp in Fig. 4 occurs when the beam transitions from monostable to bistable at 1.76 K, denoted in Fig. 5 as ΔTcr. Since the current analytical formulation does not account for initial beam curvature, the single equilibrium of the monostable beam occurs at zero transverse beam displacement, independent of temperature. As the temperature is increased beyond this value of thermomechanical buckling, the potential well of the beam is seen in Fig. 5(a) to change from possessing one minimum (in the monostable case) to two minima (in the bistable case). This is a change in the local topology at the bottom of the potential well and corresponds to a pitchfork bifurcation. For snap-through oscillations of the beam to be induced, the input mechanical and thermal energies must be sufficiently great for the beam to transition between the two local minima of potential energy by exceeding the potential energy difference between wells in Fig. 5(a).

The curvature of each potential energy component with respect to the transverse beam displacement is presented in Fig. 5(b). The mechanical component of the potential energy curvature is the mostly parabolic surface with grid lines, while the thermal component is the plane originating at zero curvature when ΔT is zero. The negative of the thermal component of potential energy curvature is plotted in order to compare the relative sizes of each component. In Fig. 5(b), regions where the mechanical curvature is greater than the thermal curvature indicate locations where the curvature of the potential well in Fig. 5(a) is always positive. Conversely, when the thermal curvature exceeds the mechanical curvature, the system may exhibit multiple states of stable equilibrium by having a nonconvex region of potential energy. The fact that the thermal component of curvature changes with increasing ΔT while the mechanical component remains constant illustrates that the mechanical strain energy is independent of temperature change, and clarifies that the thermal aspects of the potential energy depicted in Fig. 5(a) lead to the observed changes in the dynamic behavior of the beam. Additionally, it is seen that at ΔTcr, the total curvature of the potential energy at the single static equilibrium is zero since the curvatures of the mechanical and thermal components of potential energy are equal and opposite at this point. This produces a stable static equilibrium exhibiting critical slowing down that becomes unstable for increasing ΔT as the pitchfork bifurcation occurs [39]. Beyond ΔTcr, the intersection of the mechanical and thermal curvatures in Fig. 5(b) denotes inflection points of the potential energy surface in Fig. 5(a). The stable equilibria of the bistable beam exist for transverse beam displacements larger than the transverse beam displacements of these inflection points, since for ΔT larger than ΔTcr, the inflection points separate the stable and unstable equilibria.

Influence of Temperature Variation on Mechanical and Thermomechanical Impedances.

Figure 6 presents the (a) mechanical and (b) thermomechanical impedances corresponding to the results of Fig. 4 to examine the influence of change in maximum transverse temperature elevation ΔT. The impedances in Fig. 6 are plotted parametrically where the changing parameter is ΔT. While both mechanical and thermomechanical impedances exhibit similar trends, the main qualitative difference is that for increasing ΔT starting at 0 K, the thermomechanical impedance results begin at the origin of the complex impedance plane because the beam is at the reference temperature Tref at this point. The cusp discussed in Sec. 5.1 is evident in the linear and nonlinear results presented in Fig. 6, denoted by triangles. By the nonsmooth character of the cusp, a clear transition to a qualitatively different dynamic regime is suggested. Furthermore, the snap-through regime denoted by the circular marker in Fig. 6 is in close proximity to the intrawell impedance results near the bifurcation designated by the star markers in Fig. 6 and close to the analytical impedance results after the cusp occurs.

The proximity of interwell and intrawell impedances suggests that the hardening of the prebuckled thermal resonance is the origin of the separation that occurs between the intrawell regimes in the impedance plots of Fig. 3, which was inferred in Sec. 5.1. The snap-through regime denoted by gray arrows in Fig. 6 terminates when the reactances of mechanical and thermomechanical impedances vanish. This behavior qualitatively agrees with trends in the snap-through dynamics of the frequency sweeps presented in Sec. 4. Thus, vanishing reactance is seen to predict the cessation of snap-through dynamics for both increasing excitation frequency ω and temperature elevation ΔT.

Since the base acceleration amplitude and excitation frequency are held constant, the trends exhibited in the mechanical impedance of Fig. 6(a) are due only to variations in the steady-state velocity amplitude. Velocity changes are a direct result of the variation in temperature, which is shown in Fig. 5 to produce large differences in the shape of the potential energy well of the beam. In turn, the shape of the potential energy at static equilibrium determines how dynamic energy is transferred through the system and governs the sustainability of snap-through dynamics. The same trend of vanishing reactance is exhibited in the mechanical and thermomechanical impedances of Fig. 6 due to the influence that elevated temperature exerts on the velocity amplitude. Consequently, vanishing reactance may be used to forecast the termination of snap-through regimes. Thus, monitoring either force or temperature may provide such predictive capacity in applications with nonconstant temperature. The choice between a force or temperature measuring device would depend upon the ease of implementation in a particular application.

Conclusions

This research formulates an analytical model for a flat, thermomechanically loaded beam. Analytical predictions to the thermal and structural responses of the beam are obtained by equivalent linearization techniques with a multimode expansion in mechanical and thermal coordinates. Comparisons with comprehensive experimental results verify the efficacy of the analytical approach to characterize thermomechanical coupling influences on the nonlinear dynamics of a thermally buckled, harmonically excited beam. The analytical formulation also facilitates the study of temperature influences on structural displacement amplitudes corresponding to each dynamic regime of nonlinear behavior. An examination of mechanical and thermomechanical impedances indicates that bifurcations between intrawell dynamics correspond to jumps across zero reactance for distributed parameter systems. Additionally, the saddle-node bifurcation that terminates snap-through regimes is seen to exhibit vanishing mechanical and thermomechanical reactances for variation in excitation frequency and in temperature. In addition, fundamental and continuous connections are discovered in the thermomechanical impedance planes between nonlinear low amplitude dynamics of the postbuckled beam and the energetic snap-through vibration that are otherwise discontinuous by studying displacement amplitudes alone. The analytical methodology articulated here may be extended to encompass a broader class of nonlinear structural systems subjected to thermomechanical coupling. By way of direct connections to experimentally measurable impedance values, this research forms a future tool that may be used to forecast extreme dynamics of structures operating in combined thermal and mechanical loading environments.

Funding Data

  • Air Force Research Laboratory (Summer Faculty Fellowship).

  • Haythornthwaite Foundation (Young Investigator Award).

  • U.S. Department of Defense (Science, Mathematics, and Research for Transformation (SMART) Scholarship-for-Service Progam).

References

References
1.
Barnes
,
J. R.
,
Stephenson
,
R. J.
,
Welland
,
M. E.
,
Gerber
,
C.
, and
Gimzewski
,
J. K.
,
1994
, “
Photothermal Spectroscopy With Femtojoule Sensitivity Using a Micromechanical Device
,”
Nature
,
372
(
6501
), pp.
79
81
.
2.
Brinson
,
L. C.
,
1993
, “
One-Dimensional Constitutive Behavior of Shape Memory Alloys: Thermomechanical Derivation With Non-Constant Material Functions and Redefined Martensite Internal Variable
,”
J. Intell. Mater. Syst. Struct.
,
4
(
2
), pp.
229
242
.
3.
Leyva
,
I. A.
,
2017
, “
The Relentless Pursuit of Hypersonic Flight
,”
Phys. Today
,
70
(
11
), pp.
31
36
.
4.
Galea
,
S. C. P.
, and
White
,
R. G.
,
1993
, “
The Effect of Temperature on the Natural Frequencies and Acoustically Induced Strains in CFRP Plates
,”
J. Sound Vib.
,
164
(
3
), pp.
399
424
.
5.
Thornton
,
E. A.
,
1993
, “
Thermal Buckling of Plates and Shells
,”
ASME Appl. Mech. Rev.
,
46
(
10
), pp.
485
506
.
6.
Noda
,
N.
,
1991
, “
Thermal Stresses in Materials With Temperature-Dependent Properties
,”
ASME Appl. Mech. Rev.
,
44
(
9
), pp.
383
397
.
7.
Lemlich
,
R.
,
1955
, “
Effect of Vibration on Natural Convective Heat Transfer
,”
Ind. Eng. Chem.
,
47
(
6
), pp.
1175
1180
.
8.
Lemlich
,
R.
, and
Rao
,
M. A.
,
1965
, “
The Effect of Transverse Vibration on Free Convection From a Horizontal Cylinder
,”
Int. J. Heat Mass Transfer
,
8
(
1
), pp.
27
33
.
9.
Culler
,
A. J.
, and
McNamara
,
J. J.
,
2011
, “
Impact of Fluid-Thermal-Structural Coupling on Response Prediction of Hypersonic Skin Panels
,”
AIAA J.
,
49
(
11
), pp.
2393
2406
.
10.
Culler
,
A. J.
, and
McNamara
,
J. J.
,
2010
, “
Studies on Fluid-Thermal-Structural Coupling for Aerothermoelasticity in Hypersonic Flow
,”
AIAA J.
,
48
(
8
), pp.
1721
1738
.
11.
Miller
,
B. A.
, and
McNamara
,
J. J.
,
2015
, “
Time-Marching Considerations for Response Prediction of Structures in Hypersonic Flows
,”
AIAA J.
,
53
(
10
), pp.
3028
3038
.
12.
Miller
,
B. A.
,
McNamara
,
J. J.
,
Spottswood
,
S. M.
, and
Culler
,
A. J.
,
2011
, “
The Impact of Flow Induced Loads on Snap-Through Behavior of Acoustically Excited, Thermally Buckled Panels
,”
J. Sound Vib.
,
330
(
23
), pp.
5736
5752
.
13.
Daneshjo
,
K.
, and
Ramezani
,
M.
,
2004
, “
Classical Coupled Thermoelasticity in Laminated Composite Plates Based on Third-Order Shear Deformation Theory
,”
Compos. Struct.
,
64
(
3–4
), pp.
369
375
.
14.
Carrera
,
E.
,
Boscolo
,
M.
, and
Robaldo
,
A.
,
2007
, “
Hierarchic Multilayered Plate Elements for Coupled Multifield Problems of Piezoelectric Adaptive Structures: Formulation and Numerical Assessment
,”
Arch. Comput. Methods Eng.
,
14
(
4
), pp.
383
430
.
15.
Matney
,
A.
,
Mignolet
,
M. P.
,
Culler
,
A. J.
,
McNamara
,
J. J.
, and
Spottswood
,
S. M.
,
2015
, “
Panel Response Prediction Through Reduced Order Models With Application to Hypersonic Aircraft
,”
AIAA
paper No. AIAA 2015-1630.
16.
Perez
,
R.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2011
, “
Nonlinear Reduced-Order Models for Thermoelastodynamic Response of Isotropic and Functionally Graded Panels
,”
AIAA J.
,
49
(
3
), pp.
630
641
.
17.
Settimi
,
V.
,
Saetta
,
E.
, and
Rega
,
G.
,
2017
, “
Local and Global Nonlinear Dynamics of Thermomechanically Coupled Composite Plates in Passive Thermal Regime
,”
Nonlinear Dyn.
(epub).
18.
Pal
,
M. C.
,
1970
, “
Large Amplitude Free Vibration of Circular Plates Subjected to Aerodynamic Heating
,”
Int. J. Solids Struct.
,
6
(
3
), pp.
301
313
.
19.
Lee
,
J.
,
1993
, “
Large-Amplitude Plate Vibration in an Elevated Thermal Environment
,”
ASME Appl. Mech. Rev.
,
46
(
11S
), pp.
S242
S254
.
20.
Chang
,
W. P.
, and
Wan
,
S. M.
,
1986
, “
Thermomechanically Coupled Non-Linear Vibration of Plates
,”
Int. J. Non-Linear Mech.
,
21
(
5
), pp.
375
389
.
21.
Chang
,
W. P.
, and
Jen
,
S.
,
1986
, “
Nonlinear Free Vibration of Heated Orthotropic Rectangular Plates
,”
Int. J. Solids Struct.
,
22
(
3
), pp.
267
281
.
22.
Shu
,
X.
,
Zhang
,
X.
, and
Zhang
,
J.
,
2000
, “
Thermoelastic Free Vibrations of Clamped Circular Plate
,”
Appl. Math. Mech.
,
21
(
6
), pp.
715
724
.
23.
Harne
,
R. L.
, and
Goodpaster
,
B. A.
,
2018
, “
Impedance Measures in Analysis and Characterization of Multistable Structures Subjected to Harmonic Excitation
,”
Mech. Syst. Signal Process.
,
98
, pp.
78
90
.
24.
Kovacic
,
I.
,
Brennan
,
M. J.
, and
Lineton
,
B.
,
2008
, “
On the Resonance Response of an Asymmetric Duffing Oscillator
,”
Int. J. Non-Linear Mech.
,
43
(
9
), pp.
858
867
.
25.
Amabili
,
M.
,
2006
, “
Theory and Experiments for Large-Amplitude Vibrations of Rectangular Plates With Geometric Imperfections
,”
J. Sound Vib.
,
291
(
3–5
), pp.
539
565
.
26.
Yamaki
,
N.
,
Otomo
,
K.
, and
Chiba
,
M.
,
1981
, “
Non-Linear Vibrations of a Clamped Circular Plate With Initial Deflection and Initial Edge Displacement—Part 1: Theory
,”
J. Sound Vib.
,
79
(
1
), pp.
23
42
.
27.
Yamaki
,
N.
,
Otomo
,
K.
, and
Chiba
,
M.
,
1981
, “
Non-Linear Vibrations of a Clamped Circular Plate With Initial Deflection and Initial Edge Displacement—Part II: Experiment
,”
J. Sound Vib.
,
79
(
1
), pp.
43
59
.
28.
Harne
,
R. L.
, and
Wang
,
K. W.
,
2017
,
Harnessing Bistable Structural Dynamics: For Vibration Control, Energy Harvesting and Sensing
,
Wiley
,
Chichester, UK
.
29.
Hodges
,
D. H.
,
1984
, “
Proper Definition of Curvature in Nonlinear Beam Kinematics
,”
AIAA J.
,
22
(
12
), pp.
1825
1827
.
30.
Seif
,
M.
,
Main
,
J.
,
Weigand
,
J.
,
Sadek
,
F.
,
Choe
,
L.
,
Zhang
,
C.
,
Gross
,
J.
,
Luecke
,
W.
, and
McColskey
,
D.
,
2016
, “
Temperature-Dependent Material Modeling for Structural Steels: Formulation and Application
,” National Institute of Standards and Technology, Gaithersburg, MD, NIST Tech. Note No.
1907
.https://nvlpubs.nist.gov/nistpubs/TechnicalNotes/NIST.TN.1907.pdf
31.
Tauchert
,
T. R.
,
1991
, “
Thermally Induced Flexure, Buckling, and Vibration of Plates
,”
ASME Appl. Mech. Rev.
,
44
(
8
), pp.
347
360
.
32.
Nowacki
,
W.
,
1975
,
Dynamic Problems of Thermoelasticity
,
Noordhoff International Publishing
,
Leyden, The Netherlands
.
33.
Hetnarski
,
R. B.
, and
Eslami
,
M. R.
,
2009
,
Thermal Stresses—Advanced Theory and Applications
,
Springer Science+Business Media
, New York.
34.
Lorencini
,
E.
,
1970
, “
Heat Conduction With a Temperature-Dependent Thermal Conductivity Coefficient
,”
J. Eng. Phys.
,
19
, pp.
1548
1554
.https://link.springer.com/article/10.1007/BF00825355#citeas
35.
Beslin
,
O.
, and
Nicolas
,
J.
,
1997
, “
A Hierarchical Functions Set for Predicting Very High Order Plate Bending Modes With Any Boundary Conditions
,”
J. Sound Vib.
,
202
(
5
), pp.
633
655
.
37.
Goodpaster
,
B. A.
, and
Harne
,
R. L.
,
2018
, “
An Experimental Characterization of the Impedance and Spectral Content of Multistable Structural Responses During Dynamic Bifurcations
,”
ASME J. Vib. Acoust.
,
140
(5), p. 051009.
38.
Szemplińska-Stupnicka
,
W.
,
1988
, “
Bifurcations of Harmonic Solution Leading to Chaotic Motion in the Softening Type Duffing's Oscillator
,”
Int. J. Non-Linear Mech.
,
23
(
4
), pp.
257
277
.
39.
Strogatz
,
S. H.
,
1994
,
Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering
,
Addison-Wesley
,
Reading, MA
.