A finite element implementation of the stress gradient theory

In this contribution, a ﬁnite element implementation of the stress gradient theory is proposed. The implementation relies on a reformulation of the governing set of partial diﬀerential equations in terms of one primary tensor-valued ﬁeld variable of third order, the so-called generalised displacement ﬁeld. Whereas the volumetric part of the generalised displacement ﬁeld is closely related to the classic displacement ﬁeld, the deviatoric part can be interpreted in terms of microdisplacements. The associated weak formulation moreover stipulates boundary conditions in terms of the normal projection of the generalised displacement ﬁeld or of the (complete) stress tensor. A detailed study of representative boundary value problems of stress gradient elasticity shows the applicability of the proposed formulation. In particular, the ﬁnite element implementation is validated based on the analytical solutions for a cylindrical bar under tension and torsion derived by means of Bessel functions. In both tension and torsion cases, a smaller is softer size eﬀect is evidenced in striking contrast to the corresponding strain gradient elasticity solutions.

Stress gradient elasticity is a generalised continuum theory that has been proposed very recently as a counterpart of Mindlin's well-known strain gradient elasticity model [25,14,26]. Stress gradient elasticity should neither be confused with the Aifantis gradient elasticity model which is a special case of isotropic strain gradient elasticity [29,13], nor with models based on Laplacians of stress and strain which have been recently recognised as special cases of micromorphic elasticity [10,16,7]. The essential feature of the stress gradient continuum is that each material point is endowed with a third rank tensor of additional degrees of freedom complementing its usual displacement. These new degrees of freedom, called microdisplacements in [14], are conjugate to the stress gradient tensor in the generalised work of internal forces. The microstructural origin of the microdisplacements was recently highlighted based on an original homogenisation approach in [18]. In contrast, the strain gradient model is solely based on the classic displacement degrees of freedom. The resolution of boundary value problems for stress gradient elastic bodies requires the definition of suitable boundary conditions. A debate commenced after the discovery of the stress gradient theory on the possibility of prescribing all stress components at a boundary which emerges in the stress gradient framework [27]. This is in contrast to classic Cauchy mechanics where the traction vector only can be controlled. Proper extensions of Korn's inequality show that the set of boundary conditions defined in [14] leads to a wellposed boundary value problem, with associated existence and uniqueness theorems [31]. They confirm the possibility of applying extended Neumann conditions in the form of fully prescribed stress components at a boundary.
Applications of stress gradient elasticity have been rather limited up to now. They include prediction of size effects in bending with a comparison with strain gradient and micromorphic approaches [28], and particle size effects in composite materials [33]. The latter reference contains extensions of the Eshelby and Hashin-Shtrikman homogenisation approaches to heterogeneous stress gradient media in a simplified case of isotropic elasticity.
The objective of the present work is to propose the first finite element implementation of stress gradient linearised elasticity theory. The choice of appropriate nodal degrees of freedom and matrix form of the variational formulation of the boundary value problem will be discussed. Validation examples with respect to analytical solutions are provided for simple tension and torsion in a simplified case of isotropic stress gradient elasticity. Free surface boundary conditions will be shown to play an essential role in the modification of the stress distributions compared to classic elasticity. The observed size effects with respect to structure size and the intrinsic length scale of the model will be discussed in detail.
In particular, this contribution is organised as follows: after a brief recapitulation of the fundamentals of the stress gradient theory in Section 2, a finite element formulation is discussed in detail in Section 3. By making use of the latter, representative boundary value problems are studied in Section 4. Moreover, analytical solutions are derived for validation purposes. The findings are summarised and concluding remarks are given in Section 5.

Notation
Let α, β, γ, δ, ζ, η denote arbitrary first order tensors, let the standard dyadic product be indicated by α ⊗ β and let the standard single tensor contraction be given by α · β. With these definitions at hand, double and triple tensor contractions are understood in the sense [α ⊗ β] : [γ ⊗ δ] = [α · γ] [β · δ] and [α ⊗ β ⊗ γ] ∴ [δ ⊗ ζ ⊗ η] = [α · δ] [β · ζ] [γ · η]. Moreover, the symmetric part of a second order tensor Ξ is given by with e • denoting Cartesian basis vectors. In analogy with second order tensors, the additive sphericaldeviatoric decomposition of third order tensors T that are symmetric with respect to the first two indices is introduced as In particular, the spherical part of a symmetric third order tensor in a three-dimensional setting is given by where the special dyadic product ⊗ was used, see also [14]. In addition to the second order identity tensor with δ ij denoting the Kronecker-delta, the fourth order symmetric identity tensor the sixth order symmetric identity tensor and the sixth order symmetric and deviatoric identity tensor with are introduced, see also [33]. Furthermore, (right-)gradient and (right-)divergence operators are denoted as ∇ (•) and ∇ · (•), respectively.

