Architectural effect on 3D elastic properties and anisotropy of cubic lattice structures

• A large range of cubic structures was generated by using a continuum topology modelling, connectivites and properties of database can be directly plot in 2D color surface maps. • Evolution of elastic mechanical properties by the variation of relative density obtained by homogenization procedure was fitted by power law. • Power laws parameters represented the influence of topologies on the observed property. They were determined for cubic elastic constants and anisotropy. • Power laws parameters combined to computed connectivities led to the determination of different mechanical be-


H I G H L I G H T S
• A large range of cubic structures was generated by using a continuum topology modelling, connectivites and properties of database can be directly plot in 2D color surface maps. • Evolution of elastic mechanical properties by the variation of relative density obtained by homogenization procedure was fitted by power law. • Power laws parameters represented the influence of topologies on the observed property. They were determined for cubic elastic constants and anisotropy. • Power laws parameters combined to computed connectivities led to the determination of different mechanical behavior families.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
With the development of additive manufacturing processes, it is now possible to produce scaffold structures with an increasing level of tailoring. By changing the architecture, it is possible to decrease the Materials and Design 182 (2019) 108059 mass of the part while preserving a large range of accessible stiffness properties [1][2][3][4][5]. These structures are used to fill a specific volume in parts having a complex shape produced by additive manufacturing (Fig. 1a) [3]. Such complex geometries usually appear with multiaxial loadings [6,7]. On a macroscopic scale, the stress field may be seen as a combination of local elementary loadings: for instance, the part may be under pure shear load in one location, pure compression in another, and some complex loadings in between [6]. This issue cannot be reduced to uniaxial elastic properties solely, but it requires a global 3D description. In Fig. 1a, the case of a hip implant is illustrated, with external forces on the femoral neck. This results in a complex 3D stress field in the stem, with possibly the need for several optimal architectures depending on the location. On the other hand, at the microscopic scale, it is well known that a second stress field emerges from the porous character of lattices [1,8,9]. The effective stress field inside the lattice beams ( Fig. 1c) is no more relevant for mechanical design of the macroscopic part, because the local maxima interfere with the understanding of the smooth variations of the macroscopic stress (Fig. 1a) [6,10]. Therefore, there is a need for an intermediate description of the elastic properties of this complex struts network (Fig. 1b) at a mesoscale. The optimal properties are selected at a macroscopic scale while considering the microscopic features of the lattice, thanks to this intermediate element. This is made possible by using homogenization procedures to represent the scaffold through its effective elastic properties. These procedures give also the possibility to evaluate the anisotropy of a given architecture.
Most of existing works are dedicated to the study of one specific, or few topologies, like the work on the octet-truss structure [8,12,15], or the gyroid structure [18]. However, only few studies proposed to deal with a large range of architectures. We should mention the results of Xu et al., Dong et al. or Maskery et al. [16,17] working with six or seven different topologies. One of the most extended studies is the work of Vigliotti and Pasini [13] considering 12 topologies. The paper of Favre et al. [2] used 66 topologies that are however limited by a 1D description of elastic behavior.
There are several homogenization strategies and most of them can be applied to the case of strut-based lattice. Different schemes lead to determination of elastic properties based on analytical formulations, as in the works of [8,9,[12][13][14][15], while some others use numerical methods [16][17][18][19][20]. Stevens proposed a straightforward and easily implementable technique [21,22]. This procedure was successfully used in other studies with the same kind of microstructures [16,18]. This method may be considered as an upper bound for the elastic constants determined due to the kinematic boundary conditions. This work proposes to apply the homogenization procedure was based on the work of Stevens on a large set of architectures based on the works of Favre. A specific focus was done on the comparative study of the architectures concerning their elastic properties and their anisotropy. In consistence with the previous parametric studies [2,10], a large panel of topologies was swept by massive computing to generate a database of readily available solutions. Within this database, it is possible to identify different behavior groups (shear friendly or tensile friendly behavior). The reverse use of this database is a powerful tool for mechanical design to select the optimal structure for local targeted anisotropic properties. Moreover, the group combining both shear and tensile friendly behaviors could be considered as near-isotropic. This generic group may be used for general replacement of bulk isotropic materials, where no preferential orientation of the properties is required, while preserving the capability to tailor the stiffness [16].
The goal of this paper is to systematically estimate and compare stiffness constants and anisotropy of a large panel of topologies, and for a large range of relative densities. The optimal properties for specific standard mechanical loadings are determined and discussed by numerical models and experiments. First, the method for homogenization and experimental validation is presented. Second part focuses on the connectivity of cubic lattice structures. Then, the results collected in term of elastic and anisotropic properties are plotted as a function of the architecture in color maps, and as a function of density. We assess the validity of the model on a specific experimental case. Finally, the results are discussed in relation to the connectivity, and main families of behaviors are identified.

Lattice structures generation
The method to generate a database of architectures with a continuous variation of two parameters was implemented from crystallographic rules [2]. Using this method allows to sweep a range of architectures in a continuous way, while taking the advantage of symmetries to decrease the complexity of the network geometrical description. Main idea is to use a set of symmetry operations associated to the m3m point group that combines different space groups associated to the conventional cubic Bravais lattices (BCC (body centered cubic), FCC (face centered cubic) as well as the PC (primitive cubic)). These three Bravais lattices, all have a well-known lattice structure equivalent: BCC was related to the diag-structure or body centered [23], FCC to the well-known octet-truss of Desphande [8,12] and PC was associated to a simple cube [6]. One of the main advantages of this procedure is to enable a complete description of topologies by using only two parameters. This procedure allows to plot any kind of results into a 2D color map, where the 2D plane axis represent the (x, y) parametric space (see Fig. 1) and color stands for the variation of the different studied characteristics (connectivity, stiffness constant …).
All these parameters represent the information necessary to place the first beam in the space, other ones were further generated automatically by symmetry systems. Fig. 2 summarizes the procedure for a Python script generating lattices by inflating a 1D wireframe to a 3D volumic mesh, later imported to the Abaqus software. The two parameters mentioned earlier correspond to x and y coordinates of point B (Fig. 2a), while the z coordinate was maintained at mid of the cell size to preserve the continuity of the unit cell. This procedure to generate cubic lattices structures is well discussed and largely explained in a previous work, and the reader willing to further learn about this technique can refer to [2]. Each cell had a global size of 1x1x1 mm, and the radius of beams is defined as a fraction of this length.
In this study, the representative 2D space of the possible (x, y) couples was discretized by 36 topologies, each of which was produced with 4 different radii, in order to explore the 3D space of (x, y, r) for 144 different cases. All of these structures were meshed using Tetgen software [24]. The elements were set with a maximal face area of 0.005 mm 2 and a maximal volume of 0.01 mm 3 . It resulted in a mean number of elements per section between 12 and 20. These parameters were selected after a series of tests, in order to determine the convergence of the mesh for some topologies. Tetgen file format was converted to INP file format by a home-made Python script. Relative density was computed by adding up the volume of each tetrahedron on the 3D meshes. For this study the calculations were run on a single unit-cell. It is common to perform FE computation on a single cell as some publications show that for strongly periodic microstructure, a representative elementary volume limited to one unit-cell was acceptable [25,26].

Homogenization procedure
A straightforward and easily implementable homogenization routine, proposed by Steven, was used [21,22]. Stiffness matrix C ijkl is a fourth-order matrix, linking the macroscopic strain tensor ε kl to the macroscopic stress tensor σ ij in the Hooke's law σ ij = C ijkl ε kl . Due to the symmetry of σ ij , ε kl and C ijkl , C ijkl can be reduced from 81 to 21 constants. Six different boundary conditions are needed for the identification of constants. For each condition, one component of the strain tensor is equal to 1, while the five others are null. It results in a stress tensor which is obtained from the reaction forces on appropriated faces. This stress tensor provides directly the quantification of a row of C ijkl matrix, as indicated in Eq. (1): Boundary conditions are listed in Table 1. These conditions consist in three normal uniaxial strains ε x , ε y , ε z , and three shear strains γ xy , γ xz , γ yz , which respectively correspond to ε 11 , ε 22 , ε 33 , ε 12 , ε 13 , ε 23 . It was shown in other studies that this procedure gives well accurate prediction results for the macroscopic mechanical properties for lattice structures [16,18], but also for other kind of materials such as composite materials [22,27]. The homogenization procedure was implemented in Python language, and finite element analysis were performed in Abaqus 6.14, in a fully automatized procedure from the generation of structures to the post-treatment of FE analysis. For this study, an isotropic bulk material was used with a Young Modulus of 1 GPa and a Poisson coefficient set at 0.4. In a previous paper, authors have shown that the variations of stiffness with the relative density was material independent if a suitable normalization is applied [2]. The idea is to use this kind of procedure to explore the dependence of elastic constants with density. This dependence is analyzed in the Results section.
Once the stiffness matrix is resolved, it becomes possible to quantify the anisotropy of the lattice. The Zener ratio [28,29] defined in Eq. (2) was used to estimate the anisotropic behavior: This ratio was computed for each combination of topology and radius values to determine its dependence with the relative density. If this ratio is near to unity [16,29], then the structure can be seen as isotropic. It will be shown in the next section how this parameter can help to predict the ability of a topology to deform under tensile, shear, or both loading conditions.

Experimental validation
The predicted stiffness and anisotropy were assessed by experimental mechanical tests. The area of investigation was restricted to the study of one topology, to preserve a reasonable experimental time. The idea was mainly to confirm the predictability of the numerical models, and the systematic experimental test of all the topologies is not required. Therefore, the experimental Young and shear modulus were determined by uniaxial compression tests, and torsion tests on the architecture (0.25, 0.25, 0.5). Experimental measurement was compared to numerical results in the Result section "Experimental validation".
All samples were built by using selective laser melting technology (SLM) with a SLM 280 HL system (SLM Solutions). The powder material was Ti-6Al-4 V with spherical morphology and a particle size distribution centered on 45 μm. The part was built under Argon protective atmosphere with a laser of 200 W power, a scanning rate of 1650 mm·s −1 , and a hatching distance of 80 μm. In this study, we choose to fix the set of parameters to reduce the variability of the produced samples. These parameters correspond to the recommended values provided by the SLM manufacturer. Relative density was determined by relative weighting of porous sample as compared to the bulk materials of the same global dimensions.
Uniaxial quasi-static compressive tests were carried on a Zwick-Roell machine, with a load cell of 100KN and a strain rate of 1.10 −3 s −1 . Displacement measurements were carried out using an extensometer multiXtens. Samples were designed to fit with the recommendation of the norm ISO 13314 ( Fig. 3.b). They consisted of 5 repetitions of a unit-cell in the 3 directions of an orthogonal basis. The Unit-cell corresponded to a cube with a width of 3 mm. Three levels of beam radii (150, 200 and 225 μm) and bulk material were tested. It Table 1 Boundary conditions for the six different loading case.

Active strain component
Boundary condition leads to three theoretical relative densities of 21.4%, 37.0%, 44.4%, and 100%. For each density value, test was repeated four times in order to determine the incertitude of the measured Young modulus. Stiffness was determined as described by the norm ISO 13314. For the sake of simplicity, the determined quasi-elastic gradient is assimilated to the effective Young modulus.
Torsion quasi-static tests were carried out on a Zwick-Roell hydraulic machine, with a maximal couple of 250 N·m −1 at an angular rate of 0.5°·s −1 . The central area of the sample was designed in accordance with the ISO 13314 norm: it consists of a circular section of 15 mm with a height of 18 mm where unit-cells have the same dimension as for tensile testing. Two bulk cylindrical parts are added in the extremities for fixation in the jaws. For the purpose of mismatch area reduction between porous and bulk areas, a density gradient was used. It also increases the probability of failure in the central area. It consists of two slices of 3 mm of height where the struts radius is progressively increased compared to the bulk area (see Fig. 3.d). The same relative densities as for tensile tests were used. For each density value, test was repeated three times for the determination of shear moduli and its associated uncertainty values.
For both tensile and torsion tests, the uncertainty values are considered to be twice the value of the standard deviation around the average of Young modulus and shear modulus.

Results and discussion
The elastic properties can be deduced from specific effects of the connectivity, the density and the spatial organization of the lattices. The results are analyzed to separate methodically these effects on the elastic constants and on the anisotropy. First, the connectivity of the lattices was determined, as it is an important factor to understand the variations of stiffness [30,31]. Then, the dependence of elastic constants on density is analyzed for different architectures. The effect of the topology on this dependence is identified. The contribution of architecture and connectivity to the anisotropy of the lattice structures is determined. Then, the topologies are classified in groups depending on their deformation modes. Finally, the last section is an experimental validation of a particular case of interest, including discussion of divergence between theory and experiment.

Connectivity of cubic lattice structure
Connectivity in the crystallography field concerns the number of struts that are connected to the same node [30]. Since the work of Maxwell [32], Calladine [33], as well as Desphande [30], connectivity of a lattice is a direct indicator of framework rigidity. For a 3D framework, the necessary and sufficient condition to obtain a rigid structure is Z = 12 (where Z is the connectivity). When Z N 12 the system is redundant in term of rigidity. Whereas for Z b 12, the framework is underconstraint and presents a bending dominant behavior [30]. These works provide main guidelines to design of a rigid framework, but they suffers from some limitations. The first limitation was the notion of pin-jointed framework, which diverges from real situation of additively manufactured lattice structures having rigid joints [30,31]. The other limitation concerns the notion of similarly situated nodes, meaning that all nodes in the frameworks should have the same connectivity [30]. Actually, this is not the case, and a large number of lattices have several kind of nodes connectivity, like some auxetic structures [34] or the rhombic dodecahedron unit-cell [35]. Some examples of this kind of topologies are shown in Fig. 4.b. In this study, connectivity was calculated by the weighted average of each node connectivity. The weight of each connectivity is computed from the fraction of nodes, but also considers the division of nodes located at the boundaries (on the edges, vertices or faces).
One important point to note is that for any structure, there is always one node at the corner of the cell (see Figs. 1 and 2.b) and at least one secondary node present on faces, edges or inside the structure, as readers can see in Fig. 2

.b.
A 2D color map in Fig. 4a illustrates the calculated average connectivity. First, most of the topologies seem to have a very low connectivity (dark blue area in Fig. 4a). In specific areas of this map, such as the diagonal x + y = 0.5, the connectivity jumps to larger values. Some of these topologies fulfil the rigidity criterion Z = 12 [30]. Table 2 summarizes the essential data from the connectivity map (Fig. 4b). Connectivity Z 1 corresponds to the outer nodes shared with neighboring cells, and comes from point A in Fig. 2. The connectivity Z 2 is the one of inner nodes belonging to a single cell, and comes from point B in Fig. 2. The connectivity Z is the weighted average of Z 1 and Z 2 .
The category A corresponds to the structures with just one kind of node connectivity. It corresponds respectively to primitive cubic (just ConnecƟvity Y parameter X parameter Fig. 4. Connectivity as a function of x and y parameters (a) and locally generated topology for the different (x, y) couple (b) areas correspond respectively to structure groups in Table 2 (red: similarly situated, green: faces, blue: cross-shape, purple: unconnected-diag, yellow: connected diag, copper: other, dark purple dash line: strut length = ffiffiffi 2 p =2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) preserving edges of a cube as strut), octet-truss (or face centered cubic) and diag-structure (or cubic centered). The categories B and C correspond to the map boundaries in Fig. 2.b. They both present a connectivity similar to the A-type cubic primitive structure due to the sharing of lower connectivity inner nodes, but they have an outer node with a higher connectivity. The category C of "Cross-shaped" structure presents a larger connectivity due to the lower number of second order connectivity nodes (only 6, one per face). This set of structures also give a continuous way to switch from octet-truss to diag-structure. The areas D and F correspond to the so-called hexatruss structures, with a single V-shaped bond between the corners of the cube for D category, and two V-shaped bonds for F category. Group D along the diagonal of the (x,y) space (where x = y) corresponds to a lower connectivity, this is due to a large number of nodes with a connectivity equal to two, not balanced by a higher connectivity of corner nodes. The same remarks hold for the connectivity of the group F laying in the central area of the map. Group E of Strongly Connected Hexatruss is an interesting special case: the mean of (x,y) coordinates is equal to a quarter, resulting in a higher average connectivity due to new bonds between inner nodes of adjacent cells. This higher connectivity happens due to the creation of new strut with neighbor unit-cells: for these specific configurations, the distance between nodes of neighbor cells corresponds to the strut length, resulting in new bonds between nodes. In the G group, another special condition is met when the beam length is equal to ffiffiffi 2 p =2: some new connections appear also with the neighbor cells and increase the connectivity. In this study, two structures exhibit this specific configuration: the octet-truss FCC structure, and (x = 0.3; y = 0.4).
Of course, the analysis of the connectivity by reducing it to its average value only, is a significant simplification. Further consideration of its distribution may reveal a higher degree of details. This analysis is a first attempt to point out specific families of architectures with different levels of connectivity by just considering their topologies. Some families are distinguished due to a specific connectivity (A, B, C in Table 2) and others are rather explained by geometrical condition (F, G in Table 2). The areas D and F (unconnected hexatruss) exhibit the lowest connectivity due to the presence of an important number of nodes with a connectivity equal to 2, not balanced by the high connectivity of corner nodes.

Stiffness constant variation by the density decrease
Stiffness varies mainly with density. Therefore, several levels of density were studied to determine this variation. The relative density was changed by modifying the beam radius, and the resulting stiffness was obtained from FEM model. It allows to compare the topologies based on their stiffness dependence on density.
Calculations cover four levels of density for each topology. Then, a power law is used to fit this set of 4 points: where K is the studied characteristic (here the specific stiffness constant or the anisotropy ratio, but it may apply to other properties), ρ corresponds to the relative density between 0 and 1, and m* is the parameter of the power law. In this study, we choose a power law rather than the conventional Gibson-Ashby model [1], because it is a single parameter law, and it enables the display of all results for one property in a single chart. The parameter m* is determined by the slope measured on a log/ log plot of properties versus relative density. m* is a direct indicator of the decrease of a given property with the density. When equal to one, it corresponds to the mixture rule; when m* N 1 the loss of the properties is greater than the gain in weight, that may be comparable to an antagonistic effect; when m* b 1 then the loss of properties is greater than the gain, and it can be understood as a synergic effect. In this study, all properties have shown an antagonistic effect. This observation is obvious from a mechanical point of view: stress concentration and local deformation fields result in a decrease of stiffness greater than solely the effect of material loss. For all the structures, computed stiffness matrix shows at least cubic properties. This means that the complete 3D description of the constitutive elastic law of designed topologies can be described by three independent constants (C 11 , C 12 and C 44 ): This observation is consistent with the m3m point group properties. At the end of this section, we will study the evolution of these constants with density, in order to determine which kind of topology may be more interesting in terms of tensile/compressive or shear loadings.
The evaluation of the elastic constants with density holds only for a constant topology, assuming the cubic symmetry to be preserved along the deformation path. Therefore, it holds only for the linear part of the elastic behavior. Any effect of the topology change or the density change induced by deformation should be investigated by more complete theories of elasticity of materials, with voids such as micro-dilatation models [41]. The reader interested to extend the current results to large strain elasticity should refer to dedicated works on the continuum micro-dilatation approach applied to lattice structures [42]. In this case, it is expected that Eq. (3) does not hold anymore for large deformation, and its power coefficient should be continuously updated with the topology change along the deformation.
-C 11 and C 12 constant and their related m 11 and m 12 parameters Cubic materials present the particularity (as isotropic materials) to possess a tensile and shear behavior governed by independent constants. Tensile/compressive stiffness of materials only depends on C 11 and C 12 stiffness constants. They are directly associated to the calculation of equivalent Young modulus and Poisson ratio for cubic materials by the relationships [29]: Then C 11 is the parameter governing the stiffness for a uniaxial loading, and C 12 is related to the transversal response induced by the displacement.
The stiffness constants depend on both the density and the topology [1,2]. First, the effect of density is studied, regardless of the variation of the topology. Then the effect of architecture on the deformation mode is examined. An optimal structure for tensile/compressive loading is determined. Fig. 5a illustrates the variation of the specific elastic constant C 11 ¼ C 11 =C 0 11 , where C 11 0 corresponds to the associated C 11 constant of the bulk isotropic material, and Fig. 5b illustrates the variation of C 12 normalized in the same way. It is shown that both constants C 11 and C 12 vary with density following a global power law. The first power exponent is m 11 global = 2.35, and the second is m 12 global = 3. It is also obvious here, that the evolution of these constants confirm an antagonist effect, where the loss of stiffness is even higher than the weight gain. This decreasing trend is more important for C 12 : it indicates that the variation of the transverse component of strain may be more sensitive to the density than to the magnitude of longitudinal strain. If we make the analogy with isotropic materials (and especially for ceramics materials), it is well-known that Poisson ratio is highly related to relative density. It is interesting to show here that this variation exists for all the range of densities.
If we now focus on C 11 , we can observe that all computed topologies tend to fit with the same power law, regardless of architecture variations. At first glance, the reader may think that there is a negligible effect of architecture on C 11 . However, it has to be analyzed carefully while considering the range of applications of lattice structure. With an aim of mass reduction, lattice structures are generally used at low density (ρb0:3). If we draw the attention on this particular field, we can see in Fig. 1.a that the relative constant C 11 exhibits nearly 60% of variation with topology at ρ ¼ 0:3, and N80% of variation when the relative density is below 10%. For large density values, the density plays a first order role after the topology. At lower densities range, this is no more the case, and stiffness constant shows large variation with topology for fixed values of density. It indicates that architecture has a secondary order effect but tends to be more important at low densities.
For C 12 the global trend seems to be the same. It shows a suitable fit with a power law, but with a more important dispersion around the m 12global value. It was possible to describe properly the variation range by using a deviation around m 12 global equal to one (dashed copper lines in Fig. 5). This higher variability of C 12 constant is due to the high sensibility of the Poisson ratio to density. It was leading to larger variations attributed to the architecture change. This high sensitivity to architecture is observed in a larger range of relative densities as compared toC 11 stiffness constant. This variation may be explained by the results from a previous study. In the previous paper, authors have shown that some topologies exhibit auxetic effect (ν b 0) at low densities [2]. This observation may also explain the greater decrease of the C 12 stiffness constant with density.
We will now focus our attention on the effect of architecture on the evolution of these stiffness constants. The m 11 and m 12 parameters were found for each topology, which allow to consider the effect of architecture on the studied stiffness constants. Results were ploted on 2D color maps, and we will attempt here to analyze these maps with regard to previous results on connectivity, in order to explain the different variations.
As expected based on the first part of analysis, the variation of m 12 is larger than m 11 , as it can be seen on the color bars of Fig. 6: m 11 varies from about 1 to 3, and m 12 from 1 to 5.5. We can first observe that all the topologies situated at the corners of the map represent a local minima, which may signify that a single value of node connectivity tend to ensure a lesser decrease of the two stiffness constants C 11 and C 12 with decreasing density. We can also observe that for topologies in the area "connected Hexatruss" (see areas E, G in Table 2) present m 11 and m 12 parameters near to the average of all computed m parameters per constant stiffness (m 11 av = 2,09 and m 12 av = 2,68). These two areas are also associated with an higher connectivity, and it may be the cause for the steep decrease of stiffness constants at lower densities. We can see here that maximas of m 11 and m 11 are situated in areas with the smallest connectivity in map Fig. 2a (corresponding to the group F from Table 2). It means that the more nodes with a small connectivity are present in the structure, the higher is the collapse of properties with decreasing density. This phenomenon is especially pronounced for auxetic structures. This kind of materials often present re-entrant topology, as in the works of Yang et al. [34] or the works of Bückmann et al. [36] on cubic auxetic materials. These cases present two different kinds of node connectivity. It suggest that nodes with lower connectivity have a larger displacement, and they are responsible for the re-entrant effect of auxetic materials.
If we now want to analyze these maps in terms of lattice performance, it should be kept in mind that there are many possible paths in the (x, y, ρ) space to pass from a bulk isotropic material to a near zero density. The choice of relevant architectures has to be done so as to promote a rapid decrease of density while preserving mechanical properties as far as possible [1]. Therefore the best choice for the mechanical engineer would be a m parameter as small as possible. For the first map ( Fig. 6.a) the slowest decrease of properties was obtained for cubic primitive (x = 0, y = 0). It is quite intuitive: in the case of an uniaxial loading, most of the struts of this topology are directed along the loading axis, and the specific stiffness is obviously maximized. Neighbors of cubic primitive (left low corner) were globally a degraded version of cubic primitive, and it explains the smooth decrease of m 11 .
For m 12 map very high values are met in a very localized area around the coordinate x = 0.15; y = 0.25. The use of power law in Eq. (3) exclude the possibility for negative values of C 12 . However, these structures were reported to be auxetic, and it results in a negative Poisson coefficient [2]. According to Eq. (4), this is associated to negative or null values of C 12 . This explains the sudden rise of m 12 for these topologies. The occurrence of auxetic structures has been extensively discussed in previous work [2]. These conclusions can be found again using Fig. 6.b: by some appropriate data processing of m 12 , one can calculate C 12 and determine the Poisson ratio from Eq. (3), making it possible to detect the occurrence of auxetic structures.
-C 44 constant and m 44 parameters In this section, we will focus our attention on the C 44 stiffness constant, which is directly related to the shear modulus in cubic materials through the relation [29]: In Fig. 7, we can see that relative C 44 is well fitted with a power law, with an average m 44 global = 2.6. Here, we can show that m 12 global b m 44 global b m 11 global , which is conform to the observed variation of C 11 , C 12 and C 44 in the work of Li et al. [18]. Unlike the C 11 constant, there is a large variation of the m exponent for fixed values of relative density. Indeed, below only 60% of relative density, we observe a variation up to 36%, and this variation continuously increases to reach N80% for an infill below 30%.
As we concluded for C 12 , the variation of C 44 mainly depends on the density, and only in a second order on the architecture. In other words, architecture seems to play a more important role for the shear modulus than for Young modulus. It also implies that some topologies exhibit a particularly strong decrease of shear modulus compared to the others. Fig. 8 shows the different calculated m 44 parameters for each computed topology. The global variations are very similar to m 11 and m 12 maps, but the trend is reversed. The global minima were transferred from coordinate (0, 0) to the coordinates (0.5, 0.5) which correspond to diag-structure. The (0.5, 0.5) BCC structure has an optimal shear stiffness because the struts are oriented in the principal stress direction, it is the complementary of primitive cube (0, 0) in compression. Crossshaped topologies (see Table 2) also exhibit interesting shear modulus variation. Due to the cross-shape morphology of each face, this kind of topologies may be considered as degraded versions of the diagstructure, in the same way as the neighbors of cubic primitive were presented as "degraded versions" of the optimal C 11 in Fig. 6.a. It is consistent with the results in literature studies on cross-shaped-topologies, X parameter X parameter m 11 parameter m 12 parameter b/ a/ Fig. 6. a. m 11 and b. m 12 parameter for the power law which respectively link C 11 /C 11 0 and C 12 /C 12 0 to the relative density. such as the works of Dirrenberger et al. [37], or Doyoyo et al. [38] As seen in the previous maps (Fig. 6), the area that presented higher connectivity (as connected hexatruss in domains E and G) tends to have the m 44 parameter slightly lower than the average m 44 av = 2.6.
It is another confirmation that connectivity tends to improve the properties in term of shear modulus. As for the stiffness constant, the worst candidate for mass reduction were found in the area of low connectivity. Double-V Hexatruss in zone F (see Fig. 2.b) presents a large number of "V-shaped" struts assemblies along the shear direction. During the loading, it is obvious that these assemblies extend like springs and yield to low shear properties corresponding to high m 44 in Fig. 8.

Anisotropy of cubic lattices
Lattice structures are used for lightweight applications to fill up a volume that is not required to be bulk in the mechanical specifications [3]. Within a complex part, an elementary volume may be subjected to complex loading (tensile, shear …) with significant stress and strain gradients [6]. To control the deformation mode of the macroscopic part, it is of high importance to tailor the anisotropy of elastic properties at the mesoscopic scale. The anisotropy of different topologies is examined now to achieve this point. This study also gives the possibility to compare all topologies by considering simultaneously every stiffness constant that governs the elastic continuum law (Eq. (2)).
For the purpose of this study, the anisotropic ratio (or Zener ratio) is named A; A = 1 represents isotropic materials [29]. When it is N1, C 44 is predominant followed by C 11 and C 12 , in which case the materials could be seen as "shear friendly". When it is below 1, C 11 and C 12 together are predominant over C 44 . That does not directly imply that the material is "tensile friendly" due to the presence of two stiffness constants, but it may simply mean that stiffness constant related to the tensile/compressive behavior is more noticeable than C 44 constant.
As in the previous analysis, we propose here to study the evolution of this ratio with a bulk isotropic material as a reference, and 36 different architectures made of this material. To do so, we introduce the Z parameter which corresponds to the parameter binding the Zener ratio to the relative density by a power law: Three specific situations can be identified: i) Z = 0, corresponds to isotropic materials (log (A) = 0 and A is equal to 1 for all densities), ii) Z b 0, associated to shear friendly materials, and iii) Z N 0, for materials with a predominant tensile stiffness constant. Fig. 9 shows the evolution of calculated Z parameters. This map clearly tends to be separated into three domains. The first one, on the upper right corner, concentrated all the Z b 0 which are directly consistent with the result of C 44 and confirms the presence of shear friendly topologies family around the diag-structure. It also indicated that this shearing preference is compensated by a strong anisotropy of properties. The next domain of this map concerns the areas B, D, F in Table 2, corresponding to "Face topologies" and "Hexatruss". They all have the Z value above 1, in this case possibly indicating a "tensile friendly behavior". The exact behavior of these structures will be discussed just afterwards, once the global trends will be understood. The last domain concerns Connected Hexatruss areas E and G of Table 2. These two areas seem to present Z parameters which tend to conserve the isotropy of properties due to the specific architecture of the matter. Two topologies are particularly preserving A value close to 1: topologies x = 0.25; y = 0.25 with Z = 1.104; and x = 0.2; y = 0.3 with Z = 1.057. These values have to be compared to the octet-truss with Z = 1.636 which is sometimes considered as a reference case for isotropic properties [39,40], therefore it corresponds to a difference of N55%. It was very un-expected to detect these two isotropic structures in this domain, especially when we consider their high level of symmetry. These isotropic structures are resulting from a good compromise between shear and tensile stiffness constants. This in turn, comes directly from the high connectivity in the areas E and G.

