A benchmark strain gradient elasticity solution in two-dimensions for verifying computational approaches by means of the finite element method

Abstract

In elasticity, microstructure-related deviations may be modeled by strain gradient elasticity. For so-called metamaterials, different implementations are possible for solving strain gradient elasticity problems numerically. Analytical solutions of simple problems are used to verify the numerical approach. We demonstrate such a case in a two-dimensional continuum as a benchmark case for computations. As strain gradient enforces higher regularity conditions in displacements, in the finite element method (FEM), the use of standard elements is often seen as inadequate. For such piecewise or elementwise continuous elements, we examine a possible remedy to correctly simulate strain gradient elasticity problems by implementing two techniques. First, we enforce continuity of displacement gradient across elements; second, we employ a mixed finite element method where displacement and its gradient are solved both as unknowns. The results show the pros and cons of each numerical technique. All methods converge monotonically, but the mixed method is more reliable than the other one.

Keywords

Strain gradient elasticity finite element method benchmark verification analytical solution

1. Introduction

Accurate modeling of materials at microscopic length scale is necessary for design and engineering of structures with a detailed and complex substructure. The first-order (first-gradient) continuum mechanics depends on displacement and its first gradient; however, it lacks accuracy at the scales which are close to the length scale of the substructure [1], especially in capturing those phenomena which are dependent on geometric length, known as size effect. A generalization of the first-order elasticity theory is an approach for overcoming this issue. In generalized continuum mechanics, the governing equation depends on displacement, first gradient of displacement, and second gradient of displacement.

The idea of generalized mechanics dates back to the beginning of 20th century, see [2,3] for a history. Modern theories were introduced later starting around six decades ago [4 –7]. Since then, many generalized mechanics theories have been proposed, and they can be considered as specific cases of a unified theory [8].

Generalized mechanics has been widely investigated in the literature. It has been implemented for problems of elasticity [9 –13]; plasticity [14 –19]; damage modeling [20 –25]; modeling metamaterials [26 –28] such as pantographic structures [29 –31], network materials [32], viscoelastic truss structures [33], bi-pantographic structures [34], second gradient fluids [35]; gradient-enhanced homogenization [36 –39]; micropolar continua [40]; fracture mechanics [41]; biomechanics [42 –44]; and anisotropic systems [45]. Parameter determination of generalized mechanics models has been studied for static and dynamic regimes in Shekarchizadeh et al. [46, 47], respectively.

Including higher gradients of displacement in the equations results in partial differential equations of higher order. A reliable numerical computation of such equations requires suitable techniques and element type selection that ensure the monotonous convergence. For this purpose, different numerical approaches have been proposed for the strain gradient theories such as isogeometric analysis [47 –50], $C^{1}$ continuous elements [52,53], and mixed finite element formulation [54 –56]. It is beneficial to verify the computations by analytical solutions. The analytical solutions for some example problems in the generalized continua are presented, for example, in previous studies [57 –61].

In this paper, a two-dimensional problem in the framework of the strain gradient elasticity theory is solved numerically, by means of a finite element method (FEM) implemented using open-source FEniCS libraries. Different techniques are used for the implementation of the numerical code including using mixed finite element formulation, using Lagrange multipliers for imposing the boundary conditions, and enforcing continuity of first gradient across elements. The computations are compared with analytical solutions. Different implementations are compared to each other regarding convergence, computation time, and robustness.

This paper is structured in the following way. First, in Section 2, the strain gradient theory is explained, the weak form is generated to be used in the analytical solution, and the constitutive equations are presented. Then, in Section 3, an analytical solution is presented for a plate under simple shear with two different sets of boundary conditions. Next, in Section 4, the numerical implementation and weak form generation of FEM and Mixed FEM are sketched in detail. Finally, the results and error analyses of two implementations of the strain gradient elasticity theory are shown in Section 5.

2. Strain gradient elasticity

2.1. Weak form generation

One of the well-known generalizations of the conventional continuum mechanics is the strain gradient elasticity theory. In this theory, the energy density depends not only on strain, $ε_{ij}$ , but also on its gradient, $ε_{ij, k}$ , where

\begin{matrix} ε_{ij} = \frac{1}{2} (u_{i, j} + u_{j, i}), \\ ε_{ij, k} = \frac{1}{2} (u_{i, jk} + u_{j, ik}), \end{matrix}

(1)

where $u_{i}$ is the displacement field and the comma indicates differentiation in space, $X \in Ω \subset {I R}^{3}$ , in Cartesian coordinates. Therefore, we have

\begin{matrix} u_{i, j} = \frac{\partial u_{i}}{\partial X_{j}}, \\ u_{i, jk} = {(u_{i, j})}_{, k} = \frac{\partial}{\partial X_{k}} (\frac{\partial u_{i}}{\partial X_{j}}) = \frac{\partial^{2} u_{i}}{\partial X_{k} \partial X_{j}} . \end{matrix}

(2)

In equation (1), the geometric nonlinearities are neglected since the strain measure has to be linear for deriving an analytical solution in the following. We begin with the general formulation based on a scalar mathematical construct called action $A$ . We postulate the action in the time interval $T$ as

A = \int_{T} (\int_{Ω} L d V + \int_{\partial Ω} W_{s} d A + \int_{\partial \partial Ω} W_{e} d L) d t,

(3)

over the domain, $Ω$ , with its boundary, $\partial Ω$ , and the set $\partial \partial Ω$ of the edges of the boundary surface. In a three-dimensional problem, the terms $W_{s}$ and $W_{e}$ are the work done on the boundary surface elements, $d A$ , and line (edge) elements, $d L$ , respectively. Herein, we neglect the term $W_{e}$ for simplicity [62,63]. The values of $W_{s}$ are known, and as this term is applied on the boundary, for the same accuracy in an expansion up to the second order in space derivative, we set $W_{s} = W_{s} (u_{i}, u_{i, j})$ , i.e., depending only on the displacement and its first gradient but not on the second gradient.

In equation (3), $L$ is an existing Lagrangean density describing the underlying system. In the strain gradient elasticity theory for an elastic material, as discussed in Abali [64], we consider the Lagrangean density as

L = \frac{1}{2} ρ_{0} ({\overset{\cdot}{u}}_{i} {\overset{\cdot}{u}}_{i} + d^{2} {\overset{\cdot}{u}}_{i, j} {\overset{\cdot}{u}}_{i, j} + τ^{2} {\overset{\cdot\cdot}{u}}_{i} {\overset{\cdot\cdot}{u}}_{i}) - w + ρ_{0} f_{i} u_{i},

(4)

where the first three terms indicate the kinetic energy (inertial terms). The term $w$ is the stored (deformation) energy density in $J / m^{3}$ . For homogeneous materials, the stored energy density reads

w = w (u_{i, j}, u_{i, jk}) .

(5)

The last term in equation (4) denotes the potential energy with the specific body force, $f$ , in N/kg. Based on the principle of least action, for reversible systems, we are looking for solutions such that the variation of the action functional vanishes for any arbitrary test function $δ u$

δ A = 0 \forall δ u .

(6)

The variation of the action, $δ A$ , is derived by employing Taylor’s expansion with respect to a suitable perturbation parameter, as explained in Abali et al. [65]. By inserting $A$ from equation (3) in equation (6), and neglecting the inertial terms, body forces, and boundary terms acting on edges, according to the least action principle, the integral form for the second-order strain gradient elasticity for a domain $Ω$ reads

- \int_{Ω} (\frac{\partial w}{\partial u_{i, j}} δ u_{i, j} + \frac{\partial w}{\partial u_{i, jk}} δ u_{i, jk}) d V + \int_{\partial Ω} (\frac{\partial W_{s}}{\partial u_{i}} δ u_{i} + \frac{\partial W_{s}}{\partial u_{i, j}} δ u_{i, j}) d A = 0 .

(7)

After applying the product rule and the divergence (Gauss’s) theorem multiple times (for details see Appendix 1), we obtain

\begin{matrix} \int_{Ω} {(- \frac{\partial w}{\partial u_{i, j}} + {(\frac{\partial w}{\partial u_{i, jk}})}_{, k})}_{, j} δ u_{i} d V + \int_{\partial Ω} [n_{j} (\frac{\partial w}{\partial u_{i, j}} - {(\frac{\partial w}{\partial u_{i, jk}})}_{, k}) - \frac{\partial W_{s}}{\partial u_{i}}] δ u_{i} d A \\ + \int_{\partial Ω} (n_{k} \frac{\partial w}{\partial u_{i, jk}} - \frac{\partial W_{s}}{\partial u_{i, j}}) δ u_{i, j} d A = 0, \end{matrix}

(8)

since, as already remarked, the inertial terms, body forces, and boundary terms acting on edges are neglected. Here, $n$ is the outward unit normal to the boundary $\partial Ω$ . We define $W_{s}$ on Neumann boundaries as

W_{s} = t_{i} u_{i} + m_{ij} u_{i, j},

(9)

where $t$ , in N/m² (traction), and $m$ , in N m/m² (double traction), satisfy the following conditions