Stress gradient theory
The governing equations of the stress gradient theory as derived in [14] based on variational principles and on the method of virtual power will briefly be recapitulated in this section. At the outset of the theory, a generalised (volume specific) stress energy density function of the form w * (σ, R) is postulated which depends on the symmetric small deformation stress tensor σ and on the third order generalised stress tensor R. The associated minimisation problem of the (generalised) complementary energy functional on the domain B is moreover subjected to the extended set of static admissibility conditions with ρ denoting the mass density and f the (mass specific) body force density. Based on the symmetry of the stress tensor and on the observation that the volumetric part of the stress gradient is determined by (9a) since holds, the third order generalised stress tensor R is assumed to be symmetric on the first two indices and deviatoric. This entails the assumption that the deviatoric part of the stress gradient contributes to the stress energy density function w * (σ, R), in addition to the stress tensor. In order to derive the primal variational formulation of the stress gradient theory, the static admissibility conditions (9a) and (9b) are multiplied by the kinematic field variables δu and δΦ, respectively. In analogy with R, the third order tensor field δΦ is assumed to be symmetric and deviatoric. By integrating the ensuing equations over the domain B with boundary ∂B and with outward unit normal vector n, and by applying the divergence theorem one arrives at Equation (11) can be interpreted as a virtual work balance for the stress gradient theory, with the virtual work of internal forces P (i) (δu, δΦ), the virtual work of volume distributed forces P (v) (δu) and with the virtual work of contact forces P (c) (δu, δΦ). With regard to the virtual work of internal forces P (i) (δu, δΦ), energetic dualities are observed between the stresses σ and the generalised strain tensor and between the generalised stress tensor R and the kinematic variable Φ, which is henceforth referred to as the micro-displacement tensor. The latter observation motivates the introduction of a generalised strain energy density function w (e, Φ) such that The elasticity potentials w * (σ, R) and w (e, Φ) are related to each other by a Legendre(-Fenchel) transform. Moreover, the specific form of the virtual work of contact forces P (c) (δu, δΦ) stipulates Dirichlet-type boundary conditions in terms of the generalised displacements [δu ⊗ n] sym + δΦ · n and Neumann-type boundary conditions in terms of the complete stress tensor σ. A detailed discussion on the boundary conditions of the stress gradient theory is additionally provided in [14].

Finite element formulation
This section focuses on the finite element implementation of the stress gradient theory. In particular, a reformulation of the generalised virtual work equation (11) is discussed in Section 3.1 and the corresponding discretised weak form is derived in Section 3.2.

A different view on the virtual work equation
The derivation of the finite element implementation of the stress gradient theory is based on a reformulation of the generalised virtual work equation (11) in terms of one primary field variable Ψ , hereafter referred to as the generalised displacement field. The generalised displacement field is a tensor field of third order that is symmetric with respect to the first two indices and that was first introduced in the existence and uniqueness proofs of the stress gradient theory presented in [31]. It is constructed so that its deviatoric part is identical with the micro-displacement field, i.e.
whereas its spherical part is defined in terms of the (classic) displacement field, namely, The other way around, the displacement field can conveniently be reconstructed from the generalised displacement field via and it can be observed that holds, see also [31]. By making use of (12) and by invoking relation one moreover arrives at the representation of the generalised strain tensor in terms of the generalised displacements Finally, inserting (14), (16) and (17) into (11) yields the virtual work balance which is particularly useful for the finite element implementation, as it implies that either the (components) of the stress tensor or the corresponding normal projections of the generalised displacement field can independently be prescribed. In this sense, the two-field problem in terms of the displacement field u and the micro-displacement field Φ is reformulated in terms of one primary field quantity, i.e. in terms of the generalised displacement field Ψ .

