# On Finite Element Computations of Contact Problems in Micropolar Elasticity

Authors: Victor A. Eremeyev; Andrzej Skrzat; Feliks Stachowicz
Journal: Advances in Materials Science and Engineering
 (2016)

---

## Abstract

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.


---

## Body

## 1. 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–3] proposed by Cosserat brothers more than hundred years ago found applications for modelling such materials as porous solids [4, 5], bones [6–8], masonries [9, 10], beam lattices [11], and other composite materials; see [1, 12–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–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–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–31]. For rods and beams there are some phenomena as warping of cross-section, instabilities, and sensibility to imperfections are also important; see, for example, [32–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 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–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 Section2 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.


## 2. Basic Equations of the Micropolar Elasticity

Following [1–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 ui and the field of rotations θi. 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 (1)tij,i+fj=0,mij,i+ejmntmn+cj=0,where tij and mij are the stress and couple stress tensors, respectively, eijk is the Levi-Civita third-order tensor, and fj and ci are external forces and couples. Unlike in classical (Cauchy-type) continua, the tensors tij and mij are not symmetric. Equations (1) constitute the local balance of momentum and moment of momentum, respectively.The static and kinematic boundary conditions have the following form:(2)nitijAt=φj,nimijAt=ηj,uiAu=ui0,θiAu=θi0,where ni is the unit vector of external normal to the boundary A=At∪Au, φj and ηj are external forces and couples, and ui0 and θi0 are given on Au 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(3)tij=λεkkδij+μ+κεij+μεji,mij=ακkkδij+βκji+γκij,where δij is the Kronecker symbol, and λ, μ, κ, α, β, and γ are the elastic moduli. In (3) we also introduced the stretch tensor εij and the wryness tensor κij given by formulas (4)εij=uj,i-eijnθn,κij=θj,i.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(5)T=t11,t22,t33,t12,t21,t23,t32,t13,t31T,M=m11,m22,m33,m12,m21,m23,m32,m13,m31T,E=ε11,ε22,ε33,ε12,ε21,ε23,ε32,ε13,ε31T,K=κ11,κ22,κ33,κ12,κ21,κ23,κ32,κ13,κ31T,we represent the constitutive equations in the following unified form:(6)σM=CεM,C=A00B,σM=TM,εM=EKwith 18×18 symmetric stiffness matrix [C], where A and B are three-diagonal symmetric 9×9 matrices given by (7)A=λ+2μ+κλλ000000λλ+2μ+κλ000000λλλ+2μ+κ000000000μ+κμ0000000μμ+κ000000000μ+κμ0000000μμ+κ000000000μ+κμ0000000μμ+κ,B=α+β+γαα000000αα+β+γα000000ααα+β+γ000000000γβ0000000βγ000000000γβ0000000βγ000000000γβ0000000βγ.Unlike classical elasticity where the stiffness matrix in the Voigt notations is a symmetric 6×6 matrix here we have larger matrices. The values of used micropolar elastic moduli can be found from direct experiments provided in [4, 5, 8, 14] or using homogenization technique [6, 7, 9–11, 51]. 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–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.Table 1
Micropolar material parameters.
Description Symbol Formula Shear modulus [MPa] G G = 2 μ + κ 2  Poisson’s ratio [—] ν ν = λ 2 λ + 2 μ + κ  Coupling number [—] N N 2 = κ 2 ( μ + κ )  Characteristic length (torsion) [m] l t l t = β + γ 2 μ + κ Characteristic length (bending) [m] l b l b = γ 2 μ + κ  Polar ratio [—] Ψ Ψ = β + γ α + β + γ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 N=0. Indeed, if the coupling modulus κ=0 the stress tensor becomes symmetric and takes the form (8)tij=λuk,kδij+μui,j+uj,i.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,θj=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.


## 3. 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:subroutine uelmat(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,1
props,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,2
dtime,kstep,kinc,jelem,params,ndload,jdltyp,adlmag,3
predef,npredf,lflags,mlvarx,ddlmag,mdload,pnewdt,4
jprops,njpro,period,materiallib)The most important parameters areamatrx (element stiffness matrix) andrhs (right-hand side vector of residuals necessary to check the convergence of computation). In the most simple case of static linear analysis therhs 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(9)Niξ,η,ζ=181+ξξi1+ηηi1+ζζi,where ξi, ηi, and ζi are the nodal coordinates in the local coordinate system. For the Eringen strain measure the matrix of the shape functions derivatives is (10)ε11ε22ε33ε12ε21ε23ε32ε13ε31κ11κ22κ33κ12κ21κ23κ32κ13κ31=f1,100000⋯f8,1000000f1,20000⋯0f8,2000000f1,3000⋯00f8,30000f1,1000f1⋯0f8,1000f8f1,20000-f1⋯f8,20000-f800f1,2f100⋯00f8,2f8000f1,30-f100⋯0f8,30-f80000f1,10-f10⋯00f8,10-f80f1,3000f10⋯f8,3000f80000f1,100⋯000f8,1000000f1,20⋯0000f8,2000000f1,3⋯00000f8,3000f1,200⋯000f8,2000000f1,10⋯0000f8,100000f1,30⋯0000f8,3000000f1,2⋯00000f8,2000f1,300⋯000f8,30000000f1,1⋯00000f8,1u11u21u31θ11θ21θ31⋮u18u28u38θ18θ28θ38,where uij and θij are the ith displacement or microrotation components in the jth node, respectively.In the isoparametric element the shape functions (9) are applied for the mapping of local coordinates ξ, η, ζ to the coordinates x, y, z of the global Cartesian coordinate system(11)x=∑Niξ,η,ζxi;y=∑Niξ,η,ζyi;z=∑Niξ,η,ζzi.Calculation of derivativesfi,j requires application of chain rule; that is,(12)Ni,xNi,yNi,z=ξ,xη,xζ,xξ,yη,yζ,yξ,zη,zζ,zNi,ξNi,ηNi,ζor (13)Ni,xNi,yNi,z=J-1Ni,ξNi,ηNi,ζ,where J is the Jacobian matrix (14)J=x,ξy,ξz,xx,ηy,ηy,ηx,ζz,ζz,ζ,which can be found from (9) and (11).The stiffness matrix of developed element is(15)K=∫VDTCDdV,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 x, y, z(4)
Find the matrix of shape functions derivativesD (the relation between strain components and nodal displacements and microrotations)(5)
Find the constitutive matricC (the relation between stress and strain tensors)(6)
ComputeDTCD (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


## 4. 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 (θi=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 p2y=x2+z2, where p is a parameter, and x, y, and z 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.Figure 1
Fine element mesh in the contact problem.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 G=104 MPa, Poisson’s ratio ν=0.4, characteristic length for torsion lt=0.62 mm, characteristic length for bending lb=0.33 mm, coupling number N=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 lb and lt 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 Figure1 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.Figure 2
Von Mises stress [Pa]: commercial software, classical theory of elasticity.In Figure3 the distribution of couple stress myy in the vicinity of the contact zone is shown; here y 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 myy 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 size-effect, 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 micro- and nanoscales the influence of micropolar properties can be important.Figure 3
The distribution of couple stressmyy.


## 5. 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.



---

*Source: 9675604-2016-12-18.xml*