Making use of symmetries in the three-dimensional elastic inverse homogenization problem

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/330132331 Making use of Symmetries in the Three-Dimensional Elastic Inverse Homogenization Problem ArticleinInternational Journal for Multiscale Computational Engineering · January 2019 DOI: 10.1615/IntJMultCompEng.2019029111 CITATIONS 2 READS 1,251 5 authors, including: Some of the authors of this publication are also working on these related projects: Fractional quantum Hall effect View project Material design - Metamaterials View project Carlos G. Mendez Universidad Nacional del Litoral 22 PUBLICATIONS105 CITATIONS SEE PROFILE Juan Manuel Podesta National University of the Northeast 7 PUBLICATIONS25 CITATIONS SEE PROFILE Sebastian Toro Universidad Nacional del Litoral 16 PUBLICATIONS131 CITATIONS SEE PROFILE Alfredo E. Huespe Universidad Nacional del Litoral 111 PUBLICATIONS2,347 CITATIONS SEE PROFILE All content following this page was uploaded by Alfredo E. Huespe on 06 September 2019. The user has requested enhancement of the downloaded file. Making use of Symmetries in the Three-Dimensional Elastic Inverse Homogenization Problem C. M´endez1, J.M. Podest´a1, S. Toro1, A.E. Huespe1,2∗ , J. Oliver2,3 1CIMEC-UNL-CONICET, Predio Conicet Dr Alberto Cassano, CP 3000 Santa Fe, Argentina 2Centre Internacional de Metodes Numerics en Enyinyeria (CIMNE),Campus Nord UPC. 3E.T.S dEnginyers de Camins, Canals i Ports, Technical University of Catalonia (Barcelona Tech) Campus Nord UPC, M`odul C-1, c/ Jordi Girona 1-3, 08034, Barcelona, Spain Abstract The objective of this paper is the design of three-dimensional elastic metamaterials with periodic microarchitectures. The microarchitectures of these materials are attained by following an inverse design technique jointly with an homogenization-based topology optimization algorithm. In this context, we have particularly studied the connection between the symmetry of the material layout at the microscale of 3D periodic composites and the symmetry of the effective elastic properties. We have analyzed some possible Bravais lattices and space groups, which are typically associated with crystallography, to study the way in which the symmetries of these geometrical objects can be usefully used for the microarchitecture design of 3D elastic metamaterial. Following a previous work of the authors for two-dimensional problems, we suggest adopting the design domain of the topology optimization problem coincident with the Wigner-Seitz cells of specific Bravais lattices having the same point group to that of the target elasticity tensor. The numerical assessment described in this papers aims at the design of an extreme material. The solutions obtained with this procedure show that different composite microarchitectures emerge depending on the cell shape selection. Keywords: elastic symmetry; three-dimensional homogenization-based topology optimization; Wigner-Seitz 3D cells, synthesis of elastic microtructures. ∗Corresponding author. E-mail address: ahuesp[email protected] (A.E. Huespe). 1 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) 1 Introduction 1.1 Objective and problem motivation In this paper, we present a methodology aiming at the microarchitecture synthesis of three-dimensional periodic composites whose homogenized elastic properties are similar to stipulated effective elasticities. The design of composites satisfying this requirement is a well-known problem in the literature. It has been clearly explained and described in the book by [1] and references cited therein. A number of contributions addressed to solve this problem have been posteriorly published which should be mentioned in the following. In particular, typical examples where this problem naturally arise, within the context of a larger design problem, can be seen in the papers by [2] and [3]. To reach this objective, we appeal to the symmetries characterizing the material layouts of composites at the microstructural level as well as to the symmetry characterizing the stipulated effective elasticities. The symmetry notion is ubiquitous in the nature. In general, this notion results in an extremely important physical property. Particularly, some features of material responses can be predicted without resorting to experimentation or complicated calculations by only appealing to the geometrical symmetries at different length scales of the material structure. 1.2 Microarchitecture design using symmetric topologies The microarchitecture synthesis of a periodic composite with prescribed target effective elastic properties can be though as an inverse design problem that is formulated by means of an homogenization-based topology optimization problem, such as proposed in the pioneering work by [4] and widely reported in the posterior literature ([5], [6], [7], [8], [9]). Nice microstructures have been obtained by [10] with this formulation and more recently by [11] in the context of three-dimensional microstructure designs. In the present paper, we follow this approach aiming at the three-dimensional microarchitecture design of elastic metamaterials realized as biphasic composites. A particular issue arising in this type of inverse homogenization-based design approach is the selection of the domain shape where the topology optimization problem is solved. In our opinion, this domain shape would result in a good selection if, as an outcome of the approach, the same domain is also a unit cell of the designed composite. There are some aspects to be considered for selecting adequate cell shapes. For example, [7] reported that the cell size influences the obtained topology solutions. But, which is more important for the approach taken in the present contribution is 2 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) the fact that a significant number of material layout topologies could be hidden or unreachable for some frequently used cell shapes, typically cubes or right rectangular prisms in 3D ([12]). Furthermore, It has been shown that the realization of certain classes of composites, such as the Vigdergauz microstructures or the microstructures proposed by [10], could be promoted by enforcing some kind of material layout symmetry. These issues have been particularly studied in the 2D elastic material design context by the authors in previous contributions, see [3] and [13]. In view of these observations and following similar arguments to the ones given by [13], here, we propose to employ concepts taken from crystallography to define the cell shape. These concepts are intimately related to the symmetry properties of the crystal structures and the elastic target tensors. Thus, three-dimensional crystal symmetry properties, defined according to their point groups and space groups, are the primary information guiding the selection of the unit cell shape that is considered for solving the topology optimization algorithm. Our proposal consists of taking the Wigner-Seits cell of a Bravais lattice that is compatible with the point group of the target elasticity tensor. Also, the material layout at the microcell could be enforced to satisfy the space group symmetry compatible with the same elasticity tensor. By adopting this approach, it can be guarantee that the homogenized elastic properties of the designed composite will have the same, or higher, symmetry than the target ones. This methodology is a generalization to three-dimensional problems of the one reported in [13] for twodimensional cases. By using this approach in the three-dimensional problems, we notice a challenge that does not exist in the two-dimensional ones; there is not any crystal system space group guaranteeing the realization of a periodic composite with isotropic effective elastic properties. Thus, we test our proposal by designing the microarchitecture of an isotropic extreme material with different Wigner-Seitz cells. In this test, the topology optimization algorithm does not explicitly impose the isotropy constraint. The obtained microstructure results are evaluated with a criterion measuring the proximity to an effective isotropic response. A brief description of the paper is next given. Section 2 presents a brief summary of the microarchitecture design methodology that we follow. This methodology can be seen as an inverse homogenization technique formulated as a topology optimization problem. Section 3 and Appendix A give a short description of crystal symmetries and their connection with the symmetries of the material elastic properties. Evidently, a full explanation regarding the crystal structure symmetries largely exceed the scope of the present paper; thus, the interested reader should consult the specific literature for additional information on this topic. Section 4 describes the main contribution of the paper. We explain the proposed procedure guiding the selection of the cell shape and the possible space group symmetry to be imposed in the microarchitecture topology design. 3 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) A numerical assessment making use of this procedure is presented in Section 5. We show a three-dimensionl design problem of an extreme isotropic composite material with maximum shear and bulk modulus. We seek the solution by only testing a cubic crystal system and its three related Bravais lattices. Finally, Section 6 presents the conclusions and final remarks of this paper. The topology optimization algorithm is briefly described in Appendix B. 2 Homogenization-based topology optimization problem Material design via inverse homogenization refers to the problem of finding the material configuration at the microscale of a periodic composite whose effective elasticity tensor is identical to a target elasticity tensor. This problem has been formulated as an homogenization-based topology optimization problem in the pioneering works of [4], [14] and [10]. In the following years, several works using this approach and emplying a variety of topology optimization algorithms have been reported. Only to mention a few, we cite the works by [5], [11] who describe threedimensionl elastic metamaterial designs using the SIMP optimization algorithm, and [15], who have used the bi-directional evolutionary structural optimization (BESO) algorithm. Two interesting recent reviews about this methodology can be found in [16] and [17]). Here, we follow a similar homogenization-based approach. Let us consider a structure whose material is a periodic composite constituted by two isotropic elastic phases M1and M2. We take a unit cell of this material identified by Ωµ. In this microcell, phases M1and M2occupy the domains Ω1 µand Ω2 µ, respectively, see Figure 1. The characteristic function χ(y) is defined in Ωµ. It identifies the positions where the phase M1is placed and takes the following values: χ(y) = 0∀y∈Ω2 µ 1∀y∈Ω1 µ .(1) Evidently, the homogenized elasticity tensor of the composite, Ch, depends on the geometrical configuration of the phases M1and M2in Ωµ. This dependence is made explicit by introducing the notation Ch(χ). This tensor can be evaluated in Ωµby enforcing periodic boundary conditions in displacements fluctuations. Then, standard computational techniques based on finite elements ([18]) or Fast Fourier Transform ([19]) can be used to get this goal. 4 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) Micro-cell y M1 M2 Wm Wm 2 Wm 1 M2 c=1 c=0 Figure 1: Microcell used as the design domain for the Topology Optimization Problem. The material layout is defined in Ωµ. The homogenized elastic properties Chof the composites are also computed in Ωµ 2.1 The topology optimization algorithm Next, we formulate the microarchitecture inverse design problem as a topology optimization problem. Let be given the design domain Ωµand the target effective elasticity tensor ˆ C. Let also be given the space Vχcollecting together all the characteristic functions χ(y) in Ωµwhere the domains Ω1 µand Ω2 µcan be arbitrarily changed. Then, the optimization problem is formulated as follows: min χ∈Vχ kCh(χ)−ˆ Ck such that: Vh(χ)−Vobj = 0 . (2) Where Vh= (RΩµχ dΩ)/|Ωµ|, is the volume fraction of the stiff phase M1and Vobj is the target volume fraction of the same phase. The problem (2) can be reformulated and implemented by defining a continuous function ψin Ωµwhose zero level set (ψ= 0) identifies the stiff-void material interface. The position of the zero level set is iteratively updated, using a gradientlike method. The update direction for ψis defined through the topological derivative of Ch. Hence, the problem (2) is rephrased as the optimization problem (13), in Appendix B. The algorithm for solving it has been reported by Amstutz and coworkers, see [20] and [21], and a brief summary is presented in Appendix B. To facilitate the comparative analysis performed in Section 5, we use a conventional Helmholtz filter forcing the optimization algorithm in Appendix B to provide solutions displaying approximately only one length scale of the stiff phase. This filter is implemented using equation (25), Appendix B. Additional details of this filter can be found in [22] and [23]. 2.2 Issues related to the design domain selection Two implicit variables in the problem (13) of Appendix B must be defined in advance. Their values result from a decision taken by the designer and play an im5 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) portant role to govern the complexity of the attained microarchitecture topology. They are: i) the shape of the design domain Ωµ. ii) the microstructure periodicity directions. In fact, to find the effective properties of the composite, it is necessary to impose constraints on the displacement fluctuations on the boundary of Ωµwhich are compatible with the characteristic periodicity directions of the composite. The periodic constraints in the inverse homogenization problem result from a decision taken by the designer, by arguing that the microarchitecture is periodic along pre-established directions. Due to this assumption, the full material architecture is a spatial replica, by tessellation, of the cell Ωµ. 3 Crystal symmetry properties We start this Section by briefly discussing some crystal symmetry properties. In the present context of the microarchitecture design, crystal symmetries are exploited to guide the designer in making adequate decisions about the design domain selection. Crystals can be characterized according to their specific symmetry properties. These properties are inherited from the underlying Bravais lattices and from their motifs. Two crystals sharing identical point group symmetry elements are said to belong to the same crystal system. Implicit to this classification is considered the translational and glide symmetries of the crystal motif which are taken into account to characterize the crystal symmetry properties, they constitute the space group of the crystal system. We postpone until Appendix A the discussion of further details about well-established crystal symmetry properties and the description of the notation adopted to identify point and space groups which are used in this Section1. The characteristic feature referred to the crystal structure that is here stressed is the relationship between the crystal symmetry of a given material and the symmetry of its associated effective physical properties. In particular, the elastic symmetry properties. This relationship can be stated in terms of the Neumann’s principle which establishes that: the symmetry elements of any physical property of a crystal must include the symmetry elements of the point group of the crystal. (see [26], pp.20). According to this principle, the relation between crystal systems and elastic symmetry classes is given by their compatibility with similar point group transformations. This fact is summarized in Table 1 displaying six columns that are explained 1For a detailed mathematical and physical description of the following concepts: crystal system, Bravais lattice, lattice system, point group, space group, etc., the reader may consult [24] and [25]. 6 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) in the next sub-Sections. Each row, in the first column, describes an elasticity tensor class which is associated with a specific crystal system and to more than one possible point groups and space groups. The fifth column indicates the number of Bravais lattice types that share identical point group symmetries with the elasticity class shown in column 1. They are collected together as a lattice system. These Bravais lattices are depicted in Table 5, Appendix A. 3.1 Elastic symmetry classes and crystal systems The first column in Table 1 displays the structure of the matrix representation2of the elasticity tensor Cand the identities that its coefficients have to satisfy accordingly to the symmetry of the elasticity class to which it belongs to. The established relations between the first and third columns in Table 1 are determined with the procedure proposed by [27]. This procedure uses the following sequence of operations. First, a point group is taken. Subsequently, the symmetry conditions being compatible with the point group elements are enforced to the elastic response, and the resulting format of the elasticity matrix is determined. Following this procedure with the 32 point groups, it is possible to collect together the set of point groups rendering identical matrix formats. Thus, it can be concluded that every matrix format is compatible with a specific set of point groups. The crystal symmetries can be classified by using a similar procedure, but now the cross relationships, i.e. finding common symmetry elements, are performed with the space group transformations. Using this technique, only seven crystal systems are determined. Furthermore, as a consequence of the Neumann’s principle, crystal systems and elasticity classes are related by sharing similar point groups. This association is evidenced in Table 1 between the first and second column. Therefore, an elastic class is designed with the name given to the corresponding crystal system. An exception to this rule is the isotropic elastic class which has not associated any crystal system. The coefficients of the elastic matrices are described in a coordinate system whose coordinate planes coincide with the crystal symmetry planes. We call it the natural coordinate system. Additionally, by appealing to thermodynamic stability arguments, these matrices should be positive definite. 2Kelvin’s notation is assumed for the matrix form of the elasticity tensor 7 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) Table 1: Symmetry elements (point and space groups) of Crystal Systems (CS), Elasticity Classes (EC) and Lattice Systems (LS). Each row identifies the triplet (CS-EC-LS) with point group compatibility. First column describes the matrix format of the elasticity tensor class. The symbols in column 1 indicates: “∗” a possibly non-zero component; “∗−∗” two identical components; “∗ − ¯ ∗” two equal components with opposite sign; “⊗” the shear moduli is equal to (C11 −C12) (see Ting [27]). Fifth columns indicate the number of Bravais lattices constituting the lattice system defined in the sixth column. Elasticity Crystal Point Space Nbr. of BraLattice Tensor System Groups Groups vais Lattices System Triclinic 1,¯ 1 2 1 Triclinic Monoclinic 2,m,2/m 13 2 Monoclinic Orthorhombic 222,mm2, mmm 59 4 Orthorhombic Tetragonal 4,¯ 4,4/m, 442,4mm, ¯ 42m,4/mmm 68 2 Tetragonal Trigonal 3,¯ 3, 32,3m, ¯ 3m 7 1 Rhombohedral 18 1 Hexagonal Hexagonal 6,¯ 6, 6/m,622, 6mm,¯ 6m2, 6/mmm 27 Cubic 23,m¯ 3, 432,¯ 43m, m¯ 3m 36 3 Cubic Total 7 32 230 14 7 8 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) (a) (b) (c) L U WX G K L U W X GK P GH NH N G P Figure 5: Finite element meshes of the Wigner-Seitz cells (BCC and FCC cells). a) Reduced domains. Γ is the central point of the Wigner-Seitz cell; Lthe central point of a hexagonal face; X the central point of the contiguous square face; K the middle point of the edge between two hexagonal faces; U the middle point of the edge between the hexagonal and square. a) Meshes of the reduced domains. b) Meshes of the full domains obtained from the meshes of the reduced domains by successive transformations of symmetry (reflections and rotations). This meshing technique provides a symmetric finite element mesh that preserves the higher point group symmetry of the corresponding Wigner-Seitz cell. Thus, using these symmetric meshes, different kind of material distributions, with predetermined space group symmetries, can be easily implemented in the topology optimization algorithm. using the so-generated symmetric meshes, different kind of material distributions, with predetermined space group symmetries, are much simpler to implement or impose in the topology optimization algorithm. The number of finite elements in the three cells are taken such that the volume of the elements is similar in the three cases. Thus, the number of finite element in the SC, BCC and FCC cells are proportional to 1, 0.5 and 0.25, respectively. Starting configuration of the topology optimization algorithm. Topology optimization problems aiming at microstructure design, in general, contain many local minima. This characteristic induces a strong tendency to attain different solutions, depending on the initial guess configuration ([17]). Consequently, it is important to test several starting configurations to evaluate and compare the so-obtained 15 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) solutions. In the present numerical assessment, three initial configurations are tested: a) a spherical void placed around the central point of the cell; b) a cellular-like configuration with closed walls. The walls are coincident with the cell faces and the walls are of uniform thickness; c) a truss-like configuration with bars of identical sections joining the Wigner-Seitz cell vertices. 5.2 Assessment of the cell capacity for capturing isotropic responses As already mentioned above, the isotropic elastic response of the designed composite cannot be guaranteed by only enforcing a stiff phase material layout consistent with the space group of some particular crystal system. This aspect of the threedimensionl problem introduces a marked difference with respect to the 2D case which has been analyzed in the previous contribution by the authors, see [13]. In consequence, we introduce an indicator to find how close are the homogenized properties of the designed composites to isotropic responses. Based on this indicator and without imposing implicitly the isotropy constraint into the formulation of the optimization algorithm, we assess which cell provides a better response to capture this elastic feature. This indicator is computed as follows. Given an arbitrary elasticity tensor Ch, we define an isotropic elastic tensor of comparison by following the procedure proposed by [30]. First, using the components of Ch, it is evaluated: Ciso 11 =1 5(Ch 11 +Ch 22 +Ch 33) + 2 15(Ch 12 +Ch 13 +Ch 23) + 4 15(Ch 44 +Ch 55 +Ch 66) (6) Giso =1 15(Ch 11 +Ch 22 +Ch 33)−(Ch 12 +Ch 13 +Ch 23) + 3(Ch 44 +Ch 55 +Ch 66) (7) Kiso =Ciso 11 −4 3Giso (8) and with the so-determined bulk and shear moduli, Kiso and Giso, an isotropic tensor Ciso is computed with expression (5). Our assumption is that this is the closer isotropic tensor to the original Ch. It is taken as the reference tensor to perform the following analysis. Finally, we introduce a coefficient of anisotropy χ=kCh−Cisok(9) which measures the distance between Chand Ciso. A zero value of this coefficient indicates that Chis isotropic. Contrarily, a large value of χindicates that Chis far from being isotropic. 16 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) 0.060 0.064 0.068 0.072 0.076 0.150 0.145 0.155 0.160 0.165 0.170 Effective Bulk Modulus Effective Shear Modulus Simple Cubic Body Centered Cubic Face Centered Cubic Truss Cellular Spherical void 0.00 0.01 0.02 0.03 0.04 0.05 0.06 BCC FCC SP Initial configuation 0.008 0.006 0.004 0.015 0.006 0.019 0.053 0.019 0.043 SP BCC FCC (a) (b) Anisotropy coefficient c Color reference (initial configuration) Symbol reference (lattice type) Lattice type Truss Cellular Spherical void Figure 6: a) Space Keff vs. Geff : the results obtained with Primitive (SP), Face Centered (FFC) and Body-Centered (BCC) cubic Wigner-Seitz cells are superposed to the Hashin-Strikman upper bounds displayed in green lines. The initial configurations adopted in each case are distinguished with blue, red and yellow colors; b) Anisotropy coefficients χfor the nine tested case. 5.3 Analysis of results Figure 6-a plots the results for the nine different cases in the space Keff vs. Geff . They are compared with the upper Hashin-Strikman bounds, also plotted in the same Figure. The results have been computed using the three Wigner-Seitz cells: Primitive Cubic (SP), Face Centered Cubic (FFC) and Body-Centered Cubic (BCC) and the three initial configurations. The symbols identify the cell types and the colors identify the initial configurations. The corresponding coefficients of anisotropy, defined by expression (9), are plotted in Figure 6-b. Independently of the adopted initial configurations, we note from this plot that the BCC cell provides microarchitecture topologies whose effective elastic properties tend to be more isotropic respect to the solutions provided by the other two alternative cells. Figure 7 depicts the material distributions obtained with these cells and with the spherical void initial configurations. Two views of each solutions are displayed; the full cell solutions are shown in Figures b), d) and f) and the cells cut with middle planes are shown in Figures a), c) and e). Note that the material distributions assimilate to hollow topologies. Also, note the highly symmetric pattern of the material layouts obtained with the algorithm even without imposing any symmetry constraint. Table 4 shows the components of the target effective elasticity tensor computed with expression (5). These components are compared with the homogenized tensor 17 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) x z y BCC (Spherical Void) x z y FCC (Spherical Void) SP (Spherical Void) (a) (c) (e) (b) (d) (f) Figure 7: Microarchitectures with maximun bulk and shear moduli: a) and b) correspond to BCC cells (cutted with a middle plane (a) and full cell solutions (b), respectively); c) and d) correspond to FCC cells (cutted with a middle plane (c) and full cell solutions (d), respectively); e) and f) correspond to SP cells (cutted with a middle plane (e) and full cell solutions (f), respectively). 18 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) Table 4: Computed effective elasticity tensors C11 C22 C33 C44 C55 C66 C12 C13 C23 Target 0.2750 0.2750 0.2750 0.1578 0.1578 0.1578 0.1172 0.1172 0.1172 SP Cellular 0.2732 0.2732 0.2732 0.1398 0.1394 0.1392 0.1177 0.1177 0.1179 BCC Cellular 0.2559 0.2571 0.2552 0.1272 0.1274 0.1276 0.1246 0.1243 0.1244 FCC cellular 0.2570 0.2572 0.2569 0.1286 0.1274 0.1284 0.1242 0.1239 0.1241 values obtained with the SP, BCC and FCC cells and cellular-like initial configurations. We observe from Figure 6-a that, apparently, the SP solution display closer values to the target ones. However, the identity ˆ C33 = ( ˆ C11 −ˆ C12), which must be satisfied by isotropic tensors, is more tightly verified with the BCC and FCC solutions if compared with the SP solution. This behavior, which is not transparent from results in Figure 6-a, is confirmed with those depicted in Figure 6-b. The stiff phase volume fractions in these solutions are accurate to the third figure (f= 0.338). 6 Conclusions New contributions for synthesizing three-dimensionl microarchitectures of elastic composites are proposed. These contributions are addressed to enlarge the range of attainable microarchitectures in the framework of the inverse design problems formulated as homogenization-based topology optimization algorithms Based on the symmetry of the target elastic properties, we propose a procedure for selecting a spatial three-dimensionl domain where the topology optimization algorithm is solved. Furthermore, this procedure also provides a route to enforce the symmetry of the material layout within this spatial domain compatible with predefined crystal spatial groups. The so-proposed rules are derived from concepts widely developed in crystallography. This contribution is a generalization of a procedure that has been previously presented by the authors in 2D problems. The procedure has been tested by synthesizing an isotropic elastic material with prescribed maximum effective shear and bulk moduli. The results, which have been obtained without implicitly imposing a constraint of overall isotropic response, show that the BCC cells provide the tighter isotropic solution if compared with alternative cubic cells, no matter the initial configuration adopted for the topology optimization algorithm. Furthermore, by considering that the volume ratio between the SC and 19 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) BCC cells is two; then, the size of the problem with the BCC cell is half of the problem size using the SC cell for similar resolutions. An additional advantage of the proposed methodology is that very simple and precise numerical techniques, widely reported in the literature, can be used to generate the three-dimensionl Wigner-Seitz cells. Algorithms for determining the WignerSeitz cells of arbitrary distributions of points in the space, particularly the Bravais lattice atoms, are well developed in Computational Mechanics and can be easily used to determine the Wigner-Seitz cells of the 14 Bravais lattice types. Additionally, the technique based on meshing a reduced domain, and subsequently generate the mesh of the full cell through space group symmetry operations, facilitates the implementation of different types of symmetry constraints on the material layout. Acknowledgment The authors acknowledge the financial support from CONICET and ANPCyT (grants PICT 2014-3372 and 2016-2673). Authors A.E. Huespe and J. Oliver acknowledge the funding received from the Spanish Ministry of Economy and Competitiveness through the research grant DPI2017-85521-P, Project “Computational design of Acoustic and Mechanical Metamaterials” (METAMAT). A APPENDIX I: Structure and Symmetry of Crystals A.1 Bravais lattices Crystals are periodic structures, with periodicity along three linearly independent directions. This property is formalized by introducing the concept of a lattice. Lattices are defined by three vectors a,band c, the primitive vectors of the lattice, which form a basis in R3. The set of abstract points, or atoms L:= {la+ma+nc|l, m, n ∈Z}(10) constitute a Bravais lattice. In the three-dimensional space, according to the relationship between the vectors of the basis, we can distinguish 14 types of Bravais lattices which are characterized in Table 5. Notice that they are collected in a set of lattice systems. A.1.1 Unit Cell The space R3can be subdivided into cells of finite volume having all the same shape. They are called unit cells if these volumes cover all the space without overlapping 20 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) Table 5: Catalog of the 14 Bravais lattices classified according to their lattice system Lattice System Point Group Primitive Base-Centered Body-Centered Face-Centered Triclinic ¯ 1 P¯ 1 Monoclinic 2/m P2/m C2/m Orthorhombic mmm Pmmm Cmmm Immm Fmmm Tetragonal 4/mmm P4/mmm I4/mmm Rhombohedral ¯ 3m R¯ 3m Hexagonal 6/mmm P6/mmm Cubic m¯ 3m Pm¯ 3mIm¯ 3mFm¯ 3m 21 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) when are translated by the vectors of the lattice L. The primitive unit cell is a standard cell that consists of the parallelepiped generated by the primitive vectors of the lattice (third column in Table 5). Although all the 14 Bravais lattices have a primitive unit cell constructed in this way, in many cases, it is more convenient to use a conventional cell in order to better visualize the structure of the lattices (columns 4 to 6 in Table 5). These conventional cells are not necessarily unit cell because they could contain more than one lattice point and consequently their volumes are larger than the volumes of the primitive unit cell. Another standard construction for the unit cells is the Wigner-Seitz cell. This cell consists of those points of R3that are closer to a given lattice point than to any other point of L(see Table 3 for cubic lattices). Two remarkable features of the Wigner-Seitz cells are noted: i) its volume is the minimum one that a unit cell can have since, by construction, this cell contains only one lattice atom; ii) it preserves the same point group symmetry of the corresponding Bravais lattice. Primitive unit cells, in general, do not have this property. A.2 Crystallographic point groups Point symmetries are symmetry operations that leave at least one point fixed. There are five types of point symmetries in 3D. They are the inversion (¯ 1), rotations (2, 3, 4 or 6), reflection (m), rotation-inversion (¯ 2 = m,¯ 3, ¯ 4 or ¯ 6) and the identity (1). Not all rotations are allowed because they must be compatible with the discrete translation symmetry of the crystals, that is the reason because they are called crystallographic point groups. After applying the crystallographic restriction theorem, a total of 32 point groups are all the possible symmetries that a given crystal in 3D can have. They are enumerated in the third column of Table 1. By “crystal” we mean a lattice with a base or motif. Not all the point groups are compatible with all lattices, but many of them are compatible with the same group of lattices. This fact allows a classification of point groups in crystal systems, which can be seen in the second column of Table 1. The lattices themselves (or crystals without a base) have their own point group, which can be seen in the second column of Table 5. Lattices with the same point group are grouped in lattice systems and are shown in the first column of the table. Unit cells are not unequivocally defined and they can have practically any symmetry. However, as we mentioned above, Wigner-Seitz cells preserve the same symmetry as the associated Bravais lattices, and therefore, they have the same point groups shown in Table 5. As with the lattice, we are considering here Wigner-Seitz cells without a motif in its interior. 22 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) A.3 Space groups On one hand, the discrete translation symmetry of the crystal restricts the infinite possible point groups to only 32, but, on the other hand, it adds two new kind of symmetries. One type corresponds to the glide planes (denoted by a,b,c,nor d), which consist of a reflection followed by a translation. The other type corresponds to the screw axis (e.g., 32, 41, 63), which represent a rotation followed by a translation. Combining these two symmetry types with the 32 point groups, 270 space groups can be classified. In the forth column of Table 1, we can see the amount of them compatible with each crystal and lattice system. Besides the letters and numbers corresponding to the symmetries, in the nomenclature of space groups also appear a capital letter that helps to identify the compatible Bravais lattice (e.g., Pfor principal, Ffor face centered, Ifor body centered). In Table 5, the space groups of the 14 Bravais lattices can be seen below the images, and in Table 2, the 36 space groups compatible with the 3 cubic lattices are enumerated. B APPENDIX II: Solving the topology optimization problem with Topological Derivative Algorithm We summarize in this Appendix the topology optimization algorithm that is used for solving the numerical test presented in Section 5. The algorithm is a well-established technique reported in the papers by [20] and by [21]. It is a level-set method (LSM) with sensitivity computed through the topological derivative, see [31] and [32]. This technique has been implemented in a 3D code using an Augmented-Lagrangian scheme reported in [33]. Let us introduce a smooth level set-funtion defined in the microcell Ωµ,ψ∈ C0(Ωµ), satisfying ψ(y) =    <0∀y∈Ω2 µ >0∀y∈Ω1 µ 0 in the interfaces ,(11) then, the characteristic functions χ(y) in Ωµ, given by expression (1), can be redefined as follows: χ(ψ) = 0∀ψ≤0 1∀ψ > 0.(12) 23 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,) and the problem (2) is rephrased as: min ψ∈C0kCh(ψ)−ˆ Ck such that: Vh(ψ)−Vobj = 0 . (13) By making use of an Augmented Lagrangian technique, see [33], the problem (13) is rewritten as follows: max λmin ψT(ψ, λ),(14) with: T(ψ, λ) = kCh(ψ)−ˆ Ck+λ(Vh−Vobj) + α 2Vh−Vobj2(15) where Vh= (RΩµχ(ψ)dΩ)/|Ωµ|is the volume fraction of the hard phase, λis the Lagrange multiplier and αis the penalty parameter of the augmented term. The algorithm for solving the problem (14) has a loop where αis hold fixed and λis modified iteratively. The minimum of Tis searched with a descent direction algorithm. For problem (14), the descent direction is estimated with the term DψT(ψ, λ) = (Ch−ˆ C) : DψCh kCh−ˆ Ck+λ1+α(Vh−Vobj)1(16) where DψChis the topological derivative of the homogenized elasticity tensor and is given by the expressions: DψCh ij =εi µ·P·εj µ;P=m1(m2(1⊗1)+2I) ; (17) where Pis the polarization tensor, εi µand εj µare the micro-strain solutions of the homogenization problem solved with the i-th and j-th canonical macro-strains (i, j = 1,2, ..., 6), the symbols 1and Irepresent the second and fourth order unit tensor, respectively, and the coefficient m1and m2are: m1=15µδµ(ν−1) 15µ(1 −ν)+2δµ(5ν−4); m2=δλ[15µλ(1 −ν)+2λδµ(5ν−4)] −2δµ(λδµ−5µνδλ) 5δµ[3µλ(1 −ν)−3µνδλ−λδµ(1 −2ν)] ; (18) with δλ=λ−λ0;δµ=µ−µ0and (λ;µ) being the Lam`e parameters of the base (stiff) material and (λ0;µ0) being the Lam`e parameters of the material introduced as a spherical perturbation. The Poisson ratio of the base material is ν. Additional description and properties of this tensor can be found in [32], where it is called the Elastic Moment Tensor (EMT). 24 International Journal for Multiscale Computational Engineering, 17(3):261–280 (2019), Special Issue: Computational Multi-scale Modeling and Design of New Engineering Materials" Guest Editors: M. Pietrzyk, T.Burczynski, X. Oliver, A. Huespe,)