# A new fully coupled two-scales modelling for mechanical problems involving microstructure: The 95/5 technique

Djamel Missoum-Benziane, David Ryckelynck, Francisco Chinesta

## ► To cite this version:

Djamel Missoum-Benziane, David Ryckelynck, Francisco Chinesta. A new fully coupled two-scales modelling for mechanical problems involving microstructure: The 95/5 technique. Computer Methods in Applied Mechanics and Engineering, 2007, 196 (21-24), pp.2325 - 2337. 10.1016/j.cma.2006.10.013. hal-01004911

HAL Id: hal-01004911

https://hal.science/hal-01004911

Submitted on 2 Apr 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

Public Domain
# A new fully coupled two-scales modelling for mechanical problems involving microstructure: The 95/5 technique

D. Missoum-Benziane, D. Ryckelynck, F. Chinesta

LMSP, UMR 8106, CNRS-ENSAM-ESEM, 151 boulevard de l'Hôpital, F-75013 Paris

In this work, we propose a new fully coupled two-scales modelling for mechanical problems involving microstructure which makes use of a meshless constrained natural element approximation for accounting the large and very localized variations in the nodal density. This approach is an appealing choice for treating problems involving a microstructure (requiring fine local descriptions) whose evolution can be accurately described in the macroscopic scale using coarse approximations. This approach leads to a model where more than 95% of the nodes are in less than 5% of the considered domain volume.

The technique proposed in this paper allows considering simultaneously both scales, the one related to the microstructure description and the one related to its evolution (the macroscopic one). Moreover, in contrast to the vast majority of homogenization techniques, the presented approach allows an accurate description of the boundary conditions, because the microscopic domains can be located on the domain boundary.

*Keywords:* Multi-scale problems; Computational homogenization; Model reduction; Constrained natural element method; Numerical modelling

## 1. Introduction

Nowadays, experimental methods allow to describe many thermo-mechanical transformations at the microscopic scale. The main challenge is how to transfer these microscopic descriptions to the macroscopic scales related to structures or processes. Thus, for example, in the framework of the cutting process simulation, we are specially interested on the tool wear prediction. To perform this kind of prediction we need macroscopic models involving microscopic description of the thermo-mechanical transformation near the contact area between the tool and the work piece. Thus, micro-macro models allowing an accurate description of microscopic boundary conditions for both the tool and the work-piece are required.

Classical computational homogenization techniques fail in the description of the microscopic boundary conditions. To circumvent this difficulty, a full microscopic modelling of the whole mechanical system can eventually be considered, making use of parallel computing as well as domain decomposition [9,11,7]. One possibility for avoiding a full microscopic modelling lies in the use of the Arlequin approach [2]. This method consists of superposing microscopic and macroscopic descriptions in different regions of the mechanical system, that can include the domain boundary. The compatibility between the microscopic and the macroscopic displacement fields is imposed in a weak form. Other possibility lies in the use of the finite element superposition method (S-FEM) [17] which consider a displacement superposition inside each representative volume element and enforces the displacement continuity on the microscopic domain boundary.

In this work, we propose to restrict the fine microscopic descriptions to some regions in the system, whose size is very reduced in comparison to the whole domain, but it
<figure><img src="image_1.png" /><figcaption>Fig. 1. Nodal distribution in the microscopic patterns as well as in the complementary domain.</figcaption></figure>

is large enough to represent the local microstructure. In contrast to the vast majority of multilevel approaches, the one proposed in this work uses a single description of the displacement field. The cloud of nodes has higher density on the microscopic patterns $\Omega_i$, which are surrounded by a complementary domain $\Xi$ with a low density of nodes $((\cup \Omega_i) \cup \Xi = \Omega)$. This situation is depicted in Fig. 1.

In these kind of models about 95% of the nodes are located in about the 5% of the domain volume, which justifies the 95/5 appellation. The very high and localized gradient of the nodal density in the microscopic domain boundaries neighborhood is accurately accounted thanks to the meshless character of the natural element approximation considered for the interpolation of the different fields [16]. Microscopic domains are coupled to the complementary region through their boundaries. In order to insure the appropriate transmission between both kind of domains we need to define an interpolation inside each microscopic domain that is coupled to the complementary region through the nodes located on its boundary. In the same way the interpolation in the complementary domain is coupled with the microscopic patterns only through the nodes located on their boundaries. Thus, one needs define the interpolation in the non-convex domain resulting of removing the microscopic zones to the whole domain which must be compatible with the approximation defined inside each microscopic domain, as required in our approach. In the context of the natural element method two possibilities exist for accounting non-convex domains, the first one the so-called alpha-NEM proposed by Cueto et al. in [5] and the second one, the C-NEM proposed and used in [20–22]. In this work the last approach is considered because its boundary description simplicity. As described later, the natural element interpolation is based on the use of Voronoi cells instead of Delaunay triangles (used in the finite element framework) for defining the

fields approximation. It has been proved in some of our former works [12] that the quality of the Delaunay triangulation (dual of the Voronoi diagram) does not affect the quality of the computed solution in the NEM framework, even when one proceeds in very distorted meshes.

The main contribution of the present work is the extension of the elastic constitutive law computed in the microscopic regions to the complementary domain. For this purpose the homogenized behavior computed at each microscopic region spreads to its immediate neighborhood, and it is extended to the complementary domain by using a natural element interpolation. The complementary medium has a particular role in the multi-scale model, because it insures the connection between all the microscopic patterns. Moreover, as proved later, we can proceed by introducing a domain decomposition approach, where only a reduced number of modes defined in each microscopic region is considered.

In the next sections, we first provide the details about the interpolation of the behavior law computed in the microscopic patterns towards the macroscopic scale. Then, we validate the use of C-NEM (constrained natural element method) for highly heterogeneous nodal densities. The main ideas concerning the model reduction and the domain decomposition will be then addressed. Finally, some numerical examples will be considered for validating the proposed approach.

## 2. Extending microscopic information towards the macroscopic scale

The technique proposed in this work been very recent, only problems involving linear elasticity are here considered. Let's $\Omega$ be the domain of interest. A microscopic model is built on different zones $(\Omega_i)_{i=1...m}$. Each $\Omega_i$ is a microscopic pattern. These patterns are either surrounded
<figure><img src="image_2.png" /><figcaption>Fig. 2. Microscopic pattern and complementary domain.</figcaption></figure>

by a complementary medium $\Xi$ (see Fig. 2) or located on the domain boundary.

The constitutive law is assumed to be perfectly defined at the microscopic scale inside each one of the microscopic patterns. In the complementary domain we assume that the spatial variations of the averaged (or macroscopic) stress and the strain fields are smooth enough, as well as the evolution of the homogenized tensor (i.e. the microstructure is uniformly distributed in the domain or it evolves continuous and smoothly, and the number and location of the microscopic patterns is able to represent accurately the evolution of the homogenized behavior). In the case of arbitrary microstructure heterogeneities randomly distributed this technique cannot be applied, or it could be applied in each region where the homogenized behavior evolves smoothly. Since the technique that we propose is not in fact a multi-level method we cannot differentiate the displacement, stress and strain fields on the microscopic pattern and on the complementary domain. Thus, we only differentiate the constitutive law in the microscopic scale (where it is perfectly defined) and the macroscopic one (obtained by interpolating the ones homogenized in the microscopic domains).

