1 Introduction

As polymer materials and fiber composite materials are widely used in aviation and other fields, these materials are formed in various structures which exhibit thermoviscoelastic behaviors in complex environments, such as thermal environment^[1]. Due to the complexity of thermoviscoelastic constitutive equations, there are a few difficulties in analytically solving the related mathematical model of thermoviscoelastic structures. Therefore, some methods have been used to obtain numerical solutions directly, such as the finite element method (FEM)^[2], the finite difference method (FDM)^[3], the Ritz method^[4], and the Laplace transform method^[5].

Up to now, the studies about static and dynamic responses of thermoelastic and viscoelastic beams^[6] have made remarkable progress. In these studies, the linear and nonlinear governing equations of thermoelastic beams were mostly based on the Euler-Bernoulli/Timoshenko beam theory. Also, by adopting different numerical methods^[7-9], the effects of thermal load parameters and material parameters on the static and dynamic responses of thermoelastic beams^[10] were analyzed. As for the static and dynamic responses of viscoelastic beams, the material constitutive relations are differential or integral types^{[11, 3]}. Chen and Ding^[12] presented the steady-state responses of axially moving viscoelastic Timoshenko beams with a Kelvin model.

Compared with the mechanical problems of thermoelastic beams, the studies on thermoviscoelastic beams started later, and current studies focused on the linear and nonlinear constitutive equations of thermoviscoelastic materials and their dynamic responses at fixed temperatures^[13-14]. As for dynamic responses of the thermoviscoelastic beam, Bernardi and Copetti^[15] established the nonlinear vibration governing equations by using a variational formulation. A differential quadrature method (DQM) was used to calculate the dynamic responses of viscoelastic fiber-metal-laminated beams under the thermal shock by Fu and Tao^[16]. However, the authors only analyzed the results in a limiting state where the relaxation modulus is constant when time approaches infinity. In addition, compared with the FEM and the FDM, the DQM is very efficient, accurate and fast for solving the governing equations with very few numbers of discrete points. It has been widely used to analyze the static and dynamic responses of structures, linear and nonlinear vibration^[17-19]. However, it is seldom reported that the DQM is used to analyze the static and dynamic responses of a thermoviscoelastic Timoshenko beam. What is more, the influence of thermal loads on the thermoviscoelastic beam needs to be further studied.

In the present work, a mathematical model in terms of generalized displacement components for a thermoviscoelastic Timoshenko beam subject to thermal loads is presented based on the assumption of small geometric deformation and linear integral type constitutive relation. Then, the governing equations and the boundary conditions are discretized with the differential method in the temporal domain and the extended DQM in the spatial domain, which could clearly identify the inner grid and border grid points in the simulations. Finally, numerical examples are carried out to verify the efficiency and accuracy of the extended DQM, and the effects of material and thermal load parameters on quasi-static and dynamic responses of a simply-supported thermoviscoelastic Timoshenko beam are analyzed.

2 Problem formulation

Consider a Timoshenko beam with an equal cross-sectional area, the height H, the length L, and the density ρ, as shown in Fig. 1. There is no heat source inside the beam subject to the thermal field θ(x, y, z, t). Here, we only consider the mechanical load (i.e., the uniform distribution load q or the concentrated load P) or the thermal load parallel to the xz-plane.

According to the Timoshenko beam theory, the displacement components at any point of the beam can be represented as^[20]

(1)

where u(x, t), v(x, t), and w(x, t) are the displacements in the x-, y-, and z-directions, respectively, while ψ(x, t) and ϕ(x, t) denote the rotations around the z- and y-directions, respectively.

Under the assumption of small deformation, the strain components at any point of the beam can be written as

(2)

Consider an isotropic and homogenous linear thermoviscoelastic material. The three- dimensional integral type relations are given as^[4]

(3)

where S_ij and e_ij are the deviatoric stress tenors and the deviatoric strain tenors (here and hereinafter), G₁ (t) and G₂ (t) are the relaxation functions describing thermoviscoelastic material properties, and θ and φ are the infinitesimal temperature deviation from the base temperature T₀ and the thermal strain functions, respectively. The symbol ⊗ is the linear Boltzmann operator^[4] and can be written as follows:

(4)

where the symbol * is the convolution product, and the symbol · represents the derivative with respect to time t.

When Poisson's ratio μ is a constant and σ_yy = σ_zz ≈ 0, by applying the Laplace transform and its inverse transform to Eq. (3), it is not difficult to obtain

(5)

where G₃ (t) is the structural function, and G₃(t)=μG₁ (t), φ_ε is the thermal function, and which were introduced to comprehensively characterize the structural and material properties of a thermoviscoelastic thin plate^[4]. Substituting Eq. (2) into Eq. (5) yields

