A Generalized Version of a Low Velocity Impact between a Rigid Sphere and a Transversely Isotropic Strain-Hardening Plate Supported by a Rigid Substrate Using the Concept of Noninteger Derivatives

Abstract

A low velocity impact between a rigid sphere and transversely isotropic strain-hardening plate supported by a rigid substrate is generalized to the concept of noninteger derivatives order. A brief history of fractional derivatives order is presented. The fractional derivatives order adopted is in Caputo sense. The new equation is solved via the analytical technique, the Homotopy decomposition method (HDM). The technique is described and the numerical simulations are presented. Since it is very important to accurately predict the contact force and its time history, the three stages of the indentation process, including (1) the elastic indentation, (2) the plastic indentation, and (3) the elastic unloading stages, are investigated.

1. Introduction

The concept of noninteger order derivative has been intensively applied in many fields. It is worth nothing that the standard mathematical models of integer-order derivatives, including nonlinear models, do not work adequately in many cases. In the recent years, fractional calculus has played a very important role in various fields such as mechanics, electricity, chemistry, biology, economics, notably control theory, signal image processing, and groundwater problems; an excellent literature of this can be found in [1–9].

However, there exist a quite number of these fractional derivative definitions in the literature which range from Riemann-Liouville to Jumarie [10–17].The real problem that mathematicians face is that analytical solutions of these equations with noninteger order derivatives are usually not available. Since only limited classes of equations are solved by analytical means, numerical solution of these nonlinear partial differential equations is of practical importance. Though computer science is growing very fast, and numerical simulation is applied everywhere, nonnumerical issues will still play a large role [18–20]. In this paper a possibility of generalization of a low velocity impact between a rigid sphere and transversely isotropic strain-hardening plate supported by a rigid substrate that is generalized to the concept of noninteger derivatives order will be investigated.

There are many physical situations in which a thin plate made of strain-hardening materials resting on a rigid substrate is impacted by a rigid indenter. For example, such a phenomenon may be caused by the impact of hailstones, run way debris, or small stones on the panels of a vehicle or aircraft [21]. Although low velocity impact of a plate by a rigid indenter has been investigated by numerous researchers, the strain-hardening behaviour of the plate material has not been included in the analytical studies yet. Ollson [22] presented a one parameter nondimensional model for small mass impacts. Yigit and Christoforou [23, 24] have investigated the elastoplastic indentation phenomenon. They assumed the plate material to exhibit perfectly plastic behaviour and considered three stages for the indentation process: Hertzian elastic contact, elastic-perfectly plastic indentation, and Hertzian elastic unloading. Christoforou and Yigit [25, 26] used scaling rules for establishing a dynamic similarity between behaviours of the models and prototypes to present a model based on a linearized contact law with two nondimensional parameters that can be used for small as well as large mass impacts. In follow-up work [27], they obtained the nondimensional governing parameters of the low velocity impact response of composite plates through dimensional analysis and simple lumped-parameters models based on asymptotic solutions.

In this paper, approximated solutions for the generalized version of a low velocity impact between a rigid sphere and transversely isotropic strain-hardening plate supported by a rigid substrate will be obtained via the relatively new analytical method HDM.

The remaining of this paper is structured as follows: in Section 2, we present a brief history of the fractional derivative order and their properties. We present the basic ideal of the homotopy decomposition method for solving high order nonlinear fractional partial differential equations, its convergence and stability. We present the application of the HDM for system fractional nonlinear differential equations under investigation and numerical results in Section 4. The conclusions are then given in Section 5.

2. Brief History of Definitions and Properties

There exists a vast literature on different definitions of fractional derivatives. The most popular ones are the Riemann-Liouville and the Caputo derivatives. For Caputo, we have ^𝐶₀𝐷^𝛼_𝑥𝑓(𝑥)=1Γ(𝑛−𝛼)^𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑑^𝑛𝑓(𝑡)𝑑𝑡^𝑛𝑑𝑡,𝑛−1<𝛼≤𝑛. (1)

For the case of Riemann-Liouville we have the following definition: 𝐷^𝛼_𝑥𝑓(𝑥)=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑓(𝑡)𝑑𝑡. (2)

