Gradient-enriched finite element methodology for axisymmetric problems

Due to the axisymmetric nature of many engineering problems, bi-dimensional axisymmetric finite elements play an important role in the numerical analysis of structures, as well as advanced technology micro/nano-components and devices (nano-tubes, nano-wires, micro-/nano-pillars, micro-electrodes). In this paper, a straightforward C0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {C}^{0}$$\end{document}-continuous gradient-enriched finite element methodology is proposed for the solution of axisymmetric geometries, including both axisymmetric and non-axisymmetric loads. Considerations about the best integration rules and an exhaustive convergence study are also provided along with guidances on optimal element size. Moreover, by applying the present methodology to cylindrical bars characterised by a circumferential sharp crack, the ability of the present methodology to remove singularities from the stress field has been shown under axial, bending, and torsional loading conditions. Some preliminary results, obtained by applying the proposed methodology to notched cylindrical bars, are also presented, highlighting the accuracy of the methodology in the static and fatigue assessment of notched components, for both brittle and ductile materials. Finally, the proposed methodology has been applied to model the unit cell of the anode of Li-ion batteries showing the ability of the methodology to account for size effects.


Introduction
The origin of the finite element method (FEM) dates back to the mid-twentieth century, in the aerospace industry [13,57] where many relevant structures, especially rockets, can be considered as axisymmetric; hence the stress analysis of complex axisymmetric structures has always been of primary importance. This encouraged the development of axisymmetric finite elements in the first half of the 1960s [58]. However, even if axisymmetric finite elements have been originally developed for aerospace applications, they have found a wide application also in other fields (mechanical industry, transport, and containment of fluids), since many structural components are characterised by an axisymmetric geometry, such as pressure and containment vessels, pipes, cooling towers, tunnels, and shafts.
Speaking about aeronautical and mechanical components, then, fracture and fatigue problems are of great interest and an accurate estimation of the stress state is essential for a reliable assessment. Here gradient elasticity theories come into play; in fact, as extensively discussed and proven in the literature (see, for instance, [10,15,17,19,20,42,46]), they are able to overcome the deficiencies of classical elasticity in accurately determining the stress fields in specific problems, especially those for which it is necessary to account for microand nanoscale effects in the description of the overall macroscopic response. In particular, gradient-enriched theories are able to remove singularities from (micro-)stress and (micro-)strain fields in the neighbourhood of sharp crack tips and, in general, they have a smoothing effect on the stresses in the presence of stress concentrators, such as notches and holes. Moreover, they can interpret size effects in a robust and effective manner, in agreement with experiments.
The idea of enriching the equations of classical elasticity with higher-order spatial gradients was proposed by Cauchy in the early 1850s, with the aim of studying more accurately the behaviour of discrete lattice models [22][23][24]. However, this idea remained dormant until the beginning of the twentieth century when the Cosserat brothers [26] proposed to enrich the kinematics of the three-dimensional continuum equations with three micro-rotations and to include the couple-stresses in the equations of motion.
Nevertheless, the most significant revival of gradient theories took place in the 1960s when several efforts were spent in the extension of the Cosserat theory as well as couple stress theories. In particular, Toupin [56] proposed to improve the theory of classical elasticity by enriching the strain energy function with the full first gradient of the strain. Few years later, Mindlin suggested to define both kinetic and deformation energy density also in terms of micro-strains [43] and enriched the theory of classical elasticity by including the second gradient of the strain [44]. Classical elasticity was further extended by Green and Rivlin [38] by including all the gradients of the strain. The aforementioned research works, aiming to include a mathematically complete set of higher-order gradients into the classical elasticity equations, resulted in rather complicated theories.
Finally, about 20 years later, gradient elasticity theories saw a second significant renaissance in particular through the work of Eringen [29][30][31] and Aifantis et al. [1,10,46], leading to simpler theories containing only the higher-order terms needed to describe the phenomena of interest. In particular, in 1993 Ru and Aifantis [46] proposed a split operator allowing the solution of the original fourth-order Aifantis theory as an un-coupled sequence of two systems of second-order p.d.e., as described in Sect. 3, decreasing the continuity requirements from C 1 to C 0 , allowing a straightforward finite element implementation of the theory. Furthermore, as explained in [17], the stress-based formulation of the Ru-Aifantis theory presents significant advantages if compared to the displacement-based and strain-based formulations of the same theory. These are the main reasons why in this paper, as well as in [19], the stress-based Ru-Aifantis theory has been chosen for the finite element implementation.
In the literature, it is possible to find several applications of gradient elasticity to axisymmetric problems. In particular, strain-gradient theories such as the Mindlin theory have been successfully applied to study different axisymmetric problems such as stress concentrations at spherical inclusions and cavities [21,25], circular holes [32,33], Boussinesq problem [35,36]. Other researchers, instead, provided analytical solutions to axisymmetric problems such as thick-walled cylinder [34], annulus under internal and external pressure, and circular hole in infinite body [12] starting from the gradient elastic formulation proposed by Aifantis et al. [1,10,11].
Although gradient elasticity has been widely applied to axisymmetric problems leading to several analytical solutions for simple benchmark problems, the complex geometries of modern engineering components make analytical solution of practical problems nearly impossible and, therefore, numerical solutions become essential. However, as far as the authors are aware, no complete finite element implementation of gradient elasticity for axisymmetric problems has been proposed until now.
In this paper, a reduced-gradient version of the gradient-enriched finite element methodology presented in [17,19,53], for plane stress/strain and three-dimensional problems, has been developed for axisymmetric problems, with the aim of providing a comprehensive finite element methodology based on the Ru-Aifantis theory, applicable to any kind of axisymmetric problem. This includes even the non-trivial case of nonaxisymmetric loads, where all the circumferential variables as well as their derivatives with respect to the angular coordinate θ cannot be neglected, while displacements and forces must be expressed in terms of Fourier series. Appropriate numerical integration schemes have also been established, and comprehensive convergence studies have been carried out, allowing recommendations on optimal element size. The proposed methodology has also been used to analyse notched cylindrical bars subject to both static and fatigue loadings, showing that gradient elasticity could potentially be a very accurate and effective tool for the static and fatigue assessment of notched components.
After a brief description of the theoretical background of axisymmetric analysis and gradient elasticity theory provided, respectively, in Sects. 2 and 3, in Sect. 4 a possible C 0 -continuous finite element implementation of the Ru-Aifantis theory (with reduced gradient dependence) is proposed for axisymmetric solids subject to both axisymmetric and non-axisymmetric loads. In Sect. 5, details about the best integration rules to adopt are provided, while in Sect. 6, the proposed methodology is applied to solve the problem of an hollow cylinder subject to internal pressure and the numerical results compared to the analytical solution proposed by Gao and Park [34]. In Sects. 7 and 8, the convergence of the proposed methodology is studied for problems without and with singularities, respectively; furthermore, in Sect. 8, recommendations on optimal element size are provided. Section 9 is dedicated to applications of the proposed methodology to different problems, which highlight the main advantages of the present gradient-enriched methodology. In particular, these advantages consist of the stress smoothing in the presence of stress risers (removal of singularities from the stress fields around the tips of sharp cracks), which allows a more accurate static and fatigue assessment of components characterised by stress concentrators and, last but not least, the ability to describe size effects.

