Quasicontinuum Method Extended to Irregular Lattices

TL;DR

This paper extends the quasicontinuum method to irregular particle systems in heterogeneous materials, introducing five approaches that improve computational efficiency while maintaining accuracy in simulating crack propagation.

Contribution

The work develops five QC-inspired methods for irregular lattices, enabling efficient and accurate simulations of heterogeneous materials with complex crack behavior.

Findings

Significant reduction in computational cost.

High accuracy in crack trajectory prediction.

Close agreement with fully resolved particle models.

Abstract

The quasicontinuum (QC) method, originally proposed by Tadmor, Ortiz and Phillips in 1996, is a computational technique that can efficiently handle regular atomistic lattices by combining continuum and atomistic approaches. In the present work, the QC method is extended to irregular systems of particles that represent a heterogeneous material. The paper introduces five QC-inspired approaches that approximate a discrete model consisting of particles connected by elastic links with axial interactions. Accuracy is first checked on simple examples in two and three spatial dimensions. Computational efficiency is then assessed by performing three-dimensional simulations of an L-shaped specimen with elastic-brittle links. It is shown that the QC-inspired approaches substantially reduce the computational cost and lead to macroscopic crack trajectories and global load-displacement curves that are very similar to those obtained by a fully resolved particle model.

Figures

Tables

Table 1. Table 1: Relative stiffness error (RSE) for a microstructure with randomized cells subjected to a uniform tension test.

		Total Error				Homogenization Error
NoE per edge	NoE	A2	A3	A4	A5	A3	A4	A5
2	4	6.48%	7.34%	7.33%	6.48%	0.86%	0.86%	0.00%
3	26	6.47%	7.34%	7.33%	6.47%	0.86%	0.86%	0.00%
4	42	6.46%	7.34%	7.33%	6.46%	0.87%	0.87%	0.00%
5	62	6.46%	7.34%	7.32%	6.46%	0.88%	0.86%	0.00%
6	86	6.44%	7.34%	7.31%	6.44%	0.89%	0.87%	0.00%
7	116	6.43%	7.34%	7.31%	6.43%	0.91%	0.88%	0.00%
8	150	6.40%	7.34%	7.30%	6.41%	0.94%	0.90%	0.01%
9	188	6.38%	7.34%	7.28%	6.39%	0.96%	0.90%	0.01%
10	232	6.33%	7.34%	7.27%	6.35%	1.00%	0.94%	0.02%
15	524	6.09%	7.34%	7.11%	6.16%	1.24%	1.02%	0.07%
20	920	5.63%	7.67%	6.73%	5.80%	2.04%	1.10%	0.18%
25	1438	5.03%	11.17%	5.76%	5.35%	6.14%	0.72%	0.32%

Table 2. Table 2: Energy error indicator (EEI) [in 10 − 5 superscript 10 5 10^{-5} ] for a microstructure with randomized cells subjected to a uniform tension test.

NoE per edge	NoE	A2	A3	A4	A5
2	4	2.939	2.943	2.942	2.939
3	26	2.936	2.943	2.940	2.936
4	42	2.932	2.943	2.939	2.932
5	62	2.930	2.943	2.942	2.930
6	86	2.922	2.943	2.940	2.922
7	116	2.917	2.943	2.942	2.917
8	150	2.904	2.943	2.943	2.904
9	188	2.894	2.943	2.947	2.894
10	232	2.875	2.943	2.947	2.875
15	524	2.770	2.943	2.988	2.770
20	920	2.575	2.946	3.078	2.575

Table 3. Table 3: Relative stiffness error (RSE) for a microstructure with randomized cells subjected to shear.

		Total Error				Homogenization Error
NoE per edge	NoE	A2	A3	A4	A5	A3	A4	A5
2	4	36.36%	21.06%	21.05%	37.28%	-15.30%	-15.30%	0.92%
3	26	21.59%	8.89%	8.93%	22.50%	-12.70%	-12.66%	0.90%
4	42	17.31%	5.33%	5.32%	18.18%	-11.98%	-11.99%	0.87%
5	62	14.57%	3.01%	2.97%	15.45%	-11.56%	-11.61%	0.88%
6	86	12.95%	1.49%	1.40%	13.74%	-11.46%	-11.55%	0.78%
7	116	12.02%	0.63%	0.58%	12.71%	-11.39%	-11.44%	0.70%
8	150	11.17%	-0.11%	-0.14%	11.78%	-11.28%	-11.32%	0.61%
9	188	10.80%	-0.43%	-0.58%	11.38%	-11.23%	-11.38%	0.58%
10	232	10.30%	-0.88%	-0.98%	10.83%	-11.18%	-11.28%	0.53%
15	524	9.25%	-1.71%	-1.95%	9.69%	-10.96%	-11.20%	0.44%
20	920	8.54%	-0.99%	-2.03%	9.03%	-9.53%	-10.57%	0.49%
25	1438	7.51%	4.38%	-0.86%	8.05%	-3.13%	-8.37%	0.54%

Table 4. Table 4: Energy error indicator (EEI) [in 10 − 5 superscript 10 5 10^{-5} ] for a microstructure with randomized cells subjected to a shear test.

NoE per edge	total NoE	A2	A3	A4	A5
2	4	4.46	4.52	4.52	4.46
3	26	2.65	2.71	2.72	2.65
4	42	2.13	2.19	2.19	2.13
5	62	1.79	1.85	1.86	1.79
6	86	1.59	1.65	1.66	1.59
7	116	1.48	1.54	1.55	1.48
8	150	1.37	1.44	1.45	1.37
9	188	1.33	1.40	1.41	1.33
10	232	1.27	1.34	1.35	1.27
15	524	1.14	1.24	1.28	1.14
20	920	1.06	1.22	1.29	1.06

Table 5. Table 5: Numbers of particles, links, elements, repnodes, hanging nodes and DOFs for various QC approaches.

	A1	A2	A3-A5
NoP	38,400	38,400	6,752
NoL	324,672	324,672	39,695
NoE	0	4,457	4,457
NoRN	38,400	5,388	5,388
NoHN	0	33,012	1,364
DOFs	113,200	16,014	16,014

Table 6. Table 6: Numbers of incorrectly cracked links.

	A2	A3i	A3a	A4	A5
Incorrectly outside FSD	4	6	13	10	7
Incorrectly cracked 1000/1000	48	60	55	86	52
Incorrectly cracked 500/1000	15	4	4	24	4

Table 7. Table 7: Tensile strength errors.

	$F_{\max }$ error	$F_{\max }$ displacement error	$F_{\max }$ step number
A1	-	-	36
A2	3.2%	-4.5%	35
A3i	1.1%	0.4%	35
A3a	2.9%	-2.1%	35
A4	1.4%	-1.1%	35
A5	4.3%	-5.7%	35

Table 8. Table 8: Time consumption of QC simplifications (i.e., of the preparatory stage before the actual simulation).

	A1	A2	A3i, A3a	A4, A5
Generate interpolation mesh	-	0.09s	0.09s	0.09s
Transform mesh to particles	-	0.19s	0.19s	0.19s
Find element for all particles	-	0.92s	0.92s	0.92s
Global stiffness tensor	-	-	0.15s	-
Connectivity table	-	-	-	0.17s
Individual stiffness tensors	-	-	-	19.27s
Total simplification time	0.0s	1.20s	1.35s	20.64s

Table 9. Table 9: Computational times consumed by various QC approaches for a conjugate gradient (CG) solver with incomplete Cholesky preconditioning and symmetric compressed column matrix storage scheme.

	Time of QC	Computational time of one step	Total comp. time
	simplification	step time / solver time	of 1000 steps