t_{i} = \frac{\partial W_{s}}{\partial u_{i}} = n_{j} (\frac{\partial w}{\partial u_{i, j}} - {(\frac{\partial w}{\partial u_{i, jk}})}_{, k}),

(10)

m_{ij} = \frac{\partial W_{s}}{\partial u_{i, j}} = n_{k} \frac{\partial w}{\partial u_{i, jk}} .

(11)

By inserting equations (10) and (11) in equation (8), the following governing equation is obtained

{(\frac{\partial w}{\partial u_{i, j}} - {(\frac{\partial w}{\partial u_{i, jk}})}_{, k})}_{, j} = 0,

(12)

which is solved numerically to calculate the displacement field in the domain $Ω$ .

2.2. Constitutive laws

The stored energy density, $w$ , in equation (12), for centro-symmetric materials is expressed as

w = w (ε_{ij}, ε_{ij, k}) = \frac{1}{2} ε_{ij} C_{ijkl} ε_{kl} + \frac{1}{2} ε_{ij, k} D_{ijklmn} ε_{lm, n},

(13)

where $C_{ijkl}$ denotes the rank-4 stiffness tensor (first gradient), and $D_{ijklmn}$ indicates the rank-6 strain gradient stiffness tensor (second gradient) [66 –70]. For isotropic materials, we have

C_{ijkl} = c_{1} δ_{ij} δ_{kl} + c_{2} (δ_{ik} δ_{jl} + δ_{il} δ_{jk}),

(14)

\begin{matrix} D_{ijklmn} & = c_{3} (δ_{ij} δ_{kl} δ_{mn} + δ_{in} δ_{jk} δ_{lm} + δ_{ij} δ_{km} δ_{\ln} + δ_{ik} δ_{jn} δ_{lm}) + c_{4} δ_{ij} δ_{kn} δ_{ml} \\ + c_{5} (δ_{ik} δ_{jl} δ_{mn} + δ_{im} δ_{jk} δ_{\ln} + δ_{ik} δ_{jm} δ_{\ln} + δ_{il} δ_{jk} δ_{mn}) + c_{6} (δ_{il} δ_{jm} δ_{kn} + δ_{im} δ_{jl} δ_{kn}) \\ + c_{7} (δ_{il} δ_{jn} δ_{mk} + δ_{im} δ_{jn} δ_{lk} + δ_{in} δ_{jl} δ_{km} + δ_{in} δ_{jm} δ_{kl}), \end{matrix}

(15)

where $c_{1}, c_{2}$ are the two Lame constants, and $c_{3}, c_{4}, c_{5}, c_{6}, c_{7}$ are five additional parameters characterizing the substructure of the material; we refer to Abali et al. [71] for a derivation of this form. After inserting equations (14) and (15) in equation (13), the stored energy becomes

\begin{matrix} w = \frac{1}{2} c_{1} ε_{ii} ε_{jj} + c_{2} ε_{ij} ε_{ij} + 2 c_{3} ε_{ik, i} ε_{jj, k} + \frac{1}{2} c_{4} ε_{jj, i} ε_{kk, i} + 2 c_{5} ε_{ik, i} ε_{jk, j} \\ + c_{6} ε_{jk, i} ε_{jk, i} + 2 c_{7} ε_{jk, i} ε_{ji, k}, \end{matrix}

(16)

which is equivalent to the potential energy density expressed in Mindlin’s work [72, equation (11.3) on page 71]. There is a one-to-one relation between the material parameters used in our formulation $(c_{1}, c_{2}, c_{3}, c_{4}, c_{5}, c_{6}, c_{7})$ and in Mindlin’s formulation $(\tilde{λ}, \tilde{μ}, {\hat{a}}_{1}, {\hat{a}}_{2}, {\hat{a}}_{3}, {\hat{a}}_{4}, {\hat{a}}_{5})$ as follows:

c_{1} = \tilde{λ}, c_{2} = \tilde{μ}, c_{3} = \frac{1}{2} {\hat{a}}_{1}, c_{4} = 2 {\hat{a}}_{2}, c_{5} = \frac{1}{2} {\hat{a}}_{3}, c_{6} = {\hat{a}}_{4}, c_{7} = \frac{1}{2} {\hat{a}}_{5} .

(17)

The derivative of the stored energy density with respect to the first gradient of displacement is derived using equations (13), (1), and the chain rule as

\begin{matrix} \frac{\partial w}{\partial u_{i, j}} = & \frac{\partial w}{\partial ε_{qr}} \frac{\partial ε_{qr}}{\partial u_{i, j}} = \frac{1}{2} \frac{\partial (ε_{kl} C_{klmn} ε_{mn} + ε_{kl, m} D_{klmnop} ε_{no, p})}{\partial ε_{qr}} \frac{1}{2} \frac{\partial (u_{q, r} + u_{r, q})}{\partial u_{i, j}} \\ = & \frac{1}{4} C_{klmn} (δ_{kq} δ_{lr} ε_{mn} + ε_{kl} δ_{mq} δ_{nr}) (δ_{qi} δ_{rj} + δ_{ri} δ_{qj}) = \frac{1}{4} (C_{qrmn} ε_{mn} + ε_{kl} C_{klqr}) (δ_{qi} δ_{rj} + δ_{ri} δ_{qj}) \\ = & \frac{1}{4} (ε_{mn} (C_{ijmn} + C_{jimn}) + ε_{kl} (C_{klij} + C_{klji})) = \frac{1}{4} ε_{kl} (C_{ijkl} + C_{jikl} + C_{klij} + C_{klji}) = C_{ijkl} ε_{kl}, \end{matrix}

(18)

by considering the symmetry of the stiffness tensor $(C_{ijkl} = C_{jikl} = C_{klij} = C_{klji})$ . In the same way, we get the derivative of the stored energy density with respect to the second gradient of displacement as

\begin{matrix} \frac{\partial w}{\partial u_{i, jk}} = & \frac{\partial w}{\partial ε_{rs, t}} \frac{\partial ε_{rs, t}}{\partial u_{i, jk}} = \frac{1}{2} \frac{\partial (ε_{lm} C_{lmno} ε_{no} + ε_{lm, n} D_{lmnopq} ε_{op, q})}{\partial ε_{rs, t}} \frac{1}{2} \frac{\partial (u_{r, st} + u_{s, rt})}{\partial u_{i, jk}} \\ = & \frac{1}{4} D_{lmnopq} (δ_{lr} δ_{ms} δ_{nt} ε_{op, q} + ε_{lm, n} δ_{or} δ_{ps} δ_{qt}) (δ_{ri} δ_{sj} δ_{tk} + δ_{si} δ_{rj} δ_{tk}) \\ = & \frac{1}{4} (D_{rstopq} ε_{op, q} + ε_{lm, n} D_{lmnrst}) (δ_{ri} δ_{sj} δ_{tk} + δ_{si} δ_{rj} δ_{tk}) \\ = & \frac{1}{4} ((D_{ijkopq} + D_{jikopq}) ε_{op, q} + ε_{lm, n} (D_{lmnijk} + D_{lmnjik})) \\ = & \frac{1}{4} ε_{lm, n} (D_{ijklmn} + D_{jiklmn} + D_{lmnijk} + D_{lmnjik}) = D_{ijklmn} ε_{lm, n}, \end{matrix}

(19)

by considering the symmetry of the strain gradient stiffness tensor $(D_{lmnijk} = D_{lmnjik} = D_{ijklmn} = D_{jiklmn})$ . By inserting equation (14) in equation (18), and equation (15) in equation (19), we have

\frac{\partial w}{\partial u_{i, j}} = c_{1} δ_{ij} ε_{kk} + 2 c_{2} ε_{ij},

(20)

\begin{matrix} \frac{\partial w}{\partial u_{i, jk}} = & c_{3} (δ_{ij} ε_{km, m} + δ_{jk} ε_{mm, i} + δ_{ij} ε_{nk, n} + δ_{ik} ε_{mm, j}) + c_{4} δ_{ij} ε_{mm, k} \\ + & c_{5} (δ_{ik} ε_{jn, n} + δ_{jk} ε_{li, l} + δ_{ik} ε_{nj, n} + δ_{jk} ε_{im, m}) + 2 c_{6} ε_{ij, k} + 2 c_{7} (ε_{ik, j} + ε_{kj, i}) . \end{matrix}

(21)

Finally, by inserting equations (20) and (21) in equation (12), the governing equation is derived in terms of the material parameters and the gradients of strain as

\begin{matrix} c_{1} ε_{kk, i} + 2 c_{2} ε_{ij, j} - c_{3} (2 ε_{nk, nki} + ε_{mm, ikk} + ε_{mm, jij}) - c_{4} ε_{mm, kki} \\ - 2 c_{5} (ε_{nj, nij} + ε_{im, mkk}) - 2 c_{6} ε_{ji, kkj} - 2 c_{7} (ε_{ik, jkj} + ε_{kj, ikj}) = 0 . \end{matrix}

(22)