Axisymmetric loading
It is well known from the literature (see, e.g., [59]) that in axisymmetric solids subject to axisymmetric loading, displacements, strains, and stresses are all independent of the circumferential coordinate θ of the cylindrical coordinate system defined in Fig. 1. This particularity enables the study of such a problem by simply considering the generic plane section of the solid along its axis of revolution z (shaded in Fig. 1), subject to in-plane loading.
The displacement field is described by two components only, since the circumferential component u θ = 0 for axisymmetry reasons, which are functions of the radial (r ) and axial (z) coordinates only: (1) In what concerns the strains, they can be collected in the strain vector where (ignoring the spatial dependence for notational simplicity) Finally, in the case of linear elastic materials, the stress field is linked to the strain field through the following constitutive relation: where E and ν are the Young's modulus and Poisson's ratio, respectively.

Non-axisymmetric loading
In the case of an axisymmetric solid subject to non-axisymmetric loads, the circumferential component of the displacements u θ must be taken into account in addition to the radial and axial components. As a consequence of this, none of the strain and stress components can be considered null; in particular, the strain vector assumes the following form (ignoring the spatial dependence for notational simplicity) where in addition to Eq. (3), we have and the stress vector is now defined as However, it is still possible to solve this problem as a quasi-bidimensional problem by expressing the load and displacement components through Fourier series (see, e.g. [58,60]).