A1	-	8.2s / 5.72s (1.00 / 1.00)	2h 16m 40s (1.00)
A2	1.2s	13.66s / 0.47s (1.67 / 0.08)	3h 47m 41s (1.67)
A3i	1.35s	0.85s / 0.38s (0.10 / 0.07)	14m 11s (0.10)
A3a	1.35s	0.85s / 0.39s (0.10 / 0.07)	14m 11s (0.10)
A4	20.64s	0.85s / 0.39s (0.10 / 0.07)	14m 30s (0.11)
A5	20.64s	0.85s / 0.38s (0.10 / 0.07)	14m 30 (0.11)

Table 10. Table 10: Tensile strength errors for two different random microstructures.

	$F_{\max }$ error	$F_{\max }$ displacement	$F_{\max }$ step
	[%]	error [%]	number [%]
A1	-	-	13 / 36
A2	0.12 / 2.28	-7.91 / -6.77	11 / 15
A3i	-1.30 / 2.15	-1.13 / -0.05	32 / 17
A3a	-0.38 / 4.01	-4.61 / -2.65	32 / 18
A4	-0.77 / 0.80	-2.17 / -2.73	32 / 17
A5	0.18 / 3.31	-9.68 / -7.94	18 / 18

Equations

$$\delta W_{mac}=\delta W_{mic}$$

$$\delta W_{mac}=\mathbf{\sigma }:\delta \mathbf{\epsilon }$$

$$\delta W_{mic}=\frac{1}{V_0}{\int }_{V_0}{\mathbf{\sigma }}_{mic}:\delta {\mathbf{\epsilon }}_{mic}\mathrm{d}{V}_0$$

$$\delta W_{mic}=\frac{1}{V_0}\sum _{p=1}^{N_t}{L}_p{A}_p{\sigma }_{Np}\delta {\epsilon }_{Np}=\frac{1}{V_0}\sum _{p=1}^{N_t}{F_N}_p\delta \Delta L_p$$

$${\epsilon }_{Np}={\mathbf{n}}_p\cdot \mathbf{\epsilon }\cdot {\mathbf{n}}_p={\mathbf{N}}_p:\mathbf{\epsilon },\delta {\epsilon }_{Np}={\mathbf{n}}_p\cdot \delta \mathbf{\epsilon }\cdot {\mathbf{n}}_p={\mathbf{N}}_p:\delta \mathbf{\epsilon }$$

$$\mathbf{\sigma }:\delta \mathbf{\epsilon }=\frac{1}{V_0}\sum _{p=1}^{N_t}{L}_p{A}_p{\sigma }_{Np}{\mathbf{N}}_p:\delta \mathbf{\epsilon }$$

$$\mathbf{\sigma }=\frac{1}{V_0}\sum _{p=1}^{N_t}{L}_p{A}_p{\mathbf{N}}_p{\sigma }_{Np}$$

$${\sigma }_{Np}=E_p{\epsilon }_{Np}$$

$$\mathbf{\sigma }=\frac{1}{V_0}\sum _{p=1}^{N_t}{L}_p{A}_p{E}_p{\mathbf{N}}_p\otimes {\mathbf{N}}_p:\mathbf{\epsilon }={\mathbf{D}}_{\mathrm{e}}:\mathbf{\epsilon }$$

$${\mathbf{D}}_{\mathrm{e}}=\frac{1}{V_0}\sum _{p=1}^{N_t}{L}_p{A}_p{E}_p{\mathbf{N}}_p\otimes {\mathbf{N}}_p$$

$$\mathbf{D}^{num}=\left(\begin{array}[]{ccc}D_{1111}&D_{1122}&D_{1112}\\ D_{2211}&D_{2222}&D_{2212}\\ D_{1211}&D_{1222}&D_{1212}\end{array}\right)$$

$$D_{1122}=\frac{1}{V_0}\sum _{p=1}^{N_t}{L}_p{A}_p{E}_p{n}_{p1}{n}_{p1}{n}_{p2}{n}_{p2}=\frac{1}{V_0}\sum _{p=1}^{N_t}{L}_p{A}_p{E}_p{n}_{p1}{n}_{p2}{n}_{p1}{n}_{p2}=D_{1212}$$

$$\mathbf{D}^{iso}=\frac{E}{1-\nu^{2}}\left(\begin{array}[]{ccc}1&\nu&0\\ \nu&1&0\\ 0&0&\frac{1-\nu}{2}\end{array}\right)$$

$$e({\mathbf{D}}^{num},{\mathbf{D}}^{iso})=\sum _{I=1}^3{\left({\mathbf{v}}_I^T\left({\mathbf{D}}^{num}-{\mathbf{D}}^{iso}\right){\mathbf{v}}_I\right)}^2$$

$$E$$

$$\nu$$

$$e({\mathbf{D}}^{num},{\mathbf{D}}^{iso})=\sum _{I=1}^3\sum _{J=1}^3{W}_{IJ}{\left(D_{IJ}^{num}-D_{IJ}^{iso}\right)}^2=({\mathbf{D}}^{num}-{\mathbf{D}}^{iso})\mathbf{W}({\mathbf{D}}^{num}-{\mathbf{D}}^{iso})$$

$$\mathbf{W}=\left(\begin{array}[]{ccc}1&1&2\\ 1&1&2\\ 2&2&4\end{array}\right)$$

$$E$$

$$\nu$$

$$E=\frac{8}{81}(4D_{1111}+3D_{1122}+4D_{2222})$$

$$E=\frac{2}{9}(D_{1111}+6D_{1122}+D_{2222})$$

$$\mathbf{D}^{num}=\left(\begin{array}[]{cccccc}D_{1111}&D_{1122}&D_{1133}&D_{1123}&D_{1113}&D_{1112}\\ D_{2211}&D_{2222}&D_{2233}&D_{2223}&D_{2213}&D_{2212}\\ D_{3311}&D_{3322}&D_{3333}&D_{3323}&D_{3313}&D_{3312}\\ D_{2311}&D_{2322}&D_{2333}&D_{2323}&D_{2313}&D_{2312}\\ D_{1311}&D_{1322}&D_{1333}&D_{1323}&D_{1313}&D_{1312}\\ D_{1211}&D_{1222}&D_{1233}&D_{1223}&D_{1213}&D_{1212}\end{array}\right)$$

$$\mathbf{D}^{iso}=\frac{E}{(1+\nu)(1-2\nu)}\left(\begin{array}[]{cccccc}1-\nu&\nu&\nu&0&0&0\\ \nu&1-\nu&\nu&0&0&0\\ \nu&\nu&1-\nu&0&0&0\\ 0&0&0&1-2\nu&0&0\\ 0&0&0&0&1-2\nu&0\\ 0&0&0&0&0&1-2\nu\end{array}\right)$$

$$\mathbf{W}=\left(\begin{array}[]{cccccc}1&1&1&2&2&2\\ 1&1&1&2&2&2\\ 1&1&1&2&2&2\\ 2&2&2&4&4&4\\ 2&2&2&4&4&4\\ 2&2&2&4&4&4\end{array}\right)$$

$$E$$

$$\nu$$

$$a$$

$$b$$

$$RS{E}_{Ai}=\frac{{\sum }_k{R}_k^{(Ai)}-{\sum }_k{R}_k^{(A1)}}{{\sum }_k{R}_k^{(Ai)}}$$

$$EE{I}_{Ai}=\sqrt{\frac{1}{2}\sum _j^{NoL}{E}_j{A}_j{L}_j{\left({\epsilon }_j^{(Ai)}-{\epsilon }_j^{(A1)}\right)}^2}$$