Experimental validation through a specific case
For both porous and bulk samples, tests were conducted with a loading direction along the building direction for compression, and with a shear plane normal to the building direction for torsion tests. The tests lead to a value of 116.32 ±1.76 GPa for Young modulus and 45.20 ± 1.84 GPa for shear modulus. These values will be used as reference for the determination of relative stiffness. Z parameter X parameter Y parameter Fig. 9. Z parameter map which rely the anisotropic ratio to the relative density for each topology.
Eq. (5) for shear modulus. The red dashed line corresponds to numerical fitting with the associated equations: with E 0 and G 0 corresponding to the Young and shear modulus of the bulk material determined previously; E* and G* corresponding to the values obtained for the different tests; ρ* and ρ 0 corresponding respectively to porous samples density and bulk samples density. Globally the numerical model tends to be in very good agreement with experimental results. In term of relative density, deviations between numerical model and SLM manufactured samples are −1.20%, +0.32% and 1.56% respectively for the designed values of 24.1%, 37.0% and 44.7%. The structure (0.25, 0.25, 0.5) seems to possess a good manufacturability. This structure only possesses two strut orientations of 67°and 36°in respect to the building direction. Struts with a tilt angle of 36°are commonly accepted within the range of manufacturability. On the contrary the struts with a tilt angle of 67°are under the common criterion of 30°to the building plate. Good manufacturability comes probably from the rather low strut length in this structure, reducing the overhanging fraction and canceling out the negative effects of high tilt angle values. Readers may notice that the observed deviations are not monotonous according to the relative density. These reveals that inerrant defects associated to the manufacturing have different natures or different levels of significance, depending on the tested densities levels.
In Fig. 10a we can see that when the manufactured density is higher than the target values, then the stiffness is also higher for the two highest density levels. For density equal to 24%, the measured densities are higher than the nominal ones, surprisingly the associated stiffness is also above the trend. In term of shear modulus for the lower levels of density, the decrease of density induces a slight decrease of shear modulus. On the contrary, for the two other density levels, the shear modulus is lower than the predicted one despite a higher density than the target. These second observations can be an indicator that different types of defects are induced during manufacturing. The specific effect of each type of defect is not constant in regard to the relative density. In addition, the influence of defects on compressive and shear stiffness tend to be different.
To sum up, experimental testing proved a satisfactory reliability of the numerical model. Despite this good result, it is important to notice that the difference between numerical model and SLM manufactured samples highlights the occurrence of defect with different nature and proportion regarding the relative density. A further study dedicated to the determination and quantification of defects associated to SLM process will be necessary. It also seems necessary to perform an experimental validation on topologies having struts orientations lower than 30°in respect to the building plate.

