# 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

Authors: Abdon Atangana; O. Aden Ahmed; Necdet Bıldık
Journal: Abstract and Applied Analysis
 (2013)

---

## 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.


---

## Body

## 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 Section2, 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)Dxα0C(f(x))=1Γ(n-α)∫0x(x-t)n-α-1dnf(t)dtndt,=1Γ(n-α)∫0xn-1<α≤n.For the case of Riemann-Liouville we have the following definition:(2)Dxα(f(x))=1Γ(n-α)dndxn∫0x(x-t)n-α-1f(t)dt.Guy Jumarie proposed a simple alternative definition to the Riemann-Liouville derivative:(3)Dxα(f(x))=1Γ(n-α)dndxn∫0x(x-t)n-α-1{f(t)-f(0)}dt.For the case of Weyl we have the following definition:(4)Dxα(f(x))=1Γ(n-α)dndxn∫x∞(x-t)n-α-1f(t)dt.With the Erdelyi-Kober type we have the following definition:(5)D0,σ,ηα(f(x))=x-nσ(1σxσ-1ddx)nxσ(n+η)I0,σ,η+σn-α(f(x)),=x-nσ(1σxσ-1ddx)nxσ(n+η)σ>0.Here(6)I0,σ,η+σα(f(x))=σx-σ(η+α)Γ(α)∫0xtση+σ-1f(t)(tσ-xσ)1-αdt.With Hadamard type, we have the following definition:(7)D0α(f(x))=1Γ(n-α)(xddx)n∫0x(logxt)n-α-1f(t)dtt.With Riesz type, we have the following definition:(8)Dxα(f(x))=-12cos(απ/2)×{1Γ(α)(ddx)m×(∫-∞x(x-t)m-α-1f(t)dt+∫x∞(t-x)m-α-1f(t)dt)}.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:
(9)D0αf(x)=dn+1-kdxn+1-k∫0x(x-t)-αΓ(1-α)dkf(t)dtkdt.In an article published by Coimbra [17] in 2003, a variable-order differential operator is defined as follows:
(10)D0α(t)(f(x))=1Γ(1-α(x))∫0x(x-t)-α(t)df(t)dtdt+(f(0+)-f(0-))x-α(x)Γ(1-α(x)).


### 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 whichf(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 ofln(x) 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.


## 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 whichf(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 ofln(x) 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.


## 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 whichf(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 ofln(x) 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:(11)∂αU(x,t)∂tα=L(U(x,t))+N(U(x,t))+f(x,t),α>0.Subject to the initial condition(12)D0kU(x,0)=gk(x),(k=0,…,n-1),D0nU(x,0)=0,n=[α],
where ∂α/∂tα denotes the Caputo fractional order derivative operator, f is a known function, N is the general nonlinear fractional differential operator, and L 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 ∂α/∂tα on both sides of (11) to obtain
(13)U(x,t)=∑j=1n-1gj(x)Γ(α-j+1)tj+1Γ(α)∫0t(t-τ)α-1[L(U(x,τ))+N(U(x,τ))×∫0t(t-τ)α-1+f(x,τ){L(U(x,τ))+N(U(x,τ))}]dτ,
or in general by putting
(14)f(x,t)=∑j=1n-1gj(x)Γ(α-j+1)tj.We obtain(15)U(x,t)=T(x,t)+1Γ(α)∫0t(t-τ)α-1[L(U(x,τ))+N(U(x,τ))×∫0t(t-τ)α-1+f(x,τ){L(U(x,τ))+N(U(x,τ))}]dτ.In the homotopy decomposition method, the basic assumption is that the solutions can be written as a power series inp(16a)U(x,t,p)=∑n=0∞pnUn(x,t),(16b)U(x,t)=limp→1U(x,t,p),and the nonlinear term can be decomposed as(17)NU(x,t)=∑n=0∞pnℋn(U),
where p∈(0,1] is an embedding parameter. ℋn(U) is a polynomials that can be generated by
(18)ℋn(U0,…,Un)=1n!∂n∂pn[N(∑j=0∞pjUj(x,t))],n=0,1,2….The homotopy decomposition method is obtained by the graceful coupling of homotopy technique with Abel integral and is given by(19)∑n=0∞pnUn(x,t)-T(x,t)=pΓ(α)∫0t(t-τ)α-1[f(x,τ)+L(∑n=0∞pnUn(x,τ))+N(∑n=0∞pnUn(x,τ))]dτ.Comparing the terms of same powers ofp gives solutions of various orders with the first term
(20)U0(x,t)=T(x,t).


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