Equation (22) corresponds to the equation of motion derived in Mindlin’s work [72, equation (11.8) on page 72] (with the same one-to-one relation between the material parameters as in equation (17)). Herein, we have neglected the body forces and inertial terms as well.

3. Analytical solution

In this section, the analytical solution for a plate under simple shear is derived. Two cases of loading and boundary conditions are investigated.

Consider a two-dimensional plate under simple shear as in Figures 1 and 2 with infinite length in $x$ -direction, and with height $H$ in $y$ -direction. The bottom edge is fixed in both directions while the top edge is fixed only in vertical $(y)$ direction, and a displacement or traction is applied on the top edge in $x$ -direction.

Figure 1.

The plate (infinitely long in $x$ -direction) under simple shear with the boundary conditions stated in equation (26). The top edge is moving in a rail, therefore its rotation is not allowed.

Figure 2.

The plate (infinitely long in $x$ -direction) under simple shear with the boundary conditions stated in equation (33).

Considering the boundary conditions of this problem, only $u_{x}$ and its derivatives along $y$ -direction are non-zero. This so-called semi-analytical ansatz reduces the problem to a one-dimensional continuum. Therefore, the problem is one-dimensional with $u_{y} = 0$ , $u_{x} = u_{x} (y) = u$ . Hence, we have $u_{x, y} = u'$ , $u_{x, yy} = u ″$ , $u_{x, yyy} = u ″'$ , $u_{x, yyyy} = u ″″$ , and

ε_{ij} = \frac{u_{i, j} + u_{j, i}}{2} = \frac{1}{2} [\begin{matrix} 0 & u' \\ u' & 0 \end{matrix}] .

(23)

Using equation (22), the governing equation for this problem is formed as

c_{2} u ″ - (c_{5} + c_{6} + c_{7}) u ″″ = 0,

(24)

which is a fourth-order ordinary differential equation, and its general solution is

u = u_{x} (y) = q_{1} + q_{2} y + q_{3} \sinh (\frac{y}{r}) + q_{4} \cosh (\frac{y}{r}),

(25)

where $r = \sqrt{(c_{5} + c_{6} + c_{7}) / c_{2}}$ and $q_{1}, q_{2}, q_{3}, q_{4}$ are four integration constants. Four boundary conditions are necessary for finding $q_{1}, q_{2}, q_{3}, q_{4}$ . In the following, two sets of boundary conditions are discussed.

3.1. Case 1: displacement prescription

In the first case, as shown in Figure 1, the boundary conditions are considered as

u_{x} (0) = 0, u_{x} (H) = \hat{u}, m_{xy} (0) = 0, u' (H) = 0,

(26)

where a Dirichlet boundary condition is assumed on the top edge by applying a displacement $\hat{u}$ . Furthermore, the displacement and double traction, $m_{xy}$ (as defined in equation (11)) are set to zero on the bottom edge. The fourth boundary condition states that the displacement gradient is zero on the top edge. This term activates the strain gradient terms of the theory. In the case of neglecting this boundary condition, the problem of simple shear is obtained as known from the first-order theory. Therefore, this boundary condition is crucial to test the numerical implementation of the strain gradient theory. On the application level, a relatively rigid bar is adhered on top of the plate. In this way, rotating along $y$ is prevented, as visible in Figure 1, which leads to $u' (H) = 0$ .

For the first two boundary conditions in equation (26), we use equation (25), and we obtain

u_{x} (0) = q_{1} + q_{4} = 0,

(27)

u_{x} (H) = q_{1} + q_{2} H + q_{3} \sinh (\frac{H}{r}) + q_{4} \cosh (\frac{H}{r}) = \hat{u} .

(28)

The double traction, $m_{xy}$ , has been defined in equation (11) and is calculated in equation (21). For the current problem, the double traction has only one non-zero term as

\begin{matrix} m_{xy} = & \frac{\partial w}{\partial u_{x, yy}} = (c_{5} + c_{6} + c_{7}) u ″ \\ = (c_{5} + c_{6} + c_{7}) (\frac{q_{3}}{r^{2}} \sinh (\frac{y}{r}) + \frac{q_{4}}{r^{2}} \cosh (\frac{y}{r})) . \end{matrix}

(29)

Therefore, for the third boundary condition in equation (26), we have

m_{xy} (0) = (c_{5} + c_{6} + c_{7}) u ″ (0) = (c_{5} + c_{6} + c_{7}) \frac{q_{4}}{r^{2}} = 0,

(30)

and for the fourth boundary condition, we take the derivative of $u$ in equation (25) with respect to $y$ , and set $y = H$ ,

u' (H) = q_{2} + \frac{q_{3}}{r} \cosh (\frac{H}{r}) + \frac{q_{4}}{r} \sinh (\frac{H}{r}) = 0 .

(31)

By solving equations (27), (28), (30), and (31), as shown in Appendix 2, the four unknowns $(q_{1}, q_{2}, q_{3}, q_{4})$ are acquired

\begin{matrix} q_{1} = 0, q_{2} = \frac{\hat{u} \cosh (\frac{H}{r})}{H \cosh (\frac{H}{r}) - r \sinh (\frac{H}{r})}, \\ q_{3} = \frac{\hat{u} r}{r \sinh (\frac{H}{r}) - H \cosh (\frac{H}{r})}, q_{4} = 0 . \end{matrix}

(32)

3.2. Case 2: applying traction

In the second case, as shown in Figure 2, the boundary conditions are considered as

u_{x} (0) = 0, t_{x} (H) = \hat{t}, u' (0) = 0, m_{xy} (H) = 0,

(33)

where a Neumann boundary condition is assumed on the top edge by a given traction $\hat{t}$ (as defined in equation (10)). The double traction is set to zero on the top edge. Furthermore, the displacement and displacement gradient are set to zero on the bottom edge. The latter condition is responsible for the strain gradient terms of the theory. Indeed, it models a clamped edge as seen in Figure 2.

For the first boundary condition in equation (33), from equation (25), we obtain

u_{x} (0) = q_{1} + q_{4} = 0 .

(34)

The traction, $t_{x}$ with $n = (0, 1)$ , has been defined in equation (10), and for the current problem, it is reduced to

\begin{matrix} t_{x} = & \frac{\partial w}{\partial u_{x, y}} - \frac{\partial w}{\partial u_{x, yyy}} \\ = & c_{2} u' - (c_{5} + c_{6} + c_{7}) u ″' \\ = & c_{2} (q_{2} + \frac{q_{3}}{r} \cosh (\frac{y}{r}) + \frac{q_{4}}{r} \sinh (\frac{y}{r})) \\ - & (c_{5} + c_{6} + c_{7}) (\frac{q_{3}}{r^{3}} \cosh (\frac{y}{r}) + \frac{q_{4}}{r^{3}} \sinh (\frac{y}{r})) \end{matrix}

(35)

Therefore, for the second boundary condition in equation (33), we have

\begin{matrix} t_{x} (H) = c_{2} (q_{2} + \frac{q_{3}}{r} \cosh (\frac{H}{r}) + \frac{q_{4}}{r} \sinh (\frac{H}{r})) \\ - (c_{5} + c_{6} + c_{7}) (\frac{q_{3}}{r^{3}} \cosh (\frac{H}{r}) + \frac{q_{4}}{r^{3}} \sinh (\frac{H}{r})) = \hat{t} . \end{matrix}

(36)

The third and fourth boundary conditions in equation (33) give

m_{xy} (H) = (c_{5} + c_{6} + c_{7}) (\frac{q_{3}}{r^{2}} \sinh (\frac{H}{r}) + \frac{q_{4}}{r^{2}} \cosh (\frac{H}{r})) = 0,

(37)

and

u' (0) = q_{2} + \frac{q_{3}}{r} = 0,

(38)

respectively.

By solving equations (34), (36), (37), and (38), as shown in Appendix 3, the four unknowns $(q_{1}, q_{2}, q_{3}, q_{4})$ are obtained as

\begin{matrix} q_{1} = \frac{- r^{3} \hat{t} \sinh (\frac{H}{r})}{s}, q_{2} = \frac{r^{2} \hat{t} \cosh (\frac{H}{r})}{s}, \\ q_{3} = \frac{- r^{3} \hat{t} \cosh (\frac{H}{r})}{s}, q_{4} = \frac{r^{3} \hat{t} \sinh (\frac{H}{r})}{s}, \end{matrix}

(39)

where

s = c_{2} r^{2} (\cosh (\frac{H}{r}) - 1) + c_{5} + c_{6} + c_{7} .

(40)

4. Numerical implementation

For the numerical implementation of the strain gradient theory, we use the same simple shear problem as in the analytical solution. Two approaches are discussed in the following: (1) FEM and (2) Mixed FEM. The modeling, implementation, and post-processing steps are all carried out using open-source packages. For the FEM analysis, we utilize the FEniCS libraries (73) by following the computational framework as in Abali (74). FEniCS is a package of codes for solving partial differential equations, and it supports symbolic differentiation, which is exploited herein. The FEM code is developed in Python.