Discretisation
The virtual work balance (20) serves as the basis for the derivation of the discrete weak form of the stress gradient theory. More specifically speaking, the generalised displacement field and the test function η are discretised by using a finite element approach with shape functions N • , nodal values Ψ • and η • , and with n en denoting the number of element nodes. By approximating the domain B with a finite number of elements n el , by denoting the respective element domains as B e and by making use of the assembly operator A , (20) and (21b) give rise to the discrete vector of internal forces to the discrete vector of volume distributed forces and to the discrete vector of contact forces Following standard procedure, the resulting system of (in general non-linear) equations is brought into a residual-type format r h and solved for the list of generalised nodal displacements Ψ , namely The linear part of the corresponding Taylor-series expansion ∆r h is given in terms of the generalised stiffness matrix (26) where dead-loads have been assumed and with • 13 T indicating a transposition with respect to the first and third index.

Representative simulation results
After the specification of the material model and of the constitutive equations in Section 4.1, representative simulation results are the focus of this section. In particular, a cylindrical bar under tension is studied by using the proposed finite element formulation in Section 4.2 and by means of analytical methods in Section 4.3. Breaking the rotational symmetry of the tensile problem, the focus is on a rectangular bar under tension in Section 4.4 before the torsion problem of a cylindrical bar is eventually studied in Section 4.5. Finite element results shown in this section are based on element-wise mean values.

Material model
Motivated by the theoretical developments of the stress gradient theory presented in [14,31,33], a symmetric positive definite quadratic form for the generalised stress energy density function is assumed with C and C denoting the compliance and the generalised compliance tensor, respectively. By making use of the Legendre(-Fenchel) transform, the corresponding generalised strain energy density function takes the form with the stiffness tensor E and the generalised stiffness tensor E defined as In the scope of this work, a classic form for the stiffness tensor E in terms of the Lamé-type constants λ and µ, namely is adopted and the generalised stiffness tensor is assumed to take the form with denoting the material length scale parameter. Moreover, the Lamé-type constants are related to the generalised Young's modulus and to the generalised Poisson's ratio via With the specific form of the strain energy density function defined by (28)-(31) at hand, the evaluation of (13) yields the specific form of the stress and generalised stress tensor Table 1 Material parameters used in the simulations of the cylindrical bar depicted in Figure 1.

Cylindrical bar under tension
In this section, the tensile problem of a cylindrical bar as depicted in Figure 1 is studied, with positions and tensor indices referring either to a Cartesian (x 1 , x 2 , x 3 ) or to a cylindrical coordinates (r, θ, z) description. At the left boundary of the bar (z = 0 mm), generalised clamping boundary conditions of the form are assumed for the generalised displacement field. It is noted that these do, in general, not imply vanishing (classic) displacements since holds. At the right boundary (z = 100 mm), generalised traction boundary conditions of the form with the overbar indicating a prescribed quantity, are applied which result in a total axial force of 660 kN. Moreover, the outer surface of the cylindrical bar (r = 10 mm) is assumed to be stress-free, i.e.
In order to study the size-dependent material response, the problem dimensions defined in Figure 1 are kept fixed and the material length scale parameter is varied. In particular, the set of material parameters provided in Table 1 is used.
The simulation results are based on a finite element discretisation with linear (eight-node) Lagrangian elements, and integrals are evaluated numerically by using a Gaussian quadrature scheme with eight sampling points. Taking into account the rotational symmetry of the problem, the predicted stress distribution in the cross section at z = 50 mm is presented in Figure 2 with respect to the natural (physical) basis system induced by the cylindrical coordinates. First, it is observed that the boundary condition (37) causes all coefficients of the stress tensor to approach zero at the boundaries. This is a significant difference compared to the classic Cauchy continuum approach where only the normal projection of the stress tensor can be prescribed and where a constant stress profile σ zz ≈ 2100 N mm −2 is expected when an isotropic, linear elastic material response is assumed and when boundary effects are neglected. In contrast, σ zz takes a parabolic profile in the present case, with the extreme value at r = 0 mm and the slope of the profile being functions of the material length scale parameter , see Figures 2(a) and 3(a). The in-plane coefficients of the stress tensor σ rr and σ zz , which would take zero values for a classic Cauchy continuum theory, are depicted in Figures 2(b), 2(c), 3(b) and 3(c). They show a strong dependence on the material length scale parameter and are observed to be approx. two orders of magnitude smaller than the axial stresses. Due to the vanishing stress boundary condition at r = 10 mm, it is moreover observed that the element-wise mean values of all coefficients of the stress tensor approach zero at the boundary, with the resolution being limited by the finite element discretisation. In addition, the radial displacement as a function of the material length scale parameter is provided in Figure 2(d). The corresponding profile becomes linear for smaller values of the internal length scale.

