Analysis of the Viscoelastic Sphere Impact against a Viscoelastic Uflyand-Mindlin Plate considering the Extension of Its Middle Surface

In the present paper, the problem on impact of a viscoelastic sphere against a viscoelastic plate is consideredwith due account for the extension of plate’s middle surface and local bearing of sphere and plate’s materials via the Hertz theory. The standard linear solid models with conventional derivatives and with fractional-order derivatives are used as viscoelastic models, respectively, outside and within the contact domain. As a result of impact, transient waves (surfaces of strong discontinuity) are generated in the plate, behind the wave fronts of which up to the boundaries of the contact domain the solution is constructed in terms of one-term ray expansions due to short-time duration of the impact process. The motion of the contact zone occurs under the action of extension forces acting in the plate’s middle surface, transverse force, and the Hertzian contact force. The suggested approach allows one to find the time-dependence of the impactor’s indentation into the target and the Hertzian contact force.


Introduction
Nowadays fractional calculus is widely used in different fields of science and technology [1][2][3][4][5][6][7][8][9][10][11], including various dynamic problems of mechanics of solids and structures [12].The problems of impact response of thin-walled structures belong to the most challenging ones among the dynamic problems.The review of papers dealing with the application of fractional calculus models for solving the impact problems could be found in [12,13].
However, to the authors' knowledge, the impact response of viscoelastic plates utilizing the fractional derivative models has been treated only in a few articles.Thus, the problems of impact of an elastic rod or sphere against a viscoelastic Kirchhoff-Love plate or against a viscoelastic Uflyand-Mindlin plate have been considered in [14] and [15][16][17], respectively.In the case of the classical plate, the method of the Green function construction was used as the method of solution [14].In the case of the Uflyand-Mindlin plate, the ray method was utilized, according to which the solution behind the fronts of transient waves (surfaces of strong discontinuity) is constructed in terms of truncated power series with variable coefficients.Moreover, the contact force has been taken into account either according to the Hertz contact theory [15] or in terms of linear approximations allowing one to apply Laplace transform technique [16,17].
In all enumerated above papers, the extension of the plate's middle surface, as well as viscoelastic properties of an impactor, has not been considered.However, as experimental research of the behaviour of flexural thin-walled structures under low and middle velocity impacts shows [18], the contribution of extensional deformations into the impact process is rather essential and, therefore, it should be taken into account.Thus, the impact of elastic bodies upon elastic beams and plates with consideration for the transverse deformations and extension of a middle surface was studied in [19], and recently this approach has been considered when investigating the dynamic response of a viscoelastic Timoshenko-type beam impacted by an elastic sphere [20].

Shock and Vibration
In the present paper, we will generalize the approaches developed recently in [20,21] when investigating the dynamic response of a viscoelastic beam impacted by an elastic sphere with due account for beam's middle surface extension and by a viscoelastic sphere ignoring the middle surface extension, respectively, over the problem of theoretical analysis of the dynamic response of the viscoelastic plate by considering the influence of the middle surface extension on the process of deformation of the viscoelastic plate during the impact by a viscoelastic spherical impactor.
In both papers [20,21], the viscoelastic properties of the beam have been described by the standard linear solid model with integer time derivatives.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 finely in the decrease of the beam material viscosity in the contact zone.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.
It is known that viscoelastic properties of the material could be changed during its life-cycle due to different external effects, among them aging and impact loads [22], and/or temperature and irradiation [23].Dynamic mechanical studies of irradiated polyethylene were presented by Merrill et al. [23], wherein it was found the influence of irradiation on viscoelastic characteristics of polyethylene.Experimental study of concrete aging effect on the contact force and contact time during impact interaction of an elastic rod with a viscoelastic beam was carried out by Popov et al. [22], and a good correlation with the analytical model proposed in [24] was declared.
In the present paper, following [20,21,24] we will study the dynamic response of a viscoelastic Uflyand-Mindlin plate impacted by a viscoelastic sphere considering the extension of plate's middle surface.The standard linear solid model with conventional derivatives and with fractionalorder derivatives will be used as viscoelastic models, respectively, outside and within the contact domain.As a result of impact, transient waves (surfaces of strong discontinuity) are generated in the plate, behind the wave fronts of which up to the boundaries of the contact domain the solution is constructed in terms of one-term ray expansions due to short-time duration of the impact process.The motion of the contact zone is analyzed under the action of extension forces acting in the plate's middle surface, transverse forces, and the Hertzian contact force using the algebra of Rabotnov's fractional operators.

