Fundamental solution and its validation by numerical inverse Laplace transformation and FEM for a damped Timoshenko beam subjected to impact and moving loads

Abstract

This paper, taking the clamped boundary condition as an example, develops Su and Ma's fundamental solutions of the dynamic responses of a Timoshenko beam subjected to impact load. Based on that, a further extension regarding the general moving load case is also established. Kelvin–Voigt damping, whether proportionally or nonproportionally damped, is incorporated into the model, making it more comprehensive than the model of Su and Ma. Numerical inverse Laplace transformation is introduced to obtain the time-domain solution, where Durbin's formula and the corresponding convergence criteria are utilized in numerical experiments. Further, the real modal superposition method is applied at an analytical level to validate the numerical results by applying a proportionally damped condition. Total comparisons are made between the methods by sufficient case studies. The dynamic responses with and without damping effect are computed with wider slenderness to verify the correctness and effectiveness of the numerical results. Furthermore, parametric studies regarding the damping coefficients are performed to explore the nonproportional damping effect. The results show that the structural damping has significant influences on the dynamic behaviors and is especially stronger at small slender ratios. As the damping decreases the inherent frequencies and excites the low-frequency modal components more actively, a resonant phenomenon appears in high slenderness case when the beam experiences a low-speed moving load. Additionally, the computations in the moving load case indicate that the algorithm convergence is preferable when the number of grids exceeds 1000.

Keywords

Timoshenko beam inverse Laplace transformation nonproportional damping finite element method mode superposition approach

1. Introduction

The dynamic analysis of a beam structure is a fundamental task in engineering practice. In previous studies, Euler–Bernoulli beam theory has been widely used due to its simplified forms, obtained by neglecting the influences of transverse shear and cross-sectional inertia (Shafei and Shafei, 2016; Adair and Jaeger, 2017). Even now, owing to the satisfactory precision in computing thin beam models, Euler–Bernoulli theory is still in favor for use in multi-body dynamical models configured in slender ratio (e.g., Yu et al., 2012; Zhang et al., 2016). In contrast, for thick beams or when aiming to obtain a higher frequency spectrum, Timoshenko's beam theory has been demonstrated to be more accurate and applicable compared with experimental data (Han et al., 1999; Stephen, 2006; Labuschagne et al., 2009; Vakil et al., 2013). In engineering applications, dynamic analysis of a beam subjected to transient impact loads is an essential task for academic studies, in which higher modes may be involved during the computations, e.g., when thousands of mode terms are superposed (Xiaochen, 2008; Su and Ma, 2012) to obtain the inner responses under a step load. In such circumstances, Timoshenko's beam theory is more competitive and superior for engineers to conduct theoretical research, e.g., the damped beam with typical surface loads used in this paper.

Several mathematical and numerical approaches have been proposed and developed over the past half-century to compute the transient responses of the beam subjected to basic impact loads or moving loads (Ding et al., 2016), e.g., large-scale bridge structures. Among them, the complex mode superposition method was employed in Greco and Santini (2002) for dynamical studies of a simply supported beam with dampers at two ends and a surface moving load, but only Euler–Bernoulli functions were adopted. Azam et al. (2013), also using a modal superposition method, transformed the governing equations of a Timoshenko beam into a set of ordinary differential equations in state-space manners and solved them via a numerical technique. In Yufeng et al. (2002), analytical solutions of the Timoshenko beam subjected to an impulsive particle with four types of boundary conditions were developed on the basis of a similar mode superposition method. Achenbach and Sun (1965) investigated the implications of inherent parameters, e.g., load velocity, by using a complex Fourier transform, and similar derivations were conducted by Ortner and Wagner (1996) to address an initial boundary issue regarding a simply supported semi-infinite Timoshenko beam. Further, by using Lagrange equations, Kocatrk and Simsek (2006) obtained the dynamic response of an eccentrically prestressed viscoelastic Timoshenko beam incorporating a moving harmonic load, in which the trial functions that approximate the lateral deflections and the rotations of the cross-section are constructed in polynomial forms. In Huntley and Zinober (1981), double Laplace transformation was performed to evaluate the transient responses and the frequency responses of a semi-infinite Timoshenko beam; this work was considered an original study based on previous achievements. Furthermore, the wave splitting treatment was used by Johansson (2004) to derive the Green's operator in the Laplace domain and complicated time asymptotic expansions and contour integrations were computed further. In recent studies (Su and Ma, 2011, 2012), theoretical and numerical approaches have been utilized to obtain the transient waves of a simply supported Timoshenko beam, in which a typical ray solution, normal mode theory, and Laplace inversion are all included and the numerical inverse Laplace transformations were proved to be of superior flexibilities. The literature review presented above shows some remarkable achievements in similar or analogical academic domains, and other prominent developments are also available, as presented in detail in Su and Ma's (2012, 2011) work.

Further, the damping effect must be considered as an essential factor for practical application, wherein several approximate damping models have been proposed and applied to structural modelling, e.g., the Rayleigh damping model and the Kelvin–Voigt (K-V) damping model. Further, according to its inherent mathematical properties, damping can be generally classified into proportional situation and nonproportional situation; e.g., they are involved in the following studies. Hu et al. (2010) derived the control equation through the Lagrangian approach in consideration of the damping effect, in which the approximation strategy for solution was employed via shape functions and the transient responses of the cantilevered beam were obtained. To study the thermal effects on the overall dynamic response, Gu et al. (2015) utilized complex mode theory and accounted for the proportional Rayleigh damping model ultimately to perform possible numerical experiments. The finite element method (FEM) was used by Chen (2014) to analyze the impacts of the twist angle on inherent frequency, where the K-V damping effect was considered, whereas the explicit modal shapes and corresponding dynamic responses were not worked out. Liu et al. (2016) considered a nonclassically damped linear system with symmetric governing matrices, and because of the modal nonorthogonality, complex mode theory was introduced for equation decoupling. Sorrentino et al. (2003) also conducted an eigenvalue analysis of a nonproportional damped beam, and an expansion in state form was generated in conjunction with a transfer matrix technique. Such analyses were also performed in (Qu and Selvam, 2002; Li et al., 2016; Lzaro, 2016), wherein the governing matrices were uniformly symmetric, although the effects of nonproportional damping were included. Through this specific literature review, the following conclusions are drawn at this stage. The physical beam systems, particularly for general continuous systems, are mostly modelled in proportional damping conditions or via symmetrical damping matrices uniformly, whereas these simplified situations do not always hold in actual applications. Thus, developing a more flexible and applicable approach that can account for arbitrary damping conditions is necessary and of great significance. Because the response studies with nonproportionally damped condition are quite scarce, the investigations of this paper are also imperative.

Hence, the main objective of this paper is to develop the fundamental solutions to more extensive applications by means of Laplace transformations, in which the K-V damping model is considered and a nonproportional damping matrix is generated. Because only the undamped case with cantilevered boundary condition and concentrated force is computed in Su and Ma's research, it is regarded as insufficient and can be improved by the use of more physical factors and external conditions. Thus, this paper is also a direct extension of the previous work (Su and Ma, 2012). Based on the improved fundamental solutions with consideration of the damping effect in the concentrated loaded case, a special FEM configured algorithm is exploited for the moving load case. Specifically, the entire paper is organized as follows. In Section 2, the nondimensional governing equations are introduced. Two methods are introduced in Section 3 with entirely different deriving processes: one is characterized by Laplace transformation (denoted as the L-method), and the other is the normal mode method (denoted as the N-method). In fact, the N-method is introduced to validate the correctness of the L-method, in which a proportionally damped condition is performed for classical mode utilizations. Based on the results in Section 3, Section 4 presents a proposed FEM algorithm for the moving load case, and original derivations are incorporated in detail. Case studies are presented in Section 5, wherein Durbin's inverse Laplace transformation and the corresponding error criteria are shown to derive the time-domain solutions of the L-method. In the verifications, the results of undamped and proportionally damped cases are compared first. Then, the influences of nonproportional damping are evaluated through parametric studies. Further, the effectiveness of the finite element algorithm and its using conditions in the moving load case is measured via a series of numerical experiments. Eventually, the most relevant conclusions are drawn in Section 6.

2. Modeling of the governing equations

Figure 1 illustrates the physical elastic beam model bounded at ends and the principle of Timoshenko model is also presented. As shown, the original oxyz frame is established at end1 whose origin is situated at the neutral point and y follows the right-hand screw rule, L is the longitudinal span of the beam, the vertical step force is imposed at d along the x-axis, and H denotes the lateral height of beam cross-section. z_b and z_s, which together constitute a total lateral displacement z^T in the z direction, express the deformations from bending effect and shear effect respectively. According to Timoshenko's beam theory (Wu and Chen, 2015; Zhang et al., 2018b), the physical differential equations about z_b and z_s hold:

kGA \frac{\partial^{2} z_{s}}{\partial x^{2}} = ρ A \frac{\partial^{2} (z_{s} + z_{b})}{\partial t^{2}}

(1a)

EI \frac{\partial^{3} z_{b}}{\partial x^{3}} + kGA \frac{\partial z_{s}}{\partial x} = ρ I \frac{\partial^{3} z_{b}}{\partial x \partial t^{2}}

(1b)