Gradient elasticity theory and the Aifantis GradEla model
Gradient elasticity is a generic name given to a family of generalised continuum elasticity theories which suggest to take into account the effects of the micro-structure on the mechanical behaviour of a component, by including higher-order strain (or stress) in the constitutive equation. Such inclusion results to overcome many of the limits of classical elasticity in describing particular phenomena, amongst which are the removal of singularities from strain and stress fields at crack tips, the smoothing of these fields at stress concentrators, as well as the capture of size effects. This ability is due to the introduction of internal length parameters (the phenomenological multipliers of the higher-order gradient terms) which, in turn, appear in the solution of corresponding boundary value problems. The magnitude of these internal length parameters (proportional to the size of the underlying micro-structure of the material)-in relation to a characteristic specimen dimensioncontrols the component behaviour and allows for the descriptions of the aforementioned non-standard effects.
In [15], it is argued that Laplacian-based gradient theories are the most versatile. One of the most renowned Laplacian-based theories is the GradEla model developed by Aifantis et al. [1,10,46], which represents a very simple and effective approach, characterised by just one additional parameter, where the constitutive relations are enriched by means of the Laplacian of the strain as follows: where σ i j and ε kl are, respectively, the stress and strain tensors, C i jkl is the constitutive tensor, and is the internal length scale. The inclusion of the Laplacian is not arbitrary, but it appears naturally when relating local (micro) and average (macro) response and expanding the corresponding non-local integral in a Taylor series for appropriate forms of the influence kernel. This leads to the following fourth-order equilibrium equations: where u k is the displacement vector and b i are the body forces. Ru and Aifantis [46] developed a factorisation procedure which allows the solution of Eq. (9) as a decoupled sequence of two systems of second-order partial differential equations (p.d.e.), considerably simplifying the finite element implementation of this model. The first system of second-order p.d.e. to be solved is the standard classical elasticity equilibrium C i jkl u c k, jl + b i = 0 ( 1 0 ) from which it is possible to determine the field of the local displacements u c i that will be used as source term in the following second set of second-order p.d.e.
where u g i represents the field of the non-local (or gradient-affected) displacements. Finally, as shown in [17,39,40] (and also discussed in [3,5]), Eq. (11) can be redefined in terms of stresses as σ where σ g i j represents the field of the gradient-enhanced stresses.

Implementation aspects
The Ru-Aifantis theory described in Sect. 3 has been already implemented in a unified finite element methodology for both two-and three-dimensional problems [19], showing good convergence properties as well as the ability to remove singularities from the stress field in the neighbourhood of stress risers. While the aforementioned theory is based on the use of the full gradient of the stresses, in this paper a simplified version of the Ru-Aifantis theory, characterised by the use of a reduced gradient of the stresses, is developed and implemented in a straightforward C 0 finite element framework for the analysis of axisymmetric solids subject to both axisymmetric and non-axisymmetric loads. The choice of using a reduced gradient allows to obtain a simplified yet reliable methodology, more easily implementable in a finite element framework, without renouncing to the ability to remove singularities from the stress fields (see [37] for another example of reduced-gradient theory).

Axisymmetric loads
The index notation used in Sect. 3 is valid in Cartesian coordinates, while axisymmetric problems are usually addressed in cylindrical coordinates as in the following of the present paper. Therefore, in order to derive the finite element equations for axisymmetric problems unambiguously, we will now depart from the index notation used in the previous section, adopting a tensorial notation. Furthermore, two functionals, one for each of the two steps of the Ru-Aifantis theory, are defined as follows: for the first step and for the second, where u c and ε c are, respectively, the local displacement and local strain vectors, σ g is the non-local stress vector, C is the constitutive matrix, defined in Eq. (4), S = C −1 is the compliance matrix, b is the vector of the body forces in the domain Ω, t is the vector of the prescribed traction on the free portion n of the boundary. Imposing the stationarity of the first functional, that is, δW 1 = 0, the usual global system of standard elasticity equations is obtained: where d c is the vector of the nodal local displacements, Repeating the same procedure for the second functional, W 2 , leads to where s g is the vector of the nodal non-local stresses and N σ is just an expanded version of the shape function matrix N u , in order to accommodate all four non-local stress components, so that σ g = N σ s g . At this point, solving Eq. (15) for d c , and using the result as source term in Eq. (16), the nodal values of the gradient-enriched stresses s g can be easily calculated.

Non-axisymmetric loads
The finite element equations for the case of axisymmetric solids subject to non-axisymmetric loads can be determined by following the same process described in Sect. 4.1.
Regarding the first step of the Ru-Aifantis theory, the functional defined by Eq. (13) can be considered by taking into account that now the constitutive matrix is C defined in Eq. (7) and the displacements (including now also the circumferential component) are described by means of Fourier series as with where h is the order of the harmonic, while the superscripts s and a indicate, respectively, the symmetric and anti-symmetric components of the displacements (and of the loads) with respect to the θ = 0 axis.
The local strains are defined as follows: and substituting Eq.
where B s u,h = LΘ s h N u and B a u,h = LΘ a h N u are the two strain-displacement matrices for symmetric and antisymmetric displacements/loads, respectively.
After these considerations, the first step of the Ru-Aifantis theory still consists in Eq. (15), with the difference that now also the nodal forces are expressed in terms of Fourier series as and the stiffness matrix is defined as where k hl is the stiffness contribution of the hth and lth harmonics, given by Considering that dV = r dr dzdθ = r dθ dA, it can be easily demonstrated that which means that the symmetric and anti-symmetric terms are decoupled, as it should be, given the orthogonality of the Fourier series. Furthermore, it can also be proved that Hence, the off-diagonal terms of the global stiffness matrix K are all null and each harmonic can be treated separately, allowing application of the principle of superposition to determine the final result.
In what concerns the second step of the Ru-Aifantis theory, considering for simplicity only the symmetric components, defining the following functional: where and imposing the stationarity of W 2 , that is, δ W 2 = 0, the following solving system of equations is obtained: which can be easily numerically solved for the non-local stress s g,s h . The same system must be solved also for the anti-symmetric components (if needed). Finally, once the solutions for all the necessary harmonics are calculated, the global solution can be determined through superposition of the effects.