The plate under simple shear is modeled as a two-dimensional plate with the height $H = 0.5$ mm. In the problem description, the length of the plate is assumed to be infinite in the $x$ -direction. In the numerical models, we set a finite value for the length (three times of the height of the plate). Periodic boundary conditions are applied on the lateral edges of the plates; hence, the results will not depend on the value that is chosen for the plate’s length. The periodic boundary conditions applied on the left and right edges of the plate impose the following conditions on the domain: (1) the left and right edges of the plate are of course equal in length, (2) the nodes of these two edges have the same vertical coordinates such that there are corresponding nodes on both edges, and (3) the nodal value at the left is restricted to be the same as the corresponding nodal value at the right.

For the material of the plate, we assume Young’s modulus, $E = 400$ MPa and Poisson’s ratio, $ν = 0.49$ . The choice of Poisson’s ratio is crucial for consistency with the semi-analytical ansatz. The simple shear problem is valid under the constant volume assumption that is equivalent to a Poisson’s ratio of 0.5, but here it is not possible to use 0.5 value since it makes the first Lame parameter $(c_{1})$ to be infinite; therefore, we use 0.49 for Poisson’s ratio. The constitutive parameters $c_{1}$ and $c_{2}$ read

\begin{matrix} c_{1} = λ = \frac{E ν}{(1 + v) (1 - 2 ν)} = 6577.18 MPa \\ c_{2} = μ = \frac{E}{2 (1 + ν)} = 134.23 MPa . \end{matrix}

(41)

Moreover, the material parameters for strain gradient elasticity need to be chosen. For this purpose, herein, we utilize the granular micromechanical modeling [75 –78]. Based on the formulation derived in Barchiesi et al. [79], we have

\begin{matrix} c_{3} = c_{4} = \frac{ℓ^{2}}{112} λ \\ c_{5} = c_{7} = \frac{ℓ^{2}}{1120} (7 μ + 3 λ) \\ c_{6} = \frac{ℓ^{2}}{1120} (7 μ - 4 λ), \end{matrix}

(42)

where $ℓ$ is the characteristic length of the granular material of the substructure. We assume three values for $ℓ$ as ${0.1, 0.2, 0.3}$ mm, which lead to material parameters as compiled in Table 1. The negative values in Table 1 may raise concern about the positive definiteness of the strain energy function. However, the individual values of the material parameters are not of importance. The energy value must be indeed positive in order to have a unique solution. We refer to previous studies [66,70,80] for discussion of positive definiteness in the strain gradient theory. In the problem Case 1, we set $\hat{u} = 0.05$ mm and in the problem Case 2, we set $\hat{t} = 1.0$ MPa.

Table 1.

Strain gradient elasticity material parameters for constructing equation (15) calculated for different values of the characteristic length $ℓ$ .

$ℓ$ (mm)	$c_{3}$ (N)	$c_{4}$ (N)	$c_{5}$ (N)	$c_{6}$ (N)	$c_{7}$ (N)
0.1	0.59	0.59	0.18	−0.23	0.18
0.2	2.35	2.35	0.74	−0.91	0.74
0.3	5.29	5.29	1.66	−2.04	1.66

4.1. FEM

In the classical elasticity theory, we deal with a second-order partial differential equation, while in the strain gradient elasticity theory, we have a governing equation of fourth order [81]. As a result, for satisfying the $H^{2}$ regularity in the weak form, unknown fields need to belong to $C^{1}$ within the whole domain. One possible approach is to set the minimum regularity requirement for the shape functions to be $C^{1}$ [82,83]. Another implementation of the strain gradient theory with $C^{1}$ regularity is presented in Glüge [51].

For obtaining the weak form for implementing in the FEM code, we begin with the governing equation of strain gradient elasticity as derived in equation (7). The unknown displacement $u_{i}$ is represented using nodal values discretely in space and form functions for interpolation between nodal values. We circumvent using different notations for analytical functions and their discrete representations, since they never appear in the same equation. We use a discretization using (triangle) Lagrange elements, which generates piecewise continuous polynomials such that they are adequate for approximation in Hilbertian Sobolev space $H^{1}$ . This standard FEM elements of order $q$ spans $P_{q}$ on triangles, $T$ , in a two-dimensional continuum. The computational domain, $Ω$ , is discretized by dividing it in triangles, and this triangulation is denoted by $T$ . Therefore, we use a vector space for displacements

\begin{matrix} V_{FEM} = {{{u_{i}} \in H^{1} (Ω) : {u_{i}} |}_{T} \in P_{2} (T) \forall T \in T} . \end{matrix}

(43)

The Lagrange elements are $C^{0}$ continuous across element boundaries; however, it is necessary to have $C^{1}$ continuity everywhere in the domain. Instead of using more advanced elements [84], we add a term to the weak form to enforce the continuity of the first gradient of displacement across the elements. By adding such a term to equation (7), we obtain

\begin{matrix} - \int_{Ω} (\frac{\partial w}{\partial u_{i, j}} δ u_{i, j} + \frac{\partial w}{\partial u_{i, jk}} δ u_{i, jk}) d V + \int_{\partial Ω} (t_{i} δ u_{i} + m_{ij} δ u_{i, j}) d A \\ - Q \int_{Γ} (u_{i, j}^{+} - u_{i, j}^{-}) δ u_{i, j}^{+} d I = 0, \end{matrix}

(44)

where

Q = \sqrt{c_{2} (c_{5} + c_{6} + c_{7})},

(45)

and we have inserted $\frac{\partial W_{s}}{\partial u_{i}} = t_{i}$ and $\frac{\partial W_{s}}{\partial u_{i, j}} = m_{ij}$ from equations (10) and (11), respectively. The term $Q$ is the penalty parameter; it is chosen arbitrarily but in connection with the material parameters and in the adequate unit in a way that all the terms in equation (44) have the same unit. Since the formulation involves second gradient in space, we choose $q = 2$ such that $P_{2}$ elements are utilized.

In equation (44), $u_{i, j}^{+}$ and $u_{i, j}^{-}$ are the displacement gradients on the two sides of an interior facet, which is between pairs of adjacent cells of the mesh. Moreover, the integration in $d I$ is done on the set $Γ$ of the two-sided facet elements while $d A$ is for the one-sided exterior facets of the mesh, i.e., belonging the boundaries of the domain. Figure 3 shows a schematic of a triangular mesh on a domain $Ω$ , its boundary $\partial Ω$ , a cell element $d V$ , a one-sided exterior facet element $d A$ , and a two-sided facet (edge) element $d I$ . The formulation holds in two-dimensional or three-dimensional domains; however, for the sake of the visualization, we show the idea on a two-dimensional mesh in Figure 3. In a two-dimensional domain, which is the case herein, the cell element $d V$ is a surface and the facets $d A$ and $d I$ are edge elements.

Figure 3.

Schematic of a triangular mesh on a domain $Ω$ , its boundary $\partial Ω$ , a cell (surface) element $d V$ , a one-sided exterior facet (edge) element $d A$ , and a two-sided facet (edge) element $d I$ .

It is possible to add a displacement prescription in the numerical computations as a Dirichlet boundary condition. However, in the problems that are discussed in the analytical solution (Section 3), there are boundary conditions of “zero-displacement gradient” on top (Case 1) and bottom (Case 2) edges of the rectangular domain. It is not possible to add such a boundary condition directly in the numerical computations. Hence, we add the zero-displacement gradient condition into the weak form. We update the weak form of equation (44) as

\begin{matrix} - \int_{Ω} (\frac{\partial w}{\partial u_{i, j}} δ u_{i, j} + \frac{\partial w}{\partial u_{i, jk}} δ u_{i, jk}) d V + \int_{\partial Ω} (t_{i} δ u_{i} + m_{ij} δ u_{i, j}) d A \\ - Q \int_{Γ} (u_{i, j}^{+} - u_{i, j}^{-}) δ u_{i, j}^{+} d I + pQ \int_{\partial Ω} u_{i, j} δ u_{i, j} d A = 0, \end{matrix}

(46)

where $p = 10^{6}$ is a penalty factor. Here, the last integral of equation (46) enforces the gradient of displacement to be zero on the intended boundary (top or bottom edge of the plate). For the problem of Case 1 (Section 3.1), we set $t_{i} = 0$ and $m_{ij} = 0$ in the weak form (46) as there are no traction or double traction applied on the boundaries. Furthermore, the last integral of the weak form is calculated on the top edge of the domain $(y = H)$ . For the problem of Case 2 (Section 3.2), we set $t_{i} = \hat{t}$ and $m_{ij} = 0$ in the weak form (46) as the traction $\hat{t}$ is applied on the top edge of the domain and there is no double traction applied on the boundaries. Moreover, the last integral of the weak form is calculated on the bottom edge of the domain $(y = 0)$ .

4.2. Mixed FEM

The Mixed FEM is developed using the Hu–Washizu principle [85]. In the Mixed FEM, more than one function space is used for approximating the variables. Each function space is dedicated to one of the variables, and all the variables are solved simultaneously. In other words, the additional variables are calculated independently together with the primitive variable.