Comparison of simulation results and analytical solution in the case of vanishing Poisson's ratio
In order to derive an analytic solution for the boundary value problem depicted in Figure 1, a class of axisymmetric stress states is considered in the form where the stress components σ r = σ r r , σ θ = σ θθ and σ z = σ zz are unknown functions of r in the cylindrical coordinate system (r, θ, z). Subject to the assumptions of quasi-statics and vanishing body forces, see (9b) and (10), the generalised stress tensor is computed as with primes denoting derivatives with respect to r. By additionally invoking the static balance equation for the active stress components, i.e.
the generalised stress tensor can be rewritten as +σ θ e θ ⊗ e θ ⊗ e r + σ z e z ⊗ e z ⊗ e r .
For the computation of the third order tensor Φ of microdisplacements, it is necessary to evaluate the divergence of the generalised stress tensor where the equilibrium equation (40) has again been taken into account. The displacement vector and its gradient are taken in the form Under simple tensile/compression loading conditions, the displacement gradient contribution ∂uz ∂z =ε is prescribed. The generalised strain tensor can now be computed by making use of (12), (29)-(33) and by assuming  that and µ are constant in space, i.e.
Combining (42)-(44) leads to the differential system This system must be complemented by the balance relation (40) in order to be solved for the four unknown functions u r , σ r , σ θ , σ z . In accordance with Figure 1, vanishing stress boundary conditions are prescribed at the outer boundary of the cylindrical bar with radius r m , i.e.
It can be shown that the functions σ r = σ θ = 0 are no solutions of the system, with the considered boundary conditions, except in the very special case ν = 0. This is in contrast to simple tension in a cylindrical bar for the Cauchy continuum.
In the particular case ν = 0 the system simplifies and it is found that u r = 0, that σ r = σ θ = 0 and that the axial stress is solution of the ordinary differential equation A particular solution is the classic one σ z (r) = Eε, which, however, does not fulfil the boundary condition of vanishing stresses at r = r m . With the scaling x = r/˜ and with the notation y(x) = σ z (˜ x) being adopted, the homogeneous equation to solve is given by where primes indicate derivatives with respect to x. This is a special case of the modified Bessel's equation with α = 0, see [1]. The solution is given by a linear combination of modified Bessel functions, also denoted as hyperbolic Bessel functions, of the first and second kind: I 0 (x) and K 0 (x). The function K 0 is singular at x = 0 and cannot appear in the solution of the considered boundary value problem. Thus, the solution is of the form where σ 0 is the stress value at the centre that can be related to the stress value Eε = C 2 . Finally, this gives (52) Figure 4 shows the excellent agreement between the analytical and finite element solutions for three different length scales and vanishing Poisson's ratio. Regarding the analytical solution of the tensile problem given by (51) and (52), it is moreover observed that the solution converges to the one of a classic Cauchy continuum for → 0 + , except for an increasingly thin boundary layer that occurs due to the vanishing stress boundary condition at the outer surface. In addition to the latter results, the convergence of the finite element results towards the analytical solution upon mesh refinement is studied in Appendix A.