$$DE{I}_{Ai}=\sqrt{\sum _j^{NoN}{\Vert {\mathbf{u}}_j^{(Ai)}-{\mathbf{u}}_j^{(A1)}\Vert }^2}$$

Full text

1 Introduction

Discrete particle models use a network of particles interacting via discrete links or connections that represent a discrete microstructure of the modeled material. An advantage of this approach is that discrete models can naturally capture small-scale phenomena. Therefore, a variety of sophisticated discrete material models have been developed and applied in simulations of materials such as paper [2], textile [3, 4], fibrous materials [5, 6, 7], woven composite fabrics [8] or fiber composites [9]. Extensive effort has been invested into the formulation of a discrete model of concrete [10, 11, 12, 13].

Discrete mechanical models can accurately capture complex material response, especially localized phenomena such as damage or plastic softening. However, they suffer by two main disadvantages. Firstly, a large number of particles is needed to realistically describe the response of large-scale physically relevant models. This results in huge systems of equations, which are expensive to solve. Secondly, the process of assembling of this system is also computationally expensive because all discrete connections must be individually taken into account.

Both of the aforementioned disadvantages of discrete particle models can be removed by using simplified continuous models based on one of the conventional homogenization procedures. However, standard continuous models cannot capture localized phenomena in an objective way and require enrichments, e.g., by nonlocal and gradient terms, which are again computationally expensive. According to Bažant [14], the most powerful approach to softening damage in the multi-scale context is a discrete (lattice-particle) simulation of the mesostructure of the entire structural region in which softening damage can occur.

Another way to reduce the computational cost of discrete particle models is a combination of a simplified continuous model with an exact discrete description in the parts where it is needed. Such a combination of two different approaches entails that some hand shaking procedure is needed at the interface between the continuous and discrete domains [15]. The quasicontinuum (QC) method is a suitable technique combining the advantages of continuous models with the exact description of discrete particle models without additional coupling procedures.

The quasicontinuum method was originally proposed by Tadmor, Ortiz and Phillips [16, 17] in 1996. The original purpose of this computational technique was a simplification of large atomistic lattice models described by long-range conservative interaction potentials. Since that time, QC methods have been widely used to investigate local phenomena of atomistic models with long-range interactions [15, 18]. Recently, the application of QC methods has been successfully extended to other lattices and interaction potentials. For example, an application of the QC method to structural lattice models of fibrous materials with short-range nearest-neighbour interactions has been developed by Beex et al. for conservative [19] and non-conservative [20, 21] interaction potentials including dissipation and fiber sliding as well as for planar beam lattices [22], still applied to regular lattices only. An overview of applications and current directions of QC methods has been provided by Miller and Tadmor in [23, 24, 18] and in part IV of their book [25]. In last few years, a variational formulation of dissipative QC method has been done by Rokoš at al. [26], a goal-oriented adaptive version of QC algorithm has been introduced in [27] or a meshless QC method has been developed by Kochmann research group [28]. But the application of all mentioned QC methods is still restricted only to systems with regular geometry of particles.

In the present work, we extend the QC approach to irregular systems of particles with short-range interactions by axial forces. The main idea has been tentatively presented in a conference paper [29]. Here we proceed to a more systematic evaluation of the performance of various QC formulations applied to systems with elastic-brittle links. The proposed models are implemented in OOFEM [30, 31, 32], an open-source object-oriented simulation platform initially developed for finite element methods but extensible to other discretization methods.

The procedure that results from the QC method combines the following three ingredients:

Interpolation of particle displacements is used in the regions of low interest. Only a small subset of particles is selected to characterize the behavior of the entire system. These so-called repnodes (representative nodes) are used as nodes of an underlying triangular finite element mesh, and the displacements of other particles in the region of low interest are interpolated. In the regions of high interest, all particles are selected as repnodes, in order to provide the exact resolution of the particle model. This interpolation leads to a significant reduction of the number of degrees of freedom (DOFs) without inducing a large error in the regions of high interest. 2. 2.

A summation rule can be applied in order to eliminate the requirement of visiting all particles during assembly of the global equilibrium equations. If such a rule is not imposed, all particles need to be visited to construct the system of equations, which makes the process computationally expensive. If the summation rule is adopted, the contribution of all particles in each interpolation triangle is estimated based on sampling of the links that surround one single particle and proper scaling of their contribution. This makes the computational process faster but some problems occur on the interface between regions of high and low interest. The piecewise linear interpolation of displacements combined with the summation rule means that the deformation is considered as constant within each interpolation element in the regions of low interest, while the deformations of individual links in the regions of high interest are evaluated exactly. Consequently, forces of nonphysical character, called the ghost forces, appear on the interface [23, 33].

In our work, the summation procedure is based on homogenization of link networks contributing to the interpolation elements. Some of the links (truss elements) are selected to be processed exactly, in order to properly treat the interface between the exactly solved and interpolated domains and thus to eliminate the ghost forces. 3. 3.

Adaptivity provides suitable changes of the regions of high interest during the simulation process. A new triangulation of the interpolation mesh could be done, but this is actually not necessary because the type of region can be changed by adding repnodes before each step. A suitable change of the regions of high interest often leads to a substantial increase of accuracy and, in several specific cases, it is necessary in order to represent the correct physical behavior, e.g., in a crack propagation process.

2 Methods

2.1 Overview

The original QC approach was developed for regularly arranged crystal lattices, in which atoms interact at a longer distance (not just with immediate neighbors) and the interaction forces can be derived from suitable potentials. In regions of low interest, displacements were interpolated in a piecewise linear fashion, using a selected set of representative atoms (repatoms). In this context, imposition of an affine displacement field on the periodic crystal lattice can be interpreted as an application of the Cauchy-Born rule.

In the present paper, we focus on discrete particle systems with short-range elastic or elastic-brittle interactions. Such systems are typically used in simulations of heterogeneous materials. Particles in these systems are distributed randomly and, in contrast to atomistic systems, do not form regular lattices, but the idea of QC can still be used.

Three approaches based on this idea are introduced here and are compared with the fully resolved particle model in two dimensions, which is considered as the reference case. Accuracy is assessed in terms of displacement and strain errors. The number and position of repnodes are adapted to achieve the optimal result.

The computational procedure consists of the following steps:

• •

generation of particles and of connecting links,

• •

selection of repnodes and generation of interpolation elements,

• •

application of a simplification rule,

• •

assembly of global equations with repnode displacements as basic unknowns,

• •

solution of global equations (for nonlinear models using an incremental-iterative scheme),

• •

post-processing of results and error evaluation.

The details of individual steps are described in the following subsections.

2.2 Generation of input geometry

In the first step, the input geometry of the particle system is generated; it is specified by the position of all particles in the system and by the information which pairs of particles are connected by links. This process depends on the type of represented material.

The second step consists of repnode selection and generation of interpolation elements. There are two possible reasons why a certain particle is selected as a repnode:

All particles located in a region of high interest are selected as repnodes to represent the “exact” behaviour in this region. 2. 2.

In regions of low interest, a sufficient number of repnodes are needed to construct the approximation of the displacements of other particles. Such repnodes represent vertices of the interpolation elements. The basic triangulation is done by the T3D mesh generator [34]. Then all newly created vertices of the mesh elements are shifted to the position of the nearest particles and are labeled as repnodes.

2.3 Application of simplifying rule

Once the input geometry is given, it is possible to apply a suitable simplifying rule based on the idea of QC. In this paper, five approaches using various levels of simplification are considered.

2.3.1 Pure particle approach (A1)

This approach does not use any simplification and corresponds to the reference model. Only the particles and links defining the particle model are used as input. Repnodes and interpolation elements are not needed. Every single particle represents a node with independent DOFs (displacements) and the links are described by 1D truss elements. Consequently, all links are taken into account explicitly and contribute directly to the internal forces and to the stiffness matrix.