Global discussion around the elastic properties of cubic lattices structures
In this section, we propose to sum up all the previous results to give a general overview of the possibilities in term of 3D elastic behavior for cubic lattices structures generated using m3m point group. This overview is leading to six specific cases resulting from the previously presented maps. These domains are illustrated in Fig. 11. We detail briefly the particularity of each area regarding the topologies and connectivities, the stiffness constants decreasing parameters, and the anisotropy of properties.
Isotropic area: Structures having a good isotropy of properties at low range of density, mostly for (x = 0.25; y = 0.25 and x = 0.2; y = 0.3). It is mainly explained by the high connectivity of these structures, caused by a connection with neighbored lattice, and specific struts length. They present a suitable compromise between compressions and shears loadings, and they appear to be good candidates for multiaxial loading.
Shear friendly behavior: Topologies that were associated to diagstructure (neighbored topologies and some of cross-shaped structures): they present a slower decrease of C 44 compared to the two other constants. It is attributed to an important fraction of struts oriented in the direction of shear loading. It results in good performances in term of shear modulus balanced by a strong anisotropy of their properties.
Tensile friendly behavior: It is a very similar category to the previous one, but associated to the tensile/compressive aptitude of lattice. These topologies are degraded variations of the primitive cube, which presents a large proportion of strut axially distributed to the loading direction (or slightly misaligned). They also present a strong anisotropy of their properties.
Auxetic area: This area presents re-entrant topologies which leads to auxetic effect. The appearance of this effect may be suggested by strong decreasing of the C 12 relative stiffness constant. It is to be noted that using a power law blurs the complexity of the variations of Poisson ratio, and the exact determination of the auxetic properties is difficult to capture with m 12 alone.
Weak area: These structures present the smallest connectivity. It is leading to a strong decrease of each stiffness constant, which conducts to an important anisotropy of structures' properties. They may be seen as soft-structures, with a very low specific stiffness. It may be interesting for some problematic where very low elastic properties and large elastic strain need to be reached (for example biomimetic and biomechanics of rigid and soft biomaterials, or in the case of vibration damping applications).

