# A multi-scale computational method including contact for the analysis of damage in composite materials

Georgios A Drosopoulos, Peter Wriggers, Georgios E Stavroulakis

► To cite this version:

Georgios A Drosopoulos, Peter Wriggers, Georgios E Stavroulakis. A multi-scale computational method including contact for the analysis of damage in composite materials. Computational Materials Science, 2014, 95, pp.522 - 535. 10.1016/j.commatsci.2014.08.004. hal-01403838

HAL Id: hal-01403838

https://hal.science/hal-01403838

Submitted on 27 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A multi-scale computational method including contact for the analysis of damage in composite materials

Georgios A. Drosopoulos<sup>a</sup>, Peter Wriggers<sup>a</sup>, Georgios E. Stavroulakis<sup>b</sup>

In this article, a multi-scale computational homogenization scheme is proposed for the study of composite materials. A classical unilateral contact law has been incorporated in the microscopic level, for the investigation of the contact between the constitutive materials. The either-or decision resulting from the contact-no contact condition in the microscopic scale, makes the problem non-linear. This change in the contact state of the microscopic level, is taken into account by the proposed approach. Debonding between the matrix and the surrounding fibers and its impact on the macroscopic structure, are depicted. In addition, a change in the direction of the macroscopic load during analysis, results in a non-linear behavior due to the alteration of the microscopic contact state. The distribution of the displacement jump is influenced in this case, as well.

## 1. Introduction

In the present work a multi-scale, computational homogenization scheme is proposed for the study of composite materials. A classical unilateral contact law is applied in the microscopic interfaces to capture the contact behavior of the constitutive materials.

Several different, analytical and numerical approaches have been proposed in the past for the investigation of complex, non-linear, heterogeneous structures, like composites. Analytical/mathematical methods, like asymptotic homogenization [1], can be more accurate in the description of the micro structure, for relatively simple microscopic patterns and constitutive laws. On the other hand, numerical methods may be used for the simulation of complex microscopic geometries, over a statistically defined representative amount of material [2].

Numerical/computational homogenization can be extended to cover several non-linear effects, like contact, debonding, damage and plasticity [3]. According to numerical homogenization, a unit cell is explicitly solved and the resulting average quantities are then used for the determination of the parameters of a macroscopic constitutive law [4,5].

From another point of view, multi-scale computational homogenization incorporates a concurrent analysis of both the macro and the microstructure in a nested multi-scale approach, [6–14]. Within this method, the macroscopic constitutive behavior is determined during simulation, after solving the microscopic problem and transferring the necessary information on the macroscopic scale. This approach, which is generally called FE², offers the flexibility of simulating complex microstructural patterns, with every kind of non-linearity. Furthermore, the evolution of the microscopic structure can be taken into account, by using this method. More recently, some sophisticated efforts for investigating localization phenomena with computational homogenization tools, appear in the literature [15–20]. It is worth noticing that in the majority of these articles, a continuous damage law has been used to simulate failure in the microscopic scale. In the macroscopic scale a discontinuous law describing a macro crack is numerically obtained.

Some efforts toward coupling contact mechanics and computational homogenization also appear. In [21,22] contact mechanics is coupled with multi-scale homogenization for the study of rough surfaces while in [23] the investigation of a three body frictional system with rigid particles, embedded between a deformable elastic solid and a rigid surface, takes place. Debonding between the matrix and the surrounding inclusions of a composite material has been also studied in [24] by using contact mechanics and homogenization concepts. In [25,26] an interfacial failure model with normal and tangential brittle-elastic springs and a bi-linear

<footnote><sup>a</sup>Institute of Continuum Mechanics, Leibniz University of Hannover, Hannover, Germany</footnote>

<footnote><sup>b</sup>Institute of Computational Mechanics and Optimization, Technical University of Crete, Chania, Greece</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of the multi-scale, computational homogenization.</figcaption></figure>

cohesive zone model have been developed respectively, for the simulation of progressive debonding in the matrix-fiber interface of composite materials. In [27] a bilinear cohesive law is used to simulate failure between the fibers and the matrix interface; then, a comparison between periodic and minimal kinematic boundary conditions is given. In [28] a unilateral contact-friction and a damage evolution law have been used for the simulation of a masonry, brick-mortar interface, within a homogenization framework.

In this article an Augmented Lagrangian formulation has been used for the simulation of the contact state between the constitutive materials, in the microscopic scale of composite structures. A continuous macroscopic model is considered in the macro scale and an overall multi-scale contact computational homogenization scheme is developed. With this numerical scheme, several phenomena related to the microscopic contact conditions and its impact on the macroscopic model, have been investigated. Among them are included the influence of the jump of displacements on the macroscopic response, as well as the alteration of the macroscopic load direction during analysis resulting in a microscopic contact change and its impact on the macroscopic structural behavior.

## 2. A short introduction to computational homogenization

The approach adopted in this article is related to the concurrent analysis of the macroscopic and the microscopic structure, respectively. According to the classical formulation of the multi-scale computational homogenization [6–8,10], two nested boundary value problems are concurrently solved. The initial heterogeneous macroscopic structure is equivalent with a homogeneous one, in each Gauss point of which, a suitably defined RVE is correlated. This RVE includes every heterogeneity and non-linearity of the material.

With linear or periodic boundary conditions, a macroscopic strain is the loading of the RVE. After analysis and convergence of each RVE in every Gauss point, results concerning the average stress and the stiffness are given back to the macroscopic structure, Fig. 1. No assumption for the constitutive law of the macroscopic structure is a priori considered, thus the macroscopic constitutive behavior is numerically obtained. This is a practical solution to the major question of homogenization, namely which are the properties of the homogeneous constitutive law.

## 3. Averaging relations

According to the Hill–Mandel condition or energy averaging theorem, the macroscopic volume average of the variation of work equals to the local work variation, on the RVE [29]:

$$ \boldsymbol{\sigma}^M : \delta\boldsymbol{\epsilon}^M = \frac{1}{V_m} \int_{V_m} \boldsymbol{\sigma}^m : \delta\boldsymbol{\epsilon}^m dV_m \quad (1) $$

Among others, three widely used types of loading states, which satisfy the above condition, can be applied to the RVE: (a) prescribed linear displacements, (b) prescribed tractions, (c) periodic boundary conditions. In the present study both linear displacement and periodic boundary conditions have been used.

