On Finite Element Computations of Contact Problems in Micropolar Elasticity

Within the linear micropolar elasticity we discuss the development of new finite element and its implementation in commercial software. Here we implement the developed 8-node hybrid isoparametric element into ABAQUS and perform solutions of contact problems. We consider the contact of polymeric stamp modelled within the micropolar elasticity with an elastic substrate. The peculiarities of modelling of contact problems with a user defined finite element in ABAQUS are discussed. The provided comparison of solutions obtained within the micropolar and classical elasticity shows the influence of micropolar properties on stress concentration in the vicinity of contact area.


Introduction
Nowadays the interest grows to further development and application extended models of continuum mechanics in order to model micro-and nanostructured materials with complex inner structure.The basic idea of enhancement of classic Cauchy continuum model is to add additional fields describing additional degrees of freedom into constitutive equations or/and consider higher-order gradients of deformations.Among these generalized models there are the surface elasticity, micropolar or Cosserat continua, microstretched and micromorphic media, media with internal variables, gradient elasticity, and so forth.In particular, the micropolar model [1][2][3] proposed by Cosserat brothers more than hundred years ago found applications for modelling such materials as porous solids [4,5], bones [6][7][8], masonries [9,10], beam lattices [11], and other composite materials; see [1,[12][13][14][15] and reference therein.The micropolar elasticity possesses the description of size-effects and may be useful for description of the microstructured solids such as foam, bones, powders, and soils.In particular, the influence of micropolar properties may be important for the analysis of the stress concentration near holes and notches or in the vicinity of contact area.Within the Cosserat continuum model the translations and rotations determine the kinematics of the medium and the stress and couple stress tensors are introduced.The micropolar properties of material may be important near singularities or to describe observed experimentally size-effect [16][17][18][19].Let us also note that for such complex media even more general models of continuum mechanics such as gradient elasticity may be useful; see, for example, [20][21][22][23][24]. Application of extended models of continuum for modelling of such structures as open-cell foams, beam lattices, and pantographic systems is motivated by their complex inner structure.For example, in the classic beam theory moments play an important role and their influence is inherited by the homogenized models.Let us also note that for structured and nonhomogeneous beams used as structural elements for foam or lattice description some extended models were proposed in [25][26][27][28][29][30][31].For rods and beams there are some phenomena as warping of crosssection, instabilities, and sensibility to imperfections are also important; see, for example, [32][33][34][35][36][37][38][39][40][41].Considering all of these phenomena may lead to rather complex models of homogenized media.From this point of view the micropolar elasticity may be treated as first step towards modelling of 2 Advances in Materials Science and Engineering microstructured materials and their contact.Further extensions of the indentation problems can be performed based on the other enriched models of continuum; see, for example, strain gradients models [42,43].
Let us note that the effective solution of boundary value problems of micropolar elasticity as well as of other enhanced models requires advanced numerical code such as the finite element method.The generalized models of continua require usually more computational efforts than the classical elasticity since there exist more degrees of freedom.For the micropolar continuum we use an isogeometric analysis [44,45] as efficient FEM strategy which together with a hybrid mixed formulation was applied for generalized continua and structures; see, for example, [46][47][48][49][50].For the moment commercial FEM software gives one the possibility to use user defined elements and user defined procedure for implementation of nonstandard material models.Here we developed new finite element and implement it in ABAQUS.
The paper is organized as follows.In Section 2 we present the basic equations of the linear micropolar elasticity.The equilibrium equations, static and kinematic boundary conditions, and constitutive equations are given.Here we are restricted by isotropic case using the Voigt notation.In Section 3 we discuss the finite element modelling for the micropolar solids.Finally, in Section 4 we present the solution of contact problem for micropolar parabolic stamp and an elastic thick plate.

Basic Equations of the Micropolar Elasticity
Following [1][2][3] we recall here the basic equations of the linear micropolar elasticity of isotropic solids.The kinematic of a micropolar solid is described by two fields, that is, the field of translations   and the field of rotations   .The latter is responsible for the description of moment interactions of the material particles.Hereinafter the Latin indices take on value 1, 2, or 3 and we use the Einstein summation rule over repeating indices.The equilibrium equations take the form where where   is the unit vector of external normal to the boundary  =   ∪   ,   and   are external forces and couples, and  0  and  0  are given on   surface fields of translations and rotations, respectively.Obviously, other mixed boundary conditions can be introduced.
In what follows we are restricting ourselves by isotropic case.For a linear isotropic micropolar solid the constitutive equations are where   is the Kronecker symbol, and , , , , , and  are the elastic moduli.In (3) we also introduced the stretch tensor   and the wryness tensor   given by formulas Using the Voigt notation modified for the micropolar elasticity and introducing the stress and moment stress vectors with stretch and wryness vectors by the formulas we represent the constitutive equations in the following unified form: with 18 × 18 symmetric stiffness matrix [C], where A and B are three-diagonal symmetric 9 × 9 matrices given by ) ) ) , ) .
Let us note that the homogeneous micropolar model can be derived on the base of the passage from heterogeneous classical (Cauchy) continuum and on using passage from inhomogeneous micropolar (Cosserat) material; see [52][53][54][55].Analysis of general constitutive equations of anisotropic micropolar solids on the base of the material symmetry groups is performed in [2,56,57].
From the experimental point of view it is better to use another set of material parameters [4,5,8,14] listed in Table 1.
The boundary value problems of the micropolar isotropic elasticity contain also the boundary value problems of classical elasticity as a special case.This coincidence may be used for verification of the developed code.The first way to reduce the problem is to assume  = 0 or  = 0. Indeed, if the coupling modulus  = 0 the stress tensor becomes symmetric and takes the form Thus, (1) 1 transforms to the classic equilibrium equation of the linear elasticity which is independent of rotations whereas (1) 2 includes only rotations.Let us also note that (1) 2 has similar form to (1) 1 but with different elastic moduli.So in this case the problem is decoupled and translations can be determined independently of microrotations.
The second reduction is possible if one assumes the microrotations to be fixed,   = 0.In this case (1) 1 again coincides with classic equilibrium equation with Lame's moduli  and .In what follows we use both methods for verification and validation of the developed code.