Let's $\sigma$ and $\epsilon$ be the vector representation of the stress and strain, respectively, that is:

$$ \epsilon = \begin{pmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{12} \end{pmatrix} = \begin{pmatrix} \frac{\partial u_1}{\partial x} \\ \frac{\partial u_2}{\partial y} \\ \frac{1}{2} \left( \frac{\partial u_1}{\partial y} + \frac{\partial u_2}{\partial x} \right) \end{pmatrix} \quad (1) $$

and

$$ \sigma = \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{pmatrix} \quad (2) $$

At each point $P \in \Omega_i$ the elastic behavior is written as

$$ \sigma_{|P} = C(P)\epsilon_{|P} \quad \forall P \in \Omega_i \quad (3) $$

where $C$ is the local microscopic elasticity tensor (a $3 \times 3$ matrix in the 2D case). This tensor is perfectly defined in $\Omega$, despite the fact that we only consider its expression inside each microscopic domain.

Now, we consider the microscopic domain $\Omega_i$, whose boundary is denoted by $\partial\Omega_i$, that is subjected to the three boundary displacements depicted in Fig. 3. It is easy to prove that when we apply the displacement field $U_i^T(x \in \partial\Omega_i) = (x, 0)$ on its boundary, the strain average

<figure><img src="image_3.png" /><figcaption>Fig. 3. Homogenous deformation condition: (left) $\langle \epsilon_{xx} \rangle = 1$; (center) $\langle \epsilon_{yy} \rangle = 1$; (right) $\langle \epsilon_{xy} \rangle = 1$.</figcaption></figure>

in that microscopic domain results $\langle \epsilon_1 \rangle^T = (1, 0, 0)$. In the case of prescribing $U_{II}^T(x \in \partial\Omega_i) = (0, y)$ it results $\langle \epsilon_{II} \rangle^T = (0, 1, 0)$, and finally with $U_{III}^T(x \in \partial\Omega_i) = (y, x)$ the strain average is given by $\langle \epsilon_{III} \rangle^T = (0, 0, 1)$. The classical homogenization techniques assume that $\langle \sigma \rangle = \hat{C}_i \langle \epsilon \rangle$. It can be noticed that the average stresses $\langle \sigma_1 \rangle$, $\langle \sigma_2 \rangle$ and $\langle \sigma_3 \rangle$ related to the solution of the elastic problems defined by the three boundary conditions $U_i^T$, $U_{II}^T$ and $U_{III}^T$, respectively, define the columns of matrix $\hat{C}_i$.

Let's $\langle \sigma \rangle_i$ and $\langle \epsilon \rangle_i$ be the mean values of the stress and the strain, respectively, computed in $\Omega_i$ related to a real mechanical problem:

$$ \langle \sigma \rangle_i = \frac{1}{|\Omega_i|} \int_{\Omega_i} \sigma_{|P} dP \quad (4) $$

$$ \langle \epsilon \rangle_i = \frac{1}{|\Omega_i|} \int_{\Omega_i} \epsilon_{|P} dP \quad (5) $$

where $|\Omega_i|$ is the area (volume in 3D) of $\Omega_i$.

Obviously, in general heterogeneous media, where sometimes is not possible to define a representative volume element, and for general mechanical loads leading to averaged stresses and strains different to the ones that served to define $\hat{C}_i$, the following inequality holds:

$$ \langle \sigma \rangle_i - \hat{C}_i \langle \epsilon \rangle_i \neq 0 \quad (6) $$

The microscopic constitutive law inside the complementary domain is not considered because the nodal distribution in this region is too coarse for representing accurately the existing microstructure. In the complementary domain a homogenized constitutive law is derived, as we describe in the next paragraphs, in order to assure a mechanical well-posed problem in the whole domain. For any point $M$ inside $\Xi$, but close enough to $\Omega_i$ (Fig. 4), the estimated macroscopic law at that point involves the macroscopic stress and strain tensors $\sigma_{|M}$ and $\epsilon_{|M}$, respectively. In this case, the previous inequality also holds:
<figure><img src="image_4.png" /><figcaption>Fig. 4. Domain definition.</figcaption></figure>

$$ \sigma_{|M} - \hat{\mathbf{C}}_i \epsilon_{|M} \neq \mathbf{0} \qquad (7) $$

At convergence the following equalities are expected: $\sigma_{|M} \approx \langle \sigma \rangle_i$ and $\epsilon_{|M} \approx \langle \epsilon \rangle_i$, and Eqs. (6) and (7) being approximately equivalent, we can write:

$$ \sigma_{|M} \approx \langle \sigma \rangle_i + \hat{\mathbf{C}}_i (\epsilon_{|M} - \langle \epsilon \rangle_i) \\ \forall M \in \Xi \text{ in the neighborhood of } \Omega_i \qquad (8) $$

This constitutive law is only valid in the immediate vicinity of the microscopic patterns. Anywhere in the complementary domain a homogenized behavior is defined as we describe later.

This constitutive law takes into account the fact that in general cases the computed $\hat{\mathbf{C}}_i$ does not correspond to the exact homogenized tensor for the real loading conditions, mainly when the microscopic domain approaches the domain boundary. On the contrary, if $\hat{\mathbf{C}}_i$ represents accurately the homogenized behavior, then the non-local contributions vanish, and then the usual homogenized law is derived:

$$ \langle \sigma \rangle_i = \hat{\mathbf{C}}_i \langle \epsilon \rangle_i \Rightarrow \sigma_{|M} = \hat{\mathbf{C}}_i \epsilon_{|M} \qquad (9) $$

This result deserves some additional comments. Stress at point $M$, located in the immediate vicinity of $\Omega_i$, could be computed using the homogenized elastic tensor $\hat{\mathbf{C}}_i$, i.e. $\sigma_{|M} \approx \hat{\mathbf{C}}_i \epsilon_{|M}$. However, as $\hat{\mathbf{C}}_i$ is only exact for the 3 loadings that served to compute it (as previously described) we could improve the stress prediction by adding the quantity $\langle \sigma \rangle_i - \hat{\mathbf{C}}_i \langle \epsilon \rangle_i$ that constitutes the error that we are introducing when $\hat{\mathbf{C}}_i$ is considered in the actual loading conditions that generate the averages $\langle \sigma \rangle_i$ and $\langle \epsilon \rangle_i$. Thus, Eq. (8) is naturally derived, that represents a correction incorporated to the local use of the homogenized elastic tensor.

The simplest manner for defining the behavior at any point in the complementary domain $\Xi$ consists of using C-NEM interpolation constructed on the basis of the Voronoi diagram related to the cloud of points consisting of the center of gravity of the microscopic domains (the details concerning that interpolation will be addressed in the next section). Let's $\phi_i$ be the shape function related to the center of gravity of $\Omega_i$, which is perfectly defined in the whole domain $\Omega$. Thus, the macroscopic law results:

$$ \sigma_{|M} = \tilde{\sigma}_{|M} + \tilde{\mathbf{C}}_{|M} (\epsilon_{|M} - \tilde{\epsilon}_{|M}) \quad \forall M \in \Xi \qquad (10) $$

with:

$$ \tilde{\sigma}_{|M} = \sum_{i=1}^{i=V} \tilde{\phi}_i(M) \langle \sigma \rangle_i \qquad (11) $$

$$ \tilde{\epsilon}_{|M} = \sum_{i=1}^{i=V} \tilde{\phi}_i(M) \langle \epsilon \rangle_i \qquad (12) $$

$$ \tilde{\mathbf{C}}_{|M} = \sum_{i=1}^{i=V} \tilde{\phi}_i(M) \hat{\mathbf{C}}_i \qquad (13) $$

where $V$ denotes the number of natural neighbor microscopic domains visible from point $M$ (more details about the C-NEM approximation will be given in the next section).

Due to the interpolation property of the natural element approximation, when point $M$ approaches the microscopic domain $\Omega_i$, $\tilde{\sigma}_{|M}$, $\tilde{\epsilon}_{|M}$ and $\tilde{\mathbf{C}}_{|M}$ approach $\langle \sigma \rangle_i$, $\langle \epsilon \rangle_i$ and $\hat{\mathbf{C}}_i$, respectively, and consequently Eq. (10) reduces to Eq. (8).

The mechanical problem can be then written as: Find a regular enough displacement field in $\Omega$ verifying the essential boundary conditions

$$ U = U_c \text{ on } \partial_U \Omega \text{ with } \partial\Omega = \partial_U \Omega \cup \partial_F \Omega \qquad (14) $$

as well as the equilibrium equation

$$ \begin{aligned} & \sum_{i=1}^{i=m} \int_{\Omega_i} \mathrm{Tr}[\epsilon_{|P}(\mathbf{U}^*)\mathbf{C}(P)\epsilon_{|P}(\mathbf{U})] \mathrm{d}P \\ & + \int_{\Xi} \mathrm{Tr}[\epsilon_{|M}(\mathbf{U}^*)\tilde{\mathbf{C}}_{|M}\epsilon_{|M}(\mathbf{U})] \mathrm{d}M \\ & + \sum_{i=1}^{i=m} \frac{1}{|\Omega_i|} \int_{\Xi} \int_{\Omega_i} \tilde{\phi}_i(M) \mathrm{Tr}[\epsilon_{|M}(\mathbf{U}^*)(\mathbf{C}(P) - \tilde{\mathbf{C}}_{|M})] \mathrm{d}P \mathrm{d}M \\ & \times \epsilon_{|P}(\mathbf{U})] \mathrm{d}P \mathrm{d}M \\ & = \int_{\Omega} \mathbf{U}^* \mathbf{f} \mathrm{d}M + \int_{\partial_F \Omega} \mathbf{U}^* \mathbf{F}_c \mathrm{d}\Gamma \\ & \forall \mathbf{U}^* \text{ such that } \mathbf{U}^* = \mathbf{0} \text{ on } \partial_U \Omega \end{aligned} \qquad (15) $$