Problem Formulation and Governing Equations
Let a viscoelastic sphere of radius  and mass  move with the constant velocity  0 towards a viscoelastic plate along the -axis, which is perpendicular to the plane of the plate.Impact occurs at  = 0.As a result of transverse impact against the viscoelastic plate, transient waves, the longitudinal wave, and the wave of transverse shear propagate along the plate, the fronts of which are the surfaces of strong discontinuity.These waves are cylindrical surfaces-strips, the generators of which are parallel to the normal to the median plane, while the guides located in the middle plane are the circumferences extending with the normal velocities  () .The index  indicates the ordinal number of the wave:  = 1 for the longitudinal wave and  = 2 for the transverse wave.
A certain desired function  behind the fronts of wave surfaces Σ (1) and Σ (2) could be represented in terms of the ray series [25] where [, () ] () = [  /  ] () = , +() () − , −() () are discontinuities in the th-order time-derivatives of the desired function  on the waves surfaces Σ  ; that is, at  = / () ,  is the polar radius, the upper indices "+" and "−" denote that the values are calculated immediately ahead of and behind the wave fronts, respectively, and ( − / () ) is the unit Heaviside function.
Since the impact process is of short duration, then it is possible to restrict only by zero-order terms, that is, =/ ()  ( −   () ) . ( In order to find the discontinuities in the desired values, first it is necessary to write the governing equations describing the dynamic behaviour of the viscoelastic isotropic Uflyand-Mindlin plate in the polar coordinates (, ) where   and   are the forces acting in the plate plane in the and -directions, respectively,   and   are the bending moments,   is the transverse force,  is the plate deflection,   is the displacement along the radius,   is the rotational angle of the normal in the -direction,  is the density, ℎ is the plate thickness, μ and Ẽ are shear and Young's operators, respectively, σ is the Poisson's operator, D = Ẽℎ 3 /12(1 − σ2 ) is the operator of the cylindrical rigidity,  =  2 /12 is the shear coefficient, and overdots denote the time derivatives, and then to apply the kinematic condition of compatibility proposed in [26] for physical components of thin bodies to ( 3) and ( 5), where / is the Thomas -derivative [27] of the function preassigned on the moving surface.Note that condition (11) is valid in the case when spatial coordinates coincide with the ray coordinates [28].
In order to simplify the procedure for utilizing condition (11), let us interpret the surface of strong discontinuity as a layer of small thickness , within which the desired value  changes monotonically and continuously from the magnitude  + to the magnitude  − .Suppose that the ahead and back fronts of the shock layer arrive at a certain point with the fixed radius at the time instants  and  + Δ, respectively, where Δ is small.Inside the shock layer the following relationship is fulfilled.Writing (3) and ( 5) within the shock layer and considering (12), we obtain where V  = u  , and  = ẇ.Thus we have chosen only those equations of motions which will be utilized in further treatment.
Integrating ( 13) and ( 14) over  from  to  + Δ and then tending Δ → 0 yield Relationships ( 15) and ( 16) are called the dynamic conditions of compatibility.During the deduction of these conditions we neglect the terms   and   in (3) and   in (5), as well as the term / in (12), since the integrals from these terms over  from  to  + Δ tend to zero as Δ → 0.
In order to write ( 6) and ( 10) in discontinuities, it is necessary to decode the operators entering in these equations.

Decoding of Operators
For decoding the operators involved in ( 6) and (10), it is a need to specify two any operators, for example, the operator of volume extension-compression K and shear operator μ.As numerous experiments with volume stresses and strains show [29], for the majority of materials operator K is a constant value, that is, where  ∞ is a certain constant.
In further treatment the standard linear solid model with conventional derivatives and with fractional-order derivatives will be used as viscoelastic models, respectively, outside and within the contact domain, because during the impact process there occurs decrosslinking within the domain of the contact of the plate with the impactor, resulting in more freely displacements of molecules with respect to each other and finely in the decrease of the plate material viscosity in the contact zone.This circumstance allows us to describe the behaviour of the plate 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 plate material.
Thus for convenience of presentation, first we will decode fractional-order operators valid inside the contact domain, and then tending fractional parameter to unit in the final relationships we will obtain the corresponding formulae for the viscoelastic plate outside the contact domain in terms of the conventional viscoelasticity.
So within the contact domain of the impactor with the viscoelastic target, we would preassign this connectedness between the deviators of stress   and strain   by the standard linear solid model involving fractional derivatives [12], that is, or where   and   are the relaxation and retardation times, respectively,  (0 <  ≤ 1) is the fractional parameter,  0 is the relaxed shear modulus, is the Riemann-Liouville fractional derivative of the -order [12], Γ(1 − ) is the Gamma-function, and () is a certain function.
From (19) it follows that operator μ is given in the form or Here the operator has been introduced, which has the name of the dimensionless Rabotnov fractional operator [13].
If we consider formula which is proved by the direct check-up [13], in relationship (22) and introduce the following notation where  ∞ is the nonrelaxed shear modulus, then as a result we obtain where Now let us calculate the Young operator according to the formula We will start with the operator where We will seek the inverse operator (3 ∞ + μ) −1 in the form where   and   are yet unknown constants.Further we multiply operators ( 28) and ( 29), consider their property (3 ∞ + μ)(3 ∞ + μ) −1 = 1, and utilize the formula which is proved by the direct check-up [13].
As a result we obtain Equating to zero each bracket in (31), we are led to a set of two equations for defining   and whence it follows Finally we substitute operators ( 26) and ( 29) in ( 27) and consider relationship (30).As a result we have Considering that according to the set of (32) the first bracket in ( 34) is equal to zero, while the second bracket is as follows: and we finally obtain where ) < 1, and   =   .Now let us calculate the compliance operator J = Ẽ−1 which is defined in the form where ]  and   are yet unknown constants.To determine them, we will utilize the property JẼ = 1 and consider relationship (30).Then we obtain Equating to zero each bracket in (38) yields Solving the set of (39), we find Now we could calculate the Poisson's operator σ according to formula (17), which could be rewritten in the form Considering (36), from formula (41) we have For further treatment we should know the following operators: In order to calculate the operators in the right-hand side of ( 43) and (44), we assume that they have the following form: where ,  1 and ,  2 are yet unknown constants.
Equating the right sides of relationships ( 43), ( 45) and (44), (46), reducing the obtained expressions to the common denominator, and considering (30) yield Vanishing to zero the expressions in square brackets in (47), we determine unknown constants Now we could calculate the operator For this purpose, we substitute ( 36), (45), and ( 46) in (49) with due account for (30).As a result we obtain where Finally we would calculate operator Ẽσ(1 − σ2 ) −1 considering formulas (30) and (48) Further for the convenience we will represent the Rabotnov fractional operator in another form multiplying the numerator and denominator of (23) by  −    , where is the fractional integral.
As a result we obtain It has been shown in [13] that interpreting the righthand side part of (54) as a sum of the infinite descending geometrical progression with the denominator  = −   −  , we have or where is Rabotnov's fractional exponential function [29,30], which at  = 1 goes over into the conventional exponent, that is, All formulas which have been obtained for the viscoelastic fractional derivative standard linear solid model remain valid for the standard linear solid model with conventional derivatives.
Since out of the contact domain viscoelastic features of plate's material are described via the standard linear solid model with conventional derivatives, then it is possible to rewrite ( 6) and (10) with due account for the derived above operators putting  = 1.As a result we have In order to write relationships (59) in discontinuities, let us substitute there   / and / by − −1 V  and − −1 , respectively, according to (12), and then write the obtained relationships at the moments  and  + Δ.Thus, (60) (61) Expanding the integrals entering in ( 61) and ( 63) in Taylor series in terms of the small parameter Δ and limiting ourselves by the zero and first approximations yield where () may coincide with V  (),   (), or  −1 () +   (), and  3 =   .Subtracting (61) from ( 60) and ( 63) from (62), considering relationship (64) as well as that and tending Δ to zero, we have From expressions ( 15) and ( 65), (16), and (66) it is possible to find velocities of two types of transient waves, namely, longitudinal-flexural wave and transverse shear wave: Substituting the found velocities (67) and (68) in formulas ( 15) and ( 16), and limiting, as it has been already mentioned, by the zeroth terms of the ray expansions, we obtain   = −ℎ (1)  ∞ V  , (69) Relationships ( 69) and ( 70) by their form do not differ from the corresponding relationships for an elastic plate, with the only difference that a subindex ∞ has appeared at the velocity notation, which indicates that the velocities are expressed in terms of the nonrelaxed modulus of the viscoelastic model (see Appendix).