In this section, the analytical technique described in Section3 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:(21)∂tβα(t)+πEzR(1-vzrvrz)hmα2(t)=0,1<β≤2.Subject to the initial conditions(22)∂tα(0)=V0,α(0)=0.Here,E and v are Young’s modulus and Poisson’s ratio of the plate, respectively. α(t) is elastic indentation phase; m and V0 are the mass of the indenter and the initial velocity, respectively; h is the thickness of the plate and R is the radius of the spherical indenter [38, 39].Now following the description of the HDM, we arrive at the following equation:(23)∑n=0∞pnαn(t)-α(0)t-∂tα(0)=-pγΓ(β)∫0t(t-τ)β-1(∑n=0∞pnαn(τ))2dτ,γ=πEzR(1-vzrvrz)hm.Comparing the terms of the same power ofp we arrive at the following integral equations, which are very easier to compute:
(24)p0:α0(t)=V0tp1:α1(t)=-γΓ(β)∫0t(t-τ)β-1α02dτ⋮pn:αn(t)=-γΓ(β)∫0t(t-τ)β-1∑j=0n-1αjαn-j-1dτ,n≥2.Integrating the above we obtain the following solutions:(25)α0(t)=V0t,α1(t)=-γtβ+2V02Γ(1+β),α2(t)=2γ2t3+2βV03(1+β)(2+β)(3+β)Γ(4+2β),α3(t)=4γ3t4+3βV043(1+β)(3+β)×(-6Γ(2+β)2Γ(5+3β)+t1+βV0γΓ(7+3β)Γ(1+β)Γ(4+2β)Γ(6+4β)),α4(t)=-4γ4t5+4βV05Γ(4+β)Γ(6+3β)Γ2(1+β)Γ(4+2β)Γ(5+3β)Γ(6+4β)Γ(7+5β)×(2t1+βV0γΓ(5+3β)Γ(7+5β)-(2Γ(1+β)Γ(5+2β)+Γ(5+3β)Γ(7+5β)){2t1+βV0γΓ(5+3β)Γ(7+5β)}),α5(t)=-1×(3Γ3(1+β)Γ(2+2β)Γ2(4+2β)Γ(5+3β)×Γ(6+4β)Γ(7+5β)Γ(8+6β){3Γ3(1+β)Γ(2+2β)Γ2(4+2β)Γ(5+3β)})-1×(4πγ5t6+5βV06×(6t1+βV0γΓ(4+β)Γ(2+2β)Γ(2+2β)×Γ(4+2β)Γ(5+3β)Γ(6+3β)×(2Γ(1+β)Γ(7+4β)+Γ(7+5β))×Γ(8+5β)-Γ(β+1)×(2Γ(β+1)(3+β)Γ2(4+2β)×Γ(5+3β)Γ(7+3β)×Γ(6+4β)+3Γ(4+β)Γ(2+2β)×(4Γ(β+1)Γ(4+2β)×Γ(5+2β)Γ(6+3β)+(2Γ(4+2β)Γ(5+2β)+Γ(4+β)Γ(5+3β))×Γ(6+4β))Γ(7+4β))×Γ(8+6β)){6t1+βV0γΓ(4+β)Γ(2+2β)Γ(2+2β)}).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)α(t)=α0(t)+α1(t)+α2(t)+α3(t)+α4(t)+α5(t)+⋯.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 n-1 (n 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)F(α(t))=γα2(t).Figures1–6 present the approximate solution for R=0.008m, m=10-2, V0=5mm/s, h=0.0003, vrz=vzr=0.3, and E=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 x 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.(28)mi∂tβα(t)+2πRSy[2α(t)-α(tcr)]+PzπR(1-vrzvzr)h(α(t)-α(tcr))2,1<β≤2.Subject to the initial conditions(29)α(tcr)=αcr;∂tα(tcr)=Vcr.Here,Sy is the yield stress, Pz is the slope of the stress-strain curve in the plastic region and it may be defined as Pz=nEz, with 0≤n≤1. Therefore, n may be considered as a strain-hardening index. N=0 denotes a perfectly plastic behavior, whereas n=1 represents an elastic material behaviour. By increasing n 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
(30)a=4πRSymi-2πPzRαcr(1-vrzvzr)hmi,b=PzπRmi(1-vrzvzr)h,c=PzπRmi(1-vrzvzr)hα(tcr)2.Such that (28) can be reduced to
(31)∂tβα(t)+aα(t)+bα(t)2+c=0,1<β≤2.Employing the HDM, we obtain the following integral equations:(32)α0(t)=tVcr,α1(t)=-1Γ(β)∫tcrt(t-τ)β-1[aα0(τ)+bα02(τ)+c]dτ,α1(tcr)=∂tα1(tcr)=0,αn(t)=-1Γ(β)∫tcrt(t-τ)β-1×[aαn-1(τ)+b∑jn-1αj(τ)αn-j-1(τ)]dτ,αn(tcr)=∂tαn(tcr)=0,n≥0.Integrating the above we arrived at the following:(33)α0(t)=tVcr,α1(t)=-((t-tcr)β({(t-tcr)β({V0(2b(t-tcr)V0+a(2+β))}c(1+β)(2+β)+(t-tcr)}c(1+β)(2+β)+(t-tcr)×V0(2b(t-tcr)V0+a(2+β)){(t-tcr)β({V0(2b(t-tcr)V0+a(2+β))}c(1+β)(2+β)+(t-tcr)}))×(Γ(3+β))-1,a2(t)=(t-tcr)2βΓ(1+2β)Γ(2+2β)Γ(3+2β)Γ(4+2β)×(acΓ(2+2β)Γ(3+2β)Γ(4+2β)+tV0Γ(1+2β)×((a2+2bc(1+β))Γ(3+2β)Γ(4+2β)+2b(t-tcr)V0(3+β)Γ(2+2β)×({(a2+2bc(1+β))Γ(3+2β)}2b(t-tcr)V0Γ(3+2β)+aΓ(4+2β){(a2+2bc(1+β))Γ(3+2β)}))).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)α(t)=α0(t)+α1(t)+α2(t)+α3(t)+α4(t)+α5(t)+⋯.


### 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:(35)∂tβα(t)+πREz(1-vrzvzr)h(α2(t)-(1-n)(αm-αcr)2)=0,1<β≤2.Subject to the initial conditions(36)α(tm)=αm,∂tα(tm),
where αm and tm 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)α(tm)=α0(tm)+α1(tm)+α2(tm)+α3(tm)+α4(tm)+α5(tm)+⋯.For simplicity let:(38)a=πREzm(1-vrzvzr)h,b=πREzm(1-vrzvzr)h(n-1)(αm-αcr)2.Thus (35) is reduced to
(39)∂tβα(t)+aα2(t)+b=0,1<β≤2.Following carefully the steps involved in the HDM we obtain the following integral equations:(40)α0(t)=αmα1(t)=-1Γ(β)∫tmt(t-τ)β-1[aα02(τ)+b]dτ⋮αn(t)=-1Γ(β)∫tmt(t-τ)β-1[a∑j=0n-1αjαn-j-1]dτ,αn(tm)=∂tα(tm)=0,n≥1.Integrating the above we arrive at the following series solutions:(41)α0(t)=αm,α1(t)=-(aam2+b)(t-tm)βΓ(1+β),α2(t)=2aam(aam2+b)(t-tm)2βΓ(1+2β),a3(t)=-(a(aam2+b)(t-tm)3β×(8aam2Γ2(1+β)+(aam2+b)Γ(1+2β)){a(aam2+b)(t-tm)3β})×(Γ2(1+β)Γ(1+3β))-1,α4(t)=aam2(t-tm)4ββΓ(2β)Γ(4β)Γ2(1+β)Γ(1+2β)Γ(1+4β)×({×Γ(1+β)Γ(1+3β)Γ(1+4β){(aam2+b)Γ(1+2β)})}Γ(4β)Γ(1+2β)×((aam2+b)Γ(1+2β)×(8aam2Γ2(1+β)+(aam2+b)Γ(1+2β)+2(a2am4+b2)Γ(1+β)Γ(1+3β))+2aam2bΓ(2β)Γ(1+β)Γ(1+3β)×Γ(1+4β){(aam2+b)Γ(1+2β)})).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)α(t)=α0(t)+α1(t)+α2(t)+α3(t)+α4(t)+α5(t)+⋯.


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

