Verification of strain gradient elasticity computation by analytical solutions

As there are different computational methods for simulating problems in generalized mechanics, we present simple applications and their closed‐form solutions for verifying a numerical implementation. For such a benchmark, we utilize these analytical solutions and examine three‐dimensional numerical simulations by the finite element method (FEM) using IsoGeometric Analysis (IGA) with the aid of open source codes, called tIGAr, developed within the FEniCS platform. A study for the so‐called wedge forces and double tractions help to comprehend their roles in the displacement solution as well as examine the significance by comparing to the closed form solutions for given boundary conditions. It is found that numerical results are in a good agreement with the analytical solutions if wedge forces and double tractions are considered. It is also presented how the wedge forces become necessary in order to maintain equilibrium in strain gradient materials.

formation and propagation, may be realized by involving higher order gradients in order to regularize the solution of cracks by penalizing the displacements, which are localized above a given threshold [17,68,71,73,78,93]. Second, an accurate modeling of metamaterials is achieved in the framework of strain gradient theory [16,26,28,35,38,40,45,89]. Indeed, the behaviors of metamaterials depend on the morphology of microstructures [23,43,53,54] and show peculiar behaviors, such as size effects [3,[90][91][92] and band gaps [76,80], which are not captured by the classical Cauchy continuum. Exotic material behaviors may be governed by internal or boundary layers of higher strain gradients [51]. As an example of metamaterials we refer to materials with pantographic substructures [32,34] having non-negligible strain gradient terms in the formulation of deformation energy [24,25,41,65,88], which are found in recent numerical and experimental evidences, for example, in shearing tests [15], three-point bending tests [94], compression tests [87], or torsion tests [60]. A homogenization procedure of pantographic microstructures leading to second gradient materials is shown in [12][13][14]84]. Generally speaking, generalized mechanics at the macroscale results in as a homogenization procedure from the Cauchy-Boltzmann continuum at the micrometer length-scale (microscale) [57]. Strain gradient theory is able to sustain the so-called wedge forces [9,56,77], which are acting on corners (point forces) in 2D and edges (force distributions along lines) in 3D, while the classical Cauchy continuum leads to singularities when handling such issues. From [8,33,48,83], it becomes obvious that the reasons of inducing higher order strain gradient allows a continuum to sustain boundary conditions on vertices and edges of a body.
Strain gradient theory leads to higher order partial differential equations and requires an interpolation scheme for the finite element representation to guarantee a correspondingly higher order of continuity. For this reason, various numerical implementations of strain gradient theories are proposed in the literature, such as the mixed formulation [67,85], 1 continuous elements [42,66,95], and isogeometric analysis [39,56,81]. The concept of IGA [46] is a mesh-based numerical approach using shape functions in FEM identical as the basic functions in CAD generated from NURBS (Non-Uniform Rational B-Splines). In the context of strain gradient elasticity, such a formulation serves the continuity across the element boundaries. The NURBS interpolation requires only displacements as nodal degrees of freedom, and no derivatives of the displacement field are needed [39]. Some examples can be found in [19][20][21]82]. From a general viewpoint, for verification purposes, a numerical solution have to be compared to analytical results [79]. Some analytical derivations were presented for 2D cases [22,69,74,75,95]. Obtaining results of this kind for 3D cases is rather challenging.
In this paper, closed-form analytical solutions [70,72,74] are presented for cases of three dimensional benchmark problems for a centro-symmetric and isotropic material modeled by strain gradient theory. The balance equations and the boundary terms (tractions, double tractions, and wedge forces on the boundaries) are derived by means of the variational method. Numerical implementations have been developed based on IGA (Iso-Geometric Analysis) for these classical cases in order to demonstrate the benchmarking procedure. The aim of the paper is to verify our implementation based on an IGA library called tIGAr [49,50] provided on the FEniCS platform. Moreover, a comparative study allows us to understand the roles of wedge forces in maintaining equilibrium. The code uses open-source packages under GNU public license [44], and we make the codes publicly available in [2] in order to enable a scientific exchange. The paper is organized as follows: The variational formulation of strain gradient elasticity theory [1,30] is outlined in Section 2. The formulation of threedimensional problem is shown in Section 3. In Section 4 analytical and numerical results are compared for three different benchmark examples.