Equations of Motion of the Contact Domain and the Impactor
At  > 0 the displacement of the sphere's center () could be represented in the form where  is the impactor's indentation due to the local bearing of target and impactors' materials within the contact domain.
Then the equation of motion of the part of the plate being in contact with the sphere and the equation of the sphere have the form where () is the radius of the contact domain, and () is the contact force.Equations ( 72) and (73) could be solved with due account for formulas (69) and (70), as well as considering the following initial conditions: The relationships following from the generalized contact Hertz theory should be added to (72) and (73), where and in doing so the indices "1" and "2" refer to the viscoelastic plate and viscoelastic sphere, respectively.The operators valid within the contact domain are based on the fractional derivative standard linear solid models.Integrating (73) yields In order to express the contact force () in terms of displacement (), it is needed to decode the operator k.The procedure for decoding this operator as well as the construction of its inverse operator has been recently described in detail in Rossikhin et al. [21].Thus, it has been shown [21] that operator k could be written in the form where and unknown constants   and   could be determined from the set of equations wherein Considering (79), the contact force is defined as where  ∞ = 4 √ /3ae ∞ .

Shock and Vibration
Utilizing (83), it is possible to rewrite (78) in the form The values of   and   defined by ( 69) and ( 70), respectively, should be substituted in (72); then this equation together with (84) will give us a set of two governing equations involving only two unknown values,  and .However, the force   entering in (72) depends on the velocity V  according to (69); thus it is necessary to express it in terms of  and .
With this purpose in mind, we write the stress tensor   for a three-dimensional viscoelastic medium where   are components of the displacement vector, summation is carried out over two repeated indices, and   is the Kronecker symbol, ,  = 1, 2, 3, an index after a comma labels the derivative with respect to the coordinates  1 = ,  2 = , and  3 = .
Considering formula (26), as well as the expression for operator where ]  = (2/3)( ∞ / ∞ )]  , which follows from the formula let us rewrite (85) in the form Using the procedure applied above for deducing (65) and (66), from (88) we obtain Applying the geometric condition of compatibility where ]  are the components of the unit vector ⃗ ] normal to the wave surface, and   are the components of the unit vector ⃗  normal to the plate, to (89), yields Note that condition (90) differs from the Hadamard condition of compatibility, since (90) involves the second term taking the transverse shear deformation of the plate into account.
Multiplying (91) by   , we find Considering that for the cylindrical wave-strip propagating in the plate [V  ]]  = [V  ] and [ , ] = [ , ] in the cylindrical set of coordinates (, , ) and that  ∞ ( ∞ +  ∞ ) −1 =  ∞ (1 −  ∞ ) −1 , we are led to the following relationship: wherein the value [ , ], in its turn, is defined by the formula Then the discontinuity [V  ] could be written as and the value   due to (69) takes the form Considering formulas (70), ( 76), (83), and (96), as well as the relationship which is obtained from formula (12), if we substitute there the function  by the function , we could rewrite (72) in the form where Note that since the impact process is of short duration, then in the integrals entering in (84) and (98) Consider (100) and ( 84) and (98) are reduced to Adding formula (71) to ( 101) and (102), we obtain the closed set of three equations in three unknown functions, namely, (), (), and ().

