Additive manufacturing introduced substructure and computational determination of metamaterials parameters by means of the asymptotic homogenization

Metamaterials exhibit materials response deviation from conventional elasticity. This phenomenon is captured by the generalized elasticity as a result of extending the theory at the expense of introducing additional parameters. These parameters are linked to internal length scales. Describing on a macroscopic level, a material possessing a substructure at a microscopic length scale calls for introducing additional constitutive parameters. Therefore, in principle, an asymptotic homogenization is feasible to determine these parameters given an accurate knowledge on the substructure. Especially in additive manufacturing, known under the infill ratio, topology optimization introduces a substructure leading to higher-order terms in mechanical response. Hence, weight reduction creates a metamaterial with an accurately known substructure. Herein, we develop a computational scheme using both scales for numerically identifying metamaterials parameters. As a specific example, we apply it on a honeycomb substructure and discuss the infill ratio. Such a computational approach is applicable to a wide class substructures and makes use of open-source codes; we make it publicly available for a transparent scientific exchange.


Introduction
Mechanics of metamaterials is gaining an increased interest owing to additive manufacturing technologies allowing us to craft sophisticated structures with different length scales.For weight reduction, material is saved by introducing a substructure.Substructure-related change in materials response is already known [1,2,3], studied under different assumptions [4,5,6,7,8,9], and verified experimentally [10,11,12,13].Substructure-related change leads to metamaterials and this phenomenon is explained by theoretical arguments by assuming conventional elasticity in the smaller length scale (microscale) leading to generalized elasticity in the larger length scale (macroscale) [14,15,16,17,18].
For constructing theories, different length scales are often incorporated in science.For example, consider the microscale being simply the molecular structure or the lattice structure in a crystalline material conferring anisotropy upon the response at the macroscale [19,20,21,22].Another prominent structure-related anisotropy occurs in composite materials, where the microscale is composed of fibers and matrix.The alignment of fibers, and how different plies are stacked up, cause the anisotropy as well as values of effective parameters at the macroscale [23,24,25,26,27].Porous materials are frequently modeled as a full material with voids as given inclusions at the microscale, we refer to [28,29,30,31,32,33]. Additive manufacturing is capable of building metamaterials as demonstrated in [34,35,36,37,38,39].Also adding texture in 3-D printing introduces a substructure.Especially in metal 3-D printing technologies, the microscale itself is anisotropic [40,41,42].We emphasize that at the macroscale, in all examples above, the microscale structure is not detectable such that the materials substructure is smeared out that is called homogenization.
As applied to generalized mechanics, the use of homogenization techniques is challenging [43,44], since generalized mechanics is still evolving [45,46,47].There exist different homogenization techniques [48,49,50,51,52,53,54].In generalized mechanics [55,56], often a Representative Volume Element (RVE) is exploited as in [57,58], although the use of an RVE in generalized mechanics is difficult to justify [59,60].Yet there exist direct approaches [61,62] by computational homogenization methods [63,64,65,66,67,68] as well as techniques based on gamma-convergence [69,70].By means of asymptotic analysis [71,72,73,74] as already applied in [75,76,77,78], we decompose variables into global and local variations [79,80,81] and this separation makes possible to solve the elasticity problem analytically, leading to closed form relations between (known) parameters at the microscale and (sought after) parameters at the macroscale.This approach has been applied in one-dimensional problems for reinforced composites [82,83] and in two-dimensional continuum [84,85,86,87] mostly numerically.From extensive studies [88,89,90,91,92,93], we know that this method is adequate for determining metamaterials parameters.We briefly explain the derivation based on [94] and extend the method to the three-dimensional case by providing a numerical procedure by means of the Finite Element Method (FEM).Especially in honeycomb type infill substructure is our interest [95].The substructure introduces higher order effects as expected and we determine the parameters by a computational homogenization based on the asymptotic analysis.The code uses open-source packages under GNU public license [96] from the FEniCS project [97] and we make the code publicly available in [98] in order to increase the scientific exchange.