Variational formulation
In this section, strain gradient elasticity as in [58,59] will be revisited. Conventional continuum mechanics theories assume that stress at a material point is a function of state variables, such as strain, at the same point. This local assumption is adequate when the wavelength of a deformation field is orders in magnitude greater than the dominant micro-structural material length scale. However, when these two length scales are comparable, an extension becomes necessary, herein we use the strain gradient theory. Unlike the classical elasticity, in strain gradient theory the stored energy density depends on the strain, , and on the strain gradient, : Here denotes the displacement field and a comma means differentiation in space, , expressed in Cartesian coordinates, For strain gradient materials, the stored energy density, , depends on the first and second gradients of the displacement field = ( , ).
We compute the first variation of functional for the internal energy of the body int : such that we have where the variation of displacement δ is the test function in the finite element method to be used in simulations. After applying integration by parts, for the details we refer to the Appendix, we obtain In the third integral of the right hand side of the equation above, on the boundary surface, the gradient of test function δ , can be decomposed into a gradient within and normal to the surface [10,18,47,52]: where the operators and read and is the unit surface normal vector. Considering Equation (7), the last integral in Equation (6) becomes The boundary surface is assumed to be divisible into a finite number of smooth parts, Ω , each bounded by an edge, Ω . Using Stokes' divergence theorem by following [18] on each smooth surface, the following equation is obtained where is the unit tangent vector and belongs to the tangent space to Ω . According to the chain rule, we have Then the first integral on the right hand side of Equation (9) is written as By Δ(⋅), we denote the difference between values of the expression in the parentheses, which are calculated at different sides of a sharp edge. Note that integration domain of the third integral in Equation (10) is the boundary of a smooth surface. The whole body is composed of several smooth surfaces. By ∑ , the sharp edges of the body are summed up. We stress that every edge line enters twice since it belongs to two adjacent surface regions [18]. By inserting the latter into Equation (11), the last surface integral in Equation (6) becomes Thus the final integral form reads For a better analogy, we define stress and hyperstress as follows: It is observed from the last two integrals in Equation (14) that second gradient continua can sustain external surface double forces and sustain external line forces [10]. According to the principle of virtual work, ext is the external work done on the body. It is assumed to have the form [10,52]: The first part, , is the energy density because of the specific force, [4]. By Equations (14), (16), and (17), the so-called tractions , double tractions and wedge forces are expressed as By inserting Equation (15), Equation (18) into Equation (16) with the aid of Equation (14), we find out the so-called governing equation: F I G U R E 1 The schematic of a 3D block. A 3D block with numbered edges. (b) Unit tangential vector and unit normal vector

Constitutive laws
Strain energy density in case of centrosymmetric materials reads where and are the first and the second gradient elastic stiffness tensors, respectively. For isotropic materials, they are given by By 1 and 2 the Lamé constants are denoted. It is worth mentioning here that for a 3D case there are five additional parameters 3 , 4 , 5 , 6 , and 7 in the second gradient stiffness tensor . By using Equations (20)- (22) we can rewrite Equation (15) as

FORMULATION OF THE PROBLEM
A 3D block, as shown in Figure 1, is considered in the current section. The motivation for this selection is twofold. First, all surface boundaries are flat, consequently, the boundary conditions are simplified. Second, the presence of sharp edges would bring the wedge forces, and their importance will be demonstrated later on. We choose the length , width , and height ℎ to be equal ( = = ℎ). Each of the 12 edges is characterized by a unique number, see Figure 1a. Unit tangential vectors and unit normal vectors for the right and top surfaces are presented in Figure 1b.
For the considered domain Ω its boundary Ω is composed of flat surfaces leading to ( ) = 0 on Ω. Hence Equation (18) According to Equation (18), the so-called wedge forces are There are 12 edges in total for the 3D block. Let us take the edge number 6, which is blue in Figure 1b, as an example. The unit normal vectors and tangent vectors for the right and top surfaces are The wedge force on the edge number 6 is then calculated by Expressions for tractions, double tractions, and wedge forces on the remaining surfaces and edges respectively can be calculated analogously. Expressions for the stress tensor, the hyperstress tensor, and balance equations can be given explicitly in terms of the displacement by using Equations (23), (24), and (19) as shown in Appendix.