where the second and third term of the left hand member result from the virtual power principle in the complementary domain $\Xi$ making use of the behavior law expressed by Eq. (10).

### 3. Dealing with highly heterogenous nodal distributions

As previously argued, the nodal density evolves suddenly when ones moves from a microscopic domain to the complementary one. The nodal density inside the microscopic domain must be enough to represent accurately the microstructure, whereas in the complementary domain it must be enough to represent the structure response. Thus, usually, the nodal distribution becomes highly heterogenous leading to very distorted meshes if one proceeds in the finite element framework. For this reason, the use of a meshless discretization technique seems appropriate, because it is much less sensible to the aforementioned mesh distortions.

In the next paragraphs, we briefly touch upon the basis of the natural element method that will be used for discretizing the variational formulation given by Eq. (15).
The natural element method is an appealing choice among the different meshless methods, because essential boundary conditions can be imposed directly, without determination to the other properties (linear consistency, smoothness, ...). It was originally proposed by Traversoni [18], Sambridge et al. [14] and was widely investigated in elastostatic problems by Sukumar [16]. Moreover, it was successfully applied in fluid dynamic simulations involving polymer injection and metal forming simulations using updated Lagrangian formulations [12,13,1], in problems involving cracks [21] and in phase change problems involving moving discontinuities [22].

### 3.1. Constrained natural element interpolation

We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. For a more in-depth discussion on the Sibson interpolant and its application for solving second-order partial differential equations, the interested reader can refer to Braun and Sambridge [14] and Sukumar et al. [16]. The NEM interpolant is constructed on the underlying Voronoi diagram, the Delaunay tessellation being the topological dual of the Voronoi diagram.

For the sake of simplicity we only consider the 2D case, the 3D case being a direct extension.

Let $S = \{n_1, n_2, \dots, n_N\}$ be a set of nodes in $\mathbb{R}^2$. The Voronoi diagram is the subdivision of $\mathbb{R}^2$ into regions $T_i$ (Voronoi cells) defined by:

$$T_i = \{ \mathbf{x} \in \mathbb{R}^2 : d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j), \forall j \neq i \}, \quad \forall i \quad (16)$$

where $\mathbf{x}_i$ are the coordinates of node $n_i$ and $d(\cdot, \cdot)$ denotes the Euclidean distance.

The Sibson coordinates of $\mathbf{x}$ with respect to a natural neighbor $n_i$ (see Fig. 5) is defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area (volume in 3D) of the Voronoi cell associated with the point $\mathbf{x}$:

$$\phi_i(\mathbf{x}) = \frac{\text{Area}(\text{afghe})}{\text{Area}(\text{abcde})} \quad (17)$$