Rectangular bar under tension
Breaking the rotational symmetry of the boundary value problem discussed in Section 4.2 and Section 4.3, the tensile problem of a rectangular bar with a square cross section sketched in Figure 5 is studied in this section. In accordance with the previous sections, generalised clamping boundary conditions are assumed at the left boundary (z = 0 mm), and generalised traction boundary conditions that result in a total axial force of 840 kN are applied at the right boundary (z = 100 mm). At the remaining boundaries, homogeneous generalised Neumann boundary conditions are assumed. The material parameters are chosen in accordance with Table 1, the bar is discretised by using linear (eight-node) Lagrangian elements, and a Gaussian quadrature scheme with eight sampling points is used for numerical integration.
The finite element results are presented with respect to the Cartesian basis specified in Figure 5 and are evaluated for the cross section at z = 50 mm to reduce the influence of boundary effects. The simulation results are summarised in Figure 6 and differ significantly from the homogeneous, uniaxial stress state expected in the case of a classic Cauchy continuum with an isotropic, linearelastic material model: 1) A strongly inhomogeneous axial stress field σ 33 is observed in Figure 6(a-c). The axial stress takes its maximum value in the centre of the cross section and approaches zero at the boundaries. In accordance with the observations for the cylindrical bar in Section 4.2, the slope and the maximum value of the σ 33 profile depend on the material length scale parameter which penalises stress gradients in an energetic manner. 2) Significant in-plane stresses are observed in Figure 6(d-i). 3) The complex stress state in the cross-section is associated with a complex deformation of the cross-section which, in particular, does not remain square. The latter deformation mode becomes more pronounced with increasing values of the material length scale parameter . Moreover, no warping of the cross section's surface with unit normal vector e 3 is observed.
In contrast, the homogeneous uniaxial stress state that is expected (at a distance from the clamped boundary) in the case of a classic Cauchy continuum with an isotropic linear-elastic material model can be recovered by taking into account a different set of boundary conditions. More specifically speaking, the bar is again assumed to be clamped (in a generalised sense) at the left boundary (x 3 = 0 mm) and homogeneous traction boundary conditions are assumed at the right boundary (x 3 = 100 mm). At the remaining boundaries, generalised stress boundary conditions that are consistent with the expected uniaxial stress state are prescribed, i.e. The simulation results that were calculated with the set of boundary conditions (53)-(55), cf. Figure 5, and with = 1.0 mm are provided in Figure 7. For an applied axial load of 840 kN, the uniaxial stress state at a distance from the clamped boundary is given by as indicated by light grey colour. In addition, an inhomogeneous stress distribution due to the clamping conditions is observed in the vicinity of the left boundary.

Torsion test
The torsion problem of a cylindrical bar as depicted in Figure 8 is analysed in this section. The bar is assumed to be clamped (in a generalised sense) at the left boundary. The remaining boundary, except for the right end of the cylindrical bar where a torque is induced by means tangential forces, is assumed to be stress free. In particular, traction boundary conditions at the surface with outward unit normal vector n = e r , in the vicinity of the right boundary (z = 100 mm), of the form are applied that result in a torque of approx. 1100 N m with respect to the e z -axis. Moreover, by making use of the relation between the Cartesian and polar basis e r = cos (θ) e 1 + sin (θ) e 2 , (58a) an isotropic, linear elastic material model. In contrast, a parabolic profile that approaches zero at the centre (r = 0 mm) and at the outer radius (r = 10 mm) is observed in Figure 9, with the slope and the maximum shear stress value depending on the material length scale . In addition, the shear stress distribution in the cross section is provided in Figure 10 for the three different values of the length scale parameter studied.

Comparison with the analytic solution in torsion
The stress tensor takes the form where τ = σ θz is the unknown function to be determined. The static balance equations, in the absence of body forces, require that τ is a function of the sole variable r. Under these conditions, the stress gradient is deviatoric and one finds that where τ denotes the derivative of the stress component with respect to r. Its divergence is then computed as On the other hand, the displacement vector is described by a single function u θ (r, z) in the form from which the infinitesimal strain tensor is computed as Assuming that and µ are constant in space, the generalised strain tensor can now be expressed as In the analysed torsion case, generalised Hooke's law (33a) simplifies to By taking into account (61), (66) and (67), the differential system of two equations ∂u θ ∂r − u θ r = 0 (68a) is thus derived. It follows from (68a) that u θ (r, z) = α(z) r where α is a function of the sole variable z. Furthermore, (68b) shows that this function reduces to α(z) = α 0 z. Equation (68b) can now be written in the form A particular solution of this equation is the well-known linear field from classic elasticity τ (r) = µ α 0 r .
The remaining homogeneous ordinary differential equation to be solved is then where definition (48) was used. The change of variable r =˜ x and the definition y(x) = τ (˜ x) lead to the ordinary differential equation where primes indicate derivatives with respect to x. According to [1], the solutions are given by modified Bessel functions of order one, i.e.
The modified Bessel function of the second kind K 1 is singular at x = 0 and is therefore not acceptable for the considered bar. The solution of the torsion problem thus takes the form τ (r) = µ α 0 r + C 1 I 1 (r/˜ ) = µ α 0 r m r r m − I 1 (r/˜ ) where the constant C 1 was determined from the condition τ (r = r m ) = 0. The solution is found to be such that τ (r = 0) = 0, with the vanishing shear stress at the tube's axis being a consequence of the specific form of the analytical solution and not imposed to the governing differential equation. In accordance with the tensile problem discussed in Section 4.3, it is moreover observed that the analytical solution (74) converges towards the one of a classic Cauchy continuum in the limit˜ → 0 + , except for an increasingly thin boundary layer that occurs due to the vanishing stress boundary condition at the outer surface. In addition, the convergence of the finite element results towards the analytical solution upon mesh refinement is exemplarily shown for = 1.0 mm in Appendix A.
The next step consists in evaluating the torque resulting from the computed shear stress function. The only non-vanishing component of this vector is M e z and given by with x m = r m /˜ . The last integral can be evaluated exactly 1 by means of two successive integrations by parts as xm 0 The torque can now be related to the torsion angle per unit length, α 0 , as The torque value was prescribed in the finite element analysis of Section 4.5. The value α 0 can be computed from (77) and completely determines the shear stress field (74). A comparison of the analytical and numerical shear stress profiles is plotted in Figure 11 for three values of the intrinsic length scale. Excellent agreement is obtained, which shows that the usual linear shear stress profile is replaced by a non-monotonic distribution with a boundary layer close to the free surface. 1 The following identities for Bessel functions have been used: I 0 (x) = I 1 (x), I 1 (x) = I 0 (x) − I 1 (x)/x .

