Interior penalty discontinuous Galerkin FEMs for a gradient beam and CNTs

Applied Numerical Mathematics 144 (2019) 118–139

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Interior penalty discontinuous Galerkin FEMs for a gradient beam and CNTs K.G. Eptaimeros ∗ , C.Chr. Koutsoumaris, G.J. Tsamasphyros Division of Mechanics, School of Applied Mathematical and Physical Sciences, National Technical University of Athens, Greece

a r t i c l e

i n f o

Article history: Received 4 September 2018 Received in revised form 18 May 2019 Accepted 20 May 2019 Available online 23 May 2019 Keywords: Gradient elasticity Discontinuous Galerkin Interior penalty Euler-Bernoulli beam Carbon nanotubes Finite element

a b s t r a c t Carbon nanotubes (CNTs) are widely employed to modern engineering applications. The structural beam element is an acceptable model that has been used to simulate the CNTs’ response. At the same time the gradient elasticity theory, taking into account the size effect phenomena, can be considered a viable option to study the mechanical response of CNTs through a generalized gradient beam model. Given that the constitutive equations of this theory produce a boundary layer, the hp-version interior penalty discontinuous Galerkin finite element methods (IPDGFEMs) are developed in this work for the solution of a static, tensile or compressive, gradient elastic beam in macro- and micro-structure (CNT), respectively. An a priori error analysis of the aforementioned methods is performed and numerical simulations that verify the analytical findings are carried out. The consistency, the stability and the convergence of the IPDGFEMs are proved. A priori error estimates established also indicate an optimal convergence in the meshsize h, yet a slightly suboptimal convergence in the polynomial degree p. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The models of engineering, mathematics and physics are generally formulated by ordinary differential equations (ODEs) or partial differential equations (PDEs), supplemented by specific boundary and/or initial conditions. There are known procedures of finding the exact solution [25] and the approximate analytical solution [32,33,47,60] of specific differential equations. What is more, significantly interesting ODEs and PDEs have been numerically solved by employing the finite element methods (FEMs). The discontinuous Galerkin finite element method (DGFEM) is considered a substantial class of nonconforming methods. The first DGFEM was introduced in 1973 for the numerical solution of first-order hyperbolic problems [55]. Meanwhile, penalty methods were proposed for the numerical solution of second-order elliptic problems. The family of interior penalty discontinuous Galerkin (IPDG) methods combining the penalty methods with the DG methods was developed in the late 1970s and early 1980s [2,12,73]. At the same time the first contemporary DG method was presented for elliptic problems [5]. In fact, a number of methods belonging to this family has been suggested since then until nowadays [2–4,7,9,10,12,14,18, 19,21,22,27,35,39,48,49,56–58,62,70,73]. The effectiveness of DGFEMs for the solution of equations is justified thanks to their outstanding properties regarding the local conservation of the state variable, the good performance close to the solution’s discontinuities, the considerable

*

Corresponding author. E-mail addresses: [email protected] (K.G. Eptaimeros), [email protected] (C.Chr. Koutsoumaris).

https://doi.org/10.1016/j.apnum.2019.05.020 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

119

flexibility in the enforcement of the boundary and initial conditions and in the mesh design as well as the considerable variety of the polynomial degree in the local basis functions without introducing any interelement continuity requirements [2, 3,7,9,10,18,19,21,22,27,35,39,48,49,56–58,62,70–72]. The DGFEMs’ design is based on the variational formulation of the governing ODEs/PDEs and the properties of the suitable approximate space. The appropriate selection of the formulation and finite-dimensional space’s properties are capable of guaranteeing consistency, stability and convergence of the method. Meanwhile, the current use of nanomaterials (such as grapheme sheets, carbon nanotubes (CNTs) etc.) has led to a number of intricate problems in nanotechnology and engineering, since their macroscopic behavior often hinges on the material’s size effects. In particular, CNTs are widely used in recent and state of the art engineering applications, such as the construction of nanocomposites and biomaterials [24,36]. Regarding the investigation of CNTs’ mechanical response, the classic structural beam element is considered an acceptable tool for simulating CNTs [11,26,31,54,74]. An alternative route concerns the generalized nonlocal continuum theories, i.e. either higher-order gradient or integral, so that the scale effect can be taken into consideration [13,50,67,68]. Previous research projects have focused on the nonlocal integral form of (nano)beams [16,17,28,30,69]. Several experiments have been also performed on CNTs [63,75]. It is ergo worth investigating the CNTs’ response when their length dimension varies within limits of micrometer (10−6 m), where the gradient elasticity theory can be considered a viable option [6,40]. During the 1960s and 1970s, Mindlin [43–46] and Toupin [64] paid a tribute to the development of gradient theories by their innovative studies expanding on the microstructure’s dependence and its effect on the elastic response. Analytical solutions for problems of elasticity have been formulated and established by many researchers since then [1,20,23,34,37,38, 51–53,66]. The gradient elasticity theory is considered an expansion of the classical elasticity theory. Its elementary principle centers on the strain energy density being a positive-definite functional of the standard strain and the strain gradient of any order and subsequently the underlying microstructural effects of solids are taken into account. The microstructural magnitude of the constitutive equations is taken into consideration through the elastic energy deformation function depending on the strain tensor as well as the strain gradient and the rotation tensor, respectively. Toupin-Mindlin’s theory expands the definition of the continuum by bridging the gap between the classical theory of a continuum body and the theory of a crystallite lattice. Substantial research projects have been focused on the area of structural mechanics in gradient elasticity [1,23,34,37,38, 41,51,52,66]. Furthermore, a great number of researches has been carried out by applying FEMs for the numerical solution of the problems of the gradient elasticity theory [15,29,40,42,65]. What is more, this work revolves around the development of the hp-version IPDGFEMs, i.e. the symmetric interior penalty Galerkin (SIPG) and the nonsymmetric interior penalty Galerkin (NIPG), for the solution of a static, gradient elastic beam either in tension or compression with classic and microstructural material constants, respectively. In particular, our research objective focuses on performing a priori error analysis of the aforementioned methods as well as carrying out numerical simulations that verify the analytical findings. The consistency, the stability and the convergence of the IPDGFEMs are proven by the analytical results. A priori error estimates established also indicate an optimal convergence in the meshsize h, but a slightly suboptimal convergence in the polynomial degree p. What is more, the IPDG approximate displacements converge to the exact displacement for various discretizations and interpolation orders as well as for each length scale, g, that studied. Based on the conclusions reached, the IPDGFEMs developed are considered an effective tool for dealing with the problems of the gradient elasticity theory. Overall, this paper is structured as follows: Section 2 contains all the preliminary notions of the IPDG methods. The content of Section 3 deals with the general governing equation, representing a gradient elastic beam either in tension or compression. In Section 4, the boundary value problem (BVP) is deduced by just supplementing the governing equation with the suitable boundary conditions (BCs) of the gradient elasticity theory. Furthermore, Section 5 contains the IPDG weak formulation of the BVP and the consistency’s proof. Moreover, the content of Section 6 centers on the hp-version IPDGFEMs’ design for the BVP. The stability’s proof of these methods in the energy seminorm is presented as well. In Section 7, our research endeavor revolves around the proof of a priori error estimates in the energy seminorm. Section 8 then deals with the numerical simulations and the deducing results. Finally, Section 9 reiterates the conclusions of this research objective. 2. Preliminaries We denote by L p () the standard Lebesgue spaces, defined on an open subset  ⊂ , with norm || · || L p () (1  p  ∞) as well as by H k () the standard Hilbert-Sobolev space of index k  0 of real-valued functions with norm and seminorm || · ||k, and | · |k, , respectively. The norm of L 2 () is denoted by || · || for brevity. Let  be an open, bounded, convex domain in  with boundary  and let us consider a family of subdivisions P () of j . Specifically, P () is a partition of  into disjoint, open, convex element domains e = e such that

¯ = 

 e ∈P ()

¯ e, 

j

ei ∩ e = ∅ ∀ i = j ,

120

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

¯ e is empty or a vertex. The length of each element e is denoted by he . Let us assume that ¯ ei ∩  where the intersection  the family of subdivisions P () is regular [8]. The broken Sobolev space of composite order s on an open set , subjected to a subdivision P () of , is defined as j





H s (, P ()) = u ∈ L 2 () : u |e ∈ Hse (e ) ∀e ∈ P () , equipped with the associated broken norm || · ||s,P () . The local Sobolev space index of the element e is denoted by se and s := (se : e ∈ P ()) [21,62]. What is more, the finite-dimensional (sub)space (of the broken Sobolev space H s (, P ()) ) is defined as

  V hp = u ∈ L 2 ()| u ∈ P k (e ) ∀e ∈ P () ,

(1)

where k is a positive integer. The finite-dimensional space of all polynomials of degree less than or equal to k, defined on e , is denoted by P k (e ). Then, a nonnegative integer p e (the local polynomial index) is assigned to each e ∈ P (). The functions of V hp are discontinuous across the interior boundaries of the mesh. The union of element interiors and interior boundaries can be defined as

˜ = 

N el 

˜ =

e ,

e =1

Ni 

(2)

i ,

i =1

where N el is the number of finite elements and N i is the number of interior boundaries i . The following notation is adopted for the L 2 -inner product on element interiors and on interior boundaries

(a, b)˜ =

N el  (a, b)e ,

(a, b)˜ =

e =1

Ni  (a, b)i .

(3)

i =1

The norms on element interiors and interior boundaries are defined as

|| · ||2˜ =

N el  e =1

|| · ||2e ,

|| · ||2˜ =

Ni 

|| · ||2i .

(4)

i =1

˜ 0 = ˜ ∪ c and ˜ 1 = ˜ ∪ q , where c is the boundary of the axial displacement and q is the boundary of the Let  gradient of the axial displacement. We define the inner products

(a, b)˜ 0 = (a, b)˜ + (a, b)c ,

(t , v )˜ 1 = (t , v )˜ + (t , v )q ,