This approach fully resolves the “exact” particle model, and the corresponding results are used as a reference solution for evaluation of accuracy and efficiency of the following simplified approaches.

2.3.2 Hanging node approach (A2)

The first technique which exploits the QC idea for simplification of the full particle model is based on approximation of DOFs of those particles that have not been selected as repnodes. Such particles are called the hanging nodes because their DOFs are not independent unknowns but are “hanging” on auxiliary elements with displacements interpolated from the neighboring repnodes. Triangular or tetrahedral interpolation elements with vertices at the repnodes are used here. For each hanging node, the corresponding interpolation element is found. It is either the element in which the hanging node is located, or the nearest element if the hanging node is not located in any interpolation element, which sometimes occurs at a curved part of the physical boundary of the particle system; see Figure 1.

The displacement of each hanging node (grey) is a linear combination of the displacements of the vertices of the corresponding interpolation element (black). This means that DOFs of all hanging nodes are interpolated (or extrapolated) using the DOFs of repnodes as the primary unknowns. Linear interpolation is used.

All links (truss elements connecting particles) contribute to the structural stiffness matrix, but only the repnodes possess independent DOFs. The repnodes represent the nodes of an interpolation mesh, which consists of linear triangular (2D) or tetrahedral (3D) elements. These elements are used only for approximation of displacements of the nodes not selected as repnodes and do not provide a direct contribution to the internal force vector and the structural stiffness matrix.

In OOFEM implementation, the particles carrying DOFs (repnodes) are modeled as regular nodes. The particles with interpolated DOFs (hanging nodes) are represented by a special type of node for which the subset of interpolation elements can be specified. The nearest interpolation element is taken from this subset and not from all interpolation elements. This technique allows to distinguish overlapping elements on the opposite sides of a crack.

The set of links that contribute to the global equilibrium equations is the same as in the A1 approach. In regions of high interest, all nodes are repnodes and the contribution of these regions is the same as in the pure particle approach (A1).

2.3.3 Global homogenization approach (A3)

In regions of low interest, this approach replaces the stiffness that corresponds to the links by the stiffness of 2D triangular or 3D tetrahedral elements “filled” by a fictitious continuum material with properties obtained by homogenization of the discrete network. Thereby, a substantial number of truss elements can be removed from the assembly procedure, and the number of operations is significantly reduced.

In the A3 approach, only one (global) effective elastic stiffness tensor, common to all elements, is assembled from the contribution of all links. Such a tensor can be derived from the Hill-Mandel condition [35, 36], which requires that the virtual macroscopic work density at a point be equal to the average microscopic virtual work in a corresponding volume $V_{0}$ of the microstructure. This condition can be written as

[TABLE]

where

[TABLE]

is the virtual macroscopic work density (i.e., work per unit volume) and

[TABLE]

is the average virtual microscopic work density in the discrete microstructure. Here, $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}$ are the macroscopic stress and strain tensors, and $\boldsymbol{\sigma}_{mic}$ and $\boldsymbol{\varepsilon}_{mic}$ are the microscopic stress and strain tensors.

For a microstructure consisting of particles connected by links that transmit axial forces only, the integral in (3) can be replaced by a sum over all links, which leads to

[TABLE]

where $N_{t}$ is the number of links in volume $V_{0}$, $L_{p}$ and $A_{p}$ is the length and cross-sectional area of link number $p$, $\sigma_{Np}$ and $\varepsilon_{Np}$ is the axial stress and strain in that link, $F_{Np}=A_{p}\sigma_{Np}$ is the axial force and $\Delta L_{p}=L_{p}\varepsilon_{Np}$ is the elongation (change of length).

In analogy to the Cauchy-Born rule used in the original quasicontinuum theory for atomic lattices [37], we will use the simplifying assumption that the microscopic strains (actual strains as well as virtual ones) can be evaluated by projecting the macroscopic strain tensor. This assumption, in microplane theory referred to as the kinematic constraint [38], is written as

[TABLE]

where $\mathbf{N}_{p}=\mathbf{n}_{p}\otimes\mathbf{n}_{p}$ is a second-order tensor, introduced for convenience. Based on (2)–(5), the virtual work equality (1) can be rewritten as

[TABLE]

Since (6) should hold for all symmetric second-order tensors $\delta\boldsymbol{\varepsilon}$, and since $\mathbf{N}_{p}$ is symmetric, the macroscopic stress must be given by

[TABLE]

Formula (7) provides a rule for the evaluation of the macroscopic stress tensor $\boldsymbol{\sigma}$ from the microscopic stresses, in our case from the axial stresses in individual links, $\sigma_{Np}$. The formula is generally applicable, even to inelastic materials. In the particular case of a linear elastic material response, the constitutive behavior of links is described by Hooke’s law

[TABLE]

where $E_{p}$ is the microscopic elastic modulus of link number $p$ (often considered as the same for all links). Substituting (8) into (7) and exploiting the first part of (5), we obtain

[TABLE]

where

[TABLE]

is the fourth-order macroscopic elastic stiffness tensor.

In the A3 approach, the sum in (10) is taken over all links of the discrete model, and $V_{0}$ corresponds to the volume of the entire domain of analysis. Major and minor symmetries of the computed stiffness tensor $\mathbf{D}_{\rm e}$ are guaranteed because all these symmetries are exhibited by the fourth-order tensor $\mathbf{N}_{p}\otimes\mathbf{N}_{p}$, for each $p$. Once the global stiffness tensor is evaluated, the corresponding material parameters are assigned to all 2D and 3D elements. The nature of these parameters depends on the assumed type of material (e.g., isotropic, orthotropic, or general anisotropic). For instance, if the material is supposed to be macroscopically isotropic, the numerically evaluated stiffness tensor $\mathbf{D}_{\rm e}$ is replaced by its best approximation by an isotropic stiffness tensor characterized by two elastic constants. The details will be explained in Section 2.4.

Triangular or tetrahedral elements, which were in the A2 approach considered as interpolation elements for hanging nodes, are now used directly for evaluation of the structural stiffness matrix, based on the material stiffness tensor obtained in the homogenization process. Thus, all hanging nodes with interpolated DOFs and all truss elements connecting them can be removed from the computational model. This leads to a significant reduction of the computational cost, but the elimination of links must be done carefully.

Links connecting two hanging nodes (located both in the same element or in two different elements) can be removed because their stiffness is represented by the effective stiffness of the homogenized material assigned to the elements. Links connecting one repnode and one hanging node located in the same element can be removed, too, because their stiffness is also reflected by the effective material stiffness.

A special case occurs if a link passes through more than one element and connects one hanging node with one repnode. This can happen if the interpolation elements are too small, or on the interface between regions of low and high interest (see Fig. 2, in which the pink rectangle is a region of high interest and the light blue rectangle is a region of low interest). Such links should not be removed because their stiffness is not reflected by the homogenized material. Therefore, the involved hanging nodes are kept and the contribution of the links is taken into account explicitly, in addition to the contribution of the triangular or tetrahedral elements.

The global effective elastic stiffness tensor is determined by techniques that will be described in Section 2.4. For a general arrangement of particles and links, this tensor is considered as anisotropic. If the internal structure of the material exhibits no preferential directions and is expected to be macroscopically isotropic, the numerically evaluated effective elastic stiffness is replaced by its best isotropic approximation. To distinguish between the general anisotropic case and the special isotropic case, the corresponding versions of the A3 approach are referred to as A3a and A3i.

2.3.4 Local homogenization approaches (A4, A5)

