Novel differential quadrature element method for higher order strain gradient elasticity theory

# Novel differential quadrature element method for higher order strain gradient elasticity theory

Md Ishaquddin Corresponding author: E-mail address: ishaquddinmd@iisc.ac.in S.Gopalakrishnan E-mail address: krishnan@iisc.ac.in; Phone: +91-80-22932048
###### Abstract

In this paper, we propose a novel and efficient differential quadrature element based on Lagrange interpolation to solve a sixth order partial differential equations encountered in non-classical beam theories. These non-classical theories render displacement, slope and curvature as degrees of freedom for an Euler-Bernoulli beam. A generalize scheme is presented herein to implementation the multi-degrees degrees of freedom associated with these non-classical theories in a simplified and efficient way. The proposed element has displacement as the only degree of freedom in the domain, whereas, at the boundaries it has displacement, slope and curvature. Further, we extend this methodology and formulate two novel versions of plate element for gradient elasticity theory. In the first version, Lagrange interpolation is assumed in and directions and the second version is based on mixed interpolation, with Lagrange interpolation in direction and Hermite interpolation in direction. The procedure to compute the modified weighting coefficients by incorporating the classical and non-classical boundary conditions is explained. The efficiency of the proposed elements is demonstrated through numerical examples on static analysis of gradient elastic beams and plates for different boundary conditions.

Keywords: Differential quadrature element, gradient elasticity, sixth order pde, weighting coefficients, non-classical, Lagrange interpolation, Lagrange-Hermite mixed interpolation

Department of Aerospace Engineering,Indian Institute of Science Bengaluru 560012, India

## 1.0 Introduction

The differential quadrature method (DQM) is an efficient numerical tool for the solution of initial and boundary value problems. This technique was first introduced by Bellman et al.[1]. The major challenge in the application of DQM to structural problems is the implementation of multiple boundary conditions. To resolve this issue, many efficient and improved DQ schemes were developed in recent years. A comprehensive survey on the DQM development and the recent contributions to this field can be found in [2, 3, 4]. Bert et al. proposed a - method to solve problems in structural mechanics [5, 6]. Subsequent developments in this technique led to the application of this methodology to variety of problems [7, 8, 9, 10, 11, 12]. As this method could not be generalized and had limitations on accuracy, an alternative scheme was proposed which accounted boundary conditions during the formulation of weighting coefficients [13]. However, this method also could not be generalized and worked only for few specific types of boundary conditions and structures with constrained regular edges. Further improvement in these schemes were discussed in the articles by Du et al. [14, 15].

To address the difficulties in applying the DQ method for structures with discontinuous loading and geometry, Striz et al. [16, 17] developed a quadrature element method (QEM), however, due to the use of - technique the scope of this method was limited. Later, Wang et al. [18] and Chen et al.[12] proposed a differential quadrature element method (DQEM) which assumes the slope as an independent degree of freedom at the boundary. The main advantage of this method is only one grid point is required to represent the multiple degrees of freedom at the boundary. Further improvement in this field led to the development of a new method called the generalized differential quadrature rule (GDQR) [19, 20, 21]. Following this many researchers applied DEQM and GDQR techniques to variety of structural problems [22, 23, 24, 25, 26, 27, 28, 29].

In the above DQEM and GDQR techniques, Hermit interpolation functions were used to determine the weighting coefficients. In contrast, Wang et al.[28] employed the weighting coefficients based on Lagrange interpolation functions. The research inclination in the aforesaid publications was towards the solution of fourth order partial differential equations which governs the problems related to classical beam and plate theories. The DQ solution for the sixth and eighth order differential equations using GDQR technique with Hermite interpolation function was reported by Wu et al. [30, 31]. They have demonstrated the capability for structural and fluid mechanics problems. Recently, Wang et al. [32] proposed a new differential quadrature element based on Hermite interpolation to solve a sixth order partial differential equation governing the non-local Euler-Bernoulli beam. In their study, they have computed the frequencies for various combination of boundary conditions.

The classical continuum theories are effective for macro scale modelling of structural elements, neverthless, they lack efficiency to model the nano scale systems. These classical theories are governed by fourth order partial differential equations. To overcome this difficulty, several scale-dependent non-classical continuum theories are reported in the literature [33, 34, 35, 36, 37]. These non-classical continuum theories are enriched versions of classical continuum theories incorporating higher order terms in the constitutive relations. These higher order terms consists of stress and strain gradients accompanied with intrinsic length scale parameters which account for scale effects[38, 39, 40, 41]. One such class of gradient elasticity theory is the simplified theory by Mindlin et al.[33], with one gradient elastic modulus and two classical constants for structural applications. This simplified theory was used by many researchers to study static, dynamic and buckling behaviour of gradient elastic beams [42, 43, 44, 45, 46] and plates [47, 48, 49], by deriving the analytical solutions. Recently, Pegios et.al [50] developed a finite element model for gradient elastic Euler-Bernoulli beam and conducted static and stability analysis. The numerical solution of 2-D and 3-D gradient elastic structural problems using finite element and boundary element methods can be found in [51].