˜ 0 ) and t , v ∈ L 2 (˜ 1 ). The above inner products can be written in a compressed form as for a, b ∈ L 2 ( ab˜ = ab˜ + ab|c ,

t v ˜ = t v ˜ + t v |q ,

(5)

1

0

with associated norms|| · ||˜ and || · ||˜ . So, it will hold 1

0

||b||2˜ = ||b||2˜ + ||b||2c 0

or

||b||2˜ = 0

Ni 

||b||2i +

i =1

Nc 

||b||2r ,

(6)

r =1

and

|| v ||2˜ = || v ||2˜ + || v ||2q 1

or

|| v ||2˜ = 1

Ni 

|| v ||2i +

i =1

Nq 

|| v ||2 j ,

(7)

j =1

where N c is the number of exterior displacement boundary segments r ⊆ c and N q is the number of exterior displacement gradient boundary segments  j ⊆ q , as well. 3. Tension or compression of a gradient elastic beam The constitutive equations for the stress (or Cauchy stress) τ , the higher-order stress (or double stress) stress (or the axial stress on the yz-section in the direction of x) σ can be expressed as

τ = E  , μ = g 2 E ,x , σ = τ − μ,x = E ( − g 2 ,xx ),

μ and the total (8)

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

121

where E is the modulus of elasticity, g is a length scale representing a material length associated with the volumetric elastic strain energy,  is the axial strain of the beam either in tension or compression and the index ,x indicates differentiation with respect to x. The nonzero strain of Euler-Bernoulli beam theory is

 = u ,x ,

(9)

where u is the axial displacement. The equilibrium conditions require that the resultant of the internal forces should be equal to the axial force N on the uniform cross-section A. Thereby



σ d A = N.

(10)

A

Inserting the total stress

σ into eq. (10), it yields

A E ( − g 2 ,xx ) = N ,

(11)

where A E is the axial stiffness. The equilibrium equation of Euler-Bernoulli beam theory derives from

N ,x + ¯f = 0,

(12)

where ¯f is an axial force per unit length (an axially distributed load). We ergo deduce the following equation from eqs. (11) – (12)

( g 2 ,xx −  ),x =

¯f AE

= f.

(13)

4. Model problem Let us introduce an one-dimensional model problem following the concepts of Toupin and Mindlin [43–46,64]. Let  ⊂  be an open, bounded domain (i.e. a beam of length L),  be its boundary as well as c , q ,  R and  P be the axial displacement, the displacement gradient, the double force and the axial force boundaries, respectively. Eq. (13) is defined in  as

( g 2 ,xx −  ),x =

¯f AE

=f

in ,

(14)

where f ∈ L 2 (). Eq. (14) is then supplemented with the following BCs

u=c

on c ,

 · n = q on q ,

(15)

μ A = R on  R , σ A · n = P on  P , where n is the unit normal vector to the boundary, exterior to . The following relationships hold between the boundary’s parts:

c ∪  P = ,

c ∩  P = ∅,

q ∪  R = ,

q ∩  R = ∅.

(16)

We can rewrite eq. (14) and eq. (15) with eqs. (8) – (9) as:

( g 2 u ,xx − u ),xx =

¯f

=f

in ,

AE u=c

on c ,

u ,x · n = q

on q ,

2

A E g u ,xx = R

on  R ,

A E (u − g 2 u ,xx ),x · n = P

on  P ,

(17)

(18)

where c, q, R and P are the prescribed boundary of the axial displacement, the displacement gradient, the double force and the axial force, respectively. The BVP (17) – (18) has a unique solution u ∈ H 4 () depending on the problem’s data [41,42,65].

122

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

5. Weak formulation In this section, the weak formulation of the BVP (17) – (18) leading to the IPDGFEMs derives. We suppose that the problem’s solution u is a sufficiently smooth function. ˜ , the elements e := ke and e := le share the interior boundary i , where k and l For each interior boundary i ⊆  are indices with k > l. The jump across i and the mean value on i of u are defined by

[[u ]]i := u |∂ e ∩i − u |∂ e ∩i

and

u i :=

1 2



u |∂ e ∩i + u |∂ e ∩i .

The jump and the mean value definitions are extended to r ⊆ c and  j ⊆ q belonging to the boundary 

[[u ]]r = u |r ,

u r = u |r

and

[[u ]] j = u | j ,

u  j = u | j .

˜ The unit normal vector n = nk , pointing from element ke to le when k > l, is related to each interior boundary i ⊆  e and n = ne is selected to be the outward unit normal vector when a node belongs to the boundary . Since the method is nonconforming (i.e. the finite element space is not a finite-dimensional subspace of the Sobolev space H 2 ()), the broken Sobolev space H 4 (, P ()) will be used as a trial space. We multiply eq. (17) by a test function w ∈ H 4 (, P ()) and next integrate over 

(( g 2 u ,xx − u ),xx , w ) = ( f , w ) .

(19)

Afterwards, we split the inner products and then integrate by parts on every elemental integral, so we get N el N el   ( g 2 u ,xx , w ,xx )e + (( g 2 u ,xx − u ),x · n, w )∂ e e =1

e =1

N el N el    − ( g 2 u ,xx , w ,x · n)∂ e + (u ,x , w ,x )e = ( f , w )e , N el

e =1

e =1

(20)

e =1

where the outward normal to each element boundary is denoted by n. Splitting the inner products over the boundary on the left-hand side of eq. (20) and using the notations (3) and (5), we obtain

( g 2 u ,xx , w ,xx )˜ + (u ,x , w ,x )˜ + ( g 2 u ,xx − u ),x · nw ˜ 0 − g 2 u ,xx w ,x · n˜ 1 = ( f , w )˜ +

P AE

w | P +

R AE

w , x · n | R .

(21)

˜ 0 and ˜ 1 (˜ 0 = ˜ ∪ c , ˜ 1 = ˜ ∪ q ) are evaluated by applying the mathematical The inner products on the boundaries  ˜ (for identity JabK = aJbK + JaKb. The fluxes ( g 2 u ,xx − u ),x · n and g 2 u ,xx are continuous across the interior boundaries  4 2 2 example when the exact solution u ∈ H ()). The terms, including the jumps J( g u ,xx − u ),x K and J g u ,xx K, are therefore neglected without modifying the consistency. We then deduce

( g 2 u ,xx , w ,xx )˜ + (u ,x , w ,x )˜ + ( g 2 u ,xx ),x [[ w ]]˜ 0 −  g 2 u ,xx [[ w ,x ]]˜ 1 − u ,x [[ w ]]˜ 0 = ( f , w )˜ +

P AE

w | P +

R AE

w , x · n | R .

(22)

The following extra terms are introduced to guarantee the compatibility and symmetry/nonsymmetry of the formulation. Since the exact solution u and its gradient u ,x are continuous, the jumps Ju K and Ju ,x K vanish (i.e. Ju K = Ju ,x K = 0). The terms, including the aforementioned jumps, do not alter the problem’s consistency. If we select −θ( g 2 w ,xx ),x  + α J w K, θ w ,x  + δJ w K, θ g 2 w ,xx  + βJ w ,x K as test functions and then integrate the first two over ˜ 0 and the last one over ˜ 1 , we have

−θ[[u ]]( g 2 w ,xx ),x ˜ 0 + α [[u ]][[ w ]]˜ 0 = −θ c ( g 2 w ,xx ),x · n|c + αc c w |c ,

(23)

θ[[u ]] w ,x ˜ 0 + δ[[u ]][[ w ]]˜ 0 = θ c w ,x · n|c + δc c w |c ,

(24)

θ[[u ,x ]] g 2 w ,xx ˜ 1 + β[[u ,x ]][[ w ,x ]]˜ 1 = θ qg 2 w ,xx |q + βq qw ,x · n|q ,

(25)

where θ is the symmetrization parameter and θ ∈ {−1, 1}. In addition, the nonnegative piecewise continuous functions, ˜ 0 and ˜ 1 and referred to as stabilization parameters, are denoted by α , δ and β . defined on  Adding eqs. (22) – (25), it yields the problem’s IPDG weak formulation. The left-hand side of the deducing equation represents the bilinear form B sb (·, ·) and the right-hand side represents the linear functional L sb (·), respectively.

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

123

The bilinear form B sb (·, ·) is defined as

B sb (u , w ) := ( g 2 u ,xx , w ,xx )˜ + (u ,x , w ,x )˜ + ( g 2 u ,xx ),x [[ w ]]˜

0

− θ[[u ]]( g 2 w ,xx ),x ˜ 0 −  g 2 u ,xx [[ w ,x ]]˜ 1 + θ[[u ,x ]] g 2 w ,xx ˜ 1 − u ,x [[ w ]]˜ 0 + θ[[u ]] w ,x ˜ 0 + α [[u ]][[ w ]]˜ 0 + β[[u ,x ]][[ w ,x ]]˜ 1 + δ[[u ]][[ w ]]˜ 0 .

(26)

The bilinear form is symmetric for θ = −1, but it is nonsymmetric for θ = 1. The linear functional L sb (·) on H 4 (, P ()) is defined as

L sb ( w ) := ( f , w )˜ +

P AE

w | P +

R AE

w ,x · n| R − θ c ( g 2 w ,xx ),x · n|c

(27)

+ θ qg 2 w ,xx |q + θ c w ,x · n|c + αc c w |c + βq qw ,x · n|q + δc c w |c . The stabilization parameters α , αc , β, βq , δ, δc that will be determined below are dependent on the discretization parameters he and p e . The broken weak formulation (BWF) of the BVP (17) – (18) then reads:

∀ w ∈ H 4 (, P ()).

Find u ∈ bSs such that B sb (u , w ) = L sb ( w )

(28)

We denote by bSs the following function space

bSs = {u ∈ H 4 (, P ()) : u , u ,x · n, g 2 u ,xx , ( g 2 u ,xx − u ),x · n are continuous across i }. The energy seminorm, ||| · |||sb , is related to the bilinear form B sb (·, ·) and defined as

|||u |||sb = ||( g 2 )1/2 u ,xx ||2˜ + ||u ,x ||2˜ + ||α 1/2 [[u ]]||2˜ 0

1/2 +||β 1/2 [[u ,x ]]||2˜ + ||δ 1/2 [[u ]]||2˜ , u ∈ H 2 (, P ()). 1

(29)

0

Proposition 5.0.1. If α , αc , β , βq , δ , δc > 0, then ||| · |||sb is a seminorm on H 2 (, P ()). Proof. Let the quantity |||u |||sb vanishes for the piecewise constant functions, being the only non-trivial functions in H 2 (, P ()). There are ergo two cases: Either u is constant and nonzero over the entire domain  or u changes the axial displacement’s value at least once. 1/ 2 1/ 2 The former case corresponds to u = 0 and hence ||αc u ||c > 0 as well as ||δc u ||c > 0. On the other hand, the 1/ 2 latter case corresponds to Ju K = 0 at least one interior boundary, so that ||α Ju K||˜ > 0 and ||δ 1/2 Ju K||˜ > 0. As a result, |||u |||sb > 0 ∀u ∈ H 2 (, P ()) and the quantity |||u |||sb can not vanish for the piecewise constant functions. The other two axioms of the norm being the homogeneity and the triangle inequality are verified without difficulty. 2 We notice that ||| · |||sb is also a seminorm on H 4 (, P ()), since H 4 (, P ()) ⊂ H 2 (, P ()). The bilinear form (26) of the IPDG methods is well-defined on the finite element space V hp , yet it is not well-defined on H 2 (). This can be remedied by an alternative extension (i.e. B˜ sb (·, ·)) of the bilinear form to V hp + H 2 () [21]. However, this formulation is inconsistent. Furthermore, the fourth-order equation presented is supplemented with Dirichlet BCs, which are simpler in comparison with the BCs of the BVP (17) – (18). As the BCs of the BVP (17) – (18) are complicated, it is difficult to follow the procedure suggested in [21], define the lifting operators, the boundaries’ lifting and then conduct an error analysis. Given the aforementioned reasons and justifications, the procedure and the formulation introduced in [62] are followed in this study, although the bilinear form of the IPDG methods is indeed well-defined on the finite element space, but it is not well-defined on H 2 (). 5.1. Consistency It will be demonstrated that a strong solution to the BVP (17) – (18), smooth enough at the interior boundaries, is the solution to the problem in the BWF (28). To begin with, the weak continuity of fluxes across the interior boundaries i is shown.

124

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

Lemma 5.1.1. Let suppose that u ∈ H 4 (); we then get for any interior boundary i

[[u ]] w i = [[u ,x ]] w i = [[ g 2 u ,xx ]] w i = [[( g 2 u ,xx − u ),x ]] w i = 0 ∀ w ∈ L 2 (i ). Proof. According to Romkes et al. [59], the first two quantities are equal to zero for all w in L 2 (i ) when u ∈ H 2 (). The last equality is established as follows: Let i be an interior boundary and let e and e be the elements which ˜ e = int(e ∪ e ). After integrating by parts, we subsequently get for any w ∈ share the interior boundary i . Let also  ˜ e ) = C ∞ ( ˜ e) D( 0

(( g 2 u ,xx − u ),xx , w )˜ e = (( g 2 u ,xx − u ),x · n, w )∂ ˜ e − (( g 2 u ,xx ),x , w ,x )˜ e + (u ,x , w ,x )˜ e = −(( g 2 u ,xx ),x , w ,x )˜ e + (u ,x , w ,x )˜ e .

(30)

Then, the left-hand side inner product is also split and the integration by parts formula is applied to e and e . As a result, we deduce

(( g 2 u ,xx − u ),xx , w )˜ e = (( g 2 u ,xx − u ),xx , w )e + (( g 2 u ,xx − u ),xx , w )e = −(( g 2 u ,xx ),x , w ,x )˜ e + (u ,x , w ,x )˜ e + [[( g 2 u ,xx − u ),x ]] · nw i .

(31)

Thus, the identities (30) and (31) entail that

˜ e ). [[( g 2 u ,xx − u ),x ]] · nw i = 0 ∀ w ∈ D(

(32)

Ergo,

[[( g 2 u ,xx − u ),x ]] · nw i = 0 ∀ w ∈ D(i ). Since D (i ) is dense in L 2 (i ), it implies that

[[( g 2 u ,xx − u ),x ]] · nw i = 0 ∀ w ∈ L 2 (i ). To establish the equality J g 2 u ,xx K w i = 0, similar series of steps are used as well. We therefore reach the conclusion

[[ g 2 u ,xx ]] w ,x · ni = 0 ∀ w ∈ L 2 (i ). 2 Proposition 5.1.2. The BWF (28) of the BVP (17) – (18) is consistent in the space H 4 (), i.e. any solution u to the BVP, such that u ∈ H 4 (), also solves (28). Proof. For u ∈ bSs, we obtain from (28) as well as the expressions of B sb (·, ·) and L sb (·)

0 = B sb (u , w ) − L sb ( w )

= ( g 2 u ,xx , w ,xx )˜ + (u ,x , w ,x )˜ + ( g 2 u ,xx ),x [[ w ]]˜ 0 − θ[[u ]]( g 2 w ,xx ),x ˜ 0 −  g 2 u ,xx [[ w ,x ]]˜ 1 + θ[[u ,x ]] g 2 w ,xx ˜ 1 − u ,x [[ w ]]˜ 0 + θ[[u ]] w ,x ˜ 0 + α [[u ]][[ w ]]˜ 0 + β[[u ,x ]][[ w ,x ]]˜ 1 + δ[[u ]][[ w ]]˜ 0 R w ,x · n| R + θ c ( g 2 w ,xx ),x · n|c AE − θ qg 2 w ,xx |q − θ c w ,x · n|c − αc c w |c − βq qw ,x · n|q − δc c w |c .

− ( f , w )˜ −

P

(33)

AE

w | P −

Next, integrating by parts (u ,x , w ,x )˜ and ( g 2 u ,xx , w ,xx )˜ , we deduce

0 = (( g 2 u ,xx − u ),xx − f , w )˜ + [[(u − g 2 u ,xx ),x ]] w ˜ + [[ g 2 u ,xx ]] w ,x ˜

+ θ[[u ,x ]] g 2 w ,xx ˜ + θ[[u ]]( w − g 2 w ,xx ),x ˜ + θ(u − c )( w − g 2 w ,xx ),x · n|c + θ(u ,x · n − q) g 2 w ,xx |q   R P + g 2 u ,xx − w ,x · n| R + (u − g 2 u ,xx ),x · n − w | P AE

AE

+ α [[u ]][[ w ]]˜ + β[[u ,x ]][[ w ,x ]]˜ + δ[[u ]][[ w ]]˜ + αc (u − c ) w |c + βq (u ,x · n − q) w ,x · n|q + δc (u − c ) w |c .

(34)

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

125

Eq. (34) is identical to zero for all w, when

[[u ]] = [[u ,x ]] = [[ A E g 2 u ,xx ]] = [[ A E (u − g 2 u ,xx ),x ]] = 0 on ˜ ,

(35)

A E ( g 2 u ,xx − u ),xx − ¯f = 0

˜, in 

(36)

u=c

on c ,

(37)

u ,x · n = q

on q ,

(38)

A E g u ,xx = R

on  R ,

(39)

A E (u − g u ,xx ),x · n = P

on  P .

(40)

2

2

The continuity of the (axial) displacement, the displacement gradient, the double force and the (axial) force across the interior boundaries is ensured by eq. (35) (see Lemma (5.1.1)). Also, the enforcement of the governing differential equation on the element interiors is denoted by eq. (36). Eqs. (37) – (40) account for the enforcement of the BCs as well. Wherefore, we conclude that any solution u ∈ H 4 () to the BVP (17) – (18) is a weak discontinuous solution of (28). 2 6. IPDG finite element method The IPDGFEM is presented in this section for the problem (17) – (18):

Find u D G ∈ V hp such that B sb (u D G , w ) = L sb ( w )

∀ w ∈ V hp .

(41)

The functions α , αc , β , βq , δ , δc of B sb (·, ·) and L sb (·) will be defined in the coercivity property. When θ = −1, the IPDGFEM corresponds to the SIPG formulation having a symmetric bilinear form. Provided that the stabilization parameters are selected sufficiently large, this bilinear form will be coercive. This formulation was introduced for second-order elliptic equations [2,73] and fourth-order elliptic equations [62]. In addition, when θ = 1, the IPDGFEM corresponds to the NIPG formulation which was introduced for second-order elliptic equations [27,56–58] and fourth-order elliptic equations [48]. The approximate solution u D G ∈ V hp will be generally discontinuous, since there does not exist a continuity requirement in the finite element space. What is more, let us suppose that the strong solution u to the BVP satisfies the assumption of the smoothness u ∈ H 4 (). This ensures that u is a solution to (28) and ergo to (41). Galerkin orthogonality property is therefore defined as

B sb (u − u D G , w ) = 0 ∀ w ∈ V hp ,

(42)

where u is the solution of the problem and u D G is the IPDG approximate solution defined by the method (41). 6.1. Coercivity of bilinear form We prove that the bilinear form B sb (·, ·) of the IPDG methods (either SIPG or NIPG) is coercive and therefore the problem (41) will have a unique solution. Proposition 6.1.1. The hp-version IPDG methods (41) are stable in the energy seminorm (29), i.e. there is a positive constant m such that

B sb (u , u )  m|||u |||2sb

∀u ∈ V hp .

(43)

Proof. a) hp-version SIPG method: Substituting u for w and taking the symmetric bilinear form (26) where θ = −1, we then obtain