According to linear displacement boundary conditions, the loading in the boundaries of the RVE is given by the following relation:

$$ \mathbf{u}|_{\partial V_m} = \boldsymbol{\epsilon}^M \mathbf{x} \quad (2) $$

where a loading strain $\boldsymbol{\epsilon}^M$ is applied to the boundaries $\partial V_m$ of the RVE. With $\mathbf{x}$ is denoted the matrix with the undeformed coordinates of the boundary nodes of the RVE.

Periodic boundary conditions require periodic displacements, as well as antiperiodic tractions, in the opposite boundaries of the RVE. In particular, the displacements of the opposite boundaries are given by the following equations:

$$ \mathbf{u}_T - \mathbf{u}_B = \mathbf{u}_4 - \mathbf{u}_1 \quad (3a) $$

$$ \mathbf{u}_L - \mathbf{u}_R = \mathbf{u}_1 - \mathbf{u}_2 \quad (3b) $$

where the displacements in the top, bottom, left and right boundary are estimated by using the prescribed displacements of three corner nodes of the RVE, namely 1, 2 and 4, given by relation (2).

In order to proceed in the formulation of a homogenization scheme, the average quantities of both the microscopic strain and stress should be defined. The constitutive relation will be then numerically built, as is shown later in this article. The general averaging relations, are:

$$ \langle\boldsymbol{\epsilon}\rangle_{V_m} = \frac{1}{V_m} \int_{V_m} \boldsymbol{\epsilon}^m dV_m, \quad \langle\boldsymbol{\sigma}\rangle_{V_m} = \frac{1}{V_m} \int_{V_m} \boldsymbol{\sigma}^m dV_m \quad (4) $$

Eq. (4) can be further simplified. The volume average microscopic strain is equal to the macroscopic strain which has been applied as loading to the boundaries of the RVE:

$$ \langle\boldsymbol{\epsilon}\rangle_{V_m} = \boldsymbol{\epsilon}^M \quad (5) $$

In case prescribed displacements are applied to the RVE, the following simplified formulation for the macroscopic stresses, has been chosen [10,13]:

$$ \langle\boldsymbol{\sigma}\rangle_{V_m} = \frac{1}{V_m} \mathbf{f}\mathbf{x} = \boldsymbol{\sigma}^M \quad (6) $$

where $\mathbf{f}$ is the matrix of the resulting external forces in the undeformed coordinates of the boundary nodes $\mathbf{x}$ of the RVE, after microscopic analysis has been completed.

A similar relation is chosen for the macroscopic stress, in case periodic boundary conditions are used [10,13]:

$$ \langle\boldsymbol{\sigma}\rangle_{V_m} = \frac{1}{V_m} \mathbf{f}_p \mathbf{x}_p = \boldsymbol{\sigma}^M \quad (7) $$
where $\mathbf{f}_p$ denotes the external forces in the three corner nodes with prescribed displacements and $\mathbf{x}_p$ the undeformed coordinates of these nodes.

For the completion of the homogenization procedure, the stiffness of the macroscopic structure should be calculated. In particular, the consistent tangent stiffness must be incorporated in the Newton-Raphson incremental iterative procedure. This is obtained numerically, by requiring the relation between the variations of the average microscopic (e.g. the macroscopic) stresses and strains. In the next sections will be shown how it is obtained in the framework of computational homogenization, by using a classical Augmented Lagrangian formulation in the microscopic level, for the description of the unilateral contact problem.

# 4. Computational homogenization and contact mechanics

## 4.1. Formulation of the multi-scale scheme

In this section the proposed contact formulation, which is incorporated within the microscopic RVE, is presented. In particular, the classical multi-level computational homogenization approach will be initially given. It will follow the presentation of the numerical scheme adopted for the simulation of the contact between the fibers and the matrix in the RVE.

The formulation given in [10,13] is followed for the implementation of the multi-scale numerical scheme, when both prescribed displacements and periodic boundary conditions are applied to the boundaries of the RVE. For the prescribed displacements, after convergence of the RVE analysis, the system of equilibrium equations is partitioned:

$$ \begin{bmatrix} \mathbf{K}_{ff} & \mathbf{K}_{fs} \\ \mathbf{K}_{sf} & \mathbf{K}_{ss} \end{bmatrix} \begin{bmatrix} \delta \mathbf{u}_f \\ \delta \mathbf{u}_s \end{bmatrix} = \begin{bmatrix} 0 \\ \delta \mathbf{f}_s \end{bmatrix} \qquad (8) $$

With the letter $f$ are denoted the free and with $s$ the supported degrees of freedom. Eq. (8) lead to:

$$ \mathbf{K}^M \delta \mathbf{u}_s = \delta \mathbf{f}_s \qquad (9a) $$

$$ \mathbf{K}^M = \mathbf{K}_{ss} - \mathbf{K}_{sf} \mathbf{K}_{ff}^{-1} \mathbf{K}_{fs} \qquad (9b) $$

In case periodic boundary conditions are applied to the boundaries of the RVE, constraint Eqs. (3) should be taken into account. Lagrange multipliers [13] or transformation of the system of equations [10], can be used respectively. In this work the second approach has been followed, according to the description given in [30]. Thus, the initial system of equations is properly transformed to a new condensed system. Degrees of freedom that correspond to the bottom and left boundaries of the RVE are condensed out.

Just for completeness, are given here the main steps of the procedure. Initially, the constraint equations in the discretized domain are written again and partitioned to retained ($r$) and condensed ($c$) degrees of freedom:

$$ \mathbf{C} \delta \mathbf{u} = \mathbf{0} \Rightarrow \begin{bmatrix} \mathbf{C}_r & \mathbf{C}_c \end{bmatrix} \begin{bmatrix} \delta \mathbf{u}_r \\ \delta \mathbf{u}_c \end{bmatrix} = \mathbf{0} \qquad (10) $$

Eq. (10) leads to:

$$ \delta \mathbf{u}_c = \mathbf{C}_{rc} \delta \mathbf{u}_r \qquad (11a) $$

$$ \mathbf{C}_{rc} = -\mathbf{C}_c^{-1} \mathbf{C}_r \qquad (11b) $$

After transformation of the initial system of equations, the final, condensed system is obtained:

$$ (\mathbf{K}_{rr} + \mathbf{C}_{rc}^T \mathbf{K}_{cr} + \mathbf{K}_{rc} \mathbf{C}_{rc} + \mathbf{C}_{rc}^T \mathbf{K}_{cc} \mathbf{C}_{rc}) \delta \mathbf{u}_r = \delta \mathbf{r}_r + \mathbf{C}_{rc}^T \delta \mathbf{r}_c \qquad (12) $$

This system is solved similar to (8) and (9), by partitioning into prescribed and free degrees of freedom. Now the three corner nodes of the RVE correspond to the prescribed degrees of freedom.

According to the next step of the procedure, the consistent tangent stiffness of the macroscopic Gauss point should be obtained. This is achieved by formulating the relation between the variation of the macroscopic stress and the macroscopic strain. First, the relation of the prescribed displacements either on the whole RVE boundaries (linear displacements) or on the three corner nodes (periodic displacements) is given by:

$$ \delta \mathbf{u}_s = \mathbf{D} \delta \mathbf{\epsilon}^M \qquad (13) $$

where $\mathbf{D}$ is the matrix with the (undeformed) positions on the whole boundary nodes or on the three prescribed corner nodes of the RVE and $\delta \mathbf{\epsilon}^M$ is the variation of the macroscopic strain.

The average stress of the RVE, which is equal to the macro stress, is given by:

$$ \delta \boldsymbol{\sigma}^M = \frac{1}{V_m} \mathbf{D}^T \delta \mathbf{f}_s \qquad (14) $$

where $\delta \mathbf{f}_s$ is the external force vector in the boundaries of the RVE (total boundaries or three corner nodes). By substitution of (9a) and (13) into (14), the stress-strain relation is formulated and the consistent tangent stiffness of the macroscopic level, $\mathbf{K}^M$, is obtained:

$$ \delta \boldsymbol{\sigma}^M = \frac{1}{V_m} \mathbf{D}^T \mathbf{K}^M \delta \mathbf{u}_s = \frac{1}{V_m} \mathbf{D}^T \mathbf{K}^M \mathbf{D} \delta \mathbf{\epsilon}^M \Rightarrow \delta \boldsymbol{\sigma}^M = \frac{1}{V_m} \mathbf{D}^T \mathbf{K}^M \mathbf{D} \delta \mathbf{\epsilon}^M \qquad (15) $$

$$ \mathbf{C}^M = \frac{1}{V_m} \mathbf{D}^T \mathbf{K}^M \mathbf{D} \qquad (16) $$

It is noted that for the periodic boundary conditions, $\mathbf{K}^M$ is given by partitioning of relation (12) into supported and free degrees of freedom.

## 4.2. Definition of the contact interface between constituent materials

For the definition of the behavior of the interface among the constituent materials of the microscale, a unilateral contact law is considered. The simplest, node to node discretization has been chosen to define the contact condition of each interface between the matrix and the fiber of the composite material.

For one nodal pair i depicted in Fig. 2, the unilateral contact law is described by the following relations [31]:

<figure><img src="image_3.png" /><figcaption>Fig. 2. Node to node discretization for the interface among constituents adopted within Matlab.</figcaption></figure>
$$g_{Ni} = (\mathbf{u}_i^2 - \mathbf{u}_i^1) \cdot \mathbf{n}_i^1 + \delta a_i \ge 0 \quad (17a)$$

$$t_{Ni} \le 0 \quad (17b)$$

$$t_{Ni} g_{Ni} = 0 \quad (17c)$$

where $\mathbf{u}_i^1$ is the displacement vector of each node which belongs to a nodal pair of the interface, $\mathbf{n}_i^1$ is the unit normal vector to the master surface and $\delta a_i$ a potential initial gap. Inequality (17a) represents the non-penetration relation and (17b) the requirement that only compressive stresses (contact pressure) are allowed in the interface. Eq. (17c) is the complementarity relation which states that either separation with zero contact stress occurs in the interface or contact is realized with possibly non-zero contact stress.

For the implementation of the above contact formulation, a Matlab code is written. First, the unit normal vectors in each node of the master surface are calculated. Then, relations (17) are written for the whole discretized structure (see Section 4.3). Finally, an Augmented Lagrangian formulation is used for the solution of the contact problem in each interface of the microscale (Section 4.3).

According to the description given here, a simplified approach to evaluate debonding is chosen. Thus, debonding is only depicted by the unilateral contact law of relations (17). No further criteria have been used for debonding initiation/evolution and no tensile resistance has been considered in the interface. This is because the main purpose of this article is to propose a methodology for coupling computational homogenization and contact mechanics. The consideration of these aspects, which would contribute to a more realistic representation of the problem, is left for future work.

### 4.3. Contact computational homogenization

The main idea of this article is related to the investigation of the contact state between the matrix and the fibers, in the microscopic level of a composite material, within a multi-scale framework. To succeed in this, a classical Augmented Lagrangian formulation [31] is properly incorporated in the above mentioned multi-scale homogenization scheme. Thus, the Newton-Raphson incremental, iterative procedure is used for solving the contact problem within an interior loop; an exterior Augmentation is also used for updating the Lagrange multipliers. In the following lines details of the development of the equilibrium of the system, are given. Further details regarding Lagrange multipliers updating, can be found in [31].

An optimization approach has been chosen, for depicting the equilibrium of the system and the contact constraints [32]:

$$\text{minimize} \quad \left( \Pi = \frac{1}{2} \mathbf{u}^T \mathbf{k} \mathbf{u} - \mathbf{f}^T \mathbf{u} \right) \quad (18a)$$

$$\mathbf{g} = \mathbf{M}\mathbf{u} - \mathbf{\delta}\mathbf{a} \geq 0 \quad (18b)$$

$$t_N \le 0 \quad (18c)$$

With $\mathbf{g}$ is denoted the gap function, $\mathbf{\delta}\mathbf{a}$ stands for a potential initial gap while $t_N$ is the vector of the normal forces in the contact interface. The objective function of this scheme (minimization of the potential energy, (18a)) results in the equilibrium of the system, while inequalities (18b) and (18c) correspond to non-penetration and non tensile resistance conditions of the classical unilateral contact problem.

Optimality conditions of the above optimization scheme with the inequalities constraints are written as follows:

$$\nabla \Pi + \frac{1}{2} \epsilon_N \nabla \mathbf{g}^T \mathbf{g} + \nabla \lambda_N^T \mathbf{g} = 0 \quad (19)$$

