Mixed finite element formulations for strain-gradient elasticity problems using the FEniCS environment

Finite Elements in Analysis and Design 96 (2015) 23–40

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Mixed finite element formulations for strain-gradient elasticity problems using the FEniCS environment V. Phunpeng n, P.M. Baiz Department of Aeronautical Engineering, Imperial College London, SW7 2AZ, United Kingdom

art ic l e i nf o

a b s t r a c t

Article history: Received 4 June 2014 Received in revised form 2 November 2014 Accepted 24 November 2014 Available online 25 December 2014

This paper presents finite element implementations for the solution of gradient elasticity problems using the FEniCS project. The FEniCS Project provides a novel tool for the automated solution of partial differential equations by the finite element method. In particular, it offers a significant flexibility with regards to modelling and numerical discretization choices. The weak form of the gradient elasticity problem is derived from the Principal of Virtual Work. An equivalent mixed-type finite element formulation for the strain gradient elasticity problem is implemented in order to avoid the use of C1 continuous elements. The complete methodology and source codes (python scripts) are provided. Numerical results are presented and compared with well-known benchmark examples, demonstrating the applicability of FEniCS for such applications. & 2014 Elsevier B.V. All rights reserved.

Keywords: Mixed finite element formulation Strain-gradient elasticity FEniCS

1. Introduction Cutting edge composite material technology is becoming more dependent on nanosize particles and/or nanostructural systems. Nano scale systems are now commonly embedded into a standard matrix in order to obtain materials with more attractive properties than its original constituents. Carbon nanotubes (CNTs) are an example of nano-scale structures, for which a great deal of research has been carried out in order to investigate their properties and potential applications. Other nano scale aggregates that are being used in state of the art composite technology are Nanoclays and Graphene. Unfortunately, simulation technology in this area is not developing as fast as material science and manufacturing technology, potentially hampering the future progress of nano engineered materials. To study nanocomposites, classical (local) continuum theories tend to be inadequate, due to the fact that local theories do not account for micro-structural length scales of the material. Traditional or local elasticity theory considers that stresses at a point are a function of the strains at the same point (bonding forces between atoms are not considered). This assumption is valid for most continuum/macro scale simulations but in general it does not hold for micro/nano scales. To study the issue of size effects, theories which consider the material behaviour at a point as a function of the deformation of the surrounding have been proposed, these are usually referred as non-local or strain gradient theories. Non-local theories state that the stress at a point depends on the strain at all points in the continuum body [23,40]. There is a great deal of research dealing with non-local elastic theories [23–26,50] and gradient theories [34,40,61], as these could provide an efficient way to account for important scale effects. To study size effects in micro/nano composite materials, non-local/strain gradient theories currently represent one of the most attractive options as they provide a strong theoretical interpretation that converges to standard behaviours. For example, at the macro-scale, non-local/strain gradient theories seamlessly converge to the standard classical elasticity theory and are also consistent with continuum thermodynamics theory [44–47]. This type of theory can be extended to consider other configurations, e.g. non-singularities in dislocations and cracks [3–6], as well as dispersion and size effects in nano-objects [15,31]. The most widely used and most important computational method in engineering applications is the Finite Element Method (FEM). It has been implemented successfully in engineering problems in many areas such as solid and fluid mechanics [1,22,27,29,49,52], electromagnetics [30,37] and biomechanics [48,51], to mention just a few examples. Most engineering commercially available software packages for structural/material simulations are based in FEM (e.g. Abaqus [55]). Unfortunately very few provide support/capabilities to

n

Corresponding author. E-mail addresses: [email protected] (V. Phunpeng), [email protected] (P.M. Baiz).

http://dx.doi.org/10.1016/j.finel.2014.11.002 0168-874X/& 2014 Elsevier B.V. All rights reserved.

24

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

model the nano/micro scale behaviour of materials (implementations of non-local/strain gradient FEM formulations are more complex than for classical theories). Under certain conditions, there is a direct equivalence between non-local stress theories and strain gradient theories. From the two approaches, strain gradient theories are well known to be simpler to solve/implement. They have previously been considered in several publications, using either the mixed-type finite element formulation [7,53] or meshless methods [17,56]. The main difficulty of using FEM to solve strain gradient elasticity is the continuity requirement (C1 continuity). In this paper, C1 continuous finite element can be avoided by a mixed-type formulation. Establishing mixed formulations was also not an easy task in the past [9]. Fortunately, thanks to the FEniCS project [36], development of complex mixed finite element formulations has become easier. The FEniCS Project is a modern collection of open source software components directed at the automated solution of Partial Differential Equations (PDEs) by FEM [36,57]. The FEniCS environment provides the following key components: (i) a high-level user interface (written in Python) that imitates the mathematical formulation of a finite element discretization, (ii) sophisticated computational techniques that allow a large class of finite element discretizations to be concisely defined while providing highly efficient code [32,33,43]. A very important characteristic of the FEniCS environment is its flexibility for changing the governing equations and their method of discretization, allowing the use of highly unusual function spaces. This flexibility makes FEniCS particularly useful for areas of research in which computational experimentation with different constitutive relations, different auxiliary equations and different regimes are of importance (as in the case of strain gradient elasticity). Its automation focus also allows for the automation of tedious and error-prone tasks such as the assembly of finite element matrices and vectors, encouraging rapid prototyping and development. Overall, the user-specified equations and methods of discretization remain explicit while the generic finite element implementation details are kept under the hood. Advanced FEM implementations in FEniCS are increasingly taking place as the capabilities offered by the project become evident. For example, Abert et al. [2] presented a micromagnetic simulation code, Hoffman et al. [28] presented a framework for adaptive finite element computation of turbulent flow and fluid-structure interaction, Mortensen et al. [42] presented a framework for experiments in computational turbulence while Vynnytska et al. [60] evaluate the usability of the FEniCS Project for mantle convection simulations. The main objectives of this paper are (i) to implement in FEniCS for the first time a set of mixed finite element formulations for the solution of strain gradient elasticity, (ii) to compare the results with the references established in the literature and (iii) to evaluate the suitability of FEniCS for the solution of strain gradient elasticity problems. The first numerical example considered is a patch test in which the displacement functions are quadratic polynomials [63]. Material behaviour of a perfectly bonded bi-material subjected to a shear stress and an infinite plate with a hole subjected to a distribution load are also investigated. Numerical results are presented and compared with well-known benchmark solutions, highlighting the suitability of FEniCS for the solution of this kind of problems.

