Quasicontinuum Method Extended to Irregular Lattices
Karel Mikeš Milan Jirásek


[
    [
    
=====


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.















§ 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 <cit.>, textile <cit.>, fibrous materials <cit.>, woven composite fabrics <cit.> or fiber composites <cit.>.
Extensive effort has been invested into the formulation of a discrete model of concrete <cit.>.

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 <cit.>, 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 <cit.>. 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 <cit.> 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 <cit.>. 
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 <cit.> and non-conservative <cit.> interaction potentials including dissipation and fiber sliding as well as for planar beam lattices <cit.>, still applied to regular lattices only.
An overview of applications and current directions of QC methods has been provided by Miller and Tadmor in <cit.> and in part IV of their book <cit.>.
In last few years, a variational formulation of dissipative QC method has been done by Rokoš at al. <cit.>, a goal-oriented adaptive version of QC algorithm has been introduced in <cit.> or a meshless QC method has been developed by Kochmann research group <cit.>.
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 <cit.>.
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 <cit.>, 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.

  * 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 <cit.>.

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.

  * 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.






§ METHODS




 §.§ 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.




 §.§ 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. 

  * 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 <cit.>. Then all newly created vertices of the mesh elements are shifted to the position of the nearest particles and are labeled as repnodes.




 §.§ 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.



  §.§.§ 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.



  §.§.§ 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 <ref>. 

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).



  §.§.§ 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 <cit.>,
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






    δ W_mac = δ W_mic

where

    δ W_mac =
        σ:δε

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

    δ W_mic =
        1/V_0∫_V_0σ_mic:δε_mic Ṿ_0

is the average virtual microscopic work density in the discrete microstructure.
Here, σ and ε are the macroscopic
stress and strain tensors, 

and σ_mic and ε_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 (<ref>) can be replaced by a sum over all links, which leads to

    δ W_mic 
        = 1/V_0∑_p=1^N_t L_p A_p σ_Np δε_Np
        = 1/V_0∑_p=1^N_tF_N_p δΔ L_p

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,
σ_Np and ε_Np is the axial stress and strain in that link,
F_Np = A_pσ_Np is the axial force and Δ L_p=L_p ε_Np
is the elongation (change of length).

In analogy to the Cauchy-Born rule used in the original quasicontinuum theory for
atomic lattices <cit.>, 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 <cit.>,
is written as

    ε_Np = _p··_p = _p:,
         10mm δε_Np = _p·δ·_p = _p:δ

where _p =_p ⊗_p is a second-order tensor, introduced for convenience.
Based on (<ref>)–(<ref>), the virtual work equality (<ref>) can be rewritten as

    σ:δε =  
        1/V_0∑_p=1^N_t L_p A_p σ_Np_p:δ

Since (<ref>)
should hold for all symmetric second-order tensors δ, and since
_p is symmetric,
the macroscopic stress must be given by 

    σ=1/V_0∑_p=1^N_t L_p A_p _pσ_Np


Formula (<ref>)  provides a rule for the evaluation of the macroscopic stress
tensor σ from the microscopic stresses, in our case
from the axial stresses in individual links, σ_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 

    σ_Np = E_p ε_Np

where E_p is the microscopic elastic modulus of link number p (often considered
as the same for all links). Substituting (<ref>) into (<ref>) and exploiting
the first part of (<ref>),
we obtain

    σ=1/V_0∑_p=1^N_t L_p A_p E_p _p⊗_p:
        =:

where

    =1/V_0∑_p=1^N_t L_p A_p E_p _p⊗_p

is  the fourth-order macroscopic elastic stiffness tensor. 

In the A3 approach, the sum in (<ref>) 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  are guaranteed because all these symmetries are exhibited by the 
fourth-order tensor , 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  
is replaced by its best approximation by an isotropic stiffness tensor characterized
by two elastic constants. The details will be explained in Section <ref>.

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. <ref>, 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 <ref>.
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.






  §.§.§ 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.
  
  * 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. <ref>c.




  * Repnode – repnode





  * 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.

  * 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 (<ref>).
The stiffness contribution is equally distributed to all elements sharing this edge.


  * Hanging node – hanging node


  * 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.

  * 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 (<ref>).
The entire contribution is assigned to the corresponding element.

  * 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 (<ref>).
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.


  * Repnode – hanging node


  * 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.

  * 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 (<ref>).
The entire contribution is assigned to the present element.

  * 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.





 §.§ Homogenization




  §.§.§ Two-dimensional models


Formula (<ref>) 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 

    𝐃^num = 
        ([ D_1111 D_1122 D_1112; D_2211 D_2222 D_2212; D_1211 D_1222 D_1212 ])

where D_ijkl are components of the material stiffness tensor .
Note that, in two dimensions, only five of these components are independent.
Matrix 𝐃^num exhibits symmetry and, on top of that, D_1122=D_1212,
because formula (<ref>) leads to

    D_1122=1/V_0∑_p=1^N_t L_p A_p E_p n_p1n_p1n_p2n_p2 =
    1/V_0∑_p=1^N_t L_p A_p E_p n_p1n_p2n_p1n_p2 =
    D_1212


If the material stiffness is considered as anisotropic, matrix 𝐃^num
is used directly. This is done by the A5 approach, with local evaluation of 𝐃^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),
𝐃^num is approximated by a matrix 

    𝐃^iso = 
        E/1-ν^2([     1     ν     0;     ν     1     0;     0     0 1-ν/2 ])