and after substitutions from (18) the following equilibrium equation arises:

$$\mathbf{F}(\mathbf{u}) = 0 \Rightarrow \mathbf{k}\mathbf{u} - \mathbf{f} + \lambda_N^T \mathbf{M} + \epsilon_N \mathbf{M}^T \mathbf{M}\mathbf{u} - \epsilon_N \mathbf{M}^T \mathbf{\delta}\mathbf{a} = 0 \quad (20)$$

where $\epsilon_N$ is the penalty term and $\lambda_N$ is the Lagrange multiplier vector for inequality constraints. Linearization of the non-linear Eq. (20) with the Newton-Raphson, leads to the linearized Eq. (21):

$$(\mathbf{k} + \epsilon_N \mathbf{M}^T \mathbf{M}) \delta \mathbf{u} = - (\mathbf{k}\mathbf{u} - \mathbf{f} + \lambda_N^T \mathbf{M} + \epsilon_N \mathbf{M}^T \mathbf{M}\mathbf{u} - \epsilon_N \mathbf{M}^T \mathbf{\delta}\mathbf{a}) \quad (21)$$

which is rewritten as:

$$\mathbf{K} \delta \mathbf{u} = -\mathbf{G} \quad (22a)$$

$$\mathbf{K} = \mathbf{k} + \epsilon_N \mathbf{M}^T \mathbf{M} \quad (22b)$$

where vector $\mathbf{G}$ is the residual force vector which incorporates all the contributions from the active contact elements [31] and $\mathbf{k}$ is the initial stiffness matrix of the RVE.

In order to complete the proposed contact multi-scale computational homogenization scheme, coupling between the Augmented Lagrangian formulation and the multi-scale algorithm presented in Section 4.1, should be considered. In particular, after convergence of each RVE analysis, partition is applied to the tangent stiffness matrix $\mathbf{K}$ of the RVE (relation (22b)). In case linear displacement boundary conditions are applied to the RVE, the full system of equations is partitioned (relation (8)), while for periodic boundary conditions the condensed system (12) is partitioned, respectively. After partitioning, $\mathbf{K}^M$ matrix is obtained and the tangent stiffness of the macroscopic level is given by relation (16). Finally, the macroscopic stress is calculated from Eq. (14).

## 5. Some aspects of the proposed approach

### 5.1. Separation of scales

The concept of the multi-scale computational homogenization presented in this work, is based on the separation of scales. Thus, the length scale of the microscopic model should be significantly smaller than the characteristic size of the macroscopic structure. Due to this principle, the macroscopic deformation field can be considered as a uniform loading factor for the microscopic RVE. In addition, according to the description given before, the proposed method adopts a first order homogenization approach. Consequently, size effects are not taken into account. This means that the proposed method cannot yet depict severe deformation gradients in the macroscopic model resulting in localization phenomena. This extension is left for a future work.

### 5.2. Parallel analysis

The computational cost is a very important parameter for the multi-level computational homogenization. This is due to the consideration of a non-linear micro problem, into an overall, macroscopic, non-linear analysis. In order to reduce the computational cost, parallelization of the numerical scheme, has been adopted.

The proposed multi-level computation homogenization model is schematically presented in Table 1. In the shaded area of this Table is shown how parallel processing is applied in the numerical model. The adopted procedure is straightforward, thus each RVE analysis takes place in a separate core, in parallel. This is an advantage of multi-scale methods.

## 6. Details of the numerical model

In the framework of the present study, several models have been developed with the proposed contact computational homogenization scheme. Within these models, the macroscopic domain is a simple continuous, rectangular structure. The RVE, which is called in each macroscopic Gauss point, represents the microscopic
<table><caption>**Table 1**
Proposed numerical scheme.</caption><thead><tr><th colspan="4">Macroscopic analysis</th></tr></thead><tbody><tr><td colspan="4">Assign one RVE in each macroscopic Gauss Point</td></tr><tr><td colspan="4">First initial analysis</td></tr><tr><td colspan="4">ε<sup>M</sup>: initial test strain load for each RVE</td></tr><tr><td colspan="4">Linear displacements/Periodic boundary conditions</td></tr><tr><td colspan="4">RVE analysis</td></tr><tr><td colspan="4">Non-linear contact RVE analysis: Augmented Lagrangian formulation</td></tr><tr><td colspan="4">Output for the macroscopic model: Consistent tangent stiffness, C<sup>M</sup></td></tr><tr><td colspan="4">Return in the macroscopic model</td></tr><tr><td colspan="4">Application of macroscopic loading, BCs</td></tr><tr><td colspan="4">Apply macroscopic increments</td></tr><tr><td colspan="4">Apply macroscopic iterations</td></tr><tr><td colspan="4">Assembly the tangent stiffness matrix</td></tr><tr><td colspan="4">Solution of the macroscopic system</td></tr><tr><td colspan="4">Loop in every integration point of the macro model</td></tr><tr><td colspan="4">Parallel Processing</td></tr><tr><td></td><td>Macro - Gauss point 1</td><td>Macro - Gauss point 2</td><td>Macro - Gauss point 3</td></tr><tr><td colspan="4">Loading strain for the RVE:</td></tr><tr><td></td><td>ε<sup>M</sup><sub>1</sub></td><td>ε<sup>M</sup><sub>2</sub></td><td>ε<sup>M</sup><sub>3</sub></td></tr><tr><td colspan="4">RVE analysis</td></tr><tr><td colspan="4">Output for the macroscopic model: Consistent tangent stiffness - stress</td></tr><tr><td></td><td>C<sup>M</sup><sub>1, 1</sub></td><td>C<sup>M</sup><sub>2, 1</sub></td><td>C<sup>M</sup><sub>3, 1</sub></td></tr><tr><td colspan="4">Return in the macroscopic model</td></tr><tr><td colspan="4">Calculation of the internal force vector</td></tr><tr><td colspan="4">Convergence control</td></tr><tr><td colspan="4">If convergence is achived --&gt; Next increment</td></tr><tr><td colspan="4">If convergence is not achived --&gt; Next iteration</td></tr></tbody></table>

structure of a composite material; thus, it consists of a matrix and a number of fibers.

For the discretization of the contact between the matrix and the surrounding fibers, the simplest node to node case has been chosen in a small displacements framework, while friction is omitted. In addition, the tensile strength of the interface is regarded equal to zero and no other plasticity constitutive law has been considered