2. Classical and strain-gradient elasticity 2.1. Classical elasticity In this section, the variational formulation of the problem is introduced. The procedure is derived from the variation of strain energy in a volume V which equals the variation of work done (Principle of Virtual Work). From the domain Ω shown in Fig. 1, the solution of the displacement field ui is investigated through the known boundary conditions; displacement u^ i on the surface dΓ D and traction ti on the surface dΓ N . From the problem domain, the strain energy density U in the volume is expressed as Z Z δ U dΩ ¼ σ ij δεij dΩ ð1Þ Ω

Ω

where σij is the Cauchy stress (a symmetric tensor) and δεij is the arbitrary strain. To find the solution of the displacement field, Eq. (1) has to be rewritten in term of displacements. For this, the kinematic relation is required:

εij ¼ 12 ð∂j ui þ ∂i uj Þ

ð2Þ

Noting that 1

σ ij δεij ¼ σ ij ð∂i δuj þ ∂j δui Þ 2 σ ij σ ij ¼ ∂i δuj þ ∂j δui 2 2 σ ij σ ji ¼ ∂i δuj þ ∂i δuj 2 ¼ σ ij ∂i δuj

2

ð3Þ

Fig. 1. Linear elastic domain problem.

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

25

it is possible to write Z Z δ U dΩ ¼ σ ij ∂i δuj dΩ Ω

Z Ω

ð4Þ

Ω

Using the divergence theorem we can show that Z ∂i ½σ ij δuj  dΩ ¼ ni σ ij δuj dΓ

ð5Þ

Γ

and using the product rule, Z Z ∂i ½σ ij δuj  dΩ ¼ ½ð∂i σ ij Þδuj þ σ ij ð∂i δuj Þ dΩ Ω

Z Ω

ð6Þ

Ω

Rearranging the equations above we obtain Z Z Z Z Z Z Z Z σ ij ∂i δuj dΩ ¼ ni σ ij δuj dΓ  ð∂i σ ij Þδuj dΩ σ ij ∂i δuj dΩ þ ð∂i σ ij Þδuj dΩ ¼ ni σ ij δuj dΓ σ ij ∂i δuj dΩ þ bj δuj dΩ ¼ t j δuj dΓ Γ

Ω

Ω

Ω

Γ

Ω

Ω

Γ

ð7Þ

where bj denotes the body force per unit volume and tj is the traction on surface Γ . Therefore, the final form of the principle of virtual work is Z Z Z σ ij ∂i δuj dΩ þ bj δuj dΩ  t j δuj dΓ ¼ 0 ð8Þ Ω

Ω

Γ

Thus, the equilibrium equation in the body Ω is ∂i σ ij þ bj ¼ 0

ð9Þ

and the boundary condition on surface Γ is t j ¼ ni σ ij :

ð10Þ

2.2. Non-local and strain gradient elasticity In Mindlin's theory of gradient elasticity, the strain energy density is expressed in terms of the macroscopic displacements ui and the microscopic deformation gradient ηijk. The microscopic deformation gradient is formulated in three different forms [40]: Form I is defined as the second gradient of the macroscopic displacement, Form II as first gradient of the macroscopic strain and Form III there are two microscopic deformation effects including a gradient of macroscopic rotation and the second gradient of macroscopic displacement. The principle of virtual work can be written as a function of strain and strain gradient (including rotation gradient) as [35,38,40] W ¼ Wðεij ; ηijk Þ ¼ Wðεij ; η ij ; η ijk Þ

ð11Þ

the second-order deformation gradient ηijk are defined as [7,16,18,35,38,58,59,66], Form I: second gradient of displacement

ηijk ¼ ∂i ∂j uk

ð12Þ

Form II: first gradient of strain

ηijk ¼ ∂i εjk ¼ 12 ð∂i ∂k uj þ ∂i ∂j uk Þ

ð13Þ

Form III gradient of rotation

Θij ¼ 12 ðui;j  uj;i Þ ωi ¼ 12 ð∇  uÞi ¼  12 eijk Θjk η ij ¼ ωj;i η ijk ¼ ∂i εjk ¼ 13 ð∂j ∂k ui þ∂i ∂k uj þ ∂i ∂j uk Þ

ð14Þ

where ∂i is the gradient operator, ui is the displacement field, εij is the strain and ηijk is the strain gradient, Θij is the rotation tensor, wi is the rotation vector and eijk is the alternator matrix or permutation symbol. In this work, Strain-Gradient Elasticity based on Form-I and Form-II formulations are considered. In 1965, Mindlin [39] showed that the stress σij and double stress τijk are conjugated with the strain and the second-order deformation gradient:

σ ij ¼

∂W ; ∂εij

τijk ¼

∂W ∂ηijk

The variational of the total strain energy in a domain Ω can be written as Z δ W dΩ ¼ ðσ ij δεij þ τijk δηijk Þ dΩ

ð15Þ

Z

Ω

Ω

ð16Þ

26

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

Fig. 2. Surface normal vector on the domain.

Using the divergence theorem to reformulate the previous equation to the first variation of the strain energy density, we obtain (see Appendix A) Z Z Z δ W dΩ ¼  ð∂i σ ik  ∂i ∂j τijk Þδuk dΩ þ ½nj ðσ jk  ∂i τijk Þδuk þ ni τijk ∂j δuk  dΓ ð17Þ Ω

Ω

Γ

For the last term of the surface integral, the gradient of displacement variation on the boundary surface ∂j δuk (see Fig. 2) can be decomposed into a surface normal gradient and a surface gradient as [35]   ∂j δuk ¼ ∂j δuk  ðni ∂i δuk Þnj þ ðni ∂i δuk Þnj ¼ δij ∂i δuk ni nj ∂i δuk þ ðni ∂i δuk Þnj ¼ ðδij  ni nj Þ∂i δuk þ ni nj ∂i δuk

ð18Þ

we can also rewrite it as [35] ∂j δuk ¼ Dj δuk þ nj Dδuk

ð19Þ

where ni ¼ ith component of the unit surface normal vector DðÞ ¼ nk ∂ðÞ=∂xk surface normal gradient operator Dj ðÞ ¼ ðδjk  nj nk Þ∂ðÞ=∂xk surface gradient operator.

Z Γ