These approaches are refinements of A3 and aim at improving the following deficiencies:

In the A3 approach, the effective material stiffness takes into account all links. However, certain links that partially cross homogenized elements are treated explicitly, and thus a part of their stiffness is actually accounted for twice, which artificially increases the resulting structural stiffness. 2. 2.

The A3 approach assumes that, from the macroscopic point of view, material properties are the same at all points of the investigated domain, while in reality the local arrangement of links is variable across the domain.

The A4/A5 approaches remove these deficiencies by identifying effective properties of the homogenized material for each element separately. The A4 approach treats the material as isotropic while the A5 approach accounts for local anisotropy. For computation of local material parameters, it makes sense to consider general anisotropy even if the overall material behavior is isotropic. The reason is that, for small elements, the particular local arrangement of a few links can result in a significant deviation from isotropy.

Evaluation of local stiffness for all homogenized elements

The material stiffness for each element is evaluated only from the contributions of those parts of the links that are really located inside the element. Certain links are still explicitly treated as 1D truss elements and do not contribute to the stiffness of any homogenized element. In a loop over all links, the contribution of each link is computed according to the following set of rules, depending on the types of particles connected by the link. For better understanding of the rules, examples of links that correspond to individual cases (described below under labels 1ab, 2abc, 3ab) are provided in Fig. 2c.

**Repnode – repnode

**

1. (a)

If a link is located in the region of high interest:

The link is taken into account explicitly as a 1D truss element. It does not contribute to the stiffness of any homogenized element. 2. (b)

If a link is located in the region of low interest:

A link connecting two repnodes in the region of low interest is always located on an edge of an interpolation element. The stiffness contribution of this link is taken into account according to Equation (10). The stiffness contribution is equally distributed to all elements sharing this edge. 2. 2.

Hanging node – hanging node

1. (a)

If one or both ends of the link are located outside of all homogenized elements:

The link is taken into account as a truss element. It does not contribute to the stiffness of any homogenized element. 2. (b)

If both ends of the link are located in the same homogenized element:

The stiffness contribution of this link is taken into account according to Equation (10). The entire contribution is assigned to the corresponding element. 3. (c)

If the ends of the link are located in two different homogenized elements:

The stiffness contribution of this link is taken into account according to Equation (10). All elements intersected by the link are detected. The stiffness contribution is distributed to all the detected elements in proportion to the length of the part of the link inside each element. 3. 3.

Repnode – hanging node

1. (a)

If the hanging node is located outside of all homogenized elements:

The link is taken into account as a truss element. It does not contribute to the stiffness of any homogenized element. 2. (b)

If both ends of the link are located in the same element:

The stiffness contribution of this link is taken into account according to Equation (10). The entire contribution is assigned to the present element. 3. (c)

If the ends of the link are located in two different elements:

The link is taken into account as a truss element. It does not contribute to the stiffness of any homogenized element.

If a crack is modeled by overlapping elements, the same rules can be used. It is only necessary to take a decision on which side of the crack the link is located. Based on this decision, the stiffness contribution is assigned to one of the overlapping elements.

2.4 Homogenization

2.4.1 Two-dimensional models

Formula (10) is written in tensorial notation and provides the effective material stiffness tensor. In the actual numerical implementation, the tensor is represented by the corresponding matrix, based on the Voigt notation. For instance, in the two-dimensional setting, the resulting material stiffness matrix is given by

[TABLE]

where $D_{ijkl}$ are components of the material stiffness tensor $\mathbf{D}_{\rm e}$. Note that, in two dimensions, only five of these components are independent. Matrix $\mathbf{D}^{num}$ exhibits symmetry and, on top of that, $D_{1122}=D_{1212}$, because formula (10) leads to

[TABLE]

If the material stiffness is considered as anisotropic, matrix $\mathbf{D}^{num}$ is used directly. This is done by the A5 approach, with local evaluation of $\mathbf{D}^{num}$ for each homogenized element separately, and also by the A3 approach, if it is decided to use an anisotropic, globally evaluated stiffness (e.g., if the structure of the particle model is indeed anisotropic).

In the A4 approach, and also in the A3 approach with isotropic stiffness (referred to as A3i), $\mathbf{D}^{num}$ is approximated by a matrix

[TABLE]

which corresponds to the isotropic material stiffness under plane stress conditions, with $E$ and $\nu$ denoting the Young modulus and Poisson ratio of the homogenized material. Alternatively, one could adopt the plane strain assumptions, which would lead to different auxiliary values of $E$ and $\nu$ but to the same resulting matrix $\mathbf{D}^{iso}$.

Optimal values of parameters $E$ and $\nu$ are identified by minimizing a certain measure of the difference between matrices $\mathbf{D}^{num}$ and $\mathbf{D}^{iso}$. Several choices of such a measure are possible, but the results remain quite similar. The calculations presented here are based on the error measure defined as

[TABLE]

where $\mathbf{v}_{I}$ is the $I$-th eigenvector of matrix $\mathbf{D}^{iso}$. For this choice, the optimal parameters can be expressed explicitly as

[TABLE]

Note that the above expressions for $E$ and $\nu$ do not depend on stiffness coefficients $D_{1112}$ and $D_{2212}$, which are always zero for isotropic materials. If the numerically computed coefficients $D_{1112}$ and $D_{2212}$ are not small (compared to the other coefficients), the assumption of isotropy is not appropriate and a fully anisotropic stiffness should be used. Also note that a perfect matching between $\mathbf{D}^{num}$ and $\mathbf{D}^{iso}$ is possible only if coefficients $D_{1112}$ and $D_{2212}$ vanish and the other coefficients satisfy conditions $D_{2222}=D_{1111}$ and $D_{1122}=D_{1111}/3$. In this case, formulae (15)–(16) give $E=(8/9)D_{1111}$ and $\nu=1/3$.

An alternative error measure can be based on the standard Euclidean norm of fourth-order tensors. The tensorial expression $(D_{ijkl}^{num}-D_{ijkl}^{iso})(D_{ijkl}^{num}-D_{ijkl}^{iso})$, with sum over $i$, $j$, $k$ and $l$ implied by the summation convention, would be in the Voigt notation rewritten as

[TABLE]

where $W_{IJ}$ are suitable weight coefficients that can be arranged into the matrix

[TABLE]

For plane stress, the corresponding expression for the optimized elastic constants is

[TABLE]

For $D_{2222}=D_{1111}$ and $D_{1122}=D_{1111}/3$, we obtain again $E=(8/9)D_{1111}$ and $\nu=1/3$.

Alternatively, condition $\nu=1/3$ could be imposed directly. Minimization of (14) with Young’s modulus considered as the only fitting variable would lead to

[TABLE]

while minimization of (17) would give

[TABLE]

2.4.2 Three-dimensional models

In the three-dimensional setting, the resulting material stiffness matrix is given by

[TABLE]

Only fifteen components of the above matrix are independent because the coefficients $D_{ijkl}$ are invariant with respect to any permutation of the subscripts.

In the A3i and A4 approaches, $\mathbf{D}^{num}$ is approximated by an isotropic stiffness matrix in the form

[TABLE]

The weight coefficients used in the error measure according to formula (17) are

[TABLE]

By minimizing this error we obtain

[TABLE]

where

[TABLE]

In the special case of a perfectly isotropic matrix with $D_{3333}=D_{2222}=D_{1111}$ and $D_{2233}=D_{1133}=D_{1122}=D_{1111}/3$, the expressions can be simplified to $E=(5/6)D_{1111}$ and $\nu=1/4$.

2.5 Numerical simulation