(6)

By using Eq. (6), the shear force and bending moment are obtained as

(7)

where κ and I_y are the shear coefficient and the moment of inertia around the y-direction, respectively.

By taking a regular micro-hexahedron, the equilibrium equations can be written as

(8)

For a simply-supported beam, at x =0 and x =L, the boundary conditions are derived as

(9)

3 Solution procedure

In the following section, the solution procedures with the extended DQM in the spatial domain and the differential method in the temporal domain^[4] are used to study the static and dynamic responses of the simply-supported thermoviscoelastic Timoshenko beam.

3.1 Extended DQM in the spatial domain

Note that the traditional DQM cannot be used to solve the problem of multiple spatial variables because of the singular matrix created by the boundary conditions. However, the extended DQM^[21-22] could be used to avoid the ill-conditioning matrix by dividing the points in the weight coefficient matrix and function value matrix into two parts, i.e., the border grid points and the inner grid points. The extended DQM is based on the traditional DQM and the DQM suggested by Wang and Bert (DQWB), and can be applied to the multiple spatial variable issues and complex higher-order matrix boundary conditions^[22]. The first-order and higher-order derivatives of the function F at the given discrete point x=x_i for the independent variable x can be expressed as the following matrix forms:

(10)

where F=[F_b, F_d] ^T=[F₁, F_N, F₂, F₃, …, F_N-2, F_N-1]^T is a desired N line vector of the function values F_k, F_b =[F₁, F_N]^T is a desired 2 line vector of function values for the border grid points, and F_d =[F₂, F₃, …, F_N-2, F_N-1]^T is a desired N-2 line vector of function values for the inner grid points. Similarly, A_b, B_b, and C_b are 2 × N coefficient matrices for the border grid points, A_d, B_d, and C_d are N × (N-2) coefficient matrices for the inner grid points, and N is the grid number.

By using the extended DQM to separate the spatial variable of equilibrium equations, Eq. (8) can be written as

(11)

where Q_d = [Q_2d, Q^3d, …, Q_(N-2)d, Q_(N-1)d]^T is a desired N-2 line vector of the thermal load, is a desired N-2 line vector of the rotation ϕ, and w_d =[w_2d, w_3d, …, w_(N-2)d, w_(N-1)d]^T is a desired N-2 line vector of the displacement w.

Similarly, by using the extended DQM to separate the boundary conditions (i.e., Eq. (9)), the integral equations can be expressed as

(12)

In order to facilitate the programming and the expression, Eqs. (11) and (12) can be written in the following compact matrix forms:

(13)

(14)

where function matrix of the material properties, function matrix for border grid points, is a desired 2N line vector of the generalized displacement, is a desire 2(N-2) line vector of the generalized displacement for inner grid points, and is a desired 4 line vector of the generalized displacement, respectively. are similar, and C is a constant coefficient matrix.

3.2 Quasi-static analysis

When the inertia effect could be ignored, from Eqs. (13) and (14), the equilibrium equations of the quasi-static response of the beam can be written as

(15)

where , and P_j =[P_b, P_d]^T. For a standard linear solid material, the relaxation functions should satisfy the following conditions^[4]:

(16)

where g₀, g₁, ϕ_ε0, and ϕ_ε1 are material constants, and α and β are the reciprocals of mechanical relaxation time and thermal relaxation time, respectively. Substituting Eq. (16) into Eq. (15) and differentiating the result with respect to time t yield

(17)

It might be worth mentioning that the accuracy of numerical calculations depends on the number of discrete points and the discrete form. In this paper, the well-accepted Gauss-Chebyshev-Lobatto-type distribution of the grid point is adopted as follows^[23]:

(18)

where x_i is the spacing grid in the x -direction.

3.3 Dynamical analysis

In order to avoid the emergence of the singular matrix, the generalized displacement δ_d for the inner grid point is taken as a basic variable vector to solve the equilibrium equations. Substituting Eq. (16) into Eq. (13) and differentiating the result with respect to time t yield

(19)

Similarly, the boundary conditions of Eq. (12) can be transformed into

(20)

where the matrices M_db and M_dd are obtained from the matrix M_d, and M_bb and M_bd are obtained from the matrix M_d. Substituting Eq. (20) into Eq. (19) and eliminating the parameters of and δ_b yield

(21)

in which M₁ = -C^-1(M_dd - M_dbM_bdM_bb^-1), and +αP_b(t)).

The thermoviscoelastic beam is assumed to be in a natural state initially (t < 0 s). When t ≥ 0 s, the initial conditions are given as follows:

(22)

where are known. By substituting Eq. (22) into Eqs. (13) and (14), the initial conditions can be rewritten as

(23)

in which M₂ = C^-1M_dbM_bb^-1.

4 Results and discussion

In this section, the quasi-static and dynamical responses of a simply-supported Timoshenko beam will be discussed. First, the present results for elastic and viscoelastic beams are compared with the relevant studies^[24-27] to verify the accuracy and the convergence of the extended DQM^[21-22]. Then, the quasi-static and dynamical responses of a thermoviscoelastic beam are investigated.

4.1 Verification of the present discrete method

To verify the efficiency of the extended DQM, we discuss the deflection responses of elastic (as shown in Fig. 2) and viscoelastic (as shown in Fig. 3) Timoshenko beams subject to transverse uniform step mechanical loads, and the dynamic responses of Timoshenko beams subject to a transverse concentrated cyclic mechanical load by the FEM and the extended DQM (as shown in Fig. 4).

Figure 2 shows the mean squared error (MSE) of the maximum deflection of the elastic Timoshenko beam subject to a transverse uniform step mechanical load, which is obtained by comparing analytical solutions based on an alternative Timoshenko beam model^[24], the FEM and the extended DQM with the existing exact solution in Ref. [25]. The convergence of the FEM is obtained when the grid number for the FEM is taken as 100, which is also adopted in Fig. 4. In calculations as shown in Fig. 2, the parameters of the elastic metal beam are given as follows^[16]: L=1 m, H=0.01 m, b=0.1 m, g₀ =125 GPa, μ =0.3, and q=1 N/m. From Fig. 2, it can be observed that the numerical result obtained by the extended DQM is more accurate. In particular, with the increase in the number of non-uniformly distributed points, the numerical solution obtained by the extended DQM converges rapidly.

Figure 3 shows the MSE of the maximum deflection of a viscoelastic Timoshenko beam subject to a transverse homogenous step mechanical load. The beam is made of the southern yellow pine (SYP) or Douglas fir (DF). The MSE is obtained by comparing the solution obtained by the extended DQM with the exact solution in Ref. [26]. From Fig. 3, it can be seen that the MSE varies with different material properties. In addition, with the increase in the number of non-uniformly discrete points, the numerical solution converges rapidly and achieves a satisfactory accuracy. In the above calculations, the parameters of the viscoelastic timber beams in Fig. 3 are given as follows^[27]: L=1 m, H=0.2 m, b=0.2 m, q=5 000h(t) N/m, and h(t)=1, t>0; for an SYP beam, g₀ = 155.4 GPa, g₁ = 38.1 GPa, μ = 0.331, while for a DF beam, g₀ =183.9 GPa, g₁ = 31.64 GPa, μ = 0.347. Figure 4 shows the dynamic responses of a Timoshenko beam subject to a transverse concentrated cyclic mechanical load with the FEM^[26] and the extended DQM. From Fig. 4, it can be found that the numerical solutions solved with both methods converge rapidly with the increase in time. The parameters of Timoshenko beam and the initial conditions are given as follows^[26]: L=10 m, H=0.5 m, b=2 m, E₁ =9.8 × 10⁷ Pa, μ = 0.3, P=50 sin (πt) N (at the centre of the viscoelastic beam), w⁰ =0, and

From Figs. 2-4, it can be seen that the accuracy of the extended DQM is satisfactory. Compared with the FEM results, the DQM results can achieve the satisfactory accuracy with fewer discrete points. Particularly, for the long-time dynamic analysis of the structures, it can effectively reduce the cost of calculation. However, it is necessary to further study the applicability of the DQM in irregular regions and multidimensional spaces.

4.2 Thermoviscoelastic responses of a Timoshenko beam

Based on the above reliability analyses of the present discrete method, in the following section, we will examine the effect of material parameters and thermal load parameters on the quasi-static and dynamic responses of a thermoviscoelastic Timoshenko beam. Here, we pay more attention to the damping effect induced by the mechanical internal friction of materials. For simplicity, by taking the strategy used in Refs. [4] and [28], the thermal load is directly given as θ = (kz+b_θ)h(t), in which b_θ = 20 K. It should be mentioned that this linear temperature field conforms to the one-dimensional heat conduction equation without an inner heat source. In addition, the simply-supported boundary condition is adopted. The related parameters in Figs. 5-9 are taken as L = 1 m, b = 0.2 m, g₀ = 100 GPa, g₁ = 50 GPa, μ = 0.3, and ϕ_ε1 = 0.1 GPa.