If the point $\mathbf{x}$ coincides with the node $n_i$, i.e. $\mathbf{x} = \mathbf{x}_i$ then $\phi_i(\mathbf{x}_i) = 1$, and all the other shape functions vanish, i.e. $\phi_i(\mathbf{x}_j) = \delta_{ij}$ ($\delta_{ij}$ being the Kronecker's delta). The proper-

<figure><img src="image_5.png" /><figcaption>Fig. 5. Construction of the Sibson shape functions.</figcaption></figure>

ties of positivity, interpolation and partition of unity are then verified [16]:

$$\left\{
\begin{array}{l}
0 \le \phi_i(\mathbf{x}) \le 1 \\
\phi_i(\mathbf{x}_j) = \delta_{ij} \\
\sum_{i=1}^{n} \phi_i(\mathbf{x}) = 1
\end{array}
\right.
\quad (18)$$

where $n$ is the number of neighbor nodes of point $\mathbf{x}$.

The natural neighbor interpolation satisfies the local coordinate property [15], namely:

$$\mathbf{x} = \sum_{i=1}^{n} \phi_i(\mathbf{x})\mathbf{x}_i \quad (19)$$

which combined with Eqs. (18), implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness).

Natural neighbor shape functions are $C^\infty$ at any point except at the nodes, where they are only $C^0$, and on the boundaries of the Delaunay circles (spheres in 3D) where they are $C^1$, because of the discontinuity in the neighbor nodes across these boundaries.

Another important property of this interpolation is its ability to reproduce linear functions on the boundary of convex domains. The proof can be found in Sukumar et al. [16], that we illustrate in Fig. 5b: due to the fact that the Voronoi cells areas associated to points on the boundary become infinite, the contribution of internal points vanishes in the limit when the point approaches the convex boundary, and the shape functions associated with nodes $n_1$ and $n_2$ become linear on the segment $(n_1 - n_2)$. This is not the case when non convex boundaries are considered. This drawback will be considered later.

Consider an interpolation scheme for a vector-valued function $U(\mathbf{x}) : \Omega \subset \mathbb{R}^2 \to \mathbb{R}^2$, in the form:

$$U^h(\mathbf{x}) = \sum_{i=1}^{n} \phi_i(\mathbf{x}) U_i \quad (20)$$

where $U_i$ are the vectors defining the nodal degrees of freedom at the $n$ natural neighbors of point $\mathbf{x}$, and $\phi_i(\mathbf{x})$ are the Sibson coordinates defined in Eq. (17) associated with each node $n_i$. It can be noticed that Eq. (20) defines a local interpolation scheme which will be used to define both the trial and the test functions considered in the discretization of different variational formulations.

A recent development in the NEM, the constrained natural element method (C-NEM), was proposed in [21,22] in order to circumvent the problems induced by non convex domains. In this approach, a visibility criterion is introduced to restrict the natural neighbors (influent nodes). For this purpose the constrained Voronoi diagram is constructed, from which the shape functions can be easily computed. In this manner, linear interpolation is recovered along the boundary of non-convex domains, making possible the introduction of essential boundary conditions as well as the treatment of fixed or moving discontinuities.
In the C-NEM framework, the interpolation can be expressed by:

$$ \mathbf{U}^h(\mathbf{x}) = \sum_{i=1}^{V} \phi_i^C(\mathbf{x}) \mathbf{U}_i \qquad (21) $$

where V is the number of natural neighbors visible from point x and $\phi_i^C$ is the constrained natural neighbor shape function, which is actually the Sibson interpolant computed using the constrained Voronoi diagram [21].

## 3.2. Elastostatic formulation

In this section, we summarize the 2D small displacement elastostatic problem which will be considered later in the numerical examples. The equilibrium equation results in:

$$ \nabla \cdot \boldsymbol{\sigma} + \mathbf{b} = \mathbf{0} \quad \text{in } \Omega \qquad (22) $$

where $\Omega \in \mathbb{R}^2$ is the material domain, $\nabla$ is the divergence operator, $\boldsymbol{\sigma}$ is the Cauchy stress tensor and $\mathbf{b}$ is a body force term.

The constitutive relation (using a vector representation of both the stress and the strain fields) is given by:

$$ \boldsymbol{\sigma} = \mathbf{C}\boldsymbol{\epsilon} \qquad (23) $$

where $\boldsymbol{\sigma}$ and $\boldsymbol{\epsilon}$ are the vector form of the stress and linearized strain tensor (symmetric part of the gradient of displacements tensor), respectively, and $\mathbf{C}$ is the elastic tensor matrix form.

The essential and natural boundary conditions are given by:

$$ \begin{aligned} \mathbf{U} &= \mathbf{U}_c \quad \text{on } \partial_U \Omega \\ \boldsymbol{\sigma}\mathbf{n} &= \mathbf{F}_c \quad \text{on } \partial_F \Omega \end{aligned} \qquad (24) $$

where $\partial\Omega = \partial_U\Omega \cup \partial_F\Omega$ is the boundary of $\Omega$, $\mathbf{n}$ is the unit outward vector defined on $\partial\Omega$, and $\mathbf{U}_c$ and $\mathbf{F}_c$ are the prescribed displacements and tractions, respectively.

The variational formulation associated with the elastostatic problem results:

Find $\mathbf{U} \in H^1(\Omega)$ kinematically admissible ($\mathbf{U} = \mathbf{U}_c$ on $\partial_U \Omega$) such that

$$ \int_{\Omega} \boldsymbol{\sigma} \cdot \boldsymbol{\epsilon}^* \mathrm{d}\Omega = \int_{\Omega} \mathbf{f} \cdot \mathbf{U}^* \mathrm{d}\Omega + \int_{\partial_F \Omega} \mathbf{F}_c \cdot \mathbf{U}^* \mathrm{d}\Gamma, \quad \forall \mathbf{U}^* \in H_0^1(\Omega) \qquad (25) $$

where $H^1(\Omega)$ and $H_0^1(\Omega)$ are the usual Sobolev functional spaces.

On substituting the trial and test functions (both approximated in the C-NEM framework) in the above equation and using the arbitrariness of the field $\mathbf{U}^*$, the following linear system of equations is obtained after numerical integration (whose details will be presented in Section 3.3):

$$ \mathbf{K}\mathbf{d} = \mathbf{f}^{\text{ext}} \qquad (26) $$

where $\mathbf{d}$ is the vector containing the nodal displacements, the matrix $\mathbf{K}$ and the vector $\mathbf{f}^{\text{ext}}$ being given by:

$$ \mathbf{K} = \int_{\Omega} \mathbf{B}^{\text{T}} \mathbf{C} \mathbf{B} \mathrm{d}\Omega \qquad (27) $$

$$ \mathbf{f}^{\text{ext}} = \int_{\partial_F \Omega} \mathbf{F}_c \mathrm{d}\Gamma + \int_{\Omega} \mathbf{N}^{\text{T}} \mathbf{f} \mathrm{d}\Omega \qquad (28) $$

where $\mathbf{N}$ and $\mathbf{B}$ are the matrix containing the shape functions and their derivatives, respectively.

## 3.3. Stabilized conforming nodal integration

Recently, new numerical integration procedures have been proposed in the framework of meshless simulations. Domain integration using Gauss quadrature introduces significant numerical errors due to the following statements: (i) the integration cells do not match shape function supports; (ii) the non-polynomial character of the meshless shape functions [6]. On the other hand, direct nodal integration, using the nodes as integration points, leads to numerical instabilities [4].

In this work, we use the integration procedure proposed by Chen et al. [3]: the stabilized conforming nodal integration (SCNI), which was applied to the NEM by Gonzalez et al. [8] and Yoo et al. [19]. In this approach, a strain smoothing is performed to stabilize the nodal integration. The SCNI is based on the assumed strain method, in which a modified gradient is introduced at the integration point (node) [3]:

$$ \tilde{\nabla} \mathbf{U}^h(\mathbf{x}_i) = \frac{1}{|\omega_i|} \int_{\omega_i} \nabla \mathbf{U}^h(\mathbf{x}) \mathrm{d}\omega \qquad (29) $$

where $\mathbf{x}_i$ are the coordinates of node $n_i$.

Thus, the modified strain vector is given by:

$$ \tilde{\boldsymbol{\epsilon}}^h(\mathbf{x}_i) = \frac{1}{|\omega_i|} \int_{\omega_i} \boldsymbol{\epsilon}^h(\mathbf{x}) \mathrm{d}\omega = \frac{1}{|\omega_i|} \int_{\omega_i} \left\{ \begin{array}{l} \frac{\partial u_1^h(\mathbf{x})}{\partial x} \\ \frac{\partial u_2^h(\mathbf{x})}{\partial y} \\ \frac{\partial u_1^h(\mathbf{x})}{\partial y} + \frac{\partial u_2^h(\mathbf{x})}{\partial x} \end{array} \right\} \mathrm{d}\Omega \qquad (30) $$

where $u_1$ and $u_2$ are the components of vector $\mathbf{U}$.

Now, by applying the divergence theorem, it results:

$$ \tilde{\boldsymbol{\epsilon}}^h(\mathbf{x}_i) = \frac{1}{|\omega_i|} \int_{\partial\omega_i} \left\{ \begin{array}{l} u_1^h(\mathbf{x})v_1 \\ u_2^h(\mathbf{x})v_2 \\ u_1^h(\mathbf{x})v_2 + u_2^h(\mathbf{x})v_1 \end{array} \right\} \mathrm{d}\Gamma \qquad (31) $$

where $\omega_i$ is the Voronoi cell related to the node $n_i$, being $\partial\omega_i$ its boundary and $|\omega_i|$ its area. $v_1$ and $v_2$ are the components of the unit outward vector defined on the cell boundary at point $\mathbf{x}$.

Introducing now the natural element approximation of the displacement Eq. (21) we obtain:

$$ \tilde{\boldsymbol{\epsilon}}^h(\mathbf{x}_i) = \tilde{\mathbf{B}}_i \mathbf{d} \qquad (32) $$

which can be expressed explicitly by:
$$ \begin{aligned} \tilde{\epsilon}^h(\mathbf{x}_i) &= \begin{bmatrix} \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial x} & 0 & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial x} & 0 & \dots & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial x} & 0 \\ 0 & \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial y} & 0 & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial y} & \dots & 0 & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial y} \\ \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial y} & \frac{\partial \tilde{\phi}_1(\mathbf{x}_i)}{\partial x} & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial y} & \frac{\partial \tilde{\phi}_2(\mathbf{x}_i)}{\partial x} & \dots & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial y} & \frac{\partial \tilde{\phi}_N(\mathbf{x}_i)}{\partial x} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ u_{N1} & u_{N2} & \dots & u_{N2} & \dots & u_{N1} & u_{N2} \end{bmatrix} \begin{cases} u_{11} \\ u_{12} \\ u_{21} \\ u_{22} \\ \vdots \\ u_{N1} \\ u_{N2} \end{cases} \\ &= \tilde{\mathbf{B}}_i d \end{aligned} \quad (33) $$