B sb (u , u )  ||( g 2 )1/2 u ,xx ||2˜ + ||u ,x ||2˜ 







+ 2( g u ,xx ),x [[u ]]˜ 0 − 2  g 2 u ,xx [[u ,x ]]˜ 1 − 2 u ,x [[u ]]˜ 0 2

+ ||α 1/2 [[u ]]||2˜ + ||β 1/2 [[u ,x ]]||2˜ + ||δ 1/2 [[u ]]||2˜ . 0

1

0

(44)

˜ 0 and ˜ 1 on the right-hand side of (44) should be estimated. To complete the proof, the inner products on  The first inner product can be written by using Cauchy-Schwarz and Young inequalities as well as the inverse inequality [61] as

126

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

( g 2 u ,xx ),x [[u ]]˜ 0  ||( g 2 u ,xx ),x ||˜ 0 ||[[u ]]||˜ 0 

N el  ε1 e =1



N el  ε1 e =1

||( g 2 u ,xx ),x ||2∂ e +

2

2

c1 g 2

p e6 he3

1

||[[u ]]||2˜

2ε1

(45)

0

1

||( g 2 )1/2 u ,xx ||2e +

2ε1 α

||α 1/2 [[u ]]||2˜ , 0

where the constant c 1 is independent of he , p e and u. From (45), the first inner product on (44) can be estimated as follows