The time history curves of the mid-span deflection of the thermoviscoelastic Timoshenko beam subject to a step thermal load (i.e., h(t)=1, t>0) with different material parameters α (i.e., the reciprocal of mechanical relaxation time) and β (i.e., the reciprocal of thermal relaxation time) are shown in Fig. 5. With the increase in the material parameters α and β, the mechanical parameter α has a more significant effect on the time history. The reason is that, for many kinds of materials such as concrete, fiber, metal, and metal composites in the aviation and construction engineering fields, the values of thermal function parameters ϕ_ε0 and ϕ_ε1 are far less than the mechanical parameters g₀ and g₁ (approximately 10^-6 to 10^-3 times)^[29-30].

As shown in Fig. 6, the influence of thermal load graded parameter k, length-to-thickness ratio L/H, and the thermal function parameter ϕ_ε0 on the steady-state deflection of the quasi-static response for the thermoviscoelastic beam subject to a step thermal load (i.e., h(t)=1, t>0) is obvious. With the increase in k, L/H, and ϕ_ε0, the steady-state deflections of the thermoviscoelastic beam increase gradually. Note that the length-to-thickness ratio L/H can change not only bending moments and shearing forces, but also the temperature field distribution of the beam^[16].

The time history curves of the mid-span deflection of the thermoviscoelastic Timoshenko beam subject to a cyclic thermal load (i.e., h(t) = cos (ωt), ω=30 Hz, t>0) with different material parameters α and β are shown in Fig. 7. The initial conditions are given as follows: It is evident that the mechanical relaxation parameter α has a more significant effect on the time history when the system tends to a steady-state response. As the mechanical relaxation parameter α =0.5 s^-1, the dynamic response of the thermoviscoelastic beam approaches the steady state when time passes 60 s. However, as the mechanical relaxation parameter α =1.0 s^-1, the dynamic response approaches the steady state earlier when time passes 30 s. However, when the thermal relaxation parameter β is changing, the time history changes little.

Figure 8 shows the effect of thermal function parameters on the time history curves of the mid-span steady-state deflection of the thermoviscoelastic Timoshenko beam subject to a cyclic thermal load (i.e., h(t)= cos (ωt), ω=30 Hz, t>0). It can be seen that the vibration amplitude enhances with the increasing thermal function parameter ϕ_ε0. Because the thermal function parameter ϕ_ε0 directly affects the value of the thermal load.

The amplitude-frequency curves of the simply-supported thermoviscoelastic Timoshenko beam subject to a cyclic thermal load (i.e., h(t)= cos (ωt), t>0) are shown in Fig. 9. The Timoshenko beam resonates at the external excitation frequency ω = 16.45 Hz. What is more, the thermal load graded parameter k and the length-to-thickness ratio L/H have great effects on the dynamic response of the beam. With the increase in k and L/H, the amplitudes increase gradually.

It should be pointed out that the heat conduction is neglected in the present analysis. As we know, this is applicable for macro-structures with isotropic and homogenous material properties^{[4, 28]}. However, the expedient should be revalued for the functionally graded materials^[31] or microscale structures^[32]. The interplay between the mechanical internal friction of the material and the thermoelastic damping induced by the coupling of the heat conduction and the mechanical motion in a microscale structure is an interesting topic, which deserves to be studied in the future.

5 Conclusions

Based on the Timoshenko beam theory and the integral-type constitutive relation, the mathematical model for the quasi-static and dynamic responses of a thermoviscoelastic Timoshenko beam is formulated. By combining the differential method in the temporal domain and the extended DQM in the spatial domain, a simple discrete scheme is applied to the integro-partial-differential initial-boundary-value problem, and its efficiency is validated by the contrast analysis with the existing analytical and numerical solutions.

Numerical results show that the geometrical parameters, material parameters, and thermal load parameters have great effects on the quasi-static and dynamic responses of a thermoviscoelastic Timoshenko beam. With the increase in the length-to-thickness ratio L/H and the thermal load graded parameter k, the quasi-static responses and steady-state dynamic response amplitudes of the thermoviscoelastic Timoshenko beam enhance. In addition, for the quasi-static and dynamic responses of a thermoviscoelastic Timoshenko beam subject to different types of thermal loads (i.e., step and cyclic thermal loads), the thermal function parameter ϕ_ε0, introduced in this paper to comprehensively characterize the structural and material properties of a thermoviscoelastic thin plate, has a great effect on the deflection. However, compared with the thermal relaxation time, the initial vibrational responses of the beam are more sensitive to the mechanical relaxation time of the thermoviscoelastic material. These conclusions can be used as reference for beam design under complex thermal environments in the aviation and construction engineering fields.