Asymptotic homogenization
We begin with the assertion that the deformation energy at the microscale is equivalent to the deformation energy at the macroscale, for the same domain, Ω, occupied by the continuum body.We use the standard continuum mechanics notation with dV meaning the infinitesimal volume element, expressed in Cartesian coordinates as dV = dx dy dz.There is only one coordinate system used for both scales.We use a material frame, so the location of material particles is denoted by X = (X 1 , X 2 , X 3 ) = (x, y, z).Furthermore, we use "m" and "M" denoting microscale and macroscale, respectively.The domains for both scales are equivalent, large enough for allowing homogenization and small enough such that the substructure has a significant effect at the macroscale.We emphasize that a large enough domain-analogously macroscale with a large enough length scale-converges to the classical elasticity approach.
Since we model an elastic body, the deformation energy depends solely on space derivatives of displacements.At each length scale, there exists one displacement field, u m i , u M i .We stress that displacements and their derivatives are different such that the energy density is different in each position.Nevertheless, for the whole body, the total energy is equivalent at both scales.This assertion is the key axiom in nearly all homogenization theories based on the intuition that the energy applied on the body is the same although we observe a different displacement recorded by a 10 MP camera via Digital Image Correlation (DIC) compared to a displacement field captured under a microscope.
We simplify the analysis by assuming that the system at the microscale is composed of linear elastic material(s) such that the energy is quadratic in displacement gradients given by strains, ε m , with a known stiffness tensor, C m , as follows: where the comma denotes a space derivative in X and we understand Einstein's summation convention over repeated indices.For the sake of simplicity, we henceforth use linearized strain measure, and the usual (minor) symmetries of the stiffness, The system at the microscale possesses different materials.For the substructure, for example in an additively manufactured porous structure, we model the structure itself with its stiffness tensor and the voids with a nearly zero stiffness.In other words, the material is heterogeneous at the microscale.At the macroscale, the system is assumed to be homogeneous and to obey materially and geometrically linear strain gradient elasticity modeled by the following deformation energy density: with analogous symmetries We stress that G M = 0 if the macroscale is of centrosymmetric substructure and D M = 0 leads to conventional elasticity with substructure related anisotropy without higher order (strain gradient) terms.
First, we introduce a so-called geometric center: and, assuming enough continuity, approximate the macroscale displacement by the Taylor expansion around the value at the geometric center by truncating after quadratic terms.The choice of quadratic terms is justified by the nonlocality of the theory, in other words, we aim for the strain gradient theory incorporating second derivatives.All higher terms than the second derivative will be neglected.The expansion of displacement gradients reads is a vector evaluated at c X, its gradient vanishes leading to Second, we introduce spatial averaging for displacement gradients by using the latter expansions and the fact that terms evaluated at c X are constant within the domain since integration is additive and we have inserted Eq. ( 6).Thus, we obtain Third, we use the spatial averaged values in the expansions ( 7) and ( 8) Obviously, we circumvent using any spatial averaging techniques [99,100,101].Finally, we insert the latter into the energy definition and take out spatial averaged terms out of the integral By following the asymptotic homogenization method, we use a so-called homothetic ratio, , for a separation of length scales and introduce the local coordinates, Therefore, the macroscale relations in Eq. ( 12) become With the assumption that the displacement field is a smooth function at the macroscale and yperiodic in local coordinates, the mean local fluctuations vanish within the chosen domain, Ω.In other words, the effective property at the macroscale is constant representing the "oscillatory" property at the microscale.The difference between the effective (macroscale) and oscillatory (microscale) property is the fluctuation to vanish.In this regard, we decompose the microscale displacement where n u(X, y) (n = 0, 1, 2) are y-periodic.In other words, the chosen domain, Ω, acts as a Representative Volume Element (RVE) within that we seek the effective property.
We use the well-known least action principle for solving the displacement by starting off with the Lagrange function, ρf i u m i − w m , where the gravitational specific (per mass) force, f i , and the mass density, ρ, are given.For finding the variation of the action functional by the arbitrary test functions, δu, we perform an integration by part where the domain boundaries, ∂Ω, are identical to those from neighboring RVEs.Since the normal vectors, n, of neighboring surfaces, dA, are opposite, all boundaries vanish 0 =δ Derivative of the microscale displacement from Eq. ( 17) after inserting Eq. ( 15) reads Inserting the latter in Eq. ( 18) and once more using the chain rule in combination with Eq. ( 15), we obtain where separation of coefficients multiplied by the same order in and setting every term zero-since and 2 terms are independent-results in Since C m ijkl depends on y, for example consider two distinct materials at the microscale, from the first relation, we immediately conclude that 0 u i = 0 u i (X).By using this dependency, we introduce the multiplicative decomposition with the unknown tensors ϕ abc and ψ abcd .The latter decomposition is a general procedure in tensor calculus and the unknown tensors, ϕ, ψ, have no underlying assumptions.As a consequence, for u m , we have the following expression: with the first term-the sole term depending only on X, all the other terms depend on y as well-corresponding to the macroscale displacement, By using Eq. ( 24) in Eq. ( 23), we obtain the displacement gradient, and, after inserting Eq. ( 16), we acquire since we incorporate up to the second gradients in Eq. (7).By using M abcij = y c L abij + N abcij we calculate the energy at the microscale with Immediately we observe by comparing with Eq. ( 13), where Īkn = Therefore, C M , D M , G M are determined once ϕ and ψ are calculated by using the substructure.For these variables, we will obtain corresponding field equations in the following.
By inserting Eq. ( 23) in Eq. ( 21) 2 and using Analogously, by exploiting Eq. ( 21) 3 and inserting the latter, we acquire since we incorporate only up to the second derivative in Eq. ( 7).
In the case of the macroscale, with the least action principle by means of the Lagrange function, ρf i u M i − w M , we obtain after using integration by parts twice and letting the domain boundaries vanish 0 =δ since the stiffness tensors are constant at the macroscale, as well as we incorporate only up to the second derivative in Eq. (7).By using this relation in Eq. ( 32), we get By solving Eq.( 31) and Eq. ( 35) 2 , we calculate ϕ and ψ.