( g 2 u ,xx ),x [[u ]]˜ 0 

N el  ε1 e =1

2

c1 g 2

p e6 he3

1

||( g 2 )1/2 u ,xx ||2e +

2ε1 α

||α 1/2 [[u ]]||2˜ . 0

(46)

Similar arguments are used to estimate the rest of the inner products on the right-hand side of (44). We therefore arrive at the conclusion N el

 ε2 2 p e2

2

||( g 2 )1/2 u ,xx ||2e + c2 g

 g u ,xx [[u ,x ]]˜ 1  e =1

2

he

N el

 ε3 p e2

||u ,x ||2e + c3

u ,x [[u ]]˜ 0 

2

e =1

he

1 2ε3 δ

1

||β 1/2 [[u ,x ]]||2˜ ,

2ε2 β

1

||δ 1/2 [[u ]]||2˜ ,

(47)

(48)

0

where the constants c 2 and c 3 are independent of he , p e and u. After these series of steps, we gather the inequalities (46) – (48) and insert them into the right-hand side of (44). Ergo, we have

B sb (u , u ) 

N el  

1 − ε1 c 1 g

2

e =1

+

N el  

1 − ε3 c 3

e =1



+ 1−

1

ε2 β



p e6 he3 p e2 he

− ε2 c 2 g

2

p e2

||( g 2 )1/2 u ,xx ||2e

he

 1 ||u ,x ||2e + 1 − ||α 1/2 [[u ]]||2˜

ε1 α



||β 1/2 [[u ,x ]]||2˜ + 1 − 1

1

ε3 δ

0

||δ 1/2 [[u ]]||2˜

0

 m|||u |||2sb .

(49)

The minimum of the terms enclosed into the parentheses on the right-hand side of (49) is denoted by the constant m. To be more specific, (43) can be proved for m = 1/2 provided that

ε1 |e =

he3 4c 1 g 2 p e6

ε2 |e =

,

he 4c 2 g 2 p e2

and

ε3 |e =

he 2c 3 p e2

,

where we get

α = αc =

8c 1 g 2 p e6 he3

,

β = βq =

8c 2 g 2 p e2 he

and

δ = δc =

4c 3 p e2 he

.

b) hp-version NIPG method: Substituting u for w and taking the nonsymmetric bilinear form (26) where θ = 1, we then get

B sb (u , u ) = |||u |||2sb .

2

Wherefore, the bilinear form B sb (·, ·) is coercive on the finite-dimensional space V hp and the problem (41) ergo has a unique solution. 7. Error analysis An error analysis is performed in this section for IPDGFEMs (41). Our main concern is to prove hp-version a priori error estimates in the energy seminorm, ||| · |||sb , for the aforementioned methods. According to the consistency and the stability properties having proven above, the convergence of the methods can be shown.

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

127

7.1. Error estimates in the energy seminorm Convergence 7.1.1. Let p be a linear projection operator from H s (, P ()) to the finite element space V hp . The global error u − u D G can be then split as follows:

u − u D G = (u − p u ) + ( p u − u D G ) ≡ η + ξ.

(50)

Hence, we get by using the triangle inequality

|||u − u D G |||sb  |||η|||sb + |||ξ |||sb .

(51)

The part of the error in the finite element space is denoted by ξ = p u − u D G , where ξ ∈ V

hp

.

A bound on |||ξ |||sb in terms of norms of η will be provided by the following error analysis. An estimate on |||u − u D G |||sb with respect to different norms of η will be therefore obtained. The error analysis is ergo completed by quantifying the norms of η with regard to Sobolev seminorms of the solution u and the discretization parameters. For that reason, the Lemma 8 of [62] will be recalled. Theorem 7.1.2. Let us suppose that  is a bounded domain in  and P () is a regular partition of  into elements e . Let p = ( p e : e ∈ P (), p e ∈ ℵ, p e  3) be any polynomial degree vector of bounded local variation as well. The positive real piecewise constant functions α , αc , β , βq , δ and δc in each element are defined by

α = αc =

C α g 2 p e6 he3

, β = βq =

C β g 2 p e2 he

and

δ = δc =

C δ p e2 he

.

Let us also suppose that the bilinear form B sb (·, ·) is coercive (see Proposition 6.1.1). Provided that the solution u of the problem (28) belongs to the broken Sobolev space H t (, P ()), t = (t e : e ∈ P (), t e  4), the solution u D G ∈ V hp of the problem (41) will satisfy the next error bound

|||u − u D G |||2sb  C

N el 2s −4  he e 2t e −7

e =1 p e

||u ||t2e ,e ,

(52)

where 2  se  min( p e + 1, t e ) as well as C is a constant solely depends on the space dimension and t = maxe ∈P () t e . Proof. We will first estimate ξ . For that purpose, we employ the coercivity (43), the error decomposition (50) and Galerkin orthogonality (42) yielding

m|||ξ |||2sb  B sb (ξ, ξ )  | B sb (η, ξ )|.

(53)

We continue by writing the right-hand side of (53) as

m|||ξ |||2sb  |( g 2 η,xx , ξ,xx )˜ | + |(η,x , ξ,x )˜ |

+ |( g 2 η,xx ),x [[ξ ]]˜ 0 | + |[[η]]( g 2 ξ,xx ),x ˜ 0 | + | g 2 η,xx [[ξ,x ]]˜ 1 | + |[[η,x ]] g 2 ξ,xx ˜ 1 | + |η,x [[ξ ]]˜ 0 | + |[[η]]ξ,x ˜ 0 | + |α [[η]][[ξ ]]˜ 0 | + |β[[η,x ]][[ξ,x ]]˜ 1 | + |δ[[η]][[ξ ]]˜ 0 |.

(54)

To provide a bound on |||ξ |||sb with regard to appropriate norms of η , the inner products on the right-hand side of (54) should be estimated. Cauchy-Schwarz inequality and then the energy seminorm’s definition (29) are applied, so the first inner product on the right-hand side of (54) gives

|( g 2 η,xx , ξ,xx )˜ |  ||( g 2 )1/2 η,xx ||˜ ||( g 2 )1/2 ξ,xx ||˜  |||η|||sb |||ξ |||sb .

(55)

Based on (55), the first inner product of (54) can be ergo bounded as

|( g 2 η,xx , ξ,xx )˜ |  |||η|||sb |||ξ |||sb .

(56)

The second inner product of (54) can be correspondingly estimated as

|(η,x , ξ,x )˜ |  |||η|||sb |||ξ |||sb .

(57)

128

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

To estimate the stabilizing terms on the right-hand side of (54), similar series of steps are followed. Employing CauchySchwarz inequality and next the definition of the energy seminorm (29) to the first stabilizing term, we deduce

|α [[η]][[ξ ]]˜ 0 |  ||α 1/2 [[η]]||˜ 0 ||α 1/2 [[ξ ]]||˜ 0  |||η|||sb |||ξ |||sb .

(58)

From (58), the first stabilizing term of (54) can be bounded as follows

|α [[η]][[ξ ]]˜ 0 |  |||η|||sb |||ξ |||sb .

(59)

Also, the rest of stabilizing terms of (54) can be similarly bounded as

|β[[η,x ]][[ξ,x ]]˜ 1 |  |||η|||sb |||ξ |||sb ,

(60)

|δ[[η]][[ξ ]]˜ 0 |  |||η|||sb |||ξ |||sb .

We estimate the inner products of the mean value of η and the jump of ξ on the right-hand side of (54). Invoking Cauchy-Schwarz inequality and then the energy seminorm’s definition (29), it yields

|( g 2 η,xx ),x [[ξ ]]˜ 0 |  ||α −1/2 ( g 2 η,xx ),x ||˜ 0 ||α 1/2 [[ξ ]]||˜ 0  ||α −1/2 ( g 2 η,xx ),x ||˜ 0 |||ξ |||sb .

(61)

From (61), this inner product of (54) can be subsequently bounded as

|( g 2 η,xx ),x [[ξ ]]˜ 0 |  ||α −1/2 ( g 2 η,xx ),x ||˜ 0 |||ξ |||sb .

(62)

Furthermore, the remaining inner products of the corresponding form of (54) are estimated by using similar arguments. Thus, we deduce

| g 2 η,xx [[ξ,x ]]˜ 1 |  ||β −1/2  g 2 η,xx ||˜ 1 |||ξ |||sb ,

(63)

|η,x [[ξ ]]˜ 0 |  ||δ −1/2 η,x ||˜ 0 |||ξ |||sb .

To bound |||ξ |||sb in terms of norms of η , a last step is to estimate the rest of inner products on the right-hand side of (54). As in the latter case and since ξ ∈ V hp , we get by recalling the Cauchy-Schwarz inequality and afterwards the inverse inequality [61]

|[[η]]( g 2 ξ,xx ),x ˜ 0 |  ||α 1/2 [[η]]||˜ 0 ||α −1/2 ( g 2 ξ,xx ),x ||˜ 0 ⎛ ⎞1/2 N el   ||α 1/2 [[η]]||˜ 0 ⎝ ||α −1/2 ( g 2 ξ,xx ),x ||2∂ e ⎠ e =1

 ||α 1/2 [[η]]||˜ 0

⎛ ⎞1/2 N el  p e6 ⎝ c1 ||α −1/2 g 2 ξ,xx ||2 ⎠ e =1

he3

e