Numerical integration
For the numerical solution of Eqs. (15), (16), and (28), the Gauss quadrature rule has been used. As summarised in Table 1, in what concerns Eq. (15) (first step), since it coincides with the p.d.e. of classical elasticity, the standard number of integration points has been considered, in particular the bi-quadratic serendipity quadrilaterals have been under-integrated.
For the numerical solution of Eqs. (16) (second step with axisymmetric loads) and (28) (second step with non-axisymmetric loads), instead, higher integration rules are needed, except for the bi-linear quadrilateral elements, if = 0 (as reported in Table 1). On the other hand, if = 0, the quadrature rules used for the first step can be adopted for all the elements without any rank deficiency problem, similar to the 2D plane strain/stress and 3D problems, as described in [19]. The analysed problem consists of a thick-walled hollow cylinder of inner radius a = 0.5 m, outer radius b = 1.5 m, and length L = 8 m, subject to an internal pressure p i = 10 MPa (Fig. 2). The material is characterised by a Young's modulus E = 1000 MPa, Poisson's ratio ν = 0.25, and the length scale has been set as = 0.1 m. Concerning the boundary conditions, for the first step (Eq. 15) the cylinder is simply supported in the axial direction at both ends ( Fig. 2), while for the second step (Eq. 16) homogeneous natural boundary conditions are taken throughout, since, as also described in [15], this choice is the most widely accepted amongst the scientific community when dealing with gradient elasticity. For symmetry of the cylinder with respect to a plane normal to the z-axis, only half of the cross section has been modelled using 16 × 64 bi-linear and bi-quadratic quadrilateral elements and the double of linear and quadratic triangular elements (note that for triangular elements, the meshes are obtained by subdividing each quadrilateral in two triangles). Since the displacements determined according to the proposed methodology correspond to the classical elastic displacements, the numerical solution for the displacements has been compared with the following analytical solution of classical elasticity known from the literature [54]: Regarding the stresses, instead, the gradient-enriched stress fields numerically obtained by applying the developed methodology have been compared to the analytical solutions proposed by Gao and Park [34]. In Fig. 3, the displacement and stress fields, obtained by using bi-quadratic quadrilateral elements, are compared with the relative analytical solutions. Similar results have also been obtained with all the other implemented elements.
From Fig. 3, it can be observed that the numerical estimation of the radial displacements u r perfectly matches the correspondent analytical counterpart. Regarding the stress components, instead, the numerical solutions obtained by applying the proposed methodology show some differences when compared with the analytical solutions proposed by Gao and Park [34]. In particular, it is possible to observe that both the numerical σ g θθ and σ g zz stress profiles present qualitative and quantitative differences if compared to their analytical counterparts. This is due to the fact that while the methodology presented in this paper is stress-based, the analytical solutions proposed by Gao and Park [34] are obtained through a displacement-based formulation and therefore the stress components are calculated following different procedures. Furthermore, the displacement-based formulation used by Gao and Park [34] requires just two higher-order boundary conditions (in particular, they set homogeneous natural boundary conditions of the radial stress at the inner and outer surfaces) while the stressbased methodology proposed in this paper requires more higher-order boundary conditions (homogeneous natural boundary conditions of all the stress components at the inner and outer surfaces have been chosen). Therefore, the aforementioned differences are also due to a different choice of the boundary conditions. Finally, it can be easily explained why σ g zz obtained by applying the proposed methodology should be constant. To do so, let us consider first the exact stress solutions of classical elasticity known from the literature [54]: In particular, it can be observed that σ zz is constant along the radial coordinate, and therefore, due to the way the proposed methodology post-processes the stress fields, the longitudinal stress component is not affected by any gradient enrichment, leading to σ g zz = σ zz = const. Only quantitative differences can be seen, instead, in the radial stress component, due to the fact that, while in the proposed methodology the stress components are uncoupled, in the analytical solution proposed by Gao and Park [34], the stress components are coupled.
Furthermore, the analytical solution proposed by Gao and Park [34] is obtained through a full gradient methodology, while, as already explained in Sect. 4, the methodology proposed in the present paper is based on the use of a reduced gradient.
Due to the aforementioned differences (although justified), the analytical solution proposed by Gao and Park cannot be used as benchmark solution for the convergence study performed in the following section where the reference solution will be approximated through Richardson extrapolation [45].