Method of solution
We sum up the methodology proposed herein.Consider a metamaterial with a given substructure at the microscale, y.Modeling the substructure with the given C m by means of the finite element method leads to a numerical solution of ϕ by satisfying Eq.( 31): By using the solution, from Eqs. ( 28), ( 29), we determine The macroscale stiffness tensor, C M , is used in Eq.( 35) 2 in order to acquire ψ by fulfilling With this solution, we construct and determine The outcome is determining the components of C M tensor of rank four, G M tensor of rank five, and D M tensor of rank six.
In particular, for the numerical solution of Eq. ( 36) as well as Eq. ( 38), we follow the standard procedure of the finite element method [102] and utilize a finite dimensional Hilbertian Sobolev space for trial functions.The same space is used for the test functions as well, called the Galerkin procedure.The triangulation of the structure in y is established by using tetrahedrons, and we solve the discrete problem by minimizing the weak form.In order to get the weak forms, Eqs. ( 36), (38) are multiplied by arbitrary test functions of their ranks for reducing to a scalar integrated over the volume of the structure, Ω.For fulfilling the y periodicity, all boundaries are modeled as periodic boundaries by tying the nodes on corresponding surfaces.In other words, for a cube from left to right along X 1 -axis, each node, say, on the left surface has to have the same displacement as its counterpart with the same X 2 , X 3 coordinates on the right surface.Hence, technically, all boundaries are of Dirichlet type and the test functions vanish on all boundaries, for an alternative approach of weak periodicity, we refer to [103].We use herein a strong coupling with the same mesh on corresponding boundaries, since we use the RVE only at the level of parameter determination.
All the implementation is carried out in the FEniCS platform, we refer to [104] for an introduction with examples.The weak form is obtained after integrating by parts, we stress that the periodic boundary condition causes that boundary integrals vanish.Moreover, we omit distinguishing between the functions and their discrete representations, since they never occur in the same equation.
In order to calculate ϕ and ψ, by utilizing Eq.( 31) and Eq. ( 35) 2 , we obtain the following weak forms: are solved separately by setting a,b,c indices.This fact is of importance so we write out explicitly, how it is meant to do.Because of the minor symmetry, C M ijkl = C M ijlk , we know that L abkl = L bakl and ϕ abi = ϕ bai such that we solve six weak forms in order to obtain ϕ 11i , ϕ 22i , ϕ 33i , ϕ 23i , ϕ 13i , ϕ 12i , respectively.We use these values in Eq. ( 37).This method is admissible under the assumption that for each ab in Voigt's notation indices, ϕ components are per se independent.Also the use in Eq. ( 37) is justified since we obtain 21 components of the stiffness tensor as follows: Of course, depending on the substructure, it may be the case that some of ϕ components are equivalent; however, this symmetry is metamaterial specific.In the same manner, from Eq. ( 41), we use ψ abci = ψ baci and for i = 1 we solve for i = 3 In this way, we solve for ψ 11c1 . . .ψ 12c3 separately and use them to obtain G M and D M by means of Eq. ( 40).

