Strain gradient elasticity within the symmetric BEM formulation

A BSTRACT . The symmetric Galerkin Boundary Element Method is used to address a class of strain gradient elastic materials featured by a free energy function of the (classical) strain and of its (first) gradient. With respect to the classical elasticity, additional response variables intervene, such as the normal derivative of the displacements on the boundary, and the work-coniugate double tractions. The fundamental solutions - featuring a fourth order partial differential equations (PDEs) system - exhibit singularities which in 2D may be of the order 4 1/ r . New techniques are developed, which allow the elimination of most of the latter singularities. The present paper has to be intended as a research communication wherein some results, being elaborated within a more general paper [1], are reported.


INTRODUCTION
fter the pioneering work of Mindlin [2], theories of strain gradient elasticity have become very popular, particularly within the domain of nano-technologies, that is, for problems where the ratio surface/volume tends to become very large and there is a need to introduce at least one internal length.However the model introduced by Mindlin and then improved by Mindlin et al. [3] and Wu [4] leads to an excessive number of material coefficients, which at the best for isotropic materials reduce to the number of seven.In the early 1990's, Aifantis [5] introduced a signified material model of strain gradient elasticity which requires only three material coefficients, that is, the Lame' constants and one length scale parameter.The latter model was then developed further following the so-called Form II format given by Mindlin et al. [3], that is a theory centered on the existence of a free energy function of the (classical) strain and of its first gradient, which leads to the generation of symmetric stress fields (see Askes et al. [6] for historical details about the latter formulations and its applications).Formulations in the boundary element method based on the strain gradient elasticity were pioneered by Polyzos et al. [7], Karlis et al. [8], who provide a collocation BEM formulation where the simplified constitutive equation by Aifantis [5], has been adopted.In latter papers only the fundamental solutions used in the collocation approach to BEM are provided.In the case of the symmetric formulation of the BEM, Somigliana Identities (SIs) for the tractions and for the double tractions are also needed.These new SIs are necessary in order to get, through the process of modeling and weighing, a solving equation system having symmetric operators.In Polizzotto et al. [1] all the set of fundamental solutions is derived starting from the displacement fundamental solution given in [7,8].

A
The symmetric formulation is motivated by the high efficiency achieved within classical elasticity by the method in [9] with regard to the techniques used to eliminate the singularities of the fundamental solutions, the evaluation of the coefficients of the solving system and the computational procedures characterized by great implementation simplicity.This has already led to the birth of the computer code Karnak sGbem [10] operating in the classical elasticity.The objective of this paper is to experiment new techniques and procedures that, applied in the context of strain gradient elastic materials, may permit one to obtain the related solving system.

NOTATION
s a rule, a compact notation is used, with vectors and tensors denoted by bold-face symbols.The dot  , the colon products : and :: indicate the simple and higher index contraction.The following notations are also utilized: the gradient operator: grad ( ) the Laplacian operator : the component of the normal derivation of u is denoted as the quantities with bar superscripts indicate know quantities.Other symbols will be specified in the text at their first appearance.

BASIC RELATIONS IN 2D
he class of strain gradient elastic materials herein considered is featured by the following strain elastic energy 2 1 : : that is a function of the symmetric 2nd order strains u and of the strain gradient (1)   ε ε, where  is the internal length which relates the microstructure with macrostructure.Eq.(1) provides the stress tensors work-coniugate of ( 0 )  ε and (1)  ε respectively where E is the classic isotropic elasticity tensor.By the principle of the virtual work, the indefinite equilibrium equation and the following boundary conditions prove to be where the total stresses σ , the tractions t , the double tractions r and the force at the corner P are so defined ( )( ) : In eq.(4b), the symbol ( ) ( ) denotes the tangent gradient on the boundary and ( ) s K    n is the curvature on the boundary having normal n ; in eq.(4d) s denotes a unit vector tangent to the two boundary portions convergent in the corner C and the brackets    indicate that the enclosed quantity is the difference between the related values taken on two sides of the corner C. For a more exhaustive understanding of the jump    to see [2].
The constitutive equation relating the total stress σ to the strain is A T THE FUNDAMENTAL SOLUTIONS he introduction of eq.( 5) in (3a) valid in infinity domain   allows to express the latter in terms of displacements only [7,8] and to determine the fundamental solution of the displacements that for 2D solids proves to be (1 ) where  are Bessel functions of the second kind and of order 0 and 2, respectively, and   r ξ x is the distance between effect ξ and cause x points.One can note that the singularity of the fundamental solution uu G does not depend on r .Indeed, by performing the limit  r 0 one obtains (1 ) where  is the Eulero constant and the fundamental solution for isotropic gradient elasticity shows singularity of type ln( )  .In eq.( 6) the limit of uu G for 0   gives the classic isotropic elasticity solution (Kelvin) with singularity of type ln( ) r .
Table 1: Fundamental solutions for strain gradient elasticity and relative singularities.
By the fundamental solution uu G in (6), taking in account eqs.(4b,c,d) and using the well-known procedure given in [11], it is possible, by exploiting the known properties of symmetry of the fundamental solutions, to derive the entire tableau provided in Table 1.The fundamental solutions hk G of Tab. 1 are characterized by two subscripts: the first indicates the effect at ξ , i.e. displacement for h u  , traction for h t  , displacement normal derivative for h g  , double traction for h r  , corner  , a unit imposed strain for k   .To consider the way in which the fundamental solutions included in the single    and double        brackets to see [2].While the fundamental solutions represent the response in a point ξ of the unlimited domain caused by unit kinematical or mechanical action imposed at x , the response to the distributed actions is provided by the SIs that in [1] are obtained by a generalization of the Betti theorem for the gradient elastic materials.