Herein, we define three function spaces in the Mixed FEM formulation, which involves ${u_{i}, g_{ij}, M_{ij}}$ . The first one (Space 1) is a vector space for the primitive variable of the problem, the displacement $u_{i}$ , for which we choose $q = 2$ and $P_{2}$ (quadratic) shape functions. We use a distinct notation for the space derivatives in order to clearly approximate the space derivative of displacement. As the second space in the mixed formulation (Space 2), we introduce $g_{ij} = u_{i, j}$ as a variable to be identified in the computations. We use $P_{1}$ (linear) shape functions for the unknown $g_{ij}$ . For the spaces of displacement and displacement gradient, we use standard polynomial FEM elements in FEniCS.

Then, we impose $u_{i, j} = g_{ij}$ in the computations. This condition ensures that the approximated $g_{ij}$ in Space 2 is equal to the gradient of the approximated displacement in Space 1. For enforcing such an identity condition, as we will explain in the following, we define a tensor space of Lagrange multipliers, $M_{ij}$ . We choose $L^{2}$ elements for $M_{ij}$ with $P_{0}$ such that they are simply constants within elements with a jump across boundaries. We choose the family of discontinuous Galerkin elements in FEniCS which are suitable for this purpose. We set the degree of element to zero as we need constant functions.

At the end, we have the three unknown functions ${u_{i}, g_{ij}, M_{ij}}$ . For our two-dimensional problem, in each node, there are 2 unknowns for the components of $u_{i}$ , 4 unknowns for $g_{ij}$ , and 4 unknowns for $M_{ij}$ , i.e., 10 unknowns in total. These 10 unknowns are all independent in the formulation. In a nutshell, we construct the space for the Mixed FEM formulation as

\begin{matrix} V_{Mixed FEM} = {{{u_{i}} \in H^{1} (Ω) : {u_{i}} |}_{T} \in P_{2} (T) \forall T \in T, \\ \land {g_{ij}} \in H^{1} (Ω) : {g_{ij}} \in P_{1} \forall T \in T, \\ \land {M_{ij}} \in L^{2} (Ω) : {M_{ij}} \in P_{0} \forall T \in T}, \end{matrix}

(47)

within the domain of the continuum body, $Ω$ , with its closure, $\partial Ω$ , where boundary values are given on Dirichlet boundaries, $\partial Ω_{D}$ .

In the Mixed FEM formulation, we create one single mesh in the domain and each of the function spaces is constructed on the same mesh. Depending on the dimension, degree of shape function, and element type, each space has a different total number of degrees of freedom in the domain. The quadratic vector space of $u_{i}$ has 21,960, the linear tensor space of $g_{ij}$ has 11,160, and the scalar tensor space of $M_{ij}$ has 21,600 degrees of freedom. In total, the mixed space has 54,720 degrees of freedom. The test functions, $δ u_{i}, δ g_{ij}, δ M_{ij}$ , are chosen from the same space as the $u_{i}, g_{ij}, M_{ij}$ , respectively.

For generating the weak form for the Mixed FEM approach, we begin with the governing equation of equation (7). From equations (18) and (19), we have

\frac{\partial w}{\partial u_{i, j}} = C_{ijkl} ε_{kl}, \frac{\partial w}{\partial u_{i, jk}} = D_{ijklmn} ε_{lm, n} .

(48)

By inserting $\frac{\partial w}{\partial u_{i, j}}$ and $\frac{\partial w}{\partial u_{i, jk}}$ in equation (7), and also setting $\frac{\partial W_{s}}{\partial u_{i}} = t_{i}$ and $\frac{\partial W_{s}}{\partial u_{i, j}} = m_{ij}$ from equations (10) and (11), respectively, we obtain

\int_{Ω} (C_{ijkl} ε_{kl} δ u_{i, j} + D_{ijklmn} ε_{lm, n} δ g_{ij, k}) d V - \int_{\partial Ω} (t_{i} δ u_{i} + m_{ij} δ u_{i, j}) d A = 0,

(49)

where we have replaced the $δ u_{i, jk}$ with $δ g_{ij, k}$ since it is the test function related to the second gradient of displacement. For imposing the identity $u_{i, j} = g_{ij}$ (the equality of the approximated $g_{ij}$ in Space 2 and the gradient of the approximated displacement in Space 1), we define a residual as

R = \int_{Ω} M_{ij} (g_{ij} - u_{i, j}) d V .

(50)

Taking the variation of $R$ gives

\begin{matrix} δ R = \frac{\partial R}{\partial u_{i}} δ u_{i} + \frac{\partial R}{\partial u_{i, j}} δ u_{i, j} + \frac{\partial R}{\partial g_{ij}} δ g_{ij} + \frac{\partial R}{\partial M_{ij}} δ M_{ij} \\ = \int_{Ω} (M_{ij} (δ g_{ij} - δ u_{i, j}) + (g_{ij} - u_{i, j}) δ M_{ij}) d V, \end{matrix}

(51)

which should vanish in the domain $Ω$ as we aim at minimizing the residual $R$ . The complete weak form of the Mixed FEM approach is generated by adding the variation of $R$ from equation (51) to equation (49) as

\begin{matrix} \int_{Ω} (C_{ijkl} ε_{kl} δ u_{i, j} + D_{ijklmn} ε_{lm, n} δ g_{ij, k} \\ + M_{ij} (δ g_{ij} - δ u_{i, j}) + (g_{ij} - u_{i, j}) δ M_{ij}) d V \\ - \int_{\partial Ω} (t_{i} δ u_{i} + m_{ij} δ u_{i, j}) d A = 0 . \end{matrix}

(52)

For the problem of Case 1 (Section 3.1), we set $t_{i} = 0$ and $m_{ij} = 0$ in the weak form (52) as there are no traction or double traction applied on the boundaries. For the problem of Case 2 (Section 3.2), we set $t_{i} = \hat{t}$ and $m_{ij} = 0$ in the weak form (52) as the traction $\hat{t}$ is applied on the top edge of the domain and there is no double traction applied on the boundaries. The zero gradient boundary conditions of the problem Cases 1 and 2 are applied as Dirichlet boundary conditions.

5. Results and discussion

In this section, the results of the computations are presented and verified by comparing them with the analytical solutions. Furthermore, the convergence analyses of the formulations are reported.

As discussed in Section 4, two numerical implementations are used for the computations: FEM and Mixed FEM. Table 2 summarizes the details of the spaces and the meshes of the two implementations. As shown in Table 2, the total degrees of freedom in the two implementations are of the same order since we want to compare their accuracy.

Table 2.

Details of the spaces and meshes of the FEM and Mixed FEM implementations.

	Space	Variable	Space type	Element type	Shape function degree	Unknowns per node	Degrees of freedom
FEM	–	Displacement $(u_{i})$	Vector	Standard FEM	2	2	55872
Mixed FEM	1	Displacement $(u_{i})$	Vector	Standard FEM	2	2	21960
	2	Displacement gradient $(g_{ij})$	Tensor	Standard FEM	1	4	11160
	3	Lagrange multipliers $(M_{ij})$	Tensor	Discontinuous Galerkin	0	4	21600
	Total	–	–	–	–	10	54720

FEM: finite element method.

5.1. Results of Case 1

