Plane strain and plane stress elasticity under fractional continuum mechanics

In this paper, the application of non-local fractional continuum model for plane strain and plane stress elasticity is presented. The kinematics and stress concepts are discussed, and governing equations in terms of displacements for both plane problems are defined. The numerical implementation utilising generalised finite difference method is shown in detail. Three cases are solved to indicate the role of order of fractional continua and length scale: biaxial tension, pure shear and complex state. Classical (local) solution is obtained as a special case.

derivative to a space variable; therefore, the formulation is at once non-local in spatial space. Among other similar concepts [4,7,10,12,23,41], the presented formulation has the following crucial advantages [35,38]: (i) the proposed new formulation has clear physical interpretation and is developed by analogy to general framework of CCM; (ii) we deal with finite deformations (in contrast to [4,7] where small deformations are considered only); (iii) contrary to previous works, e.g. [4,7,12], the generalised fractional measures of the deformation, e.g. fractional deformation gradients or fractional strains, have the same physical dimensions as classical ones (thus, their classical interpretation remains unchanged); (iv) characteristic length scale of the particular material is defined explicitly (as in classical non-local models [2,9,14,16,17] what is more robust than implicit techniques [30,34,42]); (v) objectivity requirements are proved; (vi) and finally, the discussed concept bases on the fractional material and spatial line elements in contrast to [12] where the whole formulation bases on fractional motion (which can be important because in more general formulations, displacement field may not exist cf. [27] Box 3.1 p. 57). Finally, because the overall concept bases on the Riesz-Caputo (RC) fractional derivative, we call it Riesz-Caputo fractional continuum.
The paper is structured as follows. Section 2 deals with RC fractional continuum description. In Sect. 3, starting from general fractional small-strain elasticity, governing equations for plane strain and stress problems are defined. Section 4 deals with detailed description of numerical implementation and the solution of three non-local plane strain elasticity problems, namely biaxial tension, pure shear and complex state.

