3D strain gradient elasticity : Variational formulations, isogeometric analysis and model

This article investigates the theoretical and numerical analysis as well as applications of the three-dimensional theory of first strain gradient elasticity. The corresponding continuous and discrete variational formulations are established with error estimates stemming from continuity and coercivity within a Sobolev space framework. An implementation of the corresponding isogeometric Ritz–Galerkin method is provided within the open-source software package GeoPDEs. A thorough numerical convergence analysis is accomplished for confirming the theoretical error estimates and for verifying the software implementation. Lastly, a set of model comparisons is presented for revealing and demonstrating some essential model peculiarities: (1) the 1D Timoshenko beam model is essentially closer to the 3D model than the corresponding Euler–Bernoulli beam model; (2) the 3D model and the 1D beam models agree on the strong size effect typical for microstructural and microarchitectural beam structures; (3) stress singularities of reentrant corners disappear in strain gradient elasticity. The computational homogenization methodologies applied in the examples for microarchitectural beams are shown to possess disadvantages that future research should focus on.


Introduction
During the development of modern science and technology, continuum mechanics has turned out to be a powerful physico-mathematical framework of theories and models describing the mechanics and thermodynamics of solids and fluids.According to its underlying homogenization assumptions, it deals with the motion and stress of continuous media by ignoring the particle nature of matters.Regardless of ignoring the atomistic or molecular particles, even larger-scale microstructures, the heart of continuum mechanics are the five fundamental principles of physics: conservation of mass, balance of linear momentum, balance of angular momentum, energy conservation, and the second law of thermodynamics on entropy.These principles enable to uncover mathematical relationships between displacements, stresses, strains, velocities, and temperatures for the study of the mechanical and thermodynamical behaviour of continuous media.
The classical continuum mechanics, or Cauchy's continuum mechanics, have been successfully used in a wide range of applications of civil, mechanical, materials and geoengineering, in particular.Although Cauchy's continuum mechanics was intended to the macro-scale world, it has been applied even for nano-scale structures, e.g., for graphene in [1].However, due to the decreasing size of the Representative Volume Element (RVE), versus the micro-and nano-scale building blocks (i.e., grains, unit cells or atoms) of the material of a small-scale structure, the classical continuum mechanics hypotheses should not be valid anymore.At these smaller scales, the internal degrees of freedom of the microstructure, ignored in the classical continuum mechanics, may have a significant role in the thermomechanical behaviour of the structure.
The so-called size dependency related to material microstructure, or structural microarchitecture, has been observed in many experiments.For instance, the bending rigidity of micro-scale structures has shown to increase significantly due to size effects [2][3][4].Regarding dynamics, theoretical and experimental studies have shown that material microstructure may have fundamental effects on elastic wave propagation [5][6][7].Stress singularities in elasticity, in the presence of an applied point or line load or at a crack tip [8][9][10], are a noticeable class of nonphysical concepts appearing in the classical continuum mechanics; in this article we investigate the singularity at a reentrant corner as an example.Generalized continuum mechanics has been an attempt to capture these mechanical behaviours by extending Cauchy's continuum mechanics without losing its generality, applicability and relative simplicity.The Cosserat brothers are usually considered as the first pioneers of such theories in the beginning of the 19th century, although a trace can be followed even back to Piola [11].According to the Cosserats theory, material points are infinitesimal and rigid, and each point has rotational degrees of freedom in addition to the standard displacement degrees of freedom.In this type of formulation, moment stresses appear supplement to the classical stresses; hence, this has lead to the term couple stress theory [12,13].After a recess period in the development of continuum mechanics, in the 1960s, a new renaissance dawned.In a couple of decades, scientists introduced strain gradient theories [14,15], micromorphic theories [16,17], continua with surface energy [18,19], multipolar theories [20], and non-local theories [21].
The core of this article relies on the first strain gradient elasticity theory proposed by Mindlin [14].In the classical continuum mechanics, the deformation energy comprises only the first gradient of displacement, strain.However, in strain gradient elasticity, as the name indicates, higher displacement gradients, strain gradients, are incorporated to bring more richness to the continuum description.Strain gradient elasticity -often called second gradient elasticity as it incorporates the first and second gradient of displacement -has already been adopted for a wide variety of applications to capture size and non-local effects in small-scale structures, or microstructured materials, in general.In strain gradient elasticity, the strain tensor (being a second-order tensor) and its gradient (being a third-order tensor) are related to the (second-order) stress tensor and the (third-order) double stress tensor, respectively.In a material with a centrosymmetric microstructure, the corresponding fourth-and sixth-order elastic stiffness tensors postulate the constitutive equations.Cowin and Mehrabadi [22] and Auffray et al. [23,24] have profoundly studied the properties of these elastic stiffness tensors, leading to symmetry group classifications.
One of the exotic applications of strain gradient elasticity is the homogenization of periodic elastic media [25][26][27][28][29][30][31][32][33][34][35].The latest advances in additive manufacturing have paved the road for the design and fabrication of lattice structures at different length scales.3D printing technology offers the potential to manipulate and enhance materials and structures at different scales.By modifying the unit cell of a periodic lattice, the corresponding lattice structure attains extraordinary mechanical properties such as high stiffness, strength and toughness or vibration band gaps, while reducing the material density [36][37][38][39]31].Furthermore, it also provides possibilities for energy absorption and porosity enhancement of structures [40].Considering all of the geometrical details of a lattice unit cell, and consequently of the whole body in consideration, often results in large computational models.Theoretical or computational homogenization techniques, however, reduce the model size significantly.Accordingly, besides the standard homogenization methods, there has been many methodologies specific for the generalized homogenization of microarchitectural structures under development [41,25,42,[26][27][28][29][30][32][33][34][35] -some developed for the three-dimensional framework (or the two-dimensional framework of plane elasticity), others dedicated to the structural dimension reduction models such as bars, beams, plates and shells.In the present contribution, for obtaining the constitutive parameters for the full (3D) strain gradient elasticity tensors we apply the method of [43,44] and compare the resulting 3D model to the corresponding 1D beam models and a related simplified 3D model.
Literature on strain gradient elasticity is ample when it comes to the corresponding dimensionally reduced structural models, as can be seen in the reviews for one-dimensional beams and two-dimensional plates provided in [45][46][47].Similarly, there are plenty of studies available for the problems of plane elasticity, for both theoretical (see, e.g., [48,49]) and numerical studies (e.g., [50][51][52]), whereas for shells only a few articles can be found; see the list in [53].After all, although the original work of Mindlin was written for the general case of three-dimensional bodies, only very few later studies have focused on this framework [10,[54][55][56].Connected to the fundamental structural models, one target of the three-dimensional framework of the present article is to confirm the correctness of the so-called stiffening effect crucial for applications and characteristic to the dimension reduction models of several generalized continuum theories, as shown in the model comparison for beams and plates [57,46].With respect to this original feature -stemming from the insertion of the standard kinematical assumptions into the full three-dimensional formulation -different models can be found in literature [45], mainly due to the fact that the stiffening feature is missing from some of the first models (e.g., [58]), but the validity of the effect has still been questioned quite recently [59,60].Regarding the stiffening effect, the following has already been addressed in the literature.First, if the gradient terms of the thickness direction are included in the mathematical derivation of a dimension reduction model of beams, plates or shells, the model exhibits a stiffening effect, analogous to the different strain gradient models and modified couple stress models [57,46], whereas including only the gradient terms in the direction of the mid-axis or mid-surface essentially removes the stiffening effect.Second, some planar lattice configurations have been shown to follow the dimension reduction models incorporating the stiffening terms [29,30].The present contribution, by providing first a conforming and verified discretization for the 3D model, focuses on comparing the 1D beam models and the 3D model: the comparison reveals that the 1D (Euler-Bernoulli and Timoshenko) beam models are, indeed, in accordance with the three-dimensional model if and only if the gradient terms of the thickness direction are included in the dimension reduction -which shows that the stiffening effect is not a peculiarity of the structural models, or a result of erroneous model derivations, but shared by the 3D model as well.
One of the main obstacles in the development and applicability of the strain gradient elasticity is solving the resulting high-order differential equations.The H 2 -formulation of first strain gradient elasticity problems requires a Ritz-Galerkin method of C 1 -continuous shape functions for conformity, whereas most of the existing commercial or academic FEM-packages provide solely continuity of order C 0 -particularly for three-dimensional problems in which the complexity of the C 1 -continuous discrete system escalates: either the number of degrees of freedom or the interaction between neighbouring nodes starts to increase rapidly.As a matter of fact, in the framework of C 0 -continuity, one can formulate the problem by applying mixed formulations [61], Lagrange multipliers [62], penalty parameters [63], discontinuous Galerkin formulation [64], or other nonconforming methods such as the Morley triangle [65][66][67] and its generalization to other dimensions [68].However, the main disadvantages of the abovementioned options are the relatively low accuracy and slow convergence.Within the C 1 -continuous finite element method, the most commonly known elements -but for two-dimensional problems -are perhaps the triangular Argyris [69], Bell [70], Hsieh-Clough-Tocher triangle [71], and Powell-Sabin [72] elements.Alternatively, in the context of isogeometric analysis (IGA), H 2 -conformity can be achieved by the standard C 1continuous NURBS-shape functions also in the context of three dimensions.For the two-dimensional problems of strain gradient elasticity, investigations on the performance and applications of the different C 1 -continuous methods can be found in [50] (FEM) and with convergence analyses in the proper norms in [51] (FEM) [52,73] (IGA).For the corresponding three-dimensional problems, a standard C 0 -continuous tetrahedral finite element implementation utilizing Lagrange multipliers has been investigated in [10] and a C 1 -continuous hexahedral finite element has been proposed and applied in [54][55][56].Convergence analysis has a small role in these studies and the convergence orders remain open, in particular.
On these grounds, the targeted novelties of this article on the three-dimensional version of strain gradient elasticity and isogeometric analysis are the following; (1)-(4) from the point of view of numerics, (5)-( 6) from the point of view of mechanics: (1) variational formulations within a Sobolev space framework of order H 2 (Ω ), Ω ∈ R 3 ; (2) solvability via continuity and coercivity; (3) error estimates for conforming Ritz-Galerkin formulations; (4) an isogeometric implementation with a verification and confirmation for the theoretical error estimates; (5) model comparisons between the 3D solid model and the corresponding 1D beam models approving the size effects peculiar for strain gradient models; (6) demonstrations for model peculiarities within the three-dimensional framework, including a computational homogenization for a cellular structure and an example about the removal of a stress singularity.Our implementation has been accomplished into the GeoPDEs-package written (in Octave) for solving partial differential equations with isogeometric analysis [74].
The outlines of this article are the following: In Section 2, we briefly recall the theory of first strain gradient elasticity.Section 3 is devoted to the continuous and discrete variational formulations.In Section 4, we concisely recall the basics of isogeometric analysis and address the position of our implementation within the structure of the GeoPDEs software.In Section 5, the implementation is verified, the theoretical error estimates confirmed and a few demonstrative applications of the three-dimensional formulation are presented.Finally, we conclude this contribution in the last section.

Three-dimensional theory of strain gradient elasticity
In this section, we recall the linear theory of first strain gradient elasticity for three-dimensional continuum bodies.The related formulations are presented with the so-called index notation.
Let us consider Hamilton's principle for an independent variation δu i in the form where T and W are, respectively, kinetic and strain energy density of a body Ω ∈ R with ρ donating the mass of the material per unit volume, J i jkl standing for the inertial size effect parameters, and dot () over displacement u i is the partial time derivative.The strain energy density can be written as where ε i j denotes the standard linear strain tensor.For a rotation-independent strain measurement, the Green strain tensor is defined as with u i, j representing the gradient of the displacement.The nonlinear part of the strain can be neglected in the case of small displacements giving the linear strain tensor and its gradient In Eq. ( 3), σ i, j and µ i jk stand for the stress and double stress tensor, accordingly.They are the energy conjugates of the strain and the gradient of strain: where C i jkl and A i jklmn are the fourth-order elastic stiffness tensor and the sixth-order elastic gradient stiffness tensor, respectively.In the triclinic anisotropic case, they contain 21 and 171 independent parameters, accordingly.However, for an isotropic case, there will be only two independent parameters for the elastic stiffness tensor and five parameters for the gradient-elastic stiffness tensor; [22][23][24] cover all the symmetry classes of the classical and strain gradient elasticity.Eventually, we formulate the variational form of each energy expression for Eq. ( 1).The variation of the total kinetic energy is given as By substituting Eq. ( 6) into Eq.( 3), the variation of the total strain energy expression takes the form and the variation of the work done by external forces is given as where f i denotes a body force vector, T i represents a traction force on the boundary surface Γ , φ i stands for a body double force tensor, and τ i denotes a double traction force.
One can obtain the equation of motion in the strong form by dropping the time integration in Eq. ( 1) and by applying some standard calculus of variations [14]: where D p is the surface gradient operator defined by

Variational formulations and error estimates
In this section, the three-dimensional problems of strain gradient elasticity are formulated in a variational form within a Sobolev space setting.Within this context, the solvability of the problem is proved by the continuity and coercivity of the corresponding bilinear form and the continuity of the load functional.These results provide, in particular, a basis for many variational discretization methods such as the isogeometric Ritz-Galerkin method adopted in the present study.
According to the previous section, let us consider a solid body Ω ∈ R 3 subjected to a static body force f i , a traction force T i , a double body force φ i j and a double traction force τ i .In what follows, we will use the notation H s (Ω ) for a real Sobolev space of order s, consisting of square-integrable real-valued functions defined on Ω with square-integrable weak derivatives up to order s.The corresponding Sobolev norm and seminorm are denoted by ∥ • ∥ s and | • | s , respectively.
According to ( 7)-( 9), the variational formulation of the static case of problem ( 10)-( 16) reads as follows: where the bilinear form a : U × V → R, and the load functional l : V → R, respectively, are defined as The trial function set of the kinematically admissible displacements consists of functions satisfying the essential boundary conditions with the given classical Dirichlet data for the displacements on Γ D c (see (11)) and the given non-classical Dirichlet data for the displacement derivatives on Γ D g (see (12)), while the test function space V consists of functions satisfying the corresponding homogeneous Dirichlet boundary conditions.It is assumed, for notational simplicity, that all components are fixed on Γ D c and Γ D g .
The energy norm of the problem induced by the bilinear form is equivalent to the H 2 -norm whenever U = V , which can be seen in the Appendix of [52] providing proofs for the plane elasticity versions of the propositions below.In addition, the symmetry of the bilinear form is clearly guaranteed: For positive definite material tensors C and A, the continuity and coercivity of the bilinear form a are stated here for a domain having a fully clamped boundary Γ = Γ D c ∩ Γ D g =: Γ D (chosen in Section 5 for the numerical verification as well), for simplicity.
Proposition 2. Let us assume that Ω is convex with a Lipschitz continuous boundary Γ = Γ D , and U = V .For positive definite tensors C and A, there exists a positive constant α such that It should be noticed that these propositions have been proved in the Appendix of [52] for the two-dimensional version of strain gradient elasticity under the following simplifying condition (sometimes called the condition of weak non-locality): A i jklmn = g 2 δ mn C i jkl in which C i jkl refer to the constitutive parameters of classical elasticity, defined for isotropic materials by the Lame parameters (or Young's modulus and Poisson's ratio denoted, respectively, by E and ν in what follows), whereas g stands for a length scale parameter.Under this condition in case of isotropic materials, , where c comes from the Poincaré-Friedrichs inequality.The proofs are, however, naturally extendable to the more general cases of homogeneous materials with positive definite material tensors A, including the present framework of three-dimensional strain gradient elasticity.It is worth emphasizing that the expression of the constant α above indicates that, in any case, the constant of the positivity of the higher-order constitutive tensor is crucial particularly for coercivity.
Finally, according to Riesz Representation Theorem the problem has a unique solution, as long as the load functional l is assumed to be continuous -as in the case of a standard loading f ∈ [L 2 (Ω )] 3 chosen below for the numerical confirmation in Section 5: Proposition 3. Let us assume that Ω is convex with a Lipschitz continuous boundary Γ = Γ D , C and A are positive definite and U = V .For a continuous load functional, Problem 3 has a unique solution in V .
A conforming Ritz-Galerkin method for the problem is formulated in a standard way: In isogeometric analysis, the trial function space for a displacement field u h i is defined by where it holds that S h ⊂ H 2 (Ω ) providing C s−1 (Ω )-continuity within a patch Ω .Approximation spaces are obtained by setting With the conformity of the method, the continuity and coercivity of the continuous problem are inherited by the discrete method implying error estimates which can be proved in a standard way by imitating the steps for proving Cea's lemma [76]: This qualitative error estimate gives, with C 1 -continuous discretizations, more quantitative estimates for the error by utilizing the approximation properties of NURBS or other spline functions [77,78] (as the corresponding classical ones for polynomes [79]): Corollary.With the assumptions of Proposition 4, and by assuming that the exact solution of the Problem smooth enough, i.e., u ∈ [H p+1 (Ω )] 3 , the following error estimate holds: where c denotes an interpolation constant stemming from the choice for the discrete space.
The numerical results of Section 5 for isogeometric methods with C p−1 -continuous NURBS basis functions of order p ≥ 2 confirm these theoretical results providing the convergence order O(h p−1 ).For lower order norms, additional powers of h are obtained -as confirmed by the numerical results of Section 5.

Isogeometric analysis within GeoPDEs
The theory of the isogeometric analysis was introduced by Hughes et al. [80] to fill the gap between computeraided design and analysis.In a typical process of design, designers create CAD (Computer-Aided Design) files for geometries and thenceforth the files are translated into a corresponding FEA (Finite Element Analysis) model by reconstructing and remeshing the geometries via polynomial approximations.The heart of isogeometric analysis is to adopt CAD functions, i.e., splines such as non-uniform rational B-splines (NURBS), for representing the FEA geometry (now faithful to the CAD geometry) and further approximating the field variables of the engineering analysis in an isoparametric way.
Most of the CAD geometries are represented with NURBS.For instance, a single-patch solid geometry is defined by x : [0, 1] 3  → Ω with C i, j,k representing the vector of control point coordinates.The 3D tensor product of the NURBS shape functions R p,q,r i, j,k is given by where w i, j,k stands for the weight coefficients and N p i , M q i , L r i are the B-spline shape functions corresponding to each parametric direction ξ , η, ζ .The 1D B-spline basis functions are constructed from a given open knot vector ] by the Cox-de Boor recursion formula: with p denoting the polynomial order and n representing the number of shape functions.Hence, NURBS enable the exact CAD-representation of the geometry and furthermore they provide C s−1 -continuity (s = min( p, q, r )) within a patch.The latter property facilitates the implementation of a conforming Ritz-Galerkin method of C s−1 -continuous shape functions for high-order partial differential equations.The corresponding approximation of an unknown field variable takes the form Our implementation has been accomplished into the GeoPDEs-package written (in Octave) for solving partial differential equations with isogeometric analysis [74] (see Appendix for details).In order to implement a solver obeying the strain gradient elasticity theory, the necessary modifications are confined to the assembly of the stiffness matrix, the right-hand side entity, and the boundary conditions.In the stiffness matrix module, the main task is to implement the operators of Eq. ( 19) by calling the first and second derivatives of the shape functions from the discrete space object and to loop over the elements of the mesh object.The same procedure applies to the right-hand side entity.

Numerical convergence analysis and model comparisons
In this section, the theoretical error estimates are first confirmed and the numerical implementation is verified by studying the convergence properties of the isogeometric method.The convergence properties of the method are studied with respect to the L 2 -, H 1 -and H 2 -norms.
Second, by utilizing the verified discretization, the present 3D model and the corresponding 1D beam models are compared with each other: the transversal displacement of the central axis and a normalized bending rigidity are reported for solid (3D) beams and compared to the exact solutions of the corresponding generalized (1D) Timoshenko and Euler-Bernoulli beam models.The constitutive parameters for the model comparisons are obtained by assuming isotropic elasticity augmented with one single length scale parameter.For a set of example structures having a microarchitecture, anisotropic versions of the constitutive parameters, both C i jkl and A i jklmn , are obtained via computational homogenization.Last but not least, another peculiarity of the strain gradient formulation, smoothening of stress singularities, is demonstrated for an L-shaped three-dimensional body.

Numerical convergence analysis -fully clamped cubes
In the first verification benchmark, we adopt a simplified version of the sixth-order gradient-elastic stiffness tensor A i jklmn : the components of the tensor are defined by a multiplication of the classical fourth-order elastic stiffness tensor C i jkl and the square of an intrinsic length scale parameter g in the form Let us consider a solid cubic structure Ω with the following material and geometrical properties: The body is under the influence of a given triaxial body load f , and the homogeneous Dirichlet boundary conditions are prescribed on all of the boundaries of the domain: In particular, the boundary condition type (fully clamped) is the same as used as an assumption for the theoretical results (for simplicity, as usual).Anyway, providing an exact solution for a problem of 3D strain gradient elasticity is very challenging, hence a numerical solution provided by a higher degree of NURBS functions serves as a reference solution for the error estimates: our reference solution is obtained with p = 4.
The convergence results and the corresponding theoretical convergence rates are plotted in Fig. 1.Figs. 1(a), 1(c) and 1(e) demonstrate the error of displacement, respectively, in the L 2 -norm, H 1 -and H 2 -seminorms with order p = 2, while Figs.1(b), 1(d) and 1(f) include the error in the L 2 -norm, H 1 -and H 2 -seminorms with order p = 3. Apart from the L 2 -error with order p = 2, the convergence rates of the error measures are in accordance with the theoretical convergence rates: O(h p−1 ) in the H 2 -seminorm; O(h p ) in the H 1 -seminorm; O(h p+1 ) in the L 2 -norm.
In the second verification benchmark, the values of the stiffness tensors (C i jkl and A i jklmn ) are obtained by applying a computational homogenization methodology to an (anisotropic) octet-truss unit cell: cell size 4 × 4 × 4, strut thicknesses 8 (horizontal) and 4 (inclined); effective density 0.756 (for illustrations, see [81,82] studying the extension and bending of 3D-printed octet-truss beams).The computational homogenization method for computing the fourth-order stiffness tensor C i jkl and the sixth-order stiffness tensor A i jklmn of the unit cell is based on equating the strain energies of the detailed (full-field, fine-grain) 3D model of the unit cell and the corresponding homogenized cell under proper boundary conditions (see [43] and a simplification in [44]).For an orthotropic material, 9 independent C i jkl -components are required together with 51 independent A i jklmn -components [23].
The corresponding convergence curves for the H 1 -and H 2 -seminorms with order p = 2 are shown in Fig. 2 showing practically identical, optimal, convergence rates as the first benchmark.These results confirm the theoretical error estimates and verify the implementation.

3D solid models versus 1D beam models
In this subsection, we compare the present 3D model of strain gradient elasticity to the corresponding 1D beam models in the contexts of uniaxial bending.32)): error for displacement in different norms (blue dashed lines) with the corresponding expected theoretical convergence rates (red solid lines).