where the components of $\tilde{\mathbf{B}}_i$ are defined by:

$$ \frac{\partial \tilde{\phi}_j(\mathbf{x}_i)}{\partial x} = \frac{1}{|\omega_i|} \int_{\partial \omega_i} \phi_j(\mathbf{x}) v_1(\mathbf{x}) d\Gamma \quad (34) $$

$$ \frac{\partial \tilde{\phi}_j(\mathbf{x}_i)}{\partial y} = \frac{1}{|\omega_i|} \int_{\partial \omega_i} \phi_j(\mathbf{x}) v_2(\mathbf{x}) d\Gamma \quad (35) $$

Obviously, most of the components in the matrix $\tilde{\mathbf{B}}_i$ are zero due to the compact support of the shape functions. Introducing the C-NEM approximation, the smoothed strain and the nodal integration, the global stiffness matrix is obtained by assembling the contribution of each node $n_i$:

$$ \mathbf{K} = \sum_{i} \tilde{\mathbf{K}}_i = \sum_{i} |\omega_i| \tilde{\mathbf{B}}_i^T \tilde{\mathbf{C}} \tilde{\mathbf{B}}_i \quad (36) $$

and then, the stress in each cell can be deduced from:

$$ \tilde{\sigma}_i^h = \mathbf{C} \tilde{\mathbf{B}}_i d \quad (37) $$

which results also constant inside each Voronoi cell.

The use of the divergence theorem avoids the computation of the shape functions derivatives, requiring only the evaluation of the shape functions on the boundary of the Voronoi cells.

### 3.4. Evaluating the C-NEM ability to operate in highly heterogeneous nodal distributions

To check the C-NEM ability to operate in highly heterogeneous nodal distributions we start considering the mechanical problem reported in [10] and illustrated in Fig. 6(right) which consists of a homogeneous material subjected to the loading conditions defined in Eq. (38):

$$ \left\{ \begin{array}{l} F_x(x = 0) = \frac{2y}{L-1} \\ F_y(x = 0) = \frac{y}{L(1-y/L)} \\ F_y(x = L) = \frac{y}{L(1-y/L)} \end{array} \right. \quad (38) $$

For discretizing this problem three different nodal distributions have been considered as depicted in Fig. 6 where each microscopic domain is represented with a single point. The exact solution of this problems consists of:

$$ \left\{ \begin{array}{l} \boldsymbol{\sigma}_{11}^{\text{ex}} = \left(\frac{2y}{L} - 1\right) \left(\frac{x}{L-1}\right) \\ \boldsymbol{\sigma}_{22}^{\text{ex}} = 0 \\ \boldsymbol{\sigma}_{12}^{\text{ex}} = \frac{y}{L} \left(\frac{1-y}{L}\right) \end{array} \right. \quad (39) $$

that allows the evaluation of the relative error in $\Omega$:

$$ \eta^2 = \frac{\int_{\Omega} (\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\text{ex}}) \mathbf{C}^{-1} (\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\text{ex}}) dM}{\int_{\Omega} \boldsymbol{\sigma}^{\text{ex}} \mathbf{C}^{-1} \boldsymbol{\sigma}^{\text{ex}} dM} \quad (40) $$

For this purpose, the three clouds of nodes are considered. The second and the third one are deduced from the coarsest one by adding new microscopic patterns on the middle of the edge of each Delaunay triangle.

According to the numerical results reported in Fig. 7 the rate of convergence for both discretization techniques (C-NEM and FEM) is, as expected, of first order, being better the accuracy of the C-NEM solution. The finite element simulations have been performed using the Delaunay triangulation that serves to compute the C-NEM approximation (thus both techniques operate with the same cloud of nodes) where the usual piecewise linear finite element interpolation is defined.

It must be noticed that the exact solution cannot be exactly reproduced using linear finite element interpolations, or the C-NEM interpolation that only has linear consistency. This fact justifies the inevitable amount of error associated to the computed solution (which decreases with the mesh size) despite the homogeneous behavior here considered.