in which E and G represent Young's modulus and the shear modulus of the beam, I is the moment inertia of the cross-section along the y-axis, ρ is the volume density of the beam structure, k is the so-called shear coefficient (Jensen, 1966), giving the ratio of averaged shear strain to the strain at the centroid of cross-section, A denotes the cross-sectional area of the beam, ϕ and γ are the structural angles due to bending and shearing from the derivatives of z_b and z_s mathematically, and t is the general dimensional time.

Further, an essential hypothesis that the lateral displacement z^T is identical within the entire cross-section should be mentioned to validate the further damping model Zhang et al., 2017). Hence, by applying the specific K-V damping model (Zhao et al., 2005; Chen, 2014), the general damping effects due to the structural deformations can be evaluated by

σ = E ɛ + C_{b} \frac{\partial ɛ}{\partial t}, τ = G γ + C_{s} \frac{\partial γ}{\partial t}

(2)

where σ and ε denote the stretching stress and the strain flow within the beam, and τ describes the shear stress. C_b and C_s are the notations of viscous damping coefficients related to inner stretching and shearing motions respectively. A linear relationship is obtained for the cross-sectional rotation from Timoshenko theory, where h is introduced as an integral variable for H along the lateral span:

ɛ = \frac{\partial ϕ}{\partial x} h

(3)

Thus, the resultant bending moment M and the shear force V of the cross-section are available by following integrations (B is span of the cross-section in y-direction):

M = \int_{- H / 2}^{H / 2} σ Bhdh = EI \frac{\partial ϕ}{\partial x} + C_{b} I \frac{\partial^{2} ϕ}{\partial x \partial t}

(4a)

V = \int_{- H / 2}^{H / 2} τ Bdh = kGA γ + C_{s} kA \frac{\partial γ}{\partial t}

(4b)

Since ϕ and γ are dependent on z_b and z_s, $z_{b}$ and $z_{s}$ are chosen as the independent variables. Thus, the governing equations of the system can be derived through Newtonian superpositions in equation (1) and equation (4). Herein, $q = [\begin{matrix} z_{b} & z_{s} \end{matrix}]^{T}$ is introduced for simplification, whereby we obtain

\begin{matrix} [\begin{matrix} ρ A & ρ A \\ ρ I \frac{\partial}{\partial x} & 0 \end{matrix}] q^{. .} + [\begin{matrix} 0 & - C_{s} kA \frac{\partial^{2}}{\partial x^{2}} \\ - C_{b} I \frac{\partial^{3}}{\partial x^{3}} & - C_{s} kA \frac{\partial}{\partial x} \end{matrix}] \\ \overset{\cdot}{q} + [\begin{matrix} 0 & - kGA \frac{\partial^{2}}{\partial x^{2}} \\ - EI \frac{\partial^{3}}{\partial x^{3}} & - kGA \frac{\partial}{\partial x} \end{matrix}] q = F \end{matrix}

(5)

in which an external load

F

in vector manner is incorporated which consists of forces and moments for a general situation. For the specific impact load shown in Figure 1,

F

holds for

F = [\begin{matrix} f (x, t) & 0 \end{matrix}]^{T}

exactly;

{()}^{.}

used in equation (5) denotes the time derivative. Nondimensionalizations about system variables are performed for computational brevity by introducing

\begin{matrix} X = \frac{x}{L}, Z_{b} = \frac{z_{b}}{L}, Z_{s} = \frac{z_{s}}{L}, Z^{T} = \frac{z^{T}}{L}, \\ \bar{q} = \frac{q}{L} = [\begin{matrix} Z_{b} & Z_{s} \end{matrix}]^{T}, λ_{0} = \sqrt{\frac{kG}{E}}, L_{r} = \frac{L}{r} \end{matrix}

(r is the equivalent inertial radius of the profile gyration, defined by

I = A r^{2}

w_{0} = \sqrt{\frac{\sqrt{kGE}}{ρ L^{2}}}, C_{s}^{*} = \frac{C_{s} k}{w_{0} ρ L^{2}}, C_{b}^{*} = \frac{C_{b}}{w_{0} ρ L^{2}}, ζ = w_{0} t

, and thus the inherent parameters are refined and equation (5) is rewritten as

\begin{matrix} [\begin{matrix} 1 & 1 \\ \frac{\partial}{\partial x} & 0 \end{matrix}] \overset{··}{\bar{q}} + [\begin{matrix} 0 & - C_{s}^{*} \frac{\partial^{2}}{\partial x^{2}} \\ - C_{b}^{*} \frac{\partial^{3}}{\partial x^{3}} & - C_{s}^{*} {L_{r}}^{2} \frac{\partial}{\partial x} \end{matrix}] \\ \overset{\cdot}{\bar{q}} + [\begin{matrix} 0 & - λ_{0} \frac{\partial^{2}}{\partial x^{2}} \\ - \frac{1}{λ_{0}} \frac{\partial^{3}}{\partial x^{3}} & - λ_{0} {L_{r}}^{2} \frac{\partial}{\partial x} \end{matrix}] \bar{q} = F^{*} \end{matrix}

(6)

3. Fundamental solutions under impact load

3.1. L-method

First, the original derivations of the L-method are shown in this part by considering the extra damping effect.

3.1.1. Laplace transformation

Laplace transformation is applied to equation (6) over ζ via a substitutive parameter p initially. By supposing ${\bar{q}}_{L} = {[\begin{array}{l} Z_{b} (p) & Z_{s} (p) \end{array}]}_{L}^{T} \exp^{λ X}$ (herein, subscript $()_{L}$ , or sequential superscript $()^{L}$ , means the Laplace domain quantity. λ is a characteristic value of a fundamental spatial shape function of the beam), equation (6) is transformed to

\begin{matrix} [\begin{matrix} p^{2} & p^{2} - C_{s}^{*} p λ^{2} - λ_{0} λ^{2} \\ - p^{2} + C_{b}^{*} p λ^{2} + \frac{λ^{2}}{λ_{0}} & C_{s}^{*} {L_{r}}^{2} p + λ_{0} {L_{r}}^{2} \end{matrix}] \bar{q_{L}} = F_{L}^{*} \end{matrix}

(7)

In regard to the homogeneous form of equation (7), the presence of nontrivial solutions must be satisfied, which means

\det \begin{matrix} [\begin{matrix} p^{2} & p^{2} - C_{s}^{*} p λ^{2} - λ_{0} λ^{2} \\ - p^{2} + C_{b}^{*} p λ^{2} + \frac{λ^{2}}{λ_{0}} & C_{s}^{*} {L_{r}}^{2} p + λ_{0} {L_{r}}^{2} \end{matrix}] = 0 \end{matrix}

(8)

The mathematical solutions of this character equation hold:

λ^{2} = \frac{- b \pm \sqrt{b^{2} - 4 ac}}{2 a}

(9)

accompanied with three intermediate coefficients, a, b, and c:

\begin{matrix} a = (C_{b}^{*} p + \frac{1}{λ_{0}}) (C_{s}^{*} p + λ_{0}), \\ b = - [(C_{s}^{*} p + λ_{0}) + (C_{b}^{*} p + \frac{1}{λ_{0}})] p^{2}, \\ c = p^{2} (C_{s}^{*} {L_{r}}^{2} p + λ_{0} {L_{r}}^{2} + p^{2}) . \end{matrix}

Hence, the general spatial solutions of equation (7) can be explicitly denoted as ( $λ_{1} = (\frac{- b + \sqrt{b^{2} - 4 ac}}{2 a})^{1 / 2}, λ_{2} = (\frac{- b - \sqrt{b^{2} - 4 ac}}{2 a})^{1 / 2}$ ):

\begin{matrix} Z_{b}^{L} (p, X) = u_{1 +} (p) e^{- λ_{1} X} + u_{1 -} (p) e^{λ_{1} X} \\ + u_{2 +} (p) e^{- λ_{2} X} + u_{2 -} (p) e^{λ_{2} X} \end{matrix}

(10a)

\begin{matrix} Z_{s}^{L} (p, X) = v_{1 +} (p) e^{- λ_{1} X} + v_{1 -} (p) e^{λ_{1} X} \\ + v_{2 +} (p) e^{- λ_{2} X} + v_{2 -} (p) e^{λ_{2} X} \end{matrix}

(10b)

After substituting equation (10) into equation (7), the relationships among coefficients $u_{i \pm}$ and $v_{i \pm}$ are demonstrated satisfying

K_{i} = \frac{v_{i +}}{u_{i +}} = \frac{v_{i -}}{u_{i -}} = \frac{p^{2} - C_{b}^{*} p {λ_{i}}^{2} - \frac{{λ_{i}}^{2}}{λ_{0}}}{C_{s}^{*} {L_{r}}^{2} p + λ_{0} {L_{r}}^{2}} (i = 1 or 2)

(11)

For further derivations,

\begin{array}{l} Z_{L}^{T} = Z_{b}^{L} + Z_{s}^{L}, φ_{L} = Z_{b}^{L}'^{,} V_{L} = \frac{k G A}{w_{0}^{2} r ρ A L} \\ Z_{s}^{L}'^{=} λ_{0} L_{r} Z_{s}^{L}'^{,} M_{L} = \frac{E I}{w_{0}^{2} r ρ I} Z_{b}^{L^{''}} = \frac{L_{r}}{λ_{0}} Z_{b}^{L^{''}} \end{array}

denoting the total displacements, bending angle, shear force, and bending moment in the Laplace domain, are introduced, in agreement with previous nondimensionalizations. Specifically, these inner physical quantities are available at present because of the general solutions in equation (10). We write them in a compact form as follows:

b = Rc

(12)

where

b = (Z_{L}^{T}, ϕ_{L}, V_{L}, M_{L})^{T}

c = (u_{1 +}, u_{2 +}, u_{1 -}, u_{2 -})^{T}

R (p, X) = [R_{ij}] |_{4 \times 4}

and the phase-related spatial matrix

R (p, X)

is also represented by

\begin{matrix} R_{11} = (1 + K_{1}) e^{- λ_{1} X}, R_{12} = (1 + K_{2}) e^{- λ_{2} X}, \\ R_{13} = (1 + K_{1}) e^{λ_{1} X}, R_{14} = (1 + K_{2}) e^{λ_{2} X}; \\ R_{21} = - λ_{1} e^{- λ_{1} X}, R_{22} = - λ_{2} e^{- λ_{2} X}, \\ R_{23} = λ_{1} e^{λ_{1} X}, R_{24} = λ_{2} e^{λ_{2} X}; \\ R_{31} = - λ_{0} L_{r} K_{1} λ_{1} e^{- λ_{1} X}, R_{32} = - λ_{0} L_{r} K_{2} λ_{2} e^{- λ_{2} X}, \\ R_{33} = λ_{0} L_{r} K_{1} λ_{1} e^{λ_{1} X}, R_{34} = λ_{0} L_{r} K_{2} λ_{2} e^{λ_{2} X}; \\ R_{41} = \frac{L_{r}}{λ_{0}} {λ_{1}}^{2} e^{- λ_{1} X}, R_{42} = \frac{L_{r}}{λ_{0}} {λ_{2}}^{2} e^{- λ_{2} X}, \\ R_{43} = \frac{L_{r}}{λ_{0}} {λ_{1}}^{2} e^{λ_{1} X}, R_{44} = \frac{L_{r}}{λ_{0}} {λ_{2}}^{2} e^{λ_{2} X} . \end{matrix}

Su and Ma (2011) demonstrated that the unknown c was determined by the boundary conditions at end1 and end2. For horizontal comparison, we apply clamped conditions to both end1 and end2 as an example that differs from the situations of Su and Ma. The solving principles for derivations in other forms of boundaries are also analogous, and the boundary conditions at the two ends in immovable states physically hold:

Z_{L}^{T} |_{X = 0} = Z_{L}^{T} |_{X = 1} = 0, ϕ_{L} |_{X = 0} = ϕ_{L} |_{X = 1} = 0

(13)

In accordance with equation (12), we transform equation (13) to a matrix format $R_{b} c = t$ , i.e.,

\begin{matrix} [\begin{matrix} \begin{matrix} R_{11} (0) & R_{12} (0) \\ R_{21} (0) & R_{22} (0) \end{matrix} & \begin{matrix} R_{13} (0) & R_{14} (0) \\ R_{23} (0) & R_{24} (0) \end{matrix} \\ \begin{matrix} R_{11} (1) & R_{12} (1) \\ R_{21} (1) & R_{22} (1) \end{matrix} & \begin{matrix} R_{13} (1) & R_{14} (1) \\ R_{23} (1) & R_{24} (1) \end{matrix} \end{matrix}] c = [\begin{matrix} \begin{matrix} Z_{L}^{T} (0) \\ ϕ_{L} (0) \end{matrix} \\ \begin{matrix} Z_{L}^{T} (1) \\ ϕ_{L} (1) \end{matrix} \end{matrix}] \end{matrix}

(14)

According to Su and Ma's (2012) theory, c can be derived by

c = {R_{b}}^{- 1} t = (I - r)^{- 1} s

(15)

where a group of matrix transformations shown in sequential formulas are accompanied as well. This series of treatments are also available in Lee and Ma's (2000) original studies:

\begin{matrix} & r = - D^{- 1} (L + U) \\ R_{b} = D + L + U = D (I - r) \\ s = D^{- 1} t \end{matrix}

where

\begin{matrix} D = [\begin{matrix} \begin{matrix} R_{11} (0) & R_{12} (0) \\ R_{21} (0) & R_{22} (0) \end{matrix} & \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} \\ \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} & \begin{matrix} R_{13} (1) & R_{14} (1) \\ R_{23} (1) & R_{24} (1) \end{matrix} \end{matrix}] \\ L = [\begin{matrix} \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} & \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} \\ \begin{matrix} R_{11} (1) & R_{12} (1) \\ R_{21} (1) & R_{22} (1) \end{matrix} & \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} \end{matrix}] \\ U = [\begin{matrix} \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} & \begin{matrix} R_{13} (0) & R_{14} (0) \\ R_{23} (0) & R_{24} (0) \end{matrix} \\ \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} & \begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix} \end{matrix}] \end{matrix}