Implementation of Micropolar Elasticity into FEM Software
Efficient solution of large boundary value problems requires application of an advanced software.In this research ABAQUS commercial program has been extended by the implementation of the user element (UEL) to solve the micropolar elasticity problems.Special code written in Fortran is linked with ABAQUS software allowing the user to practice all the ABAQUS features without paying an attention to their numerical implementation.From the practical point of view the most important features are creation of sophisticated geometry, application of loads and boundary conditions, applications of constrains and contact conditions, generation of 3D meshes, using the material library, and, the most important feature, being very effective solver.
During the solver execution UEL procedure is called twice for each Gaussian point in every element.In the first call the element stiffness should be provided by UEL procedure.Very often user element procedure requires calling UELMAT code (user material procedure) necessary to obtain the relation between stress and strain increments.Another call of UEL procedure is necessary to compute residual forces, element nodal forces resulting from element stresses, which is essential in the convergence monitoring during solving nonlinear problems.Designing of own finite elements is an ambiguous task recommended to advanced users only.It should be mentioned that UEL procedures should be very carefully tested and validated.There are some important disadvantages of using user elements in ABAQUS program.First of all, ABAQUS does not recognize the shape of the element; the element is represented by the set of nodes only.For example, 2D four nodes' element can be defined as a truss structure, as a frame made up of beams, or as a quadrilateral flat element.In each of the mentioned cases different type of loads can be applied; for example, the pressure can be applied to quadrilateral element only, the bending moment can be applied to beam, and so forth.The visualization of obtained results is another matter.The visualization is not possible without detailed description of an element topology; the set of nodes does not provide this information.That is why only nodal displacements of the user elements can be displayed in ABAQUS postprocessor.Displaying strains and stresses requires developing own graphical programs or using uncomfortable techniques with "ghost" meshes made of ordinary ABAQUS elements constrained to meshes consisting of user elements.
The typical call of UEL procedure is presented below: The most important parameters are amatrx (element stiffness matrix) and rhs (right-hand side vector of residuals necessary to check the convergence of computation).In the most simple case of static linear analysis the rhs parameter can be omitted because the convergence (equilibrium equations at nodes) is achieved in the first iteration.The meaning of other parameters is described in ABAQUS documentation.UEL procedure can be slightly simplified if it contains constitutive properties; the call of UELMAT procedure is not necessary, therefore.
The whole procedure is typical for the isoparametric finite element formulation in 3D problems with the exception of considered degrees of freedom (in micropolar elasticity there are three displacement components and additionally three microrotations) and strain and stress measures (strain and stress tensors contain more components and are not symmetric).In UEL implementation of micropolar elasticity the same shape functions are used for displacement and microrotations where   ,   , and   are the nodal coordinates in the local coordinate system.For the Eringen strain measure the matrix of the shape functions derivatives is  In the isoparametric element the shape functions ( 9) are applied for the mapping of local coordinates , ,  to the coordinates , ,  of the global Cartesian coordinate system Calculation of derivatives  , requires application of chain rule; that is, or where J is the Jacobian matrix which can be found from ( 9) and (11).
The stiffness matrix of developed element is where the constitutive stiffness matrix C is presented in the previous section.The integration in ( 15) is performed as the loop over the Gaussian points.The 2 × 2 × 2 Gaussian quadrature is used.
Short description of micropolar elasticity 8-node isoparametric element implementation consists of several steps.The user element procedure based on ( 9) and ( 11)-( 15) is shortly summarized below.
(i) Calculation of stiffness matrix of 8-node isoparametric element for micropolar elasticity.In the loop over the Gaussian points (the are 8 Gaussian points) (1) Find the shape functions at each Gaussian point and their derivatives with respect to natural coordinates (2) Find the Jacobian matrix ( 14); numerically compute its inverse and determinant (3) Find the shape functions derivatives (13) with respect to Cartesian coordinates , ,  (4) Find the matrix of shape functions derivatives D (the relation between strain components and nodal displacements and microrotations) (5) Find the constitutive matric C (the relation between stress and strain tensors) (6) Compute D  CD (for the 2-point Gaussian quadrature all weights are equal to ones) and add it as the contribution to the stiffness matrix which uses the determinant of Jacobian matrix as the multiplier (ii) Additional computations (1) Find the element nodal forces resulting from element stresses and subtract them from the external nodal forces in order to compute residuals (necessary to check the rate of the solution convergence) (2) Make other computations; for example, update the strain energy