ANALYTICAL AND NUMERICAL SOLUTIONS AND COMPARATIVE STUDIES OF A SECOND GRADIENT MODEL FOR A 3D BLOCK
Let us consider the deformation modes of a block corresponding to particular components of the strain gradient tensor as shown in Table 1 [6,55]. We emphasize that for underlined indices no summation convention is applied. There are four distinct modes of deformations: Extension, torsion, non-conventional bending, and trapezoid [55,72,75]. These four modes correspond to specific non-zero strain gradient components, which will result in strain gradient effects of materials. For example, a trapezoid loading deforms the block to be trapezoid-shaped with non-zero strain gradient components 121 as shown in [55]. In this work, the extension, torsion, and non-conventional bending deformation modes are studied by considering specific displacement field solution corresponding to these non-zero strain gradient components as illustrated in Figure 2 Table 2. The values of strain gradient related parameters may be possibly given as shown in Table 2 [9]. For defining the body forces, gravitational forces are implemented by the specific force =10 N/kg. Three cases of different boundary conditions are investigated throughout this section.

Extension
The following displacement fields leading to a non-zero 222 are considered: This solution of displacement fields can be achieved by the following body forces: Due to the kinematic restrictions, the wedge forces are not imposed in this case. The imposed kinematical constraints and boundary conditions are also shown in Figure 3 for a cut slice. Plots for numerically obtained total displacement are presented in Figure 4. In Figure 5, perfect overlap is observed between the analytical and numerical results. A difference is shown for numerical solutions if double tractions are not taken into account.

Torsion
Lets now consider a torsion problem with the following displacement fields: where is a small constant with the physical dimension of the inverse of a length ( is set to be equal to 0.1 m −1 ). The analytical solution is achieved by the following body forces:  (36) Figure 6 presents the imposed kinematical constraints and boundary conditions for a cut slice. Comparing the exact solution with the numerical outcomes on three different cuts (line1, line2, line3), see Figure 7, we observe that the plot for simulation incorporating wedge forces demonstrate an almost perfect overlapping with the plot of the analytical solution. Figure 8 indicates differences for total displacement plots for different boundary conditions. It is well known that a sin-

Non-conventional bending
The so-called non-conventional bending problem is considered. The displacement field equations are conceived as where is a small constant with the physical dimension of the inverse of a length ( is set to be equal to 0.1 m −1 ). The solution of the displacement field is acquired by the following body forces: The following wedge forces calculated by Equation (27) The imposed kinematical constraints and boundary conditions are presented in Figure 9 for a cut slice. A comparison between the analytical solutions and the numerical ones along a line cut (line 1) are shown in Figure 10. Obviously, when the double tractions and wedge forces are both imposed, the numerical results are in good agreement with the analytical solutions as also observed in Figure 11.
The strain energies are calculated and presented here for the analytical and numerical solutions as shown in Figure 12.

CONCLUSIONS
In this paper, numerical simulations based on IGA are performed in order to be compared with an analytical solution obtained using the inverse method. Three exemplary cases of different boundary conditions were considered in 3D: Extension, torsion, and non-conventional bending. Numerical simulations were performed with and without taking the socalled wedge forces and double tractions into account. It was shown that the numerical results are in good agreement with the analytical solutions if the wedge forces and double tractions are considered. Besides that, we presented comparisons between the numerical simulations. Such a comparison helps to reveal and comprehend the roles of the wedge forces and double tractions in the solutions. The numerical implementation was conducted by means of an open source tool called tIGAr which is based on FEniCS library.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL. The weak form and its numerical implementation According to Equations (4), (5), (16), and (17), the weak form is shown as

O R C I D
where ( ) is equal to , and is the unit surface normal vector. For example, the of the right surface in Figure 1(b) is equal to 1 .
In this work, the wedge forces are constants and uniformly distributed on each edge with the unit of N/m. The total force acting on each edge is equal to in N. in m is equal to the length of the corresponding edge. For example, in Section 4.2, the wedge forces on edge 2 is 1 . The total force acting on edge 2 is 1 . Instead of implementing the wedge forces directly, a traction with the amplitude of 1 2 in Pa is applied on a small region (area of the region is 2 ) near edge 2 as shown in Figure A.1. For a detailed introduction of IGA we refer to [46]. In tIGAr, a self-contained implementation of single-patch explicit B-splines is realized by using a module called tIGAr.BSplines [50]. In this work, polynomial degree of the basis function has been chosen as 2 in order to acquire 1 -continuity for the formulation [39,81].

A.2
The derivation of Equation (6) The derivation of Equation (6) is shown here. Equation (6)  Therefore, the second term in the left hand side in Equation (6)
(A. 15) The other components can be calculated in the same manner. Equation (19) together with Equations (23)-(24) gives us the system of partial differential equations. For example, the balance equation in 1 direction is (A. 16)