Figure 4 shows the deformed state of the plate that has been under prescribed displacement on the top edge from FEM analysis. The effect of the zero-displacement gradient condition on the top edge $(u' (H) = 0)$ is visible in this plot. We see that the lateral edges of the plate are fully vertical in the very beginning of their upper part. In order to compare the numerical and analytical results, we consider the displacement of the right edge of the plate. The displacement $(u_{x})$ of the right edge calculated by FEM and Mixed FEM is plotted along $y$ in Figures 5 and 6, respectively, for three values of the characteristic length, $ℓ$ , which results in different material parameters, $c_{3}, c_{4}, c_{5}, c_{6}, and c_{7}$ . Figures 5 and 6 show that both FEM and Mixed FEM formulations are matching in a very satisfactory way with the analytical solution.

Figure 4.

The deformed state of the plate under prescribed displacement on the top edge (Case 1) from FEM analysis.

Figure 5.

Displacement of the right edge of the plate, under prescribed displacement on the top edge (Case 1), along $y$ : numerical computation (FEM) versus analytical solution for three different values of $ℓ$ : (a) $ℓ = 0.1$ mm, (b) $ℓ = 0.2$ mm, and (c) $ℓ = 0.3$ mm.

Figure 6.

Displacement of the right edge of the plate, under prescribed displacement on the top edge (Case 1), along $y$ : numerical computation (Mixed FEM) vs analytical solution for three different values of $ℓ$ : (a) $ℓ = 0.1$ mm, (b) $ℓ = 0.2$ mm, and (c) $ℓ = 0.3$ mm.

To see the role of strain gradient terms in systems with different sizes, we simulate the problem of Case 1 for three plates using the Mixed FEM formulation. The height of the three plates are $H = {0.2, 2.0, 20.0}$ mm. We set the characteristic length, $ℓ = 0.1$ mm, and prescribed displacement, $\hat{u} = 0.05$ mm, for the three simulations. Figure 7 shows the displacement of right edge along $y$ normalized with respect to plate height. We see, in Figure 7, that the $H = 20.0$ mm case is behaving linearly compared to the smaller heights. The reason is simply the increased ratio of the characteristic length of the geometry, $H$ , by the characteristic length of the material, $ℓ$ . As expected, for higher $H / ℓ$ values, strain gradient terms are less important. In fact, the strain gradient theory is necessary between $H / ℓ$ larger than 1 and smaller than a certain value which can be estimated using the study presented herein.

Figure 7.

Displacement of the right edge of the plate, under prescribed displacement on the top edge (Case 1) from Mixed FEM analysis for different values of plate height, $H$ , and $ℓ = 0.1$ mm.

Convergence analyses are carried out for the FEM and Mixed FEM implementations to ensure the mesh independency of the results. The error is defined as

Error = \int_{Ω} | u_{i}^{comp} - u_{i}^{ana} |^{2} d V

(53)

where $u_{i}^{comp}$ is the displacement from the numerical analysis and $u_{i}^{ana}$ is the displacement from the analytical solution. In Figure 8(a) and (b), the monotonously decreasing trend of the error on the $\log - \log$ plots is seen for FEM and Mixed FEM implementations, respectively.

Figure 8.

Mesh convergence analysis of FEM and Mixed FEM approaches in log-log scale: (a) FEM and (b) Mixed FEM.

5.2. Results of Case 2

In the problem of Case 2, a traction is applied on the top edge of the plate. Figure 9 depicts the deformation of the plate from Mixed FEM analysis. Here, we see the effect of the zero-displacement gradient condition on the bottom edge $(u' (0) = 0)$ . As a result, the lateral edges of the plate are vertical in the very lower part. For Case 2, the FEM formulation does not produce reliable results and it is not verified by the analytical solution. On the contrary, the Mixed FEM formulation yields highly matching results with the analytical solutions. Figure 10 shows the displacement $(u_{x})$ of the right edge calculated by Mixed FEM for three values of the characteristic length, $ℓ$ , which results in different values for the material parameters, $c_{3}, c_{4}, c_{5}, c_{6}, and c_{7}$ .

Figure 9.

The deformed state of the plate under applied traction on the top edge (Case 2) from Mixed FEM analysis (Deformation scale: 25:1).

Figure 10.

Displacement of the right edge of the plate, under applied traction on the top edge (Case 2), along $y$ : numerical computation (Mixed FEM) vs analytical solution for three different values of $ℓ$ . (a) $ℓ = 0.1$ mm, (b) $ℓ = 0.2$ mm, and (c) $ℓ = 0.3$ mm.

6. Conclusion

Two implementations of the strain gradient elasticity theory based on the FEM are verified. The governing equation of the strain gradient elasticity is modeled using FEM and Mixed FEM, in which three distinct function spaces for the different variables are introduced, allowing the simultaneous determination of the different unknowns on the three distinct spaces. For a simple shear problem, the computations are verified with an analytical solution. The Mixed FEM proved to be reliable in the considered boundary conditions (applying traction and displacement) while the FEM approach succeeds in predicting the displacement only in one of the cases of boundary conditions.

In the FEM implementation approach, in every node, there exist 2 degrees of freedom, while in the Mixed FEM, the number of degrees of freedom is 10 per node; however, for FEM computations, a finer mesh is necessary to give acceptable results. In a nutshell, the Mixed FEM formulation proves to be more robust and reliable for solving problems in the framework of the strain gradient elasticity theory. The codes used for the computations are developed under the GNU Public license [86] made publicly available in Abali [87] for encouraging a transparent scientific exchange.

Footnotes

Appendix 1

For deriving equation (8):

(54)

\int_{Ω} {(\frac{\partial w}{\partial u_{i, j}} δ u_{i})}_{, j} d V = \int_{Ω} \frac{\partial w}{\partial u_{i, j}} δ u_{i, j} d V + \int_{Ω} {(\frac{\partial w}{\partial u_{i, j}})}_{, j} δ u_{i} d V

(55)

\int_{Ω} (\frac{\partial w}{\partial u_{i, j}} δ u_{i}) d V = \int_{\partial Ω} n_{j} \frac{\partial w}{\partial u_{i, j}} δ u_{i} d A

Insert equation (55) in equation (54)

(56)

\int_{Ω} \frac{\partial w}{\partial u_{i, j}} δ u_{i, j} d V = \int_{\partial Ω} n_{j} \frac{\partial w}{\partial u_{i, j}} δ u_{i} d A - \int_{Ω} {(\frac{\partial w}{\partial u_{i, j}})}_{, j} δ u_{i} d V

(57)

\int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}} δ u_{i, j})}_{, k} d V = \int_{Ω} \frac{\partial w}{\partial u_{i, jk}} δ u_{i, jk} d V + \int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i, j} d V

(58)

\int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}} δ u_{i, j})}_{, k} d V = \int_{\partial Ω} n_{k} \frac{\partial w}{\partial u_{i, jk}} δ u_{i, j} d A

Insert equation (58) in equation (57)

(59)

\int_{Ω} \frac{\partial w}{\partial u_{i, jk}} δ u_{i, jk} d V = \int_{\partial Ω} n_{k} \frac{\partial w}{\partial u_{i, jk}} δ u_{i, j} d A - \int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i, j} d V

(60)

\int_{Ω} {({(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i})}_{, j} d V = \int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i, j} d V + \int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}})}_{, kj} δ u_{i} d V

(61)

\int_{Ω} {({(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i})}_{, j} d V = \int_{\partial Ω} n_{j} {(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i} d A

Insert equation (61) in equation (60)

(62)

\int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i, j} d V = \int_{\partial Ω} n_{j} {(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i} d A - \int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}})}_{, kj} δ u_{i} d V

Insert equation (62) in equation (59)

(63)

\begin{matrix} \int_{Ω} \frac{\partial w}{\partial u_{i, jk}} δ u_{i, jk} d V = \int_{\partial Ω} n_{k} \frac{\partial w}{\partial u_{i, jk}} δ u_{i, j} d A \\ - \int_{\partial Ω} n_{j} {(\frac{\partial w}{\partial u_{i, jk}})}_{, k} δ u_{i} d A + \int_{Ω} {(\frac{\partial w}{\partial u_{i, jk}})}_{, kj} δ u_{i} d V \end{matrix}

Finally, by inserting equations (56) and (63) in equation (7) and rearranging, we obtain equation (8).

Appendix 2

The parameters $q_{1}, q_{2}, q_{3}, and q_{4}$ in equation (32) are calculated using equations (27), (28), (30), and (31). From equation (30)

(64)

(c_{5} + c_{6} + c_{7}) \frac{q_{4}}{r^{2}} = 0 \Rightarrow q_{4} = 0 .

By inserting $q_{4}$ in equation (27), we have

(65)

q_{1} + q_{4} = 0 \Rightarrow q_{1} = 0 .

Inserting $q_{4}$ in equation (31) gives

(66)

q_{2} + \frac{q_{3}}{r} \cosh (\frac{H}{r}) = 0 \Rightarrow q_{2} = - \frac{q_{3}}{r} \cosh (\frac{H}{r})

and inserting $q_{1}, q_{2},$ and $q_{4}$ in equation (27) gives

(67)

q_{3} = \frac{\hat{u} r}{r \sinh (\frac{H}{r}) - H \cosh (\frac{H}{r})},

hence, from equation (66)

(68)

q_{2} = \frac{\hat{u} \cosh (\frac{H}{r})}{H \cosh (\frac{H}{r}) - r \sinh (\frac{H}{r})} .

Appendix 3

The parameters $q_{1}, q_{2}, q_{3}, and q_{4}$ in equation (39) are calculated using equations (34), (36), (37), and (38) as explained in the following. From equation (37), we find $q_{3}$ in terms of $q_{4}$ as

(69)

q_{3} = - q_{4} \frac{\cosh (\frac{H}{r})}{\sinh (\frac{H}{r})} .

From equation (38), we find $q_{2}$ in terms of $q_{3}$ and consequently in terms of $q_{4}$ as

(70)

q_{2} = \frac{- q_{3}}{r} = \frac{q_{4}}{r} \frac{\cosh (\frac{H}{r})}{\sinh (\frac{H}{r})} .

Now, we insert $q_{3}$ and $q_{2}$ , respectively, from equations (69) and (70) into equation (36) and we get

(71)

\begin{matrix} \frac{q_{4}}{r} c_{2} (\frac{\cosh (\frac{H}{r}) - \overset{2}{\cosh} (\frac{H}{r}) + \overset{2}{\sinh} (\frac{H}{r})}{\sinh (\frac{H}{r})}) \\ + \frac{q_{4}}{r^{3}} (c_{5} + c_{6} + c_{7}) (\frac{\overset{2}{\cosh} (\frac{H}{r}) - \overset{2}{\sinh} (\frac{H}{r})}{\sinh (\frac{H}{r})}) = \hat{t}, \end{matrix}