From Eq. (19), the last term of the surface integral in Eq. (17) becomes [35] Z Z ni τijk ∂j δuk dΓ ¼ ni τijk Dj ∂uk dΓ þ ni τijk nj D∂uk dΓ Γ

Γ

ð20Þ

The boundary surface Γ is bounded by edge Υ . Using Stoke's surface divergence theorem on a surface Γ , we obtain [35] I Z Z Dj ðni τijk δuk Þ dΓ ¼ ni kj τijk δuk dΥ þ ðDp np Þni nj τijk δuk dΓ

ð21Þ

Then surface integral in Eq. (17) becomes [35] Z I Z ½nj ðσ jk  ∂i τijk;i Þ þni nj τijk ðDl nl Þ  Dj ðni τijk Þδuk dΓ þ ni nj τijk Dδuk dΓ þ Σ m Δðni ki τijk Þδuk dΥ

ð22Þ

Γ

Υ

Γ

Γ

Z Ω

Γ

Υ

So that, the final form of the principle of virtual work [53] is I Z Z ½σ ij δεij þ τijk δηijk  dΩ ¼ bk δuk dΩ þ ½t k δuk þ r k Dδuk  dΓ þ Σ m pk δuk dΥ Ω

Γ

where

σ ij εij τijk ηijk bk ui δ ui tk rk pk

stress strain high-order stress strain gradient body force per unit volume displacement arbitrary displacement increment traction per unit area double traction per unit area line load along sharp edge

Υ

ð23Þ

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

27

and the equilibrium equation in the body is thus [35,53,66]: ∂i ðσ ik ∂j τjik Þ  bk ¼ 0

ð24Þ

and the boundary conditions on Γ and along Υ are (as shown in Fig. 1) [53].

 Traction t k ¼ ni ðσ ik  ∂j τjik Þ þ ni nj τijk ðDp np Þ  Dj ðni τijk Þ ð25Þ  Double traction r k ¼ ni nj τijk

ð26Þ

 Line load pk ¼ Δðni ki τijk Þ

ð27Þ

According to the above equations, gradient elasticity leads to a fourth-order PDE in terms of displacements that can be solved for example by a Hermitian finite element implementation (C1 continuity). Another format of gradient elasticity FEM implementations consists of second-order PDEs in terms of displacements and additional variables that can be solved by mixed formulations (e.g. strains and/or stresses could be approximated independently together with displacements [7,16]). Due to the complications involved in developing general C1 elements (these are easy for 1D but not for 2D or 3D), the present work will specifically focus on developing Gradient Finite Element formulations based on a Mixed Variational Approach.

3. Mixed finite element formulations Mixed finite element methods are formulations in which two or more finite element spaces are used to approximate different variables. Not only the primary variable (such as displacements in classical elasticity) but also secondary and/or tertiary variables would be approximated simultaneously (such as strains/stresses and their derivatives). These more complex implementations are sometimes called two-field or three-field mixed finite elements. There are many books and published works available for two-field or three-field mixed formulations in classical/local elasticity [14,12,65]. Establishing mixed formulations was also not an easy task in the past. Fortunately, thanks to the FEniCS project, development of complex mixed finite element formulations has become easier. In the following sections, a mixed formulation together with its computer code implementation will be shown, to highlight the feasibility of FEniCS for such applications. In the present formulations, the displacement gradients (also called relaxed strains, ψij) and higher order stresses (also called Lagrange multipliers, ρjk) together with the displacements (ui) are the unknown variables of the problem (each will be approximated independently). In order to avoid the use of C1 elements, the principle of virtual work will be rewritten. To derive the weak form of the principle of virtual work, kinematic constrains are introduced. 3.1. Kinematic constrains

ψ ij ¼ uj;i

ð28Þ

an the Kinematic relation of arbitrary variations of δu and δψ is

δψ ij ¼ δuj;i :

ð29Þ

3.2. Weak form

Z Ω

Making use of the divergence theorem, the strain energy density Eq. (23) in the domain Ω can be rewritten as [53] Z Z ½σ ij δεij þ τijk δηijk þ τijk;i ðδuk;j  δψ jk Þ dΩ ¼ bk δuk dΩ þ ½t k δψ k þ nj r k δuk;j  dΓ Ω

Γ

ð30Þ

When applying a mixed finite element formulation, we have to enforce the first variational form of the constraints by introducing Lagrange multipliers, defined as

ρjk ¼  τijk;i

ð31Þ

The constraints are enforced by writing Eq. (29) as a variational form, using an arbitrary variation of the Lagrange multipliers, δτijk;i , giving a set of extra conditions [53] Z ðψ jk  uk;j Þδτijk;i dΩ ¼ 0 ðno sum on j and kÞ ð32Þ Ω

Z Ω

After rewriting Eqs. (30) and (32), the modified virtual work becomes Z Z ðσ ij δεij  ρjk δuk;j Þ dΩ ¼ bk δuk dΩ þ t k δψ k dΓ Ω

Γ

28

Z Z

Ω Ω

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

ðτijk δηijk þ ρjk δψ jk Þ dΩ ¼ ðψ jk  uk;j Þδρjk dΩ ¼ 0

Z Γ

nj r k δuk;j dΓ

ðno sum on j and kÞ:

ð33Þ

3.3. Stability condition of mixed approximation When two or more finite element function spaces are used in mixed finite element methods to evaluate separate variables, the mixed function spaces have to be chosen carefully. For mixed finite element formulation based on strain-gradient elasticity, the necessary conditions for convergence derived by Zervos et al. [62] is sN GP Z nψ

ð34Þ

and a sufficient conditions for convergence is nu 4 nψ þ sN GP

if nψ r sN GP o nu

ð35Þ

where nu and nψ are number of degrees of freedom representing displacements (u) and relaxed strains (ψ), respectively. NGP denotes the number of Gauss points, s is the number of the components of ψ (s ¼4 for 2D problems and s¼ 9 for 3D problems). In addition, theoretically, other ways to automate the examination of stability are the inf-sup or LBB condition. Since Babuska (1973) and Brezzi(1974) presented the stability theory for mixed finite element discretizations, many papers have been published considering the appropriate discretization for mixed formulations [8–14,19–21].