Results and discussion
By virtue of 3-D printers, it is possible to manufacture complex structures with voids inside.Voids result in a porous structure at the microscale.We stress that the voids are introduced on purpose and we assume that the microscale material is full.For example in Fused Deposition Modeling (FDM), the filaments are made of non-porous material and the porosity is caused by design.This layer-by-layer manufacturing technique is coded by a software called slicer.Slicer converts the structure from the CAD design into a G-code providing the motion of the nozzle laying the melt material, i.e. print the material as a thick viscous fluid located at the given positions.For the purpose of weight reduction, all slicer softwares introduce an infill ratio, exchanging the full material with a pre-configured periodic lattice structure.Decreasing the infill ratio increases the porosity at the macroscale.One such typical honeycomb structure is a hexagonal lattice configuration as seen in Fig. 1, the CAD is utilized in Salome, the open-source integration platform for numerical simulation.The full material is replaced with this configuration, for which we compute the higher order terms for any homothetic ratio, , with the assumption that the linear isotropic material at the microscale might be linear anisotropic strain gradient at the macroscale.For the particular RVE as seen in Fig. 1, the homothetic ratio is unity, i.e. the infill ratio is around 50% meaning that the half of the space is filled with the (orange) material.The homothetic ratio is inversely related to the infill ratio, for decreasing the infill ratio increases, where = 0 reads 100% infill ratio meaning that the material is full and no substructure emerges.Obviously, for 100% infill ratio, the higher order terms, G M , D M vanish in Eq. (40).
By using the RVE, the mesh is generated in Salome by using NetGen and Mephisto algorithms as seen in Fig. 2. We emphasis that the periodic boundary conditions need corresponding meshes on the "neighboring" surfaces.An example is demonstrated in Fig. 3, where along the X 1 = X axis, the boundary surfaces are visible.All nodes on both surfaces have the same X 2 = Y and X 3 = Z coordinates such that the degrees of freedom on each node are set equivalent to the corresponding node on the neighboring surface.As the periodic boundaries reflect the given solution, they are Dirichlet boundary conditions, which means that the macroscale and microscale solutions match along the boundaries as well.Although this condition is not a priori set into the formulation, the use of RVE enforces matching boundaries.From the computational point of view, using Dirichlet   3: Demonstration of the periodic boundary conditions along X-axis, the same mesh is used such that the Y and Z coordinates are matching for nodes to be defined as the same degree of freedom boundary conditions on all surfaces, makes the problem well-defined.Hence, there are no emerging numerical problems, where we used multifrontal massively parallel sparse direct solver (mumps) for solving the weak forms and Gaussian quadrature for integration.
As usual, we write out the stiffness tensor in Voigt's notation with A, B standing for combination of two indices in the order: 11, 22, 33, 23, 13, 12 such that the rank four tensor, C M ijkl , is represented in a matrix notation, where obviously the major symmetry holds true, C M AB = C M BA , although this identity is not explicitly stated in the notation.Analogously we use α, β for three indices in the order: 111, 221,331,231,131,121,112,222,332,232,132,122,113,223,333,233,133,123 in order to be able to represent higher order terms in a matrix form as well.Specifically, for G M ijklm we have and for D M ijklmn we obtain  Computed for an RVE of 240 mm × 277.12 mm × 20 mm along X, Y , Z axes, respectively, made of an isotropic material with the Young's modulus of 110 GPa and Poisson's ratio of 0.35, we demonstrate the results in Voigt-like notation introduced above.For the stiffness tensor, we obtain 16 10 9 0 0 0 10 11 7 0 0 0 9 7 43 0 0 0 0 0 0 8 0 0 0 0 0 0 8 0 0 0 0 0 0 3 where we round off 0.1 GPa in all components.For the higher order terms, results depend on the arbitrary infill ratio set by the homothetic ratio , as follows: with ±0.1 kN/mm accuracy as well as with 0.1 TN accuracy, where 1 TN =10 12 N.A general sensitivity analysis of higher order terms is inadequate, in other words, comparison between the displacement altering because of G M and D M components is impossible.The structure dependence on the homothetic ratio as well as loading and boundary conditions affect the sensitivity.Therefore, we have written out all terms with their own accuracy and circumvent ourselves from reducing the complexity of the outcome.
Since the topology is hexagonal, centro-symmetry is lacking such that G M tensor of rank 5 fails to vanish.All components D M ×33××× regarding the second gradient along Z-axis are zero due to the chosen geometry.Obviously, the periodic boundaries along Z-axis create hollow hexagonal tubes without "porosity."Such a porous structure is indeed the case in XY -plane.Therefore, out of XY -plane the homogenization introduces a weakened structure, visible as C M 3333 being less than the half of the Young's modulus of the material itself; however, no higher order terms occur.
It is challenging to directly relate the homothetic ratio to the physical length scale and further studies are necessary in order to justify this study's parameters.