Smaller is softer
The torsion stiffness J can be defined from relation (77). It is a function of the wire radius and of the characteristic length of the stress gradient elastic continuum. Specifically speaking, holds, with J 0 denoting the usual torsion stiffness of a cylindrical bar A tensile stiffness of the cylindrical bar can also be defined from the solution found in Section 4.3. The stress distribution was given by (51) and (52). Alternatively, it can be expressed in terms of the stress value Eε, whereε is the prescribed axial displacement gradient, namely The tensile stiffness of the bar is defined as the ratio of the applied force F and the prescribed relative displacement δ = Lε, with L denoting the length of the bar. In the case of a Cauchy continuum, it is well-known that the reference stiffness is For the stress gradient continuum, the total force is computed as It follows that with the relative tensile stiffness defined as The tensile and torsion relative stiffnesses are plotted in Figure 12. The curves clearly reveal a smaller is softer tendency, meaning that the stiffness decreases when the intrinsic length scale increases. It is particularly remarkable that the bar stiffness tends towards zero when the length scale increases or when the bar radius decreases. This effect stems from the existence of a boundary layer with vanishing stress whose thickness becomes dominant for small bar radii or large intrinsic length scales. The limit case of torsion for large radii or small internal length scales is also interesting. The usual stiffness J 0 is retrieved in this limit case in spite of the presence of a stress gradient. The reason is that, at the limit, the shear stress is a linear function implying that the divergence of the stress gradient and of microdisplacements (see (63)) vanishes, leaving the usual strain tensor unaffected.
A comparison can be drawn with solutions according to existing strain gradient elasticity models, including Mindlin's original theory [25] or Aifantis simplified model [13]. No size effect is expected in tension in strain gradient elasticity at least for long enough bars fulfilling Saint-Venant's conditions. The softening effect obtained in torsion according to stress gradient elasticity is in striking contrast to the predictions of strain gradient elasticity which are notoriously associated with a smaller is stiffer effect, see [6]. Solutions of the torsion problem of strain gradient elastic bars have been provided in [2,32,20]. A torsion stiffening effect is predicted by strain gradient elasticity when the bar radius has the same order of magnitude as the intrinsic length scale and below. The linear stress profile is still valid in a strain gradient elastic circular cylindrical bar because it fulfils all required higher order boundary conditions. This is in contrast to the stress gradient solution worked out in the previous section. However, the torsion stiffness is increased by a contribution of the higher order elasticity moduli [19,20], at least if the bar is long enough to fulfil the Saint-Venant conditions, see [23] which focuses on special boundary conditions at the bar's ends. This relative stiffness enhancement is proportional to 1 + a x −2 m , where x m is the ratio of the bar radius divided by one strain gradient elastic length scale, and where a is a factor depending on the specific strain gradient model considered. This implies that the relative stiffness approaches infinity for vanishingly small radii or very large intrinsic length scales. It is recalled that the stress gradient elasticity theory predicts a vanishing relative stiffness under these conditions. Finite element solutions of the torsion problem of strain gradient elastic bars were presented in [5], confirming the absence of warping for circular bars and computing the warping of elliptical sections.