in the proposed scheme. This would probably give more realistic results. However, the main goal of this article is the examination of the contact behavior, consequently, only a pure contact problem is analyzed in the microstructure.

For both the macro, and the micro structure, plane stress, full integration finite elements have been used. The Newton-Raphson incremental iterative procedure has been also used for the solution of the microscopic problem. This problem is non-linear, due to the either-or, contact-no contact condition between the matrix and the fibers. For the numerical implementation of the whole described scheme, a Matlab code has been developed.

In every RVE model, the same fraction between the matrix and the fibers is considered, thus 62% matrix and 38% fibers. The material properties adopted in the RVE, are chosen from the literature [24] and presented here:

Matrix: E<sub>m</sub> = 68.89 Gpa, n<sub>m</sub> = 0.33.

Fibers: E<sub>f</sub> = 413.04 Gpa, n<sub>f</sub> = 0.20.

It is worth noticing, that according to the proposed research no material law is initially considered in the macroscopic scale of the problem. The material law and the relevant stress-strain relation, is built only numerically, in the course of the computational homogenization procedure.

The geometry of the macroscopic structure, the loading and the boundary conditions which are initially considered, are shown in Fig. 3(a). A second, tensile load has been also applied to the macroscopic structure, according to Fig. 3(b). In Fig. 4 the mesh of the microscopic models consisting of one, nine and twenty inclusions, is presented.

## 7. Initial results

For the investigation of the influence of the RVE mesh density, two types of analyses were conducted for the models with one, four and nine micro fibers and the macroscopic load shown in Fig. 3(a): one with denser mesh and another with less dense mesh. The denser mesh of the models with one and nine fibers is shown in Fig. 4. The mesh size of the macroscopic structure is considered constant, and equal to 25 elements, Fig. 3. According to Fig. 5, the force-displacement diagrams of the macro model, obtained for each case, are quite close indicating that within the proposed approach, there is a small sensitivity of the RVE mesh density on the results.

Another short investigation related to the influence of the macroscopic mesh on the results, was conducted. Three macro models with 25, 100 and 225 finite elements were tested, for the

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Geometry of the macro model with a mesh of 25 elements and two load cases.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. RVE mesh for increasing fibers (a) 1 fiber – mesh of 665 nodes (1330 d.o.f.), (b) 9 fibers – 2346 nodes (4692 d.o.f.), (c) 20 fibers – 3592 nodes (7184 d.o.f.).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Force-displacement diagrams of the macro structure, after consideration of a dense and a less dense mesh size for the RVE.</figcaption></figure>

microscopic model with one fiber. A concentrated macroscopic
vertical loading (instead of a distributed one) is applied to the
upper-right corner node of the macroscopic structure. The
corresponding force-displacement diagrams presented in Fig. 6
demonstrate a divergence in the response, indicating that the

macroscopic mesh size can influence the results. However as the
macroscopic mesh size is increased, the difference in the response
is reduced. To reduce computational time, the macroscopic model
with 25 elements has been considered in the following simulations
of this work.

To study the effect of the boundary conditions which are
applied to the RVE, on the macroscopic behavior of the structure,
both linear and periodic boundary conditions have been tested.
According to Fig. 7 only a small difference in the corresponding
force-displacement diagrams, is obtained.

Finally, an overview of the behavior of the composite material is
presented in the following lines. In Fig. 8 the macro structure and
the respective RVEs in two selected Gauss points are given, for the
model with four fibers, periodic boundary conditions in the RVE
and the macroscopic load shown in Fig. 3(a). According to this
Figure, debonding appears in the upper left part of the macro
model, however, it does not appear in the bottom left part of it.
Later in this article, more detailed results about debonding will
be presented.

The proposed multi-level computational homogenization
model can be used for bigger and more realistic microscopic
formulations. Thus, the same macroscopic structure and a

<figure><img src="image_7.png" /><figcaption>Fig. 6. Force-displacement diagrams of the macro structure, for increasing macroscopic mesh density.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Force-displacement diagrams of the macro structure for linear displacement and periodic boundary conditions.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Vertical displacements of the macroscopic model (m) and RVE with 4 fibers (a) 1st, (b) 7th and (c) 10th time step of the macro load.</figcaption></figure>

microscopic model with twenty fibers have been also considered, Fig. 9. In this microscopic formulation, the inclusions are irregularly distributed into the matrix of the composite material, contrary to the previous cases (see Fig. 4).

To complete the presentation of the initial results, the tensile macroscopic load shown in Fig. 3(b) is examined. Fig. 10 shows the macro structure and the respective RVEs in two selected Gauss points, for the model with four fibers and periodic boundary conditions in the RVE.

## 8. Non-linear behavior due to changing of the macro load direction

If the direction of the macroscopic loading is changed during the finite element analysis, the related force-displacement diagrams are expected to be non-linear, due to the alteration of the microscopic contact state, from contact to no contact, or the opposite. The proposed numerical scheme can depict this non-linear

behavior on the macroscopic model, due to the presence of the uni-lateral contact effects in the micro level.

In order to demonstrate how the non-linear nature of the microscopic contact problem affects the macroscopic structure, a different macroscopic loading has been tested, Fig. 11. The load shown in Fig. 11(a), which is applied in a first load step, imposes a total compression to the macro model and does not cause debonding in the RVE. In a second loading step, the load shown in Fig. 11(b) tends to change the contact status in the RVE and causes debonding in the matrix-fiber interface.

From the results of this investigation some interesting conclusions can be drawn. First, the corresponding force-displacement diagram shown in Fig. 12 is no more a straight line. Instead, it has an inclination due to the changing of the contact state in the micro structure, from no-opening to opening. In addition, the diagram indicates an improved response, in comparison with the case of a single load (where only the load shown in Fig. 11(b) exists). This is attributed to the fact that the stiffness of the structure has been increased, due to this compressive force which is initially
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** RVE with 20 fibers corresponding to macro element 1 in (a) 4th, (b) 7th, (c) 10th time step of the macro load and to macro element 21 in (d) 4th, (e) 7th, (f) 10th time step of the macro load.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Horizontal displacements for the macroscopic model (m) with tensile load and RVE with 4 fibers (a) 1st and (b) 10th time step of the macro load.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** An alternative combined loading on the macroscopic structure in two steps. (a) Loading that does not cause debonding and (b) loading that causes debonding.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Force-displacement diagrams of the macroscopic structure for a combined (in two steps) and a single loading.</figcaption></figure>

