Impact response of a Timoshenko-type viscoelastic beam considering the extension of its middle surface

In the present paper, the problem of low-velocity impact of an elastic sphere against a viscoelastic Timoshenko-type beam is studied considering the extension of its middle surface. The viscoelastic features of the beam out of the contact domain are governed by the standard linear solid model with derivatives of integer order, while within the contact domain the fractional derivative standard linear solid model is utilized, in so doing rheological constants of the material in both models are the same. However the presence of the additional parameter, i.e. fractional parameter which could vary from zero to unit, allows one to vary beam’s viscosity, since the structure of the beam’s material within this zone may be damaged, resulting in the decrease of the beam material viscosity in the contact zone. Consideration for transient waves (surfaces of strong discontinuity) propagating in the target out of the contact zone via the theory of discontinuities and determination of the desired values behind the surfaces of discontinuities upto the contact domain with the help of ray series, as well as the utilization of the Hertz theory in the contact zone allow one to obtain a set of two integro-differential equations, which govern the desired values, namely: the local bearing of the target and impactor’s materials and the displacement of the beam within the contact domain.

where [v] is the discontinuity in the velocity of longitudinal displacement on the plane transient wave propagating along the beam during the impact process with the velocity G 1 = √ E/ρ, E is Young's modulus, ρ is the density, ν is Poisson's ratio, h b is the beam thickness, and α is the local bearing the target and impactor's materials. The above formula is valid for an elastic thin plate, but it is invalid for a beam.
However, the correct formula for the beam (which will be used further in the present paper) was derived 50 years ago in the classical work by Landau and Lifshitz (1965) In the same page of (p. 345 in Vershinin 2014), the equation of motion of the contact zone was presented as where N is the longitudinal force, w and Q are transverse displacement and force, respectively, P(t) is the contact force, z is the coordinate directed along the beam axis, a dot denotes the time-derivative, a is the radius of the contact zone, in so doing the author of (Vershinin 2014) has considered that plane transient waves (surfaces of discontinuity) propagate from the contact zone during the process of impact. However, if the contact domain is a circular disk with a volume hπ a 2 , then the waves travelling from a circular contact zone are diverging circles.
In order the waves propagating from the contact zone to be the plane waves, the contact domain should be a rectangular parallelepiped with the volume 2aF, where F is the cross sectional area of the beam, and the equation should have the form In the same page of (p. 345 in Vershinin 2014), the final equation was presented as K is the shear coefficient, and µ is the shear modulus. This equation is also incorrect.
The correct equation neglecting the inertia of the contact zone (this important assumption was not mentioned in Vershinin (2014) ) has the form with the correct coefficient e = Based on the aforesaid, the further numerical treatment presented in Vershinin (2014) occurs to be invalid. [ In the present paper, we present not only the correct equation for the impact response of the elastic Timoshenko beam, but the deduction and analysis of equations describing the behaviour of a viscoelastic Timoshenko-type beam impacted by an elastic sphere are given considering the damage of the target material within the contact domain.

Problem formulation
Let an elastic sphere with the radius R and mass m move along the y−axis with a constant velocity V 0 towards a Timoshenko-type viscoelastic homogeneous isotropic beam of the width h ( Fig. 1), viscoelastic features of which are described by the standard linear solid model with conventional derivatives. The dynamic behaviour of such a beam with due account for extension of its middle surface is described by the following set of equations: where M, Q, and N are the bending moment, the shear and longitudinal forces, respectively, u and w are longitudinal and transverse displacements, respectively, ψ is the angle of rotation of the cross section around the z-axis, v =u, W =ẇ, � =ψ, F and I are the cross-sectional area and the moment of inertia with respect to the x-axis, respectively, ρ is the density, K is the shear coefficient dependent on beam's geometrical dimensions and the form of its cross section, and an overdot denotes the time derivative.
In equations (2) and (4), the operator corresponding to the Young modulus has the form where Z(t) is a desired function, E ∞ and E 0 are the non-relaxed (instantaneous modulus of elasticity, or the glassy modulus) and relaxed elastic (prolonged modulus of elasticity, or the rubbery modulus) moduli which are connected with the relaxation time τ ε and retardation time τ σ by the following relationship: In equation (3), the operator corresponding to the shear modulus has the form where µ ∞ is the nonrelaxed magnitude of the shear modulus, and n and t ε are yet unknown constants. Shitikova (2013, 2015) have shown that if the operator E is assumed to be assigned and the operator of the triaxial extension-compression K, according to experimental data (Rabotnov 1966), is considered to be time-independent, i.e. K = K ∞ , where K ∞ is the non-relaxed module, then Poisson's coefficient becomes a time-dependent operator as well as the Lamé constants and µ take the form of the time-dependent operators and µ (11), respectively, where ν ∞ and ∞ are the non-relaxed magnitudes of the corresponding operators, and The impact occurs at t = 0 (Fig. 1). When t > 0, the displacement of the sphere's center y could be represented in the form where α is the quasi-static bearing of impactor and target's materials which is connected with the contact force by the following formula according to the generalized Hertzian law: where In formula (16) indices 1 and 2 refer to the viscoelastic beam and elastic sphere, respectively, in so doing the operators E 1 and ν 1 , which act within the contact domain, differ from operators (5), (6), and (12) valid within the other parts of the target, namely: where all rheological constants for the fractional parameter γ remain the same as for γ = 1, (Rabotnov 1948) which at γ = 1 goes over into the ordinary exponent, and operator ∋ γ (τ i ) transforms into operator ∋ * 1 (τ i ). When γ → 0, the function ∋ γ (t/τ i ) tends to the Dirac delta-function δ(t).
This distinction is connected with the fact that during the impact process there occurs decrosslinking within the domain of the contact of the beam with the impactor, resulting in more freely displacements of molecules with respect to each other, and finally in the decrease of the beam material viscosity in the contact zone (Popov et al. 2015). This circumstance allows one to describe the behaviour of the beam material within the contact domain by the standard linear solid model involving fractional derivatives, since variation in the fractional parameter (the order of the fractional derivative) enables one to control the viscosity of the beam material from its initial value at γ = 1 to its vanishing at γ = 0. Thus, the substitution of operators (5), (6), and (12) with operators (17)-(19), respectively, is quite reasonable. Now the equation of motion of the sphere could be written in the form where P(t) is defined by formula (15), while the equation of motion of the contact zone, which is considered to be rigid and is restricted by the planes z = ±a (Fig. 1) is written as Equations (23) and (25) are subjected to the following initial conditions Under the assumptions made above with respect to the contact domain, transient longitudinal and transverse waves (surfaces of strong discontinuity) propagate from the boundary of the contact zone during the impact. A certain desired function Z(z, t) behind the front of the wave surface is represented in terms of the ray series (Rossikhin  and Shitikova 1995) where Z, (k) = Z, + (k) −Z, − (k) = ∂ k Z/∂t k are the discontinuities in the k-th order derivatives with respect to time t of the desired function Z(z, t) on the wave surface, the upper signs + and − denote that the given value is calculated immediately ahead of and behind the wave front, respectively, the index α labels the ordinal number of the wave, namely: α = 1 for the longitudinal wave, and α = 2 for the transverse wave, H(t) is the Heaviside function, and G (α) is the normal velocity of the surface of discontinuity.
To determine coefficients of the ray series (27), it is necessary to differentiate the governing Eqs. (1)-(4) k times with respect to time, take their difference on the different sides of the wave surface and apply the condition of compatibility for the k + 1-order discontinuities of the function Z, which has the following form (Rossikhin and Shitikova 1995): where d/dt is the complete time-derivative of the function Z, (k) (z, t) on the moving surface of discontinuity.
Since the process of impact is a transient process, then, firstly, it is possible to limit ourselves by zeroth terms of the ray series (28), and secondly, to neglect the waves reflected from the end face of the beam considering that they reach the contact zone after impactor's rebound from the beam.
Further we shall interpret a shock wave in the beam (surface of strong discontinuity) as a layer of small thickness δ, the head front of which arrives at a certain point M with the coordinate z at the moment of time t, while the back front of the shock layer reaches this point at the moment t + t. The desired values Z(z, t) at the point M, such as velocity, generalized forces and deformations, during the time increment t change monotonically and uninterruptedly from the magnitude Z − to the magnitude Z + , in so doing within the layer, according to the condition of compatibility (28), the relationship is fulfilled, which is the more accurate the smaller the value of t.
Substituting the derivatives ∂N /∂z, ∂Q/∂z, and ∂M/∂z in (1) with −G −1Ṅ , −G −1Q , and −G −1Ṁ , respectively, integrating then the resulting equations from t to t + t, and tending t to zero, we find Substituting in (2)-(4) the derivatives ∂u/∂z, ∂w/∂z, and ∂ψ/∂z by the expressions −G −1 v, −G −1 W , and −G −1 , respectively, and writing them at the moments t and t + t, we obtain Expanding the integrals entering in (32), (34), and (36) into the Taylor series with respect to the small parameter t and limiting ourselves by the zeroth and first approximations, we have Subtracting (31), (33), and (35), respectively, from (32), (34), and (36) with due account for (37)-(39), and tending t to zero, we find From relationships (30) and (40) it is possible to find the velocities of two types of transient waves: longitudinal-flexural wave and shear wave Substituting the found velocities (41) and (42) in formulae (30) and limiting, as it has been already mentioned, by the zeroth terms of the ray series, we have Note that relationships (43) differ nothing from those for an elastic beam, since at the moment of impact a viscoelastic medium behaves as an elastic medium with the unrelaxed elastic modulus. Now it is necessary to substitute the values of N and Q defined by (43) in (25). However the governing set of two equations, (23) and (25), should involve only two unknown values, α and w, while the force N entering in (25) depends on the velocity v, as it follows from (43), and therefore v should be expressed in terms of α and w.
For this purpose we write the relationship for the stress tensor in a viscoelastic medium where summation is carried out over two repeated indices, an index after a comma labels the derivative with respect to the corresponding coordinate, σ ij and u i are the stress tensor and displacement vector components, respectively, x = x 1 , y = x 2 , z = x 3 , and δ ij is the Kronecker symbol (i, j = 1, 2, 3).
Using the procedure applied above to deduce formulae (40), from relationship (44) we obtain Multiplying (45) sequentially by k i k j and s i s j , and neglecting the press of layers within the front of the surface of strong discontinuity in the direction of the vectors k and s, i.e., considering that we have whence it follows the equality Considering the generalized geometric conditions of compatibility (Rossikhin and Shitikova 2007b) relationships (47) and (48)

could be rewritten in the form
Since in further treatment one-term ray expansions will be used, then (52) and (53) take the form or (44) . Rossikhin et al. SpringerPlus (2016) 5:206 Note that formula (54) has been presented in Landau and Lifshitz (1965) for elastic beams.
where k 1 = 4 3 √ Rd −1 . A solution of the set of two equations, involving differential equation (60) and integrodifferential equation (62), allows one to find the time-dependence of the values W and α.
Note that since the impact process is of short duration, then in the integrals entering in (62) Equations (60) and (62) are subjected to the initial conditions (26).

Solution of the problem in the case of neglecting the extension of the beam's middle surface
If the extension of the middle surface is excluded from consideration, then the set of the governing equations (60) and (62) with due account for (63) is reduced to the following: where The solution of (65) could be constructed in the form where the function C(t) could be found using the method of variation of an arbitrary constant, resulting in Substituting (67) in (62), we arrive at the governing integro-differential equation subjected to the initial conditions A solution of (68) could be found using the method of successive approximations. Thus, a particular solution of (68) neglecting its terms in square brackets has the form where Relationships (76) and (77) allow one to estimate for the limiting cases the time of impactor's rebound from the target t reb and the time t max at which the indentation attains its maximal magnitude α max , i.e., at γ = 0 = e 1 t 1 + e 2 t 2 .
(78) (78)-(83) shows that the increase in the parameter γ from 0 to 1 results in the increase of the duration of contact between the impactor and the viscoelastic target, and this increase rises with the increase in the defects of moduli e 1 and e 2 and with the decrease in the relaxation times t 1 and t 2 . Moreover, within the variation of the parameter γ from 0 to 1, the maximal magnitude of the value α, as well as the time, at which the indentation attains its maximum, increase. All enumerated peculiarities are governed by the increase in viscosity of the material, from which the viscoelastic beam is made of, with the increase of the fractional parameter γ.

Elastic target
If the term α 1/2Ẇ is omitted in (58) due to its small magnitude, i.e. the inertia of the contact domain is neglected, then this equation for an elastic beam is reduced to Substituting (84) in (59), we obtain the equation for determining α(t) or, after the substitution A =α, To find the solution of (86), it is necessary to use the initial condition The governing Eq. (85) differs from the corresponding equation presented in Vershinin (2014) not only by its coefficients (due to the fact that an incorrect formula was used for transverse deformation) but by its structure as well, namely: the multiplier α was missed out in the second term, though the inertia of the contact zone was not neglected in the cited paper (Vershinin 2014).
Equation (86) could be rewritten in the form If we neglect the term eα with respect to g (α is small) in (88), then it is reduced to where The solution of (89) can be constructed analytically in terms of a series where From (90) it is seen that all coefficients a i (i ≥ 5) and b i (i ≥ 4) are expressed in terms of the coefficients a 1 and a 2 , which are defined by three different processes being caused by the shock interaction. The coefficient a 1 is responsible for the dynamic processes arising in the beam during the propagation of shear wave, but the coefficient a 2 answers for the quasistatic process occurring at local bearing of the material due to Hertz's theory and for the dynamic processes arising in the beam during the propagation of the longitudinal wave. When g = 0, which is realized at an infinitely large speed of the shear wave propagation, the solution (90) for small α goes over into the quasi-static solution obtained by Timoshenko (1934) for the Bernoulli-Euler beam.