4. Computational implementation (patch test) In this section the proposed methodology will be demonstrated while solving the patch test problem. The overall steps of the methodology are define PDE (weak form), define domains, create mesh and define function spaces, apply boundary conditions, perform computation and visualization of meshes and results. Patch tests are well known numerical examples that help to assess the convergence properties of FEM implementations. The proposed strain gradient mixed finite element implementation and discretizations will be tested with a patch test in which the displacement functions are quadratic polynomials [54,63] (traditional patch tests consider constant strains and stresses). The dimension of the patch is 0.24  0.12. The boundary conditions are the displacements ui and relaxed strains ψij at the boundary nodes. The inner nodes are investigated and the results are compared with analytical solutions [41]. 4.1. Problem definition According to Zhao et al. [64] and Soh and Wanji [54] the displacement functions of the patch test should be a quadratic polynomial that satisfies the equilibrium equations. Following analytical solutions of the patch test [54], the displacement functions are expressed as 10

uðx; yÞ ¼ ∑ ϕi ai i¼1

Table 1 Solution of the patch test problem. i

1

2

3

4

5

6

7

8

9

10

ϕ ψ

1 0

0 1

y x

x  νy

 νx y

y=2 x=2

 xy

ðy2 þ νx2 Þ=2  xy

xy

ðx2 þ νy2 Þ=2

ðx2  y2 Þ=2

ðy2  x2 Þ=2 xy

Fig. 3. Mesh of patch test for 3-node element.

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40 10

vðx; yÞ ¼ ∑ ψ i bi

29

ð36Þ

i¼1

where ai and bi are assumed as be equal to i ¼ 1; 2…, and ϕi and ψi are defined as shown in Table 1. For the boundary conditions of the patch test, the nodal displacements and relaxed strains on the boundaries are given by Eq. (36). The material properties are assumed as E ¼1000, ν ¼ 0:25 and nonlocal parameter le ¼0.1. 4.2. FEniCS implementation In this section the FEniCS implementation will be explained. For a detailed documentation of the FEniCS project, refer to the FEniCS manual [57]. First, the Dolfin library is called (DOLFIN is a C þ þ /Python library that functions as the main user interface of FEniCS):

To deal with multi-dimensional arrays and/or matrices, the Numpy library is used:

The mesh was created in ABAQUS [55] (although any other mesh generator could be employed). The input file from ABAQUS (.inp) was used to obtain the xml file called ‘PatchTest01.xml’ (see Appendix B for the complete file).

To define function spaces on the domain, different element types should be specified on a cell (triangular shape in 2D or tetrahedral in 3D). Using a 2D implementation (it is straight forward to consider a 3D implementation in FEniCS), the degrees of freedom for each variable are given as follows (displacements u, relaxed strains ψ and Lagrange multiplier ρ):     ð37Þ u ¼ fu1 u2 g21  Nu 22nu u^ 2nu 1 

Ψ ¼ ψ 11 ψ 12 ψ 21 ψ 22 

 41

     N ψ 44nψ ψ^ 4nΨ 1

ð38Þ

     N ρ 44nρ ρ^ 4nρ 1 ð39Þ  u  ψ   ρ       u ψ ρ where N , N and N are the shape functions, u^ , ψ^ . ρ^ are nodal values (nodal degrees of freedom) and n , n and n are the number of nodes per element. For displacements and relaxed strains, Lagrange polynomials of degree 2 (CG2) and 1 (CG1) will be used respectively (these basis belong to the H1 Sobolev space). For Lagrange multipliers, discontinuous Lagrangian elements of degree 0 (DG0) are employed (L2 Sobolev space). The graphical representation of these basis functions is given in Fig. 4. The total degrees of freedom for a triangular cell of the discretization described above are as follows:    1 2 3 4 5 6 1 2 3 4 5 6 u^ ¼ u1 u1 u1 u1 u1 u1 u2 u2 u2 u2 u2 u2 ð40Þ

ρ ¼ ρ11 ρ12 ρ21 ρ22



41

   1  ψ^ ¼ ψ 11 ψ 211 ψ 311 ψ 112 ψ 212 ψ 312 ψ 121 ψ 221 ψ 321 ψ 122 ψ 222 ψ 322

ð41Þ

   1 1 1 1  ρ^ ¼ ρ11 ρ12 ρ21 ρ22

ð42Þ

The functional spaces (basis functions) are defined as follows:

Fig. 4. Nodal parameters. (a) Nodal displacements. (b) Nodal relaxed strains. (c) Nodal Lagrange multipliers.

30

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

As part of a symmetric variational weak form, trial and test function are defined using the same spaces,

Using the mesh provided, FEniCS is able to define the entire domain of interest. Boundary definition is also performed using the mesh provided. In the present case (patch test), displacements and strains will be defined along all the boundary, therefore the entire boundary of the patch test will be marked as boundary part 0:

Loading in FEniCS can be defined along the boundary or in the domain, in the present case no actual tractions or body forces are applied (only displacements and relaxed strains are applied), therefore

According to Eq. (36), the boundary conditions of the patch test are defined in FEniCS as follows: uðx; yÞ ¼  4x2 þ 2xy þ 11x=4 þ 9y2 þ6y þ 1 vðx; yÞ ¼ 8x2 þ 2xy þ 4y  29y2 =8 þ2

ψ 11 ðx; yÞ ¼ 2y  8x þ 11=4 ψ 12 ðx; yÞ ¼ 16x þ 2y ψ 21 ðx; yÞ ¼ 2x þ 18y þ 6 ψ 22 ðx; yÞ ¼ 2x  29y=4 þ 4

Material properties and constants definitions in FEniCS are also straight forward:

ð43Þ

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

31

The next step in the implementation methodology involves the definition of the governing equations. First we start with the kinematic relations,

Next, the actual variational formulation is coded as follows (from Eq. (33)):

From the above portion of the script, “a” represents the bilinear form and “L” represents the linear form, which are automatically identified and extracted by FEniCS from the variational formulation. To solve the variational form with the boundary conditions, the bilinear form, linear form and known variables are assembled as a system of equations as follows:

32

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

Finally, result can be visualised in ParaView (export results to .pvd file),

The above simple implementation is robust and general as it will be shown in the following sections. It also makes evident the suitability of FEniCS to solve these complex finite element implementations.

4.3. Numerical results The numerical results and analytical solutions [41] of the patch test are listed in Table 2 (nodal displacements and relaxed strains of node E as shown in Fig. 3). The results show that triangular cells using three different discretizations (CG2CG1DG0, CG2CG2DG0 or CG3CG2DG0) can be used to efficiently interpolate displacements (u), relaxed strains (ψ) and Lagrange multipliers (ρ).