and using the formula for hyperbolic functions, $\overset{2}{\cosh} (x) - \overset{2}{\sinh} (x) = 1$ , we have

(72)

q_{4} (\frac{r^{2} c_{2} (\cosh (\frac{H}{r}) - 1) + c_{5} + c_{6} + c_{7}}{r^{3} \sinh (\frac{H}{r})}) = \hat{t},

therefore

(73)

q_{4} = \frac{r^{3} \hat{t} \sinh (\frac{H}{r})}{r^{2} c_{2} (\cosh (\frac{H}{r}) - 1) + c_{5} + c_{6} + c_{7}} .

By inserting $q_{4}$ in equations (69) and (70), we reach to the expression of $q_{3}$ and $q_{2}$ as

(74)

q_{3} = \frac{- r^{3} \hat{t} \cosh (\frac{H}{r})}{r^{2} c_{2} (\cosh (\frac{H}{r}) - 1) + c_{5} + c_{6} + c_{7}},

and

(75)

q_{2} = \frac{r^{2} \hat{t} \cosh (\frac{H}{r})}{r^{2} c_{2} (\cosh (\frac{H}{r}) - 1) + c_{5} + c_{6} + c_{7}} .

Equation (34) gives the $q_{1}$ as

(76)

q_{1} = \frac{- r^{3} \hat{t} \sinh (\frac{H}{r})}{r^{2} c_{2} (\cosh (\frac{H}{r}) - 1) + c_{5} + c_{6} + c_{7}} .

Acknowledgements

The work illustrated in this paper is part of N. Shekarchizadeh’s thesis for the PhD course in “Mathematical Models for Engineering, Electromagnetism and Nanosciences” at Sapienza University of Rome.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Navid Shekarchizadeh

Bilen Emek Abali

Alberto Maria Bersani

References

Müller

. The experimental evidence for higher gradient theories. In: Bertram

Forest

(eds) Mechanics of strain gradient materials (CISM International Centre for Mechanical Sciences), vol. 600. Cham: Springer International Publishing, 2020, pp. 1–18.

dell Isola

Andreaus

Placidi

. At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: an underestimated and still topical contribution of gabrio piola. Math Mech Solids 2015; 20(8): 887–928.

dell Isola

Della Corte

Giorgio

. Higher-gradient continua: the legacy of piola, mindlin, sedov and toupin and some future research perspectives. Math Mech Solids 2017; 22(4): 852–872.

Mindlin

Tiersten

. Effects of couple-stresses in linear elasticity. Arch Ration Mech Anal 1962; 11(1): 415–448.

Mindlin

. Second gradient of strain and surface-tension in linear elasticity. Int J Solids Struct 1965; 1(4): 417–438.

Toupin

. Theories of elasticity with couple-stress. Arch Ration Mech Anal 1964; 17(2): 85–112.

Eringen

. Linear theory of micropolar elasticity. J Math Mech 1966; 15: 909–923.

Neff

Ghiba

I-D

Madeo

, et al. A unifying perspective: the relaxed linear micromorphic continuum. Contin Mech Thermodyn 2014; 26(5): 639–681.

Kaiser

Forest

Menzel

. A finite element implementation of the stress gradient theory. Meccanica 2021; 56(5): 1109–1128.

10.

Rosi

Placidi

Auffray

. On the validity range of strain-gradient elasticity: a mixed static-dynamic identification procedure. Eur J Mech A Solids 2018; 69: 179–191.

11.

Mousavi

Reddy

Romanoff

. Analysis of anisotropic gradient elastic shear deformable plates. Acta Mech 2016; 227(12): 3639–3656.

12.

Tran

Niiranen

. A geometrically nonlinear Euler–Bernoulli beam model within strain gradient elasticity with isogeometric analysis and lattice structure applications. Math Mech Complex Syst 2020; 8(4): 345–371.

13.

Khakalo

Niiranen

. Form II of Mindlin’s second strain gradient theory of elasticity with a simplification: for materials and structures from nano-to macro-scales. Eur J Mech A Solids 2018; 71: 292–319.

14.

Jebahi

Forest

. Scalar-based strain gradient plasticity theory to model size-dependent kinematic hardening effects. Contin Mech Thermodyn 2021; 33(4): 1223–1245.

15.

Scherer

Besson

Forest

, et al. Strain gradient crystal plasticity with evolving length scale: application to voided irradiated materials. Eur J Mech A Solids 2019; 77: 103768.

16.

Wulfinghoff

Alipour

Rezaei

, et al. Generalized interface models with damage in gradient plasticity. PAMM 2016; 16(1): 411–412.

17.

van Beers

McShane

Kouznetsova

, et al. Grain boundary interface mechanics in strain gradient crystal plasticity. J Mech Phys Solids 2013; 61(12): 2659–2679.

18.

Khakalo

Laukkanen

. Strain gradient elasto-plasticity model: 3d isogeometric implementation and applications to cellular structures. Comput Methods Appl Mech Eng 2022; 388: 114225.

19.

Amar

Chiricotto

Giacomelli

, et al. Mass-constrained minimization of a one-homogeneous functional arising in strain-gradient plasticity. J Math Anal Appl 2013; 397(1): 381–401.

20.

Placidi

Misra

Barchiesi

. Two-dimensional strain gradient damage modeling: a variational approach. Z Angew Math Phys 2018; 69(3): 56.

21.

Yang

Misra

. Higher-order stress-strain theory for damage modeling implemented in an element-free Galerkin formulation. Comput Model Eng Sci 2010; 64(1): 1–36.

22.

Sharma

Peerlings

RHJ

Shanthraj

, et al. An FFT-based spectral solver for interface decohesion modelling using a gradient damage approach. Comput Mech 2020; 65(4): 925–939.

23.

Geers

Borst

Brekelmans

, et al. Validation and internal length scale determination for a gradient damage model: application to short glass-fibre-reinforced polypropylene. Int J Solids Struct 1999; 36(17): 2557–2583.

24.

Abali

Klunker

Barchiesi

, et al. A novel phase-field approach to brittle damage mechanics of gradient metamaterials combining action formalism and history variable. J Appl Math Mech/Z Angew Math Mech 2021; 101(9): e202000289.

25.

Nguyen

Niiranen

. A second strain gradient damage model with a numerical implementation for quasi-brittle materials with micro-architectures. Math Mech Solids 2020; 25(3): 515–546.

26.

Seppecher

Alibert

J-J

dell’Isola

. Linear elastic trusses leading to continua with exotic mechanical interactions. J Phys: Conf Ser 2011; 319: 012018.

27.

Turco

. How the properties of pantographic elementary lattices determine the properties of pantographic metamaterials. In: Abali

Altenbach

dell’Isola

, et al. (eds) New achievements in continuum mechanics and thermodynamics: a tribute to Wolfgang H. Müller (Advanced structured materials), vol. 108. Cham: Springer International Publishing, 2019, pp. 489–506.

28.

Placidi

Greco

Bucci

, et al. A second gradient formulation for a 2D fabric sheet with inextensible fibres. Z Angew Math Phys 2016; 67(5): 114.

29.

Vangelatos

Yildizdag

Giorgio

, et al. Investigating the mechanical response of microscale pantographic structures fabricated by multiphoton lithography. Extreme Mech Lett 2021; 43: 101202.

30.

Ciallella

Pasquali

Golaszewski

, et al. A rate-independent internal friction to describe the hysteretic behavior of pantographic structures under cyclic loads. Mech Res Commun 2021; 116: 103761.

31.

Spagnuolo

Andreaus

Misra

, et al. Mesoscale modeling and experimental analyses for pantographic cells: effect of hinge deformation. Mech Mater 2021; 160: 103924.

32.

Rahali

Eremeyev

Ganghoffer

. Surface effects of network materials based on strain gradient homogenized media. Math Mech Solids 2020; 25(2): 389–406.

33.

Glaesener

Bastek

J-H

Gonon

, et al. Viscoelastic truss metamaterials as time-dependent generalized continua. J Mech Phys Solids 2021; 156: 104569.

34.

Barchiesi

Eugster

dell’Isola

, et al. Large in-plane elastic deformations of bi-pantographic fabrics: asymptotic homogenization and experimental validation. Math Mech Solids 2020; 25(3): 739–767.

35.

Rosi

Giorgio

Eremeyev

. Propagation of linear compression waves through plane interfacial layers and mass adsorption in second gradient fluids. J Appl Math Mech/Z Angew Math Mech 2013; 93(12): 914–927.

36.

Geers

MGD

Kouznetsova

Brekelmans

WAM

. Gradient-enhanced computational homogenization for the micro-macro scale transition. J Phys IV France 2001; 115: 145–152.

37.

Ganghoffer

Reda

. A variational approach of homogenization of heterogeneous materials towards second gradient continua. Mech Mater 2021; 158: 103743.

38.

Barboura

. Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures. Int J Solids Struct 2018; 136: 60–76.

39.

Yang

Abali