Guy Jumarie proposed a simple alternative definition to the Riemann-Liouville derivative: 𝐷^𝛼_𝑥𝑓(𝑥)=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛𝑥₀(𝑥−𝑡)^{𝑛−𝛼−1}𝑓(𝑡)−𝑓(0)𝑑𝑡. (3)

For the case of Weyl we have the following definition: 𝐷^𝛼_𝑥𝑓(𝑥)=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛∞_𝑥(𝑥−𝑡)^{𝑛−𝛼−1}𝑓(𝑡)𝑑𝑡. (4)

With the Erdelyi-Kober type we have the following definition: 𝐷^𝛼_0,𝜎,𝜂𝑓(𝑥)=𝑥^−𝑛𝜎1𝜎𝑥^𝜎−1𝑑𝑑𝑥^𝑛𝑥^{𝜎(𝑛+𝜂)}𝐼^𝑛−𝛼_{0,𝜎,𝜂+𝜎}𝑓(𝑥),𝜎>0. (5)

Here 𝐼^𝛼_{0,𝜎,𝜂+𝜎}𝑓(𝑥)=𝜎𝑥^{−𝜎(𝜂+𝛼)}Γ(𝛼)^𝑥₀𝑡^{𝜎𝜂+𝜎−1}𝑓(𝑡)(𝑡^𝜎−𝑥^𝜎)^1−𝛼𝑑𝑡. (6)

With Hadamard type, we have the following definition: 𝐷^𝛼₀𝑓(𝑥)=1Γ()𝑥𝑑𝑛−𝛼𝑑𝑥^𝑛𝑥₀𝑥log𝑡^{𝑛−𝛼−1}𝑓(𝑡)𝑑𝑡𝑡. (7)

With Riesz type, we have the following definition: 𝐷^𝛼_𝑥𝑓(𝑥1)=−()×12cos𝛼𝜋/2Γ(𝛼)𝑑𝑑𝑥^𝑚×^𝑥_−∞(𝑥−𝑡)^{𝑚−𝛼−1}𝑓(𝑡)𝑑𝑡+^∞_𝑥(𝑡−𝑥)^{𝑚−𝛼−1}𝑓(𝑡)𝑑𝑡. (8)

We will not mention the Grunward-Letnikov type here because it is in series form [28]. This is not more suitable for analytical purpose.

In 1998, Davison and Essex [16] published a paper which provides a variation to the Riemann-Liouville definition suitable for conventional initial value problems within the realm of fractional calculus [28]. The definition is as follows: 𝐷^𝛼₀𝑓(𝑥)=𝑑^{𝑛+1−𝑘}𝑑𝑥^{𝑛+1−𝑘𝑥}₀(𝑥−𝑡)^−𝛼Γ()𝑑1−𝛼^𝑘𝑓(𝑡)𝑑𝑡^𝑘𝑑𝑡. (9)

In an article published by Coimbra [17] in 2003, a variable-order differential operator is defined as follows: 𝐷₀^𝛼(𝑡)𝑓(𝑥)=1Γ(1−𝛼(𝑥))^𝑥₀(𝑥−𝑡)^{−𝛼(𝑡)}𝑑𝑓(𝑡)+𝑓0𝑑𝑡𝑑𝑡₊−𝑓(0⁻)𝑥^{−𝛼(𝑥)}Γ((𝑥.1−𝛼)) (10)

2.1. Advantages and Disadvantages

2.1.1. Advantages [28]

It is very important to point out that all these fractional derivative order definitions have their advantages and disadvantages; here we will include Caputo, variational order, Riemann-Liouville Jumarie, and Weyl [28]. We will examine first the variational order differential operator. Anomalous diffusion phenomena are extensively observed in physics, chemistry, and biology fields [19, 29]. To characterize anomalous diffusion phenomena, constant-order fractional diffusion equations are introduced and have received tremendous success. However, it has been found that the constant-order fractional diffusion equations are not capable of characterizing some complex diffusion processes, for instance, diffusion process in inhomogeneous or heterogeneous medium [30]. In addition, when we consider diffusion process in porous medium, if the medium structure or external field changes with time, in this situation, the constant-order fractional diffusion equation model cannot be used to well characterize such phenomenon [31, 32]. Still in some biology diffusion processes, the concentration of particles will determine the diffusion pattern [33, 34]. To solve the above problems, the variable-order (VO) fractional diffusion equation models have been suggested for use [34].