5. Numerical examples In this section, the previous FEniCS implementation (patch test) will be modified in order to solve other well known benchmark problems in the area of non-local and strain gradient elasticity. The modification will only involve domain/boundary definition (extracted from the mesh file) and boundary conditions (displacement and loading defined as described by the problem of interest).

Table 2 Numerical results of the patch test at point E (x¼ 0.04, y¼ 0.02). Results

u

v

∂u=∂x

∂v=∂x

∂u=∂y

∂v=∂y

Exact [41] Zhao et al. [64] Soh and Wanji [54] Present work (using CG2CG1DG0) Present work (using CG2CG2DG0) Present work (using CG3CG2DG0)

1.2288 1.2288 1.2288 1.2288 1.2288 1.2288

2.0929 2.0929 2.0929 2.0929 2.0929 2.0929

2.4700 2.4700 2.4700 2.4700 2.4700 2.4700

0.6800 0.6800 0.6800 0.6800 0.6800 0.6800

6.4400 6.4400 6.4400 6.4400 6.4400 6.4400

3.9350 3.9350 3.9350 3.9350 3.9350 3.9350

Fig. 5. A bimaterial under uniform shear and sketch of the mesh.

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

33

5.1. Boundary layer problem Two materials bonded perfectly are subjected to a shear stress σ 1 21 as shown in Fig. 5. The shear modulus and non-local parameter of the two materials are assumed as [53,56]

μ1 ¼ 2μ2 and l1 ¼ l2 ¼ l

ð44Þ

First, the domain is defined via a 2D mesh of triangular cells. In this problem, it is also necessary to define the different material properties for different parts of the domain. The mesh and the material properties were created in ABAQUS [55]. The input file from ABAQUS (.inp) was used to obtain the xml files for mesh and material properties (these are called ‘10  10meshI.xml’, ‘10  10meshI_E.xml’ and ‘10  10meshI_nu.xml’). For the material properties, a ‘MeshFunction’ over cells is used to represent subdomains of materials.

To create functions of material properties (E and nu), a function space is implemented to collect all elemental properties (discrete function). The material properties in the xml files were read from the mesh functions and are turned to be functions of material properties (E and nu):

When all other material properties are defined, the functions of E and nu are taken into account:

Contrary to the patch test problem, the present case involves the application of different boundary conditions in different parts of the surface, therefore lets define different regions of the boundary:

34

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

In order to avoid rigid body motion, it is also necessary to define points to fix the domain. In this case, the bottom left corner and the bottom right corner are used:

Traction, displacement and relaxed strain boundary conditions are now defined. From Fig. 5, at the bottom left corner we set u ¼ v ¼ 0 and at the bottom right corner u ¼0. The relaxed strains φ11 and φ22 are zero along the left and the right of the boundary.

Fig. 6. Numerical results of boundary layer analysis. (a) Strain ϵ12 using CG2CG1DG0. (b) The comparison of the normalized strain on various discreatizations.

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

35

Results for this problem are provided in Fig. 6 and Table 3, from where it is clear that the agreement between FEniCS solutions and analytical solutions is excellent.

5.2. Stress concentration problem The final example considered is an infinite plate with a circular hole as shown in Fig. 7. The plate is subjected to a uniformly distributed load p, and in this case we will study the influence of size effect in the solution (stress concentration around the hole). The results are compared with analytical and numerical results of published works [41,53,54,64]. First, the domain of the problem is created in Abaqus. The geometry and mesh are saved as an input file (.inp) and then these are converted to the necessary xml file (mesh file). Different meshes with different mesh ratios were considered to investigate the convergence. In the present example, due to symmetric considerations only a quarter of the model will be modelled:

Table 3 Numerical results of the boundary layer analysis at x1 ¼ 0 and l ¼ 0.1. x2 =l

10.9898 9.5309 8.1667 6.8909 5.6979 4.5822 3.5389 2.5633 1.6510 0.7978 0.0000  0.7978  1.6510  2.5633  3.5389  4.5822  5.6979  6.8909  8.1667  9.5309  10.9898

Analytical solutions [53]

CG2CG1DG0

CG2CG2DG0

CG3CG2DG0

ϵ12 =ϵ 12

ϵ12 =ϵ 12 (% error)

ϵ12 =ϵ 12 (% error)

ϵ12 =ϵ 12 (% error)

0.6667 0.6667 0.6666 0.6664 0.6659 0.6644 0.6602 0.6495 0.6240 0.5666 0.4444 0.3834 0.3547 0.3419 0.3366 0.3345 0.3337 0.3334 0.3334 0.3333 0.3333

0.6669 (0.04) 0.6667 (0.00) 0.6668 (0.04) 0.6664 (0.00) 0.6660 (0.01) 0.6639 (0.07) 0.6591 (0.18) 0.6463 (0.49) 0.6163 (1.24) 0.5628 (0.68) 0.4523 (1.76) 0.3917 (2.18) 0.3589 (1.20) 0.3435 (0.47) 0.3372 (0.20) 0.3345 (0.01) 0.3338 (0.02) 0.3333 (0.05) 0.3333 (0.01) 0.3332 (0.04) 0.3333 (0.01)

0.6665 (0.02) 0.6665 (0.02) 0.6665 (0.02) 0.6663 (0.02) 0.6657 (0.04) 0.6638 (0.08) 0.6589 (0.20) 0.6466 (0.46) 0.6184 (0.91) 0.5605 (1.08) 0.4583 (3.11) 0.3910 (1.98) 0.3582 (1.00) 0.3434 (0.44) 0.3371 (0.17) 0.3347 (0.05) 0.3337 (0.01) 0.3334 (0.00) 0.3334 (0.00) 0.3334 (0.00) 0.3334 (0.01)

0.6665 (0.02) 0.6665 (0.02) 0.6665 (0.02) 0.6664 (0.01) 0.6659 (0.00) 0.6646 (0.03) 0.6609 (0.10) 0.6510 (0.22) 0.6273 (0.52) 0.5715 (0.87) 0.4499 (1.23) 0.3796 (0.99) 0.3535 (0.34) 0.3411 (0.24) 0.3362 (0.12) 0.3343 (0.05) 0.3336 (0.02) 0.3334 (0.00) 0.3334 (0.01) 0.3334 (0.01) 0.3334 (0.02)