According to equation (15), the objective is transformed to solving the unknown s, which was defined as a source function of boundary loading in wave propagating analyses.

3.1.2. Response under impact load

Considering the detailed loading situation shown in Figure 2, we can see that an interior impact load $F_{0} H (t)$ located at $X = d^{*}$ is applied to the beam where the Heaviside function H(t) is introduced for the impact load. The Laplace-domain responses of the system excited by this load, according to Su and Ma (2011)'s wave propagating model, can be expressed as

b = R (I - r)^{- 1} s^{*} - {Rs}_{0}^{*}

(16)

in which b is composed of a general solution

s^{*}

and a particular solution

s_{0}^{*}

, specifically

s_{0}^{*} (X, p) = (s_{+} (d) 0)^{T}, if 0 < X < d^{*};

(17a)

s_{0}^{*} (X, p) = (0 s_{-} (d))^{T}, if d^{*} < X < 0;

(17b)

s^{*} (X, p) = (s_{+} (d) s_{-} (d))^{T}, for all X

(17c)

where

s_{+}

and

s_{-}

, as the most fundamental quantities, are still undetermined. Therefore, a basic semi-infinite beam model, as depicted in Figure 2, is constructed to perform the sequential derivations. Note that, according to the equivalence principle, the step force is equally decomposed into two

(F_{0} H (t)) / 2

parts, with one for the left segment and the other for the right segment. Based on a previous derivation process (Su and Ma, 2012), we develop and obtain the dimensionless solutions of this semi-infinite beam model that considers the damping effect. The exact solutions are

\tilde{Z_{b}^{L}} = \tilde{u_{1 +}} (p) e^{- λ_{1} \tilde{X}} + \tilde{u_{2 +}} (p) e^{- λ_{2} \tilde{X}}

(18a)

\tilde{Z_{s}^{L}} = K_{1} \tilde{u_{1 +}} (p) e^{- λ_{1} \tilde{X}} + K_{2} \tilde{u_{2 +}} (p) e^{- λ_{2} \tilde{X}}

(18b)

Here, $\tilde{()}$ denotes the intermediate quantities of the semi-infinite beam model. Symmetric boundary conditions at $X = d^{*}$ (equivalent to $\tilde{X} = 0$ ) are considered to obtain the final solutions, which includes a nonrotation condition and a nonshear condition of the cross-section as a result of the symmetry of the two segments. After substituting the Laplace-domain forms of these physical requirements into equation (18), we derive

- K_{1} \tilde{u_{1 +}} λ_{1} - K_{2} \tilde{u_{2 +}} λ_{2} = \frac{{F_{0}}^{*}}{2 p}

(19a)

- \tilde{u_{1 +}} λ_{1} - \tilde{u_{2 +}} λ_{2} = 0

(19b)

in which

F_{0}

is nondimensionalized through

{F_{0}}^{*} = \frac{F_{0}}{{w_{0}}^{2} ρ AL}

. Thus, the explicit solutions of the semi-infinite beam model are

\tilde{Z_{b}^{L}} = \frac{{F_{0}}^{*} (e^{- λ_{1} \tilde{X}} + e^{λ_{1} \tilde{X}})}{2 p λ_{1} (K_{2} - K_{1})} - \frac{{F_{0}}^{*} (e^{- λ_{2} \tilde{X}} + e^{λ_{2} \tilde{X}})}{2 p λ_{2} (K_{2} - K_{1})}

(20a)

\tilde{Z_{s}^{L}} = \frac{K_{1} {F_{0}}^{*} (e^{- λ_{1} \tilde{X}} + e^{λ_{1} \tilde{X}})}{2 p λ_{1} (K_{2} - K_{1})} - \frac{K_{2} {F_{0}}^{*} (e^{- λ_{2} \tilde{X}} + e^{λ_{2} \tilde{X}})}{2 p λ_{2} (K_{2} - K_{1})}

(20b)

which are very crucial in the entire solving process, and after translating

\tilde{X}

to the original X coordinate, equation (20) is rewritten as