To obtain more information about the error distribution, we can also plot a local error related to each element (in the FEM context) or to each Voronoi cell (in the C-NEM framework) according to:

$$ \text{Err}_c^{\text{FEM}} = (\boldsymbol{\sigma}_c^{\text{FEM}} - \boldsymbol{\sigma}_c^{\text{ex}}) \mathbf{C}_c^{-1} (\boldsymbol{\sigma}_c^{\text{FEM}} - \boldsymbol{\sigma}_c^{\text{ex}}) \quad (41) $$

<figure><img src="image_6.png" /><figcaption>Fig. 6. Different discretizations of the test problem (each microscopic domain is represented by a single point that in fact contains a certain number of nodes).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 7. CNEM and FEM discretization errors.</figcaption></figure>

where $\sigma_c^{\text{FEM}}$ is the computed solution using the finite element method in the element e (as a linear interpolation of the displacement field is considered, the resulting stress and strain fields are constant inside the elements).

In the C-NEM context we have:

$$ \text{Err}_{c}^{\text{C-NEM}} = (\boldsymbol{\sigma}_{c}^{\text{C-NEM}} - \boldsymbol{\sigma}_{c}^{\text{ex}})\mathbf{C}_{c}^{-1}(\boldsymbol{\sigma}_{c}^{\text{C-NEM}} - \boldsymbol{\sigma}_{c}^{\text{ex}}) \quad (42) $$

where $\boldsymbol{\sigma}_c^{\text{C-NEM}}$ is the computed solution using the constrained natural element method in the Voronoi cell c (as we are using the stabilized conforming nodal integration the resulting stress and strain fields are constant inside the Voronoi cells).

The just defined local error is not affected by the volume of the element or the cell in order to be used at both scales simultaneously. We can notice in Fig. 8 an uniform

C-NEM error distribution, in contrast to the error distribution associated with the finite element resolution.

### 4. Model reduction and domain decomposition

In order to consider high performance computing making use of parallel platforms, we propose in this section an efficient domain decomposition algorithm. For this purpose we differentiate the degrees of freedom (displacements) related to nodes inside the microscopic domains and the ones related to the complementary domain (which include the nodes located on the microscopic domains boundaries).

As the resolution takes place simultaneously at both scales we will propose in the next paragraphs a model reduction operating in the microscopic domains that allows a significant reduction in the number of degrees of freedom involved in the global discretization. Moreover, we propose a strategy that extract the significant functions that must be introduced to update the microscopic approximation bases.

Let's $q_M$ be the vector containing the degrees of freedom related to the complementary domain (including those related to the microscopic domains boundaries). Let's $q_{\mu i}$ be the degrees of freedom related to the microscopic domain $\Omega_i$ (excluding its boundary). The discrete linear system can be expressed in the matrix form:

$$ \left[ \begin{matrix} \underline{K}_{MM} & \underline{K}_{M\mu 1} & \cdots & \underline{K}_{M\mu i} & \cdots \\ \underline{K}_{M\mu 1}^{\text{T}} & \underline{K}_{\mu 1\mu 1} & \cdots & 0 & \cdots \\ \vdots & \vdots & \ddots & \vdots & \ddots \\ \underline{K}_{M\mu i}^{\text{T}} & 0 & \cdots & \underline{K}_{\mu i\mu i} & \cdots \\ \vdots & \vdots & \ddots & \vdots & \ddots \end{matrix} \begin{matrix} q_M \\ q_{\mu 1} \\ \vdots \\ q_{\mu i} \\ \vdots \end{matrix} \right] = \begin{matrix} \underline{F}_M \\ \underline{F}_{\mu 1} \\ \vdots \\ \underline{F}_{\mu i} \\ \vdots \end{matrix} \quad (43) $$

<figure><img src="image_8.png" /><figcaption>Fig. 8. Error distribution related to the C-NEM and FEM solutions.</figcaption></figure>
We assume that this linear system is solved by using a particular iteration algorithm, and that at iteration $n$ the microscopic field in $\Omega_i$ is approximated using the approximation functions whose nodal expression constitute the columns of matrix $\underline{\phi}_{\underline{\mu}\underline{i}}^{(n)}$, i.e. at iteration $n$ one can express:

$$ \underline{q}_{\underline{\mu}\underline{i}}^{(n)} = \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)} \underline{a}_{\underline{\mu}\underline{i}}^{(n)} \quad (44) $$

being $\underline{a}_{\underline{\mu}\underline{i}}^{(n)}$ the coefficients of the reduced basis approximation at that iteration.

Introducing the reduced form of the microscopic fields into the linear system (43) it results:

$$ \underline{K} \begin{Bmatrix} \underline{q}_M^{(n)} \\ \underline{a}_{\underline{\mu}1}^{(n)} \\ \vdots \\ \underline{a}_{\underline{\mu}\underline{i}}^{(n)} \\ \vdots \end{Bmatrix} = \begin{Bmatrix} \underline{F}_M \\ \underline{\phi}_{\underline{\mu}1}^{(n)T} \underline{F}_{\underline{\mu}1} \\ \vdots \\ \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)T} \underline{F}_{\underline{\mu}\underline{i}} \\ \vdots \end{Bmatrix} \quad (45) $$

where

$$ \underline{K} = \begin{bmatrix} \underline{K}_{MM} & \underline{K}_{M\mu1}\underline{\phi}_{\underline{\mu}1}^{(n)} & \cdots & \underline{K}_{M\mu i}\underline{\phi}_{\underline{\mu}i}^{(n)} & \cdots \\ \underline{\phi}_{\underline{\mu}1}^{(n)T} \underline{K}_{M\mu1}^{T} & \underline{\phi}_{\underline{\mu}1}^{(n)T} \underline{K}_{M\mu1\underline{\mu}1} & \underline{\phi}_{\underline{\mu}1}^{(n)} & 0 & \cdots \\ \vdots & \vdots & \ddots & \vdots & \ddots \\ \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)T} \underline{K}_{M\mu i}^{T} & 0 & \cdots & \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)T} \underline{K}_{M\mu\underline{i}\underline{\mu}\underline{i}} & \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)} \\ \vdots & \vdots & \ddots & \vdots & \ddots \\ \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)T} \underline{K}_{M\mu i}^{T} & 0 & \cdots & \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)T} \underline{K}_{M\mu\underline{i}\underline{\mu}\underline{i}} & \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)} \end{bmatrix} \quad (46) $$

that can be rewritten in the following compact form:

$$ \begin{bmatrix} \underline{K}_{MM} & \underline{K}_{MR}^{(n)} \\ \underline{K}_{MR}^{(n)T} & \underline{K}_{RR}^{(n)} \end{bmatrix} \begin{Bmatrix} \underline{q}_M^{(n)} \\ \underline{a}^{(n)} \end{Bmatrix} = \begin{Bmatrix} \underline{F}_M \\ \underline{F}_R^{(n)} \end{Bmatrix} \quad (47) $$