The governing equation under investigation here is given as follows:(21)∂tβα(t)+πEzR(1-vzrvrz)hmα2(t)=0,1<β≤2.Subject to the initial conditions(22)∂tα(0)=V0,α(0)=0.Here,E and v are Young’s modulus and Poisson’s ratio of the plate, respectively. α(t) is elastic indentation phase; m and V0 are the mass of the indenter and the initial velocity, respectively; h is the thickness of the plate and R is the radius of the spherical indenter [38, 39].Now following the description of the HDM, we arrive at the following equation:(23)∑n=0∞pnαn(t)-α(0)t-∂tα(0)=-pγΓ(β)∫0t(t-τ)β-1(∑n=0∞pnαn(τ))2dτ,γ=πEzR(1-vzrvrz)hm.Comparing the terms of the same power ofp we arrive at the following integral equations, which are very easier to compute:
(24)p0:α0(t)=V0tp1:α1(t)=-γΓ(β)∫0t(t-τ)β-1α02dτ⋮pn:αn(t)=-γΓ(β)∫0t(t-τ)β-1∑j=0n-1αjαn-j-1dτ,n≥2.Integrating the above we obtain the following solutions:(25)α0(t)=V0t,α1(t)=-γtβ+2V02Γ(1+β),α2(t)=2γ2t3+2βV03(1+β)(2+β)(3+β)Γ(4+2β),α3(t)=4γ3t4+3βV043(1+β)(3+β)×(-6Γ(2+β)2Γ(5+3β)+t1+βV0γΓ(7+3β)Γ(1+β)Γ(4+2β)Γ(6+4β)),α4(t)=-4γ4t5+4βV05Γ(4+β)Γ(6+3β)Γ2(1+β)Γ(4+2β)Γ(5+3β)Γ(6+4β)Γ(7+5β)×(2t1+βV0γΓ(5+3β)Γ(7+5β)-(2Γ(1+β)Γ(5+2β)+Γ(5+3β)Γ(7+5β)){2t1+βV0γΓ(5+3β)Γ(7+5β)}),α5(t)=-1×(3Γ3(1+β)Γ(2+2β)Γ2(4+2β)Γ(5+3β)×Γ(6+4β)Γ(7+5β)Γ(8+6β){3Γ3(1+β)Γ(2+2β)Γ2(4+2β)Γ(5+3β)})-1×(4πγ5t6+5βV06×(6t1+βV0γΓ(4+β)Γ(2+2β)Γ(2+2β)×Γ(4+2β)Γ(5+3β)Γ(6+3β)×(2Γ(1+β)Γ(7+4β)+Γ(7+5β))×Γ(8+5β)-Γ(β+1)×(2Γ(β+1)(3+β)Γ2(4+2β)×Γ(5+3β)Γ(7+3β)×Γ(6+4β)+3Γ(4+β)Γ(2+2β)×(4Γ(β+1)Γ(4+2β)×Γ(5+2β)Γ(6+3β)+(2Γ(4+2β)Γ(5+2β)+Γ(4+β)Γ(5+3β))×Γ(6+4β))Γ(7+4β))×Γ(8+6β)){6t1+βV0γΓ(4+β)Γ(2+2β)Γ(2+2β)}).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)α(t)=α0(t)+α1(t)+α2(t)+α3(t)+α4(t)+α5(t)+⋯.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 n-1 (n 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)F(α(t))=γα2(t).Figures1–6 present the approximate solution for R=0.008m, m=10-2, V0=5mm/s, h=0.0003, vrz=vzr=0.3, and E=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 x 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.(28)mi∂tβα(t)+2πRSy[2α(t)-α(tcr)]+PzπR(1-vrzvzr)h(α(t)-α(tcr))2,1<β≤2.Subject to the initial conditions(29)α(tcr)=αcr;∂tα(tcr)=Vcr.Here,Sy is the yield stress, Pz is the slope of the stress-strain curve in the plastic region and it may be defined as Pz=nEz, with 0≤n≤1. Therefore, n may be considered as a strain-hardening index. N=0 denotes a perfectly plastic behavior, whereas n=1 represents an elastic material behaviour. By increasing n 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
(30)a=4πRSymi-2πPzRαcr(1-vrzvzr)hmi,b=PzπRmi(1-vrzvzr)h,c=PzπRmi(1-vrzvzr)hα(tcr)2.Such that (28) can be reduced to
(31)∂tβα(t)+aα(t)+bα(t)2+c=0,1<β≤2.Employing the HDM, we obtain the following integral equations:(32)α0(t)=tVcr,α1(t)=-1Γ(β)∫tcrt(t-τ)β-1[aα0(τ)+bα02(τ)+c]dτ,α1(tcr)=∂tα1(tcr)=0,αn(t)=-1Γ(β)∫tcrt(t-τ)β-1×[aαn-1(τ)+b∑jn-1αj(τ)αn-j-1(τ)]dτ,αn(tcr)=∂tαn(tcr)=0,n≥0.Integrating the above we arrived at the following:(33)α0(t)=tVcr,α1(t)=-((t-tcr)β({(t-tcr)β({V0(2b(t-tcr)V0+a(2+β))}c(1+β)(2+β)+(t-tcr)}c(1+β)(2+β)+(t-tcr)×V0(2b(t-tcr)V0+a(2+β)){(t-tcr)β({V0(2b(t-tcr)V0+a(2+β))}c(1+β)(2+β)+(t-tcr)}))×(Γ(3+β))-1,a2(t)=(t-tcr)2βΓ(1+2β)Γ(2+2β)Γ(3+2β)Γ(4+2β)×(acΓ(2+2β)Γ(3+2β)Γ(4+2β)+tV0Γ(1+2β)×((a2+2bc(1+β))Γ(3+2β)Γ(4+2β)+2b(t-tcr)V0(3+β)Γ(2+2β)×({(a2+2bc(1+β))Γ(3+2β)}2b(t-tcr)V0Γ(3+2β)+aΓ(4+2β){(a2+2bc(1+β))Γ(3+2β)}))).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)α(t)=α0(t)+α1(t)+α2(t)+α3(t)+α4(t)+α5(t)+⋯.


## 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:(35)∂tβα(t)+πREz(1-vrzvzr)h(α2(t)-(1-n)(αm-αcr)2)=0,1<β≤2.Subject to the initial conditions(36)α(tm)=αm,∂tα(tm),
where αm and tm 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)α(tm)=α0(tm)+α1(tm)+α2(tm)+α3(tm)+α4(tm)+α5(tm)+⋯.For simplicity let:(38)a=πREzm(1-vrzvzr)h,b=πREzm(1-vrzvzr)h(n-1)(αm-αcr)2.Thus (35) is reduced to
(39)∂tβα(t)+aα2(t)+b=0,1<β≤2.Following carefully the steps involved in the HDM we obtain the following integral equations:(40)α0(t)=αmα1(t)=-1Γ(β)∫tmt(t-τ)β-1[aα02(τ)+b]dτ⋮αn(t)=-1Γ(β)∫tmt(t-τ)β-1[a∑j=0n-1αjαn-j-1]dτ,αn(tm)=∂tα(tm)=0,n≥1.Integrating the above we arrive at the following series solutions:(41)α0(t)=αm,α1(t)=-(aam2+b)(t-tm)βΓ(1+β),α2(t)=2aam(aam2+b)(t-tm)2βΓ(1+2β),a3(t)=-(a(aam2+b)(t-tm)3β×(8aam2Γ2(1+β)+(aam2+b)Γ(1+2β)){a(aam2+b)(t-tm)3β})×(Γ2(1+β)Γ(1+3β))-1,α4(t)=aam2(t-tm)4ββΓ(2β)Γ(4β)Γ2(1+β)Γ(1+2β)Γ(1+4β)×({×Γ(1+β)Γ(1+3β)Γ(1+4β){(aam2+b)Γ(1+2β)})}Γ(4β)Γ(1+2β)×((aam2+b)Γ(1+2β)×(8aam2Γ2(1+β)+(aam2+b)Γ(1+2β)+2(a2am4+b2)Γ(1+β)Γ(1+3β))+2aam2bΓ(2β)Γ(1+β)Γ(1+3β)×Γ(1+4β){(aam2+b)Γ(1+2β)})).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)α(t)=α0(t)+α1(t)+α2(t)+α3(t)+α4(t)+α5(t)+⋯.


## 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.



---

*Source: 671321-2013-03-31.xml*