applied to the macroscopic model (first load of Fig. 11(a), then load of Fig. 11(b)).

Furthermore, this simple example demonstrates how the macroscopic loading conditions may potentially influence the behavior of the structure when debonding in the microstructure occurs, as well as how debonding itself can influence the macroscopic behavior.

## 9. Consideration of the displacement jump

In [2] is given that the average microscopic strain of the RVE is equal with the loading strain of relation (2) only when the material of the RVE is perfectly bonded. Thus, when debonding occurs the displacement jump should be taken into account. The proposed contact multi-scale homogenization algorithm considers the jump of displacements, which arises between the matrix and the surrounding fibers, as is described by the following relation:

$$ \epsilon^M = \frac{1}{V_m} \int_{V_m} \epsilon^m dV_m + \frac{1}{2V_m} \int_{\partial V_1 \cap \partial V_2} (\|\mathbf{u}\| \otimes \mathbf{n} + \mathbf{n} \otimes \|\mathbf{u}\|) dA \quad (23) $$

where $\epsilon^M$ and $V_m$ are the loading strain and the volume of the RVE, $\partial V_1$, $\partial V_2$ the common boundaries between the matrix and the fibers and $\|\mathbf{u}\|$ the displacement jump which arises after RVE simulation. According to the adopted procedure, the macroscopic strain $\epsilon^M$, which is applied to the RVE, incorporates both the (element by element) integration of strains and the jump (relation (23)).

An example of a single RVE with one fiber is presented here, to show this detail. First, one strain loading that causes debonding, is applied to the RVE. When the RVE simulation is completed, the average strain of the material is calculated, by using the strain integration, in the right part of Eq. (23). With this element by element integration, the average strain except the areas with the gaps, is estimated. Then, the contribution of the displacement jump is calculated by using the second part of the right term of relation (23). The summation of the integration of strains and the jump contribution, results in the strain loading that was initially applied to the RVE.

In a second example, a strain loading which causes an almost zero debonding is applied to the RVE and the above mentioned procedure is repeated. It is confirmed that the loading strain is approximately equal to the element by element integration of strains of the deformed RVE, as the jump contribution is close to zero.

In particular, for a loading strain $\epsilon^M = [0.03 \ 0.06 \ 0.09]^T$ a displacement jump arises after RVE simulation, Fig. 13(a). The element by element integration of strains after RVE simulation, is estimated equal to $[0.006437 \ 0.020688 \ 0.049126]^T$ and the jump of displacements equal to $[0.023526 \ 0.039193 \ 0.040718]^T$. By adding these two vectors, is obtained the vector $[0.029963 \ 0.059881 \ 0.089844]^T$, which is almost equal to the initial loading strain $\epsilon^M$.

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Single RVE analysis (a) with displacement jump and (b) with almost zero displacement jump after the RVE simulation.</figcaption></figure>

For a loading strain $\epsilon^M = [-0.03 \ -0.06 \ 0.09]^T$ the displacement jump after the RVE simulation is very small, close to zero, Fig. 13(b). The element by element integration of strains after RVE simulation, is estimated equal to $[-0.030865 \ -0.060155 \ 0.085260]^T$, thus very close to the loading strain. A quite small jump of displacements equal to $[0.000906 \ 0.000235 \ 0.004859]^T$ is also obtained. By adding these two vectors, it is received the vector $[-0.029959 \ -0.059920 \ 0.090119]^T$, which is almost equal to the initial loading strain $\epsilon^M$, as well.

## 10. Investigation of the debonding in the RVE

As a continuation of the previous section, the influence of the displacement jump resulting in debonding between the constitutive materials, is examined here. A new overall multi-scale model has been developed, where debonding is not permitted between the matrix and the fibers, meaning that a tie connection has been applied between them. It is noted that this is a linear model since the contact-no contact condition has been removed from the micro formulation. Instead, only a fixed contact (tie) restriction has been used.

According to the results, the behavior of the macro structure is significantly improved, in comparison with the case that opening between the constitutive materials is permitted. In Fig. 14 the macroscopic force-displacement diagrams obtained from the two models, from an RVE with one fiber and the macroscopic load shown in Fig. 3(a), are presented. It is shown that debonding causes a reduction in the stiffness of the material resulting in a more flexible behavior, contrary to the model with a tie connection between the matrix and the fibers (e.g. no debonding), in the RVE. Similar diagrams are obtained for the models with more microscopic fibers.

One of the main advantages of the proposed numerical scheme is that it can offer an insight of the debonding in each RVE, corresponding to each integration point of the macroscopic model. Consequently, a measurement of debonding, which can be considered as

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Force-displacement diagrams in case debonding is permitted, or is not permitted in the RVE.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** External macroscopic force-debonding area (sum in the four integration points of each macroscopic element) for the model with four micro fibers.</figcaption></figure>

a debonding map in the macroscopic model, can be given. This is a reflection of the microscopic behavior, on the macroscopic structure.

In particular, the evaluation of debonding throughout the multi-scale analysis, can offer more details about the distribution of it in the macroscopic domain, as well as the areas of the macro structure where this phenomenon becomes more severe. In diagrams shown in Fig. 15 the summation of the debonding corresponding to the four integration points of some macroscopic elements is given for the model with four micro fibers and the external loading

shown in Fig. 3(a). According to this Figure, element 21 of the macroscopic structure suffers the biggest debonding, while element 22 and element 16 follow (the enumeration of the macroscopic elements is given in Fig. 3).

Based on the previously mentioned results, the distribution of debonding in the macroscopic structure is given in Fig. 16. These values appear in each RVE, corresponding to each macroscopic integration point. The map of debonding on the macro structure can be used for the evaluation of this phenomenon, during the incremental-iterative procedure of the multi-level analysis. Fig. 17 shows the distribution of debonding for the second (tensile) load given in Fig. 3(b).

Furthermore, for the alternative macroscopic loading of Fig. 11, the corresponding force-debonding diagrams are shown in Fig. 18. According to this Figure, the diagrams are not linear. The change of the contact status in the microscopic level from contact to no-contact, which is also depicted in the non-linear diagram of Fig. 12, results in a non-linear response: initially no debonding appears; as the macroscopic external load is gradually increased, debonding appears in some RVEs and the corresponding diagrams become non-linear.