\begin{matrix} Z_{b}^{L} = \frac{{F_{0}}^{*} (e^{- λ_{1} (X - d^{*})} + e^{λ_{1} (X - d^{*})})}{2 p λ_{1} (K_{2} - K_{1})} \\ - \frac{{F_{0}}^{*} (e^{- λ_{2} (X - d^{*})} + e^{λ_{2} (X - d^{*})})}{2 p λ_{2} (K_{2} - K_{1})} \end{matrix}

(21a)

\begin{matrix} Z_{s}^{L} = \frac{K_{1} {F_{0}}^{*} (e^{- λ_{1} (X - d^{*})} + e^{λ_{1} (X - d^{*})})}{2 p λ_{1} (K_{2} - K_{1})} \\ - \frac{K_{2} {F_{0}}^{*} (e^{- λ_{2} (X - d^{*})} + e^{λ_{2} (X - d^{*})})}{2 p λ_{2} (K_{2} - K_{1})} \end{matrix}

(21b)

Further, according to the relationships of equation (21) and equation (17), the explicit $s_{+}$ and $s_{-}$ are derived ultimately as

s_{+} (d) = (u_{1 +}, u_{2 +})

(22a)

s_{-} (d) = (u_{1 -}, u_{2 -})

(22b)

where

\begin{matrix} & u_{1 +} = - \frac{{F_{0}}^{*} e^{λ_{1} d^{*}}}{2 p λ_{1} (K_{1} - K_{2})}, u_{1 -} = - \frac{{F_{0}}^{*} e^{- λ_{1} d^{*}}}{2 p λ_{1} (K_{1} - K_{2})}; \\ u_{2 +} = \frac{{F_{0}}^{*} e^{λ_{2} d^{*}}}{2 p λ_{2} (K_{1} - K_{2})}, u_{2 -} = \frac{{F_{0}}^{*} e^{- λ_{2} d^{*}}}{2 p λ_{2} (K_{1} - K_{2})} . \end{matrix}

Based on the basic solutions in equation (22), the exact $R (p, X)$ of equation (12) is also obtainable. It is also encouraging that current results possess analytical levels as there is no approximation within the derivations. Because extra damping effects, different boundary conditions and additional nondimensionalizations are considered, the situations involved here are more complicated compared with those of Su and Ma (2012). Overall, the current derivations have no conflict with that of Su and Ma.

3.2. N-method