Convergence study
The convergence behaviour of the implemented elements in the case of both axisymmetric and nonaxisymmetric loads has been studied. For this purpose, the L2-norm error defined as where σ e and σ c are, respectively, the extrapolated and calculated values of the stresses, has been plotted against the number of degrees of freedom (nDoF). From now on, the reference solutions have been approximated through Richardson extrapolation [45].
As already explained in [19], the stresses are determined as primary variable, by solving the second step of the Ru-Aifantis theory (either Eq. (16) or (28)) and not as secondary variables, like in standard finite element methodologies based on classical elasticity. This means that, although the error conserves its proportionality with the nDoF, the theoretical convergence rate for the stresses is different from the one determined in classical elasticity. In particular, Helmholtz equations such as Eqs. (16) and (28) are characterised by the following proportionality [41]: where n is the polynomial order.
Equation (32) clearly shows that the e 2 -nDoF curve in a bi-logarithmic system of axis is represented by a straight line, whose slope represents the convergence rate p of the numerical solution to the extrapolated solution. In particular, the theoretical convergence rate is p = −1 for linear elements and p = −1.5 for quadratic elements.
In what concerns the case of axisymmetric load, the thick-walled hollow cylinder problem presented in Sect. 5 has been considered. The domain has been modelled by using all the implemented elements, starting from a 4 × 16 mesh and refining the mesh, by doubling the number of elements in both directions 4 times, up to a 64 × 256 mesh. The length parameter has been set as = 0.1 m.
Regarding the case of axisymmetric solids subject to non-axisymmetric loads, two different problems have been analysed. The first one consists of a cylindrical bar of radius R = 1 m and length L = 8 m, subject to a bending moment M = 10 6 Nm (Fig. 4, left), while the second one consists of the same cylinder subject now to a torsional moment T = 10 6 Nm (Fig. 4, right). In both cases, Young's modulus E = 1000 MPa, Poisson's ratio ν = 0.25, and the length scale has been set as = 0.1 m.
For symmetry of the cylinder with respect to a plane normal to the z-axis, in both cases, only half of the cross section has been modelled (Fig. 4). Also, in this case all the implemented elements have been used to discretise the domain, starting from a 4 × 16 mesh and performing the same mesh refinement described in the case of axisymmetric loads. Concerning the boundary conditions, in both problems, for the first step (Eq. 15 with K defined by Eq. 22), homogeneous essential boundary conditions have been imposed, so that u z = 0 (bending problem) and u z = u θ = 0 (torsional problem) along the axis of symmetry (Fig. 4), while for the second step (Eq. 28) homogeneous natural boundary conditions have been chosen, as in the previous case.
In Fig. 5, the L2-norm error is plotted against the nDoF in a bi-logarithmic system for the three analysed problems, from which it is possible to appreciate that all the implemented elements produce numerical solutions characterised by convergence rates in good agreement with the theoretical predictions when applied to any kind of axisymmetric problem.