Fig. 7. The square plate with a circular hole and mesh for the quadrant.

36

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

Boundaries of the quarter model will be now marked (defined) for application of boundary conditions,

Fig. 8. Convergence of the stress concentration using CG2CG1DG0 (ν ¼ 0:0 and l ¼ 0.01).

Fig. 9. Convergence of the stress concentration using CG2CG1DG0 (ν ¼ 0:0 and l ¼0.1).

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

37

Finally, loading and boundary conditions are defined:

In the present example as there is only one material definition, the rest of the script follows exactly the same steps as the patch test example. Results are presented in Figs. 8 and 9 and Tables 4 and 5. The comparison of the stress concentration factor between the analytical solution [41], the two-dimensional simulations of quadrilateral elements [53,54,64] and the three-dimensional simulations [66] are presented in Table 6. The results show that the internal length scale l affects the stress concentration factor significantly and the stress concentration factors obtained with CG2CG1DG0 elements are in excellent agreement with the analytical solutions (Fig. 10).

Table 4 Convergence of stress concentration factor (ν ¼ 0:0) based on strain-gradient formulation II using the CG2CG1DG0 elements and the mesh ratio of 60. Meshes

No. of nodes

No. of elements

Analytical [41]

Present work (Form II)

Error (%)

30  30h1A 30  30h1B 30  30h1C 30  30h1D 30  30h1E 30  30h1F

121 441 961 1681 2601 3721

200 800 1800 3200 5000 7200

2.9984 2.9984 2.9984 2.9984 2.9984 2.9984

3.0991 3.0348 3.0184 3.0115 3.0079 3.0056

3.66 1.21 0.67 0.44 0.32 0.24

Table 5 Convergence of stress concentration factor (ν ¼ 0:0) based on Strain-gradient formulation I using the CG2CG1DG0 elements and the mesh ratio of 60. Meshes

No. of nodes

No. of elements

Analytical [41]

Present work (Form I)

Error (%)

30  30h1A 30  30h1B 30  30h1C 30  30h1D 30  30h1E 30  30h1F

121 441 961 1681 2601 3721

200 800 1800 3200 5000 7200

2.9984 2.9984 2.9984 2.9984 2.9984 2.9984

3.1013 3.0384 3.0230 3.0169 3.0139 3.0122

3.43 1.33 0.82 0.62 0.52 0.46

Table 6 Comparison of stress concentration factor (ν ¼ 0:0) with different non-local parameters with analytical solutions [41]. l/a

Analytical [41]

Numerical [53,54,64]

Zybell [66]

Present work (Form II)

Present work (Form I)

0.01 0.1 1.0

2.9984 2.8779 1.8888

3.003 (0.15%) 2.902 (0.84%) 1.912 (1.23%)

3.010 (0.39%) 2.889 (0.38%) 1.901 (0.65%)

3.0056 (0.24%) 2.8855 (0.26%) 1.8837 (0.27%)

3.0122 (0.46%) 2.8898 (0.41%) 1.8939 (0.27%)

Fig. 10. Meshes for an infinite plate and hole. (a) 30  30h1A. (b) 30  30h1B. (c) 30  30h1C. (d) 30  30h1D.

38

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

6. Conclusion This paper has demonstrated the capabilities of the FEniCS Project to solve mixed finite element formulation for the study of strain gradient elasticity. The numerical results using 2D triangular cells show that discretizations CG2CG1DG0, CG2CG2DG0 and CG3CG2DG0 are in excellent agreement with analytical solutions and other numerical solutions. As expected, the higher the degree of the discretization the better the solution agreement (but at a higher computational cost). As demonstrated in the examples, FEniCS capabilities enabled easy implementation of the variational problems considered. The rich form language provided is helpful in this regard as it allows complex functionals to be quickly formulated and computed. It is expected that the scripts provided in this paper will help the scientific community to develop more general and advanced mechanical models for accurate simulations of nano/micro engineered composites.

Acknowledgements We would like to acknowledge the support of the Thai Government scholarship programme for funding, without which this assignment would not have been possible. We would also like to extend our sincerest thanks to Dr. Jack Hale for his advice and fruitful discussions.

Appendix A A.1. Derivation of strain energy density Using the divergence theorem, the variation of the total strain energy Eq. (16) can be reformulated. The first term was already shown in the strain energy density of classical elasticity (Eqs. (1)–(7)), while the high order term can be rewritten as Z Z τijk δηijk dΩ ¼ τijk ∂i ∂j δuk dΩ ð45Þ Ω

Z Ω

Z Ω

Z Ω

Z Ω

Ω

Using the divergence theorem it is possible to write Z ∂i ½τ ijk ∂j δuk  dΩ ¼ ni τijk ∂j δuk dΓ

ð46Þ

It is also possible to write (using the product rule), Z ∂i ½τ ijk ∂j δuk  dΩ ¼ ð½∂i τijk ∂j δuk þ τijk ∂i ∂j δuk Þ dΩ

ð47Þ

Rearranging Eqs. (46) and (47), we can obtain Z Z τijk δηijk dΩ ¼ ni τijk ∂j δuk dΓ  ð∂i τijk Þ∂j δuk dΩ

ð48Þ

It is necessary to reapply the procedure to the last term of Eq. (48). Using the divergence theorem, Z ∂j ½ð∂i τijk Þδuk  dΩ ¼ nj ð∂i τijk Þδuk dΓ

ð49Þ

Γ

Ω

Γ

Ω

Γ

and the product rule, Z Z ∂j ½ð∂i τijk Þδuk  dΩ ¼ ½ð∂i ∂j τijk Þδuk þ ð∂i τijk Þ∂j δuk  dΩ Ω

Z Ω

Z

ð50Þ

Ω

Rearranging Eqs. (49) and (50) we obtain Z Z ð∂i τijk Þ∂j δuk dΩ ¼ nj ð∂i τijk Þδuk dΓ  ð∂i ∂j τijk Þδuk dΩ Γ

Ω

Therefore, it possible to write Z  Z Z τijk δηijk dΩ ¼ ni τijk ∂j δuk dΓ  nj ð∂i τijk Þδuk dΓ  ð∂i ∂j τijk Þδuk dΩ Ω ZΓ Z Γ Z Ω ¼ ni τijk ∂j δuk dΓ  nj ð∂i τijk Þδuk dΓ þ ð∂i ∂j τijk Þδuk dΩ Γ