which corresponds to the isotropic material stiffness under plane stress conditions,
with E and ν 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 ν but to the same resulting matrix 𝐃^iso.

Optimal values of parameters E and ν are identified by minimizing
a certain measure of the difference between matrices 𝐃^num and 𝐃^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

    e(𝐃^num,𝐃^iso)=	
        ∑_I=1^3(𝐯_I^T( 𝐃^num -𝐃^iso)𝐯_I)^2

where 𝐯_I is the I-th eigenvector of matrix 𝐃^iso.






For this choice, the optimal parameters can be expressed explicitly as





    E   =   4(D_1111+2D_1122+D_2222)(4D_1111-7D_1122+4D_2222)/33D_1111+6D_1122+33D_2222
    ν   =   D_1111+62D_1122+D_2222/33D_1111+6D_1122+33D_2222

Note that the above expressions for E and ν 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 𝐃^num and 𝐃^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 
(<ref>)–(<ref>)
give E=(8/9)D_1111 and ν=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

    e(𝐃^num,𝐃^iso)=	
        ∑_I=1^3∑_J=1^3 W_IJ( D_IJ^num -D_IJ^iso)^2
         = (𝐃^num-𝐃^iso)𝐖(𝐃^num-𝐃^iso)

where W_IJ are suitable weight coefficients 

that can be arranged into the matrix

    𝐖 =   ([ 1 1 2; 1 1 2; 2 2 4 ])

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

    E   =   (D_1111+2D_1122+D_2222)(D_1111+14D_1122+D_2222)/2(3D_1111+12D_1122+3D_2222)
    
    	ν   =   2(D_1111-D_1122+D_2222)/3D_1111+12D_1122+3D_2222

For D_2222=D_1111 and D_1122=D_1111/3, we obtain again 
E=(8/9)D_1111 and ν=1/3. 

Alternatively, condition ν=1/3 could be imposed directly. Minimization of 
(<ref>) with Young's modulus considered as the only fitting variable would lead to

    E=8/81(4D_1111+3D_1122+4D_2222)

while minimization of (<ref>) would give

    E=2/9(D_1111+6D_1122+D_2222)




  §.§.§ Three-dimensional models


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

    𝐃^num = 
        ([ 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 ])

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, 𝐃^num is approximated by an isotropic stiffness matrix in the form 

    𝐃^iso = 
        E/(1+ν)(1-2ν)([  1-ν    ν    ν    0    0    0;    ν  1-ν    ν    0    0    0;    ν    ν  1-ν    0    0    0;    0    0    0 1-2ν    0    0;    0    0    0    0 1-2ν    0;    0    0    0    0    0 1-2ν ])

The weight coefficients used in the error measure according to formula (<ref>) are

    𝐖 = 
        ([ 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 ])

By minimizing this error we obtain

    E   =   (a+2b)(a+11b)/15a+39b
    ν   =   2a+b/5a+13b

where

    a   =    D_1111 + D_2222 + D_3333
    
    	b   =    D_2233 + D_1133 + D_1122

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 ν=1/4.



 §.§ Numerical simulation


Approaches A1–A5 described above have been implemented 
into the OOFEM open-source code <cit.>.

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.





 §.§ 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

    RSE_Ai = ∑_k R_k^(Ai) - ∑_k R_k^(A1)/∑_k R_k^(Ai)

where ∑_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).







  * The energy error indicator (EEI) is defined as

    EEI_Ai = √(1/2∑_j^NoL E_jA_jL_j (ε^(Ai)_j - ε^(A1)_j)^2)

where ε^(Ai)_j is the axial strain at link j, and the sum is taken over all links (NoL denotes the number of links).

  * The total displacement error indicator (DEI) is defined as

    DEI_Ai = √(∑_j^NoN‖𝐮^(Ai)_j - 𝐮^(A1)_j‖^2)