Beam bending -parameter studies
For investigating the validity and sensitivity of two different 1D beam models of strain gradient elasticity, the bending behaviour of the 3D model (numerical solution) and the corresponding 1D beam models (analytical solutions) are compared to each other for a cantilever beam under a uniform transversal body load with the following material, geometrical and boundary conditions: In particular, by (35) we adopt the same material assumptions as in the first benchmark above (isotropic elasticity augmented with one length scale parameter) and accomplish parameter studies with respect to the length scale parameter g and the thickness-to-length ratio t/L.For six beams, with different combinations of length (L = 3, 5, 10) and the material length scale parameter (g = 0.1, 0.9), the deflection of the neutral axis is plotted in Fig. 3.For the sake of comparison, also the exact deflections of the corresponding strain gradient Euler-Bernoulli beam theory (gEB) and strain gradient Timoshenko beam theory (gT) are plotted in the figures.It is worth noting that the thickness-to-length ratios of the beams are t/L = 1/3, 1/5, 1/10, meaning very thick or thick cantilever beams.
An obvious observation is that the deflection of the 3D strain gradient model (g3D, blue circles) follows the Timoshenko beam model (gT, red solid line) in all of the cases -regardless of the thickness-to-length ratio or gradient contribution.Furthermore, the accuracy of the Euler-Bernoulli solution (gEB, green dashed line) decreases as the thickness-to-length ratio decreases -similarly to classical elasticity.Besides, the accuracy of the Euler-Bernoulli model decreases as the length scale parameter increases -ending up to a very clear deviation from the other two models (50% of the maximum deflection) in case of a very thick beam (t/L = 1/3) and relatively large gradient contribution (g = 0.9).