Approaches A1–A5 described above have been implemented into the OOFEM open-source code [30, 31, 32]. OOFEM computes displacements of repnodes and hanging nodes and strains and stresses in truss elements and homogenized planar or spatial elements. Afterwards, some post-processing procedures are required to evaluate the error of each approach and to plot the computed results.

The OOFEM input for the A4/A5 approaches is almost the same as for A3, with only one exception. An A4/A5 input file contains a large number of materials with different parameters, which are then assigned to individual elements. Numerical tests show that computations with A4/A5 are slightly slower than with A3 because approaches A4/A5 deal with more materials. However, this difference is almost negligible. The process of identifying material parameters is visibly slower for A4/A5 and the difference depends on the type of homogenization. Even if the time of the initial set-up process is taken into account, simulations based on A3–A5 turn out to be several times faster than those based on A1 and A2.

2.6 Error measures

Accuracy of the simplifying approaches A2–A5 is evaluated by comparing the results to the exact approach, A1. The following error measures are used for that purpose:

The relative stiffness error (RSE) is measured via the relative reaction error, defined as

[TABLE]

where $\sum_{k}R_{k}^{(Ai)}$ is the sum of reactions in the loading direction at all nodes $k$ with a prescribed nonzero displacement (the simulations are performed under direct displacement control). 2. 2.

The energy error indicator (EEI) is defined as

[TABLE]

where $\varepsilon^{(Ai)}_{j}$ is the axial strain at link $j$, and the sum is taken over all links ($NoL$ denotes the number of links). 3. 3.

The total displacement error indicator (DEI) is defined as

[TABLE]

where $\mathbf{u}^{(Ai)}_{j}$ is the displacement vector at node $j$ and $NoN$ is the number of nodes.

3 Simple tests of elastic response

3.1 Two-dimensional periodic lattices

3.1.1 Lattice geometry and properties

For verification of the implemented methods, the first tests are run for regular lattices composed of periodically repeated cells. The properties of the cells are adjusted such that the resulting macroscopic behavior be isotropic.

The microstructure is generated by periodically repeating a square-shaped basic cell of size $L_{x}\times L_{y}$ with crossed diagonals; see Figure 3. The moduli of individual links in horizontal, vertical and diagonal directions are denoted as $E_{x}$, $E_{y}$ and $E_{d}$, respectively. For simplicity, the cross-sectional area $A$ of all links within a cell is considered to be the same. When multiple cells are combined into a rectangular pattern, the horizontal and vertical links on the inter-cell boundaries are merged and their resulting area is doubled.

Due to periodicity, evaluation of the homogenized stiffness can be based on formula (10) applied to one single cell and combined with formula (11). The resulting stiffness matrix of a homogenized two-dimensional continuum is

[TABLE]

where $t$ is the out-of-plane thickness, $L_{d}=\sqrt{L_{x}^{2}+L_{y}^{2}}$ is the length of the diagonal link, and $\alpha=\arctan(L_{y}/L_{x})$ is an angle characterizing the inclination of diagonals; see Fig. 3. In general, the stiffness matrix given by (33) corresponds to an orthotropic material.

By comparing (33) with formula (13) for the stiffness matrix of an isotropic material under plane stress conditions, we find that the periodic cell leads to macroscopic isotropy if the following conditions are satisfied:

[TABLE]

Condition (36) implies that, in the case of isotropy, the macroscopic Poisson ratio is restricted to $\nu=1/3$. From conditions (34)–(35) combined with relations $L_{x}=L_{d}\cos\alpha$ and $L_{y}=L_{d}\sin\alpha$ we obtain

[TABLE]

Finally, we can substitute $\nu=1/3$ and (37)–(38) into condition

[TABLE]

and link the diagonal stiffness

[TABLE]

to geometrical parameters of the cell and to the macroscopic elastic modulus.

For a cell of a given geometry and for a prescribed macroscopic elastic modulus, the characteristics of individual links in the cell can be obtained from (40) and (37)–(38). To ensure that all moduli are positive, angle $\alpha$ must be between $\pi/6$ and $\pi/3$, which means that the ratio $L_{x}:L_{y}$ must be between $1:\sqrt{3}$ and $\sqrt{3}:1$.

3.1.2 Direct tension test

The presented QC approaches are first subjected to simple tests in direct tension and shear on a regular lattice shown in Fig. 4a, with boundary conditions according to Fig. 5. The lattice is composed of $50\times 50$ periodic isotropic cells. No region of high interest is defined in these elementary tests.

For direct tension (Fig. 5a), the displacements of individual particles computed using the A1 approach (fully resolved particle model) exactly correspond to a linear displacement field that would arise in a homogeneous continuum. All horizontal links are stretched in the same way and transmit the same axial forces, and similar statements apply to the vertical links and to the diagonal links. The purpose of this test is to check whether the simplified approaches A2–A5 lead to the exact results. Indeed, this is the case, even on irregular meshes shown in Fig. 6, provided that the link stiffnesses are tuned up such that the macroscopic behavior is isotropic. The test is analogous to patch tests of finite elements, because it demonstrates that a solution with a uniform strain field is captured exactly by the numerical method. Neither interpolation nor homogenization errors arise and all approaches pass the patch test.

If the macroscopic behavior of the lattice is not isotropic, approaches A2 and A5 still lead to the exact results, while approach A3 does so only if the globally homogenized material stiffness is considered as anisotropic. For A4, the locally homogenized material stiffness is by definition isotropic, which induces a homogenization error.

The fact that the global homogenization procedure gives exactly the same stiffness matrix $\mathbf{D}^{num}$ given by (33) as the homogenization of one cell is obvious, because the whole domain is composed of an integer number of identical cells. On the other hand, it is not immediately clear that the same stiffness matrix is obtained by local homogenization over an arbitrary triangular element, including cases when edges of the element are not aligned with lattice directions (horizontal, vertical and diagonal).

A graphical proof of this interesting property is sketched in Fig. 7, which shows a typical triangular element and the underlying regular lattice. The key point is that the sum of the lengths of intersections of the triangle with vertical links, $\sum_{p}L_{y,p}$, is exactly equal to the triangle area, $A_{e}$, divided by the horizontal spacing between the horizontal links, $L_{x}$. As shown in Fig. 7, this sum multiplied by $L_{x}$ directly corresponds to evaluation of the triangle area by numerical integration using the trapezoidal rule. Since the function to be integrated is piecewise linear and one of the integration points is always located at the point where the slope changes (i.e., at the projection of one vertex), the numerical quadrature is exact, which means that

[TABLE]

An analogous statement holds for horizontal links, for ascending diagonal links, and for descending diagonal links. Consequently, the relative proportions of horizontal, vertical, ascending diagonal and descending diagonal links in each single triangular element are the same as the overall relative proportions, and formula (10) always leads to the same material stiffness as formula (33), valid for one periodic cell. This would not be the case if the vertices of the triangle were placed at arbitrary locations and not at grid points.

3.1.3 Shear test

For shear, the boundary conditions shown in Fig. 5b do not lead to a uniform strain field (note that the top and bottom parts of the boundary are considered as traction-free and their displacements are not prescribed). The relative stiffness errors evaluated according to formula (30) for meshes with different numbers of elements (NoE) are graphically presented in Fig. 9a. The error decreases fast with increasing number of elements, which indicates that the error is due to interpolation. The homogenization procedure, even when performed for each triangular element separately, does not introduce any error, for the same reasons as explained in the previous subsection on direct tension (see Fig. 7).

An exception is the finest interpolation mesh with 26 elements per edge, for which the elements are not much bigger than the lattice cells. As illustrated in Fig. 8, for such fine meshes there exist links (marked in red) that connect a repnode with a particle located in an interpolation element not connected to that repnode. Such links are accounted for explicitly and are excluded from homogenization, which leads to disturbances. The A5 approach with locally evaluated anisotropic material stiffness is then superior to the A3 and A4 approaches, which assume isotropy.