With the Jumarie definition which is actually the modified Riemann-Liouville fractional derivative, an arbitrary continuous function needs not to be differentiable; the fractional derivative of a constant is equal to zero and more importantly it removes singularity at the origin for all functions for which 𝑓(0)=constant, for instant, the exponentials functions and Mittag-Leffler functions [28].

With the Riemann-Liouville fractional derivative, an arbitrary function needs not to be continuous at the origin and it needs not to be differentiable.

One of the great advantages of the Caputo fractional derivative is that it allows traditional initial and boundary conditions to be included in the formulation of the problem [5, 12]. In addition its derivative for a constant is zero.

It is customary in groundwater investigations to choose a point on the centreline of the pumped borehole as a reference for the observations and therefore neither the drawdown nor its derivatives will vanish at the origin, as required [13]. In such situations where the distribution of the piezometric head in the aquifer is a decreasing function of the distance from the borehole, the problem may be circumvented by rather using the complementary, or Weyl, fractional order derivative [13].

2.1.2. Disadvantages [28]

Although these fractional order derivatives display great advantages, however, they are not applicable in all the situations. We will begin with the Liouville-Riemann type.

The Riemann-Liouville derivative has certain disadvantages when trying to model real-world phenomena with fractional differential equations [28]. The Riemann-Liouville derivative of a constant is not zero. In addition, if an arbitrary function is a constant at the origin, its fractional derivation has a singularity at the origin for instant exponential and Mittag-Leffler functions. Theses disadvantages reduce the field of application of the Riemann-Liouville fractional derivative.

Caputo’s derivative demands higher conditions of regularity for differentiability: to compute the fractional derivative of a function in the Caputo sense, we must first calculate its derivative. Caputo derivatives are defined only for differentiable functions while functions that have no first-order derivative might have fractional derivatives of all orders less than one in the Riemann-Liouville sense.

With the Jumarie fractional derivative, if the function is not continuous at the origin, the fractional derivative will not exist, for instance, what will be the fractional derivative of ln(𝑥) and many other ones [28].

Variational order differential operator cannot easily be handled analytically. Numerical approach is some time needs to deal with the problem under investigation.

Although Weyl fractional derivative found its place in groundwater investigation, it is still displaying a significant disadvantage; because the integral defining these Weyl derivatives is improper, greater restrictions must be placed on a function [28]. For instance, the Weyl derivative of a constant is not defined. On the other hand general theorems about Weyl derivatives are often more difficult to formulate and prove than are corresponding theorems for Riemann-Liouville derivatives.

3. Method Description [35, 36]

To illustrate the basic idea of this method, we consider a general nonlinear nonhomogeneous fractional partial differential equation with initial conditions of the following form: 𝜕^𝛼𝑈(𝑥,𝑡)𝜕𝑡^𝛼=𝐿(𝑈(𝑥,𝑡))+𝑁(𝑈(𝑥,𝑡))+𝑓(𝑥,𝑡),𝛼>0. (11)