Closure
In this contribution, a finite element implementation of the stress gradient theory was proposed and studied in detail. The implementation relied on a reformulation of the governing set of field equations (of the stress gradient theory) in terms of one primary, tensor-valued field variable of third order. Since the latter field variable was based on the classic displacement field and on a microdisplacement field-type quantity, it was referred to as a generalised displacement field. It was then shown that a reformulation of the weak form of the balance equations in terms of the generalised displacement field is particularly suitable for a finite element-based implementation, especial with regard to (the implementation of) boundary conditions. More specifically speaking, a peculiarity of the stress gradient theory as opposed to the classic Cauchy continuum was that the complete stress tensor could be prescribed at the boundaries. With the finite element implementation at hand, representative boundary value problems were studied in detail. In particular, a cylindrical bar under tension was analysed for which an analytical solution by means of Bessel functions could be derived for validation purposes in the case of a circular cross section. Breaking the rotational symmetry of the bar, the study of a rectangular bar under tension revealed a significantly different deformation behaviour as compared to a classic Cauchy continuum. These results strongly differ from strain gradient elasticity which does not predict any size effect for a bar in tension, at least in the sense of Saint-Venant. Finally, the (unusual) shear stress distribution in a cylindrical bar under torsion was studied by means of numerical and analytical methods. Note that the analytical solution of the simple tension problem was provided explicitly for a circular bar in the special case of vanishing Poisson's ratio. In contrast, the solution presented for torsion loading in isotropic stress gradient elasticity is general.
The stress gradient continuum represents a generalisation of a classic Cauchy continuum. Amongst others, it allows for generalised boundary conditions in terms of the complete stress tensor and naturally accounts for a material length scale. The simulation examples provided in this work illustrate a remarkable feature of stress gradient elasticity which predicts that smaller is softer in the presence of stress free boundaries. This is in contrast to smaller is stiffer effects according to strain gradient elasticity. Free boundary layer effects can be accounted for by means of second strain gradient elasticity and essentially lead to an increase of the tensile stiffness of a bar [9,22]. Smaller is softer effects have been reported in materials science for the apparent Young's modulus of nanowires in [36] in relation to free surface effects, see [9] and references quoted therein. This trend was also observed in the plasticity of bulk metallic glasses according to [17]. It may also be relevant in the context of boundary layer effects in elastic composite materials, see the bending-gradient theory for thick heterogeneous plates discussed in [30]. Additionally, both smaller is stiffer as well as smaller is softer effects have been reported in bending experiments on human cortical bone samples [35,8]. These seemingly contradictory results prompted the analyses presented in [34], which suggests by employing simple analytical laminate beam theories and finite elementbased studies of two-dimensional beam-type samples with periodic heterogeneities that the latter experimental findings can be explained by the particular properties of the material surface. From a modeling point of view, recent research efforts moreover focus on the combination of gradient and nonlocal elasticity theories to account for both smaller is softer as well as smaller is stiffer effects that may, possibly, be related to different material length scales, [3]. Applications of nonlocal strain gradient theories to simulate beam-and shelltype nanoscale structures are discussed in, for example, [11,24], and it remains to be shown how these theories are related to the stress gradient theory that serves as the basis for the present contribution.
In addition, the simplified isotropic elasticity law used in the present work involving a single additional intrinsic length scale should be extended to more general sixth order tensors of stress gradient elasticity, including isotropy and anisotropy [14,4]. Furthermore, more complex geometries and loadings will be investigated in future works so as to reveal the full potential of the stress gradient theory. Another challenging issue to be considered is the introduction of nonlinear aspects of the stress gradient theory into the proposed finite element approach, regarding either the finite deformation framework established in [15] or plasticity theory [12]. The latter would allow for a comparison between stress gradient plasticity and strain gradient plasticity predictions for a wealth of physical situations which were explored in the past only in the case of strain gradient approaches [21,6].

A Convergence studies
The convergence upon mesh refinement of the finite element results discussed in Section 4.3 and Section 4.5 is briefly studied in this section. In particular, the convergence of the element-wise mean values of the axial stress σ zz for the tensile problem and of the in-plane shear stress σ θz for the torsion problem against the analytical solution is exemplarily shown for = 1.0 mm in Figure 13. The simulation results for the different finite element meshes are additionally provided in Figure 14 and Figure 15 in terms of surface plots.

Declarations
Funding: Open Access funding provided by Projekt DEAL.

Conflicts of interest/Competing interests:
The authors declare that they have no conflict of interest.
Availability of data and material: Research data and material will be provided by the authors upon request.