Parabolic Stamp Indentation
Using the presented above finite element we analyzed few 3D static problems for solids with certain singularities such as notch, hole, or small contact area for the contact problem of two solids [17].For comparison of solutions with the linear elasticity we used the fact that the micropolar elasticity provides the same solutions as classic elasticity if we assume specific form of microrotations (  = 0) or for specific set of elastic moduli that is without coupling ( = 0).Similar results obtained by the commercial software using the classical theory of elasticity and the results acquired by use of user element while reduced number of material data is used prove that the solution obtained by the UEL procedure is reliable.In this research the Hertzian contact between a parabolic stamp and a half space is considered.For classical Hertz theory of an elastic contact we refer, for example, to [58, pp. 102-14], where formulas for stress distribution are also given.By parabolic stamp we mean a solid with surface in the form of an elliptic paraboloid that is a surface given by equation  2  =  2 +  2 , where  is a parameter, and , , and  are the Cartesian coordinates, respectively; see Figure 1.The main problem arising in this type of analysis is that ABAQUS program does not recognize the volumes and faces of user elements (user element is represented by a set of nodes).The only possible contact formulation is the contact between the surface (ordinary finite elements implemented in ABAQUS) and the set of nodes (UEL elements), therefore.The solution is less precise than the solution of the same problem considering the contact of two surfaces although reasonable results can be achieved for sufficiently dense meshes.
The contact of a parabolic stamp made of polystyrene foam with the flat elastic plate is considered.The foam is modelled as a micropolar material.For the foam we use following material data [5]: shear modulus  = 104 MPa, Poisson's ratio ] = 0.4, characteristic length for torsion   = 0.62 mm, characteristic length for bending   = 0.33 mm, coupling number  = 0.04, and polar ratio Ψ = 1.5.The average cell size of the foam is about 1 mm; it correlates with characteristic lengths.As in the case of homogenization of regular lattice structures [11], the foam cell size plays a role of a characteristic length parameter, and   and   depend on it.For the plate we consider the classic elastic material with the same shear modulus and Poisson's ratio as for foam.
In Figure 1 the finite element mesh is presented which is dense enough in the contact zone in order to obtain the results confirming Hertz's formulas.The distribution of von Mises stress obtained by ABAQUS within classical linear theory of elasticity is shown in Figure 2. Very similar results (distribution and magnitude) are obtained using developed micropolar user elements when microrotations are fixed.This demonstrates excellent convergence of the distribution of stress and stress magnitudes using the developed element to the classical solutions.So the developed UEL procedure is reliable.
In Figure 3 the distribution of couple stress   in the vicinity of the contact zone is shown; here  is the vertical axis.Couple stress appears only in the stamp (modelled as the micropolar medium) near the contact zone.As expected the distribution of   is antisymmetric.The magnitude of couple stress (all components) is several orders smaller than the magnitude of Cauchy stresses.Thus, the influence of the micropolar properties is restricted only to zones with singularities such as sudden change in normal stresses in contact area and free surface.On the other hand, the micropolar elasticity gives the possibility to capture sizeeffect, observed for microstructured solids such as foams as well as in the case of nanoindentation.All solved benchmark tests for  ̸ = 0 considered here have shown that the difference between classical elasticity and micropolar one is essential in an area of size of the characteristic length.In such areas the stresses are not symmetric and there is a redistribution of the deformation energy among translational and rotational degrees of freedom; one can expect rotations of microstructural elements such as foam struts in foams or grains in granular media.So for local behavior such as in the case of contact problem with small area and deformations near other types of stress concentrators and at the microand nanoscales the influence of micropolar properties can be important.

Conclusions
We discussed finite element approach adopted to the linear micropolar elasticity in order to model microstructured solids such as porous materials and beam lattices.The new 8-node hybrid micropolar isoparametric element and its implementation in ABAQUS are presented.Here we analyzed the contact problem between two elastic solids.Comparison of solutions based on classical and micropolar elasticity is carefully discussed.Numerical tests have shown that couple stress appears almost in the vicinity of contact zone.

Figure 1 :Figure 2 :
Figure 1: Fine element mesh in the contact problem.

Figure 3 :
Figure 3: The distribution of couple stress   .
and   are the stress and couple stress tensors, respectively,   is the Levi-Civita third-order tensor, and   and   are external forces and couples.