Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2024.49.0020 Acta Polytechnica CTU Proceedings 49:20–25, 2024 © 2024 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague EFFECT OF GEOMETRY ON HOMOGENISED PROPERTIES OF SELECTED AUXETIC METAMATERIALS Nataša Jošková∗, Martin Doškář Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague, Czech Republic ∗ corresponding author: joskonat@cvut.cz Abstract. Our study investigates the influence of geometrical parameters of two types of auxetic metamaterials on their effective properties. In particular, we focus on three-dimensional lattice structures, which we represent with discrete beam models of their respective Periodic Unit Cells (PUCs). Limiting the scope of the study to linear elasticity, we compute the effective response of PUCs by plugging the kinematic ansatz of the first-order numerical homogenisation into the strain energy expression arising from the Direct Stiffness Method and minimising the energy with respect to the periodic fluctuation field. The obtained effective stiffness matrices are post-processed to arrive at elastic parameters as Poisson’s ratios coefficients, that are reported in different directions with respect to the key geometrical parameters. Keywords: Auxetic metamaterial, first-order homogenisation, periodic unit cell, effective Poisson’s ratio. 1. Introduction Metamaterials are artificial materials with properties beyond those of materials found in nature. These prop- erties are mainly determined by their microstructure rather than by the chemical or physical parameters of the bulk constituents from which they are made [1]. Due to technological advances in recent decades, com- plex microstructures of these metamaterials can be produced by manufacturing techniques such as 3D and even 4D printing [2], optical lithography [3], or elec- trospinning. Mathematical modelling is then needed for efficient design of metamaterials by circumventing lenghty experimental search for their optimal design. Our study of the influence of geometry on the effec- tive Poisson’s ratio focuses on two variants of a three- dimensional auxetic metamaterial (cubic and hexago- nal) proposed by Bückmann et al. [3] and shown in Figure 1. Both designs exhibit a periodic microstruc- ture allowing us to investigate only the response of a Periodic Unit Cell (PUC) as their Representative Volume Element (RVE). 2. Geometry of investigated metamaterial The microstructure of metamaterial is composed of arranged bow-tie structures, the geometry of which is controlled by the angle δ located between the diagonal beam of the central bow-tie structure in the upper left quadrant and the yz plane passing through its initial node, see insets on the right-hand side in Figure 2. The 3D metamaterial is created by rotating the central bow-tie structure located in the xz plane around its vertical centre beams, as shown in Figure 2. The cubic microstructure is obtained by rotation by 90°, from (a). Cubic variant. (b). Hexagonal variant. Figure 1. Two investigated auxetic metamaterials. Insets on the right-hand side show a top view of the microstructures. 20 https://doi.org/10.14311/APP.2024.49.0020 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 49/2024 Effect of geometry on homogenised properties of auxetic metamaterial (a). Cubic microstructure. (b). Hexagonal prismatic microstructure. Figure 2. Periodic Unit Cell of cubic and regular hexagonal prismatic microstructure along with their corresponding side view including angle δ, which con- trols the geometry of the metamaterial. which the central structure is formed (dark blue). This central structure is complemented by a similar one (cyan) shifted by half a period in yz plane. To create the PUC of regular hexagonal prism with equal height and long base diagonal length, we perform 60° and 120° rotations of the bow-tie structure (purple). The height H of both PUCs will be considered as a single unit height. Here, we model PUC with discrete beams; we know the position of each node and the orientation of the beams that connect them. For this study, we assume a circular beam cross-section with diameter d = 0.1H. Consequently, cross-sectional characteristics follow as: A = π 4 · d2, I = π 64 · d4, J = 2I = π 32 · d4, (1) where A is the cross-section area, I is the second moment of inertia, and J is the polar moment of inertia. The volume of PUC is for the cubic variant Vc = H3 and Vh = 3 √ 3 8 · H3 for the hexagonal one. 3. Direct stiffness method To compute a mechanical response of the PUC model, we use a linear discrete beam model that can be described by a linear relation: F = Ku, (2) where the vector F contains the forces and moments ap- plied on all nodes, K is the stiffness matrix – assembled from individual submatrices Ki for each beam – and u is the displacement vector that successively contains subvectors of displacements ui, vi, wi and rotations φx,i, φy,i, φz,i of individual nodes i. 3.1. Local stiffness matrix The local stiffness matrix Kl i for the ith beam follows from the Bernoulli-Euler beam theory [4, 5] and has size 12 × 12, due to the 6 unknowns located at the beginning (index b) and the end (index e) node of the beam. Individual parts that contribute to the stiffness matrix can be divided into 4 submatrices pertinent to: (1) membrane behaviour:[ X l b X l e ] = EA Li [ 1 −1 −1 1 ] [ ul b ul e ] , (3) (2) bending in xy plane:  Y l b M l z,b Y l e M l z,e  = 2EIz Li  6 L2 i 3 Li − 6 L2 i 3 Li 3 Li 2 − 3 Li 1 − 6 L2 i − 3 Li 6 L2 i − 3 Li 3 Li 1 − 3 Li 2   vl b φl z,b vl e φl z,e  , (4) (3) bending in xz plane:  Zl b M l y,b Zl e M l y,e  = 2EIy Li  6 L2 i − 3 Li − 6 L2 i − 3 Li − 3 Li 2 3 Li 1 − 6 L2 i 3 Li 6 L2 i 3 Li − 3 Li 1 3 Li 2   wl b φl y,b wl e φl y,e  , (5) (4) torsion:[ M l x,b M l x,e ] = GJ Li [ 1 −1 −1 1 ] [ φl x,b φl x,e ] . (6) In the equations above, Li refers to the length of a corresponding beam, E denotes Young’s modulus, G stands for the shear modulus, and moments of inertia Iy, Iz are equal to I due to the circular cross- section of the beam. The local stiffness matrix Kl i is obtained by combining the submatrices from Equa- tions (3)–(6), each contributing to its specific degrees of freedom (DOFs). 21 Nataša Jošková, Martin Doškář Acta Polytechnica CTU Proceedings 3.2. Global stiffness matrix The beams constituting the PUC’s microstructure have different orientations. Hence, it is necessary to transform local displacements, rotations and end forces from Equations (3)–(6) into a global coordinate system. In our study, we took the approach of building a rotation matrix Ri for the ith beam with Euler angles. The beam’s initial position is established by align- ing its local coordinate system with the global one. Afterwards, we execute an extrinsic rotation around the global coordinate axes until the desired position is achieved. Rotation matrices: Rx = 1 0 0 0 cos α − sin α 0 sin α cos α  , (7) Ry =  cos β 0 sin β 0 1 0 − sin β 0 cos β  , (8) Rz = cos γ − sin γ 0 sin γ cos γ 0 0 0 1  , (9) determine rotation around global x, y, z axes by an- gles α, β, γ, respectively, while each angle stands for rotation from the latest beam’s position. To obtain the rotation matrix Ri, we perform matrix multiplication: Ri = Rz · Ry · Rx, (10) where the sequence of the elements is based on order in which the beam is rotated, starting with rotation around the global x axis. The nodal displacements and rotations are trans- formed on each side of the beam equally, so we can create a transformation matrix Ti by arranging the rotation matrix Ri on its diagonal: Ti =  Ri 0 0 0 0 Ri 0 0 0 0 Ri 0 0 0 0 Ri  . (11) The local stiffness matrix is then transformed into a global coordinate system using the relation: Kg i = TT i Kl iTi, (12) from which we obtain the global stiffness matrix Kg i for each beam and localise them to the stiffness matrix K using Boolean localisation matrices Li: K = ∑ i LT i Kg i Li. (13) 4. Homogenisation The homogenisation process substitutes a heteroge- neous PUC at the microscopic level with a corre- sponding macroscopic constitutive model, allowing us to study the effective behaviour of the PUC when treated as a material considering its microstructure. In our study, we are using the first-order numerical homogenisation to obtain effective metamaterial prop- erties. 4.1. Displacement decomposition In the first-order homogenisation, the total displace- ment field u⃗(x⃗) is assumed in the form: u⃗(x⃗) = u⃗ E(x⃗) + u⃗ ∗(x⃗), (14) where u⃗ E denotes the macroscopic and u⃗∗ is the fluc- tuation part of the displacement field caused by the heterogeneity of the metamaterial [6]. The macro- scopic part u⃗ E of the displacement field corresponds to a situation under which an entire cell composed of homogeneous material would be subjected to a con- stant macroscopic strain tensor E: E = Exx Exy Exz Eyx Eyy Eyz Ezx Ezy Ezz  , (15) which results in a displacement field uE given as: u⃗ E(xi) = E · xi = u E i v E i w E i  . (16) For the first-order homogenisation, it is further assumed that the volumetric average of the gradient of the entire displacement field u⃗(x⃗) corresponds to the prescribed macroscopic deformation E. By applying the symmetric gradient operator ∇s = 1 2 (∇ + ∇T) and averaging the result over a unit cell Ω, we obtain: E = 1 |Ω| ∫ Ω ∇su⃗(x⃗) dx⃗ = 1 |Ω| ∫ Ω ∇su⃗ E(x⃗) + ∇su⃗ ∗(x⃗) dx⃗, (17) where Ω represents the microscale domain of interest, being the metamaterial’s macroscopic point [7]. The macroscopic part of the deformation u⃗ E(x⃗) in Equation (16) is defined such that: E = 1 |Ω| ∫ Ω ∇su⃗ E(x⃗)dx⃗. (18) Consequently, the fluctuation part u⃗ ∗(x⃗) of the dis- placement field u⃗(x⃗), must have a zero volumetric average gradient, i.e.: 1 |Ω| ∫ Ω ∇su⃗ ∗(x⃗)dx⃗ = 0. (19) 22 vol. 49/2024 Effect of geometry on homogenised properties of auxetic metamaterial For our discrete beam model, Equation (16) in a ma- trix form reads as: uE i =  xi 0 0 0 1 2 zi 1 2 yi 0 yi 0 1 2 zi 0 1 2 xi 0 0 zi 1 2 yi 1 2 xi 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  E = QE i E, (20) and couples nodal DOFs with the macroscopic de- formation E, which is the vectorial representation of a symmetric second-order tensor: E = [ Ex Ey Ez Γyz Γxz Γxy ]T . (21) Now we can write the original degrees of freedom u depending on the macroscopic deformation E and the fluctuation unknowns u∗. Let’s define the extended displacement vector: û = [ E u∗ ] , (22) then we can express degrees of freedom of our discrete beam model as: u = [ QE I ] [ E u∗ ] = Q û, (23) where I is the square identity matrix and QE is com- posed of the blocks QE i corresponding to the expres- sion (20). 4.2. Periodic boundary conditions To satisfy the constraint (19), we introduce Periodic Boundary Conditions (PBC), which is a natural model assumption for materials with periodic microstruc- tures. Let’s denote Π(x) the mapping from the source part Γs of the boundary Γ onto its periodic image. The fluctuation DOFs at a periodic point, denoted as u∗(Π(x)), are then equivalent to the corresponding DOFs u∗(x) at the boundary of the unit cell Γs: u∗(Π(x)) = u∗(x) ∀ x ∈ Γs. (24) To this end, we established a new vector a, contain- ing unknown periodic fluctuation DOFs, which maps to u∗ through the Boolean matrix P∗: u∗ = P∗a. (25) This step significantly reduces the number of un- knowns in our equations. To prevent the PUC from moving as a whole unit during deformation, we prescribe zero fluctuation dis- placements at the node in the centre of the top edge of the PUC, illustrated with a black node in Figure 3, while rotations remain free. Due to the periodicity, the fluctuation displacements must also vanish at the node in the centre of the bottom face of the PUC. Figure 3. Fixed fluctuation displacements at nodes in the centre of the top and bottom face of the PUC, marked down with the black dots. Similarly to Equation (23), the matrix P̂ connects the macroscopic deformation E and the fluctuation unknowns a, to the extended DOFs û: û = [ E u∗ ] = [ I P∗] [ E a ] = P̂ â. (26) By connecting the macroscopic and fluctuation parts of the displacement field from Equation (26) and considering Equation (23) we get an expression for full-field displacement DOFs u based on unknown â: u = QP̂ â. (27) 4.3. Energy minimisation For every macroscopic deformation E, there is a certain state into which the cell deforms, because it naturally attempts to reach the state with the lowest energy. The energy E of a discrete beam model, which represents our PUC composed of beams and nodes, can be written as: E = 1 2uTKu. (28) Similarly, for a linear elastic material of volume V that is subjected to a uniform deformation E at the macroscopic level, the relation for the energy can be expressed as: EM = V 1 2ETDhomE, (29) where Dhom is the material stiffness matrix of the desired homogenised metamaterial and vector E con- tains the individual macroscopic components of the deformation tensor in the vectorial form introduced in Equation (21). Plugging the unknowns from Equation (27) into Equation (28) yields an expression for PUC’s energy based on macroscopic deformations E and fluctuation unknowns a: E(E, a) = 1 2 [ E a ]T [ K̂EE K̂Ea K̂aE K̂aa ] [ E a ] = 1 2 âTK̂ â, (30) where stiffness matrix K̂, and its four submatrices K̂ij pertinent to E and a, follow from: K̂ = P̂TQTK Q P̂. (31) 23 Nataša Jošková, Martin Doškář Acta Polytechnica CTU Proceedings Since we are interested in the response of the ho- mogenised PUC to a prescribed macroscopic defor- mation E and not in displacement of its individual nodes, we express a with respect to the macroscopic deformation E of the cell. To this end, we keep the macroscopic deformation E fixed and determine the fluctuation displacements a as the solution ã(E) that minimizes the energy E(E, a) for the given deformation E as: ã(E) = argmin a∈Rn 1 2(ETK̂EEE + ETK̂Eaa + aTK̂aEE + aTK̂aaa), (32) where n represents the number of unknown fluctuation DOFs. As a result of the matrix K̂ being both sym- metric and positive definite, the quadratic form (30) attains a global minimum at the point of its zero gradient: ∇aE(E, a) ∣∣ a=̃a(E) = KaEE + Kaaã(E) = 0. (33) This provides us with the expression for the minimizer: ã(E) = −K−1 aa KaEE. (34) Substituting the expression above into Equation (30) yields an energy dependent entirely on the macroscopic deformations E: Ẽ(E) = E(E, ã(E)), (35) Ẽ(E) = 1 2 ET(K̂EE − K̂EaK̂−1 aa K̂aE) E = 1 2 ETKeffE. (36) Comparing the expression (36) with the formula for the energy of a homogeneous material of volume V = |Ω| subjected to a constant deformation E in the form in Equation (29), we arrive at the relation for the homogenised material stiffness of the auxetic metamaterial as: Dhom = 1 V Keff. (37) 4.4. Effective Poisson’s ratio The procedure introduced in the previous sections yields the whole effective stiffness matrix. However, comparing and discussing the entire stiffness matrix is cumbersome. Here, we focus on the homogenised Poisson’s effect as it is the primal objective of the auxetic metamaterial. Because our structures are symmetric in three mu- tually perpendicular directions, we can expect the overall orthotropic response. Consequently, there are three different values of the Poisson’s ratio depending on the direction of the prescribed relative deformation. Clearly, the Poisson’s ratios ν in the xy and xz planes, νxy and νxz, will be the same due to the symmetry of PUC. Furthermore, νyz and νzy should be equal for the same reason. We will refer to the stress in the metamaterial Σ to distinguish the macroscopic level from the microscopic one. The stress-strain relation is linear because we work in the range of Hooke’s law, so macroscopic stress can be written as: Σ = DhomE, (38) which can be broken down into individual components as follows: Σx Σy Σz Σyz Σxz Σxy  =  Dxx Dxy Dxz 0 0 0 Dyx Dyy Dyz 0 0 0 Dzx Dzy Dzz 0 0 0 0 0 0 Gyz 0 0 0 0 0 0 Gxz 0 0 0 0 0 0 Gxy   Ex Ey Ez Γyz Γxz Γxy  . (39) To determine the Poisson’s ratio νij , we perform a virtual uniaxial tension/compression experiment in which we prescribe the macroscopic strain in the ith di- rection and compute the strain in the jth direction for the requirement of zero macroscopic stress in the jth and kth direction. This experiment results in the following relation: ẼEi j = DikDjk − DijDkk DjjDkk − D2 jk Ei. (40) Using (40) we can express Poisson’s ratio for any di- rection, following its definition as a negative ratio of the derived and prescribed macroscopic strain. Con- sequently, the effective Poisson’s ratio of an auxetic metamaterial is given by: νij = − ẼEi j Ei = DikDjk − DijDkk D2 jk − DjjDkk . (41) 5. Results We parameterised the two microstructural geome- tries (cubic and hexagonal) from Section 2, recall Figures 1 and 2, with an angle δ ∈ (0°, 45°⟩. This range was chosen to avoid beams’ overlaps. The nu- merical results comply with our assumptions of equal Poisson’s ratio values in following directions: νxy = νxz, νyx = νzx, νyz = νzy, (42) see also overlapping lines in Figure 4. We observed auxetic behaviour in the whole range of δ with cu- bic PUC, while the metamaterial with hexagonal PUC exhibits pure auxetic properties only for an- gles δ ∈ ⟨10.56°, 45°⟩. For the lower values of δ, the metamaterial is auxetic only in the xy and xz planes, with the highest Poisson’s ratio value of νyz = 0.33, which leads to the lateral contraction during stretching in the yz plane. Poisson’s ratios νyx and νzx are the most influenced by the metamaterial’s geometry and attain their min- imum, within the investigated range of δ, from which 24 vol. 49/2024 Effect of geometry on homogenised properties of auxetic metamaterial 0 5 10 15 20 25 30 35 40 45 -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 xy xz yx zx yz zy (a). Cubic microstructure. 0 5 10 15 20 25 30 35 40 45 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 xy xz yx zx yz zy (b). Hexagonal prismatic microstructure. Figure 4. Poisson’s ratio as a function of angle δ for cubic and regular hexagonal prismatic microstructure. their value starts increasing and slowly approaching remaining Poisson’s ratios. Specifically, for the cubic PUC, the global minimum νyx = −1.98 is achieved by the geometry of angle δ = 13.78°. Hexagonal PUC exhibits its minimal Poisson’s ratio of value νyx = −1.55 by the angle δ = 15.92°. The remaining pairs of Poisson’s ratios, i.e. νxy, νxz and νyz, νzy, on the other hand, tend to decrease with the increasing angle δ across the entire range. Additionally, we observe certain values of δ, where some Poisson’s ratios coincide. In the cubic PUC, νxy and νyz attain the same value of −0.05 for δ = 3.32°. The hexagonal PUC exhibits two such values of δ, δ = 16.96° and δ = 36.98°, for which Pois- son’s ratios are equal to νxy = νyz = −0.26 and −0.63, respectively. In addition, the hexagonal PUC have equal Pois- son’s ratios νxy and νyx of value νxy = −0.82 at the angle δ = 44.71°. 6. Conclusion We investigated the influence of the geometry con- trolled by the angle δ in the bow-tie part of the mi- crostructure on the effective Poisson’s ratios of two three-dimensional metamaterials. Their PUCs were modelled with discrete beam elements, and the effec- tive metamaterial properties were determined using Direct Stiffness Method and the first-order numerical homogenisation, resulting in a connection between their microstructure and macroscopic behaviour. Given the symmetries of both investigated mi- crostructural geometries, we obtain three distinct val- ues of Poisson’s ratios as functions of angle δ. Two of these values are monotonously decreasing with in- creasing angle δ, while the remaining value (same for νyx and νzx) exhibits a minimum within the studied range of ν. The cubic PUC features the minimum value νyx = −1.98 for δ = 13.78°, while the hexagonal PUC attains the minimum νyx = −1.55 for δ = 15.92°. In conclusion, only the cubic PUC delivers auxetic behaviour in all directions for all investigated values of δ. The hexagonal PUC shares the same trait only for δ ≥ 10.56°. Acknowledgements This work was supported by the Grant Agency of the Czech Technical University in Prague, grant No. SGS23/032/OHK1/1T/11. References [1] E. Barchiesi, M. Spagnuolo, L. Placidi. Mechanical metamaterials: A state of the art. Mathematics and Mechanics of Solids 24(1):212–234, 2019. https://doi.org/10.1177/1081286517735695 [2] X. Zhou, L. Ren, Z. Song, et al. Advances in 3D/4D printing of mechanical metamaterials: From manufacturing to applications. Composites Part B: Engineering 254:110585, 2023. https: //doi.org/10.1016/j.compositesb.2023.110585 [3] T. Bückmann, N. Stenger, M. Kadic, et al. Tailored 3D mechanical metamaterials made by dip-in direct-laser-writing optical lithography. Advanced Materials 24(20):2710–2714, 2012. https://doi.org/10.1002/adma.201200584 [4] A. Kassimali. Matrix analysis of structures. Cengage Learning, Stamford, Australia, 2nd edn., 2012. ISBN 978-1-111-42620-0. [5] W. McGuire, R. H. Gallagher, R. D. Ziemian. Matrix structural analysis. John Wiley, New York, 2nd edn., 2000. ISBN 9781507585139. [6] J. C. Michel, H. Moulinec, P. Suquet. Effective properties of composite materials with periodic microstructure: A computational approach. Computer Methods in Applied Mechanics and Engineering 172(1–4):109–143, 1999. https://doi.org/10.1016/S0045-7825(98)00227-8 [7] M. Doškář, J. Novák. A jigsaw puzzle framework for homogenization of high porosity foams. Computers & Structures 166:33–41, 2016. https://doi.org/10.1016/j.compstruc.2016.01.003 25 https://doi.org/10.1177/1081286517735695 https://doi.org/10.1016/j.compositesb.2023.110585 https://doi.org/10.1016/j.compositesb.2023.110585 https://doi.org/10.1002/adma.201200584 https://doi.org/10.1016/S0045-7825(98)00227-8 https://doi.org/10.1016/j.compstruc.2016.01.003 Acta Polytechnica CTU Proceedings 49:20–25, 2024 1 Introduction 2 Geometry of investigated metamaterial 3 Direct stiffness method 3.1 Local stiffness matrix 3.2 Global stiffness matrix 4 Homogenisation 4.1 Displacement decomposition 4.2 Periodic boundary conditions 4.3 Energy minimisation 4.4 Effective Poisson's ratio 5 Results 6 Conclusion Acknowledgements References