Additive manufacturing introduced substructure and computational determination of metamaterials parameters by means of the asymptotic homogenization
AAdditive manufacturing introducedsubstructure and computational determination of metamaterials parameters by means of the asymptotic homogenization
Bilen Emek Abali , ∗ Emilio Barchiesi Technische Universit¨at Berlin, Institute of Mechanics, MS 2,Einsteinufer 5, 10587 Berlin, Germany Uppsala University, Division of Applied Mechanics,Department of Materials Science and Engineering, Box 534, SE-751 21 Uppsala, Sweden International Research Center on Mathematics and Mechanics of Complex Systems,Universit`a degli Studi dell’Aquila,Via Giovanni Gronchi 18 - Zona industriale di Pile 67100, L’Aquila, ItalyDedicated to Prof. Holm Altenbach on the occasion of his 65 th birthday Abstract
Metamaterials exhibit materials response deviation from conventional elasticity. This phe-nomenon is captured by the generalized elasticity as a result of extending the theory at theexpense of introducing additional parameters. These parameters are linked to internal lengthscales. Describing on a macroscopic level a material possessing a substructure at a microscopiclength scale calls for introducing additional constitutive parameters. Therefore, in principle,an asymptotic homogenization is feasible to determine these parameters given an accurateknowledge on the substructure. Especially in additive manufacturing, known under the infillratio, topology optimization introduces a substructure leading to higher order terms in me-chanical response. Hence, weight reduction creates a metamaterial with an accurately knownsubstructure. Herein, we develop a computational scheme using both scales for numericallyidentifying metamaterials parameters. As a specific example we apply it on a honeycombsubstructure and discuss the infill ratio. Such a computational approach is applicable to awide class substructures and makes use of open-source codes; we make it publicly available fora transparent scientific exchange.
Keywords:
Metamaterials, Homogenization, Generalized mechanics, Finite Element Method (FEM) ∗ Corresponding author, ORCID: 0000-0002-8735-6071, email: [email protected] a r X i v : . [ c s . C E ] O c t Introduction
Mechanics of metamaterials is gaining an increased interest owing to additive manufacturing tech-nologies allowing us to craft sophisticated structures with different length scales. For weightreduction, material is saved by introducing a substructure. Substructure-related change in ma-terials 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 metamaterialsand this phenomenon is explained by theoretical arguments by assuming conventional elasticityin 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 crystallinematerial conferring anisotropy upon the response at the macroscale [19, 20, 21, 22]. Anotherprominent structure-related anisotropy occurs in composite materials, where the microscale iscomposed 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 themicroscale, we refer to [28, 29, 30, 31, 32, 33]. Additive manufacturing is capable of buildingmetamaterials as demonstrated in [34, 35, 36, 37, 38, 39]. Also adding texture in 3-D printingintroduces a substructure. Especially in metal 3-D printing technologies, the microscale itselfis anisotropic [40, 41, 42]. We emphasize that at the macroscale, in all examples above, themicroscale structure is not detectable such that the materials substructure is smeared out that iscalled 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 tech-niques [48, 49, 50, 51, 52, 53, 54]. In generalized mechanics [55, 56], often a Representative VolumeElement (RVE) is exploited as in [57, 58], although the use of an RVE in generalized mechanics isdifficult to justify [59, 60]. Yet there exist direct approaches [61, 62] by computational homogeniza-tion methods [63, 64, 65, 66, 67, 68] as well as techniques based on gamma-convergence [69, 70]. Bymeans of asymptotic analysis [71, 72, 73, 74] as already applied in [75, 76, 77, 78], we decomposevariables into global and local variations [79, 80, 81] and this separation makes possible to solvethe elasticity problem analytically, leading to closed form relations between (known) parameters atthe microscale and (sought after) parameters at the macroscale. This approach has been appliedin 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 thatthis method is adequate for determining metamaterials parameters. We briefly explain the deriva-tion based on [94] and extend the method to the three-dimensional case by providing a numericalprocedure by means of the Finite Element Method (FEM). Especially in honeycomb type infillsubstructure is our interest [95]. The substructure introduces higher order effects as expected andwe determine the parameters by a computational homogenization based on the asymptotic analy-sis. 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.2
Asymptotic homogenization
We begin with the assertion that the deformation energy at the microscale is equivalent to thedeformation energy at the macroscale, (cid:90) Ω w m d V = (cid:90) Ω w M d V , (1)for the same domain, Ω, occupied by the continuum body. We use the standard continuum mechan-ics notation with d V meaning the infinitesimal volume element, expressed in Cartesian coordinatesas d V = d x d y d z . There is only one coordinate system used for both scales. We use a materialframe, so the location of material particles is denoted by X = ( X , X , X ) = ( x, y, z ). Further-more, we use “m” and “M” denoting microscale and macroscale, respectively. The domains forboth scales are equivalent, large enough for allowing homogenization and small enough such thatthe substructure has a significant effect at the macroscale. We emphasize that a large enoughdomain—analogously macroscale with a large enough length scale—converges to the classical elas-ticity approach.Since we model an elastic body, the deformation energy depends solely on space derivatives ofdisplacements. At each length scale, there exists one displacement field, u m i , u M i . We stress thatdisplacements and their derivatives are different such that the energy density is different in eachposition. Nevertheless, for the whole body, the total energy is equivalent at both scales. Thisassertion is the key axiom in nearly all homogenization theories based on the intuition that theenergy applied on the body is the same although we observe a different displacement recorded bya 10 MP camera via Digital Image Correlation (DIC) compared to a displacement field capturedunder a microscope.We simplify the analysis by assuming that the system at the microscale is composed of linear elasticmaterial(s) such that the energy is quadratic in displacement gradients given by strains, ε m , witha known stiffness tensor, C m , as follows: w m = 12 C m ijkl ε m ij ε m kl , ε m ij = 12 ( u m i,j + u m j,i + u m k,i u m k,j ) , (2)where the comma denotes a space derivative in X and we understand Einstein ’s summationconvention over repeated indices. For the sake of simplicity, we henceforth use linearized strainmeasure, ε m ij = 12 ( u m i,j + u m j,i ) , (3)and the usual (minor) symmetries of the stiffness, C m ijkl = C m jikl = C m ijlk , we obtain w m = 12 C m ijkl u m i,j u m k,l . (4)The system at the microscale possesses different materials. For the substructure, for example inan additively manufactured porous structure, we model the structure itself with its stiffness tensorand the voids with a nearly zero stiffness. In other words, the material is heterogeneous at themicroscale. At the macroscale, the system is assumed to be homogeneous and to obey materiallyand geometrically linear strain gradient elasticity modeled by the following deformation energydensity: w M = 12 C M ijkl u M i,j u M k,l + 12 D M ijklmn u M i,jk u M l,mn + G M ijklm u M i,j u M k,lm , (5)with analogous symmetries C M ijkl = C M jikl = C M ijlk as well as D M ijklmn = D M ijklmn = D M ikjlmn = D M lmnijk and G M ijklm = G M jiklm = G M ijkml . We stress that G M = 0 if the macroscale is of centro-symmetric substructure and D M = 0 leads to conventional elasticity with substructure relatedanisotropy without higher order (strain gradient) terms.3irst, we introduce a so-called geometric center: c X = 1 V (cid:90) Ω X d V , (6)and, assuming enough continuity, approximate the macroscale displacement by the
Taylor ex-pansion around the value at the geometric center by truncating after quadratic terms. The choiceof quadratic terms is justified by the nonlocality of the theory, in other words, we aim for the straingradient theory incorporating second derivatives. All higher terms than the second derivative willbe neglected. The expansion of displacement gradients reads u M i ( X ) = u M i (cid:12)(cid:12)(cid:12) c X + u M i,j (cid:12)(cid:12)(cid:12) c X ( X j − c X j ) + 12 u M i,jk (cid:12)(cid:12)(cid:12) c X ( X j − c X j )( X k − c X k ) . (7)Since u M i (cid:12)(cid:12)(cid:12) c X is a vector evaluated at c X , its gradient vanishes leading to u M i,l ( X ) = u M i,j (cid:12)(cid:12)(cid:12) c X δ jl + 12 u M i,jk (cid:12)(cid:12)(cid:12) c X ( δ jl ( X k − c X k ) + ( X j − c X j ) δ kl ) , = u M i,l (cid:12)(cid:12)(cid:12) c X + u M i,lk (cid:12)(cid:12)(cid:12) c X ( X k − c X k ) ,u M i,lm ( X ) = u M i,lk (cid:12)(cid:12)(cid:12) c X δ km = u M i,lm (cid:12)(cid:12)(cid:12) c X . (8)Second, we introduce spatial averaging for displacement gradients by using the latter expansionsand the fact that terms evaluated at c X are constant within the domain (cid:104) u M i,j (cid:105) = 1 V (cid:90) Ω u M i,j d V = u M i,j (cid:12)(cid:12)(cid:12) c X + u M i,jk (cid:12)(cid:12)(cid:12) c X ¯ I k , (cid:104) u M i,jk (cid:105) = 1 V (cid:90) Ω u M i,jk d V = u M i,jk (cid:12)(cid:12)(cid:12) c X , (9)with ¯ I k = 1 V (cid:90) Ω ( X k − c X k ) d V = 1 V (cid:90) Ω X k d V − V (cid:90) Ω c X k d V = 0 , (10)since integration is additive and we have inserted Eq. (6). Thus, we obtain (cid:104) u M i,j (cid:105) = u M i,j (cid:12)(cid:12)(cid:12) c X , (cid:104) u M i,jk (cid:105) = u M i,jk (cid:12)(cid:12)(cid:12) c X . (11)Third, we use the spatial averaged values in the expansions (7) and (8) u M i ( X ) = u M i (cid:12)(cid:12)(cid:12) c X + (cid:104) u M i,j (cid:105) ( X j − c X j ) + 12 (cid:104) u M i,jk (cid:105) ( X j − c X j )( X k − c X k ) ,u M i,j ( X ) = (cid:104) u M i,j (cid:105) + (cid:104) u M i,jk (cid:105) ( X k − c X k ) ,u M i,jk ( X ) = (cid:104) u M i,jk (cid:105) . (12)Obviously, we circumvent using any spatial averaging techniques [99, 100, 101]. Finally, we insert4he latter into the energy definition and take out spatial averaged terms out of the integral (cid:90) Ω w M d V = (cid:90) Ω (cid:18) C M ijlm u M i,j u M l,m + 12 D M ijklmn u M i,jk u M l,mn + G M ijklmn u M i,j u M k,lm (cid:19) d V = 12 C M ijlm (cid:90) Ω u M i,j u M l,m d V + 12 D M ijklmn (cid:90) Ω u M i,jk u M l,mn d V + G M ijklm (cid:90) Ω u M i,j u M k,lm d V = 12 C M ijlm (cid:90) Ω (cid:16) (cid:104) u M i,j (cid:105) + (cid:104) u M i,jk (cid:105) ( X k − c X k ) (cid:17)(cid:16) (cid:104) u M l,m (cid:105) + (cid:104) u M l,mn (cid:105) ( X n − c X n ) (cid:17) d V ++ 12 D M ijklmn (cid:90) Ω (cid:104) u M i,jk (cid:105)(cid:104) u M l,mn (cid:105) d V + G M ijlmn (cid:90) Ω (cid:16) (cid:104) u M i,j (cid:105) + (cid:104) u M i,jk (cid:105) ( X k − c X k ) (cid:17) (cid:104) u M l,mn (cid:105) d V = V (cid:18) C M ijlm (cid:104) u M i,j (cid:105)(cid:104) u M l,m (cid:105) + (cid:0) C M ijlm ¯ I kn + D M ijklmn + 2 G M ijlmn ( X k − c X k ) (cid:1) (cid:104) u M i,jk (cid:105)(cid:104) u M l,mn (cid:105) ++ 2 G M ijlmn (cid:104) u M i,j (cid:105)(cid:104) u M l,mn (cid:105) (cid:19) , (13)by using ¯ I kn = 1 V (cid:90) Ω ( X k − c X k )( X n − c X n ) d V . (14)By following the asymptotic homogenization method, we use a so-called homothetic ratio, (cid:15) , for aseparation of length scales and introduce the local coordinates, y j = 1 (cid:15) ( X j − c X j ) . (15)Therefore, the macroscale relations in Eq. (12) become u M i ( X ) = u M i (cid:12)(cid:12)(cid:12) c X + (cid:15)y j (cid:104) u M i,j (cid:105) + 12 (cid:15) y j y k (cid:104) u M i,jk (cid:105) ,u M i,j ( X ) = (cid:104) u M i,j (cid:105) + (cid:15)y k (cid:104) u M i,jk (cid:105) ,u M i,jk ( X ) = (cid:104) u M i,jk (cid:105) . (16)With the assumption that the displacement field is a smooth function at the macroscale and y -periodic in local coordinates, the mean local fluctuations vanish within the chosen domain, Ω. Inother 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 microscaledisplacement u m ( X ) = u ( X , y ) + (cid:15) u ( X , y ) + (cid:15) u ( X , y ) + O ( (cid:15) ) , (17)where n u ( X , y ) ( n = 0, 1, 2) are y -periodic. In other words, the chosen domain, Ω, acts as aRepresentative 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 themass density, ρ , are given. For finding the variation of the action functional by the arbitrary testfunctions, δ u , we perform an integration by part where the domain boundaries, ∂ Ω, are identicalto those from neighboring RVEs. Since the normal vectors, n , of neighboring surfaces, d A , are5pposite, all boundaries vanish0 = δ (cid:90) Ω (cid:0) ρf i u m i − w m (cid:1) d V , (cid:90) Ω (cid:0) ρf i δ u m i − C m ijkl u m k,l δ u m i,j (cid:1) d V , (cid:90) Ω (cid:16) ρf i + (cid:0) C m ijkl u m k,l (cid:1) ,j (cid:17) δ u m i d V − (cid:90) ∂ Ω C m ijkl u m k,l n j δ u m i d A , ρf i + (cid:0) C m ijkl u m k,l (cid:1) ,j . (18)Derivative of the microscale displacement from Eq. (17) after inserting Eq. (15) reads u m i,j = (cid:16) u i ( X , y ) + (cid:15) u i ( X , y ) + (cid:15) u i ( X , y ) + O ( (cid:15) ) (cid:17) ,j = u i,j + (cid:15) u i,j + (cid:15) u i,j + δ kj (cid:15) ∂∂y k (cid:16) u i + (cid:15) u i + (cid:15) u i (cid:17) + O ( (cid:15) )= u i,j + ∂ u i ∂y j (cid:15) + (cid:15) u i,j + ∂ u i ∂y j + (cid:15) u i,j + (cid:15) ∂ u i ∂y j + O ( (cid:15) ) . (19)Inserting the latter in Eq. (18) and once more using the chain rule in combination with Eq. (15),we obtain ρf i + (cid:32) C m ijkl (cid:18) u k,l + ∂ u k ∂y l (cid:15) + (cid:15) u k,l + ∂ u k ∂y l + (cid:15) u k,l + (cid:15) ∂ u k ∂y l (cid:19)(cid:33) ,j ++ 1 (cid:15) ∂∂y j (cid:32) C m ijkl (cid:18) u k,l + ∂ u k ∂y l (cid:15) + (cid:15) u k,l + ∂ u k ∂y l + (cid:15) u k,l + (cid:15) ∂ u k ∂y l (cid:19)(cid:33) = 0 (20)where separation of coefficients multiplied by the same order in (cid:15) and setting every term zero—since (cid:15) and (cid:15) terms are independent—results in1 (cid:15) ∂∂y j (cid:16) C m ijkl ∂ u k ∂y l (cid:17) = 0 , (cid:15) (cid:32)(cid:18) C m ijkl ∂ u k ∂y l (cid:19) ,j + ∂∂y j (cid:0) C m ijkl u k,l (cid:1) + ∂∂y j (cid:16) C m ijkl ∂ u k ∂y l (cid:17)(cid:33) = 0 ,ρf i + (cid:0) C m ijkl u k,l (cid:1) ,j + (cid:16) C m ijkl ∂ u k ∂y l (cid:17) ,j + ∂∂y j (cid:0) C m ijkl u k,l (cid:1) + ∂∂y j (cid:16) C m ijkl ∂ u k ∂y l (cid:17) = 0 ,(cid:15) (cid:32)(cid:0) C m ijkl u k,l (cid:1) ,j + (cid:16) C m ijkl ∂ u k ∂y l (cid:17) ,j + ∂∂y j (cid:0) C m ijkl u k,l (cid:1)(cid:33) = 0 ,(cid:15) (cid:0) C m ijkl u k,l (cid:1) ,j = 0 . (21)Since C m ijkl depends on y , for example consider two distinct materials at the microscale, from thefirst relation, we immediately conclude that u i = u i ( X ). By using this dependency, we introducethe multiplicative decomposition u i = u a,b ( X ) ϕ abi ( y ) , u i = u a,bc ( X ) ψ abci ( y ) , (22)with the unknown tensors ϕ abc and ψ abcd . The latter decomposition is a general procedure in tensorcalculus and the unknown tensors, ϕ , ψ , have no underlying assumptions. As a consequence, for u m , we have the following expression: u m i = u i ( X ) + (cid:15) u a,b ( X ) ϕ abi ( y ) + (cid:15) u a,bc ( X ) ψ abci ( y ) + O ( (cid:15) ) , (23)6ith the first term—the sole term depending only on X , all the other terms depend on y aswell—corresponding to the macroscale displacement, u M = u ( X ) . (24)By using Eq. (24) in Eq. (23), we obtain the displacement gradient, u m i,j = (cid:16) u M i + (cid:15)u M a,b ϕ abi + (cid:15) u M a,bc ψ abci (cid:17) ,j + O ( (cid:15) )= u M i,j + ∂ϕ abi ∂y j u M a,b + (cid:15)ϕ abi u M a,bj + (cid:15) ∂ψ abci ∂y j u M a,bc + (cid:15) ψ abci u M a,bcj + O ( (cid:15) )= (cid:16) δ ia δ jb + ∂ϕ abi ∂y j (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) L abij u M a,b + (cid:15)u M a,bc (cid:16) ϕ abi δ jc + ∂ψ abci ∂y j (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) N abcij + (cid:15) ψ abci u M a,bcj + O ( (cid:15) ) , (25)and, after inserting Eq. (16), we acquire u m i,j = L abij (cid:104) u M a,b (cid:105) + (cid:15) (cid:104) u M a,bc (cid:105) y c L abij + (cid:15) (cid:104) u M a,bc (cid:105) N abcij , (26)since we incorporate up to the second gradients in Eq. (7). By using M abcij = y c L abij + N abcij wecalculate the energy at the microscale (cid:90) Ω w m d V = 12 (cid:90) Ω P (cid:16) C m ijkl L abij L cdkl (cid:104) u M a,b (cid:105)(cid:104) u M c,d (cid:105) + 2 (cid:15)C m ijkl L abij M cdekl (cid:104) u M a,b (cid:105)(cid:104) u M c,de (cid:105) ++ (cid:15) C m ijkl M abcij M defkl (cid:104) u M a,bc (cid:105)(cid:104) u M d,ef (cid:105) (cid:17) d V = V (cid:16) ¯ C abcd (cid:104) u M a,b (cid:105)(cid:104) u M c,d (cid:105) + 2 ¯ G abcde (cid:104) u M a,b (cid:105)(cid:104) u M c,de (cid:105) + ¯ D abcdef (cid:104) u M a,bc (cid:105)(cid:104) u M d,ef (cid:105) (cid:17) . (27)with ¯ C abcd = 1 V (cid:90) Ω C m ijkl L abij L cdkl d V , ¯ G abcde = (cid:15)V (cid:90) Ω C m ijkl L abij M cdekl d V , ¯ D abcdef = (cid:15) V (cid:90) Ω C m ijkl M abcij M defkl d V . (28)Immediately we observe by comparing with Eq. (13), C M ijlm = ¯ C ijlm ,G M ijlmn = ¯ G abcde ,C M ijlm ¯ I kn + D M ijklmn + 2 (cid:15)y k G M ijlmn = ¯ D ijklmn , (29)where ¯ I kn = (cid:90) Ω P ( X k − c X k )( X n − c X n ) d V = (cid:15) (cid:90) Ω P y k y n d V . (30)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.7y inserting Eq. (23) in Eq. (21) and using u = u ( X ), we obtain ∂∂y j (cid:0) C m ijkl u k,l (cid:1) + ∂∂y j (cid:16) C m ijkl ∂ u a,b ϕ abk ∂y l (cid:17) =0 ,∂C m ijkl ∂y j δ ak δ bl u a,b + ∂∂y j (cid:16) C m ijkl ∂ϕ abk ∂y l (cid:17) u a,b =0 ,∂∂y j (cid:32) C m ijkl (cid:16) δ ak δ bl + ∂ϕ abk ∂y l (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) L abkl (cid:33) =0 . (31)Analogously, by exploiting Eq. (21) and inserting the latter, we acquire ρf i + (cid:0) C m ijkl u k,l (cid:1) ,j + (cid:16) C m ijkl ∂ u a,b ϕ abk ∂y l (cid:17) ,j + ∂∂y j (cid:0) C m ijkl u a,bl ϕ abk (cid:1) + ∂∂y j (cid:16) C m ijkl ∂ u a,bc ψ abck ∂y l (cid:17) =0 ,ρf i + C m ijkl u k,lj + C m ijkl u a,bj ∂ϕ abk ∂y l + ∂∂y j (cid:0) C m ijkl ϕ abk (cid:1) u a,bl + ∂∂y j (cid:16) C m ijkl ∂ψ abck ∂y l (cid:17) u a,bc =0 ,ρf i + C m ickl u a,bc (cid:32) δ ak δ bl + ∂ϕ abk ∂y l (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) L abkl + u a,bc ∂∂y j (cid:32) C m ijkl (cid:16) ϕ abk δ cl + ∂ψ abck ∂y l (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) N abckl (cid:33) =0(32)Equations (21) , are identically fulfilled (cid:0) C m ijkl u k,l (cid:1) ,j + (cid:16) C m ijkl ∂ u k ∂y l (cid:17) ,j + ∂∂y j (cid:16) C m ijkl u k,l (cid:17) =0 ,C m ijkl u a,blj ϕ abk + C m ijkl ∂ u a,bcj ψ abck ∂y l + ∂∂y j (cid:16) C m ijkl u a,bcl ψ abck (cid:17) =0 ,C m ijkl u k,lj =0 ,C m ijkl u a,bclj ψ abck =0 , (33)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 boundariesvanish0 = δ (cid:90) Ω (cid:0) ρf i u M i − w M (cid:1) d V , (cid:90) Ω (cid:16) ρf i δ u M i − C M ijkl u M k,l δ u M i,j − D M ijklmn u M l,mn δ u M i,jk − G M ijklm δ u M i,j u M k,lm − G M ijklm u M i,j δ u M k,lm (cid:17) d V , ρf i + C M ijkl u M k,lj − D M ijklmn u M l,mnjk + G M ijklm u M k,lmj − G M kjilm u M k,jlm , ρf i + C M ijkl u M k,lj , (34)since the stiffness tensors are constant at the macroscale, as well as we incorporate only up to thesecond derivative in Eq. (7). By using this relation in Eq. (32), we get − C M icab u M a,bc + C m ickl u a,bc L abkl + u a,bc ∂∂y j (cid:32) C m ijkl N abckl (cid:33) =0 , − C M icab + C m ickl L abkl + ∂∂y j (cid:16) C m ijkl N abckl (cid:17) =0 . (35)8y solving Eq.(31) and Eq. (35) , we calculate ϕ and ψ . We sum up the methodology proposed herein. Consider a metamaterial with a given substructureat the microscale, y . Modeling the substructure with the given C m by means of the finite elementmethod leads to a numerical solution of ϕ by satisfying Eq.(31): ∂∂y j (cid:32) C m ijkl L abkl (cid:33) = 0 , L abkl = δ ak δ bl + ∂ϕ abk ∂y l . (36)By using the solution, from Eqs. (28), (29), we determine C M abcd = ¯ C abcd = 1 V (cid:90) Ω C m ijkl L abij L cdkl d V . (37)The macroscale stiffness tensor, C M , is used in Eq.(35) in order to acquire ψ by fulfilling − C M icab + C m ickl L abkl + ∂∂y j (cid:16) C m ijkl N abckl (cid:17) = 0 , N abckl = ϕ abk δ cl + ∂ψ abck ∂y l . (38)With this solution, we construct M abcij = y c L abij + N abcij , ¯ I kn = (cid:15) (cid:90) Ω P y k y n d V . (39)and determine G M abcde = ¯ G abcde = (cid:15)V (cid:90) Ω C m ijkl L abij M cdekl d V , ¯ D abcdef = (cid:15) V (cid:90) Ω C m ijkl M abcij M defkl d V ,D M ijklmn = ¯ D ijklmn − C M ijlm ¯ I kn − (cid:15)y k G M ijlmn . (40)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 standardprocedure of the finite element method [102] and utilize a finite dimensional Hilbert ian
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 solvethe 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 overthe volume of the structure, Ω. For fulfilling the y periodicity, all boundaries are modeled asperiodic boundaries by tying the nodes on corresponding surfaces. In other words, for a cube fromleft to right along X -axis, each node, say, on the left surface has to have the same displacementas its counterpart with the same X , X coordinates on the right surface. Hence, technically,all boundaries are of Dirichlet type and the test functions vanish on all boundaries, for analternative approach of weak periodicity, we refer to [103]. We use herein a strong coupling withthe same mesh on corresponding boundaries, since we use the RVE only at the level of parameterdetermination.All the implementation is carried out in the FEniCS platform, we refer to [104] for an introductionwith examples. The weak form is obtained after integrating by parts, we stress that the periodic9oundary condition causes that boundary integrals vanish. Moreover, we omit distinguishing be-tween 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) , we obtain the following weakforms: (cid:90) Ω C m ijkl L abkl ∂ δ ϕ abi ∂y j d V = 0 , (cid:90) Ω (cid:18) − C M icab δ ψ abci + C m ickl L abkl δ ψ abci − (cid:16) C m ijkl N abckl (cid:17) ∂ δ ψ abci ∂y j (cid:19) d V = 0 , (41)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 (cid:90) Ω C m ijkl L kl ∂ δ ϕ i ∂y j d V = 0 , (cid:90) Ω C m ijkl L kl ∂ δ ϕ i ∂y j d V = 0 , (cid:90) Ω C m ijkl L kl ∂ δ ϕ i ∂y j d V = 0 , (cid:90) Ω C m ijkl L kl ∂ δ ϕ i ∂y j d V = 0 , (cid:90) Ω C m ijkl L kl ∂ δ ϕ i ∂y j d V = 0 , (cid:90) Ω C m ijkl L kl ∂ δ ϕ i ∂y j d V = 0 , (42)in order to obtain ϕ i , ϕ i , ϕ i , ϕ i , ϕ i , ϕ i , 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 21components of the stiffness tensor as follows: C M1111 = 1 V (cid:90) Ω C m ijkl L ij L kl d V , L kl = δ k δ l + ∂ϕ k ∂y l ,C M1122 = 1 V (cid:90) Ω C m ijkl L ij L kl d V , L kl = δ k δ l + ∂ϕ k ∂y l ,. . .C M1212 = 1 V (cid:90) Ω C m ijkl L ij L kl d V , L kl = δ k δ l + ∂ϕ k ∂y l . (43)Of course, depending on the substructure, it may be the case that some of ϕ components areequivalent; however, this symmetry is metamaterial specific. In the same manner, from Eq. (41),we use ψ abci = ψ baci and for i = 1 we solve (cid:90) Ω (cid:18) − C M1 c δ ψ c + C m1 ckl L kl δ ψ c − (cid:16) C m1 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 , (cid:90) Ω (cid:18) − C M1 c δ ψ c + C m1 ckl L kl δ ψ c − (cid:16) C m1 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 , (cid:90) Ω (cid:18) − C M1 c δ ψ c + C m1 ckl L kl δ ψ c − (cid:16) C m1 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 , (cid:90) Ω (cid:18) − C M1 c δ ψ c + C m1 ckl L kl δ ψ c − (cid:16) C m1 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 , (cid:90) Ω (cid:18) − C M1 c δ ψ c + C m1 ckl L kl δ ψ c − (cid:16) C m1 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 , (cid:90) Ω (cid:18) − C M1 c δ ψ c + C m1 ckl L kl δ ψ c − (cid:16) C m1 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 , (44)10or i = 2 (cid:90) Ω (cid:18) − C M2 c δ ψ c + C m2 ckl L kl δ ψ c − (cid:16) C m2 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 ,. . . (cid:90) Ω (cid:18) − C M2 c δ ψ c + C m2 ckl L kl δ ψ c − (cid:16) C m2 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 , (45)for i = 3 (cid:90) Ω (cid:18) − C M3 c δ ψ c + C m3 ckl L kl δ ψ c − (cid:16) C m3 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 ,. . . (cid:90) Ω (cid:18) − C M3 c δ ψ c + C m3 ckl L kl δ ψ c − (cid:16) C m3 jkl N ckl (cid:17) ∂ δ ψ c ∂y j (cid:19) d V = 0 . (46)In this way, we solve for ψ c . . . ψ c separately and use them to obtain G M and D M by meansof Eq. (40). By virtue of 3-D printers, it is possible to manufacture complex structures with voids inside. Voidsresult in a porous structure at the microscale. We stress that the voids are introduced on purposeand 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. Thislayer-by-layer manufacturing technique is coded by a software called slicer. Slicer converts thestructure from the CAD design into a G-code providing the motion of the nozzle laying the meltmaterial, i.e. print the material as a thick viscous fluid located at the given positions. For thepurpose of weight reduction, all slicer softwares introduce an infill ratio, exchanging the full materialwith a pre-configured periodic lattice structure. Decreasing the infill ratio increases the porosityat the macroscale. One such typical honeycomb structure is a hexagonal lattice configuration asseen in Fig.1, the CAD is utilized in Salome, the open-source integration platform for numericalsimulation. The full material is replaced with this configuration, for which we compute the higherorder terms for any homothetic ratio, (cid:15) , with the assumption that the linear isotropic material atthe microscale might be linear anisotropic strain gradient at the macroscale. For the particularRVE as seen in Fig.1, the homothetic ratio is unity, i.e. the infill ratio is around 50% meaningthat the half of the space is filled with the (orange) material. The homothetic ratio is inverselyrelated to the infill ratio, for decreasing (cid:15) the infill ratio increases, where (cid:15) = 0 reads 100% infillratio meaning that the material is full and no substructure emerges. Obviously, for 100% infillratio, 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 asseen in Fig. 2. We emphasis that the periodic boundary conditions need corresponding meshes onthe “neighboring” surfaces. An example is demonstrated in Fig. 3, where along the X = X axis,the boundary surfaces are visible. All nodes on both surfaces have the same X = Y and X = Z coordinates such that the degrees of freedom on each node are set equivalent to the correspondingnode 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 matchalong the boundaries as well. Although this condition is not a priori set into the formulation, theuse of RVE enforces matching boundaries. From the computational point of view, using
Dirichlet X -axis, the same mesh is usedsuch that the Y and Z coordinates are matching for nodes to be defined as the samedegree of freedomboundary conditions on all surfaces, makes the problem well-defined. Hence, there are no emergingnumerical problems, where we used multifrontal massively parallel sparse direct solver (mumps)for solving the weak forms and Gauss ian quadrature for integration.As usual, we write out the stiffness tensor in
Voigt ’s notation with A , B standing for combinationof two indices in the order: 11, 22, 33, 23, 13, 12 such that the rank four tensor, C M ijkl , is representedin a matrix notation, C M AB = C M1111 C M1122 C M1133 C M1123 C M1113 C M1112 C M2211 C M2222 C M2233 C M2223 C M2213 C M2212 C M3311 C M3322 C M3333 C M3323 C M3313 C M3312 C M2311 C M2322 C M2333 C M2323 C M2313 C M2312 C M1311 C M1322 C M1333 C M1323 C M1313 C M1312 C M1211 C M1222 C M1233 C M1223 C M1213 C M1212 , (47)where obviously the major symmetry holds true, C M AB = C M BA , although this identity is not ex-plicitly 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 torepresent higher order terms in a matrix form as well. Specifically, for G M ijklm we have G M Aα = G M11111 G M11221 G M11331 G M11231 G M11131 G M11121 G M11112 G M11222 G M11332 G M11232 G M11132 G M11122 G M11113 G M11223 G M11333 G M11233 G M11133 G M11123 G M22111 G M22221 G M22331 G M22231 G M22131 G M22121 G M22112 G M22222 G M22332 G M22232 G M22132 G M22122 G M22113 G M22223 G M22333 G M22233 G M22133 G M22123 G M33111 G M33221 G M33331 G M33231 G M33131 G M33121 G M33112 G M33222 G M33332 G M33232 G M33132 G M33122 G M33113 G M33223 G M33333 G M33233 G M33133 G M33123 G M23111 G M23221 G M23331 G M23231 G M23131 G M23121 G M23112 G M23222 G M23332 G M23232 G M23132 G M23122 G M23113 G M23223 G M23333 G M23233 G M23133 G M23123 G M13111 G M13221 G M13331 G M13231 G M13131 G M13121 G M13112 G M13222 G M13332 G M13232 G M13132 G M13122 G M13113 G M13223 G M13333 G M13233 G M13133 G M13123 G M12111 G M12221 G M12331 G M12231 G M12131 G M12121 G M12112 G M12222 G M12332 G M12232 G M12132 G M12122 G M12113 G M12223 G M12333 G M12233 G M12133 G M12123 , (48)13nd for D M ijklmn we obtain D M αβ = D M111111 D M111221 D M111331 D M111231 D M111131 D M111121 D M111112 D M111222 D M111332 D M111232 D M111132 D M111122 D M111113 D M111223 D M111333 D M111233 D M111133 D M111123 D M221111 D M221221 D M221331 D M221231 D M221131 D M221121 D M221112 D M221222 D M221332 D M221232 D M221132 D M221122 D M221113 D M221223 D M221333 D M221233 D M221133 D M221123 D M331111 D M331221 D M331331 D M331231 D M331131 D M331121 D M331112 D M331222 D M331332 D M331232 D M331132 D M331122 D M331113 D M331223 D M331333 D M331233 D M331133 D M331123 D M231111 D M231221 D M231331 D M231231 D M231131 D M231121 D M231112 D M231222 D M231332 D M231232 D M231132 D M231122 D M231113 D M231223 D M231333 D M231233 D M231133 D M231123 D M131111 D M131221 D M131331 D M131231 D M131131 D M131121 D M131112 D M131222 D M131332 D M131232 D M131132 D M131122 D M131113 D M131223 D M131333 D M131233 D M131133 D M131123 D M121111 D M121221 D M121331 D M121231 D M121131 D M121121 D M121112 D M121222 D M121332 D M121232 D M121132 D M121122 D M121113 D M121223 D M121333 D M121233 D M121133 D M121123 D M112111 D M112221 D M112331 D M112231 D M112131 D M112121 D M112112 D M112222 D M112332 D M112232 D M112132 D M112122 D M112113 D M112223 D M112333 D M112233 D M112133 D M112123 D M222111 D M222221 D M222331 D M222231 D M222131 D M222121 D M222112 D M222222 D M222332 D M222232 D M222132 D M222122 D M222113 D M222223 D M222333 D M222233 D M222133 D M222123 D M332111 D M332221 D M332331 D M332231 D M332131 D M332121 D M332112 D M332222 D M332332 D M332232 D M332132 D M332122 D M332113 D M332223 D M332333 D M332233 D M332133 D M332123 D M232111 D M232221 D M232331 D M232231 D M232131 D M232121 D M232112 D M232222 D M232332 D M232232 D M232132 D M232122 D M232113 D M232223 D M232333 D M232233 D M232133 D M232123 D M132111 D M132221 D M132331 D M132231 D M132131 D M132121 D M132112 D M132222 D M132332 D M132232 D M132132 D M132122 D M132113 D M132223 D M132333 D M132233 D M132133 D M132123 D M122111 D M122221 D M122331 D M122231 D M122131 D M122121 D M122112 D M122222 D M122332 D M122232 D M122132 D M122122 D M122113 D M122223 D M122333 D M122233 D M122133 D M122123 D M113111 D M113221 D M113331 D M113231 D M113131 D M113121 D M113112 D M113222 D M113332 D M113232 D M113132 D M113122 D M113113 D M113223 D M113333 D M113233 D M113133 D M113123 D M223111 D M223221 D M223331 D M223231 D M223131 D M223121 D M223112 D M223222 D M223332 D M223232 D M223132 D M223122 D M223113 D M223223 D M223333 D M223233 D M223133 D M223123 D M333111 D M333221 D M333331 D M333231 D M333131 D M333121 D M333112 D M333222 D M333332 D M333232 D M333132 D M333122 D M333113 D M333223 D M333333 D M333233 D M333133 D M333123 D M233111 D M233221 D M233331 D M233231 D M233131 D M233121 D M233112 D M233222 D M233332 D M233232 D M233132 D M233122 D M233113 D M233223 D M233333 D M233233 D M233133 D M233123 D M133111 D M133221 D M133331 D M133231 D M133131 D M133121 D M133112 D M133222 D M133332 D M133232 D M133132 D M133122 D M133113 D M133223 D M133333 D M133233 D M133133 D M133123 D M123111 D M123221 D M123331 D M123231 D M123131 D M123121 D M123112 D M123222 D M123332 D M123232 D M123132 D M123122 D M123113 D M123223 D M123333 D M123233 D M123133 D M123123 , (49)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 structureby means of the approach explained in Eqs. (36)-(40).Computed for an RVE of 240 mm × .
12 mm ×
20 mm along X , Y , Z axes, respectively, madeof an isotropic material with the Young ’s modulus of 110 GPa and
Poisson ’s ratio of 0 .
35, wedemonstrate the results in
Voigt -like notation introduced above. For the stiffness tensor, weobtain C M AB =
16 10 9 0 0 010 11 7 0 0 09 7 43 0 0 00 0 0 8 0 00 0 0 0 8 00 0 0 0 0 3
GPa , (50)where we round off 0 . (cid:15) , as follows: G M Aα = (cid:15) − − −
21 84 26 − − −
40 4 − −
18 7 11 944 − − −
24 96 30 − − −
46 3 − −
15 6 9 740 − − −
16 63 20 − − −
30 19 − −
83 35 51 410 0 0 0 0 0 3 − −
15 6 9 7 −
14 65 21 − − − − −
16 7 10 8 0 0 0 0 0 0 37 − − − − −
11 11 − − kN/mm , (51)with ± . D M αβ = (cid:15) − − −
58 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 − − −
47 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 − − −
275 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 −
48 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 −
53 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 −
16 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 − − −
77 0 0 0 0 0 0 0 0 00 0 0 0 0 0 − − −
63 0 0 0 0 0 0 0 0 00 0 0 0 0 0 − − −
366 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 −
64 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 −
70 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 −
22 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 TN , (52)with 0 . N. A general sensitivity analysis of higher order terms isinadequate, in other words, comparison between the displacement altering because of G M and D M (cid:15) as well as loadingand boundary conditions affect the sensitivity. Therefore, we have written out all terms with theirown 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 tovanish. All components D M × ××× regarding the second gradient along Z -axis are zero due to thechosen geometry. Obviously, the periodic boundaries along Z -axis create hollow hexagonal tubeswithout “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 M3333 being less thanthe 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 furtherstudies are necessary in order to justify this study’s parameters.
Generalized mechanics has been already studied in 1950s as a purely academic research. Additivemanufacturing opens the door for crafting structures with substructures (microscale), called infills,leading to different length scales performing simultaneously at the macroscale, thus, making thegeneralized elasticity necessary for accurate modeling. Involving strains, conventional elasticitynecessitates 21 material parameters. Generalized elasticity with strain gradients introduces addi-tional to the 21 (different) parameters in C M rank 4 tensor, another 108 parameters in G M rank5, and 171 parameters in D M rank 6 tensors. Asymptotic analysis results in micro-macro-scalerelations that we briefly yet thoroughly demonstrated in this work. Finally, a new methodologyis proposed for using the substructure and determining all the parameters in generalized elasticityby using computations based on the finite element method (FEM). In order to present the methodon a particular case of hexagonal honeycomb substructure, open-source codes based numericalimplementation is established under GNU public license [96], the code is available in [98] in orderto allow a transparent scientific exchange. Acknowledgements
B. E. Abali’s work was partly funded by a grant from the Daimler and Benz Foundation.
References [1] A. Eringen & E. Suhubi (1964)
Nonlinear theory of simple micro-elastic solids , InternationalJournal of Engineering Science, :pp. 189–203[2] R. Mindlin (1964) Micro-structure in linear elasticity , Archive for Rational Mechanics andAnalysis, :pp. 51–78[3] A. Eringen (1968)
Mechanics of micromorphic continua , E. Kr¨oner (Ed.)
Mechanics of Gen-eralized Continua , pp. 18–35, Springer-Verlag, Berlin[4] P. Steinmann (1994)
A micropolar theory of finite deformation and finite rotation multiplica-tive elastoplasticity , International Journal of Solids and Structures, :pp. 1063–1084155] V. A. Eremeyev, L. P. Lebedev, & H. Altenbach (2012)
Foundations of micropolar mechanics ,Springer Science & Business Media[6] C. Polizzotto (2013)
A second strain gradient elasticity theory with second velocity gradientinertia–part i: Constitutive equations and quasi-static behavior , International Journal ofSolids and Structures, :pp. 3749–3765[7] C. Polizzotto (2013)
A second strain gradient elasticity theory with second velocity gradi-ent inertia–part ii: Dynamic behavior , International Journal of Solids and Structures, :pp. 3766–3777[8] E. A. Ivanova & E. N. Vilchevskaya (2016)
Micropolar continuum in spatial description ,Continuum Mechanics and Thermodynamics, :pp. 1759–1780[9] B. E. Abali (2018)
Revealing the physical insight of a length scale parameter in metamaterialsby exploring the variational formulation , Continuum Mechanics and Thermodynamics, :pp. 885–894[10] F. dell’Isola, P. Seppecher, M. Spagnuolo, E. Barchiesi, F. Hild, T. Lekszycki, I. Giorgio,L. Placidi, U. Andreaus, M. Cuomo, S. R. Eugster, A. Pfaff, K. Hoschke, R. Langkemper,E. Turco, R. Sarikaya, A. Misra, M. De Angelo, F. D’Annibale, A. Bouterf, X. Pinelli,A. Misra, B. Desmorat, M. Pawlikowski, C. Dupuy, D. Scerrato, P. Peyre, M. Laudato,L. Manzari, P. G¨oransson, C. Hesch, S. Hesch, P. Franciosi, J. Dirrenberger, F. Maurin,Z. Vangelatos, C. Grigoropoulos, V. Melissinaki, M. Farsari, W. Muller, B. E. Abali,C. Liebold, G. Ganzosch, P. Harrison, R. Drobnicki, L. Igumnov, F. Alzahrani, & T. Hayat(2019)
Advances in pantographic structures: design, manufacturing, models, experimentsand image analyses , Continuum Mechanics and Thermodynamics, :pp. 1231–1282[11] E. Barchiesi, M. Spagnuolo, & L. Placidi (2019)
Mechanical metamaterials: a state of theart , Mathematics and Mechanics of Solids, :pp. 212–234[12] F. dell’Isola, E. Turco, A. Misra, Z. Vangelatos, C. Grigoropoulos, V. Melissinaki, &M. Farsari (2019)
Force–displacement relationship in micro-metric pantographs: Experi-ments and numerical simulations , Comptes Rendus M´ecanique, :pp. 397–405[13] W. H. M¨uller (2020)
The experimental evidence for higher gradient theories , Mechanics ofStrain Gradient Materials , pp. 1–18, Springer[14] C. Pideri & P. Seppecher (1997)
A second gradient material resulting from the homogenizationof an heterogeneous linear elastic medium , Continuum Mechanics and Thermodynamics, :pp. 241–257[15] V. P. Smyshlyaev & K. Cherednichenko (2000)
On rigorous derivation of strain gradienteffects in the overall behaviour of periodic heterogeneous media , Journal of the Mechanicsand Physics of Solids, :pp. 1325–1357[16] P. Seppecher, J.-J. Alibert, & F. dell’Isola (2011)
Linear elastic trusses leading to continuawith exotic mechanical interactions , Journal of Physics: Conference Series , vol. 319, p.012018, IOP Publishing[17] H. Abdoul-Anziz & P. Seppecher (2018)
Strain gradient and generalized continua obtained byhomogenizing frame lattices , Mathematics and mechanics of complex systems, :pp.213–250[18] K. K. Mandadapu, B. E. Abali, & P. Papadopoulos (2018)
On the polar nature and invarianceproperties of a thermomechanical theory for continuum-on-continuum homogenization ,arXiv preprint arXiv:1808.02540[19] A. Reuss (1929)
Berechnung der Fließgrenze von Mischkristallen auf grund der Plas-tizit¨atsbedingung f¨ur Einkristalle. , ZAMM-Journal of Applied Mathematics and Mechan-ics/Zeitschrift f¨ur Angewandte Mathematik und Mechanik, :pp. 49–58[20] Z. Hashin & S. Shtrikman (1962)
On some variational principles in anisotropic and nonho-mogeneous elasticity , Journal of the Mechanics and Physics of Solids, :pp. 335–342[21] B. Shafiro & M. Kachanov (2000)
Anisotropic effective conductivity of materials with non-randomly oriented inclusions of diverse ellipsoidal shapes , Journal of applied physics, :pp. 8561–8569[22] R. Lebensohn, Y. Liu, & P. P. Castaneda (2004)
On the accuracy of the self-consistent pproximation for polycrystals: comparison with full-field numerical simulations , ActaMaterialia, :pp. 5347–5361[23] V. Levin (1976) Determination of composite material elastic and thermoelastic constants ,Mechanics of Solids, :pp. 119–126[24] J. Willis (1977)
Bounds and self-consistent estimates for the overall properties of anisotropiccomposites , Journal of the Mechanics and Physics of Solids, :pp. 185–202[25] V. Kushnevsky, O. Morachkovsky, & H. Altenbach (1998)
Identification of effective propertiesof particle reinforced composite materials , Computational mechanics, :pp. 317–325[26] R. Sburlati, R. Cianci, & M. Kashtalyan (2018)
Hashin’s bounds for elastic properties ofparticle-reinforced composites with graded interphase , International Journal of Solids andStructures, :pp. 224–235[27] N. Shekarchizadeh, R. Jafari Nedoushan, T. Dastan, & H. Hasani (2020)
Experimental andnumerical study on stiffness and damage of glass/epoxy biaxial weft-knitted reinforcedcomposites , Journal of Reinforced Plastics and Composites, p. 0731684420938446[28] J. D. Eshelby (1957)
The determination of the elastic field of an ellipsoidal inclusion, andrelated problems , Proceedings of the Royal Society of London. Series A, Mathematicaland Physical Sciences, pp. 376–396[29] T. Mori & K. Tanaka (1973)
Average stress in matrix and average elastic energy of materialswith misfitting inclusions , Acta Metallurgica, :pp. 571–574[30] S. Kanaun & L. Kudryavtseva (1986)
Spherically layered inclusions in a homogeneous elasticmedium , Journal of Applied Mathematics and Mechanics, :pp. 483–491[31] Z. Hashin (1991)
The spherical inclusion with imperfect interface , Journal of applied Me-chanics, :pp. 444–449[32] L. Nazarenko (1996)
Elastic properties of materials with ellipsoidal pores , International Ap-plied Mechanics, :pp. 46–52[33] L. Dormieux, D. Kondo, & F.-J. Ulm (2006)
Microporomechanics , John Wiley & Sons[34] D. M. Kochmann & G. N. Venturini (2013)
Homogenized mechanical properties of auxeticcomposite materials in finite-strain elasticity , Smart Materials and Structures, :p.084004[35] L. Placidi, L. Greco, S. Bucci, E. Turco, & N. L. Rizzi (2016)
A second gradient formulationfor a 2d fabric sheet with inextensible fibres , Zeitschrift f¨ur angewandte Mathematik undPhysik, :p. 114[36] E. Turco, M. Golaszewski, I. Giorgio, & F. D’Annibale (2017)
Pantographic lattices withnon-orthogonal fibres: Experiments and their numerical simulations , Composites Part B:Engineering, :pp. 1–14[37] Y. Solyaev, S. Lurie, & A. Ustenko (2018)
Numerical modeling of a composite auxetic meta-materials using micro-dilatation theory , Continuum Mechanics and Thermodynamics, pp.1–9[38] G. Ganzosch, K. Hoschke, T. Lekszycki, I. Giorgio, E. Turco, & W. H. M¨uller (2018) , Technische Mechanik, :pp. 233–245[39] H. Yang, G. Ganzosch, I. Giorgio, & B. E. Abali (2018)
Material characterization and com-putations of a polymeric metamaterial with a pantographic substructure , Zeitschrift f¨urangewandte Mathematik und Physik, :p. 105[40] L. Hitzler, M. Merkel, W. Hall, & A. ¨Ochsner (2018)
A review of metal fabricated with laser-and powder-bed based additive manufacturing techniques: process, nomenclature, materi-als, achievable properties, and its utilization in the medical sector , Advanced EngineeringMaterials, :p. 1700658[41] X. Wang, J. A. Mu˜niz-Lerma, M. A. Shandiz, O. Sanchez-Mata, & M. Brochu (2019)
Crystallographic-orientation-dependent tensile behaviours of stainless steel 316l fabricatedby laser powder bed fusion , Materials Science and Engineering: A, :p. 138395[42] L. Hitzler, J. Hirsch, J. Tomas, M. Merkel, W. Hall, & A. ¨Ochsner (2019)
In-plane anisotropyof selective laser-melted stainless steel: The importance of the rotation angle increment nd the limitation window , Proceedings of the Institution of Mechanical Engineers, PartL: Journal of Materials: Design and Applications, :pp. 1419–1428[43] U. M¨uhlich, L. Zybell, & M. Kuna (2012) Estimation of material properties for linear elasticstrain gradient effective media , European Journal of Mechanics-A/Solids, :pp. 117–130[44] M. Laudato (2020)
Nonlinear phenomena in granular solids: Modeling and experiments , Developments and Novel Approaches in Nonlinear Solid Body Mechanics , pp. 179–189,Springer[45] I. Giorgio, M. De Angelo, E. Turco, & A. Misra (2019)
A Biot–Cosserat two-dimensionalelastic nonlinear model for a micromorphic medium , Continuum Mechanics and Thermo-dynamics, pp. 1–13[46] H. Altenbach, W. H. M¨uller, & B. E. Abali (eds.) (2019)
Higher Gradient Materials andRelated Generalized Continua , vol. 120 of
Advanced Structured Materials, (256 pages) ,Springer, Cham[47] W. M¨uller, W. Rickert, & E. Vilchevskaya (2020)
Thence the moment of momentum , ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f¨ur Angewandte Mathematikund Mechanik, :p. e202000117[48] C. Pideri & P. Seppecher (1997)
A homogenization result for elastic material reinforcedperiodically with high rigidity elastic fibres , Comptes Rendus de l’Academie des SciencesSeries IIB Mechanics Physics Chemistry Astronomy, :pp. 475–481[49] S. Forest, R. Dendievel, & G. R. Canova (1999)
Estimating the overall properties of heteroge-neous cosserat materials , Modelling and Simulation in Materials Science and Engineering, :p. 829[50] V. Kouznetsova, M. G. Geers, & W. M. Brekelmans (2002)
Multi-scale constitutive mod-elling of heterogeneous materials with a gradient-enhanced computational homogenizationscheme , International Journal for Numerical Methods in Engineering, :pp. 1235–1260[51] W. Pietraszkiewicz & V. Eremeyev (2009)
On natural strain measures of the non-linearmicropolar continuum , International Journal of Solids and Structures, :pp. 774–787[52] I. Giorgio, N. Rizzi, & E. Turco (2017)
Continuum modelling of pantographic sheets forout-of-plane bifurcation and vibrational analysis , Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, :p. 20170636[53] F. dell’Isola, I. Giorgio, M. Pawlikowski, & N. Rizzi (2016)
Large deformations of planarextensible beams and pantographic lattices: heuristic homogenization, experimental andnumerical examples of equilibrium , Proceedings of the Royal Society A: Mathematical,Physical and Engineering Sciences, :p. 20150790[54] G. Maurice, J.-F. Ganghoffer, & Y. Rahali (2019)
Second gradient homogenization of mul-tilayered composites based on the method of oscillating functions , Mathematics and Me-chanics of Solids, :pp. 2197–2230[55] J. Li (2011)
A micromechanics-based strain gradient damage model for fracture prediction ofbrittle materials–part i: Homogenization methodology and constitutive relations , Interna-tional Journal of Solids and Structures, :pp. 3336–3345[56] H. Askes & E. C. Aifantis (2011)
Gradient elasticity in statics and dynamics: an overview offormulations, length scale identification procedures, finite element implementations andnew results , International Journal of Solids and Structures, :pp. 1962–1990[57] J. R. Willis (1981)
Variational and related methods for the overall properties of composites , Advances in applied mechanics , vol. 21, pp. 1–78, Elsevier[58] P. Franciosi, M. Spagnuolo, & O. U. Salman (2018)
Mean green operators of deformable fibernetworks embedded in a compliant matrix and property estimates , Continuum Mechanicsand Thermodynamics, pp. 1–32[59] T.-H. Tran, V. Monchiet, & G. Bonnet (2012)
A micromechanics-based approach for thederivation of constitutive elastic coefficients of strain-gradient media , International Jour-18al of Solids and Structures, :pp. 783–792[60] Y. Rahali, I. Giorgio, J. Ganghoffer, & F. dell’Isola (2015)
Homogenization `a la piola producessecond gradient continuum models for linear pantographic lattices , International Journalof Engineering Science, :pp. 148–172[61] H. Reda, I. Goda, J. Ganghoffer, G. L’Hostis, & H. Lakiss (2017) Dynamical analysis of ho-mogenized second gradient anisotropic media for textile composite structures and analysisof size effects , Composite Structures, :pp. 540–551[62] J.-F. Ganghoffer, G. Maurice, & Y. Rahali (2019)
Determination of closed form expressionsof the second-gradient elastic moduli of multi-layer composites using the periodic unfoldingmethod , Mathematics and Mechanics of Solids, :pp. 1475–1502[63] Y. Rahali, M. Assidi, I. Goda, A. Zghal, & J.-F. Ganghoffer (2016)
Computation of theeffective mechanical properties including nonclassical moduli of 2.5 d and 3d interlocks bymicromechanical approaches , Composites Part B: Engineering, :pp. 194–212[64] K. ElNady, I. Goda, & J.-F. Ganghoffer (2016) Computation of the effective nonlinear me-chanical response of lattice materials considering geometrical nonlinearities , Computa-tional Mechanics, :pp. 957–979[65] Y. Rahali, F. Dos Reis, & J.-F. Ganghoffer (2017)
Multiscale homogenization schemes forthe construction of second-order grade anisotropic continuum media of architectured ma-terials , International Journal for Multiscale Computational Engineering, [66] B. E. Abali, H. Yang, & P. Papadopoulos (2019)
A computational approach for determi-nation of parameters in generalized mechanics , H. Altenbach, W. H. M¨uller, & B. E.Abali (Eds.)
Higher Gradient Materials and Related Generalized Continua , AdvancedStructured Materials, vol. 120, chap. 1, pp. 1–18, Springer, Cham[67] B. E. Abali & H. Yang (2020)
Parameter determination of metamaterials in generalized me-chanics as a result of computational homogenization , D. A. Indeitsev & A. M. Krivtsov(Eds.)
Advanced Problems in Mechanics. APM 2019 , Lecture Notes in Mechanical Engi-neering, chap. 2, pp. 22–31, Springer, Cham[68] A. Skrzat & V. A. Eremeyev (2020)
On the effective properties of foams in the framework ofthe couple stress theory , Continuum Mechanics and Thermodynamics, pp. 1–23[69] J.-J. Alibert, P. Seppecher, & F. dell’Isola (2003)
Truss modular beams with deformationenergy depending on higher displacement gradients , Mathematics and Mechanics of Solids, :pp. 51–73[70] J. Alibert & A. Della Corte (2015)
Second-gradient continua as homogenized limit of panto-graphic microstructured plates: a rigorous proof , Zeitschrift f¨ur angewandte Mathematikund Physik, :pp. 2855–2870[71] A. Bensoussan, J.-L. Lions, & G. Papanicolaou (1978)
Asymptotic Analysis for PeriodicStructures , North-Holland, Amsterdam[72] S. J. Hollister & N. Kikuchi (1992)
A comparison of homogenization and standard mechanicsanalyses for periodic porous composites , Computational Mechanics, :pp. 73–95[73] P. W. Chung, K. K. Tamma, & R. R. Namburu (2001)
Asymptotic expansion homogenizationfor heterogeneous media: computational issues and applications , Composites Part A:Applied Science and Manufacturing, :pp. 1291–1301[74] I. Temizer (2012)
On the asymptotic expansion treatment of two-scale finite thermoelasticity ,International Journal of Engineering Science, :pp. 74–84[75] S. Forest, F. Pradel, & K. Sab (2001) Asymptotic analysis of heterogeneous cosserat media ,International Journal of Solids and Structures, :pp. 4585–4608[76] V. A. Eremeyev (2016)
On effective properties of materials at the nano-and microscalesconsidering surface effects , Acta Mechanica, :pp. 29–42[77] J.-F. Ganghoffer, I. Goda, A. A. Novotny, R. Rahouadj, & J. Sokolowski (2018)
Homogenizedcouple stress model of optimal auxetic microstructures computed by topology optimiza-tion , ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f¨ur AngewandteMathematik und Mechanik, :pp. 696–717[78] E. Turco (2019)
How the properties of pantographic elementary lattices determine the prop- rties of pantographic metamaterials , B. Abali, H. Altenbach, F. dell’Isola, V. Eremeyev,& A. ¨Ochsner (Eds.) New Achievements in Continuum Mechanics and Thermodynamics ,vol. 108 of
Advanced Structured Materials , pp. 489–506, Springer[79] M. Peszynska & R. E. Showalter (2007)
Multiscale elliptic-parabolic systems for flow andtransport. , Electronic Journal of Differential Equations (EJDE), :pp. 1–30[80] J. Pinho-da Cruz, J. Oliveira, & F. Teixeira-Dias (2009)
Asymptotic homogenisation in linearelasticity. part i: Mathematical formulation and finite element modelling , ComputationalMaterials Science, :pp. 1073–1080[81] Y. Efendiev & T. Y. Hou (2009)
Multiscale finite element methods: theory and applications ,vol. 4, Springer Science & Business Media[82] C. Boutin (1996)
Microstructural effects in elastic composites , International Journal of Solidsand Structures, :pp. 1023–105[83] E. Barchiesi, F. dell’Isola, M. Laudato, L. Placidi, & P. Seppecher (2018)
A 1d continuummodel for beams with pantographic microstructure: Asymptotic micro-macro identificationand numerical results , Advances in Mechanics of Microstructured Media and Structures ,pp. 43–74, Springer[84] A. Bacigalupo (2014)
Second-order homogenization of periodic materials based on asymptoticapproximation of the strain energy: formulation and validity limits , Meccanica, :pp.1407–1425[85] C. Boutin, I. Giorgio, L. Placidi, et al. (2017)
Linear pantographic sheets: Asymptotic micro-macro models identification , Mathematics and Mechanics of Complex Systems, :pp.127–162[86] L. Placidi, U. Andreaus, A. Della Corte, & T. Lekszycki (2015)
Gedanken experiments for thedetermination of two-dimensional linear second gradient elasticity coefficients , Zeitschriftf¨ur angewandte Mathematik und Physik, :pp. 3699–3725[87] M. Cuomo, L. Contrafatto, & L. Greco (2014)
A variational model based on isogeometricinterpolation for the analysis of cracked bodies , International Journal of Engineering Sci-ence, :pp. 173–188[88] R. Peerlings & N. Fleck (2004) Computational evaluation of strain gradient elasticity con-stants , International Journal for Multiscale Computational Engineering, [89] J. Li (2011)
Establishment of strain gradient constitutive relations by homogenization ,Comptes Rendus M´ecanique, :pp. 235–244[90] J. Li & X.-B. Zhang (2013)
A numerical approach for the establishment of strain gradientconstitutive relations in periodic heterogeneous materials , European Journal of Mechanics-A/Solids, :pp. 70–85[91] S. Barboura & J. Li (2018) Establishment of strain gradient constitutive relations by usingasymptotic analysis and the finite element method for complex periodic microstructures ,International Journal of Solids and Structures, :pp. 60–76[92] M. M. Ameen, R. Peerlings, & M. Geers (2018)
A quantitative assessment of the scale sepa-ration limits of classical and higher-order asymptotic homogenization , European Journalof Mechanics-A/Solids, :pp. 89–100[93] A. Porubov & E. Grekova (2020) On nonlinear modeling of an acoustic metamaterial , Me-chanics Research Communications, :p. 103464[94] H. Yang, B. E. Abali, D. Timofeev, & W. H. M¨uller (2019)
Determination of metamate-rial parameters by means of a homogenization approach based on asymptotic analysis ,Continuum Mechanics and Thermodynamics, pp. 1–20[95] T. Tancogne-Dejean, N. Karathanasopoulos, & D. Mohr (2019)
Stiffness and strength ofhexachiral honeycomb-like metamaterials , Journal of Applied Mechanics, [96] Gnu Public (2007),
Gnu general public license
Fenics
Supply code for computations , http://bilenemek.abali.org/[99] S. Whitaker (1967)
Diffusion and dispersion in porous media , AIChE Journal, :pp.2020–427[100] J. C. Slattery (1967)
Flow of viscoelastic fluids through porous media , AIChE Journal, :pp. 1066–1071[101] W. G. Gray & P. Lee (1977)
On the theorems for local volume averaging of multiphasesystems , International Journal of Multiphase Flow, :pp. 333–340[102] T. I. Zohdi (2018)
Finite Element Primer for Beginners , Springer[103] F. Larsson, K. Runesson, S. Saroukhani, & R. Vafadari (2011)
Computational homogenizationbased on a weak format of micro-periodicity for RVE-problems , Computer Methods inApplied Mechanics and Engineering, :pp. 11–26[104] B. E. Abali (2017)