where 𝐮^(Ai)_j is the displacement vector at node j and NoN is the number of nodes.






§ SIMPLE TESTS OF ELASTIC RESPONSE




 §.§ Two-dimensional periodic lattices






  §.§.§ 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× L_y with
crossed diagonals; see Figure <ref>. 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 (<ref>) applied
to one single cell and combined with formula (<ref>).
The resulting stiffness matrix of a homogenized two-dimensional continuum is 

    𝐃^num   =
       2A/L_xL_yt([ E_xL_x+E_dL_dcos^4α  E_dL_dcos^2αsin^2α                   0;  E_dL_dcos^2αsin^2α E_yL_y+E_dL_dsin^4α                   0;                   0                   0  E_dL_dcos^2αsin^2α ])

where t is the out-of-plane thickness,
L_d=√(L_x^2+L_y^2) is the length of the diagonal link,
and α=arctan(L_y/L_x) is an angle characterizing the inclination of diagonals; see Fig. <ref>. In general, the stiffness matrix
given by (<ref>) corresponds to an orthotropic material.

By comparing (<ref>) with formula (<ref>) 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:

    E_xL_x+E_dL_dcos^4α   =    E_yL_y+E_dL_dsin^4α
    
    	E_dL_dcos^2αsin^2α   =   ν(E_yL_y+E_dL_dsin^4α)
    ν   =   1-ν/2




Condition (<ref>) implies that, in the case of isotropy, the macroscopic Poisson ratio 
is restricted to ν=1/3. 
From conditions (<ref>)–(<ref>)
combined with relations L_x=L_dcosα and L_y=L_dsinα we obtain

    E_x    =    E_d cosα (3-4cos^2α) 
    
        	E_y    =    E_d sinα (3-4sin^2α)

Finally, we can substitute ν=1/3 and (<ref>)–(<ref>) into condition

    E/1-ν^2 = 2A/L_xL_yt  ( E_xL_x+E_dL_dcos^4α)

and link the diagonal stiffness 

    E_d = 3tL_dE/8Asin 2α

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
(<ref>) and (<ref>)–(<ref>).
To ensure that all moduli are positive, angle α must be between  π/6 and
π/3, which means that the ratio L_x:L_y must be between 1:√(3) and √(3):1. 





  §.§.§ 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. <ref>a,
with boundary conditions according to Fig. <ref>. 
The lattice is composed of 50× 50 periodic isotropic cells.
No region of high interest is defined in these elementary tests. 



For direct tension (Fig. <ref>a), 
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. <ref>, 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 
𝐃^num given by (<ref>) 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. <ref>,
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, ∑_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. <ref>, 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

    A_e = L_x ∑_p L_y,p

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 (<ref>) always leads to the same material stiffness
as formula (<ref>), 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. 



  §.§.§ Shear test




For shear, the boundary conditions shown in Fig. <ref>b 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 (<ref>) for
meshes with different numbers of elements (NoE) are 

graphically presented in Fig. <ref>a.
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. <ref>).





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. <ref>, 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. <ref>b, 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.









































 §.§ Two-dimensional lattices with randomized cells




  §.§.§ Lattice geometry and properties


In this series of tests, 
the microstructure is obtained by randomization of a periodic microstructure composed
of 50× 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. <ref>b.

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










  §.§.§ 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 <ref>
and graphically presented in Fig. <ref>.
The energy errors (EEI) are listed in Table <ref>.
Isotropic parameters for A3 and A4 have been determined according to (<ref>).
For comparison, the RSEs are plotted again in Fig. <ref>, this time with parameters A3 and A4 determined according to (<ref>).
The isotropic stiffness is overestimated for both matrix norms.
Nevertheless, norm (<ref>) 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 <ref> (see Fig. <ref>). 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 <ref> 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. 








												





												





  §.§.§ Shear test



For the shear test of a randomized lattice, the relative stiffness error is listed
in Table <ref> and graphically presented in Fig. <ref> with A3 and A4 according to the first matrix norm,
and again in Fig. <ref> 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 <ref>.
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. 






											
																	
											





												








 §.§ 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 (<ref>) 
and (<ref>) 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.




 §.§ Three-dimensional patch tests


To check the performance in three dimensions, basic patch tests are performed on cube samples composed of 24× 24× 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 (<ref>) only.
Same as in the 2D case, no region of high interest is defined in these tests.



  §.§.§ 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. <ref>.
The energy errors (EEI) are plotted in Fig. <ref> and
the total displacement error indicator (DEI) in Fig. <ref>.

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.










  §.§.§ Shear test