In terms of the stiffness error (RSE), the A2 approach (hanging nodes) gives in the present example somewhat higher accuracy than the other approaches, even though the difference is not dramatic. As shown in Fig. 9b, the energy error indicator (EEI) has the same value for all approaches (A2–A5) and decreases to very low levels as the interpolation mesh is refined.

3.2 Two-dimensional lattices with randomized cells

3.2.1 Lattice geometry and properties

In this series of tests, the microstructure is obtained by randomization of a periodic microstructure composed of $50\times 50$ cells. The randomization is achieved by random shifts of particle positions. The maximal shift values in $x$ and $y$ directions are $0.3L_{x}$ and $0.3L_{y}$, respectively. The nodes located on edges are shifted along the edges only. The final result is evaluated as an average of computations with five different randomized microstructures. An example of a randomized microstructure is plotted in Fig. 4b.

Patch tests have been performed for a number of triangular meshes with various sizes of elements; see Fig. 6.

3.2.2 Direct tension test

For the direct tension test, the relative stiffness errors (RSE) corresponding to different meshes and to approaches A2–A5 are listed in Table 1 and graphically presented in Fig. 10. The energy errors (EEI) are listed in Table 2. Isotropic parameters for A3 and A4 have been determined according to (14). For comparison, the RSEs are plotted again in Fig. 11, this time with parameters A3 and A4 determined according to (17). The isotropic stiffness is overestimated for both matrix norms. Nevertheless, norm (17) provides a significantly higher error and seems to be unacceptable for application to uniform tension. The error induced by the A3 approach is independent of the refinement of the interpolation mesh, except for very fine meshes. The reason is that the A3 approach deals with the same homogenized stiffness matrix in all elements, and if the whole sample is replaced by a homogeneous continuum, the prescribed boundary conditions induce uniform strain across the sample. Such a state of uniform strain is then captured by all meshes because the underlying finite elements pass the standard patch test. Of course, the reference solution obtained using the fully resolved particle model is somewhat different, because the microstructure is not completely regular. The error of the A3 approach is thus nonzero; it does not depend on the mesh, with the exception of very fine meshes, for which some links are considered explicitly, as already discussed in Section 3.1.3 (see Fig. 8). For approaches A4 and A5, the total error slowly decreases as the mesh is refined. As expected, A4 is seen to be more accurate than A3, and A5 is still more accurate.

The last three columns in Table 1 present the homogenization errors for approaches A3–A5. The homogenization error is defined as the difference between the total error of the given approach and the total error of the A2 approach, which does not use any homogenization. The homogenization errors are seen to be below 1%, and for the A5 approach to be virtually nil, with the exception of very fine meshes. Even for coarser meshes, the homogenization error slightly increases with mesh refinement. For approaches A3 and A4, this can be attributed to the anisotropic character of the local arrangements of links. The deviation from isotropy is more pronounced on fine interpolation meshes.

3.2.3 Shear test

For the shear test of a randomized lattice, the relative stiffness error is listed in Table 3 and graphically presented in Fig. 12 with A3 and A4 according to the first matrix norm, and again in Fig. 13 with A3 and A4 according to the second matrix norm. For this shear test, the results of comparison of different matrix norms is totally opposite to the results for tension. The second matrix norm predicts correct results comparable with other approaches, whereas the stiffness predicted with the first norm is significantly underestimated and the total stiffness error can even become negative, which indicates that the response is too compliant. The energy error is listed in Table 4. It is seen that A5 gives almost the same total error as A2, and the error decreases with refinement of the interpolation mesh. On the other hand, approaches A3 and A4 give seemingly a smaller relative stiffness error, which can even become negative. Their energy error is larger than for A2 and A5.

Most of the error is due to interpolation. The homogenization error is smaller, and for the A5 approach it is negligible.

3.3 Two-dimensional tests – conclusions

In partial summary, the tensile test of a regular lattice serves as a patch test that verifies that the implementation is correct because all the approaches reproduce the exact solution with no error. The shear test of a regular lattice shows that, upon refinement, the error tends to zero.

For the tensile test of a randomized lattice, the error in stiffness caused by interpolation remains above 5% even on very fine meshes. For the shear test, it remains above 7.5%. This is the intrinsic error that needs to be accepted. Isotropic homogenization in some cases leads to an increase of compliance, which counteracts the increase of stiffness due to interpolation. The resulting total error in stiffness is thus in certain cases near zero or even negative, but the energy error always remains positive.

Differences in the performance of homogenization based on error measures (14) and (17) in different patch tests indicate that isotropic homogenization of anisotropic materials can be dangerous. Therefore, it is better to avoid approaches A3i and A4 if the homogenized microstructure is strongly anisotropic.

3.4 Three-dimensional patch tests

To check the performance in three dimensions, basic patch tests are performed on cube samples composed of $24\times 24\times 24$ cells. All presented QC approaches are again applied to direct tension and shear. The initial periodic 3D microstructure is randomized in the same way as in 2D. The final result is evaluated as an average of computations with five different randomized microstructures. Parameters of isotropic stiffness (for the A3i and A4 approaches) are obtained by using matrix norm (17) only. Same as in the 2D case, no region of high interest is defined in these tests.

3.4.1 Direct tension test

For the three-dimensional direct tension test, the relative stiffness errors (RSE) corresponding to different meshes and to approaches A2–A5 are plotted in Fig. 15. The energy errors (EEI) are plotted in Fig. 16 and the total displacement error indicator (DEI) in Fig. 17.

In terms of the RSE, all approaches exhibit the same behavior as in two dimensions. The only difference is that the homogenized isotropic stiffness is underestimated instead of overestimated. The EEI of A3i and A3a remains constant for all mesh sizes (except the finest mesh). This confirms the fact that, for a uniform displacement field, homogenization used by these approaches is independent of the refinement of the interpolation mesh. On the other side, the EEI of A4 increases with mesh refinement because smaller elements are statistically more anisotropic. The order of DEI reflects the quality of homogenized procedures used in all QC approaches.

3.4.2 Shear test

For the three-dimensional shear test, the RSE, EEI and DEI are plotted in Fig. 15, Fig. 16 and Fig. 17, respectively. The performance of all approaches is quite similar. The response of A3 and A4 is too compliant due to homogenization. Therefore these approaches appear to be more accurate than A2 and A5 in terms of RSE but not in terms of EEI and DEI.

3.4.3 Three-dimensional tests – conclusions

Three-dimensional patch tests have shown that local anisotropic homogenization (i.e., the A5 approach) provides results with an almost zero homogenization error, but an intrinsic stiffness error around 15% due to interpolation must be taken into account.

4 Failure simulations

To assess the efficiency and accuracy of QC-based approaches in applications that involve inelastic material response, failure of an L-shaped specimen is simulated in three dimensions. The behavior of links that connect particles is assumed to be elastic-perfectly brittle, with link breakage occurring at a critical level of tensile strain.

The L-shaped specimen shown in Fig. 21 has dimensions $300\times 300\times 100$ mm. It is fixed at the bottom section and loaded by prescribed vertical and horizontal displacements imposed at the left end section. As a result, the non-convex corner is opened, the upper part of the specimen is bent and the horizontal part is twisted.

The microstructure is generated with a density of 20 nodes along the short edge and randomized. This procedure results in 38,400 particles (with 113,200 unknown DOFs) connected by 324,672 links. Material parameters are considered to be the same for all links.