In Fig. 19 the distribution of debonding is given, for this alternative loading of the macroscopic structure. Contrary to Figs. 16 and 17, is shown that the distribution of debonding changes, during analysis. Small values close to zero initially appear, due to the total compression of the macroscopic structure, Fig. 11(a). As the

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Distribution of debonding (mm²) in the macroscopic structure for (a) the 1st, (b) the 5th and (c) the 10th time step (shear macroscopic load, Fig. 3(a)).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Distribution of debonding (mm²) in the macroscopic structure for (a) the 1st, (b) the 5th and (c) the 10th time step (tensile macroscopic load, Fig. 3(b)).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** External macroscopic force-debonding area (sum in the four integration points of each macroscopic element) for the model with four micro fibers and the alternative loading.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Distribution of debonding (mm²) in the macroscopic structure for the combined loading in (a) the 1st, (b) the 3rd and (c) the 10th time step.</figcaption></figure>

loading which causes debonding is gradually applied (Fig. 11(b)), debonding is increased and debonding distribution in the macro geometry is changed. The final distribution is similar to the one obtained from the model with the non-combined loading, Figs. 16(c) and 19(c) respectively.

It is worth noticing that the above mentioned non-linear behavior is not attributed to a specific assumption, but is gradually obtained as analysis of the macroscopic model and the detailed RVE with unilateral contact interfaces, is evolved.

## 11. Comparison between the proposed model and direct numerical simulation analysis

For the validation of the proposed method and the results presented in the previous sections, a direct heterogeneous macroscopic finite element model is developed with the commercial software Abaqus. The same approach has been adopted in several scientific publications [10,16,17] for creating a reference solution. Comparison with the results obtained from the multi-scale computational homogenization method is considered (a) for the initial macroscopic shear loading, (b) for the tensile loading and (c) for the combined load case with compressive and shear loading resulting in a non-linear force-displacement diagram. Since the main goal of the article is to present the methodology for developing a numerical scheme which incorporates classical contact mechanics into multi-scale homogenization and not to solve a specific material science problem, the mentioned "numerical experiment" is considered as an adequate validation of the proposed approach, thus no further experimental results have been used.

The multi-scale model which is considered for the comparison, is the model with the matrix and one fiber in each integration point of

the macroscopic structure, according to Fig. 4(a). It is compared with a direct numerical simulation (DNS) model shown in Fig. 20. Unilateral contact law in each interface between the matrix and the fiber, together with plane stress, full integration elements have been used in the DNS analysis. The dimensions and the boundary conditions of the DNS model are the same with the homogeneous macroscopic model considered in the multi-scale analysis, given in Fig. 3.

Before presenting the results of the comparison, some details regarding the DNS analysis in Abaqus as well as some conceptual similarities - differences between the approach proposed in this

<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Mesh of the direct numerical simulation model.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Force-displacement diagrams obtained from the proposed multi-scale method and the DNS analysis for (a) shear, (b) tensile, and (c) combined compressive and shear macroscopic loading.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Multi-scale model and shear loading: Vertical displacements of the macroscopic model (m) and RVE in three Gauss points in the last time step.</figcaption></figure>

article and the DNS analysis are given. A surface to surface discretization is used to simulate contact between the matrix and each of the fibers in the DNS analysis. Concerning the requirements in the normal and the tangential direction of the interfaces, unilateral contact stands for the description of the opening/contact condition in the normal direction. Zero tensile resistance is considered in the matrix-fiber interface, thus no traction – separation

law is used. In the tangential direction some soft springs were used, to avoid rotation of the fibers.

On the proposed contact multi-scale homogenization scheme a node to node discretization was used. To eliminate rotation of the fibers, some soft springs were used. Similar to DNS analysis, no traction – separation law has been used, indicating zero tensile resistance in each interface.
<figure><img src="image_24.png" /><figcaption>**Fig. 23.** DNS Analysis and shear loading: Vertical displacements of the macroscopic model (m) in the last time step (a scale factor two has been used).</figcaption></figure>

One important difference between the two methods is that in the framework of the multi-scale approach, a significantly smaller number of contact interfaces than the DNS model, appears. For instance, for the model with the matrix and one fiber in the RVE, one unilateral contact interface appears. On the contrary, 100 (10 × 10) interfaces are introduced in the respective DNS model (see also Figs. 4(a) and 20). For the case of 20 fibers in the RVE, a DNS model with 40000 (200 × 200) interfaces would be needed for the simulation in Abaqus. Due to the large number of the contact interfaces, numerical singularities would possibly arise, indicating difficulties to obtain a solution to this problem.

<figure><img src="image_25.png" /><figcaption>**Fig. 25.** DNS Analysis and tensile loading: Horizontal displacements of the macroscopic model (m) in the last time step (a scale factor five has been used).</figcaption></figure>

For the three mentioned macroscopic load cases, the respective force-displacement diagrams given in Fig. 21 show a good comparison between the DNS model and the proposed approach.

For the shear macroscopic loading, the state of the contact condition between the matrix and the fiber is similar in both models, according to Figs. 22 and 23. No opening occurs in the bottom and left part of the macroscopic model while opening appears in the top and left part of it. A smaller opening appears in the center of structure.

When a tensile macroscopic loading is applied to the structure, some debonding appears in the whole macroscopic domain. A similar image is obtained from both the multi-scale model and the DNS analysis, according to Figs. 24 and 25.

<figure><img src="image_26.png" /><figcaption>**Fig. 24.** Multi-scale model and tensile loading: Horizontal displacements of the macroscopic model (m) and RVE in three Gauss points in the last time step.</figcaption></figure>
## 12. Conclusions

In the present article, a multi-scale computational homogenization approach is proposed for the study of composite materials. A classical unilateral contact model has been incorporated in the overall numerical scheme, for the study of the microscopic structure. Thus, an Augmented Lagrangian formulation is used for the evaluation of the contact state between the matrix and the surrounding fibers, in each RVE. Several microscopic models are then tested.

The presence of the contact law in the RVE, affects the behavior of the macroscopic structure. Depending on the macro load, the contact state in the RVE can be changed from opening to no-opening and vice versa, resulting in a non-linear macroscopic force-displacement diagram. In addition, is shown that a proper macroscopic load combination can result in the increase of the response, because the contact status in the micro scale is properly affected by the macro load.