Conclusions
Generalized mechanics has been already studied in 1950s as a purely academic research.Additive manufacturing opens the door for crafting structures with substructures (microscale), called infills, leading to different length scales performing simultaneously at the macroscale, thus, making the generalized elasticity necessary for accurate modeling.Involving strains, conventional elasticity necessitates 21 material parameters.Generalized elasticity with strain gradients introduces additional to the 21 (different) parameters in C M rank 4 tensor, another 108 parameters in G M rank 5, and 171 parameters in D M rank 6 tensors.Asymptotic analysis results in micro-macro-scale relations that we briefly yet thoroughly demonstrated in this work.Finally, a new methodology is proposed for using the substructure and determining all the parameters in generalized elasticity by using computations based on the finite element method (FEM).In order to present the method on a particular case of hexagonal honeycomb substructure, open-source codes based numerical implementation is established under GNU public license [96], the code is available in [98] in order to allow a transparent scientific exchange.

Figure 1 :Figure 2 :
Figure 1: Honeycomb structure in Salome and a possible representative volume element (RVE)shown opaque within the transparent structure, orange denotes the 3-D printed material and gray is void (air) modeled with a significantly low modulus

Figure
Figure 3: Demonstration of the periodic boundary conditions along X-axis, the same mesh is used such that the Y and Z coordinates are matching for nodes to be defined as the same degree of freedom where the symmetry holds true, D M αβ = D M αβ .Therefore, we determine 21 components for C M AB , 108 components for G M Aα , and 171 components for D M αβ in this work for the honeycomb structure by means of the approach explained in Eqs.(36)-(40).