Convergence in the presence of singularities and recommendations on optimum element size
As is well known [1,3,5,10,46], one of the most important abilities of gradient elasticity is the removal of singularities from the strain and stress fields such as those obtained at the tips of sharp cracks when classical elasticity is used. See also the analyses in [15,17,19], as well as a more recent discussion in [6,7]. Since the presence of cracks reduces drastically the convergence rate of standard finite element methodologies, very fine meshes are required in the neighbourhood of the stress riser to ensure that convergence of the solution is reached, with significant increasing in computational cost. Hence, the analysis of the convergence behaviour of the proposed gradient-enriched finite element methodology is of great interest. Furthermore, an accurate error estimation has allowed the formulation of recommendations on optimal element size also for axisymmetric finite elements that together with the recommendations provided in [19] for plane and three-dimensional finite elements constitutes a complete and useful meshing guideline for gradient-enriched finite elements.
To analyse the aforementioned aspects, three different problems (Fig. 6) consisting of a cylindrical bar characterised by the presence of a circumferential crack and subject, respectively, to axial load F = 10 6 N (axisymmetric), bending moment M = 10 6 Nm, and torque T = 10 6 Nm (non-axisymmetric loads) have been studied. The cylinder has a radius R = 1 m, a length L = 8 m, and the crack is 0.25 m deep. Material properties, boundary conditions, and meshes are taken as in Sect. 7.
As in Sect. 7, the convergence rate is determined as the slope of the straight line obtained by plotting in a bi-logarithmic graph the L2-norm error, defined by Eq. (31), versus the nDoF (Fig. 7). From the literature [59], it is known that in the presence of singularities, the error on classical stresses is proportional to the nDoF as where λ = 0.5 for a nearly closed crack, which leads to a theoretical convergence rate p = −0.25 for both linear and quadratic elements.  Fig. 7, it is possible to observe that, when singularities are involved, both linear and quadratic elements show approximatively the same convergence rate (in accordance with Eq. (33)). Furthermore, the proposed gradient-enriched methodology produces a significant improvement in terms of convergence rate; in particular, the solutions of the three analysed problems are all characterised by a convergence rate almost three times higher than the correspondent theoretical value (typical of standard finite element methodology based on classical elasticity). As already explained in [19], this improvement is due to two main factors: -removal of singularities from the stress fields; -gradient-enriched stresses are primary variables, instead of secondary variables (as in standard classical elasticity-based finite element methodologies). Analysing then the error affecting the numerical solutions (determined through Eq. (31)), it is possible to provide recommendations on optimal element size. Table 2 summarises the ratios between the element size and the length scale that ensure an error of about 5% or lower. The recommendations summarised in Table 2 highlight the fact that using the proposed methodology, even in the presence of singularities, it is possible to obtain accurate solutions using a fairly coarse mesh, reducing consistently the computational efforts.

Removal of singularities
As already mentioned in Sect. 1, one of the features of the proposed methodology is the ability to remove singularities from the stress field, such as those emerging near the tip of sharp cracks. This ability can be easily shown by considering the problems presented in Sect. 8 and comparing the stress profiles obtained by applying the proposed methodology with those produced by classical elasticity. In Fig. 8, in fact it can be seen that while the classical elastic stress fields are characterised by an unbounded peak in correspondence with the crack tip, the gradient-enriched ones converge to a unique finite solution upon mesh refinement (the stress profiles presented in Fig. 8 were obtained by using bi-quadratic quadrilateral elements, similar results are produced by all the other implemented elements).