In the simulation, the prescribed displacement is increased until the critical value of tensile strain is reached in the most stretched link. Subsequently, the link is removed and the loading is imposed again. Repeating this process results in a series of link breakages that define the macroscopic crack trajectory.

4.1 Cylindrical fully resolved domain

For the purpose of this test, a cylindrical fully resolved domain (FSD) is specified around the non-convex corner; see Fig. 21. The FSD radius is set to 50 mm, i.e., to one half of the length of the shortest specimen edge. This FSD occupies 11.8 % of the volume of the entire domain. Basic characteristics of models used by simplified QC approaches are listed in Table 5.

The simulation proceeds until 1000 links are cracked. Relative errors evaluated by comparing the results of simplified approaches to the full particle model (approach A1) are depicted in Fig. 22 and 23. Crack opening error and the energy error indicator (EEI) defined in (31) are depicted in Fig. 24. The ranking of individual QC approaches according to EEI in the first step reflects the quality of homogenization methods used by these approaches. On the other side, in the last step, the ranking of approaches is different and depends on which links are cracked.

The capability of all simplified approaches A2-A5 to correctly predict the link that will break next (i.e., the link with maximum strain) depends on the current microstructure. The numbers of incorrectly cracked links during the failure process are compared in Table 6. The numbers in the line denoted as $500/1000$ indicate how many of the first 500 cracked links for the given approach are not found in the list of the first 1000 cracked links in the exact solution (approach A1). For simplified approaches, the precise sequence of cracked links is not always the same but just a small number of cracked links are predicted incorrectly. Even though a few incorrectly cracked links appear, the simplified approaches are able to predict a correct overall crack trajectory, provided that the FSD is selected properly.

Another important quantity is the maximum value of loading force $F_{\rm max}$ observed in the force-displacement diagram. The results are shown in Table 7. For all simplified approaches, the maximum force is predicted in the 35th step, while the full particle model (approach A1) gives the maximum force in the 36th step. However, relative errors in $F_{\rm max}$ as well as relative errors in the corresponding displacement are just a few percent.

Computational times consumed by individual components of various QC-based approaches are summarized in Table 8. Two most demanding procedures are the assignment of interpolation elements to all particles and the assembly of individual stiffness tensors from all elements. Searching for the interpolation element is done independently for each particle. The distribution of link stiffnesses to individual tensors is also done independently for each link. Therefore, parallelization of the loop over particles or links can easily be envisaged with an almost ideal expected speed-up.

The computational times of one step and of the whole simulation are shown in Table 9. The QC approaches are able to reduce the computational time of one step more than ten times. Even if the computational time needed for the initial simplification is high, the total simulation time is significantly reduced if the total number of steps is huge.

4.2 Fully resolved domain for crack propagation

The cylindrical FSD used in the previous section is not convenient for simulations of a long crack trajectory. Even if cracking of individual links outside the FSD is implemented, propagation of the crack outside the FSD is not accurate. To be able to predict the correct crack trajectory, it is necessary to change the FSD adaptively or set up a priori an FSD in the region where crack propagation is expected. For that purpose a new wedge-shaped FSD is selected; see Fig. 25. For this FSD, cracking of 2000 links is computed and cracked links for approaches A1 and A5 are depicted in Figs. 26–27 for two random realizations of the underlying particle model. The list of cracked links for the A5 approach is not the same as for the exact approach A1, but the macroscopic shape of the crack is quite similar to the exact solution.

The force-displacement diagram corresponding to a crack growth simulation in a microscopically elastic-brittle material exhibits high oscillations; see the scattered gray lines in Figs. 28–29. Interpretation and comparison of such force-displacement diagrams is facilitated by smoothing the results. The smoothing procedure replaces the values of force and displacement in each characteristic point by their weighted averages around this point. Different numbers of neighboring points have been used, ranging from $\pm$5 to $\pm$100. Three types of weight functions have been considered, namely constant, linear (closer points have a stronger influence) and bell-shaped polynomial (in the form often used by nonlocal material models) given by

[TABLE]

where $s$ is one plus the number of points between the averaged point and the point for which the weight is evaluated. Parameter $R$ equals to one plus the number of considered neighboring points, and $\langle\ldots\rangle$ are Macauley brackets denoting the positive part. The results for different numbers of neighboring points are compared for constant and linear weight functions in Fig. 28 and 29, respectively. The results obtained with the bell-shaped weight function are very similar to the results with linear weights. Constant weights lead to sharper shapes of the final diagrams and, for the same number of used neighboring points, the oscillations are smaller in comparison with linear or polynomial weights. All three variants of smoothing reflect the character of the original diagram and reduce oscillations.

Diagrams for all approaches smoothed with constant weights are compared in Fig. 30 and 31 for 10 and 50 neighboring points. The same diagrams smoothed with linear weights are compared in Fig. 32 and 33.

In accordance with patch tests, the initial elastic response of the simplified approaches is stiffer than the exact solution. However, the shape of the softening branch is captured by all approaches very well.

4.3 Large-scale computations

The maximum size of the problem to be solved is limited by the size of available memory. Typically, today’s office PCs come with 4GB of RAM. For such a computer, the maximum size of the L-shaped specimen that can be solved by OOFEM depends on the fineness of the microstructure. For a direct solver with a symmetric skyline matrix storage format, the maximum fineness of microstructure that can be solved using the pure particle model (approach A1) is 27 particles along the shortest edge, which leads to 282,852 DOFs and 820,846 links. The more powerful conjugate gradient iterative solver with a dynamically growing compressed column matrix storage format is able to solve examples with 46 particles along the shortest edge and with 1,424,068 DOFs and 4,190,040 links. Larger examples cannot be run on a single standard PC due to lack of memory.

QC-based approaches can handle particle models with microstructure density up to 89 particles along the shortest edge using the same cheap PC with 4GB of RAM. The full particle model has 3,493,161 particles, 31,007,592 links and 10,439,878 DOFs. The QC approach uses a fully resolved domain with 433,265 repnodes and 3,694,490 links that cover 11.8 % of the solved domain. In the remaining part of the domain, an interpolation mesh consisting of 4,457 elements and 1,098 mesh nodes is used, which results into a total of 3,698,947 elements and 1,302,939 DOFs. In this simplified solution, the link with the maximum strain is predicted correctly and the relative error in maximum strain is just a few percent.

Modern supercomputers would allow to solve problems with a finer microstructure but such computations can be quite expensive. Furthermore, even for supercomputers, a finite limit on the size of the problem always exists and QC-based approaches can make even larger systems solvable.

5 Conclusions and future work

The presented example has demonstrated that QC-based methods can lead to a substantial reduction of the computational cost. The error induced by this reduction can be kept within acceptable limits by suitably setting the region of high interest (fully resolved domain, FSD).

Macroscopic properties associated with the global stiffness are naturally affected by a certain error induced by interpolation. Local phenomena such as cracking are well captured by sufficiently large FSDs.

Approaches A3i and A4 based on isotropic homogenization tend to underestimate the global stiffness if the material is significantly anisotropic. In such cases, these approaches seemingly appear to be more accurate in certain examples but accuracy and convergence of this approach are not guaranteed. On the other hand, the A5 approach based on local anisotropic homogenization seems to be very powerful. The homogenization error of A5 is negligible and this approach provides almost as accurate results as A2 while running substantially faster.

Both elastic and simple inelastic material models have been presented. So far, the inelastic behavior has been considered to have the form of brittle failure on the microscopic level. Future work will deal with optimization of efficiency and extensions to elastoplasticity and to softening material response, e.g., to damage-based models.

Acknowledgement

Financial support received from the Czech Science Foundation (GAČR project No. 14-00420S) is gratefully acknowledged.