Conclusions
In this study, a large database of cubic lattices structures has been generated to produce a large panel of topologies. The following relevant conclusions were drawn: • We proposed here to determine lattices connectivity based on a weighted average of the nodes network. Some specific families of topologies with different levels of connectivity were identified and correlated to the mechanical properties. • The use of a homogenization procedure allowed us to determine the different elastic constants of the homogenized material. The modeling work has revealed that the stiffness matrix always preserves the symmetry of the point group that was used to generate the structures. It was shown that the dependence of elastic constants to the relative density follows a power law. The domains corresponding to a low value of m 11 , m 12 , m 44 on color maps are directly correlated to the domains of high connectivity. The optimal structure maximizing the stiffness for the lowest density is the one minimizing the m parameters for a given loading. It is interesting to notice that the optimal topology varies with the nature of the loading: for the pure tensile strain, the x = 0.0 and y = 0.0 (PC) is optimal, while for pure shear strain the x = 0 0.5 and y = 0.5 (BCC) is optimal. Elastic constants are determined within the framework of small strain, with a constant cubic architecture. Any deviation from this work hypothesis would require using continuum models such as micro-dilatation. The estimation of elastic constants holds for the initial linear stage of elastic deformation. • A parameter quantifying the variation of anisotropic ratio (Zener ratio) was determined from the elastic constants. It was shown that some topologies preserve an attractive isotropy for interestingly low densities below 16%. • The experimental validation of the numerical model on a specific case shows a good agreement with the theoretical result. However, some discrepancy was observed, and it is attributed to the geometric deviation in the manufactured topology.
• The cross-correlation of all color maps including connectivity, elastic constants variations, and anisotropy variation highlight the existence of different behavior families: isotropic, tensile friendly, shear friendly, auxetic, and weak.