⎛ ⎞1/2 N el  c 1  ||α 1/2 [[η]]||˜ 0 ⎝ ||( g 2 )1/2 ξ,xx ||2e ⎠ , e =1

Cα

(64)

where the constant c 1 is independent of ξ , he and p e . Moreover, we select c 1 /C α  1 without loss of generality. Thereby, the use of energy seminorm’s definition (29) gives

⎛ ⎞1/2 N el  c 1 ||α 1/2 [[η]]||˜ 0 ⎝ ||( g 2 )1/2 ξ,xx ||2e ⎠  ||α 1/2 [[η]]||˜ 0 ||( g 2 )1/2 ξ,xx ||˜ e =1

Cα

 |||η|||sb |||ξ |||sb .

(65)

From (64) – (65), this inner product of (54) can be bounded as follows

|[[η]]( g 2 ξ,xx ),x ˜ 0 |  |||η|||sb |||ξ |||sb .

(66)

What is more, the rest of inner products of the corresponding form of (54) are estimated by following the aforementioned procedure in a similar manner. As a result, we have

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

129

|[[η,x ]] g 2 ξ,xx ˜ 1 |  |||η|||sb |||ξ |||sb ,

(67)

|[[η]]ξ,x ˜ 0 |  |||η|||sb |||ξ |||sb .

Thereafter we insert the inequalities (56) – (57), (59) – (60), (62) – (63), (66) – (67) into the right-hand side of (54). It then yields

  |||ξ |||sb  C |||η|||sb + ||α −1/2 (η,xx ),x ||˜ 0 + ||β −1/2 η,xx ||˜ 1 + ||δ −1/2 η,x ||˜ 0 .

(68)

By combining the mathematical expression (51) with (68), we deduce

  |||u − u D G |||2sb  C |||η|||2sb + ||α −1/2 (η,xx ),x ||2˜ + ||β −1/2 η,xx ||2˜ + ||δ −1/2 η,x ||2˜ . 1

0

0

(69)

A bound on |||u − u D G |||sb regarding norms of η has been therefore obtained. To complete the proof, the terms on the right-hand side of (69) should be estimated. We notice that η ∈ / V hp . Using the energy seminorm’s definition (29), the first term of (69) yields

|||η|||2sb  C

N el   e =1

 ||η,xx ||2e + ||η,x ||2e + ||α 1/2 [[η]]||2˜ + ||β 1/2 [[η,x ]]||2˜ + ||δ 1/2 [[η]]||2˜ . 1

0

0

(70)

By recalling Lemma 8 of [62] for the first two norms of (70), we obtain

||η,x ||e  C

s −1

hee

t −1 p ee

||u ||te ,e ,

||η,xx ||e  C

s −2

hee

t −2

p ee

||u ||te ,e .

Next, we estimate the term of the stabilization parameter

||α 1/2 [[η]]||2˜  C 0

N el  p e6

h3 e =1 e

(71)

α of (70). It is therefore deduced

||η||2∂ e .

(72)

Now, employing Lemma 8 of [62] to (72), we get

C

N el  p e6 e =1

he3

2

||η||∂ e  C

N el 2s −4  he e e =1

2t e −7

pe

||u ||t2e ,e .

(73)

From (72) – (73), the factor with the stabilization parameter

||α 1/2 [[η]]||2˜  C 0

N el 2s −4  he e 2t e −7

e =1 p e

α of (70) can be bounded as follows

||u ||t2e ,e .

(74)

The remaining terms with the stabilization parameters β and δ of (70) can be similarly bounded as follows

||β 1/2 [[η,x ]]||2˜  C 1

||δ

1/2

2

[[η]]||˜  C 0

N el 2s −4  he e 2t e −5 e =1 p e

||u ||t2e ,e , (75)

N el 2s −2  he e

||u ||t2e ,e . 2t e −3 p e =1 e

After that the insertion of the inequalities (71) and (74) – (75) into the right-hand side of (70) yields

|||η

|||2sb

C

 N el 2s −4  he e e =1

2t −4 pe e

+

2se −4

he

2t −5 pe e

+

2se −4

he



2t −7 pe e

 he2se −4 ||u ||t2e ,e . 2t e −7 p e =1 e

||u ||t2e ,e

+C

 N el 2s −2  he e e =1

2t −2 pe e

+

2se −2

he

2t −3 pe e

 ||u ||t2e ,e

N el

C

(76)

The remaining factors on the right-hand side of (69) will be estimated as well. The term with the stabilization parameter

α can be written as ||α −1/2 (η,xx ),x ||2˜  C 0

N el  h3 e

p6 e =1 e

||(η,xx ),x ||2∂ e .

(77)

130

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

Now, using Lemma 8 of [62] to (77), we have

C

N el  h3

e ||( ,xx ),x ||2∂ e 6 p e =1 e

η

C

N el 2s −4  he e 2t e −1

e =1 p e

||u ||t2e ,e .

(78)

From (77) – (78), the term with the stabilization parameter

||α −1/2 (η,xx ),x ||2˜  C 0

N el 2s −4  he e 2t e −1

e =1 p e

α of (69) can be bounded as follows

||u ||t2e ,e .

(79)

The rest of the terms with the stabilization parameters β and δ of (69) can be similarly estimated as

||β −1/2 η,xx ||2˜  C 1

||δ −1/2 η,x ||2˜  C 0

N el 2s −4  he e 2t e −3

e =1 p e

N el 2s −2  he e 2t e −1

e =1 p e

||u ||t2e ,e , (80)

||u ||t2e ,e .

Inserting the inequalities (76), (79) – (80) into the right-hand side of (69), it gives

|||u

− u D G |||2sb

C

N el 2s −2  he e

||u ||t2e ,e 2t e −1 p e e =1

+C

 N el 2s −4  he e e =1

 he2se −4 ||u ||t2e ,e . 2t e −7 p e e =1

2t e −1

pe

+

2se −4

he

2t e −3

pe

+

2se −4

he

2t e −7

pe

 ||u ||t2e ,e

N el

C

2

(81)

It is worth mentioning that the resulting a priori error estimate of the hp-version IPDG methods (either SIPG or NIPG) is optimal in h, yet it is p-suboptimal by 3/2 orders of p. Similar conclusions have been also drawn for the biharmonic equation [62]. 8. Numerical validation In this section, the hp-version IPDGFEMs are numerically tested for a gradient elastic beam in tension or compression with classic and microstructural material constants. Thus, a specific BVP appearing a boundary layer is introduced. A convergence study is then conducted for different discretizations and interpolation orders and subsequently a comparison is made between the numerical findings and the analytical results of the error analysis. This section is ended by comparing the numerical solution with the exact solution for three different values of the length scale g, such as 0.01L, 0.03L and 0.05L. 8.1. Model problem The beam is fixed on its left boundary and either an axial tensile load or an axial compressive load, P¯ , acts on its right boundary. It is assumed that the beam’s length is L. An axially distributed load, ¯f , is also applied along the beam. The problem can be therefore formulated as

g 2 u (4) (x) − u (2) (x) = u (0) = 0,

¯f (x) AE

u ( L ) = ε1 ,

= f (x) in  = (0, L ),

R (0) = R¯ ,

(82)

P ( L ) = P¯ .

(83)

In the applications followed, the positive sign in the prescribed boundary of the gradient of the axial displacement and in the prescribed boundary of the axial force implies that the beam is subjected to an axial tensile load, whereas the negative sign implies that the beam is subjected to an axial compressive load, respectively. Two different types of constants are selected for the beam: (i) L = 1 m, g = 0.01L, ¯f = 0, A = 3.1416 × 10−2 m2 , E = 120 × 109 Pa, (ii) L = 1 m, g = 0.01L, ¯f = 0, A = 3.1416 × 10−2 m2 , E = 120 × 109 Pa,

ε1 = 0.6, R¯ = 0 and P¯ = 103 N. ε1 = −0.4, R¯ = 0 and P¯ = −106 N.

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

131

Table 1 Numerical errors and convergence rates for the hp-version SIPG method under uniform mesh refinement (beam (i)). Method SIPG, p = 3

SIPG, p = 4

SIPG, p = 5

SIPG, p = 6

N el

|||u − u DG |||sb

4 8 16 32 64

7.2933 × 10−2

4.6786 × 10−2 2.3384 × 10−2 1.0404 × 10−2 3.4955 × 10−3

6.4047 × 10−1 1.0006 1.1684 1.5735

6.0922 × 10−4 1.6408 × 10−4 2.5168 × 10−5 2.3839 × 10−6

1.0827 1.8925 2.7048 3.4002

4 8 16 32 64

5.2299 × 10−2 2.4249 × 10−2 9.3818 × 10−3 2.5513 × 10−3 4.5000 × 10−4

1.1088 1.3700 1.8786 2.5032

7.8521 × 10−4 2.4207 × 10−4 4.0348 × 10−5 3.4526 × 10−6 1.7322 × 10−7

1.6976 2.5849 3.5467 4.3170

4 8 16 32 64

3.2812 × 10−2 1.2013 × 10−2 3.2304 × 10−3 4.8425 × 10−4 4.3799 × 10−5

1.4497 1.8947 2.7379 3.4668

4.2991 × 10−4 9.4034 × 10−5 9.4934 × 10−6 4.1778 × 10−7 1.0025 × 10−8

2.1928 3.3082 4.5061 5.3811

4 8 16 32 64

1.9586 × 10−2 5.5899 × 10−3 9.3460 × 10−4 7.4268 × 10−5 3.4069 × 10−6

1.8089 2.5804 3.6535 4.4462

2.3001 × 10−4 3.4878 × 10−5 1.9855 × 10−6 4.3819 × 10−8 5.2008 × 10−10

2.7213 4.1348 5.5018 6.3967

k1

||u − u DG ||

k2

1.2903 × 10−3

Table 2 Numerical errors and convergence rates for the hp-version NIPG method under uniform mesh refinement (beam (i)). Method

N el

|||u − u DG |||sb

k1

||u − u DG ||

k2

NIPG, p = 3

4 8 16 32 64