We note that $\underline{K}_{MM}$ representing the macroscopic contribution to the stiffness matrix remains unchanged during the iterative resolution. So, it seems pertinent to proceed combining a block Gaussian elimination with a Cholesky factorization of $\underline{K}_{MM}$ (instead of computing directly $\underline{K}_{MM}^{-1}$).

If we choose as principal unknowns the reduced state variables $\underline{a}^{(n)}$ the previous linear system reduces to:

$$ (\underline{K}_{RR}^{(n)} - \underline{K}_{MR}^{(n)T} \underline{K}_{MM}^{-1} \underline{K}_{MR}^{(n)}) \underline{a}^{(n)} = \underline{F}_R^{(n)} - \underline{K}_{MR}^{(n)T} \underline{K}_{MM}^{-1} \underline{F}_M \quad (48) $$

Then, using $\underline{a}^{(n)}$ we can compute the nodal displacements $\underline{q}_M$ from:

$$ \underline{q}_M = \underline{K}_{MM}^{-1} (\underline{F}_M - \underline{K}_{MR}^{(n)} \underline{a}^{(n)}) \quad (49) $$

To update the reduced approximation bases that are used to approximate the solution in the microscopic domains we compute $\delta\underline{q}_{\underline{\mu}\underline{i}}^{(n+\frac{1}{2})}$ that represents a correction in the nodal degrees of freedom related to the nodes in the microscopic domains:

$$ \underline{K}_{\underline{\mu}\underline{i}\underline{\mu}\underline{i}} \delta\underline{q}_{\underline{\mu}\underline{i}}^{(n+\frac{1}{2})} = \underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)} \quad (50) $$

with

$$ \underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)} = \underline{F}_{\underline{\mu}\underline{i}} - \underline{K}_{\underline{\mu}\underline{i}\underline{\mu}\underline{i}} \underline{\phi}_{\underline{\mu}\underline{i}}^{(n)} \underline{a}_{\underline{\mu}\underline{i}}^{(n)} - \underline{K}_{M\underline{\mu}\underline{i}\underline{i}}^{T} \underline{q}_{M} \quad (51) $$

The local correction $\delta\underline{q}_{\underline{\mu}\underline{i}}^{(n+\frac{1}{2})}$ and the first three Krylov's subspaces generated by $\underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}: \underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}, \underline{K}_{\underline{\mu}\underline{i}\underline{\mu}\underline{i}}\underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}, \underline{K}_{\underline{\mu}\underline{i}\underline{\mu}\underline{i}}^{2}\underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}$ are used for defining the updated reduced basis $\underline{\phi}_{\underline{\mu}\underline{i}}^{(n+1)}$. The iteration procedure continues until reaching convergence, i.e. a residual norm $\|\underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}\|$, $\forall i$ small enough. This approach can be easily parallelized.