Subject to the initial condition 𝐷^𝑘₀𝑈(𝑥,0)=𝑔_𝑘(𝑥),(),𝐷𝑘=0,…,𝑛−1_𝑛0𝑈(𝑥,0)=0,𝑛=[𝛼], (12) where 𝜕^𝛼/𝜕𝑡^𝛼 denotes the Caputo fractional order derivative operator, 𝑓 is a known function, 𝑁 is the general nonlinear fractional differential operator, and 𝐿 represents a linear fractional differential operator. The method first step here is to transform the fractional partial differential equation to the fractional partial integral equation by applying the inverse operator 𝜕^𝛼/𝜕𝑡^𝛼 on both sides of (11) to obtain 𝑈(𝑥,𝑡)=^𝑛−1_𝑗=1𝑔_𝑗(𝑥)Γ𝑡𝛼−𝑗+1^𝑗+1Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1[𝐿(𝑈((𝑈((𝑥,𝜏))+𝑁𝑥,𝜏))+𝑓𝑥,𝜏)𝑑𝜏, (13) or in general by putting 𝑓(𝑥,𝑡)=^𝑛−1_𝑗=1𝑔_𝑗(𝑥)Γ𝑡𝛼−𝑗+1^𝑗. (14)

We obtain 𝑈()()+1𝑥,𝑡=𝑇𝑥,𝑡Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1[𝐿(𝑈((𝑈((𝑥,𝜏))+𝑁𝑥,𝜏))+𝑓𝑥,𝜏)𝑑𝜏. (15)

In the homotopy decomposition method, the basic assumption is that the solutions can be written as a power series in 𝑝𝑈𝑥,𝑡,𝑝=^∞_𝑛=0𝑝^𝑛𝑈_𝑛(𝑥,𝑡),𝑈(𝑥,𝑡)=lim_𝑝→1𝑈𝑥,𝑡,𝑝, (16a)(16b)

and the nonlinear term can be decomposed as 𝑁𝑈(𝑥,𝑡)=^∞_𝑛=0𝑝^𝑛ℋ_𝑛(𝑈), (17) where 𝑝∈(0,1] is an embedding parameter. ℋ_𝑛(𝑈) is a polynomials that can be generated by ℋ_𝑛𝑈₀,…,𝑈_𝑛=1𝜕𝑛!𝜕𝑝^𝑛𝑛𝑁^∞_𝑗=0𝑝^𝑗𝑈_𝑗(𝑥,𝑡),𝑛=0,1,2…. (18)

The homotopy decomposition method is obtained by the graceful coupling of homotopy technique with Abel integral and is given by ^∞_𝑛=0𝑝^𝑛𝑈_𝑛()()=𝑝𝑥,𝑡−𝑇𝑥,𝑡Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1𝑓(𝑥,𝜏)+𝐿^∞_𝑛=0𝑝^𝑛𝑈_𝑛(𝑥,𝜏)+𝑁^∞_𝑛=0𝑝^𝑛𝑈_𝑛(𝑥,𝜏)𝑑𝜏. (19)

Comparing the terms of same powers of 𝑝 gives solutions of various orders with the first term 𝑈₀(𝑥,𝑡)=𝑇(𝑥,𝑡). (20)

4. Application of the Method to Solve the Governing Differential Equations

In this section, the analytical technique described in Section 3 is employed to obtain the solutions of the governing differential equations in each of the mentioned three contact stages. The derivation of this equation can be found in [37].

4.1. Solution of the Governing Differential Equation in the Elastic Indentation Phase

The governing equation under investigation here is given as follows: 𝜕^𝛽_𝑡𝛼(𝑡)+𝜋𝐸_𝑧𝑅1−𝑣_𝑧𝑟𝑣_𝑟𝑧𝛼ℎ𝑚²(𝑡)=0,1<𝛽≤2. (21)

Subject to the initial conditions 𝜕_𝑡𝛼(0)=𝑉₀,𝛼(0)=0. (22)

Here, 𝐸 and 𝑣 are Young’s modulus and Poisson’s ratio of the plate, respectively. 𝛼(𝑡) is elastic indentation phase; 𝑚 and 𝑉₀ are the mass of the indenter and the initial velocity, respectively; ℎ is the thickness of the plate and 𝑅 is the radius of the spherical indenter [38, 39].

Now following the description of the HDM, we arrive at the following equation: ^∞_𝑛=0𝑝^𝑛𝛼_𝑛(𝑡)−𝛼(0)𝑡−𝜕_𝑡𝛼(0)=−𝑝𝛾Γ𝛽^𝑡₀(𝑡−𝜏)^𝛽−1∞_𝑛=0𝑝^𝑛𝛼_𝑛(𝜏)²𝑑𝜏,𝛾=𝜋𝐸_𝑧𝑅1−𝑣_𝑧𝑟𝑣_𝑟𝑧.ℎ𝑚 (23)

Comparing the terms of the same power of 𝑝 we arrive at the following integral equations, which are very easier to compute: 𝑝⁰∶𝛼₀(𝑡)=𝑉₀𝑡𝑝¹∶𝛼₁(𝑡)𝛾=−Γ𝛽^𝑡₀(𝑡−𝜏)^𝛽−1𝛼²₀⋮𝑝𝑑𝜏_𝑛∶𝛼_𝑛(𝑡)𝛾=−Γ𝛽^𝑡₀(𝑡−𝜏)^{𝛽−1𝑛−1}_𝑗=0𝛼_𝑗𝛼_{𝑛−𝑗−1}𝑑𝜏,𝑛≥2. (24)

Integrating the above we obtain the following solutions: 𝛼₀(𝑡)=𝑉₀𝛼𝑡,₁(𝑡)=−𝛾𝑡^𝛽+2𝑉²₀Γ,𝛼1+𝛽₂(𝑡)=2𝛾²𝑡^3+2𝛽𝑉³₀1+𝛽2+𝛽3+𝛽Γ,𝛼4+2𝛽₃(𝑡)=4𝛾³𝑡^4+3𝛽𝑉⁴₀3×−1+𝛽3+𝛽6Γ2+𝛽²Γ+𝑡5+3𝛽^1+𝛽𝑉₀𝛾Γ7+3𝛽ΓΓΓ,𝛼1+𝛽4+2𝛽6+4𝛽₄(𝑡)=−4𝛾⁴𝑡^5+4𝛽𝑉⁵₀Γ4+𝛽Γ6+3𝛽Γ²ΓΓΓΓ×1+𝛽4+2𝛽5+3𝛽6+4𝛽7+5𝛽2𝑡^1+𝛽𝑉₀Γ−ΓΓ,𝛼𝛾Γ5+3𝛽7+5𝛽2Γ1+𝛽5+2𝛽+Γ5+3𝛽7+5𝛽₅(𝑡)=−1×3Γ³1+𝛽Γ2+2𝛽Γ²ΓΓΓ4+2𝛽5+3𝛽×Γ6+4𝛽7+5𝛽8+6𝛽⁻¹×4√𝜋𝛾⁵𝑡^6+5𝛽𝑉⁶₀×6𝑡^1+𝛽𝑉₀ΓΓΓΓ×Γ×Γ𝛾Γ4+𝛽2+2𝛽2+2𝛽×Γ4+2𝛽5+3𝛽6+3𝛽2Γ1+𝛽7+4𝛽+Γ7+5𝛽×Γ8+5𝛽−Γ𝛽+12Γ𝛽+13+𝛽²ΓΓ×ΓΓ+ΓΓΓ.4+2𝛽×Γ5+3𝛽7+3𝛽×Γ6+4𝛽+3Γ4+𝛽2+2𝛽4Γ𝛽+14+2𝛽×Γ5+2𝛽6+3𝛽2Γ4+2𝛽5+2𝛽+Γ4+𝛽5+3𝛽×Γ6+4𝛽7+4𝛽×Γ8+6𝛽 (25)

In the same manner one can obtain the rest of the components. But in this case, few terms were computed and the asymptotic solution is given by 𝛼(𝑡)=𝛼₀(𝑡)+𝛼₁(𝑡)+𝛼₂(𝑡)+𝛼₃(𝑡)+𝛼₄(𝑡)+𝛼₅(𝑡)+⋯. (26)

Remark 1. Equation (21) was solved in [37] via the homotopy perturbation method for =2. In the HPM, the initial guess or first component of the series solution may not be unique, whereas with the HDM the first component is uniquely defined as the Taylor series expansion of order 𝑛−1 (𝑛 is the order of the partial differential equation). This is one of the advantages that the HDM has over HPM.

The contact force in the elastic indentation phase may be interpreted in terms of the indentation value [37] 𝐹(𝛼(𝑡))=𝛾𝛼²(𝑡). (27)

Figures 1–6 present the approximate solution for 𝑅=0.008m, 𝑚=10⁻², 𝑉₀=5mm/s, ℎ=0.0003, 𝑣_𝑟𝑧=𝑣_𝑧𝑟=0.3, and 𝐸=75GPa. The approximate solutions of main problem have been depicted in Figures 1, 2, 3, 4, 5, and 6 which plotted according to different 𝛽 values as function of time for a fixed 𝑥 and as function of space and time.

Figure 1 

Approximate solution (26) of the governing differential equation in the elastic indentation phase for 𝛽=1.9.

Figure 2 

Approximate solution (26) of the governing differential equation in the elastic indentation phase for 𝛽=2.

Figure 3 

Approximate solution of the contact force in the elastic indentation phase (27) with 𝛽=1.9.

Figure 4 

Approximate solution of the contact force in the elastic indentation phase (27) with 𝛽=2.

Figure 5 

Surface showing the approximate solution of the governing differential equation in the elastic indentation phase equation (21) for 𝛽=1.9.

Figure 6 

Approximate solution of the contact force in the elastic indentation phase equation (27) for 𝛽=2.

4.2. Solution of the Governing Differential Equation in the Plastic Indentation Phase

The governing equation under investigation here is given as follows. 𝑚_𝑖𝜕^𝛽_𝑡𝛼(𝑡)+2𝜋𝑅𝑆_𝑦2𝛼(𝑡)−𝛼𝑡_𝑐𝑟+𝑃_𝑧𝜋𝑅1−𝑣_𝑟𝑧𝑣_𝑧𝑟ℎ𝛼(𝑡)−𝛼𝑡_𝑐𝑟²,1<𝛽≤2. (28)

Subject to the initial conditions 𝛼𝑡_𝑐𝑟=𝛼_𝑐𝑟;𝜕_𝑡𝛼𝑡_𝑐𝑟=𝑉_𝑐𝑟. (29)

Here, 𝑆_𝑦 is the yield stress, 𝑃_𝑧 is the slope of the stress-strain curve in the plastic region and it may be defined as 𝑃_𝑧=𝑛𝐸_𝑧, with 0≤𝑛≤1. Therefore, 𝑛 may be considered as a strain-hardening index. 𝑁=0 denotes a perfectly plastic behavior, whereas 𝑛=1 represents an elastic material behaviour. By increasing 𝑛 from 0 to 1, behaviour of the material approaches elastic behaviour. In addition, initial conditions of this phase or the initial velocity correspond to the values attained at the critical indentation at the end of the elastic indentation stage based on (26). For simplicity let 𝑎=4𝜋𝑅𝑆_𝑦𝑚_𝑖−2𝜋𝑃_𝑧𝑅𝛼_𝑐𝑟1−𝑣_𝑟𝑧𝑣_𝑧𝑟ℎ𝑚_𝑖𝑃,𝑏=^𝑧𝜋𝑅𝑚_𝑖1−𝑣_𝑟𝑧𝑣_𝑧𝑟ℎ,𝑃𝑐=^𝑧𝜋𝑅𝑚_𝑖1−𝑣_𝑟𝑧𝑣_𝑧𝑟ℎ𝛼𝑡_𝑐𝑟². (30)

Such that (28) can be reduced to 𝜕^𝛽_𝑡𝛼(𝑡)+𝑎𝛼(𝑡)+𝑏𝛼(𝑡)²+𝑐=0,1<𝛽≤2. (31)

Employing the HDM, we obtain the following integral equations: 𝛼₀(𝑡)=𝑡𝑉_𝑐𝑟,𝛼₁(𝑡)1=−Γ𝛽^𝑡_𝑡𝑐𝑟(𝑡−𝜏)^𝛽−1𝑎𝛼₀(𝜏)+𝑏𝛼²₀(𝜏)𝛼+𝑐𝑑𝜏,₁𝑡_𝑐𝑟=𝜕_𝑡𝛼₁𝑡_𝑐𝑟𝛼=0,_𝑛(𝑡)1=−Γ𝛽^𝑡_𝑡𝑐𝑟(𝑡−𝜏)^𝛽−1×𝑎𝛼_𝑛−1(𝜏)+𝑏^𝑛−1_𝑗𝛼_𝑗(𝜏)𝛼_{𝑛−𝑗−1}(𝜏)𝛼𝑑𝜏,_𝑛𝑡_𝑐𝑟=𝜕_𝑡𝛼_𝑛𝑡_𝑐𝑟=0,𝑛≥0. (32)

Integrating the above we arrived at the following: 𝛼₀(𝑡)=𝑡𝑉_𝑐𝑟,𝛼₁(𝑡)=−𝑡−𝑡_𝑐𝑟^𝛽𝑐1+𝛽2+𝛽+𝑡−𝑡_𝑐𝑟×𝑉₀2𝑏𝑡−𝑡_𝑐𝑟𝑉₀×Γ+𝑎2+𝛽3+𝛽⁻¹,𝑎₂(𝑡)=𝑡−𝑡_𝑐𝑟^2𝛽ΓΓΓΓ×ΓΓ1+2𝛽2+2𝛽3+2𝛽4+2𝛽𝑎𝑐Γ2+2𝛽3+2𝛽4+2𝛽+𝑡𝑉₀Γ×𝑎1+2𝛽₂+2𝑏𝑐1+𝛽Γ3+2𝛽Γ4+2𝛽+2𝑏𝑡−𝑡_𝑐𝑟𝑉₀Γ×3+𝛽2+2𝛽2𝑏𝑡−𝑡_𝑐𝑟𝑉₀Γ.3+2𝛽+𝑎Γ4+2𝛽 (33)

Using the package Mathematica, in the same manner one can obtain the rest of the components. But in this case, few terms were computed and the asymptotic solution is given by 𝛼(𝑡)=𝛼₀(𝑡)+𝛼₁(𝑡)+𝛼₂(𝑡)+𝛼₃(𝑡)+𝛼₄(𝑡)+𝛼₅(𝑡)+⋯. (34)

4.3. Solution of the Governing Differential Equation of the Unloading Phase

The governing equation of motion of the indenter mass in the unloading phase under investigation here is given as follows: 𝜕^𝛽_𝑡𝛼(𝑡)+𝜋𝑅𝐸_𝑧1−𝑣_𝑟𝑧𝑣_𝑧𝑟ℎ𝛼²(𝑡)−(1−𝑛)𝛼_𝑚−𝛼_𝑐𝑟²=0,1<𝛽≤2. (35)

Subject to the initial conditions 𝛼𝑡_𝑚=𝛼_𝑚,𝜕_𝑡𝛼𝑡_𝑚, (36) where 𝛼_𝑚 and 𝑡_𝑚 are the maximum indentation value and its relevant occurrence time, respectively. At the maximum indentation time, the velocity of the indenter becomes zero. Therefore, the values corresponding to this time may be used as initial conditions for the unloading stage [37].

Initial conditions of this phase may be obtained from solutions of the previous stage at the time of the maximum indentation. The velocity of the indenter at the time instant that it attains its maximum indentation is zero. Therefore, time of the maximum indentation may be determined by differentiating (34), with respect to time and setting the resulting equation equal to zero. Solving this equation, the time of the maximum indentation is obtained. Substituting this time into (34) yields the value of the maximum indentation as 𝛼𝑡_𝑚=𝛼₀𝑡_𝑚+𝛼₁𝑡_𝑚+𝛼₂𝑡_𝑚+𝛼₃𝑡_𝑚+𝛼₄𝑡_𝑚+𝛼₅𝑡_𝑚+⋯. (37)

For simplicity let: 𝑎=𝜋𝑅𝐸_𝑧𝑚1−𝑣_𝑟𝑧𝑣_𝑧𝑟ℎ,𝑏=𝜋𝑅𝐸_𝑧𝑚1−𝑣_𝑟𝑧𝑣_𝑧𝑟ℎ(𝑛−1)𝛼_𝑚−𝛼_𝑐𝑟². (38)

Thus (35) is reduced to 𝜕^𝛽_𝑡𝛼(𝑡)+𝑎𝛼²(𝑡)+𝑏=0,1<𝛽≤2. (39)

Following carefully the steps involved in the HDM we obtain the following integral equations: 𝛼₀(𝑡)=𝛼_𝑚𝛼₁(𝑡)1=−Γ𝛽^𝑡_𝑡𝑚(𝑡−𝜏)^𝛽−1𝑎𝛼²₀(𝜏)⋮𝛼+𝑏𝑑𝜏_𝑛(𝑡)1=−Γ𝛽^𝑡_𝑡𝑚(𝑡−𝜏)^𝛽−1𝑎^𝑛−1_𝑗=0𝛼_𝑗𝛼_{𝑛−𝑗−1}𝛼𝑑𝜏,_𝑛𝑡_𝑚=𝜕_𝑡𝛼𝑡_𝑚=0,𝑛≥1. (40)

Integrating the above we arrive at the following series solutions: 𝛼₀(𝑡)=𝛼_𝑚,𝛼₁(𝑡)=−𝑎𝑎²_𝑚+𝑏𝑡−𝑡_𝑚^𝛽Γ,𝛼1+𝛽₂(𝑡)=2𝑎𝑎_𝑚𝑎𝑎²_𝑚+𝑏𝑡−𝑡_𝑚^2𝛽Γ,𝑎1+2𝛽₃(𝑡)=−𝑎𝑎𝑎²_𝑚+𝑏𝑡−𝑡_𝑚^3𝛽×8𝑎𝑎²_𝑚Γ²1+𝛽+𝑎𝑎²_𝑚Γ×Γ+𝑏1+2𝛽₂1+𝛽Γ1+3𝛽⁻¹,𝛼₄(𝑡)=𝑎𝑎²_𝑚𝑡−𝑡_𝑚^4𝛽𝛽Γ2𝛽Γ4𝛽Γ²ΓΓ×ΓΓ×1+𝛽1+2𝛽1+4𝛽4𝛽1+2𝛽𝑎𝑎²_𝑚Γ×+𝑏1+2𝛽8𝑎𝑎²_𝑚Γ²1+𝛽+𝑎𝑎²_𝑚Γ𝑎+𝑏1+2𝛽+2²𝑎⁴_𝑚+𝑏²Γ1+𝛽Γ1+3𝛽+2𝑎𝑎²_𝑚ΓΓ.𝑏Γ2𝛽1+𝛽1+3𝛽×Γ1+4𝛽 (41)

Using the package Mathematica, in the same manner one can obtain the rest of the components. But in this case, few terms were computed and the asymptotic solution is given by 𝛼(𝑡)=𝛼₀(𝑡)+𝛼₁(𝑡)+𝛼₂(𝑡)+𝛼₃(𝑡)+𝛼₄(𝑡)+𝛼₅(𝑡)+⋯. (42)

5. Conclusion and Discussion

Low velocity impact between a rigid sphere and a transversely isotropic strain-hardening plate supported by a rigid substrate was extended to the concept of noninteger derivatives. The governing equations of the elastic indentation were obtained by Yigit and Christoforou [23, 24]. The contact was assumed to be elastic, and the stresses through the thickness were assumed to be constant. The stress expressions are only valid when no permanent deformation results due to the impact. The experimental evidence reported by Poe Jr. and Illg [39] and Poe Jr. [40] confirms the maximum value of the transverse. Normal stress has the dominant influence on the failure of a plate subjected to impact loads. The third phase is assumed to be an elastic one again.

A brief history of the fractional derivative orders was presented. Advantages and disadvantages of each definition were presented. The new equations were solved approximately using the relatively new analytical technique, the homotopy decomposition methods. The numerical simulations showed that the approximate solutions are continuous and increase functions of the fractional derivative orders. The method used to derive approximate solution is very efficient, easier to implement, and less time consuming. The HDM is a promising method for solving nonlinear fractional partial differential equations.

Conflict of Interests

The authors declare that they have no conflict interests.

Authors’ Contribution

A. Atangana and A. Ahmed made the first draft and N. Bıldık corrected and improved the final version. All the authors read and approved the final draft.

Acknowledgment

The authors would like to thank the referee for some valuable comments and helpful suggestions.