7.2821 × 10−2 4.6630 × 10−2 2.3287 × 10−2 1.0359e × 10−2 3.4833 × 10−3

6.4310 × 10−1 1.0017 1.1686 1.5724

1.3463 × 10−3 6.1618 × 10−4 1.6487 × 10−4 2.6755 × 10−5 3.2684 × 10−6

1.1275 1.9020 2.6234 3.0332

4 8 16 32 64

5.2235 × 10−2 2.4237 × 10−2 9.3493 × 10−3 2.5405 × 10−3 4.4839 × 10−4

1.1078 1.3743 1.8797 2.5023

7.9982 × 10−4 2.4681 × 10−4 4.3014 × 10−5 4.0172 × 10−6 2.2853 × 10−7

1.6963 2.5205 3.4205 4.1358

4 8 16 32 64

3.2883 × 10−2 1.2027 × 10−2 3.2182 × 10−3 4.8227 × 10−4 4.3643 × 10−5

1.4511 1.9019 2.7384 3.4660

4.4036 × 10−4 9.8293 × 10−5 1.0298 × 10−5 4.7992 × 10−7 1.5111 × 10−8

2.1635 3.2548 4.4234 4.9891

4 8 16 32 64

1.9712 × 10−2 5.5932 × 10−3 9.3089 × 10−4 7.3988 × 10−5 3.3962 × 10−6

1.8173 2.5870 3.6532 4.4453

2.3842 × 10−4 3.6953 × 10−5 2.1426 × 10−6 4.9047 × 10−8 6.7284 × 10−10

2.6898 4.1082 5.4490 6.1878

NIPG, p = 4

NIPG, p = 5

NIPG, p = 6

Thereafter three different types of constants are selected for the CNT: (iii) L = 10−6 m, g = 0.03L, ¯f = 0, A = 3.1416 × 10−18 m2 , E = 1040.6750 × 109 Pa,

ε1 = 10−4 , R¯ = 0 and P¯ = 10−9 N.

(iv) L = 10−7 m, g = 0.02L, ¯f = 0, A = 3.1416 × 10−18 m2 , E = 1040.6750 × 109 Pa,

ε1 = 10−4 , R¯ = 0 and P¯ = 10−9 N.

(v) L = 10−6 m, g = 0.05L, ¯f = 0, A = 3.1416 × 10−18 m2 , E = 1040.6750 × 109 Pa,

ε1 = −10−4 , R¯ = 0 and P¯ = −10−9 N.

The selected values of the stabilization constants are: C α = C β = 12, C δ = 6.

132

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

Table 3 Numerical errors and convergence rates for the hp-version SIPG method under p-refinement (beam (i)). Method SIPG, N el = 4

SIPG, N el = 8

SIPG, N el = 16

SIPG, N el = 32

SIPG, N el = 64

p

|||u − u DG |||sb

3 4 5 6

7.2933 × 10−2

5.2299 × 10−2 3.2812 × 10−2 1.9586 × 10−2

1.1560 2.0892 2.8299

7.8521 × 10−4 4.2991 × 10−4 2.3001 × 10−4

1.7265 2.6995 3.4305

3 4 5 6

4.6786 × 10−2 2.4249 × 10−2 1.2013 × 10−2 5.5899 × 10−3

2.2845 3.1480 4.1958

6.0922 × 10−4 2.4207 × 10−4 9.4034 × 10−5 3.4878 × 10−5

3.2081 4.2376 5.4398

3 4 5 6

2.3384 × 10−2 9.3818 × 10−3 3.2304 × 10−3 9.3460 × 10−4

3.1746 4.7779 6.8025

1.6408 × 10−4 4.0348 × 10−5 9.4934 × 10−6 1.9855 × 10−6

4.8764 6.4843 8.5823

3 4 5 6

1.0404 × 10−2 2.5513 × 10−3 4.8425 × 10−4 7.4268 × 10−5

4.8858 7.4470 10.2840

2.5168 × 10−5 3.4526 × 10−6 4.1778 × 10−7 4.3819 × 10−8

6.9050 9.4644 12.3680

3 4 5 6

3.4955 × 10−3 4.5000 × 10−4 4.3799 × 10−5 3.4069 × 10−6

7.1259 10.4400 14.0070

2.3839 × 10−6 1.7322 × 10−7 1.0025 × 10−8 5.2008 × 10−10

9.1140 12.7700 16.2290

k3

||u − u DG ||

k4

1.2903 × 10−3

Table 4 Numerical errors and convergence rates for the hp-version NIPG method under p-refinement (beam (i)). Method

p

|||u − u DG |||sb

k3

||u − u DG ||

k4

NIPG, N el = 4

3 4 5 6

7.2821 × 10−2 5.2235 × 10−2 3.2883 × 10−2 1.9712 × 10−2

1.1549 2.0739 2.8068

1.3463 × 10−3 7.9982 × 10−4 4.4036 × 10−4 2.3842 × 10−4

1.8100 2.6745 3.3651

3 4 5 6

4.6630 × 10−2 2.4237 × 10−2 1.2027 × 10−2 5.5932 × 10−3

2.2746 3.1401 4.1992

6.1618 × 10−4 2.4681 × 10−4 9.8293 × 10−5 3.6953 × 10−5

3.1803 4.1259 5.3659

3 4 5 6

2.3287 × 10−2 9.3493 × 10−3 3.2182 × 10−3 9.3089 × 10−4

3.1722 4.7793 6.8036

1.6487 × 10−4 4.3014 × 10−5 1.0298 × 10−5 2.1426 × 10−6

4.6705 6.4067 8.6106

3 4 5 6

1.0359 × 10−2 2.5405 × 10−3 4.8227 × 10−4 7.3988 × 10−5

4.8856 7.4464 10.2820

2.6755 × 10−5 4.0172 × 10−6 4.7992 × 10−7 4.9047 × 10−8

6.5910 9.5218 12.5100

3 4 5 6

3.4833 × 10−3 4.4839 × 10−4 4.3643 × 10−5 3.3962 × 10−6

7.1262 10.4400 14.0050

3.2684 × 10−6 2.2853 × 10−7 1.5111 × 10−8 6.7284 × 10−10

9.2477 12.1720 17.0670

NIPG, N el = 8

NIPG, N el = 16

NIPG, N el = 32

NIPG, N el = 64

8.2. Convergence study A priori error estimates stated in the Theorem 7.1.2 are illustrated by a series of numerical simulations. The asymptotic behavior of the errors of the IPDG methods (41) (i.e. the SIPG and the NIPG) in a sequence of successively finer meshes is examined for various values of the number of elements, N el , and the polynomial degree, p. Provided that the selected values of the stabilization constants are large enough, a priori error estimates will be valid for the SIPG [62]. Let us assume that the solution is smooth and discontinuous piecewise polynomials are employed. A uniform mesh having been successively refined is used and the parameter, N el = L /h, indicates the mesh refinement (i.e. h-refinement). On the other hand, the parameter, p (i.e. the degree of the approximation polynomial), indicates the p-refinement. The error is

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

133

Fig. 1. Convergence of the hp-version SIPG and NIPG in the energy seminorm of the error under h-refinement ((a) beam (i), (c) beam (ii)). Convergence of the hp-version SIPG and NIPG in the L 2 -norm of the error under h-refinement ((b) beam (i), (d) beam (ii)).

defined as e = u − u D G . The convergence rates of the IPDG methods in the energy seminorm and L 2 -norm are denoted by k1 , k2 , k3 and k4 , respectively. k1 and k2 are determined by applying the following formulas for h

k1 =

1 ln 2



ln

|||eh |||sb |||eh/2 |||sb



,

k2 =



1 ln 2

ln

||eh || , ||eh/2 ||

(84)

while k3 and k4 are determined by applying the next formulas for p

k3 =

1 ln( p + 1) − ln p



ln

|||e p |||sb |||e p +1 |||sb



, k4 =

1 ln( p + 1) − ln p

 ln

||e p || . ||e p +1 ||

(85)

Tables 1 – 4 summarize the convergence rates k1 , k2 , k3 and k4 for the gradient elastic beam (i). In Tables 1 – 2, a comparison is made between the energy seminorm and L 2 -norm of the error in the approximation to u with respect to the number of finite elements, N el , in a sequence of uniform subdivisions for polynomial degree 3  p  6. We notice that |||u − u D G |||sb and ||u − u D G || converge to zero, as N el increases and the mesh is subsequently refined for each fixed p. The error energy seminorm also exceeds the error L 2 -norm for each p and N el . Furthermore, SIPG and NIPG methods exhibit the optimal convergence O (h p −1 ) in the energy seminorm and the optimal convergence O (h p +1 ) in the L 2 -norm, as the mesh is refined for each fixed p. Tables 3 – 4 present the energy seminorm and L 2 -norm of the error in the approximation to u with respect to the polynomial degree, p, in a sequence of increasing degree of the polynomial approximation for fixed N el = 4, 8, 16, 32, 64,

134

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

Fig. 2. Convergence of the hp-version SIPG and NIPG in the energy seminorm of the error under h-refinement. (a) CNT (iii), (b) CNT (iv), (c) CNT (v).