For the three-dimensional  shear test, the RSE, EEI and DEI are plotted in Fig. <ref>,
Fig. <ref> and Fig. <ref>, 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.









  §.§.§ 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.







§ 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. <ref> has
dimensions 300× 300× 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.



 §.§ 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. <ref>. 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 <ref>.

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. <ref> and <ref>. Crack opening error and the energy error indicator (EEI) defined in (<ref>) are depicted in Fig. <ref>.
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 <ref>. 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_ max observed in the force-displacement diagram. The results are shown in Table <ref>. 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_ 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 <ref>. 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 <ref>. 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.


							
































































				
	

							

	


	






 §.§ 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. <ref>. For this FSD, cracking of 2000 links is computed and cracked links for approaches A1 and A5 are depicted in Figs. <ref>–<ref>
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. <ref>–<ref>.
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 ±5 to ±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

    w(s) = ⟨ 1 - s^2/R^2⟩^2

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 ⟨…⟩ 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. <ref> and <ref>, 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. <ref> and <ref> for 10 and 50 neighboring points.
The same diagrams smoothed with linear weights are compared in Fig. <ref> and <ref>.

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.






















	


























 §.§ 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.




§ 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.

99

=0.0ex
=0.0ex
=0.0ex
=0.0ex





Liu10 J. Liu, Z. Chen, and K. Li,
“A 2-d lattice model for simulating the failure of paper”,
Theoretical and Applied Fracture Mechanics, 54(1), 1-10, 2010.

BeexVer13 L. Beex, C. Verberne, and R. Peerlings,
“Experimental identification of a lattice model for woven fabrics: Application to electronic textile”,
Composites Part A: Applied Science and Manufacturing, 48, 82-92, 2013.

BeePee15 L. Beex, R. Peerlings, K. van Os and M. Geers,
“The mechanical reliability of an electronic textile investigated using the virtual-power-based quasicontinuum method”,
Mechanics of Materials, 80, 52–66, 2015.

WilBeex13 D. Wilbrink, L. Beex, and R. Peerlings,
“A discrete network model for bond failure and frictional sliding in fibrous materials”,
International Journal of Solids and Structures, 50(9), 1354-1363, 2013.

RidGon10 A. Ridruejo, C. González and J. LLorca,
“Damage micromechanisms and notch sensitivity of glass-fiber non-woven felts: An experimental and numerical study”,
Journal of the Mechanics and Physics of Solids, 58, 1628-1645, 2010.

KulaUesa12 A. Kulachenko and T. Uesaka,
“Direct simulations of fiber network deformation and failure”
Mechanics of Materials, 51, 1-14, 2012.

PengCao04 X.Q. Peng and J. Cao,
“A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics”,
Composites: Part A: Applied science and manufacturing, 36, 859-874, 2005.

Baz90 Z. P. Bažant, M. R. Tabbara, M. T. Kazemi, and G. Pijaudier-Cabot,
“Random particle model for fracture of aggregate or fiber composites”,
Journal of Engineering Mechanics, 116(8), 1686–1705, 1990.

Jin16 C. Jin, N. Buratti, M. Stacchini, M. Savoia and G. Cusatis,
“Lattice discrete particle modeling of fiber reinforced concrete: Experiments and simulations”,
European Journal of Mechanics-A/Solids, 57, 85-107, 2016.

CusPel11 G. Cusatis, D. Pelessone, and A. Mencarelli,
“Lattice discrete particle model (ldpm) for failure behavior of concrete i: Theory”,
Cement and Concrete Composites, 33(9), 881-890, 2011.

LilMie03 G. Lilliu and J. van Mier,
“3d lattice type fracture model for concrete”,
Engineering Fracture Mechanics, 70(7-8), 927-941, 2003.

LiuDen07 J. Liu, S. Deng, J. Zhang, and N. Liang. Lattice type of fracture model for concrete. Theoretical and Applied Fracture Mechanics, 48(3):269 – 284, 2007.

Baz10 Z. P. Bažant.
“Can multiscale-multiphysics methods predict softening damage and structural failure?”,
International Journal for Multiscale Computational Engineering, 8(1), 2010.


CurMil03 W. A. Curtin and R. E. Miller,
“Atomistic/continuum coupling in computational materials science”,
Modelling and Simulation in Materials Science and Engineering, 11, R33-R68, 2003.

TadPhi96 E. Tadmor, R. Phillips, and M. Ortiz,
“Mixed atomistic and continuum models of deformation in solids”,
Langmuir, 12, 4529-4534, 1996.