Solutions by means of the normal mode approach are also generated in this part for comparison. Because the K-V damping model is incorporated, the difficulties in solving the modal responses are increased. Generally, the traditional mode superposition method can be divided into real mode theory and complex mode theory. Although real mode theory is more concise and easier to manipulate, its utilization should ensure that, e.g., the system matrices possess proportional damping condition $K M^{- 1} C = C M^{- 1} K$ (Caughey and O'Kelly, 1965; Adhikari, 2006) (here, $K$ denotes the stiffness matrix of the entire system, $M$ denotes the corresponding inertial matrix, and $C$ describes the damping matrix); otherwise, the complex-domain derivation cannot be avoided. As the N-method is used to verify the correctness of the L-method, it is preferable to transform the derivations into the real domain. Because C_b and C_s are still undetermined, a perfect proportional damping condition is constructed; this treatment also avoids the lengthy derivations in complex mode theory.

3.2.1. Configuration of proportional damping

In accordance with equation (7), λ is reserved as the general spatial solution. We consider the homogeneous form of equation (6) at first, thus

\begin{matrix} [\begin{matrix} 1 & 1 \\ 0 & λ \end{matrix}] [\begin{matrix} {Z_{s}}^{\cdot \cdot} \\ {Z_{b}}^{\cdot \cdot} \end{matrix}] + [\begin{matrix} - C_{s}^{*} λ^{2} & 0 \\ - C_{s}^{*} {L_{r}}^{2} λ & - C_{b}^{*} λ^{3} \end{matrix}] [\begin{matrix} {Z_{s}}^{\cdot} \\ {Z_{b}}^{\cdot} \end{matrix}] \\ + [\begin{matrix} - λ_{0} λ^{2} & 0 \\ - λ_{0} {L_{r}}^{2} λ & - \frac{λ^{3}}{λ_{0}} \end{matrix}] [\begin{matrix} Z_{s} \\ Z_{b} \end{matrix}] = 0 \end{matrix}

(23)

To construct a proportionally damped condition, equation (23) is compacted to $M [\overset{··}{Z_{s}} \overset{··}{Z_{b}}]^{T} + C [\overset{\cdot}{Z_{s}} \overset{\cdot}{Z_{b}}]^{T} + K [Z_{s} Z_{b}]^{T} = 0$ , and we derive (Caughey and O'Kelly, 1961, 1965)

K M^{- 1} C = [\begin{matrix} C_{s}^{*} λ_{0} λ^{2} (λ^{2} - {L_{r}}^{2}) & - C_{b}^{*} λ_{0} λ^{4} \\ C_{s}^{*} {L_{r}}^{2} λ (λ_{0} λ^{2} - λ_{0} {L_{r}}^{2} + \frac{λ^{2}}{λ_{0}}) & C_{b}^{*} λ^{3} (- λ_{0} {L_{r}}^{2} + \frac{λ^{2}}{λ_{0}}) \end{matrix}]

(24a)

C M^{- 1} K = [\begin{matrix} λ_{0} C_{s}^{*} λ^{2} (λ^{2} - {L_{r}}^{2}) & - C_{s}^{*} \frac{λ^{4}}{λ_{0}} \\ C_{s}^{*} {L_{r}}^{2} λ (λ_{0} λ^{2} - λ_{0} {L_{r}}^{2} + \frac{C_{b}^{*}}{C_{s}^{*}} λ_{0} λ^{2}) & C_{b}^{*} λ^{3} (- \frac{C_{s}^{*}}{C_{b}^{*}} \frac{1}{λ_{0}} {L_{r}}^{2} + \frac{λ^{2}}{λ_{0}}) \end{matrix}]

(24b)

In fact, the calculated values of C_b and C_s in the K-V damping model originate mainly from experimental tests or empirical data; e.g., in Zhang et al. (2017), approximate damping coefficients were employed to conduct the numerical computations. Zhao et al. (2005) also analyzed the stability of the beam system via mathematic discussions by taking the damping coefficients as arbitrary nonnegative parameters. Herein, note that if $\frac{C_{s}^{*}}{C_{b}^{*}} = {λ_{0}}^{2}$ is set up, then the inner damping will satisfy a proportional damping standard $K M^{- 1} C = C M^{- 1} K$ . This original relational expression was proposed by Caughey and O'Kelly (1965) to guarantee that the system possesses a complete set of real orthonormal eigenvectors, i.e., the classical normal modes. Adhikari (2006) developed the necessary conditions about the classical normal modes for the systems with singular mass matrices, in which a viscously damped linear system with nonnegative definite system matrices was considered.

Aiming to validate the L-method, such particular treatment is primarily mathematical reasonable, so the complex domain derivations is avoided successfully. In fact, Chen (2014) set C_b and C_s to identical values in actual computations and analyses. Capsoni et al. (2013) also discussed the critical conditions of C_b and C_s via decoupling the shearing and bending schemes, in which the relations of C_b and C_s were demonstrated to be interdependent. Gu et al. (2015) derived a solution in the structural aspect by decoupling the governing equations and considering the K-V damping model at the same time, and a proportional Rayleigh damping model, as shown in their equation (58), was introduced for each mode in the eventual response analysis. Caughey and O'Kelly (1963) also highlighted that the usual treatment of linearly damped parameter systems assumed that the system equations may be transformed into a symmetrical set of equations. This assumption is justified in passive systems, which is beneficial to the inherent modal analyses presented in this part. Hence, the proportional damping condition constructed here has no conflict with the physical mechanism, and thus the real-domain derivations are performed in the following sections. In addition, the parametric studies regarding the two damping coefficients are to be conducted in case studies to explore their influence laws. Once the system owns the proportional damping matrix, the derivations of dynamic responses under loading state can be greatly shortened. Further, the classical normal modes conditions are also met for derivations (Caughey and O'Kelly, 1965). The original notations $C_{s}^{*}$ and $C_{b}^{*}$ are still effective in the sequential derivations, and the inherent relationship $\frac{C_{s}^{*}}{C_{b}^{*}} = {λ_{0}}^{2}$ holds between them in the N-method.

The solving process of the N-method is resolved into two parts: the undamped case and the damped case.

3.2.2. Undamped case

According to the principle of variable separation, the general nondamping solution of equation (23) is $[\binom{Z_{s}}{\begin{matrix} Z_{b} \end{matrix}}] = [\binom{C_{1}}{\begin{matrix} C_{2} \end{matrix}}] e^{iw ζ}$ . And on account of the presence of nontrivial $[\binom{C_{1}}{\begin{matrix} C_{2} \end{matrix}}]$ , the characteristic equation is generated as ( $Ω = w^{2}$ )

λ^{4} + Ω λ_{0} + \frac{1}{λ_{0}} λ^{2} + Ω^{2} - λ_{0} {L_{r}}^{2} Ω = 0

(25)

whose roots are mathematically expressed as

λ^{2} = - \frac{1}{2} Ω (λ_{0} + \frac{1}{λ_{0}}) [1 \pm \sqrt{\frac{(λ_{0} - \frac{1}{λ_{0}})^{2} + 4 \frac{λ_{0} {L_{r}}^{2}}{Ω}}{(λ_{0} + \frac{1}{λ_{0}})^{2}}}]

(26)

Equation (26) constitutes the fundamental modal solutions in the undamped case, which can be generally expressed through

Z_{i} = a_{1} e^{λ_{1} X} + a_{2} e^{- λ_{1} X} + a_{3} e^{λ_{2} X} + a_{4} e^{- λ_{2} X}

(27)

In accordance with the L-method, only clamped boundaries are concerned and the solving processes of other types of constraints, e.g., free end, simply-supported end, etc., are also analogous. In agreement with equation (13), the clamped constraints in time domain hold:

\begin{matrix} Z^{T} (0, ζ) |_{X = 0} = Z^{T} (0, ζ) |_{X = 1} = 0, \\ ϕ (0, ζ) |_{X = 0} = ϕ (0, ζ) |_{X = 1} = 0 \end{matrix}

(28)

If $Ω < λ_{0} {L_{r}}^{2}$ , the solution configurations of equation (26) consist of two pairs of roots in real aspect and imaginary aspect respectively, which are denoted as $\pm η$ and $\pm α$ . Correspondingly, the frequency equation and the characteristic vectors are

| \begin{matrix} - \frac{λ_{0} η^{2}}{Ω} \sinh (η) - \frac{K_{11}}{K_{12}} \frac{λ_{0} α^{2}}{Ω} sin (α) & \frac{λ_{0} η^{2}}{Ω} (cos (α) - \cosh (η)) \\ K_{11} (η \cosh (η) - α cos (α)) & K_{11} η \sinh (η) - K_{12} \frac{η^{2}}{α^{2}} α sin (α) \end{matrix} | = 0

(29)

and

\begin{matrix} [\begin{matrix} Z_{sn} \\ Z_{bn} \end{matrix}] = [\begin{matrix} \sinh (η X) \\ K_{11} \sinh (η X) \end{matrix}] + \frac{A_{2}}{A_{1}} [\begin{matrix} \cosh (η X) \\ K_{11} \cosh (η X) \end{matrix}] \\ + \frac{A_{3}}{A_{1}} [\begin{matrix} sin (α X) \\ K_{12} sin (α X) \end{matrix}] + \frac{A_{4}}{A_{2}} \frac{A_{2}}{A_{1}} [\begin{matrix} cos (α X) \\ K_{12} cos (α X) \end{matrix}] \end{matrix}

(30)

in which

\begin{matrix} & K_{11} = - \frac{λ_{0} η^{2} + Ω}{Ω}, K_{12} = \frac{λ_{0} α^{2} - Ω}{Ω}, \frac{A_{4}}{A_{2}} = \frac{η^{2}}{α^{2}}, \\ \frac{A_{3}}{A_{1}} = - \frac{K_{11}}{K_{12}}, \frac{A_{2}}{A_{1}} = - \frac{\sinh (η) + \frac{K_{11}}{K_{12}} \frac{α^{2}}{η^{2}} sin (α)}{\cosh (η) - cos (α)} \end{matrix}

If $Ω = λ_{0} {L_{r}}^{2}$ , the solution configurations of equation (26) are composed of two multiple roots 0 and a pair of imaginary roots $\pm α$ . The corresponding frequency equation and the characteristic vectors are

| \begin{matrix} sin (α) & cos (α) - 1 \\ K_{21} cos (α) - K_{21} & - K_{21} sin (α) - α \end{matrix} | = 0

(31)

and

\begin{matrix} [\begin{matrix} Z_{sn} \\ Z_{bn} \end{matrix}] = \frac{A_{1}}{A_{3}} \frac{A_{3}}{A_{4}} [\begin{matrix} X \\ - X \end{matrix}] + \frac{A_{2}}{A_{4}} [\begin{matrix} 1 - \frac{Ω}{2 λ_{0}} X^{2} \\ \frac{Ω}{2 λ_{0}} X^{2} \end{matrix}] \\ + \frac{A_{3}}{A_{4}} [\begin{matrix} sin (α X) \\ K_{21} sin (α X) \end{matrix}] + [\begin{matrix} cos (α X) \\ K_{21} cos (α X) \end{matrix}] \end{matrix}

(32)

in which

\begin{matrix} K_{21} = \frac{λ_{0} α^{2} - Ω}{Ω}, \frac{A_{2}}{A_{4}} = - \frac{λ_{0} α^{2}}{Ω}, \\ \frac{A_{1}}{A_{3}} = α K_{21}, \frac{A_{3}}{A_{4}} = \frac{1 - cos (α)}{sin (α)} \end{matrix}

If $Ω > λ_{0} {L_{r}}^{2}$ , the solution configurations of equation (26) are composed of two pairs of imaginary roots $\pm α$ and $\pm β$ , whose frequency equation and characteristic vectors hold:

| \begin{matrix} sin (β) - \frac{α}{β} \frac{K_{31}}{K_{32}} sin (α) & cos (β) - cos (α) \\ β K_{31} (cos (β) - cos (α)) & - β K_{31} cos (β) + α K_{32} sin (α) \frac{β^{2}}{α^{2}} \end{matrix} | = 0

(33)

and

\begin{matrix} [\begin{matrix} Z_{sn} \\ Z_{bn} \end{matrix}] = [\begin{matrix} sin (β X) \\ K_{31} sin (β X) \end{matrix}] + \frac{A_{2}}{A_{1}} [\begin{matrix} cos (β X) \\ K_{31} cos (β X) \end{matrix}] \\ + \frac{A_{3}}{A_{1}} [\begin{matrix} sin (α X) \\ K_{32} sin (α X) \end{matrix}] + \frac{A_{4}}{A_{2}} \frac{A_{2}}{A_{1}} [\begin{matrix} cos (α X) \\ K_{32} cos (α X) \end{matrix}] \end{matrix}

(34)

in which

\begin{matrix} K_{31} = \frac{λ_{0} β^{2} - Ω}{Ω}, K_{32} = \frac{λ_{0} α^{2} - Ω}{Ω}, \frac{A_{4}}{A_{2}} = - \frac{β^{2}}{α^{2}}, \\ \frac{A_{3}}{A_{1}} = - \frac{β}{α} \frac{K_{31}}{K_{32}}, \frac{A_{2}}{A_{1}} = \frac{sin (β) - \frac{α}{β} \frac{K_{31}}{K_{32}} sin (α)}{cos (α) - cos (β)} \end{matrix}

Equations (29)–(34) show the complete modal solutions of the undamped system, which are also consistent (Yufeng et al., 2002; Su and Ma, 2011; Zhang et al., 2018a) with same or similar boundary conditions. Herein, all three cases are of both mathematical and applicable meanings in regard to the mode analyses.

3.2.3. Damped case

Similar to the undamped case, we introduce the superposition-form solutions for damped condition at first (subscript $()_{d}$ denotes the damped case here):

[\begin{matrix} Z_{s} \\ Z_{b} \end{matrix}]_{d} = \sum_{n = 1}^{\propto} [\begin{matrix} Z_{sn} (X) \\ Z_{bn} (X) \end{matrix}]_{d} Q_{n} (ζ)

(35)

where

Q_{n} (ζ)

represents the modal response of mode n (positive integers n and sequential m denote the modal numbers). Then, substituting equation (35) into equation (23) and in consideration of the external loads shown in equation (6), premultiplying

[\begin{matrix} Z_{sm} & Z_{bm} \end{matrix}]_{d}

, and finally integrating the consequence through operator

\int_{0}^{1} d X

and utilizing the classical modal orthogonality in the meanwhile (Zhang et al., 2017), we are able to transform equation (23) to a general second-order ordinary equation as follows:

{Q_{m}}^{\cdot \cdot} + 2 κ_{m} w_{m} {Q_{m}}^{\cdot} + Ω_{m} Q_{m} = F_{m}^{*}

(36)

in which

Ω_{m} = {w_{m}}^{2} = \frac{- \int_{0}^{1} (Z_{sm}^{2} + Z_{bm} Z_{sm} + Z_{bm} Z_{bm}') d X}{\int_{0}^{1} (λ_{0} Z_{sm} Z_{sm}'' + λ_{0} {L_{r}}^{2} Z_{bm} Z_{sm}' + \frac{1}{λ_{0}} Z_{bm} Z_{bm}^{(3)}) d X}

κ_{m} = \frac{- \int_{0}^{1} (C_{s}^{*} Z_{sm} Z_{sm}'' + C_{s}^{*} {L_{r}}^{2} Z_{bm} Z_{sm}' + C_{b}^{*} Z_{bm} Z_{bm}^{(3)}) d X}{2 w_{m} \int_{0}^{1} (λ_{0} Z_{sm} Z_{sm}'' + λ_{0} {L_{r}}^{2} Z_{bm} Z_{sm}' + \frac{1}{λ_{0}} Z_{bm} Z_{bm}^{(3)}) d X}

As presented in Zhang et al.'s (2017, 2018a) work, the explicit solution of equation (36) can be obtained according to the discussions about damping classifications. Herein, we only compute the dynamic responses from an initial equilibrium circumstance of the system.

If $κ_{m} = 0$ ,

{Q_{m} (ζ) = \frac{1}{w_{m}} \int}_{0}^{ζ} F_{m}^{*} (τ) sin [w_{n}^{i} (ζ - τ)] d τ

(37)

If $0 < κ_{m} < 1$ ,

{Q_{m} (ζ) = \frac{1}{w_{m}^{d}} \int}_{0}^{ζ} F_{m}^{*} (τ) e^{- κ_{m} w_{m} (ζ - τ)} sin [w_{m}^{d} (ζ - τ)] d τ

(38)

If $κ_{m} = 1$ ,

{Q_{m} (ζ) = \int}_{0}^{ζ} F_{m}^{*} (τ) (ζ - τ) e^{- w_{m} (ζ - τ)} d τ

(39)

If $κ_{m} > 1$ ,

Q_{m} (ζ) = \frac{\int_{0}^{ζ} F_{m}^{*} (τ) [\begin{matrix} e^{(- κ_{m} + \sqrt{{κ_{m}}^{2} - 1}) w_{m} (ζ - τ)} \\ + e^{(- κ_{m} - \sqrt{{κ_{m}}^{2} - 1}) w_{m} (ζ - τ)} \end{matrix}] d τ}{2 w_{m} \sqrt{{κ_{m}}^{2} - 1}}

(40)

in which the explicit modal impact load at

d^{*}

F_{m}^{*} = - {F_{0}}^{*} H (ζ) Z_{sm} (X = d^{*})

(41)

4. Expansion to moving load case

Based on the fundamental solutions in impact load case, it is possible and meaningful to expand the L-method to wider applications in engineering, namely the moving load case conducted in this part. When a moving load acts on the beam, the spatial variable X is added to $F_{L}^{*}$ and complicate the analytical derivations of L-method. Under such circumstances, a numerical strategy is put forward here originally to solve the moving load issue. On account of the achievements in Section 3.1.2, we exploit a FEM based technique for the moving load case.

As illustrated in Figure 3, the entire beam is divided into N successive grids longitudinally via N + 1 nodes. Notations (j) and (j) are introduced for each segment and each node respectively. Actually, the intervals within the N + 1 nodes can be configured uniformly or unequally, and one can also refine the divisions in the vicinity of the most concerned regions. Correspondingly, an equivalent treatment is performed for the moving load by discretizing it into a series of immovable concentrated forces. As shown in Figure 4, $F_{0} H_{j} (ζ)$ of segment (j) is enabled when moving load moves within interval $X_{[j - 1]} < X \leq X_{[j]}$ , i.e.,

H_{j} (ζ) = H δ (ζ - ζ_{j - 1}) - H δ (ζ - ζ_{j}), at X = d_{j}^{*}

(42)

where

ζ_{[j - 1]} = \frac{w_{0} x_{j - 1}}{v}

and

ζ_{[j]} = \frac{w_{0} x_{j}}{v}

are the nondimensional coordinates of any two contiguous nodes. It is noted that the applied position of each force unit is set to the midpoint of segment uniformly, namely

d_{j}^{*} = \frac{X_{[j]} + X_{[j - 1]}}{2}

. Specifically, the moving load is transformed to the superposition form of

F_{0} H_{j} (ζ)

f (x, t) = \sum_{j = 1}^{N} F_{0} H_{j} (t)

(43)

According to the principle of the Laplace method, transformation $L [f (t - t_{0})] = e^{- s t_{0}} F (s)$ holds on the basis of $L [f (t)] = F (s)$ . Thus, the responses excited by segment (j), consistent with equation (21), can be expressed as

Z_{b}^{L} j_{=} (e^{- p ζ_{[j - 1]}} - e^{- p ζ_{[j]}}) [\frac{F_{0}^{*} (e^{- λ_{1} (X - d_{j}^{*})} + e^{λ_{1} (X - d_{j}^{*})})}{2 p λ_{1} (K_{2} - K_{1})} - \frac{F_{0}^{*} (e^{- λ_{2} (X - d_{j}^{*})} + e^{λ_{2} (X - d_{j}^{*})})}{2 p λ_{2} (K_{2} - K_{1})}]

(44a)

\begin{matrix} Z_{s}^{L} j_{=} (e^{- p ζ_{[j - 1]}} - e^{- p ζ_{[j]}}) \\ [\frac{K_{1} F_{0}^{*} (e^{- λ_{1} (X - d_{j}^{*})} + e^{λ_{1} (X - d_{j}^{*})})}{2 p λ_{1} (K_{2} - K_{1})} - \frac{K_{2} F_{0}^{*} (e^{- λ_{2} (X - d_{j}^{*})} + e^{λ_{2} (X - d_{j}^{*})})}{2 p λ_{2} (K_{2} - K_{1})}] \end{matrix}

(44b)

Herein, explicit $u_{1 +}^{j}, u_{1 -}^{j}, u_{2 +}^{j}, u_{2 -}^{j}$ of segment (j) in the moving load case are also given.

u_{1 +}^{j} = - (e^{- p ζ_{[j - 1]}} - e^{- p ζ_{[j]}}) \frac{{F_{0}}^{*} e^{λ_{1} d_{j}^{*}}}{2 p λ_{1} (K_{1} - K_{2})}

(45a)

u_{1 -}^{j} = - (e^{- p ζ_{[j - 1]}} - e^{- p ζ_{[j]}}) \frac{{F_{0}}^{*} e^{- λ_{1} d_{j}^{*}}}{2 p λ_{1} (K_{1} - K_{2})}

(45b)

u_{2 +}^{j} = (e^{- p ζ_{[j - 1]}} - e^{- p ζ_{[j]}}) \frac{{F_{0}}^{*} e^{λ_{2} d_{j}^{*}}}{2 p λ_{2} (K_{1} - K_{2})}

(45c)

u_{2 -}^{j} = (e^{- p ζ_{[j - 1]}} - e^{- p ζ_{[j]}}) \frac{{F_{0}}^{*} e^{- λ_{2} d_{j}^{*}}}{2 p λ_{2} (K_{1} - K_{2})}

(45d)

After substituting equations (44) and (45) into equation (16) and further superposing the series together, we derive the eventual numerical solution. In matrix format, it has the form:

b = R (I - r)^{- 1} \sum_{j = 1}^{N} s_{j}^{*} - R \sum_{j = 1}^{N} s_{j 0}^{*}

(46)

which is in numerical level with N-nodes approximations.

5. Case studies

Numerical experiments are conducted in this part to verify the correctness of the derivations, and thorough comparisons are made between the L-method and the N-method. Because the damping effect is considered a key factor, contrastive analyses regarding the undamped and damped systems are also conducted. As a result, parametric impacts on the dynamical responses are measured for damping coefficients, taking advantage of the L-method. Further, computations are made for the moving load case, where the conditions of use of the FEM algorithm are evaluated through convergence experiments regarding the division number N. Moreover, the dynamic performances of the system in undamped and damped conditions are analyzed.

On account of the complexities of equations (16) and (46), the numerical Laplace inversion method, which was proved fairly effective by Durbin (1974) and Su and Ma (2012), is applied to obtain the final time-domain solutions. According to Durbin's formula, the inversion of a general Laplace-domain function $f (x, ζ)$ can be approximated numerically through

\begin{matrix} f (x, ζ) \approx \\ \frac{e^{ζ_{0} ζ}}{T} [\begin{matrix} - \frac{1}{2} Re {F (x, ζ_{0})} + \sum_{k = 1}^{N_{0}} Re {F (x, ζ_{0} + \frac{ik π}{T})} cos (\frac{k π ζ}{T}) - \sum_{k = 0}^{N_{0}} Im {F (x, ζ_{0} + \frac{ik π}{T})} sin (\frac{k π ζ}{T})] \end{matrix} \end{matrix}

(47)

and an error criterion for equation (47), which contributes to the parameter determinations, is also included, i.e.,

| Re {F (x, ζ + i N_{0} \frac{π}{T})} | \leq \frac{ε T}{e^{ζ T}}, | Im {F (x, ζ + i N_{0} \frac{π}{T})} | \leq \frac{ε T}{e^{ζ T}}

(48)

in which parameter T is generally dependent on the total response time ζ considered in the computations. ζ₀ and N₀, which are introduced as a Laplace coefficient and a truncation parameter, are determined by equation (48) with an ε precision.

5.1. Verifications

Verifications are conducted first to validate the effectiveness of the L-method. Because slenderness is an important factor in beam modelling (Chang and Lee, 2009), two typical slender ratios L_r = 10 and L_r = 100 are incorporated as the representatives of a short beam and a long beam. Initial data based on a real structure are derived from Zhang et al. (2017). In detail,

\begin{matrix} & L = 1800 m, A = 20 m^{2}, I = 309.1563 m^{4}, \\ E = 3.63 \times 10^{10} N m^{- 2}, G = 1.5783 \times 10^{10} N m^{- 2}, \\ ρ = 5.1282 \times 10^{3} kg m^{- 3}, C_{b} = 2 \times 10^{8} pas, \\ k = 5 / 6, and F_{0} = 2 \times 10^{4} N \end{matrix}

which are nondimensionalized through the definitions of Section 2. Further,

ε = 10^{- 6}

is applied to the parameter determinations of equation (47).

Figure 5 shows the essential modal characteristics of the system computed through N-method for deeper analyses, including the inherent frequencies in Figure 5(a) and the corresponding damping rates in Figure 5(b). We can observe that the frequencies in the case of L_r = 100 are much smaller than that of L_r = 10 case within the first 40 modes and that they approach in the vicinity of 50th mode with comparable magnitudes. Further, note that the damping ratios are apparently larger in case of L_r = 10. In other words, higher modes in a small slender ratio will reach the over-damping vibrational state earlier under external excitations. It is also proved that the critical damping for both cases appears near the condition of $Ω = λ_{0} {L_{r}}^{2}$ . Specifically, the critical damping lies at the third mode for the case of L_r = 10 and at the 49th mode for the case of L_r = 100.

Figures 6 –13 illustrate the dynamic responses of the two methods with identical computing conditions. Note that the loads are imposed on $d^{*} = 0.5$ and $d^{*} = 0.25$ , respectively, which are constructed to explore a symmetric situation and an asymmetric situation. The receivers are set at X = 0.2 and X = 0.8, following Su and Ma (2012). In actual calculations, the first 2000 modes are superposed for the N-method to guarantee the computational accuracy (Su and Ma, 2011; Zhang et al., 2018a). For the L-method, a MATLAB searching program, which is interpreted below, is exploited for ζ₀ and T to adapt the various computing conditions from Figures 6 –13. First, ζ_max is determined in advance to start the searching program of the computations. Herein, ζ_max is the maximum response time, and $5 * 2 pi / w_{1}$ is utilized as the actual experimental ζ_max. Next, the searching regions regarding the time are expanded to [ζ_max, $10 \times ζ_{max}$ ] to derive a pair of preferable ζ₀ and T, which is further divided into 1000 equal subintervals. To each subinterval, the errors in equation (47) are evaluated and recorded once equation (47) is satisfied. Eventually, the smallest error within the database is chosen and corresponds to the optimal ζ₀ and T.

As depicted in Figures 6 and 7 or Figures 10 and 11(a), the responses at X = 0.2 and X = 0.8 under symmetric $d^{*} = 0.5$ loading condition demonstrate the L-method is fairly stable because the plots, whether in undamped condition or damped condition, are absolutely identical. High agreements are also observed in contrastive plots of the two methods, where the accordant fluctuating tendencies of the responses in different case conditions validate the effectiveness of the numerical solutions. Furthermore, apparent discrepancies appear between the methods, especially in the case of L_r = 10, e.g., Figures 8 and 9, the separations from the plots are quite prominent. Further, regular phenomena of case L_r = 10 illustrated in Figure 8(b) to Figure 13(b) are detected in damped condition, where the peaks and troughs from the responses of L-method are either advanced or postponed compared with the standard N-method. In contrast, these phase differences in the case of L_r = 100, as proved by Figure 10(b) to Figure 13(b), almost disappear. This imperfection regarding the L-method reveals the numerical algorithm of equation (47) is not positively unexceptionable for arbitrary small to large L_r because of its high dependence on ζ₀ and T. Zhang et al. (2018a), by means of numerical experiments, concluded that Durbin's formula may endure poor arithmetic flexibility, which is the greatest limitation of its computational effectiveness. Furthermore, the L-method clearly behaves worse in Figures 8 and 9 compared with other the set of figures, as the computing conditions of Figures 8 and 9 incorporate an asymmetric load at $d^{*} = 0.25$ . This demonstrates that the precision of the L-method is related to the spatial factor, as demonstrated by the fact that equation (47) is closely dependent on x. Given that the load position is not always specific in engineering, a more advanced searching technology for ζ₀ and T are demanded and will be of great significance in future studies. In fact, these concerns and discussions are not considered by Su and Ma (2012) because they used only one pair of ζ₀ and T in their studies. Generally, although the results of L-method have not coincided with the N-method accurately, the trends of the responses, as demonstrated by the fluctuations in each plot, are basically consistent. Additionally, damping has similar impacts on the dynamical responses for both methods, but the implications are more prominent in small L_r. The displacements Z^T also converge rapidly in the low L_r condition, echoing the results in Figure 5(b) that the L_r = 10 case has much larger damping rates. In contrast, the damping in case of L_r = 100 behaves very inactively, and the responses in the undamped and damped conditions almost coincide within the given time period.

Thus, it is concluded that although the preconditions of equation (47) decrease the effectiveness of L-method, the reasonable and decent predictions of the dynamical responses and the free requirements for the damping conditions indicate that the numerical method is a good alternative and has favorable application prospects in engineering. The full analyses in Section 5.1 are also considered extensions of Su and Ma (2012).

5.2. Damping variations

Although the N-method is established at an analytical level, the solutions from normal mode theory can only handle proportional damping conditions. Note that the numerical solution that is demonstrated to be valid in Section 5.1 has no restriction on damping conditions, and parametric studies about the two damping coefficients are conducted via the L-method to evaluate their implications. Specifically, continuous variations regarding C_b and C_s, varied via proportion ratios ranging from 0 to 1 relative to the original proportionally damped conditions, are performed in the computations.

Figures 14 and 15 exhibit the results of various structural damping, in which the combination $d^{*} = 0.5$ and X = 0.2 is taken as an example. Moreover, variations regarding the bending rotation ϕ, the shearing force $\bar{V}$ and the sectional moment $\bar{M}$ are obtained by equation (12). Overall, the changes in C_b and C_s have significant implications on the dynamical performances. In particular, when C_b and C_s increase from ratio = 0 to ratio = 0.1, the variations in Z^T, $ϕ, \bar{V}$ and $\bar{M}$ are considerably large. After ratio = 0.1, the responses converge to the original proportional damping values accordingly. Note that the separated plots in intervals $[0, 0.1]$ and [0.1,1] for each quantity present common piecewise-linear trends uniformly. Additionally, L_r affects the responses to some extent, in which the curves of $C_{b} - ζ = 1000$ and $C_{s} - ζ = 1000$ , as demonstrated by Z^T of Figure 15(a), tend to coincide with each other in the case of L_r = 100. Consistent with Z^T, the plots of other quantities, i.e., ϕ, V and M, present similar tendencies. Overall, the consequences computed in different slenderness indicate that the nonproportional damping has accordant mechanical actions on the system.

5.3. Moving load case

The dynamical responses of the beam subjected to moving load are computed in this section for the algorithmic validation, and the uniform node divisions are applied to equation (46) for the numerical experiments and different node numbers $N = 100, 500, 1000, 10000$ are adopted to investigate the application conditions of the FEM. A low-frequency constant $v_{0} = 0.001 m / s$ , which provides a longer action time along the longitudinal span, is employed in equation (42) to detect the response variations (Zhang et al., 2018a). Herein, only the proportionally damped condition is considered to explore the computational effectiveness and efficiency.

Figures 16 and 17 illustrate the results of moving load cases that are also characterized by undamped cases and damped cases respectively. Note that the time course the process when the load passes on the beam surface is equivalent to the longitudinal span along the axis from spatial X = 0 to X = 1. From the response amplitudes, it is concluded that the plots gradually converge to a stable result with the increment of N. The plots in cases of N = 1000 and N = 10000 almost coincide each other within the whole span range precisely, whereas the responses in cases of N = 100 and N = 500 show great discrepancies that also increase with the action time and almost diverge at X = 1, e.g., Figures 16(a) and 17(a). Because the computational costs grow greatly with N and become a dominant limitation of the algorithm, N = 1000 is deemed to be large enough and to have superior precision for actual applications. Further, the results in the L_r = 10 case show that the damped responses in Figure 16(b) have no fluctuations compared with Figure 16(a), which is attributed to the higher damping ratios evaluated by Figure 5(b). The greater damping effects in small slenderness suppress the vibrational movements and decrease the response amplitudes. In contrast, the damped responses in Figure 17(b) in the case of L_r = 100 behave even more significantly compared with the undamped case in Figure 17(a), as a result of the following factors. Due to the damping effects, the inherent frequency in the L_r = 100 case is reduced to lower $w_{m}^{d} = w_{m} \sqrt{1 - {κ_{m}}^{2}}$ in advance (Zhang et al., 2017). When the low-speed moving load passes, the low-component external excitations tend to interact with more low-frequency modes of the beam, leading to the resonance phenomenon and the considerable enlarged response amplitude. Such a phenomenon is also demonstrated by Zhang et al.'s (2017, 2018a) damped cases with identical structural data. Overall, the numerical results generated by FEM can predict the essential vibrational mechanism reasonably; thus, the algorithmic effectiveness is validated, and the corresponding application condition is obtained preliminarily.

The L-method has strong application prospects due to its uniform expressions for different damping conditions. A discretization strategy in the moving load case is employed to approximate a moving load, and this is also in keeping with the traditional FE approaches, e.g., a moving load in ANSYS or ABAQUS platform is realized by a series of stationary loads (Dimitrovová and Rodrigues, 2011; Hou and Zhao, 2012). Up to now, the commercial FE tools still depend on secondary developments to simulate such a moving load, resulting in a deficiency of the algorithmic flexibility. Further, the massive node information in ANSYS or ABAQUS environment should be computed simultaneously to iterate the next increment, whereas the L-method only computes the information of concerned nodes that has potential to achieve lower computational costs. In addition, the L-method also has potential to handle some fluid-structure issues, e.g., the nonproportional damping and elastic effects of water in Zhang et al. (2017) can be inserted into the beam equations to obtain a similar fundamental solution. On the whole, the flexible L-method is an effective alternative to address the similar issues.

6. Conclusions

By applying the Timoshenko beam model and the K-V damping model, in this paper, the fundamental responses of a typical beam system subjected to impact load were derived, in which both analytical and numerical solutions were generated. The influences of structural damping were incorporated into the original derivations, which were considered an improvement to Su and Ma's work. Further, an extending solution regarding the moving load case was established by FEM. The Durbin formula was introduced into the L-method to obtain the ultimate time-domain solutions. In the case studies, the contrastive analyses based on the L-method and N-method indicated that the numerical solutions perform well in predicting the dynamical responses, whereas its computational precisions were demonstrated to be affected by the inherent parameters. The results also revealed the significant implications of the structural damping, e.g., the damping effects behave more prominently in small L_r situation. Parametric studies regarding the damping coefficients were conducted via the L-method, in which a uniform piecewise-linear law was concluded for each physical quantity. Furthermore, the computations in the moving load case demonstrate that the algorithmic effectiveness is proportional to the grid number, where refined division performs well. The damping is proved sensitive to the slenderness, and diverging responses can be detected by the L-method for higher slenderness. Overall, the L-method is more flexible due to the free damping requirement and thus has significant application prospects.

Footnotes

Acknowledgement

The authors gratefully acknowledge the fund support provided by Prof. Ming Zhu who is now a doctoral advisor at Beihang University.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported in part by the National Key R&D Program of China (grant number 2016YFB1200100).

References

Achenbach

Sun

(1965) Moving load on a flexibly supported Timoshenko beam. International Journal of Solids and Structures 1(4): 353–370.

Adair

Jaeger

(2017) A power series solution for rotating nonuniform Euler-Bernoulli cantilever beams. Journal of Vibration and Control. 24(17): 3855–3864.

Adhikari

(2006) Damping modelling using generalized proportional damping. Journal of Sound and Vibration 293(1–2): 156–170.

Azam

Mofid

Khoraskani

(2013) Dynamic response of Timoshenko beam under moving mass. Scientia Iranica 20(1): 50–56.

Capsoni

Vigan

Bani-Hani

(2013) On damping effects in Timoshenko beams. International Journal of Mechanical Sciences 73(8): 27–39.

Caughey

O'Kelly

MEJ

(1961) Effect of damping on the natural frequencies of linear dynamic systems. Journal of the Acoustical Society of America 33(11): 1458–1461.

Caughey TK and O'Kelly MEJ (1963) General theory of vibration of damped linear dynamic systems. Technical report, Dynamics Laboratory, California Institute of Technology, USA.

Caughey

O'Kelly

MEJ

(1965) Classical normal modes in damped linear dynamic systems. Journal of Applied Mechanics 32(3): 583–588.

Chang

Lee

(2009) Free vibration of a single-walled carbon nanotube containing a fluid flow using the Timoshenko beam model. Physics Letters A 373(10): 982–985.

10.

Chen

(2014) Parametric studies on bending vibration of axially-loaded twisted Timoshenko beams with locally distributed Kelvin–Voigt damping. International Journal of Mechanical Sciences 88: 61–70.

11.

Dimitrovová Z and Rodrigues AFS (2011) An enhanced moving window method: Applications to high-speed tracks. In: International conference on civil, structural and environmental engineering computing, Stirlingshire, UK 2011.

12.

Ding

Zhu

(2016) Effects of axial load and structural damping on wave propagation in periodic Timoshenko beams on elastic foundations under moving loads. Physics Letters A 380(31–32): 2335–2341.

13.

Durbin

(1974) Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate's method. Computer Journal 17(4): 371–376.

14.

Greco

Santini

(2002) Dynamic response of a flexural non-classically damped continuous beam under moving loadings. Computers and Structures 80(26): 1945–1953.

15.

Qin

Chu

(2015) Analytical analysis of the thermal effect on vibrations of a damped Timoshenko beam. Mechanical Systems and Signal Processing 60: 619–643.

16.

Han

Benaroya

Wei

(1999) Dynamics of transversely vibrating beams using four engineering theories. Journal of Sound and Vibration 225(5): 935–988.

17.

Hou

Zhao

(2012) Response analysis of saturated asphalt pavement under action of moving load. Transportation Standardization 21: 12–14.

18.

Wang

Zhang

(2010) Approximate analytical solutions and experimental analysis for transient response of constrained damping cantilever beam. Applied Mathematics and Mechanics 31(11): 1359–1370.

19.

Huntley

Zinober

ASI

(1981) Transient responses of the Timoshenko beam evaluated by the double Laplace transform. Computers and Structures 14(3–4): 319–323.

20.

Jensen

(1966) On the shear coefficient in Timoshenko's beam theory. Journal of Sound and Vibration 33(2): 621–635.

21.

Johansson

(2004) A purely directional wave splitting for the homogeneous Timoshenko beam. Journal of Sound and Vibration 275(1): 299–316.

22.

Kocatrk

Simsek

(2006) Dynamic analysis of eccentrically prestressed viscoelastic Timoshenko beams under a moving harmonic load. Computers and Structures 84(31): 2113–2127.

23.

Labuschagne

van Rensburg

NFJ

van der Merwe

(2009) Comparison of linear beam theories. Mathematical and Computer Modelling 49(1–2): 20–30.

24.

Lee

(2000) Transient elastic waves propagating in a multi-layered medium subjected to in-plane dynamic loadings. I. Theory. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 456(1998): 1355–1374.

25.

Wang

(2016) Accurate method for harmonic responses of non-classically damped systems in the middle frequency range. Journal of Vibration and Control 22(2): 426–441.

26.

Liu

Lian

Liang

et al. (2016) An improved complex multiple-support response spectrum method for the non-classically damped linear system with coupled damping. Bulletin of Earthquake Engineering 14(1): 161–184.

27.

Lzaro

(2016) Eigensolutions of non-proportionally damped systems based on continuous damping sensitivity. Journal of Sound and Vibration 363: 532–544.

28.

Ortner

Wagner

(1996) Solution of the initial-boundary value problem for the simply supported semi-infinite Timoshenko beam. Journal of Elasticity 42(3): 217–241.

29.

Selvam

(2002) Efficient method for dynamic condensation of nonclassically damped vibration systems. AIAA Journal 40(40): 368–375.

30.

Shafei

(2016) Oblique impact of multi-flexible-link systems. Journal of Vibration and Control 24(5): 904–923.

31.

Sorrentino

Marchesiello

Piombo

BAD

(2003) A new analytical technique for vibration analysis of non-proportionally damped beams. Journal of Sound and Vibration 265(4): 765–782.

32.

Stephen

(2006) The second spectrum of Timoshenko beam theory—further assessment. Journal of Sound and Vibration 292(1–2): 372–389.

33.

(2011) Theoretical analysis of transient waves in a simply-supported Timoshenko beam by ray and normal mode methods. International Journal of Solids and Structures 48(3–4): 535–552.

34.

(2012) Transient wave analysis of a cantilever Timoshenko beam subjected to impact loading by Laplace transform and normal mode methods. International Journal of Solids and Structures 49(9): 1158–1176.

35.

Vakil

Sharbati

Vakil

et al. (2013) Vibration analysis of a Timoshenko beam on a moving base. Journal of Vibration and Control 21(6): 1068–1085.

36.

Chen

(2015) Free vibration analysis of a Timoshenko beam carrying multiple spring–mass systems by using the numerical assembly technique. International Journal for Numerical Methods in Engineering 50(5): 1039–1058.

37.

Xiaochen W (2008) Numerical analysis for dynamic response of bridge under moving loads. PhD thesis, Zhejiang University of Technology, China.

38.

Wen

Shen

et al. (2012) Propagation of flexural wave in periodic beam on elastic foundations. Physics Letters A 376(4): 626–630.

39.

Yufeng

Yuansong

Dechao

et al. (2002) Elastic impact on finite Timoshenko beam. Acta Mechanica Sinica 18(3): 252–263.

40.

Zhang X, Zhu M and Liang H (2017) Dynamic analysis of the continuous fluid-structure system based on Timoshenko model and considering damping. In: 58th AIAA/ASCE/AHS/ASC structures, structural dynamics, and materials conference, Location?, Date? 2017.

41.

Zhang

Zhu

Zhe Wu

(2018a) Response analyses of Timoshenko beam based on K-V damping model. Journal of Beijing University of Aeronautics and Astronautics 44(3): 500–507.

42.

Zhang

Zhu

et al. (2018b) Dynamic analysis of Timoshenko beam with arbitrary constraints and a further optimization based on least energy principle. Mathematical Problems in Engineering 2018(6): 1–19.

43.

Zhang

Pang

Fan

(2016) Analyses of transverse vibrations of axially pretensioned viscoelastic nanobeams with small size and surface effects. Physics Letters A 380(29–30): 2294–2299.

44.

Zhao

Liu

Zhang

(2005) Stability for the Timoshenko beam system with local Kelvin–Voigt damping. Acta Mathematic Sinica 21(3): 655–666.