A CASE STUDY
s a premise, it should be noted that the theory of strain gradient should be applied to those bodies where the ratio between surface and volume in 3D and between boundary and area in 2D is high.In these cases, the body should be considered as composed of a nucleus, a crust and its boundary.The example we want to meet will show that on the basis of this theory the behavior of the cortex has different characteristics from the classical behavior of the nucleus, and this depends on the presence of new variables such as double traction r and normal derivative g of the displacements.Obviously, the body will suffer the effects due to the presence of these new variables.From computational point of view, the fundamental solutions present in the Tab. 1 show various orders of singularity, up to 4 1/ r in the column related to u , up to 3 1/ r in that related to C u ; furthermore the singularities present in the column C u present an order lower than that of u , relatively to the same effect.In order to investigate the techniques useful to remove the singularities from the blocks of coefficients, a simple application in two-dimensional space is shown, but at this stage involving the study of fundamental solutions having the maximum order of singularity equal to 2 1/ r .It has to be remembered that in the symmetric BEM the coefficients are obtained by a double integration, the first (inner) regarding the modeling of the cause through appropriate shape functions and the second (external) regarding the effect weighing through shape functions, but dual in the energetic sense.For this purpose, the example shown concerns a plate completely constrained on the boundary, subjected to the simple distortions u , C u and distortions of higher order / n    g u .This example is also used to develop the techniques of rigid motion that could be used in other applications to compute the coefficient blocks in which the techniques available are not sufficient to suppress the residual singularities.For the plate of Fig. 1a the boundary conditions are Let us proceed to write the SIs on the boundary and, following the compact notation proposed by Panzeca et al. [9], the latter equations become Let us now discretize the plate boundary into boundary elements (Fig. 1b) and introduce the modeling of the quantities on the boundary as functions of the nodal variables through suitable matrices of shape functions h N , with , , , h t r u g  ; ; ; Since the vector C u collects nodal quantities to not be modeled, the following identity has to be imposed In eqs.(12a-d) the quantities T , R , U , G have the meaning of nodal quantities and precisely: force and double traction as unknowns, displacements and normal derivatives of displacements as known quantities.To evaluate the vector G , because the discontinuity of the normal derivative at the corner C, a double node belonging respectively to the two portions of boundary afferent to C has to be considered.
Introducing the modeling (12a-d) and the relation (13) into eqs.(11a,b)and performing the weighing second Galerkin, by using the shape functions in dual form, one writes and, following the compact notation proposed by Panzeca et al. [9], the latter equations become Since some components of the vector U and of the vector ; Finally, the solving system is rewritten in compact way with obvious meaning of symbols with K symmetric flexibility matrix of the system, X vector of nodal unknowns, u L and g L load vectors due to U e G respectively, being u g   L L L the total load vector.In order to evaluate the coefficients of the equation system, the numerical techniques to remove the singularities and the rigid motion strategy have to be employed.

Observation about rigid motions
Let the plate of Fig. 2 be subjected to a) a rigid motion of translation obtained as the sum of constant displacement fields x u and y u (Fig. 2a), b) a rigid motion of rotation having wideness  around any point (Fig. 2b), c) a rigid motion of roto-translation obtained as a combination of a) and b).In all these cases the displacement and normal derivative of the displacement fields are entirely known both in the domain and on the boundary of the plate.Particularly -for the case of rigid motion of translation (Fig. 2a) , 0 , 0 with c1 and 2 c being the wideness of the displacements; -for the case of rigid motion of rotation having width f (Fig. 2b) ,c,d,e