TadOrt96 E. B. Tadmor, M. Ortiz, and R. Phillips,
“Quasicontinuum analysis of defects in solids”,
Philosophical Magazine A, 73, 1529-1563, 1996.

TadMil05 E. Tadmor and R. Miller,
“The theory and implementation of the quasicontinuum method”,
Handbook of Materials Modeling, Springer Netherlands, 663-682, 2005.

BeePee14 L. Beex, R. Peerlings, and M. Geers,
“Central summation in the quasicontinuum method”,
Journal of the Mechanics and Physics of Solids, 70, 242-261, 2014.

BeePee14a L. Beex, R. Peerlings, and M. Geers,
“A multiscale quasicontinuum method for dissipative lattice models and discrete networks”,
Journal of the Mechanics and Physics of Solids, 64, 154-169, 2014.

BeePee2014c L. Beex, R. Peerlings, and M. Geers,
“A multiscale quasicontinuum method for lattice models with bond failure and fiber sliding”,
Computer Methods in Applied Mechanics and Engineering, 269, 108-122, 2014.

beex15 L. A. A. Beex, O. Rokoš, J. Zeman and S. P. A. Bordas,
“Higher‐order quasicontinuum methods for elastic and dissipative lattice models: uniaxial deformation and pure bending”,
GAMM‐Mitteilungen, 38(2), 344-368, 2015.

MilTad09 R. E. Miller and E. B. Tadmor,
“A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods”,
Modelling and Simulation in Materials Science and Engineering, 17, 053001, 2009.

MilTad02 R. E. Miller and E. B. Tadmor, 
“The quasicontinuum method: Overview, applications and current directions”, 
Journal of Computer-Aided Materials Design, 9, 203-239, 2002.

TadMil11 E. B. Tadmor and R. E. Miller,
“Modeling Materials: Continuum, Atomistic and Multiscale Techniques”,
Cambridge University Press, 2011.

rokos16 O. Rokoš,  L. A. Beex, J. Zeman and R. H. Peerlings,
“A Variational Formulation of Dissipative Quasicontinuum Methods”,
in press... arXiv preprint arXiv:1601.00625, 2016.


MamLar15 A. Memarnahavandi, F. Larsson and K. Runesson,
“A goal-oriented adaptive procedure for the quasi-continuum method with cluster approximation”,
Computational Mechanics, 55, 617-642, 2015.

koch14 D. M. Kochmann and G. N. Venturini,
“A meshless quasicontinuum method based on local maximum-entropy interpolation”,
Modelling and Simulation in Materials Science and Engineering, 22, 2014.

MikJir15 K. Mikeš and M. Jirásek,
The quasicontinuum method extended to disordered materials,
in Proceedings of the 15th International Conference on Civil, Structural and
Environmental Engineering Computing, Civil-Comp Press, Stirlingshire, Scotland, 2015.

oofem01 B. Patzák and Z. Bittnar
“Design of object oriented finite element code.”
Advances in Engineering Software, 32(10), 759-767, 2001.

oofem12 B. Patzák and D. Rypl
“Object-oriented, parallel finite element framework with dynamic load balancing.”
Advances in Engineering Software, 47(1), 35-50, 2012.

Pat12 B. Patzák, 
“OOFEM – an object-oriented simulation tool for advanced modeling of materials and structures”, 
Acta Polytechnica, 52, 59-66, 2012.

SheMil99 V. Shenoy, R. E. Miller, E. B. Tadmor, D. Rodney, R. Phillips, and M. Ortiz,
“An adaptive finite element approach to atomic-scale mechanics—the quasicontinuum method”,
Journal of the Mechanics and Physics of Solids, 47, 611-642, 1999.

Ryp15 D. Rypl,
T3D mesh generator,
[online]. [cit. 2015-5-1]. Available at: http://ksm.fsv.cvut.cz/dr/t3d.html.










Hill63 R. Hill,
“Elastic properties of reinforced solids: some theoretical principles”
Journal of the Mechanics and Physics of Solids 11(5), 357-372, 1963.

Man72 J. Mandel,
“Plasticité classique et viscoplasticité. International Centre for Mechanical Sciences, Courses and Lectures No. 97.” Springer, Udine 1972.


SteiEliz06 P. Steinmann, A. Elizondo and R. Sunyk,
“Studies of validity of the Cauchy-Born rule by direct comparison of continuum and atomistic modelling. Modelling and Simulation in Materials”, Science and Engineering, 15(1), S271, 2006.

Ozbolt01 J. Ožbolt, Y. Li, and I. Kožar,
“Microplane model for concrete with relaxed kinematic constraint”, International Journal of Solids and Structures, 38(16), 2683-2711, 2001.




