Study of Damping of Bare and Encased Steel I-Beams Using the Thermoelastic Model

: Steel I-beams are a fundamental structural component in civil construction. They are one of the main load-bearing components in a building that must withstand both the structure and any incoming external perturbations, such as seismic events. To avoid damage to the structure, the building must be designed to dissipate the maximum amount of energy possible. One way energy can be dissipated is through internal or structural damping, of which thermoelasticity is one of the causes, especially in low-frequency harmonic excitations. The main goal of this study is to analyze the amount of damping in an I-beam generated by thermoelasticity and when encased in a Portland cement concrete layer, using a Finite Element model. It was found that, due to the geometry of the I-Beam, the damping coefﬁcient as a function of frequency has two local maxima, as opposed to the traditional single maximum in rectangular beams. Encasing an I-beam in a concrete layer decreases the overall damping. While the extra coating protects the beam, the reduction in damping leads to a lower energy dissipation rate and higher vibration amplitudes


Introduction
Seismic events are low-frequency, high-energy events that can impart heavy structural damage on structures [1][2][3][4].The energy propagating through the soil into the structure can be partially dissipated rather than converted into deformation.Several strategies exist to maximize damping, either through deliberate design of stationary structures and materials [5][6][7][8][9], or via the employ of dynamic absorber devices (active or passive) which make use of moving parts to counter-balance incoming vibrations [10][11][12].On the other hand, natural dissipation phenomena are always present in any structure.Macroscopically, they originate from the friction between bolted or riveted parts.Microscopically, the main damping phenomenon is structural damping [13][14][15].Structural damping is an umbrella term that includes all damping sources intrinsic to a solid structure and excludes damping due to viscous interaction with a fluid or direct friction with another moving part [13].
Since the source of structural damping is not unique [13], several models have attempted to predict the dynamic behavior of a structure.One of the most well-known models is the hysteretic model [13,[16][17][18][19].This model, however, causes controversy in the scientific community.The consensus is that this model can only be used in the analysis of steady-state vibrations because of the supposed non-causality of the transient response [16,18,[20][21][22][23][24].In actuality, this model has trouble modeling the transient response due to an unstable differential equation.Therefore, a more generic model is needed to cover both regimes.Fractional derivative models are also common in the study of structural damping, especially in polymeric materials [19,[25][26][27][28].However, these are empirical models which share a significant drawback: there is no single definition for fractional derivatives.A commonly overlooked alternative is the thermoelastic model.
The thermoelastic model is derived from first principles and fully couples the elastic wave equation with the heat equation.While the origins of this model can be traced back to the 1950s [29,30], because of its complexity, only with the advent of modern computers could this model be more thoroughly exploited [30][31][32][33][34][35].This model simulates damping phenomena at the micro-scale [36][37][38][39], since that effect is more pronounced in smaller structures.However, it is also noticeable in larger structures at lower frequencies [35,40].Since seismic events tend to occur in the lower band of the frequency spectrum, thermoelasticity could contribute significantly to the overall damping.This paper compares the thermoelastic damping behavior of concrete-encased steel I-beams to that of bare beams.Steel I-beams are fundamental structural elements used in buildings and, together with columns, are their main load-bearing components.Understanding the dynamic behavior of these elements, especially in what damping is concerned, can help engineers design structures that are better prepared to handle incoming dynamic loads.While there have been several studies on thermoelasticity applied to beams or plates (either "regular" size or micro/nano-sized), most sources have focused on the study of elements with rectangular cross-sections [30][31][32][33][36][37][38][39]41,42].
The case studies presented in this work focus on beams under a bending load, and the analysis uses solid elements instead of the more common beam elements [30,32,33,39].Since solid elements are more versatile, the variation of temperature across the thickness of the beam can be studied in more detail.The main disadvantage is that solid elements are computationally more expensive since the analysis requires more elements than it otherwise would in order to capture bending behavior accurately.
The article is divided into four main sections: Introduction, Materials and Methods, Results and Discussion, and Conclusions.The thermoelastic model differential equations are presented in the Materials and Methods section.This section also covers the formulation of the 3D hexahedral element and the nuances in its use in conjunction with the thermoelastic model (e.g., boundary conditions).In the Results and Discussion section, the model validity is assessed, mesh generation and convergence studies are performed, and the results of the simulations are presented.

The Thermoelastic Model
The thermoelastic model or linear thermoelasticity dynamic model results from the complete coupling of the elastic wave equation with the heat equation and is given by (1): where u is the displacement vector (u = [u x , u y , u z ] T ), θ is the temperature variation, f is the body load and q is the body heat flux.All function quantities are functions of the Euclidean coordinates and time.The constants in the model are the material properties: density (ρ), Lamé constants (µ and λ), specific heat capacity (C p ), and thermal conductivity (k).
The Lamé constants are given by: where E is the Young modulus and ν is the Poisson coefficient.Since temperature is critical to the problem, all coefficients must be specified at the same base temperature T 0 , the point around which linearization is performed in the thermoelastic model.Finally, the coefficient γ is the coupling term and is given by: where α is the linear thermal expansion coefficient.Since Equation ( 1) is a linearized equation, it is only valid for small deformations and small temperature variations: θ T 0 1.The coupling term in the elastic wave equation is γ∇θ and results from the thermal strain in elastic materials: where θ is the strain generated by temperature variation and δ is the Kronecker Delta.The total strain, is given by: where M is the mechanical strain (stress-induced strain).The coupling term in the heat equation is γT 0 (∇ • u) and is the result of the definition of entropy: where s is the specific entropy and is the total strain.A more detailed demonstration of the coupling term in the thermal equation can be seen in [29].

Energy Flow and Balance in the Thermoelastic Model
Because it couples two fields in a single system, the thermoelastic model has some particularities regarding energy flow, especially when compared with the standard elastic wave equation.Analyzing the energy flow makes it possible to understand the phenomena behind the damping predicted by this model.A complete representation of the energy flow in the system can be seen in Figure 1.The conservative energies are the potential energy (V), kinetic energy (T ) and exergy (B), which are given by Equation ( 7): where C is the elasticity tensor, u i is the displacement vector (u) in index notation, Ω is the domain of the structure and V is the volume.Please note that the potential energy in this model is not the same as the elastic potential in traditional elasticity because it uses the complete strain and not just the mechanical strain.The only non-conservative energy in the domain is the work lost to irreversible entropy generation.One of the characteristics of the Equation ( 1) is that the damping term is not explicit, much unlike all other damping models used in mechanical vibrations.Damping occurs by generating irreversible Entropy due to the existing internal heat fluxes.The irreversible entropy is generated whenever there is a heat flux in the material, and the total energy lost this way is given by the Gouy-Stodola theorem [31]: where W S is the work lost to irreversible entropy.Internal entropy generation is not the only way the system can lose energy.Energy can also be lost through the boundaries either by imposing a fixed temperature, imposed nonzero flux, or through convection.When the system is thermally open, the entropy can flow into and out of the domain due to the generated heat fluxes.These boundary fluxes are the Work (Equation (9a)) and Heat (Equation (9b)): where f i and P i are the index notations of the body force ( f ), and the loads applied at the boundary, q and q B are the body and boundary heat fluxes, respectively, and Γ is the domain of the boundaries of the structure.
In light of the First Law of Thermodynamics, the complete energy balance is given by the following equation: where ∆U is the variation of the internal energy.In this energy balance, the work done on the system is positive.
If the study is focused on harmonic analysis, then all conservative energies in a cycle are zero, and, as such, the Equation (10) reduces to: In this study, all thermal boundary conditions are Adiabatic, and the heat transferred (Q) is always zero.The dissipated energy can be determined either by Equation (9a) or Equation (8) since they have to be equal.However, the dissipated work is not a very good measure of the damping of a system since it is susceptible to the presence of resonances and, as such, is not generic or comparable enough.An alternative is to use the loss angle or an equivalent damping factor.
The loss angle is a quantity widely used in material science to characterize damping in viscoelastic materials and is the phase of a complex Young modulus definition: where E * is the dynamic modulus, E is the storage modulus (reduces to the Young Modulus in purely elastic materials), E is the loss modulus that quantifies the dissipated elastic energy, and φ E is the loss angle.
The loss angle is related to the hysteretic model and the hysteretic damping coefficient: where η is the hysteretic damping coefficient.
For the hysteretic model, the damping coefficient can be extracted from the work made by a load using the maximum elastic potential energy: where k is the stiffness of the spring, and X is the complex amplitude of the displacement.The denominator of Equation ( 14) is the maximum potential energy in a cycle.The equivalent hysteretic damping coefficient is also known in the literature as the quality factor (Q −1 TED ) (The TED under script refers to the thermoelastic damping.The Quality Factor can be defined for any damping) [40,43,44].

Finite Element Model
The thermoelastic model was implemented on an in-house Finite Element Analysis software in 2D quadrilateral and 3D hexahedral solid elements (linear and quadratic for 2D and linear for 3D).The software was written in C++ and Intel® Math Kernel Library performs the linear algebra routines.All elements support static, dynamic, and harmonic analysis.
For this study, a solid linear hexahedral 3D element was defined according to the Figure 2.For the given reference frame, the shape functions are given by Equation (15): where φ are the shape functions and u, v and w are the element internal coordinates that are in the interval [0, 1] ∈ R.
To determine the element matrices, the weak form of Equation ( 1) must be determined.The first step is to expand the differential operators and explicitly write the equation as a system of four differential equations: where u x , u y and u z are the displacements in the x, y and z directions, respectively.As represented in Equation ( 16), the system has three elasticity equations and one thermal equation.
The discretization process for the displacements and temperature variations is explained in detail in Appendix A.
After the discretization process is done, Equation (A4) can be joined with Equations (A11a), (A11b) and (A14) and the full discretized system is given by: where M, C and K, are the mass, velocity (the term damping matrix should be avoided because the damping in the thermoelastic model is not completely contained in this matrix) and stiffness matrices, given by: and f and B are the nodal body loads (forces and heat fluxes) and the boundary terms (Equation ( 26)), respectively.The matrices M, L and N are given by Equations (A5) and (A12).
It is important to note that while Equation ( 17) is written with second-order derivatives, the system is, in fact, a third-order system (this can be seen in the mass matrix, which is non-invertible).Also, Equation ( 17) can be applied to a single element or the entire mesh, assuming each matrix is correctly assembled.

Boundary Conditions
Boundary conditions (BC) are a fundamental component of boundary value differential equation problems, sometimes overlooked.While the boundary conditions for individual problems in Elasticity and Heat conduction are well-studied, their coupling generates new boundary conditions with unique properties.These conditions must be explicitly presented, particularly in the case of Neumann boundary conditions (Neumann BCs).
Neumann BCs are conditions proportional to the gradient of displacements (or temperature variation).In classical Elasticity, they are employed to incorporate any applied forces into the system at the boundaries of the structure.In heat problems, Neumann BCs are used to introduce external heat fluxes.
Using a 1D example for simplicity, a Neumann BC for the Elasticity problem is: where E is the material Young Modulus, A is the area of the boundary, P is the distributed load applied at the boundary and Γ is the boundary.
When using the thermoelastic model, according to Equation ( 5), the strain not only depends on the applied force (similarly to Equation ( 19)) but also on the temperature variation.As such, the BC in Equation ( 19) has to be updated to include the temperature variation, becoming: This new boundary condition is pivotal to the model as it addresses a seemingly paradoxical behavior inherent in thermoelasticity under adiabatic conditions: • When a body is stretched, it cools down, and when compressed, it heats up.

•
When heated, a body expands, and when cooled, it shrinks.
The first behavior is embedded in the differential equations themselves, whereas the second behavior is represented by the boundary condition in Equation (20).
The seemingly minor modification described above has profound implications for both the problem and the subsequent finite element model.Firstly, as evident from Equation ( 20), the thermoelastic model exhibits not Neumann boundary conditions (BCs) but rather Robin BCs.This distinction arises because the BC is a linear combination of the first derivative of a state (displacement, u x ) and a state (θ).
The second consequence pertains to the finite element implementation, compelling the analysis to consider all boundaries, even those without externally applied loads.
The third implication arises when multiple materials are employed on a mesh, as exemplified in the case of a concrete-encased beam.In classical elasticity and heat conduction problems, a finite element model mesh can accommodate multiple materials without special modification, provided that individual elements maintain a single material.This accommodation is feasible due to the continuity of forces inside the mesh, extending across material boundaries.In contrast, when transitioning to the thermoelastic model, only internal forces remain continuous.The final monomial in Equation (20) introduces discontinuity (This term behaves akin to a fictitious force, contingent upon the material properties and temperature variation).Consequently, in the thermoelastic model, internal boundaries need to be added at material interfaces within the mesh.For a 1D problem, this scenario is depicted in Figure 3.
As seen in Figure 3b, while the internal force ( f ) is in balance, the "thermal forces" are not, and an imbalance generates a net force.If the standard implementation of a finite element model is used, this imbalance is ignored, leading to incorrect results.

Implementation of the Boundary Conditions
The hexahedral thermoelastic solid element has six boundaries, identified in Figure 4.A series of boundary terms appeared when performing the integration by parts on Equation (A2).These terms are: where Γ are the boundaries according to Figure 4.
According to Equation (20), the boundaries terms defined in Equation ( 22) are also: where P x is the x coordinate of the distributed load applied on the surface of the boundary and n # x is the x coordinate of the normal vector of the given boundary.
Applying the discretizations in Equation (A3) and a similar one to Equation (A6), the discretized boundary terms are: where P x is the nodal distributed load values (the values sampled at the element nodes).
The matrices m are the boundary matrices and are similar to Equation (A5a), but only integrated at the assigned boundary: The matrices m nx are the boundary matrices with the x coordinate of the normal vector: As with body loads, the computational cost of integrating loads and temperature variations at boundaries is considerably reduced through the use of boundary matrices.However, in the context of Equation ( 25), it is crucial to exercise caution to prevent a substantial normal rotation.Should the problem result in a significant rotation of the boundary normal, it becomes imperative to update these matrices during the simulation.
Finally, the element matrices in Equations ( 24) and ( 25) are not assembled like the other element matrices.Instead, they are only assembled for elements at the structure boundaries and using only the boundary matrices corresponding to the exposed element sides.
With the boundaries defined, the boundary term, B, in Equation ( 17) is defined as: where N b is the number of boundaries in the mesh added to any internal boundary generated by a material interface.

Harmonic Solution
To get the solution of the harmonic problem, a Fourier transform must be applied to Equation (17): where ω is the angular frequency of the excitation and U, Θ, F, Q, P and Q b are the complex amplitudes of the displacements, temperature variation, body forces, body heat fluxes, boundary loads and boundary heat fluxes, respectively.The matrix K Θ contains the Robin BC elements of the thermoelastic BCs: Finally, to apply Dirichlet BCs, the system matrix must be decoupled accordingly.
To obtain the complex amplitudes (and the corresponding amplitudes and phases), Equation ( 27) is solved using the linear equation solver.

Energy Computation in the Finite Element Model
For a harmonic analysis, as stated by Equation (11), only the boundary energy fluxes and the work lost to irreversible Entropy need to be determined.To determine these quantities used the finite element model, Equations ( 8), (9a) and (9b) must be discretized and integrated for a complete cycle (0 ≤ t < 2π ω ) The discretized Gouy-Stodola theorem for a cycle is given by: where To compute the work and heat transfer in a cycle, the nodal loads must be determined first.Due to the Dirichlet BCs (direct imposition of values of DoFs), reaction loads (forces and heat fluxes) are generated in response.To determine these loads, the input vector (the right side of Equation ( 27)) has to be recomputed: The quantities with the "node" subscript are the nodal loads, and the system matrix is the full matrix without decoupling.
With the nodal loads, the work and heat in a cycle are given by: Note that for Equations ( 29) and (31b), if ω = 0, both equations result in an indeterminate form.This indeterminate form can be lifted by computing the limits, which always yields zero for both equations.
Finally, the equivalent hysteretic damping coefficient is determined by:

Geometries and Materials
In this work, two geometries are studied, a symmetrical ASTM A6 wide flange I-beam (W150 × 100 × 24.0, Figure 5a and Table 1) bare or completely encased (Figure 5b).The length of all beams in this study is 1.6 m, guaranteeing a length/height ratio 10.Two thicknesses of the encasing (t O ) were used: 35 mm and 17.5 mm.
Apart from one case study, where the beam is made of copper, all case studies use standard steel (for more information on the case studies, please refer to Section 3.3).For the casing, on the encased studies, it was used an average Portland cement concrete (28 days cure).The material properties can be seen in Table 2.For all properties, the base temperature, T 0 , is 293.15K.

Validation Study
During the literature review, no published results were found for damping I-beams.The thermoelastic model has been extensively validated previously; for example, Zener compared the predicted maximum quality factor for different rectangular beams [43].More recently, Zuo et al. compared the numerical simulation of a double-clamped quartz micro-beam, showing a good fit with the experimental data [41].
A study using a double-clamped rectangular beam in [40] was replicated to validate the developed element.The validation study is also compared with base benchmarks using similar analyses done by Zener [43] and Lifshitz and Roukes [44].
The geometry used in the validation is a double-clamped rectangular beam with a length (L) of 90 µm and a width (b) of 4.5 µm; the thickness (h) varies throughout the study, Figure 6.The objective of the authors in Parayil et al. [40] was to study the evolution of the quality factor (the equivalent hysteretic damping coefficient) at the first bending mode, with an increasing beam thickness by using different implementations of a Timoshenkobased thermoelastic model.The material used in the simulations was silicon with the properties in Table 3.The comparative results are shown in Figure 7. Comparative study of the evolution of the quality factor (Q −1 TED ) with the beam thickness (h) for the first bending mode, adapted from [40].

Solid Element
The present study (SLD-FEM) matches the Thimoshenko-based model (TBT-A1) results for all thicknesses.The study also fits well with the results obtained by Zener (EBT-Z) and Lifshitz (EBT-LR) for more slender beams, which is expected since both analyses rely on an Euler-Bernoulli beam model.
The mesh used in the SLD-FEM study varies since the thickness also changes.After a convergence study for each case, the final mesh has 220 elements across the length and 19 elements across the width.The number of elements across the thickness, N h , is given by the following expression: where h is the thickness of the beam.

Mesh Convergence
Only two geometries were used in this study: the bare I-beam and the encased I-beam, and, as such, only two meshes were used.The implementation of the solid element was validated against standard cases for the elastic model.
To determine the mesh parameters, a mesh convergence study was done for both meshes.Since, in this article, only harmonic simulations are used, the convergence of the first two natural frequencies was chosen as the convergence criterion.
For simplicity, the mesh parameters were defined as a function of the number of elements across the length of the beam (N L ): where L is the beam length and N t w , N h w , N w b , N w t , N h b , N h t are the number of elements across the web thickness, web height, bottom flange width, top flange width, bottom flange thickness and top flange thickness, respectively.These equations were established to have elements with a maximum proportion of 5:1 (length relative to any other dimension).
For the encased I-beam, the number of elements in the casing thickness (t O ) is given by: The mesh convergence for the bare I-beam can be seen in Table 4 and graphically in Figure 8a.The selected mesh was the last, with 240 elements across the length, Figure 8b.The second mesh convergence study was done for an encased W150 × 100 × 24.0 I-beam, with a casing thickness (t O ) of 35 mm.The convergence data can be seen in Table 5 and graphically in Figure 9a.It was considered that the mesh converged at 140 elements across the length, Figure 9b.

Case Studies
For this article, nine different case studies were used.The description of all cases can be seen in Table 6.All case studies use a beam with a 10:1 ratio to the total height of the I-beam (10 H), and all cases use a 1 N/m distributed load applied at x = 0.25 m, Figure 10.The position of the force was chosen to avoid nodes in the bending mode shapes, guaranteeing that the maximum vibration modes are excited.Regarding the BCs, all case studies feature double-clamped beams (at x = 0 and x = 1.6 m), and all boundaries are Adiabatic.
Adiabatic boundaries signify that the model enforces no heat flux at any boundary.Other potential boundary conditions include the imposition of temperature (Dirichlet), convection (Robin), and radiation (non-linear BC), each influencing the damping characteristics of the beam.The imposition of a constant temperature is unrealistic, assuming surfaces maintain a constant time-controlled temperature.Physically, this would require the application of temperature-controlled heaters on all surfaces, essentially creating a thermal clamp.Conversely, convection and radiation boundary conditions are more realistic.Convection represents heat transfer between a solid and a fluid (or vice versa) and is modeled by Newton's law of cooling.However, in terms of damping, this law has no impact as it assumes the surrounding fluid is inviscid, resulting in no energy dissipation under harmonic excitation [35].Realistic convection studies necessitate modeling the structure-fluid interaction, which is computationally expensive and minimally impactful since natural air convection is overshadowed by conduction in this study.Consequently, it is approximated by imposing adiabatic BCs.The energy lost through radiation is even smaller, particularly considering that beams in buildings are typically not directly exposed to solar radiation, and inter-surface radiation indoors is negligible.
In the case of encased beams, case studies SI018 and SI019 mirror SI009 and SI012, respectively, but employ a single theoretical volume-mixed material instead of the two material phases.This theoretical material results from a volume-weighted average of steel (7.64%) and Portland cement concrete (92.36%) properties.These studies aim to discern the influence of geometric material phase distribution on results, determining if one-dimensional beam elements suffice.If volume-mix results align with two-phase case studies, solid elements become unnecessary, and computationally less expensive beam elements are preferable.
The concrete casing thickness of 35 mm adheres to standard civil construction practices and represents the maximum thickness allowed by the Eurocode.

Harmonic Behavior of I-Beams
The thermoelastic model is more prominent in low frequencies; as such, all case studies were run in two frequency sweeps: a long-range "low" resolution sweep from 0 to 5500 rad/s with a resolution of 6 rad/s and a short-range variable high-resolution frequency sweep from 0 to 30 rad/s.The first sweep was done to understand how the model behaves in a broader range and compare it with the other models.The second sweep is more focused on representing the evolution of the equivalent damping coefficient with frequency.
The long-range frequency sweep Bode diagram for the first three case studies is shown in Figure 11.The Bode diagram shown in Figure 11 was taken on the line of the applied load (Figure 10) at the center of the beam (y direction).The three case studies show the comparison between the standard elastic model (SI001), the hysteretic model (SI002, with η = 0.02), and the thermoelastic model (SI003), Table 6.All models show nearly coincident resonant frequencies but very different behaviors in damping (more evidently in the phase).The first two resonant frequencies for the three models are shown in Table 7.
The Bode diagram, depicted in Figure 11, was captured along the line of the applied load (Figure 10) at the beam center in the y-direction.The comparison across three case studies, as illustrated in Table 6, includes the standard elastic model (SI001), the hysteretic model (SI002) with a damping ratio of η = 0.02, and the thermoelastic model (SI003).Although all models exhibit nearly coincident resonant frequencies, they manifest distinct damping behaviors, particularly evident in the phase.The first two resonant frequencies for each model are summarized in Table 7.In terms of temperature distribution, a reasonable interpretation can be made using the first mode shape of the beam, Figure 12a (the values at the axis do not have the physical meaning as the actual displacements depend on the initial conditions).The cross-section at the anti-node can be seen in Figure 12b.As expected, the beam cools down on the traction side and heats at the compression side of the beam.The temperature variation is smooth across the cross-section, except for some concentration at the interfaces between the flanges and the web.
Using Equation ( 32), the behavior of the equivalent damping coefficient for the I-beam thermoelastic case studies (SI003 and SI004) is shown on Figure 13.
The overall behavior closely aligns with anticipated thermoelastic patterns, as demonstrated in previous works [34,35,45].As anticipated, the copper beam exhibits higher damping than steel, primarily attributable to the superior thermal conductivity of copper relative to steel.
The occurrence of maximum damping corresponds to excitation frequencies aligning with the inverse of a time constant in the thermal equation [40].However, as depicted in Figure 13, two distinct maxima are observable, particularly pronounced in the copper case study.One plausible explanation for this phenomenon lies in the use of 3D solid elements on a relatively thick beam, revealing more complex behaviors compared to conventional one-dimensional beam elements found in literature [32,34,45].Additionally, the I-beam shape introduces two thicknesses-the web and flange-with comparable dimensions, bringing the first two thermal time constants into closer proximity.

Harmonic Behavior of Encased I-Beams
Two frequency sweeps were also done in the encased I-beam simulations (case studies SI009 to SI012, Table 6).However, the "low" resolution sweep was done with a different range: 0 to 2320 rad/s.The simulation of the thermoelastic model is considerably more expensive, computationally, than the elastic model.Since this work focuses on lower frequencies, a range that only captures the first mode of vibration was chosen.The first frequency sweep Bode diagram is shown in Figure 14.The first resonant frequencies for the case studies in Figure 14 are shown on Table 8.As evident from Tables 7 and 8, the resonant frequency of the thermoelastic model is consistently higher than the natural frequency of the elastic model.This observation may initially seem counterintuitive, as damping typically decreases the resonant frequency in comparison to an undamped system.However, the thermoelastic model introduces unique characteristics related to the Heat equation.
Recalling the energy balance in a thermoelastic system, illustrated in Figure 1 and governed by Equations ( 7) and ( 8), the Heat equation introduces, alongside damping, the concept of Exergy.Exergy serves as an additional energy storage mechanism, quantifying the potential of temperature imbalances.Since the squared natural frequencies of a system are proportional to the ratio between potential energies and kinetic energy, an increase in the former results in higher natural frequencies.This phenomenon arises because the Heat equation lacks a "thermal inertia" term.
The absence of "thermal inertia" has been a topic of discussion in the scientific community [46][47][48][49], primarily because Fourier's Law predicts instantaneous excitation propagation throughout the domain.While acceptable in Classical Physics, this contradicts Modern Physics, which asserts that no phenomena can travel faster than the speed of light.
In the Classical Thermodynamics context, Fourier's Law remains a viable approximation because thermal conduction phenomena are "very slow" relative to the speed of a hypothetical thermal wave, and the "thermal inertia" is likely orders of magnitude smaller compared to the conduction term, allowing it to be safely disregarded.However, when coupling the Heat equation with other equations predicting faster phenomena, such as the Wave Equation, neglecting "thermal inertia" could lead to deviations.Unfortunately, reliable values for proposed "thermal inertias" for the materials used could not be found during the research for this work.
The distribution of the temperature variation in the cross-section of the encased beam (for the case studies SI012 and SI019, at the resonance) can be seen in Figure 15.
By comparing Figure 15b,d, the impact of two material phases on temperature distributions becomes evident.The temperature profile for the volume mix closely resembles that of case study SI003 (Figure 12b), displaying a gradual transition from colder regions (under tension) to hotter regions (under compression).However, when examining the results of case study SI012, where distinct material phases are present, the temperature distribution is no longer smooth.In this instance, a noticeable disparity emerges between the temperature profiles of the two material phases.Specifically, the central I-beam exhibits a higher temperature variation than the surrounding concrete casing.This discrepancy emphasizes that the material phase distribution is a crucial factor that cannot be overlooked, rendering the use of one-dimensional beam elements inadequate.
The computation of the evolution of the equivalent damping coefficient for the encased case studies is illustrated in the plot presented in Figure 16.
The first observation one can take from Figure 16 is the overall decrease in the frequency of the maximum point.The damping values also decrease significantly in all frequencies.A better comparison between all case studies can be seen in Figure 17 by plotting the damping coefficients in a log-log scale.The decrease in the thermoelastic damping coefficient between the cased and bare I-beams can be attributed to two key characteristics of the model and geometry.Firstly, thermoelastic damping is typically higher in thin geometries where temperature variations are more pronounced, resulting in higher energy losses.The addition of a thick casing substantially increases the overall bulk of the beam, mitigating the impact of the I-beam's thin features.Secondly, heat conduction is more effective than convection (assumed to be nearly adiabatic in this study).Despite having a smaller thermal conduction value than steel, the concrete layer diffuses temperature variations across its section.This diffusion significantly diminishes temperature gradients, following the Gouy-Stodola theorem (Equation ( 8)), thereby reducing the energy lost to irreversible temperature generation.When halving the casing thickness (as observed in case study SI017), the damping decrease is less pronounced, and within a certain frequency interval, it becomes comparable to that of the bare I-beam.

Conclusions
This study primarily aimed to investigate the damping of bare and concrete-encased I-beams due to the inherent thermoelastic effect.The results obtained from simulations using harmonic Finite Element Analysis reveal the following:

•
The thermoelastic model predicts an increase in the resonant frequencies, potentially attributable to the absence of "thermal inertia" in Fourier's Law.

•
A steel I-beam exhibits a maximum damping coefficient of 3.3 × 10 −4 , and due to the cross-section geometry of the beam, the damping evolution with frequency demonstrates two local maxima.

•
The addition of a concrete casing to the I-beam reduces the maximum damping coefficient (approximately 2 × 10 −4 ), indicating that the beam dissipates less energy due to thermoelasticity compared to the bare configuration.

•
While the added layer protects the beam, the decreased damping value results in less energy dissipation and a slower reduction in vibration amplitude.

•
Given the distinct behavioral differences across the beam cross-section between the two material phases in case study SI012, traditional one-dimensional beam elements should not be employed, as they would oversimplify the problem.
This study prompts several questions for future research.Foremost is the imperative to experimentally validate the obtained results under controlled conditions.The model developed for this study can assist in designing the experiment and selecting suitable sensors and actuators.Another question to be explored in future work is the increase in resonant frequencies due to the absence of "thermal inertia" in the heat equation.The speed of the hypothetical thermal wave (or "second sound", as termed by some authors) should be studied for the materials used in this investigation to refine and update the thermoelastic model.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations and List of Symbols
The following abbreviations are used in this manuscript: Up to this point, the equations are just the rewriting of the elasticity equation.To apply the element definition, the following discretization must be done: Replacing the discretized displacements and temperature variation into Equation (A2), the discretized equation is obtained: where M, L xx , L xy , L xz and N x are the element matrices.These matrices depend only on the geometry of the element and are given by the following equations (in index notation where i are the matrix rows and j are the matrix columns): L xy ij = Ω ∂φ i ∂x (u, v, w) ∂φ j ∂y (u, v, w) J dudvdw, (A5c) where J is the Jacobian of the transformation from the element reference frame, (u, v, w), to the global cartesian reference frame, (x, y, z).To be more accurate, it should be used the absolute value of the determinant.However, for the finite element analysis to yield accurate results, all elements must be "proper" or regular.In this case, the Jacobian is always positive.The matrices M, L xx are symmetric, and L yx = L T xy and L zx = L T xz .The last integral to be handled is the one at the right side of the Equation (A4).This integral is numerically evaluated for each applied force in standard implementations of Finite Element models.While usable in the case of a static analysis, where the system is only solved once, this method can lead to very high numerical costs in the cases of multiple solutions that happen in dynamic analysis or harmonic frequency sweeps.To reduce the numerical charge, a similar discretization to the one shown in Equation (A3) can also be applied to the body force: f x (x, y, z, t) ≈ ∑ i φ i (x, y, z) f x i (t).
(A6) While this discretization is only an approximation, if the applied load is constant or linear, it will yield the same results.For other functions, the approximation is as good as the approximation of the displacements.Applying the discretization of the force, the integral becomes: where f x is the value of the body force (in the x direction) at each node of the element.
For the elasticity equations in the y and z direction, finding the weak formulation and subsequent discretization is similar to the x direction.
Multiplying both equations by the test function and integrating them in the domain, the following equations are obtained: The discretizations in Equation (A3) can be expanded to include the new discretizations in the following equations: üy (x, y, z, t) = ∑ i φ i (x, y, z) üy i (t), üz (x, y, z, t) = ∑ i φ i (x, y, z) üz i (t), f y (x, y, z, t) = ∑ i φ i (x, y, z) f y i (t), f z (x, y, z, t) = ∑ i φ i (x, y, z) f z i (t), (A10) and can be applied to Equation (A9), to obtain the discretized equations in y and z directions, resulting in Equations (A11a) and (A11b

Figure 1 .
Figure 1.Energy flow in the thermoelastic model.The wave equation originates energies marked in blue, and energies marked in orange are originated by the heat equation.

Figure 3 .
(a)  Interface between two 1D elements with the same material (marked element 1 and 2), (b) interface between two 1D elements with different materials.

Figure 4 .
Figure 4. Boundary identification in the linear solid hexahedral thermoelastic element.

4 Figure 7 .
Figure 7. Comparative study of the evolution of the quality factor (Q −1TED ) with the beam thickness (h) for the first bending mode, adapted from[40].

Figure 8 .
Figure 8.(a) Convergence of the first two natural frequencies for the W150 × 100 × 24.0 beam ("•" represents the convergence of the first mode and " " represents the convergence of the second mode), (b) mesh for the W150 I-beam with 240 elements across the length.

Figure 9 .
Figure 9. (a) Convergence of the first two natural frequencies for the encased W150 × 100 × 24.0 I-beam ("•" represents the convergence of the first mode and " " represents the convergence of the second mode), (b) mesh for the encased W150 I-beam with 140 elements across the length.

Figure 10 .
Figure10.Amplitude e position of the applied distributed load.

Figure 11 .
Figure 11.Bode diagram of the long-range frequency sweep of case studies 1 to 3, measured at the center of the distributed load (point: (0.25, 0, 0.16)).

Figure 12 .
Figure 12.(a) First mode shape of the SI003 case study, (b) mid beam cross-section showing the temperature variation across the section.

Figure 15 .Figure 16 .
Figure 15.(a) First mode shape of the case study SI012, (b) mid beam cross-section showing the temperature variation across the section for SI012, (c) first mode shape of the case study SI019, (d) mid beam cross-section showing the temperature variation across the section for SI019.

Funding:
This research received no external funding Data Availability Statement: Data are contained within the article.

Table 1 .
Dimensions of the section of a W150 × 100 × 24.0 wide flange beam (all dimensions in mm).

Table 2 .
Properties of the materials used in this study.All properties are taken at a base temperature, T 0 , of 293.15 K.

Table 6 .
Case studies used in this work.

Table 7 .
First two resonant frequencies.

Table 8 .
Resonant frequencies of the first mode for the encased beams simulations.
The demonstration and definition of the differential equation weak form and the subsequent element expressions are presented here.Starting with the x direction, multiplying the first elasticity equation by the test function and integrating it into the domain (Ω) results in: B xy and B xz are the boundary terms.These terms are studied in more detail in Section 2.4.
üy dVIntegrating by parts the second derivative term in both equations results in the following: ). ρM üy + λL T xy + µL xy u x + µL xx + (λ + 2µ)L yy + µL zz u y + λL yz + µL T yz u z +γN y θ = Mf y + B yx + B yy + B yz ,+ λL T yz + µL yz u y + µL xx + µL yy + (λ + 2µ)L zz u z +γN z θ = Mf z + B zx + B zy + B zz , (A11b)where B yx , B yy , B yz , B zx , B zy , B zz are the boundary terms.The new element matrices are given by:The last equation that needs to be discretized is the heat equation.Multiplying the test function and integrating the equations, one obtains:Integrating by parts and applying the discretizations in Equations (A3) and (A10), the discretized heat equation is given by:ρC p M θ + k(L xx + L yy + L zz )θ + γT 0 (N x ux + N y uy + N z uz ) = Mq + B θx + B θy + B θz , (A14)where B θx , B θy , B θz are the boundary terms.