Observations about the load vector as a consequence of the rigid motion
1) If the solid is subject to a rigid motion of translation (Fig. 2a) having assigned values of the displacements c1 and 2 c , one has  U 0 and  G 0 but the total load vector has to be ( ) and as a consequence the solution vector X has to be null.
2) If the solid is subject to a rigid motion of rotation (Fig. 2a) having assigned rotation vector f  φ , in the generic node i one has with i r vector distance between the instantaneous center of rotation and the node i , i n e i s being respectively the unit normal and tangent vectors to the boundary elements on which the node i lies, but the total load vector has to be because the solution vector X has to be null.
3) If the solid is subjected to a combination of the rigid motions 1) and 2), the condition u g    L L L 0 has to be always valid.

Example of coefficient calculation in presence of singularity
The greater difficulties consist in removing the singularities when the cause is focused at the corner.The simple observation that the vector u regards all the nodes of the boundary, including the vector C u at the corner, allows to eliminate the higher order singularities after the first integration.Indeed, the use of the distribution theory permits to evaluate the effect on the boundary without limit operations and removing the strong singularities through the summing of effects, whereas the remaining ones are eliminated through the second integration.The singularities present in the fundamental solutions used for the plate analysis require shape functions of type ( 0 )  C .In order to show how operate to remove singularity in the coefficients of the solving system, we examine the load vector and, suppressing for simplicity the subscript yy, through the process of weighing one obtains The following contributions bb g , bc g bC g are provided distinguishing between the singular and the regular part (sing ) ( reg ) ( reg ) The sum of the effects in terms of distributions bb g , bc g bC g eliminates the strong singularity 1/ r , while the weak singularity ( ) Log r present in the bb g distribution is eliminated in the weighing operation by integration by parts.Eqn.(31) gives the coefficient in closed form.

Numerical results
Let us analyze the two-dimensional plate of Fig. 1 having the following physical and mechanical characteristics: The plate is subjected to a linear displacement distribution u , including the nodal displacement C u , and to the linear displacement normal derivatives / n    g u , all imposed on the boundary.The displacement and displacement gradient fields imposed on the boundary have the distribution shown in Fig. 4a,b.The field response is derived by the SIs after the solution being obtained.For the example here considered, characterized by the presence of eight nodes on the boundary and four corner (Fig. 1b), the know nodal vector U has sixteen components, whereas C U defined by I 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 U H U U 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 shows eight components; the vector G collects twenty-four components, due to the presence of the double node belonging respectively to the two portions of boundary afferent to the corners,.Two load conditions are considered:  First load condition For this condition the vectors of know nodal quantities U , C U and G are     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 For this condition the vectors of known nodal quantities U , C U and G are       0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 Tractions x T , y T , double tractions x R , y R , obtained by analysis, are summarized in Tab. 2 and Tab. 3.
Fig. 5 and Fig. 6 show the fields of displacements x u , of standard strain ( 0 ) xx  , of high-order strain (1)   2 2 ( 0 )  Comments The diagrams of Fig. 5,6 show how the distribution of normal derivatives x g imposed on the boundary influence the trend of the displacement field, in particular the slope near the boundary.
-in the first load conditition the imposed value 0 x g  involves zero slope at 1 x   for the displacement field and the related field of the standard strain ( 0 ) xx  show null values at 1 x   ; the field of high-order strain show high values at 1 x   and values close to zero in the centre of the plate -in the second load condition the imposed value 1 x g  provides at 1 x   the same slope of the displacement field inside the plate; the standard strain ( 0 ) xx  show constant unit value, whereas the field of the high-order strain is null and the total strain coincides with the standard strain.This is the case in which the imposed value for x g gets the same solution valid for isotropic classic elasticity.


, stress for h   ; whereas the second subscript indicates -through a T work-conjugate rule -the cause applied at x , i.e. a unit concentrated force for k u  , a unit surface relative displacement for k t  , a unit concentrated double force for k g  , a unit higher order surface distortion for k r  , a unit corner force for

Figure 1 :
Figure 1: Plate: a) geometry and constraints; b) discretization into boundary elements and linear modelling of the boundary quantities.
10b) In the SIs:  the bar superscripts characterize known quantities;  the apex CPV indicates that the corresponding integral is evaluated as Cauchy Principal Value to which the related free term 1 ( ) 2  is added;  C u is a vector containing only the known displacement of the corner; are, respectively, the displacement and normal derivative of displacement distribution on the boundary, computed as the difference between the response to the values of displacement imposed on the two portions of boundary afferent to the corner.Introducing the latter SIs given in eqs.(10a,b) into the boundary conditions eqs.(9a,b) one has: 1

Figure 4 :
Figure 4: Distribution on the boundary a) x u and b) x g .

Figure 5 :
Figure 5: Field mechanical characteristics regarding the first load condition.

Figure 6 :
Figure 6: Field mechanical characteristics regarding the second load condition.
it is opportune to introduce the following relation

Table 2 :
Results for the first load condition.

Table 3 :
Results for the second load condition.