In this paper, we propose for the first time a novel differential quadrature beam element based on Lagrange interpolation to solve a sixth order partial differential equation associated with gradient elastic Euler-Bernoulli beam theory. Further, we extend this methodology and formulate two novel versions of plate element for gradient elastic Kirchhoff plate theory. In the first version, the Lagrange interpolations are used in both and direction, and later, mixed-interpolations are used with Lagrange interpolations in direction and continuous Hermite in direction. A novel way to impose the classical and non-classical boundary conditions for gradient elastic beam and plate elements are presented. A new procedure to compute the higher order weighting coefficients for the proposed elements are explained in detail. The efficiency and the performance of the elements are established through numerical examples.

## 1 Strain gradient elasticity theory

Mindlin’s [33] strain gradient micro-elasticity theory with two classical and one non-classical material constants is consider in the present study. The two classical material constants are Lam constants and the non-classical one is related to intrinsic bulk length . In what follows, the theoretical basis required to formulate the differential quadrature beam and plate element for gradient elasticity theory are presented. Further, the classical and non-classical boundary conditions associated with the gradient elastic Euler-Bernoulli beam and Kirchhoff plate are discussed.

### 1.1 Gradient elastic Euler-Bernoulli beam

The stress-strain relations for 1-D gradient elastic theory are defined as [52, 44]

 τ =2με+λtrεI ς =g2[2μ∇ε+λ∇(% trε)I] (1)

where, , are Lam constants. is the Laplacian operator and I is the unit tensor. , denotes Cauchy and higher order stress respectively, and () are the classical strain and its trace which are expressed in terms of displacement vector w as:

 ε=12(∇w+w∇),trε=∇w (2)

It follows from the above equations the constitutive relations for an Euler-Bernoulli gradient beam can be stated as

 τx=Eεx,ςx=g2ε′x,εx=−z∂2w(x,t)∂x2 (3)

For the above state of stress and strain the strain energy expression in terms of displacement can be written as

 Ub=12∫L0EI[(w′′)2+g2(w′′′)2]dx (4)

The potential energy of the applied load is given by

 (5)

The total potential energy of the beam is given by

 Πb=Ub+Vb (6)

where, , and are the Young’s modulus, area, moment of inertia, respectively. and are the transverse load and displacement of the beam. , and are shear force, bending moment and higher order moment acting on the beam.

Using the principle of minimum potential energy [55]:

 δΠb=δ(Ub+Vb)=0 (7)

and performing integration by parts, we get the governing equation for a gradient elastic Euler-Bernoulli beam as

 EI(w′‵′−g2w‵′′)+qb=0 (8)

and the associated boundary conditions are:

Classical :

 V =EI[w′′′−g2w‵′]=0orw=0,atx=(0,L) M =EI[w′′−g2w′‵′]=0orw′=0,atx=(0,L) (9)

Non-classical :

 ¯M =[g2EIw′′′]=0% orw′′=0,atx=(0,L) (10)

The list of classical and non-classical boundary conditions employed in the present study for a gradient elastic Euler-Bernoulli beam are as follows

Simply supported :
classical :   ,    non-classical :   at

Clamped :
classical :   ,    non-classical :   at

Next, the constitutive relations, governing equation and the associated classical and non-classical boundary conditions for a gradient Kirchhoff plate are presented.

### 1.2 Gradient elastic Kirchhoff plate

The strain-displacement relations for a Kirchhoff’s plate theory can be defined as [56]

 εxx=−z¯wxx,εyy=−z¯wyy,γxy=2εxy=−2z¯wxy (11)

where, is transverse displacement of the plate. The stress-strain relations for a gradient elastic Kirchhoff plate is given by [52, 37]:

Classical:

 τxx= E1−ν2(εxx+νεyy) τyy= E1−ν2(εyy+νεxx) (12) τxy= E1+νεxy

Non-classical:

 ςxx= g2E1−ν2∇2(εxx+νεyy) ςyy= g2E1−ν2∇2(εyy+νεxx) (13) ςxy= g2E1+ν∇2εxy