Müller

, et al. Verification of asymptotic homogenization method developed for periodic architected materials in strain gradient continuum. Int J Solids Struct 2022; 238: 111386.

40.

Ivanova

Vilchevskaya

. Micropolar continuum in spatial description. Contin Mech Thermodyn 2016; 28(6): 1759–1780.

41.

Mousavi

Aifantis

. Dislocation-based gradient elastic fracture mechanics for in-plane analysis of cracks. Int J Fract 2016; 202(1): 93–110.

42.

Giorgio

Andreaus

dell’Isola

, et al. Viscous second gradient porous materials for bones reconstructed with bio-resorbable grafts. Extreme Mech Lett 2017; 13: 141–147.

43.

Giorgio

Andreaus

Scerrato

, et al. A viscoporoelastic model of functional adaptation in bones reconstructed with bio-resorbable materials. Biomech Model Mechanobiol 2016; 15(5): 1325–1343.

44.

Giorgio

Spagnuolo

Andreaus

, et al. In-depth gaze at the astonishing mechanical behavior of bone: a review for designing bio-inspired hierarchical metamaterials. Math Mech Solids 2021; 26(7): 1074–1103.

45.

Auffray

Dirrenberger

Rosi

. A complete description of bi-dimensional anisotropic strain-gradient elasticity. Int J Solids Struct 2015; 69–70: 195–206.

46.

Shekarchizadeh

Abali

Barchiesi

, et al. Inverse analysis of metamaterials and parameter determination by means of an automatized optimization problem. J Appl Math Mech/Z Angew Math Mech 2021; 101(8): e202000277.

47.

Shekarchizadeh

Laudato

Manzari

, et al. Parameter identification of a second-gradient model for the description of pantographic structures in dynamic regime. Z Angew Math Phys 2021; 72(6): 190.

48.

Fischer

Klassen

Mergheim

, et al. Isogeometric analysis of 2D gradient elasticity. Comput Mech 2011; 47(3): 325–334.

49.

Makvandi

Reiher

Bertram

, et al. Isogeometric analysis of first and second strain gradient elasticity. Comput Mech 2018; 61(3): 351–363.

50.

Rudraraju

Van der Ven

Garikipati

. Three-dimensional isogeometric solutions to general boundary value problems of toupin’s gradient elasticity theory at finite strains. Comput Methods Appl Mech Eng 2014; 278: 705–728.

51.

Yang

Abali

Müller

. On finite element analysis in generalized mechanics. In: Indeitsev

Krivtsov

(eds) International Summer School-conference “advanced problems in mechanics.” Cham: Springer, 2019, pp. 233–245.

52.

Glüge

. A C1 incompatible mode element formulation for strain gradient elasticity. In: Altenbach

Müller

Abali

(eds) Higher gradient materials and related generalized continua (Advanced structured materials), vol. 120. Cham: Springer International Publishing, 2019, pp. 95–120.

53.

Papanicolopulos

S-A

Zervos

Vardoulakis

. A three-dimensional C1 finite element for gradient elasticity. Int J Numer Methods Eng 2009; 77(10): 1396–1415.

54.

Phunpeng

Baiz

. Mixed finite element formulations for strain-gradient elasticity problems using the FEniCS environment. Finite Elem Anal Des 2015; 96: 23–40.

55.

Shu

King

Fleck

. Finite elements for materials with strain gradient effects. Int J Numer Methods Eng 1999; 44(3): 373–391.

56.

Amanatidou

Aravas

. Mixed finite element formulations of strain-gradient elasticity problems. Comput Methods Appl Mech Eng 2002; 191(15–16): 1723–1751.

57.

Cordero

Forest

Busso

. Second strain gradient elasticity of nano-objects. J Mech Phys Solids 2016; 97: 92–124.

58.

Yang

Timofeev

Abali

, et al. Verification of strain gradient elasticity computation by analytical solutions. J Appl Math Mech/Z Angew Math Mech 2021; 101(12): e202100023.

59.

Rizzi

Khan

Ghiba

I-D

, et al. Analytical solution of the uniaxial extension problem for the relaxed micromorphic continuum and other generalized continua (including full derivations). Arch Appl Mech. Epub ahead of print 17 November 2021. DOI: 10.1007/s00419-021-02064-3.

60.

Rizzi

Hütter

Madeo

, et al. Analytical solutions of the cylindrical bending problem for the relaxed micromorphic continuum and other generalized continua. Contin Mech Thermodyn 2021; 33(4): 1505–1539.

61.

Zervos

Papanicolopulos

S-A

Vardoulakis

. Two finite-element discretizations for gradient elasticity. J Eng Mech 2009; 135(3): 203–213.

62.

Steigmann

dell’Isola

. Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching. Acta Mech Sin 2015; 31(3): 373–382.

63.

dell’Isola

Steigmann

. A two-dimensional gradient-elasticity theory for woven fabrics. J Elast 2015; 118(1): 113–125.

64.

Abali

. Revealing the physical insight of a length scale parameter in metamaterials by exploring the variational formulation. Contin Mech Thermodyn 2018; 31(4): 885–894.

65.

Abali

Müller

dell’Isola

. Theory and computation of higher gradient elasticity theories based on action principles. Arch Appl Mech 2017; 87(9): 1495–1510.

66.

Nazarenko

Glüge

Altenbach

. Positive definiteness in coupled strain gradient elasticity. Contin Mech Thermodyn 2021; 33(3): 713–725.

67.

Nazarenko

Glüge

Altenbach

. Inverse Hooke’s law and complementary strain energy in coupled strain gradient elasticity. J Appl Math Mech/Z Angew Math Mech 2021; 101(9): e202100005.

68.

Abali

Barchiesi

. Additive manufacturing introduced substructure and computational determination of metamaterials parameters by means of the asymptotic homogenization. Contin Mech Thermodyn 2021; 33: 993–1009.

69.

Lazar

. Irreducible decomposition of strain gradient tensor in isotropic strain gradient elasticity. J Appl Math Mech/Z Angew Math Mech 2016; 96(11): 1291–1305.

70.

dell’Isola

Sciarra

Vidoli

. Generalized Hooke’s law for isotropic second gradient materials. Proc R Soc A: Math Phys Eng Sci 2009; 465(2107): 2177–2196.

71.

Abali

Müller

Eremeyev

. Strain gradient elasticity with geometric nonlinearities and its computational evaluation. Mech Adv Mater Mod Process 2015; 1(1): 4.

72.

Mindlin

. Micro-structure in linear elasticity. Arch Ration Mech Anal 1964; 16(1): 51–78.

73.

Logg

Mardal

K-A

Wells

. (eds). Automated solution of differential equations by the finite element method: the FEniCS book (Lecture notes in computational science and engineering), vol. 84. Berlin: Springer, 2012.

74.

Abali

. Computational reality (Advanced structured materials), vol. 55. Singapore: Springer Nature, 2017.

75.

Chang

Misra

. Packing structure and mechanical properties of granulates. J Eng Mech 1990; 116(5): 1077–1093.

76.

Misra

Yang

. Micromechanical model for cohesive materials based upon pseudo-granular structure. Int J Solids Struct 2010; 47(21): 2970–2981.

77.

Yang

Ching

Misra

. Higher-order continuum theory applied to fracture simulation of nanoscale intergranular glassy film. J Nanomech Micromech 2011; 1(2): 60–71.

78.

Yang

Misra

. Micromechanics based second gradient continuum theory for shear band modeling in cohesive granular materials following damage elasticity. Int J Solids Struct 2012; 49(18): 2500–2514.

79.

Barchiesi

Misra

Placidi

, et al. Granular micromechanics-based identification of isotropic strain gradient parameters for elastic geometrically nonlinear deformations. J Appl Math Mech/Z Angew Math Mech 2021; 101: e202100059.

80.

Eremeyev

Lurie

Solyaev

, et al. On the well posedness of static boundary value problem within the linear dilatational strain gradient elasticity. Z Angew Math Phys 2020; 71(6): 182.

81.

Tsepoura

Papargyri-Beskou

Polyzos

, et al. Static and dynamic analysis of a gradient-elastic bar in tension. Arch Appl Mech 2002; 72(6): 483–497.

82.

Niiranen

Khakalo

Balobanov

, et al. Variational formulation and isogeometric analysis for fourth-order boundary value problems of gradient-elastic bar and plane strain/stress problems. Comput Methods Appl Mech Eng 2016; 308: 182–211.

83.

Khakalo

Niiranen

. Isogeometric analysis of higher-order gradient elasticity by user elements of a commercial finite element software. Comput Aided Des 2017; 82: 154–169.

84.

Arnold

Falk

Winther

. Finite element exterior calculus, homological techniques, and applications. Acta Numer 2006; 15: 1–155.

85.

Washizu

. Variational methods in elasticity and plasticity. 3rd ed. New York: Pergamon Press, 1982.

86.

GNU General Public License, 2007, http://www.gnu.org/licenses/gpl.html

87.

Abali

. Supply code for computations, 2021, http://bilenemek.abali.org/