Static and fatigue assessment of notched bars
The aforementioned ability to remove singularities, or more in general the stress smoothing in correspondence with stress concentrators, leads to more physically realistic results and therefore to more accurate static and fatigue assessments of components presenting stress risers. In particular, it allows the static and fatigue assessments, by considering the values of the relevant stresses directly at the crack/notch tip and not into the analysed body, as it happens in most of the existing assessment procedures, with evident simplifications in the assessment process.
In [16,48], similarities and differences between gradient elasticity and the Theory of Critical Distances (TCD) (for an in-depth description of this theory, see [52]) are investigated, showing the advantages of these two theories in the static and fatigue assessment of cracked components. Furthermore, in [48] a relation linking the two length parameters, characterising the aforementioned theories has been proposed, namely where L is the characteristic length of the TCD, defined as in [52]. The proposed methodology has been applied to notched cylindrical bars (Fig. 9), subject to both static and fully reversed cyclic loadings (load ratio R = −1), in order to show the accuracy of the proposed methodology, with defined through Eq. (34), in the determination of both static and fatigue strength of notched components (for more results, see also [20]).
The accuracy of the proposed methodology has been tested against a wide range of materials (both brittle and ductile) and notch root radii.
All the relevant data about the analysed problems are summarised in Tables 3 and 4, where σ 0 is the inherent material strength, σ th is the ultimate tensile nominal stress, Δσ 0 is the plain fatigue limit stress range, and Δσ th is the range of the nominal fatigue strength (all the nominal stresses are referred to the gross section of the specimens).
The accuracy of the proposed methodology in estimating both the static strength and the fatigue limit of notched cylindrical bars has been verified by defining the following errors: for static and fatigue problems, respectively, where σ eff and Δσ eff are, respectively, the stress (for static problems) and the stress range (for fatigue problems) numerically obtained at the notch tip. In Fig. 10, the aforementioned errors are plotted for the different analysed materials, from which it is possible to appreciate the accuracy of the proposed methodology, with errors mainly ranging between −10% and +30% (dashed lines). Considering the usual error band ±20%, it is possible to observe that the present methodology produces results with errors falling into a band of the same size, but 10% more conservative.     In previous publications by Aifantis et al. [2,4,14,27,28], it was shown that the GradEla model can conveniently capture the occurrence of elastic size effects in small one-dimensional objects and specimens with holes commonly used for advanced technology applications. More recently, its extension to describe size effects due to thermomechanical, chemomechanical, and electromechanical couplings has been reviewed in [7]. The purpose of this section is to expand this discussion for an emerging technological problem of internal stress development and capacity fade in Lithium-ion (Li-ion) rechargeable battery electrodes by elaborating on the classical elasticity treatment first used for this problem in the pioneering work of [8] (see also related chapters and references in [9]). In fact, nowadays, the majority of electronic devices, such as laptops, tablets, smart phones, and digital cameras, are powered by Li-ion batteries. Li-ion batteries have replaced their older sisters nickel-cadmium (Ni-Cd) and nickel-metal-hydride (Ni-MH) batteries due to their unique characteristics, such as lightness, high energy density, and non-toxicity [8].
In Li-ion batteries, the active sites in the anode and (to a lesser extent) in the cathode are subject to significant expansions upon lithium insertion, which leads to a severe damage of the surrounding glass/ceramic matrix and consequently to a deterioration of the electrochemical properties of the electrodes. It is then evident how critical an accurate estimation of the stresses in the matrix is for a subsequent damage evaluation.
However, the dimensions characterising this problems are usually of the order of nanometres, which makes classical elasticity not accurate enough for the estimation of the stresses, due to its incapability to capture size effects, particularly significant at the nanoscale. In this section, the proposed gradient-enriched methodology is applied to study the aforementioned problem, in order to show its ability to easily describe size effects, obtaining a more accurate size-dependent estimation of the stresses, with a small additional computational effort.
The active sites of the electrodes can be fabricated in different shapes like thin layers, spheres, or fibres. In this section, the third shape is considered and analysed.
Fibre-like electrodes can be considered as a series of unit cells formed by a cylindrical Li-ion active site of radius r = a, surrounded by a glass or ceramic matrix (inert with respect to Li) of radius r = b.
Two different problems have been analysed: in first instance, only the outer matrix has been modelled as an hollow cylinder, and after that, both the active site and the surrounding matrix have been modelled as two nested cylinder made of two different materials. The material properties are taken as E g = 75 GPa, E s = 41 GPa, ν g = 0.25, and ν s = 0.33 (see [8]), where the subscript g refers to the glass matrix and the subscript s to the active site. Furthermore, the glass matrix and the active site are characterised by an internal length g and s , respectively. Concerning the geometry dimensions, the inner and outer radii have been defined, respectively, as a = γ a 0 and b = γ b 0 , where a 0 and b 0 are, respectively, the inner and outer radii of the domain, as given in  [8], and different sizes of the annular ring have been analysed, by varying the value of the scaling coefficient γ from 0.2 to 4.0.
In both cases, due to their axisymmetry, the bi-dimensional problems shown in Fig. 11 have been modelled using bi-quadratic quadrilateral elements. Regarding the boundary conditions, for the first set of equations (15), they have been chosen as u r (a) = Δ − δ and u r (b) = 0.
In fact, it can be considered that, at the inner surface of the matrix (r = a), since active site and matrix are in contact, when the active site reaches its maximum expansion Δ, the matrix pushes back a distance δ. On the other hand, at the outer matrix-matrix boundary (r = b) of adjoining unit cells, the displacements cancel each other, considering the fibre-like unit cells tightly constrained into the outer case. For the second set of equations (16), instead, homogeneous natural boundary conditions are taken throughout. In Fig. 12 (left), the profile of the radial stress along the radius (obtained by modelling only the matrix) is plotted for the first 300 nm, in the case of γ = 1, for different values of the length scale g . In Fig. 12 (right), instead, the profile of the same stress, obtained when both active site and matrix are modelled, is plotted for different values of the two length scales; for comparison purposes, one of the stress profiles obtained by modelling the matrix only is also plotted.
In particular, from Fig. 12 (left) it can be observed that the maximum compressive radial stress, occurring at the inner surface of the matrix, decreases for increasing values of the length scale g . The same happens  when the whole unit cell is modelled, even if the maximum compressive radial stress is not experienced at the interface between active site and matrix, but at a certain distance inside the matrix, depending on the value of the adopted length scales s and g , as shown in Fig. 12 (right). Furthermore, focusing the attention on the radial stress profile in the active site, it can be observed that the increase in the value of the length scale s , representative of the material of the active site, leads to a smoothing effect on the radial stress profile itself.
It is also worth noticing that different values of s and g can be interpreted as different materials used for the active site and the matrix, respectively.
The size effect analysis has been performed by studying the stress ratio defined as where σ ref is a reference stress, assumed to be equal to the maximum compressive radial stress obtained by applying classical elasticity, while σ g rr is the non-local radial stress. This stress ratio is an indicator for the spare capacity of the battery towards mechanical failure.
In Fig. 13, the stress ratios resulting from the two analysed problems are plotted against the ratio between the internal radius a and the length scale g (the length scales are taken as g = 8 nm in both cases and s = 4 nm). It can be observed that, while the local radial stress determined through classical elasticity does not show any size effect, gradient elasticity produces, in both cases, a stress ratio that increases for decreasing dimensions of the unit cell, while for increasing size of the unit cell it exhibits a horizontal asymptote, coincident to the solution of classical elasticity, which means that the maximum compressive non-local radial stress decreases for decreasing size of the unit cell, while it tends to the value determined according to classical elasticity for increasing dimensions. The only appreciable difference in Fig. 13 is that modelling the entire unit cell, the size effects become more pronounced, if compared with the case in which only the matrix is modelled.