respectively. We note that |||u − u D G |||sb and ||u − u D G || diminish and converge to zero, when p augments and the (computational) mesh is fixed. The error energy seminorm also takes greater values than the error L 2 -norm for each N el and p. Moreover, Figs. 1a, 1c, 2 display the convergence of the hp-version SIPG and NIPG energy seminorms of the error under h-refinement, respectively. To elucidate, a comparison is made between the hp-version SIPG and NIPG energy seminorms of the error in the approximation to u with the mesh parameter N el and the polynomial degree ranges from 3 to 6 for the gradient elastic beams ((i), (ii)). In addition, we make a similar comparison for the CNTs ((iii), (iv), (v)) with 4  p  6. It is clear that hp-version SIPG and NIPG error energy seminorms present indiscernible differences for each fixed p. What is more, a conclusion drawn is that hp-version SIPG and NIPG error energy seminorms converge to zero at the optimal rate O (h p −1 ) when the mesh is refined for each fixed p. Theorem 7.1.2 is therefore verified. Furthermore, in Figs. 1b, 1d, we exhibit the convergence of the hp-version SIPG and NIPG ||u − u D G || under h-enrichment, respectively. In particular, Figs. 1b, 1d exhibit the hp-version SIPG and NIPG L 2 -norms of the error in the approximation to u with the mesh parameter N el , for 3  p  6. Optimal rates of convergence are observed, as the mesh parameter, N el , increases. To be more specific, the hp-version SIPG and NIPG ||u − u D G || converge to zero at the rate O (h p +1 ), when N el augments for each fixed p. Fig. 3 exhibits the convergence of the energy seminorms and the L 2 -norms of the error with respect to p-refinement for the gradient elastic beams ((i), (ii)), as N el is fixed for the hp-version SIPG and NIPG. In Fig. 4, a similar comparison is made for the convergence of the error energy seminorms with respect to p-refinement for the CNTs ((iii), (iv), (v)). An exponential rate of convergence is expected under p-refinement, because the problem’s solution is an analytic function. In fact, the convergence plots become straight lines on a linear-log scale, as the degree of the approximate polynomials augments (Fig. 3, Fig. 4). As a result, this states an exponential convergence in p. Theorem 7.1.2 is ergo verified.

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

135

Fig. 3. Convergence of the hp-version SIPG and NIPG in the energy seminorm of the error under p-refinement ((a) beam (i), (c) beam (ii)). Convergence of the hp-version SIPG and NIPG in the L 2 -norm of the error under p-refinement ((b) beam (i), (d) beam (ii)).

Fig. 5a – 5d respectively exhibit the exact and the IPDG approximate axial displacements (BVP (82) – (83)) for various discretizations, interpolation orders and three different values of the length scale g, such as 0.01L, 0.03L and 0.05L. Both exact and approximate displacements are normalized with their corresponding maximum values. The normalization suggests that the aforementioned figures are virtually identical when either a beam ((i), (ii)) or a CNT ((iii), (iv), (v)) is used. We notice that the IPDG approximate displacements, obtained from hp-version SIPG or NIPG method, converge to the exact displacement for each length scale g that studied. What is more, we remark that both exact and approximate displacements appear to be less flexible as the length scale, g, decreases and tends to zero. Eq. (82) leads to the equation of the classical elasticity ignoring the material microstructure when the length scale, g, tends to zero. 9. Conclusions This research objective revolves around the development of the hp-version IPDGFEMs, i.e. the SIPG and the NIPG methods, for the solution of a static gradient elastic beam either in tension or compression with classic and microstructural material constants, respectively. Our research endeavor specifically focuses on performing a priori error analysis of the aforementioned methods as well as carrying out numerical simulations that verify the analytical findings. In particular, the consistency, the stability and the convergence of the IPDGFEMs are proven by the analytical results. Besides, a priori error estimates indicate an optimal convergence in the meshsize h, but a slightly suboptimal convergence in the polynomial degree p. Moreover, the analytical findings of the IPDGFEMs (i.e. Theorem 7.1.2) are verified by the numerical simulations, regarding the cases of the above structures, for various discretizations and interpolation orders. For each fixed p, the corresponding

136

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

Fig. 4. Convergence of the hp-version SIPG and NIPG in the energy seminorm of the error under p-refinement. (a) CNT (iii), (b) CNT (iv), (c) CNT (v).

numerical results show that the error energy seminorms of the SIPG and NIPG methods converge to zero at the optimal rate O (h p −1 ) when the number of elements, N el , increases and the mesh is therefore refined. Furthermore, in beams ((i), (ii)), the numerical results exhibit that the error L 2 -norms of the SIPG and NIPG methods converge to zero at the rate O (h p +1 ), as the N el augments for each fixed p. When the polynomial degree increases and N el is fixed, however, the numerical results illustrate that the error energy seminorms and the error L 2 -norms of the SIPG as well as NIPG methods demonstrate an exponential convergence in p, respectively. The IPDG approximate displacements converge to the exact displacement for various discretizations and interpolation orders as well as for each length scale, g, that studied. Both exact and approximate displacements also appear to be less flexible as the length scale, g, decreases and tends to zero. According to the conclusions drawn, the IPDGFEMs developed are considered an effective tool for dealing with the problems of the gradient elasticity theory. Acknowledgements K.G. Eptaimeros gratefully acknowledges the support of the General Secretariat for Research and Technology Award Program, administered by the Research Grants Special Account of the National Technical University of Athens (67/103100, the title of the individual project: “Research on the mechanical behavior of carbon nanotubes”). C.Chr. Koutsoumaris gratefully acknowledges the support of the General Secretariat for Research and Technology Award Mechanics Program, administered by the Research Grants Special Account of the National Technical University of Athens (Award Programs of K.A. 63/1913, 63/1947, 63/1880, 63/1882 as decision of the General Secretariat for Research and Technology 71644/28.04.2016).

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

137

Fig. 5. Normalized exact displacement, hp-version SIPG and NIPG normalized approximate displacements for g = 0.01L , 0.03L , 0.05L. (a) p = 3, N el = 64. (b) p = 4, N el = 32. (c) p = 5, N el = 16. (d) p = 6, N el = 12.

References [1] N. Aravas, Plane-strain problems for a class of gradient elasticity models-a stress function approach, J. Elast. 104 (1–2) (2011) 45–70. [2] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (4) (1982) 742–760. [3] D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2001) 1749–1779. [4] I. Babuška, M. Suri, The optimal convergence rate of the p-version of the finite element method, SIAM J. Numer. Anal. 24 (4) (1987) 750–776. [5] G.A. Baker, Finite element methods for elliptic equations using nonconforming elements, Math. Comput. 31 (137) (1977) 45–59. ˇ ´ M. Canadija, R. Luciano, F. Marotti de Sciarra, Application of gradient elasticity to armchair carbon nanotubes: size effects and [6] R. Barretta, M. Brˇcic, constitutive parameters assessment, Eur. J. Mech. A, Solids 65 (2017) 1–13. [7] G. Becker, L. Noels, A fracture framework for Euler-Bernoulli beams based on a full discontinuous Galerkin formulation/extrinsic cohesive law combination, Int. J. Numer. Methods Eng. 85 (2011) 1227–1251. [8] D. Braess, Finite Elements. Theory, Fast Solvers and Application in Solid Mechanics, Cambridge University Press, United Kingdom, 1997. [9] E. Burman, A. Ern, I. Mozolevski, B. Stamm, The symmetric discontinuous Galerkin method does not need stabilization in 1D for polynomial orders p  2, C. R. Acad. Sci. Paris I 345 (2007) 599–602. [10] Z. Cai, X. Ye, S. Zhang, Discontinuous Galerkin finite element methods for interface problems: a priori and a posteriori error estimations, SIAM J. Numer. Anal. 49 (2011) 1761–1787. [11] C.F. Cornwell, L.T. Wille, Elastic properties of single-walled carbon nanotubes in compression, Solid State Commun. 101 (8) (1997) 555–558. [12] J. Douglas, T. Dupont Jr., Interior penalty procedures for elliptic and parabolic Galerkin methods, in: Computing Methods in Applied Sciences, Springer, 1976, pp. 207–216. [13] W.H. Duan, C.M. Wang, Y.Y. Zhang, Calibration of nonlocal scaling effect parameter for free vibration of carbon nanotubes by molecular dynamics, J. Appl. Phys. 101 (2) (2007) 024305.

138

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