Beam bending -microarchitectural structures
In this subsection, we compare the present 3D model and the corresponding 1D beam models in the context of bending rigidity.Four different strain gradient models are compared to each other: -3D strain gradient model (A3D) with anisotropic material parameters C i jkl and A i jklmn (obtained via computational homogenization) -3D strain gradient model (g3D) with anisotropic material parameters C i jkl but one length scale parameter g (A i jklmn = g 2 δ mn C i jkl ) -1D strain gradient Timoshenko beam model (gT) with Young's modulus E, shear modulus G and one length scale parameter g -1D strain gradient Euler-Bernoulli beam model (gEB) with Young's modulus E and one length scale parameter g.
In the first example, the actual material parameters C i jkl , A i jklmn , E, G and g are obtained via computational homogenization methods applied to detailed finite element models of 3D classical elasticity for a set of comparable beam structures made of the same metamaterial and having the same thickness-to-length ratio -but the number of representative volume elements (RVE), or unit cells, in the thickness direction is different.A bending size effect has been observed for such microarchitectural beams in [29] by accomplishing a set of virtual experiments for planar lattices, but now we extrude the 2D lattice configurations in the direction of the x 2 -coordinate and form cellular anisotropic 3D structures (see Fig. 4) serving as reference structures for the generalized 3D and 1D models listed Table 1 Equitriangular unit cell: Dimensions and material properties.above.The geometrical and material properties of the unit cell, presented in Table 1 for the sake of clarity, have been introduced in [29,33].
For modelling the beam-like structures as homogenized prismatic 3D solids (see Fig. 5), the classical orthotropic (essentially transversally isotropic) elasticity determined by the cellular microarchitecture (via C i jkl , see Table 2) is augmented by two different types of constitutive models of gradient elasticity: the complete set of generalized constitutive parameters (A i jklmn ); the one-parameter (g) simplification already adopted in the previous subsections (but there combined with the C i jkl -values of isotropy).For the former, the material parameters (51 independent A i jklmn -values in orthotropy [23,24], less in transversal isotropy but not listed here) are obtained by the computational homogenization method already applied to the octet-truss unit cell in Section 5.1.For the latter, the length scale parameter can be obtained from the corresponding beam models sharing the same parameter (see [29]).
A series of beam structures is labelled according to the number of equilateral triangles in the thickness direction (x 3 in Fig. 4 (left)) -N = 1, 2, 4, 8, 16, 32, 64 -as depicted in Fig. 4 (right) for the first four beams of the series.As the unit cell of the metamaterial detailed in Fig. 4 (left) consists of two triangles, the corresponding series according to the number of unit cells in the thickness direction is actually n = 1/2, 1, 3/2, 2, 4, 8, 16, giving thicknesses t = nh 3 /2 and lengths L = 2nh 1 which result in thickness-to-length ratio t/L = (h 3 / h 1 )/4 (notified to mean very thick beams for configurations close to an equilaterally triangular honeycomb lattice resulting from the dimensions given in Table 1).
For modelling the beam structures by 1D beam models (with the central line following the x 1 -direction according to Fig. 4), parameter C 11 required for the classical beam models (for uniaxial bending in the x 1 x 3 -plane) needs to be augmented with the length scale parameter g.Namely, it has been shown in [29,33] -by using fine-grained 3D models of classical elasticity as references for the model validation -that in uniaxial bending these structures follow the strain gradient beam models.The g-value specific to the metamaterial has been obtained by a series of bending tests in [29,33], and the same value is utilized here for the g3D-model through the simplification defined in (32).According to this model simplification, the cantilever deformations of the homogenized 3D strain gradient solids are depicted in Fig. 5 for the first four beams of the series.Now that the material parameters are at hand for all of the model types (A3D, g3D, gT and gEB), bending rigidity, or flexural rigidity, can be determined as the resistance of the structure to an external loading: where F and δ are, respectively, the vertical applied force and the corresponding maximum displacement of the structure.According to the classical Euler-Bernoulli beam theory (EB), the corresponding theoretical bending rigidity of a cantilever beam is given as with E standing for Young's modulus of the material and L denoting the length of the beam, whereas I denotes the second moment of inertia of the cross-sectional profile, and the latter equality is valid for rectangular cross sections of a homogeneous material with thickness t and depth b.This identity shows, in particular, that beams made of the same material (E) and having constant depth (b) and constant thickness-to-length ratio (t/L) possess the same bending rigidity.It should be notified that for the series of beam structures illustrated in Fig. 4 these conditions are valid, apart the homogeneity -and the fact that the beams are very thick and are hence not expected to follow the EB beam theory, if any: . Indeed, the corresponding bending rigidity given by the classical Timoshenko beam theory (T) -appropriate for thicker beams -is given as where the last equality is valid for rectangular cross sections of homogeneous material (E and shear modulus G) with thickness t and depth b.This form shows, in particular, that for a fixed length rigidity D T approaches D E B when thickness approaches zero, whereas for constant thickness-to-length ratios rigidities D T and D E B stay apart (according to the product of ratio 3E/G and (t/L) 2 ).The corresponding strain gradient beam theories gEB and gT can be shown to follow, respectively, the formulae [29] D g E B ≈ 3 ) , ) . ( The latter forms for D g E B and D gT show how, with a fixed g-value, these rigidities rapidly increase for small thickness values, when compared to the classical counterparts D E B and D T , respectively -this character will be next connected to the size effect observed in [29,33] for the series of cellular beams of Fig. 4 now modelled as homogenized 3D strain gradient solids as illustrated in Fig. 5. Finally, for the microarchitectural beam structures working as cantilevers the normalized bending rigidities are plotted in Fig. 6 (left) for the different modelling options: -A3D (IGA, approximation, circle markers) -g3D (IGA, approximation, diamond markers) -gEB (analytical, exact, red dashed line) -gT (analytical, exact, blue dashed line); the last two following (41).The normalization has been done with respect to the classical EB-model according to (39) (black solid line), while the classical T-model follows (40) (dashed green line) -these rigidities present the corresponding prismatic beams made of a homogeneous material with Young's modulus C 11 .These plots confirm the qualitative consistency, with respect to the stiffening effect, between the 3D-models (A3D-model, g3D-model) and the 1D-model (gT-model) for thick beams (which, in turn, has already been shown in [29,33] to agree with the 3D fine-grain reference models of the beam structures of Fig. 4).However, the homogenization methodology used for obtaining the complete A-tensor leads to an essentially stiffer contribution (A3D, red circle markers for N = 1, 2, 4, 8) than the simplified approach of one g parameter (g3D, diamond markers for N = 1, 2, 4, 8), whereas the corresponding orthotropic C-tensors (Table 2) agree, which is confirmed by the overlapping markers  for N = 16, 32, 64.Finally, these plots show that for the strain gradient Euler-Bernoulli beam theory these beam structures are too thick, as expected (whereas similar but thinner microarchitectural beam structures have been shown to follow gEB [29] when stiffening).Deflections along the central line of one beam (N = 2), plotted in Fig. 6 (right; the black solid line as a reference), show that the gT-and g3D-model are the ones closest to the 3D fine-grain reference model (black solid line), whereas the A3D-model (red circles) is too stiff and the 3D model including the homogenized C-tensor alone (yellow squares) is too flexible.
For the second example, the equitriangular unit cell of Fig. 4 (left) is copied five times in the x 2 -direction and then doubled in the x 1 -direction forming a short tubular beam depicted in Fig. 7 (left): 5h 2 × 2h 1 × h 3 (length × depth × thickness in the coordinate directions indicated by the subscripts).Via a series of similar beams (but formed by repeating the unit cell in both x 2 -and x 3 -directions), it has been shown in [33] that such beams exhibit the stiffening effect in uniaxial bending in the x 2 x 3 -plane with g = 0.88 (mm).This value is now confirmed by clamping one end of the structure (at x 2 = 0) and bending the structure as a cantilever by a shear load at the free end in the x 3 -direction and matching the maximum deflection according to the corresponding beam model.Next, in order to form a biaxial bending state, a corresponding load is simultaneously applied in the x 1 -direction.For different models, the square sum of the two deflection values (u and w in the directions x 1 and x 3 , respectively) are plotted in Fig. 8 (right).The g3D-model agrees well with the reference model (a fine-grain 3D-solid model by a commercial FEA software Comsol), whereas the A3D-model is again too stiff and the 3D model including the homogenized C-tensor alone is too flexible, as in the previous example.Using only one g-parameter (validated now according to uniaxial bending in the x 2 x 3 -plane) applies reasonably well to this problem of biaxial bending as well since for the chosen tubular structure bending in the x 2 x 3 -plane is the deformation exhibiting strong stiffening.The stiffening effect is essentially weaker in the other bending direction which, however, does not suffer too much from the slightly incorrect g-contribution it gets via the one-parameter simplification.
According to the derivations of anisotropic strain gradient beam models in [82,83], it can be deduced that in the x 2 x 3 -plane bending, for instance, the g-parameter should essentially form (via equation ( 32)) the A x 2 x 2 x 3 x 2 x 2 x 3parameter associated to the strain gradient term mainly responsible for the possible stiffening effect.In the simplified model, however, the g-parameter finds its way to the other non-zero values of the A-tensor as well (according to the orthotropic C-tensor and Eq. ( 32)), now especially to the A x 2 x 2 x 1 x 2 x 2 x 1 -parameter associated to the term mainly responsible for the possible stiffening effect in the x 2 x 1 -plane bending -now weaker for the present structure, however.
After all, we can conclude that the models based on the 3D strain gradient theory, with the aid of the implemented conforming Ritz-Galerkin method providing the numerical results, exhibit bending behaviours qualitatively agreeing with 1D strain gradient beam models.The simplified g3D-model agrees with the 1D Timoshenko beam model quantitatively as well, in uniaxial bending.It is evident, however, that for more general deformation states a homogenization methodology -more reliable and more general than the ones applied here which have turned out to either produce too stiff homogenized solids or to involve too few independent parameters -should be developed in order to obtain a complete A3D-model supporting material anisotropies via both elasticity tensors, C i jkl and A i jklmn .