Conclusions
In this paper, a (reduced-) gradient-enriched FE methodology has been presented for the solution of axisymmetric problems, subject to both axisymmetric and non-axisymmetric loads. A comprehensive discussion about the best Gauss integration rules to be used, error estimation and convergence behaviour is also included.
The methodology proposed in the present paper completes the one presented in [19] for plane stress/strain and three-dimensional problems, creating a unified FE methodology based on gradient elasticity, applicable to any sort of engineering problem.
The proposed methodology has been firstly applied to the problem of a hollow cylinder subject to a uniform internal pressure, and the results compared with the correspondent analytical solutions. This comparison has shown that, while for the displacements all the implemented elements produce results coincident with the analytical solution, for the stresses the numerical estimates differ from the analytical counterparts proposed by Gao and Park [34]. These differences can be ascribed to the fact that while the proposed methodology is a stressgradient elastic methodology, the one used by Gao and Park represents a displacement-based formulation, this leading in particular to different higher-order boundary conditions and different relationship amongst the stress components (while in the analytical solution proposed by Gao and Park the stress components are coupled, in the developed methodology they are uncoupled). Moreover, whereas Gao and Park [34] in their work use a full gradient approach, the methodology proposed in the present paper is based on the use of a reduced gradient.
The convergence behaviour of the proposed methodology has been checked in both the absence and the presence of singularities, for three different loading conditions. The results have shown that for problems without singularities, overall, the convergence rates of the numerical solutions are well in line with the correspondent theoretical predictions. In the presence of singularities, instead, the convergence rates are higher than the correspondent theoretical values, typical of traditional finite element methodologies. The cause of this faster convergence has been mainly attributed to two aspects: -the removal of singularities, characteristic of gradient elasticity theories; -the non-local stresses determined as primary variables, instead of secondary variables as in traditional finite element methodologies.
Furthermore, a guideline on optimal element size has been suggested, highlighting that relatively coarse meshes can produce sufficiently accurate solutions, with associated reduction in computational costs.
The ability of the proposed methodology to remove singularities from the stress fields has been shown by analysing three problems, consisting of a cylindrical bar characterised by a circumferential crack and subject to three different loading conditions.
The accuracy of the present methodology in the static and fatigue assessment of notched cylindrical bars has been also investigated for both brittle and ductile materials, covering a wide range of notch root radii as well. This study has shown that the proposed gradient-enriched FE methodology produces very accurate results, with errors ranging from −10% and +30%, which correspond to an error band equal to the usual ±20% but shifted 10% towards the conservative side. In addition, the proposed methodology allows the static and fatigue assessment of cracked/notched components by considering the relevant stress values directly at the tip of the stress riser (i.e. on the surface), avoiding the need to know a priori the failure location into the body, with evident advantages from an industrial point of view. Thus, it is possible to state that the proposed methodology can potentially become a powerful tool in the static and fatigue assessment of components presenting stress concentrators.
Finally, the proposed finite element methodology has been applied to evaluate the radial stress inside the glass/ceramic matrix surrounding the active sites in the electrodes of Li-ion batteries, showing the ability of gradient elasticity to describe size effects, or in other words the different mechanical behaviours of the matrix when the dimensions of the unit cell are scaled; in particular, it has been pointed out that scale effects become more significant for decreasing dimensions of the unit cell.
Thus, it is possible to conclude that, using gradient elasticity, with a limited additional computational cost with respect to classical elasticity, it is possible to take into account size effects, obtaining more accurate estimations of the stresses experienced by the matrix. Furthermore, this feature of gradient elasticity enables a more accurate choice of the material to use for the matrix (different materials mean different length scales) as well as the most appropriate dimensions of the unit cell.
In a similar way, gradient elasticity can be also employed to consider size effects for the benchmark problems presented in [9] to obtain input for the most favourable configurations (spheres vs. cylindrical fibres) and material selection for preventing fracture. Some interesting results are also expected by considering the parameters Δ and δ as being dependent on the concentration of the diffusing Li-ions. This will be discussed independently in a forthcoming paper.