[14] G. Engel, K. Garikipati, T.J.R. Hughes, M.G. Larson, L. Mazzei, R.L. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Comput. Methods Appl. Mech. Eng. 191 (2002) 3669–3750. [15] K.G. Eptaimeros, Continuous interior penalty finite element method for a 6th-order bending gradient elastic (nano)beam equation, AIP Conf. Proc. 1978 (2018) 470031, https://doi.org/10.1063/1.5044101. [16] K.G. Eptaimeros, C.C. Koutsoumaris, G.J. Tsamasphyros, Nonlocal integral approach to the dynamical response of nanobeams, Int. J. Mech. Sci. 115–116 (2016) 68–80. [17] K.G. Eptaimeros, C.C. Koutsoumaris, I.T. Dernikas, Th. Zisis, Dynamical response of an embedded nanobeam by using nonlocal integral stress models, Composites, Part B, Eng. 150 (2018) 255–268. [18] A.T. Eyck, A. Lew, Discontinuous Galerkin methods for non-linear elasticity, Int. J. Numer. Methods Eng. 67 (9) (2006) 1204–1243. [19] X. Feng, T. Lewis, M. Neilan, Discontinuous Galerkin finite element differential calculus and applications to numerical solutions of linear and nonlinear partial differential equations, J. Comput. Appl. Math. 299 (2016) 68–91. [20] H.G. Georgiadis, C.G. Grentzelou, Energy theorems and the integral in dipolar gradient elasticity, Int. J. Solids Struct. 43 (2006) 5690–5712. [21] E.H. Georgoulis, P. Houston, Discontinuous Galerkin methods for the biharmonic problem, IMA J. Numer. Anal. 29 (2009) 573–594. [22] E.H. Georgoulis, E. Hall, J.M. Melenk, On the suboptimality of the p-version interior penalty discontinuous Galerkin method, J. Sci. Comput. 42 (2010) 54–67. [23] A.E. Giannakopoulos, K. Stamoulis, Structural analysis of gradient elastic components, Int. J. Solids Struct. 44 (2007) 3440–3451. [24] R.F. Gibson, E.O. Ayorinde, Y.F. Wen, Vibrations of carbon nanotubes and their composites: a review, Compos. Sci. Technol. 67 (1) (2007) 1–28. [25] Z. Hammouch, T. Mekkaoui, P. Agarwal, Optical solitons for the Calogero-Bogoyavlenskii-Schiff equation in (2+1) dimensions with time-fractional conformable derivative, Eur. Phys. J. Plus 133 (7) (2018) 248. [26] V.M. Harik, Ranges of applicability for the continuum beam model in the mechanics of carbon nanotubes and nanorods, Solid State Commun. 120 (7–8) (2001) 331–335. [27] P. Houston, C. Schwab, E. Süli, Discontinuous hp-finite element methods for advection-diffusion problems, SIAM J. Numer. Anal. 39 (6) (2002) 2133–2163. [28] C.C. Koutsoumaris, K.G. Eptaimeros, A research into bi-Helmholtz type of nonlocal elasticity and a direct approach to Eringen’s nonlocal integral model in a finite body, Acta Mech. 229 (2018) 3629–3649. [29] C.C. Koutsoumaris, K.G. Eptaimeros, The gradient beam: a confrontation between the analytical closed type and numerical type solution, AIP Conf. Proc. 1978 (2018) 470032, https://doi.org/10.1063/1.5044102. [30] C.C. Koutsoumaris, K.G. Eptaimeros, G.J. Tsamasphyros, A different approach to Eringen’s nonlocal integral stress model with applications for beams, Int. J. Solids Struct. 112 (2017) 222–238. [31] A. Krishnan, E. Dujardin, T.W. Ebbesen, P.N. Yianilos, M.M.J. Treacy, Young’s modulus of single-walled nanotubes, Phys. Rev. B 58 (20) (1998) 14013–14019. [32] A. Kumar, S. Kumar, Residual power series method for fractional Burger types equations, Nonlinear Eng. 5 (4) (2016) 235–244. [33] A. Kumar, S. Kumar, S.P. Yan, Residual power series method for fractional diffusion equations, Fundam. Inform. 151 (1–4) (2017) 213–230. [34] D.C.C. Lam, F. Yang, A.C.M. Chong, J. Wang, P. Tong, Experiments and theory in strain gradient elasticity, J. Mech. Phys. Solids 51 (2003) 1477–1508. [35] M.G. Larson, A.J. Niklasson, Analysis of a non symmetric discontinuous Galerkin method for elliptic problems: stability and energy error estimates, SIAM J. Numer. Anal. 42 (1) (2004) 252–264. [36] K.-T. Lau, C. Gu, D. Hui, A critical review on nanotube and nanotube/nanoclay related polymer composite materials, Composites, Part B, Eng. 37 (6) (2006) 425–436. [37] M. Lazar, G.A. Maugin, A note on line forces in gradient elasticity, Mech. Res. Commun. 33 (2006) 674–680. [38] K.A. Lazopoulos, A.K. Lazopoulos, Bending and buckling of thin strain gradient elastic beams, Eur. J. Mech. A, Solids 29 (2010) 837–843. [39] J. Li, J.M. Melenk, B. Wohlmuth, J. Zou, Optimal a priori estimates for higher order finite elements for elliptic interface problems, J. Appl. Numer. Math. 60 (2010) 19–37. [40] D.M. Manias, T.K. Papathanasiou, S.I. Markolefas, E.E. Theotokoglou, Analysis of a gradient-elastic beam on Winkler foundation and applications to nano-structure modelling, Eur. J. Mech. A, Solids 56 (2016) 45–58. [41] S.I. Markolefas, D.A. Tsouvalas, G.I. Tsamasphyros, Theoretical analysis of a class of mixed, C 0 -continuity formulations for general dipolar gradient elasticity boundary value problems, Int. J. Solids Struct. 44 (2007) 546–572. [42] S.I. Markolefas, D.A. Tsouvalas, G.I. Tsamasphyros, Some C 0 -continuous mixed formulations for general dipolar linear gradient elasticity boundary value problems and the associated energy theorems, Int. J. Solids Struct. 45 (2008) 3255–3281. [43] R.D. Mindlin, Micro-structure in linear elasticity, Arch. Ration. Mech. Anal. 16 (1964) 51–78. [44] R.D. Mindlin, Second gradient of strain and surface tension in linear elasticity, Int. J. Solids Struct. 1 (4) (1965) 417–438. [45] R.D. Mindlin, N.N. Eshel, On first-gradient theories in linear elasticity, Int. J. Solids Struct. 4 (1968) 109–124. [46] R.D. Mindlin, H.F. Tiersten, Effects of couple-stresses in linear elasticity, Arch. Ration. Mech. Anal. 11 (1) (1962) 415–448. [47] V.F. Morales-Delgado, J.F. Gómez-Aguilar, K.M. Saad, M.A. Khan, P. Agarwal, Analytic solution for oxygen diffusion from capillary to tissues involving external force effects: a fractional calculus approach, Physica A 523 (2019) 48–65. [48] I. Mozolevski, E. Süli, A priori error analysis for the hp-version of the discontinuous Galerkin finite element method for the biharmonic equation, Comput. Methods Appl. Math. 3 (4) (2003) 596–607. [49] J.T. Oden, I. Babuska, C.E. Baumann, A discontinuous hp-finite element method for diffusion problems, J. Comput. Phys. 146 (2) (1998) 491–519. [50] A. Pantano, D.M. Parks, M.C. Boyce, Mechanics of deformation of single- and multi-wall carbon nanotubes, J. Mech. Phys. Solids 52 (4) (2004) 789–821. [51] S. Papargyri-Beskou, D.E. Beskos, Static analysis of gradient elastic bars, beams, plates and shells, Open Mech. J. 4 (2010) 65–73. [52] S. Papargyri-Beskou, K.G. Tsepoura, D. Polyzos, D.E. Beskos, Bending and stability analysis of gradient elastic beams, Int. J. Solids Struct. 40 (2003) 385–400. [53] C. Polizzotto, A gradient elasticity theory for second-grade materials and higher-order inertia, Int. J. Solids Struct. 49 (2012) 2121–2137. [54] P. Poncharal, Z.L. Wang, D. Ugarte, W.A. De Heer, Electrostatic deflections and electromechanical resonances of carbon nanotubes, Science 283 (5407) (1999) 1513–1516. [55] W.H. Reed, T.R. Hill, Triangular Mesh Methods for the Neutron Transport Equation, Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, Los Alamos, 1973. [56] B. Rivière, M.F. Wheeler, A discontinuous Galerkin method applied to non-linear parabolic equations, Lect. Notes Comput. Sci. Eng. 11 (1999) 231–244. [57] B. Rivière, M.F. Wheeler, V. Girault, Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems I, Comput. Geosci. 3 (1999) 337–360. [58] B. Rivière, M.F. Wheeler, V. Girault, A priori error estimates for finite element methods based on discontinuous approximation spaces for elliptic problems, SIAM J. Numer. Anal. 39 (3) (2001) 902–931. [59] A. Romkes, S. Prudhomme, J.T. Oden, A priori error analyses of stabilized discontinuous Galerkin method, Comput. Math. Appl. 46 (2003) 1289–1311. [60] K.M. Saad, O.S. Iyiola, P. Agarwal, An effective homotopy analysis method to solve the cubic isothermal auto-catalytic chemical system, AIMS Math. 3 (1) (2018) 183–194.

K.G. Eptaimeros et al. / Applied Numerical Mathematics 144 (2019) 118–139

139

[61] C. Schwab, p- and hp-Finite Element Methods. Theory and Applications to Solid and Fluid Mechanics, Numerical Mathematics and Scientific Computation, Clarendon/Oxford University Press, New York, 1998. [62] E. Süli, I. Mozolevski, hp-version interior penalty finite element DGFEMs for the biharmonic equation, Comput. Methods Appl. Mech. Eng. 196 (2007) 1851–1863. [63] E.T. Thostenson, Z. Ren, T.W. Chou, Advances in the science and technology of carbon nanotubes and their composites: a review, Compos. Sci. Technol. 61 (2001) 1899–1912. [64] R.A. Toupin, Elastic materials with couple-stress, Arch. Ration. Mech. Anal. 11 (1962) 385–414. [65] G.I. Tsamasphyros, S.I. Markolefas, D.A. Tsouvalas, Convergence and performance of the h- and p-extensions with mixed finite element C 0 -continuity formulations, for tension and buckling of a gradient elastic beam, Int. J. Solids Struct. 44 (2007) 5056–5074. [66] K.G. Tsepoura, S. Papargyri-Beskou, D. Polyzos, D.E. Beskos, Static and dynamic analysis of gradient elastic bars in tension, Arch. Appl. Mech. 72 (2002) 483–497. [67] Q. Wang, K.M. Liew, Application of nonlocal continuum mechanics to static analysis of micro- and nano-structures, Phys. Lett. A 363 (3) (2007) 236–242. [68] Q. Wang, C.M. Wang, The constitutive relation and small scale parameter of nonlocal continuum mechanics for modelling carbon nanotubes, Nanotechnology 18 (7) (2007) 075702. [69] Y.B. Wang, X.W. Zhu, H.H. Dai, Exact solutions for the static bending of Euler-Bernoulli beams using Eringen’s two phase local/nonlocal model, AIP Adv. 6 (8) (2016) 085114. [70] Z. Wang, Q. Tang, W. Guo, Y. Cheng, Sparse grid discontinuous Galerkin methods for high-dimensional elliptic equations, J. Comput. Phys. 314 (2016) 244–263. [71] L. Wei, X. Zhang, S. Kumar, A. Yildirim, A numerical study based on an implicit fully discrete local discontinuous Galerkin method for the time-fractional coupled Schrödinger system, Comput. Math. Appl. 64 (8) (2012) 2603–2615. [72] L. Wei, Y. He, A. Yildirim, S. Kumar, Numerical algorithm based on an implicit fully discrete local discontinuous Galerkin method for the time-fractional KdV-Burgers-Kuramoto equation, J. Appl. Math. Mech. (Z. Angew. Math. Mech.) 93 (1) (2013) 14–28. [73] M.F. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM J. Numer. Anal. 15 (1) (1978) 152–161. [74] E.W. Wong, P.E. Sheehan, C.M. Lieber, Nanobeam mechanics: elasticity, strength, and toughness of nanorods and nanotubes, Science 277 (5334) (1997) 1971–1975. [75] M.F. Yu, B.S. Files, S. Arepalli, R.S. Ruoff, Tensile loading of ropes of single wall carbon nanotubes and their mechanical properties, Phys. Rev. Lett. 84 (24) (2000) 5552–5555.