The influence of debonding on the structural response of the macro model has been also presented. In case debonding is not permitted in the microscopic level, the behavior of the macroscopic structure is significantly improved. This attributed to the fact that when debonding appears in the RVE, a loss of stiffness is caused in the macro level.

In the framework of the proposed approach, distribution of debonding can be represented in the macroscopic model. Contour diagrams are used, as a tool for the evaluation of debonding in the macroscopic level. Relevant diagrams showing debonding in each RVE when the external macroscopic load is gradually increased, are presented. For a different macroscopic load combination, these diagrams become non-linear: first a branch with small or zero debonding appears, then another branch with increasing debonding follows.

In the proposed method a first order homogenization approach has been followed. As a result it cannot be used for depicting localization phenomena in the macroscopic scale. The extension of the method on this type of problems as well as the incorporation of a contact constitutive formulation representing tensile resistance in the RVE interfaces, will be the subject of future investigation.

## Acknowledgements

The research project is implemented within the framework of the Action Supporting Postdoctoral Researchers of the Operational Program "Education and Lifelong Learning" (Actions Beneficiary: General Secretariat for Research and Technology), and is co-financed by the European Social Fund (ESF) and the Greek State.

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

Co-financed by Greece and the European Union

<figure><img src="image_29.png" /></figure>

## References

[1] E. Sanchez-Palencia, Non-Homogeneous Media and Vibration Theory, Lecture Notes in Physics, Springer, Berlin, 1980.

[2] T.I. Zohdi, P. Wriggers, An Introduction to Computational Micromechanics, Springer, The Netherlands, 2008.

[3] V.P. Nguyen, M. Stroeven, L.J. Sluys, J. Multisc. Model. 3 (2011) 1–42.

[4] H.E. Pettemann, S. Suresh, Int. J. Solids Struct. 37 (2000) 5447–5464.

[5] C. Dascalu, G. Bilbie, E.K. Agiasofitou, Int. J. Solids Struct. 45 (2008) 409–430.

[6] P.M. Suquet, Local and Global Aspects in the Mathematical Theory of Plasticity, Plasticity Today: Modelling, Methods and Applications, Elsevier, London, 1985.

[7] R. Smit, W. Brekelmans, H. Meijer, Comput. Methods Appl. Mech. Eng. 155 (1998) 181–192.

[8] F. Feyel, Comput. Mater. Sci. 16 (1999) 344–354.

[9] C. Miehe, J. Schröder, J. Schotte, Comput. Methods Appl. Mech. Eng. 171 (1999) 387–418.

[10] V. Kouznetsova, Computational Homogenization for the Multi-Scale Analysis of Multi-Phase Materials, PhD Thesis, Technical University Eindhoven, The Netherlands, 2002.

[11] G.A. Drosopoulos, P. Wriggers, G.E. Stavroulakis, Contact analysis in multi-scale computational homogenization, in: M. Jirásek, O. Allix, N. Moës, J. Oliver (Eds.), Proceedings of 3rd International Conference on Computational Modeling of Fracture and Failure of Materials and Structures (CFRAC), Prague, Czech Republic, 2013, p. 226.

[12] G.A. Drosopoulos, P. Wriggers, G.E. Stavroulakis, Incorporation of contact mechanics in multi-level computational homogenization for the study of composite materials, in: Proceedings of 3rd International Conference on Computational Contact Mechanics (ICCCM), Lecce, Italy, 2013.

[13] C. Miehe, A. Koch, Arch. Appl. Mech. 72 (2002) 300–317.

[14] T.J. Massart, R.H.J. Peerlings, M.G.D. Geers, Int. J. Damage Mech. 16 (2007) 199–226.

[15] T.J. Massart, R.H.J. Peerlings, M.G.D. Geers, Int. J. Numer. Meth. Eng. 69 (2007) 1022–1059.

[16] V.P. Nguyen, M. Stroeven, L.J. Sluys, Eng. Fract. Mech. 79 (2012) 78–102.

[17] T. Belytschko, S. Loehnert, J.H. Song, Int. J. Numer. Meth. Eng. 73 (2008) 869–894.

[18] C.V. Verhoosel, J.J.C. Remmers, M.A. Gutiérrez, R. de Borst, Int. J. Numer. Meth. Eng. 83 (2010) 1155–1179.

[19] E.W.C. Coenen, V.G. Kouznetsova, E. Bosco, M.G.D. Geers, Int. J. Fract. 178 (2012) 157–178.

[20] E.W.C. Coenen, V.G. Kouznetsova, M.G.D. Geers, Int. J. Numer. Meth. Eng. 90 (2012) 1–21.

[21] P. Wriggers, J. Nettingsmeier, Computational Contact Mechanics, CISM Courses and Lectures, vol. 498, SpringerWienNewYork, Italy, pp. 129–161.

[22] P. Wriggers, J. Reinelt, Comput. Methods Appl. Mech. Eng. 198 (2009) 1996–2008.

[23] I. Temizer, P. Wriggers, Comput. Methods Appl. Mech. Eng. 198 (2008) 377–396.

[24] P. Wriggers, G. Zavarise, T.I. Zohdi, Comput. Mater. Sci. 12 (1998) 39–56.

[25] A. Caporale, R. Luciano, E. Sacco, Comput. Struct. 84 (2006) 2200–2211.

[26] V. Kushch, S. Shmegera, P. Brøndsted, L. Mishnaevsky Jr., Int. J. Eng. Sci. 49 (2011) 17–29.

[27] H.M. Inglis, P.H. Geubelle, K. Matouš, Phil. Mag. 88 (2008) 2373–2397.

[28] E. Sacco, F. Lebon, Int. J. Solids Struct. 49 (2012) 3666–3680.

[29] R. Hill, J. Mech. Phys. Solids 11 (1963) 357–372.

[30] R.D. Cook, D.S. Malkus, M.E. Plesha, R.J. Witt, Concepts and Applications of Finite Element Analysis, forth ed., John Wiley and Sons, Inc, United States, 2002.

[31] P. Wriggers, Computational Contact Mechanics, Springer, The Netherlands, 2006.

[32] S.A. Georgiadis, G.E. Stavroulakis, E.K. Koltsakis, P.D. Panagiotopoulos, Comput. Struct. 34 (1990) 735–752.