The use of $\delta\underline{q}_{\underline{\mu}\underline{i}}^{(n+\frac{1}{2})}$, $\underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}$, $\underline{K}_{\underline{\mu}\underline{i}\underline{\mu}\underline{i}}\underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}$, $\underline{K}_{\underline{\mu}\underline{i}\underline{\mu}\underline{i}}^{2}\underline{R}_{\underline{\mu}\underline{i}\underline{i}}^{(n)}$ as approximation functions is quite natural, but we have proved that the introduction of the three microscopic displacement fields computed in the homogenization procedure can accelerate significantly the convergence. For proving it, we consider the same problem, using seven functions (the correction field, the three first Krylov's subspace and the three displacement fields coming from the homogenization procedure) or only the four first functions. As noticed in Fig. 9 the rate of convergence related to the first strategy is two times higher than the second one.

## 5. Numerical examples and discussion

First we validate the use of very heterogeneous nodal distributions and densities. For this purpose we consider again the problem with a known solution related to a homogeneous material, described in Section 3 (Fig. 6). The domain $\Omega$ is partitioned in some microscopic patterns $(\Omega_i)_{i=1...m}$ and the complementary domain $\Xi$. Because we are considering a homogeneous material we have for each microscopic pattern $\Omega_i: \hat{\mathbf{C}} = \mathbf{C}$. So, the weak formulation associated with Eq. (15) reduces naturally to:

$$ \underline{U} = \underline{U}_c \text{ on } \partial_{\underline{U}} \underline{\Omega} \text{ with } \partial\underline{\Omega} = \partial_{\underline{U}}\underline{\Omega} \cup \partial_{\underline{F}}\underline{\Omega} $$

and

$$ \int_{\Omega} \mathrm{Tr}[\epsilon_{|M}(\underline{U}^*) \underline{C} \epsilon_{|M}(\underline{U})] \mathrm{d}M = \int_{\Omega} \underline{U}^* \underline{f}_c \mathrm{d}M + \int_{\partial_{\underline{F}}\Omega} \underline{U}^* \underline{F}_c \mathrm{d}\Gamma \\ \forall \underline{U}^* \text{ such that } \underline{U}^* = 0 \text{ on } \partial_{\underline{U}}\underline{\Omega} \quad (52) $$

The different cloud of nodes considered in the discretization of the previous variational formulation are shown in Fig. 10, where only the center of gravity of the microscopic patterns are represented. Each microscopic pattern is described by using 600 nodes. The finest model contains more than 200,000 nodes.

Fig. 11 depicts the relative error on the whole domain $\Omega$ which is compared with the one associated to the microscopic domains $\Omega_i$. The rate of convergence in the microscopic domains is higher than the global one, and similar accuracy is obtained around 4% of error.

This error is due to the field interpolation – that only has linear consistency –, as in the example considered in
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Analysis of the reduced approximation basis.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Different microscopic domains distributions.</figcaption></figure>

Section 3.4, because in the homogeneous case Eq. (8) reduces to Eq. (9) and in consequence the computed behavior everywhere in the complementary domain according to

<figure><img src="image_12.png" /><figcaption>Fig. 11. Convergence analysis: relative error in the whole domain and in the microscopic zones.</figcaption></figure>

Eq. (12) results the one given by Eq. (9) (the C-NEM interpolation having linear consistency it can reproduce exactly the constant elasticity tensor everywhere in the complementary domain).

The second simulation concerns a heterogeneous medium, depicted in Fig. 12, composed of a series of cylindrical layers (with constant thickness) with different elastic properties. This example is chosen for validating the behavior extension defined by Eqs. (8) and (10). The layers are numbered sequentially increasing with the layer radius and denoted by $L_i$. The odd layers have an isotropic elastic behavior with a Young modulus $E_1$ and a Poisson's coefficient $v_1$, being $E_2$ and $v_2$ the corresponding values in the even layers. The exact solution can be easily obtained by axial symmetry considerations:

$$U_{\text{ex}} = U_{\text{r}}(r) \mathbf{e}_{\text{r}}$$

with:

$$U_{\text{r}}(r) = A_i r + \frac{B_i}{r} \quad r \in L_i$$

where $\mathbf{e}_{\text{r}}$ is the unit vector defining the radial direction, and the constant $A_i$ and $B_i$ are computed for accounting the transmission conditions on the layers interfaces.

The domain of interest (illustrated in Fig. 12) is simply extracted from the just described problem by imposing appropriate displacements on the domain boundary (according to the exact solution).

Fig. 13 depicts the computed $C_{11}$ component of the elastic tensor. The computed equivalent Von-Mises stress is illustrated in Fig. 14 where we can observe both the fine and the coarse scale variation of the Von-Mises stress,
<figure><img src="image_13.png" /><figcaption>Fig. 12. Heterogenous problem involving two materials with different elastic properties.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Component C_{11} of the elastic tensor distribution.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Von-Mises equivalent stress.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15. Error distribution.</figcaption></figure>

the relative error, computed from Eq. (40), being 5%. The local error is depicted in Fig. 15 which proves an uniform distribution. It must be highlighted that the computed error in the microscopic domains located on the domain boundary is similar to the one computed in the domains located inside. We notice that the maximum error is concentrated in the microscopic domains boundary neighborhood possibly due to the linearity of CNEM shape functions on the interfaces, which are not rich enough to describe the stress field.

## 6. Conclusion

In this paper, we have proposed a new fully coupled two-scales modelling for mechanical problems involving microstructure. We use a meshless constrained natural element interpolation for taking into account the large and very localized variations in the nodal density. This technique allows considering simultaneously both scales, the one related to the microstructure description and the one related to its evolution (the macroscopic one). For this purpose the domain of interest is divided in some microscopic zones representing the microstructure evolution in the whole macroscopic domain, which do not cover the entire domain, and a complementary domain. The constitutive laws are defined only in the microscopic zones at the microscopic scale. The homogenized laws accurately defined at the microscopic level are then extended, using an appropriate technique, to that complementary domain.

In contrast to the vast majority of homogenization techniques, our approach allows an accurate description of the boundary conditions, because the microscopic domains can be located on the domain boundary. Moreover, to make possible the use of parallel computing, an appropriate resolution scheme was developed. This scheme is based on a domain decomposition algorithm coupled with a model reduction that allows a significant reduction in the number of degrees of freedom involved in the global discretization and which proceeds by extracting the significant functions

that must be introduced to update the approximation bases used in the microscopic domains.

Some numerical examples have been considered, allowing the comparison between the computed and the exact solutions. According to those results we can conclude:

* The ability of C-NEM to deal with highly heterogeneous nodal distributions.

* The convergence of the proposed approach.

* The efficiency of the proposed technique for extending the behavior computed in some microscopic regions to the whole domain.

* The good quality of the results, in both scales.

* The good quality of the results in the neighborhood of the domain boundaries, that allows taking into account accurately the edge effects.

The technique proposed in this work been very recent, only problems involving linear elasticity have been considered. The extension of that strategy to non-linear problems requires its formulation in the context of standard linearization procedures. This extension constitutes a work in progress.

## References

[1] I. Alfaro, J. Yvonnet, E. Cueto, M. Doblaré, F. Chinesta, Meshless methods with application to metal forming, Comput. Methods Appl. Mech. Engrg. 195 (48–49) (2006) 6661–6675.

[2] H. Ben Dhia, Problèmes mécaniques multi Echelles: la méthode Arlequin, C.R. ACAD. SCI. I-MATH. 326 (1998) 899–904.

[3] J.S. Chen, C.T. Wu, Y. Yoon, A stabilized conforming nodal integration for Galerkin mesh-free methods, Int. J. Numer. Methods Engrg. 50 (2) (2001) 435–466.

[4] J.S. Chen, S. Yoon, C.T. Wu, Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods, Int. J. Numer. Methods Engrg. 53 (2002) 2587–2615.

[5] E. Cueto, M. Doblaré, L. Gracia, Imposing essential boundary conditions in the natural element method by means of density-scaled Alpha-shapes, Int. J. Numer. Methods Engrg. 49 (2000) 519–546.

[6] J. Dolbow, T. Belytschko, Numerical integration of the Galerkin weak form in meshfree methods, Comput. Mech. 23 (1999) 219–230.
[7] C. Farhat, F.X. Roux, A method of finite element tearing and interconnecting and its parallel solution algorithm, Int. J. Numer. Methods Engrg. 32 (1991) 1205-1227.

[8] D. Gonzalez, E. Cueto, M.A. Martinez, M. Doblaré, Numerical integration in Natural neighbor Galerkin methods, Int. J. Numer. Methods Engrg. 60 (2004) 2077-2104.

[9] P. Ladeveze, O. Loiseau, D. Dureisseix, A micro-macro and parallel computational strategy for highly heterogeneous structures, Int. J. Numer. Methods Engrg. 52 (2001) 121-138.

[10] P. Ladeveze, J.P. Pelle, La maitrise du calcul en mécanique linéaire et non-linéaire, Hermès Science Publications, Paris, 2001.

[11] P. Ladeveze, A. Nouy, On a multiscale computational strategy with time and space homogenization for structural mechanics, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3061-3087.

[12] M.A. Martinez, E. Cueto, M. Doblaré, F. Chinesta, A meshless simulation of injection processes involving short fibers composites, Int. J. Forming Process. 4 (2001) 217-236.

[13] M.A. Martinez, E. Cueto, M. Doblaré, F. Chinesta, Updated lagrangian free surface flow simulations with natural neighbour Galerkin methods, Int. J. Numer. Methods Engrg. 60 (2004) 2105-2129.

[14] M. Sambridge, J. Braun, M. McQueen, Geophysical parameterization and interpolation of irregular data using natural neighbors, Geophys. J. Int. 122 (1995) 837-857.

[15] R. Sibson, A vector identity for the Dirichlet tesselations, Math. Proc. Cambridge Philosophical Soc. 87 (1980) 151-155.

[16] N. Sukumar, B. Moran, T. Belytschko, The natural elements method in solid mechanics, Int. J. Numer. Methods Engrg. 43 (1998) 839-887.

[17] N. Takano, Y. Okuno, Three-scale finite element analysis of heterogeneous media by asymptotic homogenization and mesh superposition methods, Int. J. Solids Struct. 41 (2004) 4121-4135.

[18] L. Traversoni, Natural neighbor finite elements, International Conference on Hydraulic Engineering Software Hydrosoft Proceedings, vol. 2, Computational Mechanics Publications, England, 1994, pp. 291-297.

[19] J. Yoo, B. Moran, J.S. Chen, Stabilized conforming nodal integration in the natural-element method, Int. J. Numer. Methods Engrg. 60 (1998) 861-890.

[20] J. Yvonnet, D. Ryckelynck, P. Lorong, F. Chinesta, Interpolation naturelle sur les domaines non-convexes par l'utilisation du diagramme de Voronoi contraint: méthode des éléments c-naturels, Rev. Euro. élément. Finis. 13 (2003) 487-509.

[21] J. Yvonnet, D. Ryckelynck, P. Lorong, F. Chinesta, A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM), Int. J. Numer. Methods Engrg. 60 (2004) 1451-1474.

[22] J. Yvonnet, F. Chinesta, P. Lorong, D. Rynckelynck, The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces, Int J Thermal Sci. 44 (2005) 559-569.