Fractional deformation gradients
The concept of fractional kinematics differs from the classical kinematics by replacing the classical deformation gradient (first derivative of motion) with the fractional deformation gradient (fractional derivative of motion) [37]. As will be shown, fractional deformation gradient operates by analogy to classical deformation gradient. The important difference is that fractional deformation gradient is non-local; thus, all other related measures (e.g. strain tensor) are also non-local. It is because of fractional differential operator is defined over the interval, which is in contrast to classical differential operator which is defined in a point.
Another important aspect of the presented formulation is that through the fractional deformation gradient, we introduce fractional material and line elements. They are the foundations for all deformation measures, simultaneously living standard meaning of motion. In other words, every single point in material configuration transforms to a single point in current configuration (like in the CCM), but deformation is defined on conjugated with this point fractional line element. One can also think that this fractional line element includes through the fractional differential operator the information about specific material (e.g. steel, concrete, wood) structure, i.e. grain size or defects distribution.
Let us then follow the milestones of the FCM concept-for details please follow [35,[37][38][39]. The description is given in the Euclidean space. We refer to B as the reference configuration of the continuum body, while S denotes its current configuration. Points in B are denoted by X and in S by x. Coordinate system for B is denoted by {X A } with base E A , and for S, we have {x a } with base e a .
We generalise the classical deformation gradient and its inverse as follows whereF X andF x are fractional deformation gradients, X and x are length scales in B and S, respectively, φ defines the regular motion of the material body, while ϕ its inverse. We assume additionally that = X = x ; hence, the isotropic non-locality is considered. The fractional differential operator D α is defined as the RC fractional derivative. For clarity, the RC derivative, of function f (t) (t ∈ (a, b) ⊆ R and 0 < α < 1-when α is an integer, the usual definition of a derivative is used), reads [1,18]  where α > 0 denotes the real order of the derivative, D denotes 'derivative' (RC stands for Riesz-Caputo), a, t, b are so-called terminals, and the factor (2−α) (2) appears for objectivity reasons. The terminals a and b are in relation to the physical length scale of a particular material. In Eq.
The introduction of fractional deformation gradients results that one can think that deformation is observed in 'fractional' space. It means that we will operate on fractional spatial and material line, surface and volume elements. In this sense, we can introduce the following relations where dx and dX are the fractional spatial and material line elements, respectively, while dx and dX are classical spatial and material line elements, respectively. The situation is summarised in Fig. 1.
Finally, it is important that the length scale is not arbitrary due to objectivity requirements. The length scale is closely related to the choice of terminals in the fractional differential operator. It can be thought as putting the physical constraints on the obtained fractional generalisation of classical kinematics. In the presented formulation for arbitrary point in a body, say X A , the terminals for RC fractional derivative are a = X A − L 2 and b = L 2 + X A , where L is the length of the interval for fractional differentiation. It can be shown that the corresponding length scale is then [38] = L 2 .

Fractional strains and stresses
We define the finite fractional strains by analogy to the CCM. Thus, the difference in scalar products in actual and reference configurations allows us to define 4 concepts of finite strains cf. Fig. 1. Thus, one can define: where E is the classical Green-Lagrange strain tensor or its fractional counterpart (symmetric), e is the classical Euler-Almansi strain tensor or its fractional counterpart (symmetric), and ♦ F can be replaced with F orF X orF x or α F. Other fractional measures like left/right Cauchy-Green tensors or left/right stretch tensor can be defined by analogy to classical one by exchanging classical deformation gradient with fractional one. It is clear that when fractional deformation gradient is used related measure of deformation is non-local.
Infinitesimal fractional strain is obtained as in the classical set-up by considering relationship between fractional strains and fractional displacement gradients tensors and omitting higher order terms in obtained relations. Hence, infinitesimal fractional Cauchy strain tensor is defined as where and gradũ In above equations, U and u denote material and spatial displacement. The result given by Eq. (12) allows to observe crucial differences with one obtained in [22] (Eq. 9), [4] (Eq. 2.7) and [7] (Eq. 20). Thus first, as mentioned, Eq. (12) is obtained as a special case of general fractional finite strains. Second, in contrast to the previous results, in Eq. (12), the fractional derivative operates on the finite interval. Finally, the length scale is given explicitly and is in the relation with the interval over which the fractional differential operator acts. It is clear that fractional continua need the corresponding fractional stresses in general. However, the restrictions to be held by fractional stresses are analogous to the classical one, but fulfilled in auxiliary 'fractional' (or phenomenological) space, as mentioned. To understand this logic (without loss of generality), purely, mechanical problem is considered below.
Let us first introduce (following Fig. 1) the material (spatial) volume element dV (dv) or its fractional counterpart dṼ (dṽ), and the material (spatial) surface element dS (ds) or its fractional counterpart dS (ds)cf. [39]. Furthermore, for brevity, we propose the denotation ♦ (·) and one can replace (·) with classical quantity or fractional counterpart.
So, the fractional Cauchy (true) traction vectort (ñ) exerted on ds with outward normalñ is obtained as a transformation of classical Cauchy (true) traction vector t (n) exerted on ds with outward normal n. Thus, according to the diagram in Fig. 1, we havẽ As in CCM, we postulate the following relationship between fractional traction and fractional Cauchy stress tensort whereσ denotes fractional Cauchy stress tensor. Based on the above relations and the fact that we have the fundamental result that where σ is classical Cauchy stress tensor. The relation given by Eq. (19) should be interpreted that in the spatial picture, there is no need to postulate new definition for stresses, but fractional one is just transformation of the classical Cauchy measure to auxiliary 'fractional' space. This formal treatment is somehow equivalent to postulating more and more complicated material functions for classical description, because additional parameters coming from fractional formulation, namely α and (together with type of fractional derivative), play the role of new martial parameters. Furthermore, selection of a specific path in Fig. 1 causes that the remaining one can be considered utilising the concept of dual variables [20].
To show that the introduced fractional stress field satisfies analogous relations as in the classical formulation the balance of momentum in spatial description is considered. First, it is clear that the conservation of mass holds, hence or shortly where ρ 0 (ρ 0 ) is the reference mass density (fractional counterpart), ρ (ρ) is spatial mass density (fractional counterpart), and F is a Jacobian. The following relation holds where υ is a velocity, and f is a body force per unit mass. By applying the divergence theorem to Eq. (22), we have so, the local form is or in the absence of inertia forces Furthermore, the symmetry of ♦ σ (so σ orσ ) can be checked in classical manner utilising the balance of moment of momentum.
In this point, it is necessary to point out that transformation to material description is analogous to CCM and requires the use of a suitable deformation gradient ♦ F (defined previously). Finally, the result given by Eq. (24) gives justification and is analogical for the one postulated in [4] for similar CCM generalisation.

General governing equations
Let us consider the problem of static deformation under the assumption that material is linear elastic and small fractional deformation holds. The governing equations are (for purpose the denotation( ·) introduced in previous section is omitted) In the above equation, we have denoted: b is the body force, L e is the stiffness tensor, U and σ are parts of boundary where the displacements and the tractions are applied, respectively. Now assuming material isotropy, we have where μ, λ are Lemé constants. Hence, in view of Eqs. (27), (26) 1 and (26) 2 , the general problem of linear isotropic fractional elasticity in terms of displacements is

Plane stress and plane strain
To obtain fractional plane stress and plane strain problems, we follow classical restrictions imposed on Eq. (26). For numerical implementation purposes, both problems in terms of displacements are as follows: Plane stress Plane strain where for biaxial tension, for pure shear and for complex state In the above equations, reference displacement d = 1e − 6 m, while height and width of the plate equal H = P 1,y − P 3,y = 5 m and S = P 1,x − P 3,x = 5 m, respectively.  [24,29]. This scheme assumes the uniform spatial discretisation (uniformly distributed grid points), where additionally x = y. The number of rows and columns in the grid varies dependently to assumed and equals from 9 × 9 to 25 × 25 (or equivalently from = 0.625 m to = 0.208 m, because = x = y in the presented scheme), respectively. Example, grid for = 0.625 m is presented in Fig. 4. Such concept introduces additional fictitious nodes outside the body. It is important that we have assumed that the displacements u and v in these fictitious nodes (which appears due to non-locality [8]) are the same as the closest boundary displacements. The concept is presented in Fig. 4. In this way, the solution of a set of two partial fractional differential equations is replaced by the solution of linear equations.
The set of linear equations which govern the problem is obtained through the following logic.
Let us assume that φ is either u or v (φ should not be confused with motion Sect. 2). We propose the following approximation for derivatives at the particular point of interest (i, j) Calculation of D where f (n) (t k ) denotes classical n-th derivative at t = t k . • for the right-sided derivatives, we use: Now, knowing that α ∈ (0, 1) ⇒ n = 1, let us assume additionally that m = 2 ⇒ k = 1. For such a case, one can obtain the following approximations (cf. [38]): where Based on Eqs. (44) and (45), we can calculate a numerical approximation of particular derivatives in surrounding points Central differences are assumed for all first-order derivative approximations; however, these are defined on the basis of the points (i ± 0.5, j) or (i, j ± 0.5) for midpoints and on the basis of points (i ± 1, j) or (i, j ± 1) when the point of interest coincident with a grid point ∂φ In order to calculate necessary approximations of first-order derivatives ∂φ ∂ y , the indexes of row and column should be exchanged and x should be replaced by y.

Finally, the numerical approximations of
Assuming α = 1.0, we obtain A = 0.5, B = C = D = 0, and the above equations are reduced to the classical finite differences. For example, the ∂ ∂ x D x 1 φ i, j following second-order central form is obtained Notice that because of using midpoints (i ± 0.5 and j ± 0.5), the presented derivatives are completely defined on the basis of 21 points as presented in Fig. 5.
As the last part of the BVP approximation, let us consider discrete form for fractional strains calculation.

Results discussion
Three boundary value problems of a plane strain state are analysed (see Fig. 3). Numerical calculations were prepared utilising Python 2.7 programming language with scientific computing numpy/scipy library. For postprocessing, matplotlit library was used.
The particular attention is focused on the analysis of influence of fractional continua order (α) and length scale ( ). It will be shown that under the presented formulation, classical solution is obtained for α 1 (independently on ) and for 0 (independently on α).

Biaxial tension and pure shear
The two first problems consider basic states, namely biaxial tension and pure shear. We focus the attention to the analysis of strains (it is clear that because of linear constitutive relation, for both cases, the analogical results are obtained for the stress field). The classical solutions are uniform strains within the whole body. In Fig. 6, the influence of coefficient α and length scale on the strain difference ♦ ε − ε is presented. The value of coefficient α ranges from 0.0 to 1.0 with a step of 0.1, while the length scale is a function of a point grid density and ranges from 0.625 m for grid 9 × 9 to 0.208 m for grid 25 × 25. The strain values are taken from the centre point P c = (0, 0) (see Fig. 2).
Let us start the discussion with the results for pure shear (Fig. 6b). It can be observed that changing any variable causes fluctuations around the classical solution. Coefficient α approaching the value of 1.0 reduces the difference to the minimum, which is in agreement to an observation that the scheme becomes classical finite differences for α = 1.0-cf. Eq. (65). The convergence to classical solution can be also visible when the length scale decreases for α > 0.3. It is in agreement with expectations that going with 0, one should converge classical solution as well. When α < 0.3, strong oscillations appear and neither variation of α nor causes predictable change of the error.
In contrast to the pure shear problem, the error for biaxial tension does not depend to α or (see Fig. 6a). The observed small error, with no general tendency, is in range of 3e − 12 which in comparison with the classical strains value can be considered as a numerical approximation error. Thus, the problem represents the case, where fractional kinematics does not influence considerably the solution.
Concluding, for both biaxial tension and pure shear tests, the displacement field is strongly restricted by boundary conditions. It is the reason why independently whether fractional or classical descriptions were applied no significant difference in obtained displacement (so simultaneously strain and stress fields also) field is observed. The situation is different when complex state is considered-cf. following discussion.

Complex state
Similarly to the previously discussed problems, for complex state, the influence of coefficient α and length scale is analysed as well. However, the results also consider displacements and stresses, in both directions, and are presented instead of single point along paths-horizontal and diagonal denoted as γ − γ and β − β, respectively (see also Fig. 2). Figures 7, 8 and 9 show displacements, strains and stresses along γ − γ and β − β in dependency to α coefficient. In all cases, the influence of α can be observed. The closer α is to 1.0, the smaller the difference is from classical solution. There exists a limit for α, below which the v component of displacement begins to oscillate. The limit changes along with the length scale and is equal to c.a. 0.2 and 0.1 for = 0.625 m and ≥ 0.312 m, respectively. Smaller length scale also reduces the oscillation which for = 0.208 m are not observed before α = 0.0.
Analysing the results for various length scales, it can be observed that when 0, the isolines in the graphs became more and more vertical. The influence of α decreases, and we observe that for 0, the solution reaches reference one independently of assumed α. Such observations are in agreement with the previous ones as well as the fractional derivative definition [21,24,31].

Conclusions
Fractional calculus appears to be a powerful tool for modelling in mechanics (physics). In this paper, the application in continuum mechanics is presented. In a natural way, because of fractional differential operator definition, non-local formulation is obtained. It is important that because of the existence of many fractional derivatives, one can redefine the presented formalism in terms of them. In this sense, one can choose optimal formulation dependently, e.g. on material being described, it defines simultaneously the mainstream of future investigations. It is clear that classical solution is recovered as a special case.
In this paper, we present the application of RC fractional continuum to plane problems of elasticity. Governing equations, in terms of displacements, for both plane strain and plane stress are defined. Detailed algorithm for numerical solution is presented as well. Three plane strain problems are solved: biaxial tension, pure shear and complex state. The particular attention is focused on the analysis of influence of fractional continua order (α) and length scale ( ). It is shown that the classical (local) formulation is recovered for cases when order of fractional continua goes to one (α 1) or length scale tends to zero ( 0). Both parameters α and (as well as type of fractional differential operator) should be thought as a new material parameters. Their influence can be interpreted, in a phenomenological sense, as an introduction of some multiscale effects. In other words, instead of postulating more and more complicated material functions, one can apply fractional differentiation (instead of classical integer order one) resulting in the high flexibility of model, together with limited number of additional parameters.