Acta Polytechnica CTU Proceedings doi:10.14311/APP.2017.13.0080 Acta Polytechnica CTU Proceedings 13:80–84, 2017 © Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/app ADAPTIVE QUASICONTINUUM SIMULATION OF ELASTIC-BRITTLE DISORDERED LATTICES Karel Mikeš∗, Milan Jirásek Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: karel.mikes.1@fsv.cvut.cz Abstract. The quasicontinuum (QC) method is a computational technique that can efficiently handle atomistic lattices by combining continuum and atomistic approaches. In this work, the QC method is combined with an adaptive algorithm, to obtain correct predictions of crack trajectories in failure simulations. Numerical simulations of crack propagation in elastic-brittle disordered lattices are performed for a two-dimensional example. The obtained results are compared with the fully resolved particle model. It is shown that the adaptive QC simulation provides a significant reduction of the computational demand. At the same time, the macroscopic crack trajectories and the shape of the force-displacement diagram are very well captured. Keywords: quasicontinuum method, adaptive crack propagation, elastic-brittle lattice model. 1. Introduction Discrete particle models can effectively capture com- plex material responses, especially localized phenom- ena such as damage or plastic softening. The main dis- advantage of particle-based approaches is that a huge number of particles is needed to describe the response of large-scale physically relevant models. This results in an extensive system of equations, which is expen- sive to solve. Furthermore, the process of assembling this system is also computationally expensive because all discrete connections must be individually taken into account. A quasicontinuum (QC) based method can remove both of these disadvantages by combining continuum and atomistic approaches. In order to simulate crack propagation, the QC method needs to be combined with an adaptive algorithm that allows crack growth in arbitrary directions and initialization of new cracks. In this work, the material is represented by particles interacting via elastic-brittle links forming a disor- dered two-dimensional lattice. Only axial interaction between particles is considered and the behavior of links is assumed to be perfectly elastic-brittle, with link breakage occurring at a critical level of tensile strain. 2. QC method The QC method can efficiently handle high-resolution particle models by combining continuum and discrete approaches. This method was originally proposed in 1996 [1] for regular atomistic lattices with long-range interactions. Since that time, the QC method has been widely used and extended to applications for a variety of materials represented by regular lattices [2]. An extension of the QC method to irregular lattices has recently been developed by the authors [3]. The main idea of the QC method is to reduce the number of degrees of freedom (DOFs) and the asso- ciated computational cost without losing the exact atomistic description in regions of interest. Therefore, two types of regions in the investigated domain are considered. In regions of high interest, the pure parti- cle approach is used and all particles carry their own independent DOFs. By contrast, in regions of low interest, continuum assumptions are applied and the computational model is significantly simplified. The procedure that results from the QC method can be briefly presented in the following three steps: (1.) Interpolation (2.) Summation (3.) Adaptivity 2.1. Interpolation Interpolation of DOFs is used in regions of low interest. To simplify the full particle model, only a small subset of particles is selected to represent the entire system. These so-called repnodes serve as the nodes of a back- ground triangular mesh that is used to interpolate the DOFs of other particles in regions of low interest. On the other hand, in regions of high interest, all particles are selected as repnodes to provide the exact particle representation. This interpolation leads to a significant reduction of the number of DOFs without losing the exact particle description in regions where a high resolution is needed. 2.2. Summation The interpolation provides a significant reduction of the number of DOFs but all particles still need to be visited to construct the system of governing equations, which makes the process computationally expensive. 80 http://dx.doi.org/10.14311/APP.2017.13.0080 http://ojs.cvut.cz/ojs/index.php/app vol. 13/2017 Adaptive Quasicontinuum Simulation of Disordered Lattices The summation rule can be applied in order to elimi- nate the requirement of visiting all particles during the assembly of the system. 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 properly scaling the contribution. This makes the computational process faster but some problems can occur on the interface between regions of high and low interest. Because of the interpolation and the summation, the deformation is considered as constant within each in- terpolation 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 between the regions of low and high interest [4]. In this work, the summation procedure is realized by a homogenization of links contributing to the in- terpolation elements. To eliminate the ghost forces, some specific links are selected to be processed ex- actly, in order to capture the interface between the fully resolved and interpolated domains. 2.3. Adaptivity Adaptivity allows to adjust the regions of low and high interest during the simulation process. The type of region can be changed by adding repnodes before each step. An appropriate modification of the regions of high interest leads to a substantial increase of ac- curacy. Moreover, in specific cases such as simulation of crack propagation or damage evolution, adaptivity is necessary in order to represent the correct phys- ical behavior. After the insertion of new repnodes, a new triangulation of the interpolation mesh can be constructed, in order to improve the accuracy of interpolation and summation in areas of low interest. 3. QC-based approaches In this section, the exact approach and two approaches based on the idea of QC with different levels of sim- plification are introduced. 3.1. Fully resolved approach This approach does not use any simplification. Every single particle represents a node with independent DOFs. All links are taken into account explicitly and contribute directly to the stiffness matrix. Conse- quently, this approach provides the “exact” result, which is used as a reference solution for evaluation of accuracy and efficiency of the following, simplified approaches. 3.2. QC approach with interpolation In this approach, only the interpolation rule is used to simplify the full particle model. In the regions of high interest, all particles are selected as repnodes. By contrast, in the area of low interest, only the parti- cles forming the vertices of interpolation elements are selected as repnodes. All remaining particles are kept in the model and their DOFs are linearly interpolated by using the underlying interpolation mesh. Such particles are called hanging nodes because their DOFs are not independent but are “hanging” on appropriate repnodes. All link connections are considered explicitly in the whole domain. 3.3. QC approach with interpolation and homogenization In this approach, 2D finite elements are used not only to interpolate DOFs but also to replace the stiffness that corresponds to the microstructure. Consequently, a substantial number of links and hanging nodes can be removed from the particle model. Material properties of 2D elements are identified by homogenization of the effective elastic stiffness ten- sor De representing the microstructure of the links. Applying the equivalence of the overall virtual work expressed for the continuous material and for the microstructure represented by discrete links, the fol- lowing formula for the effective elastic stiffness tensor can be derived: De = 1 V Nt∑ i=1 LiEAini ⊗ ni ⊗ ni ⊗ ni (1) Here, ni is the unit vector specifying the direction of the given link, E is the Young modulus, A is the cross-sectional area and L is the length of the link. The summation is taken over Nt links occupying vol- ume V . The effective stiffness tensor De can be evaluated globally for all elements, or locally for each element separately. Afterwards, the material parameters of 2D elements can be identified as isotropic, orthotropic, or arbitrarily anisotropic. Different homogenization procedures for disordered lattices with normal inter- actions are described in [3]. In this work, only the most accurate homogeniza- tion rule with a local anisotropic stiffness tensor and an arbitrarily anisotropic material is used. Evaluation of the effective material stiffness tensor given by (1) is done for each element separately. The stiffness tensor of each element is obtained only from the contribu- tions of the parts of the links that are located in that particular element. 4. Adaptive algorithm Adaptivity is one of the key ingredients of the QC method. The area of high interest and the geome- try of interpolation mesh can be arbitrarily changed during the simulation to optimize the accuracy and computational costs. Furthermore, for QC simulation of materials with softening or damage, an adaptive 81 Karel Mikeš, Milan Jirásek Acta Polytechnica CTU Proceedings algorithm is required to be able to predict the exact failure mechanism. In this paper, it is assumed that cracked links tend to form a macroscopic crack. The QC adaptive algo- rithm is designed to be able to simulate the growth of existing macroscopic crack as well as the initialization of new cracks. In each simulation step, this algorithm is able to change both the area of high interest and the interpolation elements in the area of low interest. The overall algorithm of an adaptive CQ simulation is organized as follows: Algorithm 1 Adaptive QC simulation 1: load input geometry 2: generate interpolation mesh 3: calculate local stiffness tensors and material pa- rameters 4: assemble stiffness matrix 5: assemble load vector 6: repeat 7: solve one elastic step 8: break link with maximal strain 9: update area of high interest 10: update interpolation mesh 11: update local stiffness tensors and material pa- rameters 12: update stiffness matrix and load vector 13: until 14: postprocessing 4.1. Area of high interest update The change of the area of high interest is realized around newly cracked links to allow crack growth. When a new link is broken, all particles within a given updating radius are labeled as repnodes. This update is realized only if the breaking link is located not too far from existing macroscopic cracks. Therefore, a number of previously cracked links within a checking radius is evaluated and the update is realized only if a given critical number of previously cracked links is reached. Consequently, the area of high interest is updated frequently but not necessarily in all steps. This allows to account not only for crack branching but also for nucleation of new macroscopic cracks in different locations. This procedure frequently requires searching for particles within a given radius, which can be time- consuming for large sets of particles. Therefore, it is convenient to store particles in an octree structure that enables fast spatial searching. 4.2. Interpolation mesh update If the area of high interest has changed, it seems natu- ral to change the geometry of the interpolation mesh as well to provide a better shape of interpolation ele- ments at the interface between areas of low and high interest. However, this change is extremely expensive, as explained in section 4.3, and does not bring any sig- nificant improvement of accuracy. Therefore, it turns out to be more effective to use relatively small ele- ments with fixed geometry instead of changing the size of elements during the simulation. Only the elements that are completely covered by the extended area of high interest become useless and can be removed from the model. If the interpolation elements are not changed, the updating radius must be larger than the size of in- terpolation elements to ensure that the nodes of one interpolation element are not located on the opposite sides of the crack (otherwise, special finite element techniques such as XFEM must be used). 4.3. Stiffness tensors update Generally, there are three reasons why stiffness tensors must be updated. (1.) Breakage of a link (2.) Change of a region of high interest (3.) Change of the interpolation mesh The breakage of one link occurs in each simulation step. The resulting update of stiffness tensors can be simply realized by subtracting the stiffness contributions of broken links from all corresponding stiffness tensors. The change of the area of high interest is realized by adding new repnodes. Consequently, the type of ending particles of some links is changed. All links changed to type repnode–repnode and some specific links changed to type repnode–hangingnode need to be solved explicitly from this moment. All such links must be identified and their stiffness influence is sub- tracted in the same way as for broken links. The change of the geometry of interpolation mesh re- quires a recalculation of local stiffness tensors. This re- calculation is extremely demanding because all links must be visited. Even if the mesh is changed only locally near a newly created area of high interest and only some elements are affected, it is not easy to identify which links contribute to these elements and this step significantly slows down the simulation time. The recalculation of stiffness can be sped up by storing the information of element contributions during the assemblage of stiffness tensors at the beginning of the simulation. But this can be done only at the price of additional memory demands, which may become critical for large numbers of links. 5. Results To assess the efficiency and accuracy, the presented adaptive QC algorithm has been used to simulate failure in two dimensions. The L-shaped specimen with dimensions 100 × 100 mm (which may correspond to a joint of a frame) is fixed at the bottom section and loaded by prescribed vertical displacements im- posed at the right end section; see Figure 1. As a result of this loading, the non-convex corner is opened 82 vol. 13/2017 Adaptive Quasicontinuum Simulation of Disordered Lattices 0 20 40 60 80 100 0 10 20 30 40 50 60 70 80 90 100 ↑ ↑ ↑ ↑ ↑ ↑ 0 20 40 60 80 100 0 10 20 30 40 50 60 70 80 90 100 ↑ ↑ ↑ ↑ ↑ ↑ 0 20 40 60 80 100 0 10 20 30 40 50 60 70 80 90 100 ↑ ↑ ↑ ↑ ↑ ↑ Figure 1. Interpolation elements and area of high interest (black) and broken links (red) in the first (left), 100-th (middle) and the last step of simulation (right). 0 20 40 60 80 100 0 10 20 30 40 50 60 70 80 90 100 Figure 2. Broken links computed with exact particle approach A0 (black) and with QC approach with local anisotropic homogenization A3 (red). and cracking of links is expected to start in this re- gion. Therefore, a small initial area of high interest is prescribed around this corner; see Figure 1 (left). The random microstructure is generated with a den- sity that corresponds to 67 particles along the short edge, leading to a total of 19,000 particles connected by 64,064 links. Basic characteristics of models used by full particle and QC simplified simulations are listed in Table 1. Material parameters are considered to be the same for all links. As expected, cracking starts in the non-convex cor- ner. Duringcrack growth, the area of high interest is increased around newly broken links. Interpolation elements completely covered by the area of high inter- est are removed from the model, but the geometry of the remaining interpolation elements is not changed. The areas of high interest in different loading steps are depicted in Figure 1. Overall 240 simulation steps 0 1 2 3 4 0 5 10 15 20 25 30 35 40 45 Displacement [mm] F o rc e [ N ] Exact particle approach A0 QC simplified approach A3 Figure 3. Force-displacement diagrams computed with exact particle approach A0 (black) and with QC approach with local anisotropic homogenization A3 (red). A0 A1 A3 Particles 19,000 19,000 2,199 Links 64,064 64,064 2,929 Elements 0 1,542 1,542 Repnodes 19,000 1,665 1,665 Hnaging nodes 0 17,335 534 DOFs 37,799 3,129 3,129 Table 1. Numbers of particles, links, elements, repn- odes, hanging nodes and unknown DOFs for exact approach (A0), QC approaches with interpolation only (A1) and QC approach with local anisotropic homogenization (A3). are performed until the specimen is almost completely split; see Figure 1 (right). Three different levels of simplification are used and compared. • Full (exact) particle approach (A0). • QC simplified approach with interpolation only (A1). 83 Karel Mikeš, Milan Jirásek Acta Polytechnica CTU Proceedings • QC simplified approach with interpolation and local anisotropic homogenization (A3). The results obtained for both QC simplified ap- proaches are almost equivalent. Both approaches pre- dict the same broken links in the same order. Gen- erally, approaches involving homogenization tend to a stiffer response in comparison with the interpola- tion approach. In this case, the influence of homoge- nization is negligible and the maximum difference in relative stiffness is just 0.24 %. In comparison with the exact approach, the QC approaches give a slightly stiffer response. The relative stiffness error of the initial elastic branch is 11.8 %. The maximum relative stiffness error is below 18 % during almost the whole simulation. Only in the last six steps (when the specimen is almost completely split) the error increases up to 25 %. The prediction of cracked links is not perfectly exact, but more than 90 % of links are broken correctly and the difference in the final shape of the macroscopic crack is almost negligible; see Figure 2. The simplified results are also very accurate in terms of the force-displacement diagram. Despite a stiffer elastic response, the trend of the unloading branch is very well captured and the shape of the QC simpli- fied force-displacement diagram resembles the exact solution; see Figure 3. Times consumed by individual parts of the simu- lation for various approaches are listed in Table 2. The initialization procedures are realized just once at the beginning of the simulation. On the other hand, the repeated procedures are realized in each simula- tion step. For QC approach A1, with interpolation only, the realization of QC simplification is relatively fast, but the initial assemblage of the stiffness matrix and load vector is extremely slow because of a huge number of hanging nodes with interpolated DOFs. By contrast, QC approach A3, involving local homog- enization, needs a relatively long time to realize QC simplification including local effective stiffness tensors calculation. However, the initial assemblage is very fast because the number of links and hanging nodes is significantly reduced. For both QC approaches, a long extra time is needed only at the beginning of the simulation. This initial time becomes negligible in simulations with a large number of repeated steps. The most important fact is that each simulation step of the simplified approach is twice as fast for A1 and almost 19 times faster for A3 in comparison with the full particle model. 6. Conclusions This paper described a QC adaptive algorithm with automatic changes of area of high interest designed for simulation of crack propagation in disordered lattices. The presented algorithm allows not only propagation of an existing crack but also initialization of new cracks. A0 A1 A3 Times of initial procedures [s]: QC simplification - 0.83 8.08 Assemble stiffness matrix 1.09 5.81 0.15 Assemble load vector 1.92 6.36 0.21 Times of repeated procedures [s]: Solve one step 3.99 0.25 0.15 Update one step 1.08 2.31 0.12 Table 2. Time consumed by individual parts of sim- ulation for various approaches. The presented example has shown that QC based approaches in combination with an adaptive algorithm are able to capture the exact macroscopic crack trajec- tory even if the area of crack propagation is not known at the beginning of simulation. The force-displacement diagram can be predicted with high accuracy in spite of the fact that the initial elastic response of simplified models is slightly stiffer. The adaptive extension of the area of high interest has turned out to be sufficient to obtain accurate results and a new triangulation of interpolation mesh is not necessary. The application of QC approaches provides a sig- nificant reduction of the number of unknown DOFs and of links used in numerical models. This reduction leads to a significant speed-up of one simulation step, at the price of an increased cost of the initial prepa- ration (model simplification) before the first step. If both interpolation and homogenization are applied, the speed-up factor of one step is 18.8 and the whole simulation of crack propagation (including 240 steps) gets more than 15 times faster. Acknowledgements Financial support received from the Czech Technical Uni- versity (SGS project 17/043/OHK1/1T/11) and the Czech Science Foundation (GAČR projects 14-00420S and 17- 04150J) is gratefully acknowledged. References [1] E. B. Tadmor, M. Ortiz, R. Phillips. Quasicontinuum analysis of defects in solids. Philosophical magazine A 73(6):1529–1563, 1996. [2] D. M. Kochmann, J. S. Amelang. The quasicontinuum method: Theory and applications. In Multiscale Materials Modeling for Nanomechanics, chap. The quasicontinuum method: Theory and applications, pp. 159–193. Springer, 2016. [3] K. Mikeš, M. Jirásek. Quasicontinuum method extended to irregular lattices. Computers and Structures 192:50–70, 2017. [4] R. E. Miller, 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(5):053001, 2009. 84