where, , ,, are the classical Cauchy stresses and ,, denotes higher order stresses related to gradient elasticity.

The strain energy for a gradient elastic Kirchhoff plate is given by [37, 49]:

 Up=Ucl+Usg (14)

where, and are the classical and gradient elastic strain energy given by

 Ucl=12D∫∫A[¯w2xx+¯w2yy+2¯w2xy+2ν(¯wxx¯wyy−¯w2xy)]dxdy (15)
 Usg= 12g2D∫∫A[¯w2xxx+¯w2yyy+3(¯w2xyy+¯w2xxy)+2ν(¯wxyy¯wxxx+ ¯wxxy¯wyyy−¯w2xyy−¯w2xxy]dxdy (16)

here, .

The potential energy of the external load is defined as

 Wp=∫∫Aq0(x,y)¯wdxdy+∫sVn¯wds−∫sMn¯w′(n)ds−∫s¯Mn¯w′′(n)ds (17)

where is the transverse load on the plate. and are the normal and tangential to the boundary point corresponding to and coordinates axis.

Using the principle of minimum potential energy and performing integration by parts over the area, we obtain the governing equation for a gradient elastic Kirchhoff plate as

 D∇4¯w−g2D∇6¯w−q0=0 (18)

where,

 ∇4=∂4¯w∂x4+∂4¯w∂y4+2∂4¯w∂x2∂y2, ∇6=∂6¯w∂x6+∂6¯w∂y6+3∂6¯w∂x4∂y2+3∂6¯w∂x2∂y4

the associated boundary conditions for the plate with domain defined over (), (), are listed below.

Classical boundary conditions :

 or ¯w=0,atx=(0,lx)
 or ¯w=0,aty=(0,ly) (19)
 Mx=−D(∂2¯w∂x2+ν∂2¯w∂y2)+g2D[∂4¯w∂x4+ν∂4¯w∂y4+(3−ν)∂4¯w∂x2∂y2]=0 or ¯wx=0,atx=(0,lx) My=−D(∂2¯w∂y2+ν∂2¯w∂x2)+g2D[∂4¯w∂y4+ν∂4¯w∂x4+(3−ν)∂4¯w∂x2∂y2]=0 or ¯wy=0,aty=(0,ly)

Non-classical boundary conditions :

 ¯Mx=−g2D(∂3¯w∂x3+ν∂3¯w∂x∂y2)=0or¯wxx=0,atx=(0,lx) ¯My=−g2D(∂3¯w∂y3+ν∂3¯w∂y∂x2)=0or¯wyy=0,aty=(0,ly) (21)

The concentrated force at the free corner is given by

 R=2D(1−ν)∂2¯w∂x∂y−2g2D(1−ν)[∂4¯w∂x3∂y+∂4¯w∂y3∂x] (22)

where and are the length and width of the plate. ,  are the shear force, ,  are the bending moment and ,  are the higher order moment. is the concentrated force at the free corner.

The different boundary conditions employed in the present study for a gradient elastic Kirchhoff plate are:

Simply supported edge :
or      at
or      at

Clamped edge :
at
at

Free edge :
at
at

Free corner :

## 2 Differential quadrature elements for gradient elasticity theory

In this section, first, we formulate a differential quadrature element based on Lagrangian interpolation function for 1-D gradient elastic Euler-Bernoulli beam. Next, we develop two new versions of gradient plate element with different choice of interpolation functions. The grid employed in the present study is unequal Gauss–Lobatto–Chebyshev points given by

 zi=12[1−cos(i−1)πN−1] (23)

where is the number of grid points and are the coordinates of the grid. For the plate analysis is employed.

### 2.1 Differential quadrature element for gradient Euler-Bernoulli beam

The th order derivative of the deflection at location for a N-node 1-D beam element is assumed as

 wni(x,t)=N∑j=1Lnj(x)wj (24)

are Lagrangian interpolation functions in co-ordinate. The Lagrange interpolation functions can be defined as[4, 3],

 Lj(x)=β(x)β(xj)=N∏k=1(k≠j)(x−xk)(xj−xk) (25)

where

The first order derivative of the above shape functions can be written as

 Aij=L′j(xi)⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩∏Nk=1(k≠i,j)(xi−xk)/∏Nk=1(k≠j)=(xj−xk)(i≠j)∑Nk=1(k≠i)1(xi−xk) (26)

The conventional higher order weighting coefficients are computed as

 Bij=N∑k=1AikAkj,Cij=N∑k=1BikAkj,Dij=N∑k=1BikBkj,(i,j=1,2,...,N) (27)

here, , and are weighting coefficients for second, third, and fourth order derivative, respectively.

Let us consider a N-node gradient Euler-Bernoulli beam element as shown in the Figure 1.

Each interior node has displacement as the only degree of freedom, and the boundary nodes has 3 degrees of freedom , , . These extra boundary degrees of freedom related to slope and curvature are introduced in to the formulation through modifying the conventional weighting coefficients. The new displacement vector now includes the slope and curvature as additional degrees of freedom at the element boundaries as: . The modified weighting coefficient matrices accounting for slope and curvature degrees of freedom at the boundaries are derived as follows:

First order derivative matrix:

 ¯Aij=⎧⎨⎩Aij(i,j=1,2,⋯,N)0(i=1,2,⋯,N;j=N+1,⋯,N+4) (28)

Second order derivative matrix:

 ¯Bij=⎧⎨⎩Bij(i=2,3,⋯,N−1;j=1,2,⋯,N)0(i=2,3,⋯,N−1;j=N+1,⋯,N+4) (29)
 ¯Bij=N−1∑k=2AikAkj(i=1,N;j=1,2,⋯,N) (30)
 ¯Bi(N+1)=Ai1;¯Bi(N+2)=AiN(i=1,N) (31)

Third order derivative matrix:

 ¯Cij=⎧⎪⎨⎪⎩∑Nj=1∑Nk=1¯BikAkj(i=2,3,...,N−1)0(i=2,3,⋯,N−1;j=N+1,⋯,N+4) (32)
 ¯Cij=N−1∑k=2BikAkj (i=1,N;j=1,2,⋯,N) (33)
 ¯Ci(N+3)=Ai1;¯Ci(N+4)=AiN(i=1,N) (34)

Fourth order derivative matrix:

 ¯Dij=N+4∑j=1N∑k=1Bik¯Bkj(i=1,2,...,N) (35)

Fifth order derivative matrix:

 Vij=⎧⎨⎩¯Dij(i=2,3,⋯,N−1;j=1,2,⋯,N)0(i=2,3,⋯,N−1;j=N+1,⋯,N+4) (36)
 Vij=N−1∑k=2BikBkj(i=1,N;j=1,2,⋯,N) Vi(N+3)=Bi1;Vi(N+4)=BiN(i=1,N) (37)
 ¯Eij=N∑k=1AikVkj(i=1,2,...,N;j=1,2,...,N+4) (38)

Sixth order derivative matrix:

 ¯Fij=N∑k=1BikVkj(i=1,2,...,N;j=1,2,...,N+4) (39)

Here, , , , , and are first to sixth order modified weighting coefficients matrices, respectively. Using the above Equations (28)-(39), the governing differential equation (8), in terms of the differential quadrature at inner grid points is written as

 EIN+4∑j=1¯Dijwbj−g2EIN+4∑j=1¯Fijwbj=qb(xi)(i=2,3,...,N−1) (40)

The boundary forces given by Equations (1.1)-(10), in terms of differential quadrature are expressed as

Shear force:

 Vi=EIN+4∑j=1¯Cijwbj−g2EIN+4∑j=1¯Eijwbj(i=1,N) (41)

Bending moment:

 Mi=EIN+4∑j=1¯Bijwbj−g2EIN+4∑j=1¯Dijwbj(i=1,N) (42)

Higher order moment:

 ¯Mi=g2EIN+4∑j=1¯Cijwbj (i=1,N) (43)

here and correspond to the left support and right support of the beam, respectively.

Once the boundary conditions in Equations (1.1)-(10) are applied, we get the following system of equations in the matrix form as

 ⎡⎢ ⎢ ⎢⎣kbb−kbdkdb−kdd⎤⎥ ⎥ ⎥⎦⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ΔbΔd⎫⎪ ⎪ ⎪⎬⎪ ⎪ ⎪⎭=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩fbfd⎫⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪⎭ (44)

where the subscript and indicates the boundary and domain of the beam. , and , are the boundary and domain forces and displacements of the beam, respectively. Now expressing the system of equations in terms of domain dofs , we get

 [kdd−kdbk−1bbkbd]{Δd}={fd−kdbk−1bbfb} (45)

The solution of the above system of equations renders the displacements at the domain nodes of the beam element. The boundary displacements are computed from Equation (44), and forces are computed from the Equations (41)-(43).

### 2.2 Differential quadrature element for gradient elastic plates

Here, we present two versions of novel differential quadrature element for a gradient elastic Kirchhoff plate. First, the differential quadrature element based on Lagrange interpolation in and direction is formulated. Next, the differential quadrature element based on Lagrange-Hermite mixed interpolation, with Lagrangian interpolation is direction and Hermite interpolation assumed in direction is presented. Similar to the beam elements discussed in the previous section, the plate element has only displacement as degrees of freedom in the domain and at the plate edges it has 3 degrees of freedom , , or , depending upon the edge. At the corners the element has five degrees of freedom , , , and . The new displacement vector now includes the slope and curvature as additional degrees of freedom at the element boundaries given by: , where .

A differential quadrature gradient plate element for a grid is shown in the Figure 2. Here, are the number of grid points in and directions, respectively. It can be seen that the element has three degrees of freedom on each edge and five degrees of freedom at the corners. In the figure, the slope and curvature dofs with first subscript and , correspond to the edges and , respectively. Similarly, the slope and curvature dofs with second subscript and , correspond to the edges and , respectively.

#### 2.2.1 Lagrange interpolation based differential quadrature element for gradient elastic plates

The deflection for a node differential quadrature rectangular plate element is assumed as

 w(x,y,t)=Nx∑p=1Ny∑q=1Lp(x)Lq(y)~wpq(t) (46)

where is the nodal deflection vector and and are the Lagrange interpolation functions in and directions, respectively. The slope and curvature degrees of freedom at the element boundaries are accounted while computing the weighting coefficients of higher order derivatives as discussed in section 2.1. Using the 1-D Lagrange interpolation functions derived in the section 2.1, the governing differential equation (18) in the differential quadrature syntax at inner grid points can be expressed as

 D[Nx+4∑r=1¯Dxpr~wrs+2Nx+4∑q=1Ny+4∑r=1¯Bxpq¯Bysr~wqr+Ny+4∑r=1¯Dysr~wpr]− g2D[Nx+4∑r=1¯Fxpr~wrs+3Nx+4∑q=1Ny+4∑r=1¯Dxpq¯Bysr~wqr+3Nx+4∑q=1Ny+4∑r=1¯Bxpq¯Dysr~wqr+ Ny+4∑r=1¯Fysr~wpr]=¯qo(xp,ys) (p=2,3,...,Nx−1;s=2,3,...,Ny−1) (47)

The boundary forces, Equation (1.2)-(22), in terms of differential quadrature form are expressed as
Shear force:

 Vx= −D(Nx+4∑r=1¯Cxpr~wrs+(2−ν)Nx+4∑q=1Ny+4∑r=1¯Axpq¯Bysr~wqr)+ g2D[Nx+4∑r=1¯Expr~wrs+(3−ν)Nx+4∑q=1Ny+4∑r=1¯Axpq¯Dysr~wqr+3Nx+4∑q=1Ny+4∑r=1¯Cxpq¯Bysr~wqr] (p=1,Nx;s=2,3,...,Ny−1) (48)
 Vy= −D(Ny+4∑r=1¯Cysr~wpr+(2−ν)Nx+4∑q=1Ny+4∑r=1¯Bxpq¯Aysr~wqr)+ g2D[Ny+4∑r=1¯Eysr~wpr+(3−ν)Nx+4∑q=1Ny+4∑r=1¯Dxpq¯Aysr~wqr+3Nx+4∑q=1Ny+4∑r=1¯Bxpq¯Cysr~wqr] (p=2,3,...,Nx−1;s=1,Ny) (49)

Bending moment:

 Mx= −D(Nx+4∑r=1¯Bxpr~wrs+νNy+4∑r=1¯Bysr~wpr)+ g2D[Nx+4∑r=1¯Dxpr~wrs+νNy+4∑r=1¯Dysr~wpr+(3−ν)Nx+4∑q=1Ny+4∑r=1¯Bxpq¯Bysr~wqr] (p=1,Nx;s=2,3,...,Ny−1) (50)
 My= −D(Ny+4∑r=1¯Bysr~wpr+νNx+4∑r=1¯Bxpr~wrs)+ g2D[Ny+4∑r=1¯Dysr~wpr+νNx+4∑r=1¯Dxpr~wrs+(3−ν)Nx+4∑q=1Ny+4∑r=1¯Bxpq¯Bysr~wqr] (p=2,3,...,Nx−1;s=1,Ny) (51)

Higher order moment:

 ¯Mx= −g2D(Nx+4∑r=1¯Cxpr~wrs+νNx+4∑q=1Ny+4∑r=1¯Axpq¯Bysr~wqr) (p=1,Nx;s=2,3,...,Ny−1) (52)
 ¯My= −g2D(Ny+4∑r=1¯Cysr~wpr+νNx+4∑q=1Ny+4∑r=1¯Bxpq¯Aysr~wqr