Γ

Ω

Now we can substitute Eqs. (7) and (52) into the total strain energy density Eq. (16) to obtain, Z δ W dΩ ¼ ðσ ij δεij þ τijk δηijk Þ dΩ Ω ZΩ Z ¼ ni σ ij δuj dΓ  ð∂i σ ij Þδuj dΩ Γ Ω Z Z Z þ ni τijk ∂j δuk dΓ  nj ð∂i τijk Þδuk dΓ þ ð∂i ∂j τijk Þδuk dΩ Ω ZΓ ZΓ ¼  ð∂i σ ij Þδuj dΩ þ ð∂i ∂j τijk Þδuk dΩ Z

Ω

Ω

ð51Þ

ð52Þ

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

Z

Z ni τijk ∂j δuk dΓ  nj ð∂i τijk Þδuk dΓ ΓZ Γ ZΓ ¼  ð∂i σ ik Þδuk dΩ þ ð∂i ∂j τijk Þδuk dΩ ZΩ Z Ω Z þ nj σ jk δuk dΓ þ ni τijk ∂j δuk dΓ  nj ð∂i τijk Þδuk dΓ ni σ ij δuj dΓ þ

þ

39

Z

Γ

Γ

Γ

Finally we obtain Z δ W dΩ ¼  ð∂i σ ik  ∂i ∂j τijk Þδuk dΩ Ω ZΩ þ ½nj ðσ jk  ∂i τijk Þδuk þ ni τijk ∂j δuk  dΓ :

ð53Þ

Z

Γ

ð54Þ

Appendix B B.1. XML code for the patchtest The following code illustrates the XML format for the Patchtest mesh:

References [1] B. Emek Abali, et al., Technical university of berlin, institute of mechanics, chair of continuums mechanics and material theory, computational reality, 2011. URL: /http:/ www.lkm.tu-berlin.de/ComputationalReality/S. [2] C. Abert, L. Exl, F. Bruckner, D. Suess, A. Drews, magnum.fe: a micromagnetic finite-element simulation code based on fenics, J. Magn. Magn. Mater. 345 (2013) 29–35. [3] E.C. Aifantis, Update on a class of gradient theories, Mech. Mater. 35 (2003) 259–280. [4] E.C. Aifantis, Non-singular dislocation fields, Mater. Sci. Eng. 3 (2009) 012026. [5] E.C. Aifantis, On scale invariance in anisotropic plasticity, gradient plasticity and gradient elasticity, Int. J. Eng. Sci. 47 (2009) 1089–1099. [6] E.C. Aifantis, A note on gradient elasticity and nonsingular crack fields, Int. J. Eng. Sci. 20 (2011) 103–105. [7] E. Amanatidou, A. Aravas, Mixed finite element formulations if strain-gradient elasticity problems, Comput. Methods Appl. Mech. Eng. 191 (2002) 1723–1751. [8] D.N. Arnold, F. Brezzi, J. Douglas, Peers: a new mixed finite element for plane elasticity, Jpn. J. Appl. Math. 1 (1984) 347–367. [9] D.N. Arnold, F. Brezzi, J. Douglas, Mixed finite element methods for elliptic problems, Comput. Methods Appl. Mech. Eng. 82 (1990) 281–300. [10] D.N. Arnold, F. Brezzi, M. Fortin, A stable nite element for the stokes equations, Calcolo 21 (4) (1984) 337–344. [11] D.N. Arnold, R.S. Falk, A new mixed formulation for elasticity, Jpn. J. Appl. Math. 1 (1984) 347–367. [12] D.N. Arnold, R.S. Falk, R. Winther, Mixed finite element methods for linear elasticity with weakly imposed symmetry, Math. Comput. 76 (2007) 1699–1723. [13] D.N. Arnold, R. Winther, Mixed finite elements for elasticity, Numer. Math. 92 (2002) 401–419. [14] D.N. Arnold, R. Winther, Mixed Finite Elements for Elasticity in the Stress–Displacement Formulation, American Mathematical Society, 2002. [15] H. Askes, E.C. Aifantis, Gradient elasticity and flexural wave dispersion in carbon nanotubes, Phys. Rev. B 80 (2009) 195412. [16] H. Askes, E.C. Aifantis, Gradient elasticity in statics and dynamics: an overview of formulations, length scale identification procedures, finite element implementations and new results, Int. J. Solids Struct. 48 (13) (2011) 1962–1990. [17] H. Askes, E.C. Aifantis, Numerical modeling of size effects with gradient elasticity – formulation, meshless discretization and examples, Int. J. Fract. 117 (2011) 347–358. [18] H. Askes, I. Morata, E.C. Aifantis, Finite element analysis with staggered gradient elasticity, Comput. Struct. 86 (2008) 1266–1279. [19] I. Babuka, M. Suri, Locking effects in the finite element approximation of elasticity problems, Numer. Math. 62 (1992) 439–463. [20] I. Babuska, The nite element method with lagrange multipliers, Numer. Math. 20 (1973) 179–192.

40

V. Phunpeng, P.M. Baiz / Finite Elements in Analysis and Design 96 (2015) 23–40