Smoothening stress singularities -L-shaped domain
This example demonstrates how strain gradient elasticity removes, or better smoothens, the stress singularities of classical elasticity appearing at reentrant corners.Let us consider a three-dimensional L-shaped body, depicted in Fig. 9 (left), with the following material properties and imposed boundary conditions: where Γ D corresponding to the nonhomogeneous Dirichlet boundary condition consists of the faces marked with the red colour in Fig. 9.The size of the body is determined by setting a = 1.Fig. 9 (right) shows the development of stress component σ x x at point P situating in the reentrant corner as the number of elements increases (from ca 100 to ca 1000000).In classical linear elasticity (LE), the value of the stress component grows unlimitedly as depicted by the blue solid line, whereas in the realm of strain gradient linear elasticity (gLE), the value is controlled due, and according, to the intrinsic length scale parameter: the red line for g = 0.001 is still very close to the diverging classical solution, while the yellow, magenta and green lines for g = 0.01, 0.03 and 0.04, respectively, already produce convergence to bounded stress values.

Conclusion
In this article, we have analysed and demonstrated the prospects of the theory of first strain gradient elasticity in the realm of three-dimensional solids by providing a mathematical formulation, a numerical method and a software implementation with error analysis and model comparisons.Finally, we have addressed some occasions in which the classical continuum mechanics crucially fails to capture the results of virtual experiments (which could be replaced by real experiments).
On the mathematical side, we have first provided variational formulations within a Sobolev space framework H 2 with solvability results and error estimates stemming from continuity and coercivity.For numerical results, we have implemented a numerical method by following the path of isogeometric analysis and provided a numerical convergence analysis for confirming the theoretical convergence results and for verifying the implementation.The convergence rates in the H 2 -, H 1 -and L 2 -norms turn out to be optimal.
Regarding model peculiarities, we have compared the 3D model to the corresponding 1D strain gradient Euler-Bernoulli and Timoshenko beam models in the cases of (1) an isotropic cantilever beam and (2) a transversally isotropic cantilever beam in uniaxial bending.The 3D strain gradient theory has been shown to exhibit the bending size effect -the smaller, the stiffer -similarly to the 1D strain gradient beam models (the Timoshenko model in the case of the chosen very thick beams).In the beam models, the term responsible for the stiffening effect becomes explicitly visible (see [57]) but in the present 3D model it is embedded in the somewhat complex formulation.The agreement between the 3D and 1D models, with respect to this effect, confirms that the accomplishment of the dimension reduction for reaching the beam models requires taking into account the full strain gradient of the 3D formulation (not only the derivative in the direction of the beam axis), which can be seen natural: the standard kinematical assumptions are inserted into the full 3D formulation and the standard mathematical steps (integration over the beam cross section, essentially) results in a variational formulation along the beam axis, and further the corresponding Euler equations (after integration by parts).Finally, we have applied the strain gradient model for the smoothening of stress singularities in the case of a reentrant corner (of an L-shaped three-dimensional domain) by investigating the effect of an intrinsic length scale parameter possessed by the model.Applying the present 3D strain gradient elasticity model with the present numerical method to more general structures requires the application of generalized homogenization methods for obtaining more than one length scale parameter (cf. the 2D problem setting in [84] for lattice beams) -although in many applications using just the generalized dimension reduction models is a more robust, and according to the present model comparison a valid and accurate, option for many beam-like structures at least.The examples of the present study call for a general and reliable computational homogenization methodology -more general and more reliable than the ones applied above which have turned out to either produce too stiff homogenized solids or to involve too few independent parameters -in order to obtain 3D material tensors supporting material anisotropies.One interesting application in this direction would be embedding pantographic structures, for which a fairly solid theoretical understanding is already available [26,[85][86][87][88], into the present three-dimensional framework, within both linear and geometrically nonlinear regimes.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix. GeoPDEs implementation
GeoPDEs [74] is an open-source software for isogeometric analysis (see Fig. 10), written in Octave compatible with Matlab.The essential features of GeoPDEs are the modularity, dimensionless implementation, providing multipatches, and the vectorization.
A solver in GeoPDEs is initiated by obtaining initial data such as the geometry, required continuity, loadings, and boundary conditions.In the first step, a geometry structure consisting of the parametrization of the geometry and its first and second derivatives is created in the solver block.In the second step, the mesh class instantiates an object holding the necessary information for the estimation of the quadrature rules and parametrization.In the third step, the space class takes the NURBS basis functions of the geometry structure and computes the basis functions and their derivatives at given quadrature points of the mesh object.Thenceforth, the stiffness matrix and right-hand side entity are assembled, and the boundary conditions are imposed.Eventually, the unknown degrees of freedom are obtained by solving the linear system.

Proposition 1 .
Let us assume that Γ = Γ D and U = V .There exists a positive constant C such that a

Fig. 1 .
Fig. 1.A clamped unit cube with an isotropic C-tensor and one single length scale parameter g (according to Eq. (32)): error for displacement in different norms (blue dashed lines) with the corresponding expected theoretical convergence rates (red solid lines).

Fig. 2 .
Fig. 2. A clamped unit cube with anisotropic Cand A-tensors according to computational homogenization for an octet-truss unit cell: error for deflection in the H 1 -norm (left) and H 2 -norm (right) with order p = 2.

Fig. 3 .
Fig. 3. Cantilever beams with an isotropic C-tensor and one length scale parameter g: deflection distribution along the neutral axis for different length L and length scale g values with p = 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Equitriangular unit cell -beams in the x 1 -direction: normalized bending rigidity for 3D and 1D models.. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Equitriangular unit cell -a beam in the x 2 -direction: 3D fine-grain model (left); square sum of deflections (u and w) in biaxial bending along the central line coordinate (s) (right).

Fig. 9 .
Fig. 9.An L-shaped solid body: geometry (left); the development of stress component σ x x at point P in mesh refinements with different values of the length scale parameter g (right).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
[75]r a time period of t ∈ [t 1 , t 2 ].W ext stands for the work done by the external forces.The kinetic energy density takes the form[75]

Table 2
Equitriangular unit cell: Coefficients of the fourth-order stiffness tensor C i j [GPa].