Approximate Analytical Solution of the Governing Set of Equations
Now we will consider several particular cases.
as it has been shown in [21], then in (101) and (102) the terms involving the sums vanish to zero; therefore, (102) is reduced to while (101) could be retained in the form of (73).
If we neglect the inertia of the contact domain due to its smallness, then it could be rewritten as Consider relationships and introduce a new variable Equation ( 73) with due account of (105) could be rewritten in the form where (/) = α .If we ignore the deformation of the middle surface, that is, putting  = 0, then (108) is reduced to where  1 =  ∞ /.Note that (109) is the Abel equation of the second kind [31], and its solution was found in [25,32] in the form where It has been emphasized in [25,32] that coefficients  1 and  2 define two different processes being caused by the shock interaction of the impactor and the target: the coefficient  1 is responsible for the dynamic processes generating in the plate during propagation of the surfaces of the discontinuity, while the coefficient  2 describes the quasistatic processes occurring at the local bearing of the material due to Hertzian theory.If  1 → 0, what corresponds to an infinite large velocity of transverse shear wave propagation (Kirchhoff-Love plate), then solution (110) for small  goes over into the quasistatic solution [18].
If  → ∞ in (108), then we are led to the well-known equation which describes the impact interaction of colliding elastic bodies according to the Hertz theory [33].
If  ̸ = 0, then since  is a small value, (108) could be rewritten in the form Let us put  =  0  in the right-hand side of (112) and integrate it.As a result we obtain From (113) we determine the refined duration of contact where  0 reb = 2 −1 ∞ , as well as the refined time at which  attains its maximal value where  0 max =  −1 ∞ , and the maximal indentation where  0 max = (1/2) 0  −1 ∞ .

The Case
where Δ = ∑ Comparing formulas (114) and ( 119), ( 115) and (120), as well as (116) and (121), the following conclusion could be made: all characteristic values calculated at  = 1 are larger than those calculated at  = 0, and all characteristic values calculated at other magnitudes of  (0 <  < 1) locate between these two corresponding limiting values in order of increasing magnitudes of the value .