[21] F. Brezzi, On the existence, uniqueness, and approximations of saddle-point problems arising from lagrange multipliers, Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge 8 (R-2) (1974) 129–151. [22] E. Dick, Introduction to finite element methods in computational fluid dynamics, Comput. Fluid Dyn. (2009) 235–274. [23] A. Cemal Eringen, Nonlocal polar elastic continua, Int. J. Eng. Sci. 10 (1972) 1–16. [24] A. Cemal Eringen, On differential equations of non local elasticity and solutions of screw dislocation and surface waves, J. Appl. Phys. 54 (9) (1983) 1266–1279. [25] A. Cemal Eringen, Nonlocal Continuum Field Theories, Springer-Verlag, New York, 2002. [26] A. Cemal Eringen, D.G.B. Edelen, On nonlocal elasticity, Int. J. Eng. Sci. 10 (1972) 233–248. [27] A.J.M. Ferreira, MATLAB Codes for Finite Element Analysis: Solids and Structures, Springer, Dordrecht, London, 2008. [28] J. Hoffman, J. Jansson, R. Vilela de Abreu, N.C. Degirmenci, N. Jansson, M. Nazarov, K. Müller, J.H. Spühler, Unicorn Parallel adaptive finite element simulation of turbulent flow and fluid structure interaction for deforming domains and complex geometry, Comput. Fluids 80 (2013) 310–319. [29] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, 2000. [30] J.M. Jin, The Finite Element Method in Electromagnetics, 3rd edition, Wiley-IEEE Press, 2014. [31] J. Kioseoglou, I. Konstantopoulos, G. Ribarik, G.P. Dimitrakopulos, E.C. Aifantis, Non-singular dislocation and crack fields: implications to small volumes, Microsyst. Technol. 15 (2009) 117–121. [32] R.C. Kirby, Algorithm 839: fiat, a new paradigm for computing nite element basis functions, ACM Trans. Math. Softw. 30 (2004) 502–516. [33] R.C. Kirby, A. Logg, Efficient compilation of a class of variational forms, ACM Trans. Math. Softw. 33, 117 (2007), pp. 17, 20. [34] S. Kong, S. Zhou, Z. Nie, K. Wang, Static and dynamic analysis of micro beams based on strain gradient elasticity theory, Int. J. Eng. Sci. 47 (2009) 487–498. [35] D.C. Lam, F. Yang, A.C. Chong, J. Wang, P. Tong, Experiments and theory in strain gradient elasticity, J. Mech. Phys. Solids 51 (2003) 1477–1508. [36] A. Logg, K.A. Mardal, G.N. Wells, Automated Solution of Differential Equations by the Finite Element Method, Lecture Notes in Computational Science and Engineering, Springer, 2012. [37] J. Lu, D.V. Thiel, Computational and visual electromagnetics using an integrated programming language for undergraduate engineering students, IEEE Trans. Magn. 36 (2000) 1000–1003. [38] R.D. Mindlin, Micro-structure in linear elasticity, Arch. Ration. Mech. Anal. 10 (1964) 51–78. [39] R.D. Mindlin, Second gradient of strain and surface tension in linear elasticity, Int. J. Solids Struct. 1 (1965) 417–438. [40] R.D. Mindlin, N.N. Eshel, On first strain-gradient theories in linear elasticity, Int. J. Solids Struct. 4 (1968) 109–124. [41] R.D. Mindlin, Influence of couple-stresses on stress concentrations, Exp. Mech. 3 (1963) 1–7. [42] M. Mortensen, H.P. Langtangen, G.N. Wells, A fenics-based programming framework for modelling turbulent flow by the Reynolds-averaged Navierstokes equations, J. Comput. Phys. 263 (2014) 283–302. [43] K.B. Ølgaard, G.N. Wells, Optimizations for quadrature representations of nite element tensors through automated code generation, ACM Trans. Math. Softw. 37 (2010) 8:1–8:23. [44] C. Polizzotto, Stress gradient versus strain gradient constitutive models within elasticity, Int. J. Solids Struct. 45 (2008) 4820–4834. [45] C. Polizzotto, Thermodynamics-based gradientplasticity theories with an application to interface models, Int. J. Solids Struct. 45 (2008) 4820–4834. [46] C. Polizzotto, A gradient elasticity theory for second-grade materials and higher order inertia, Int. J. Solids Struct. 49 (2012) 2121–2137. [47] C. Polizzotto, A second strain gradient elasticity theory with second velocity gradient inertia. Part I. Constitutive equations and quasi-static behavior, Int. J. Solids Struct. 50 (2013) 749–765. [48] E.J. Rayfield, Finite element analysis and understanding the biomechanics and evolution of living and fossil organisms, Ann. Rev. Earth Planet. Sci. 35 (2007) 541–576. [49] J.N. Reddy, An Introduction to the Finite Element Method, McGraw-Hill, New York, 1993. [50] J.N. Reddy, S.D. Pang, Nonlocal continuum theories of beams for the analysis of carbon nanotubes, J. Appl. Phys. 103 (2008). [51] S. Schuttea, S.P.W. van den Bedemb, F. Keulenb, et al., A finite-element analysis model of orbital biomechanics, Vis. Res. 46 (2006) 1724–1731. [52] S. Shen, Finite element methods in fluid mechanics, Ann. Rev. Fluid Mech. 9 (1977) 421–445. [53] J.Y. Shu, W.E. King, N.A. Fleck, Finite elements for materials with strain gradient effects, Int. J. Numer. Methods Eng. 44 (1999) 373–391. [54] A.K. Soh, C. Wanji, Finite element formulations of strain gradient theory for microstructures and c0  1 patch test, Int. J. Numer. Methods Eng. 61 (2004) 433–454. [55] Dassault Systemes, ABAQUS 6.9 User Documentation, Internet Manual, Simulia., 2001. [56] Z. Tang, S. Shen, S.N. Atluri, Analysis of materials with strain-gradient effects: a meshless local Petrov–Galerkin (mlpg) approach, with nodal displacements only, Comput. Model. Eng. Sci. 4 (2003) 177–196. [57] The FEniCS Team, The Fenics Project, 2011. [58] R.A. Toupin, Elastic materials with couple stresses, Arch. Ration. Mech. Anal. 11 (1962) 385–414. [59] G.I. Tsamasphyros, C.D. Vrettos, A mixed finite volume formulation for the solution of gradient elasticity problems, Arch. Appl. Mech. 80 (2010) 609–627. [60] L. Vynnytska, M.E. Rognes, S.R. Clark, Benchmarking fenics for mantle convection simulations, Comput. Geosci. 50 (2013) 95–105. [61] J. Yang, Special Topics in the Theory of Piezoelectricity, Springer, USA, 2009. [62] A. Zervos, S.A. Papanicolopulos, I. Vardoulakis, Two finite-element discretizations for gradient elasticity, J. Eng. Mech. 135 (2009) 203–213. [63] J. Zhao, W.J. Chen, S.H. Lo, A refined nonconforming quadrilateral element for couple stress/strain gradient elasticity, Int. J. Numer. Methods Eng. 85 (2011) 269–288. [64] J. Zhao, C. Wanji, S.H. Lo, A refined nonconforming quadrilateral element for couple stress/strain gradient elasticity, Int. J. Numer. Methods Eng. 85 (2011) 269–288. [65] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method: The Basis, 5th edition, vol. 1, Butterworth-Heinemann, Oxford, 2000. [66] L. Zybell, U. Mühlich, M. Kuna, Z.L. Zhang, A three-dimensional finite element for gradient elasticity based on a mixed-type formulation, Comput. Mater. Sci. 52 (2012) 268–273.