Issue 
EPJ Applied Metamaterials
Volume 2, 2015
Advanced Metamaterials in Microwaves, Optics and Mechanics



Article Number  13  
Number of page(s)  14  
DOI  https://doi.org/10.1051/epjam/2015013  
Published online  26 January 2016 
https://doi.org/10.1051/epjam/2015013
Research Article
Homogenization of 1D and 2D magnetoelastic lattices
^{1}
George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, Georgia
30332, United States
^{2}
Daniel Guggenheim School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, Georgia
30332, United States
^{*} email: mschaeffer3@gatech.edu
Received:
25
September
2015
Accepted:
17
November
2015
Published online: 26 January 2016
This paper investigates the equivalent inplane mechanical properties of one dimensional (1D) and two dimensional (2D), periodic magnetoelastic lattices. A lumped parameter model describes the lattices using magnetic dipole moments in combination with axial and torsional springs. The homogenization procedure is applied to systems linearized about stable configurations, which are identified by minimizing potential energy. Simple algebraic expressions are derived for the properties of 1D structures. Results for 1D lattices show that a variety of stiffness changes are possible through reconfiguration, and that magnetization can either stiffen or soften a structure. Results for 2D hexagonal and reentrant lattices show that both reconfigurations and magnetization have drastic effects on the mechanical properties of lattice structures. Lattices can be stiffened or softened and the Poisson’s ratio can be tuned. Furthermore for certain hexagonal lattices the sign of Poisson’s ratio can change by varying the lattice magnetization. In some cases presented, analytical and numerically estimated equivalent properties are validated through numerical simulations that also illustrate the unique characteristics of the investigated configurations.
Key words: Magnetoelastic / Homogenization / Tunable properties / Reconfiguration
© M. Schaeffer and M. Ruzzene, published by EDP Sciences, 2015
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
In previous studies [1, 2], the authors have shown that periodic magnetoelastic structures can be reconfigured and tuned to change the way that plane waves travel through them. Though it is evident that such reconfiguration and tuning also affects the mechanical properties of the structures, these changes and the quantification thereof have not yet been studied in detail.
Assembly of conventional materials in novel ways creates new structures that can perform better than their constituents. Such structures are often composed of “metamaterials” and are called “metastructures”. Periodic and cellular structures are common products of such an endeavor as they can be extremely stiff for their weight, or exhibit novel wave propagation properties [3–8]. Oftentimes, it is convenient to investigate periodic metamaterial assemblies as equivalent continuua through homogenization. The process of homogenization has been investigated for discrete periodic structures by, for example Born [9], Gonella and Ruzzene [10], Maradudin et al. [11], Suiker et al. [12]. Similar methods are applied to find the equivalent properties of the discrete models discussed in this publication.
Related work has been presented [13] in which finite magnetoelastic structures show that their Poisson’s ratio can be tuned with the application of a nonuniform external magnetic field and through changes in constituent dimensions. The magnetic dipole moments in the 2D structures discussed therein are oriented in the plane of the system and the changes in Poisson’s ratio are a result of the reorientation of the system’s rigid entities under the constant application of the external magnetic field. Also related is the homogenization methodology in reference [10], which is well suited for calculating equivalent mechanical properties of periodic lattices and is applied to hexagonal and reentrant beam lattices. It is also generalized for arbitrarily complex elastic structures, but is not presented for use in systems exhibiting nonnearest neighbor interactions, such as magnetic interactions, and must therefore be extended for application to magnetoelastic lattices. In contrast to the work by Grima et al. [13], the structures discussed in this publication contain magnetic dipole moments perpendicular to the lattice plane, and mechanical properties including stiffnesses are presented. The changes Poisson’s ratio in this publication are not brought about solely by geometric orientation as in reference [13], but also through changes in interaction stiffnesses, which are caused by the included magnets. Furthermore, this publication focus specifically on the multistability and bistability of the studied lattices and how reconfiguration effects properties.
The objective of this paper is to extend the homogenization method presented in reference [10] to structures with nonnearest neighbor interactions and quantify the tunability of the equivalent mechanical properties of reconfigurable periodic magnetoelastic lattices such as those studied in reference [1]. This includes stiffness and density for 1D lattices and Poisson’s ratios and stiffnesses for 2D lattices.
The paper is organized as follows. Following this introduction the modeling approach used is discussed in Section 2, which contains the methodology for the study of periodic lattices, including the identification of equilibrium configurations (Sect. 2.4). Then the homogenization procedure applied is described in Section 3. Finally the homogenization method is applied to magnetoelastic structures to obtain solutions for 1D and 2D structures in Sections 4 and 5 respectively. Section 4 includes a comparison to numerical simulations for validation of the homogenization method. Conclusions are then summarized in Section 6.
2 Theoretical background
2.1 Lattice model
Magnetoelastic lattices are modeled as systems of permanent magnetic particles with translational and rotational degrees of freedom (DOFs). Axial and torsional springs connect the particles to form an elastic structure. The particles have a finite radius r_{p}, and feature both translational and rotational inertias. All particles are identical in terms of their inertial and magnetic properties. Similarly, all springs connecting the particles are massless and feature the same axial and torsional stiffness. The interactions between particles are governed by: (1) magnetic interactions, and (2) elastic interactions through axial and torsional springs. The energy functionals associated with the interactions above and the governing equations are described in the following subsections. Details can be found in reference [1].
2.2 Energy functionals
The behavior of the generic ith particle of a 2D lattice is described by the position of its center of mass r_{i} = x_{i}i + y_{i}j, and by a rotation angle θ_{i}. Accordingly, the vector of the generalized coordinates associated with the ith particle is given by:(1)
2.2.1 Mechanical interactions
The strain energy associated with the elastic interactions is defined in terms of axial and torsional contributions. A phenomenological, simplified model is formulated whereby the axial spring is characterized by a bilinear forcedisplacement relation, which is chosen as a simplified model for a contact condition between two particles. The contact stiffness is described as proportional to the stiffness k_{a} through a parameter α > 1. A contact distance r_{c} = 2r_{p} defines the contact conditions. When the contact conditions are satisfied the axial stiffness of the bilinear spring is αk_{a}. The strain energy associated with the mechanical components connecting particles i and j is therefore given by:(2)where k_{τ} denotes the torsional spring constant, r _{ij} = r _{j} − r _{i}, r_{0} defines the interparticle distance when the axial springs are undeformed, and H is the Heaviside unit step function.
2.2.2 Magnetic interactions
The magnetic interaction force is described by modeling the permanent magnets as magnetically rigid dipole moments [14]. Furthermore, all magnetic dipole moments are considered perpendicular to the direction of motion. The force between two dipoles is given by Yung et al. [15]. In the 1D and 2D lattices considered, the motion of the masses is constrained to the direction i and to the lattice plane i , j respectively. This simplifies the expression of the magnetic interaction force as follows:(3)where μ_{0} is the permeability of free space, m _{i} and m _{j} are the dipole moments of particles i and j, n = r _{ij}/r_{ij} is a unit vector and ψ ^{(m)} = m_{i} m_{j}3μ_{0}/(4π) contains all the magnetic constants. Evaluating the work of the magnetic force (Eq. (3)) from a reference state at r_{ij} → ∞ to the current state r_{ij} defines a magnetic potential energy function .
The magnetic force described by equation (3) in theory acts for r → ∞. This implies that the lattice is, in general, characterized by nonnearest neighbor interactions that extend to infinity. However, the magnitude of the magnetic force rapidly decays with interparticle distance, which for practical purposes allows the introduction of a radius of influence r_{∞}, beyond which magnetic interactions are neglected. As many as 3,000 particles may be included in a typical analysis of a magnetized structure.
2.3 Lattice description and unit cell energy
The considered lattices are characterized by periodic distributions of masses along the i direction in the case of the 1D lattices, and over the plane i, j for the 2D lattices. Their positions are described through a set of lattice vectors di for the 1D lattices and d_{1}, d_{2} for the 2D lattices.
A unit cell is here defined as the smallest structure which can be repeated to describe the mass, stiffness, and magnetic characteristics of the assembly. The total energy of the unit cell is given by the sum of the potential energy U and kinetic energy T and the work of the external forces W:(4)where e = 1, …, E identifies the particles belonging to a unit cell. The magnetic interactions can be considered as external forces or, since they are magnetically rigid, they can be included in U_{e}. As the lattice is periodic, this energy functional is equal for all unit cells and its expression does not need to keep track of a specific unit cell identifier, which simplifies notation. The generalized coordinates for the unit cell are therefore:(5)
The set of governing equations describing the behavior of the cell and its interactions with connected units can be derived from the energy functional through application of Lagrange’s equations. The resulting set of governing equations is:(6)where are unit cell level matrices, and are unit cell level vectors of generalized coordinates and magnetic forces.
2.4 Identification of equilibrium configurations
The energy functional in equation (4), restricted to the time constant terms, is employed for the identification of the periodic equilibrium configurations and for the analysis of their stability which is discussed in this section. The stable periodic equilibrium configurations are identified through an optimization search algorithm that seeks for local energy minima [16]. The cell interacting with a reference unit are denoted with the index n = −N, …, N in 1D lattices, and the index pair n = −N, …, N; m = −M, …, M in 2D lattices.
In order to reduce the number of possible combinations only configurations consisting of four or less particles are considered in 1D lattices, and only configurations consisting of two particles are considered in 2D lattices. This is enforced by imposing an Eparticle periodicity (where E = 1, 2, or 4) as a constraint to the optimization problem defining the equilibrium conditions. This is practically implemented by imposing periodic conditions prior to the minimization of the static energy functional, i.e. by letting:(7)with n = −N, …, N; m = −M, …, M, and where for simplicity q = q_{0,0} denotes the DOFs of the reference unit. This implies that the unknowns of the nonlinear algebraic problem seeking for a minimum of the energy functional are: the unit cell DOFs q including the coordinates x, y of the particles within a unit cell and their rotation θ, and the lattice vectors d_{1}, d_{2} in the case of 2D lattices. The resulting nonlinear system of equations is expressed as:(8)where x = [q, d_{1}, d_{2}]. The problem in equation (8) is solved through an unconstrained QuasiNewton optimization scheme as described by Rao [17].
Results of this process are shown in Figures 1 and 2 for 1D and 2D lattices respectively. All 1D configurations pictured in Figure 1 include periodic systems with 1–4 particle unit cells. All 2D configurations to be discussed are pictured in Figure 2. Additional examples of equilibrium configurations are found in reference [1].
Figure 1. Periodic equilibrium examples for 1D structures with one (a), two (b–d), and four (e–g) particle unit cells. Red and blue circles denote magnets with dipole moments perpendicular to the lattice plane pointing towards and away from the reader respectively. Black lines denote elastic connections between magnets. Reprinted from M. Schaeffer and M. Ruzzene, Wave propagation in multistable magnetoelastic lattices, Int. J. Solids Struct. 56–57 (2015) 78–95, copyright (2015), with permission from Elsevier. 
Figure 2. Periodic equilibrium examples for 2D structures with two particle unit cells. Red and blue circles denote magnets with dipole moments perpendicular to the lattice plane pointing towards and away from the reader respectively. Black lines denote elastic connections between magnets. Particles in the unit cell have been highlighted by placing yellow dots in their centers. Reprinted from M. Schaeffer and M. Ruzzene, Wave propagation in multistable magnetoelastic lattices, Int. J. Solids Struct. 56–57 (2015) 78–95, copyright (2015), with permission from Elsevier. 
3 Homogenization procedure
3.1 Linearization about equilibrium
The considered homogenization approach assumes a linear elastic continuum. This requires linearization of the governing equations under the assumption that the lattice undergoes small displacements from a stable equilibrium.
A first order Taylor series expansion is employed to linearize the forces. Thus, the gradient of the forces with respect to the unit cell DOFs defines a stiffness matrix, and the equations of motion for a unit cell can be expressed as:(9)where q_{m,n} now denotes the displacement from equilibrium of the DOFs in cell m, n, and the absence of external forcing is assumed. In equation (9), K_{m,n} describes both the linearized elastic and magnetic interactions between the considered unit cell, here identified by the subscripts m = 0, n = 0, and the (2N + 1)(2M + 1) − 1 connected neighboring cells.
3.2 Reduction of the degrees of freedom
The homogenization method, developed by Gonella and Ruzzene [10], derives an equation of motion from the lattice structure that is compared to the continuum equations of motion in order to extract the mechanical properties. Before the method is applied the lattice must be reduced to an equivalent system with only one particle per unit cell, which ensures that the spacing between the particles makes a square grid in the lattice coordinate space. This is required to achieve cancellation between the 1st order partial derivatives that would otherwise manifest in the equation of motion derived by the homogenization method and break the analogy to the continuum equation of motion. Guyan reduction [18] accomplishes this task. Through the algebraic manipulation and substitution in Guyan reduction the n_{m} DOFs of the master particles m are used to define the n_{s} DOFs of the slave particles s. Any particle in a unit cell can be chosen as the master particle as long as the same particle is chosen in each unit cell. Partitioning the equation of motion into master and slave DOFs gives:(10)where n_{dof} = n_{m} + n_{s} is the number of DOFs describing the particles within the distance r_{∞} from a reference unit cell.
Though the equilibrium geometry of the structure is periodic, the deformation is arbitrary and cannot be described by the reference unit cell alone. Therefore, equation (10) is the equation of motion for an entire structure, not just one cell as equation (9), and [q_{m}; q_{s}] contains the q_{m,n} corresponding to all m, n pairs. Having equation (10), a transformation matrix is then calculated, which, when applied, solves for the master DOFs in terms of the slave DOFs(11)where I is the identity matrix. Applying the transformation matrix produces the reduced mass and stiffness matrices , and the displacement vector .(12)
In final preparation for homogenization methodology the mass and stiffness matrices are partitioned into blocks, each of which describe the effect of cell m, n’s DOFs on the reference cell’s DOFs. For the lattices studied, equation (13) separates out three equations from the n_{m} reduced equations of motion, which correspond to the three DOFs of the master node in the reference unit cell.(13)
It is worth noting here that for a static problem Guyan reduction introduces no assumptions, and thus no error regarding static elastic properties.
3.3 Longwavelength approximation
The homogenization method follows the procedure described by Gonella and Ruzzene [10] and is summarized here. First the lattice coordinate space η_{1}, η_{2} is defined using the lattice vectors d_{1}, d_{2}:(14) (15)
Then it is assumed that the displacement of the particles is continuous in space and can be represented by a second order Taylor series in the lattice coordinate space η_{1}, η_{2}. Thus, the master DOFs of the cell md_{1} + nd_{2} away from the reference unit cell may be approximated as:(16)
Application of equation (16) to equation (13) produces:(17)where, although the absence of external forcing will be assumed, a force vector f continuous in x, y is included for generality. Equation (17) assumes that for m, n ≠ 0, meaning there is no inertial coupling between cells. Of note is the fact that the matrices in equation (17) are dependent not on nearest neighbors only but on the cells up to M or N indexes away, where M, N ≠ 1 in general, allowing the inclusion of nonnearest neighbor magnetic interactions. The matrices in equation (17) are defined:(18)
Equation (18) differs from that in reference [10] as the typo omitting the 1/2 from A′ and C′ has been fixed. In order to convert to the physical space x, y from the lattice coordinate space η_{1}, η_{2} the derivatives are converted using:(19) (20)
Note that these equations differ from those in reference [10] slightly as two small typos have been fixed. The change of basis gives the continuum equation of motion in the physical space x, y:(21)where A, …, F, F_{m} are defined in terms of A′, …, F′, F′_{m} and E as follows:(22)
As a comparison will be made with the equation of motion for a continuum the equation must be expressed per unit volume V, which for a 2D structure is:(23)where b is the outofplane thickness of the structure. The new equation of motion is:(24)where b is a vector of body forces.
From here the procedure of [10] is continued, omitting terms that decrease with the length scale of the micro structure ε ^{3}. Then the rotational DOF θ is algebraically removed (assuming no moments are applied) and the two resulting equations are:(25)where q_{x}, q_{y} and b_{x}, b_{y} are the components of q and b respectively, and a … f are constants. Equation (25) is compared with the equations of motion for a thin orthotropic material (plane stress) which are:(26)
Comparing like terms allows algebraic solution of the density ρ the Young’s moduli E_{1} and E_{2}, the Poisson’s ratios ν_{12} and ν_{21}, and the shear modulus G_{12} in terms of a … f.
4 Homogenization of 1D magnetoelastic lattices
The procedure is first applied to the 1D lattices shown in Figures 1a and 1b. The analytical expressions of the equivalent properties of both structures derived below are used to check the convergence of the approximate solutions described in the previous section.
For both lattices the same approach is used. A variable s, used to denote the pattern type, allowing specification between the two lattices mathematically. It is −1 for alternating polarities (Figure 1b) and 1 if all polarities are the same (Figure 1a). First the equation of motion for the generic particle i is derived. It is assumed that the lattice is in equilibrium, and evenly spaced with distance h between all nearestneighbor magnets. The equation of motion is:(27)where u_{i} is the displacement of particle i in the x direction, and is the equivalent linearized stiffness of the magnetic interaction with a particle distance hn away from particle i.
Assuming a continuous distribution of displacement, a second order Taylor series can be expanded about particle i for a step n, which approximates u_{i+n}. The i is dropped from the notation as its value is not important so u_{i}, u_{i+n} → u, u_{n}, yielding:(28)
For 1D lattices the unit cell volume V is:(29)where b and w are the thicknesses of the rectangular prism unit cell in the 2 directions perpendicular to the length of the lattice. Substituting equation (28) into equation (27) and dividing by the unit cell volume V simplifies to:(30)
By comparing equation (30) to the 1D wave equation(31)where ρ is the material density, it can be seen that the equivalent Young’s modulus E is:(32)
Linearizing equation (3) it can be determined that:(33)where r_{n} is the distance from particle i to particle n and s is defined above.
4.1 Young’s modulus for two particle cells
For convenience, the modulus of elasticity is normalized to the unmagnetized case (ψ ^{(m)} = 0), E_{o}. If all magnets are repelling (s = 1) (Figure 1a) the normalized modulus of elasticity is:(34)where h_{o} = h/r_{0} and ζ is the Reimann zeta function defined as [19]:(35)
Since there exist methods to accurately calculate the Reimann zeta function and it is computationally reasonable to achieve a value for it to machine precision the radius of influence r_{∞} is not used in this case and no modeling error is introduced by omitting magnetic interactions in this 1D structure. Furthermore, the normalized magnetization is:(36)
The normalized value can be interpreted as the ratio of the magnetic force at r = r_{0} to the force required to fully compress the axial spring in the absence of contact, i.e. k_{a}r_{0}. Hence, this normalized quantity provides a reference of the strength of the magnetic interaction with respect to the spring restoring force. In general, will simply denote magnitude for a structure, since attraction and repulsion of magnets depends on the specific magnet pair in a lattice. Therefore when the terms “high” and “low” are used to refer to it is the magnitude that is being discussed [1].
If the polarity of the magnets alternates (s = −1) (Figure 1b) the normalized modulus of elasticity is:(37)
Now only h_{o} must still be defined. The distance between each particle h is a stable critical point of the energy functional equation (4) when reduced to its static terms. This leads to the following polynomial problem for nearest neighbors repelling and attracting respectively:(38) (39)
The variation of the equivalent Young’s modulus with is presented in Figure 3. The figure describes the tunability of the stiffness with changes in magnetization for simple structures defined by a 1 and 2 magnet unit cell. It is observed that the presence of repulsive magnetic interactions produces a hardening of the structure (solid black curve), and the presence of attractive nearest neighbor magnetic interactions has a softening effect (dashed red curve). There is a critical magnetization for the case with attracting nearest neighbors near where the structure approaches 0 stiffness as the structure becomes unstable. This is the magnetization magnitude at which the magnetic forces overpower the elastic forces from the axial springs. Beyond this point only reconfigured versions of the structure are stable.
Figure 3. Exact analytical solution of normalized modulus of elasticity for: the structure pictured in Figure 1a (black –) and the structure pictured in Figure 1b (red – –). 
4.2 Young’s modulus for four particle cells
The larger a unit cell becomes, the more possible configurations there are. The 1D results in this paper only consider unit cells of up to four magnets for simplicity. Furthermore, not all possible magnetic polarity patterns are considered. The results from four magnet unit cells are now discussed, with two magnet unit cells considered as a special case of four magnet unit cells. The group of configurations considered here are those pictured in Figures 1b–1g. The configurations in this group can reconfigure into any other configurations in the group, provided that both can exist for a given magnetization. Any structures that have contacting particles require a minimum magnetization to hold the particles in contact against the compressed spring’s force. In addition, when a structure is magnetized highly enough certain configurations will not be possible as magnetic forces will overcome elastic forces. Therefore, there are values of that do not permit all the configurations. For the following results a contact stiffness of α = 40 and a contact distance of r_{c} = 0.5r_{0} are used. The properties of the different configurations are shown in Figures 4 and 5. For this single system can take any of the six studied configurations, which have different stiffnesses and densities. The modulus of elasticity E and density ρ are normalized by the modulus of elasticity and density calculated when , which are E_{0} and ρ_{0}. Conveniently, the normalized density is the reciprocal of the normalized equilibrium unit cell width d, i.e. d/d_{0} = (ρ/ρ_{0})^{−1}, where d_{0} is the equilibrium distance calculated when . For the simplest cases, which are shown in Figure 3, d/d_{0} = 2h/2r_{0}) = h_{0}.
Figure 4. Normalized modulus of elasticity of structures with an even number of magnets per unit cell; configuration in Figure 1b (black circles), configuration in Figure 1c (red squares), configuration in Figure 1e (red – –), configuration in Figure 1f (black –), configuration in Figure 1g (blue – ·). 
Figure 5. Normalized density of structures with an even number of magnets per unit cell; configuration in Figure 1b (black circles), configuration in Figure 1c (red squares), configuration in Figure 1d (black rings), configuration in Figure 1e (red – –), configuration in Figure 1f (black –), configuration in Figure 1g (blue – ·). 
Figure 6. Comparison between normalized modulus of elasticity calculated by homogenization and numerical simulation for the configuration in Figure 1a (blue – · and black squares respectively), the configuration in Figure 1b (red – – and black circles respectively), and the configuration in Figure 1c (black – and black diamonds respectively). 
The Young’s modulus and its variation with is presented in Figure 4. Through reconfiguration alone there are as many values available as configurations. The structures corresponding to Figure 4 all soften in response to an increase of their magnetization. This is due to nearest neighbors that are in attraction, because the linearized model of attracting magnets exhibits negative axial stiffness. If both reconfiguration and changes in are considered there is no stiffness between the values of E/E_{0} ≈ 0.2 and 2 that cannot be achieved by at least one configuration, and there are some stiffnesses which can be achieved by as many as four configurations. This shows that stiffness adaptation is not limited to a few discrete values, allowing easier integration into possible design applications. The structure pictured in Figure 1d is shown separately in Figure 7 because its stiffness is much greater than the other structures’ as the contact stiffness used is α = 40. The structure is not a full 40 times stiffer though because the attractive magnetic interactions soften the structure, even under the minimum magnetization required for it to exist. Having only contact interactions, the structure in Figure 1d is the stiffest possible configuration for a structure with alternating magnetic poles and the parameters considered here, and thus gives an upper bound. The same can be said for the density, as contacting particles provides the closest packed structure.
Figure 7. Comparison between normalized modulus of elasticity calculated by homogenization and numerical simulation for the configuration in Figure 1d (red – – and black circles respectively). 
Figure 8. Normalized moduli and Poisson’s ratio for the hexagonal lattice with repelling magnets (a, c, and e) and with attracting nearest neighbors (b, d, and f). k _{ τ }/k _{ τo } = 10^{2} (red – –), k _{ τ }/k _{ τo } = 10^{0} (black triangles), k _{ τ }/k _{ τo } = 10^{−1} (blue – ·), and k _{ τ }/k _{ τo } = 10^{−3} (black –). 
Figure 9. Poisson’s ratio for the hexagonal lattice with repelling magnets (a) and with attracting nearest neighbors (b) for values of k _{ τ } that allow changes in to change the sign of ν. k _{ τ }/k _{ τo } = 10^{−0.1} (blue circles), k _{ τ }/k _{ τo } = 10^{−0.2} (red squares), k _{ τ }/k _{ τo } = 10^{−0.3} (black –), k _{ τ }/k _{ τo } = 10^{−0.4} (red – –), k _{ τ }/k _{ τo } = 10^{−0.5} (blue – ·), and k _{τ}/k _{ τo } = 10^{−0.6} (black diamonds). 
The densities of the 2 and 4 magnet configurations are presented in Figure 5. The density is primarily effected by the configuration, but there is some change with . The density can be approximately doubled through reconfiguration. The configurations pictured in Figures 1c and 1f exhibit similar densities because both have the same ratio of contacting to noncontacting magnets.
The configurations pictured in Figures 1c and 1f are worthy of note because they would be considered equivalent in terms of E and ρ if only nearest neighbor interactions were considered. In fact, for low enough they have almost equivalent properties. But with higher values of their dissimilar properties show that nonnearestneighbor interactions have a significant effect. Even so, the density does not change significantly with increasing , which means that a structure can have little change in its global size but a significant change in its stiffness by changing only its internal configuration. For example, if reconfiguring from the configuration in Figure 1c to the one in Figure 1f would produce only a 2.7% increase in density, but a 33% decrease in stiffness.
4.3 Numerical simulations
For the purpose of validating the homogenization method, as applied to the magnetoelastic lattices, numerical simulations are completed. The equations of motion are integrated using an explicit, variablestep RungeKutta method to obtain the time history of the positions of the particles in a finite lattice. This provides a framework for performing numerical experiments. A finite sample of each configuration to be validated, 40 particles long, is compressed quasistatically over a strain of less than 0.001, in which the structure can be considered linear. The compression is performed through an applied displacement of one end of the system while the other end is held fixed. There is linear damping included as described in reference [1] to damp out vibrations, so slow enough forcing will ensure that dynamic forces are negligible. To obtain a strain rate slow enough to consider quasistatic a convergence study is performed. The force versus displacement information is used to calculate a stiffness simply by comparing the change in reaction force to the structure’s axial strain. By normalizing to the case when the simulation results can be compared with the homogenization results directly. Their comparison is shown in Figures 6 and 7, with curves representing the approximate analytical results and the shapes representing the numerical results.
5 Results: 2D lattices
For simplicity the results presented for 2D magnetoelastic structures are restricted to systems described by a two magnet unit cell. A magnetization pattern where all magnets repel (Figure 2a) is used to investigate and compare to the pattern where nearest neighbors attract (Figure 2b). The magnetization pattern with attracting nearest neighbors is used to showcase the drastic changes in properties that can be achieved by reconfiguration. The reconfiguration considered is the transition to the reentrant lattice (Figure 2c) as this is possibly the easiest way to reconfigure the magnetoelastic hexagonal lattice studied based on previous research by the authors Schaeffer and Ruzzene [1, 20]. In contrast with the hexagonal patterns, the reentrant lattice is anisotropic. Figure 10a defines the 1 and 2 directions for the reentrant lattice, which are used to identify the orthotropic properties. Furthermore, the reentrant lattice cannot exist for all the same values of k_{τ}/k_{τo} and as the hexagonal lattice. Namely, there is an upper limit on k_{τ}/k_{τo} near 0.1. To explore the tunability of the equivalent mechanical properties of the aforesaid lattices the homogenization procedure as described in Section 3 is applied to structures with varying torsional stiffnesses k_{τ} and magnetization magnitudes . As in the previous section a contact stiffness of α = 40 and a contact distance of r_{c} = 0.5r_{0} are used. The properties are normalized with respect to a reference case defining the reference properties E_{o} and G_{o}, which is the unmagnetized () hexagonal lattice with k_{τo} = 1 N/rad and k_{ao} = 1 N/m.
Figure 10. Reentrant unit cell with coordinate system 1, 2 pictured (a), and θ defined as the angle between coordinate axes 1, 2 and x, y (b). 
The first results presented are for the two considered hexagonal lattice structures (Figures 2a and 2b). The homogenization produces properties that describe an isotropic material, i.e. E_{1} = E_{2} = E, ν_{12} = ν_{21} = ν, and G_{12} = E/(2(1 + ν)) = G, which is to be expected for regular hexagonal geometry [4]. Variation in k_{τ} is seen to cause substantial changes in material stiffness so the value of k_{τ} is incremented exponentially. Indeed, as k_{τ} → 0 the quantities E, G → 0 and the lattice approaches a mechanism. In this state the structure’s global deformation is dominated by the rotational deformation of the torsional springs. In such a case magnetization has a stiffening effect on the structure, which is evident in Figures 8a–8d for the cases where k_{τ}/k_{τo} = 10^{−3} (solid black line), though it is more prominent in the lattice with only repelling magnetic interactions. The stiffening of the structure can be explained predominantly through the repulsion of the nextnearest neighbor interactions, which are repulsive in both of the lattices discussed here. As the lattices are dominated by rotational deformation, the nearest neighbor magnetic interactions do not significantly effect the global properties. The exception comes when the the effective axial stiffness is nearly 0 close to the point of instability for attracting nearest neighbors, e.g. near . Though nextnearest neighbor interactions can explain the dominant trend for both magnetization patterns, comparison of Figures 8a and 8b reveals that more distant magnetic interactions have an nonnegligible effect on the structure. Thus, efforts to include the nonnearest neighbor interactions appear to be valuable. Significant tunability is seen for low k_{τ} when the magnetic interactions dominate the structure’s properties. For example, in the case where k_{τ}/k_{τo} = 10^{−3} in Figure 8a E increases by a factor of 48 over the range of shown.
Increasing k_{τ} causes the global deformation of the lattice to be dominated by the deformation of the axial springs in the lattice. In the limit as k_{τ} → ∞ the structure will become like the one discussed in reference [21]. The axial deformation dominated structures are stiffened by repulsive nearest neighbor interactions and softened by attractive nearest neighbor interactions, due to the equivalent linearized axial stiffness of magnetic interactions being positive for repulsion and negative for attraction. The greatest torsional stiffness shown is k_{τ}/k_{τo} = 10^{2} (dashed red line). Increasing k_{τ}/k_{τo} above this does not have a significant effect on the mechanical properties. In addition, decreasing the torsional stiffness below k_{τ}/k_{τo} = 10^{−3} does not make a significant change. Therefore, the properties of the hexagonal lattices discussed are bounded by the solid black and red dashed curves in Figure 8.
An interesting property of both magnetization patterns is that the Poisson’s ratio can be either positive or negative depending on the torsional stiffness k_{τ} and the magnetization (Figures 8e and 8f). In the limit as k_{τ} → 0, when the compliance of the system is dominated by the torsional springs, it is observed that ν → 1, as expected for the analogous structure of the regular hexagonal honeycomb with slender beams [4]. In the limit as k_{τ} → ∞ it is observed that ν → −1/3, which is consistent with equations given by Grima et al. [21] for regular hexagonal lattices with the angles between bars constrained. In this second case the structure’s compliance is dominated by the axial springs. Thus, it is clear that the ratio between rotational and axial stiffness affects the sign of ν. Since magnetization can have a similar effect to changing the rotational stiffness, as discussed above, for proper values of k_{τ} it is possible to switch between positive and negative ν simply by changing . The ranges of k_{τ} that allow for such tunability of ν for both magnetization patterns discussed are shown in Figure 9.
As stated previously the magnetoelastic hexagonal lattice with attracting nearest neighbors can be reconfigured into the reentrant lattice (Figure 2c) where strong contact forces hold the structure in equilibrium against magnetic attraction. The homogenized properties for varying values of k_{τ}/k_{τo} and are presented in Figure 11 and they are compared to the hexagonal lattices that reconfigure into them. In contrast to analyses of reentrant lattices such as those by Gibson and Ashby [4] or Gonella and Ruzzene [10], where the geometry of a cell is explicitly defined, the reentrant lattices discussed in this article have their geometry defined by the equilibrium configuration associated with a chosen value of . The equilibrium length of one axial spring relative to another may differ by more than 10% depending on the choice of k_{τ}/k_{τo} and , but the equilibrium lengths are approximately equal for low k_{τ}/k_{τo} and . Therefore, the tuning of properties with changing discussed hereafter are caused in part by the change in the reentrant lattice equilibrium geometry.
Figure 11. Normalized equivalent orthotropic properties for the reentrant lattice with α = 40 and k _{ τ }/k _{ τo } = 10^{−1} (blue · –), and k _{ τ }/k _{ τo } = 10^{−3} (black –). 
Reconfiguration from the hexagonal to the reentrant lattice has a drastic effect on the stiffness of the lattice. The drastic change is seen comparing in Figures 11, 8b, 8d, and 8f. Two choices of k_{τ}/k_{τo} are used. The value k_{τ}/k_{τo} = 10^{−1} is near the maximum value of k_{τ}/k_{τo} for which the reentrant lattice can exist, and the value k_{τ}/k_{τo} = 10^{−3} is representative of the asymptotic solution as k_{τ}/k_{τo} → 0. Values for E_{1} and E_{2} (Figures 11a and 11b) in the reentrant lattices are generally much higher than E in hexagonal lattices of equal k_{τ}/k_{τo}. The internal structure of the reentrant lattice makes strain along the one and two directions produce axially dominant local deformation, regardless of the torsional stiffness. Furthermore, the stiff contact interactions make E_{1} stiffer than values of E associated with axial deformation dominant hexagonal lattices with the same and magnetization pattern. In the two direction, the increased density from the reconfiguration makes values of E_{2} stiffer than associated values of E for axial deformation dominant hexagonal lattices. For the k_{τ}/k_{τo} = 10^{−3} case the Young’s moduli of the hexagonal configuration are dominated by magnetic interactions, so since the minimum magnetization for the existence of the reentrant configuration is low in this case the Young’s moduli can be increased by more than two orders of magnitude simply through structural reconfiguration to the reentrant lattice. Even for higher magnetizations there is over a 10fold increase in stiffness due to reconfiguration. For k_{τ}/k_{τo} = 10^{−1} the structure is still significantly stiffened by reconfiguration, but not as drastically as in the previous case. The k_{τ}/k_{τo} = 10^{−1} also has a more limited range of in which it can exist in the reentrant configuration.
As in the hexagonal lattice the shear modulus G_{12} in the reentrant lattice is dependent on k_{τ}. When k_{τ}/k_{τo} = 10^{−3} reconfiguration can account for over order of magnitude increase in G_{12}, but when k_{τ}/k_{τo} = 10^{−1} the shear moduli for the hex and reentrant lattices are more comparable. The Poisson’s ratios also change significantly with reconfiguration, with the reentrant lattices exhibiting lower Poisson’s ratios than corresponding hexagonal lattices in general. Although the reentrant lattice structure is generally known for being auxetic [4] this magnetoelastic example has positive ν_{12} and ν_{21} because of the stiff contact interactions that prevent the structure from folding in on itself. However, sufficiently soft contact interactions do allow negative ν_{12} and ν_{21}.
In the reentrant lattice the shear moduli are small compared to the axial stiffnesses, and this produces significant variation of stiffness with respect to direction. Figure 12 shows the variation of the equivalent properties in the x, y coordinate system, which is redefined here to be oriented by θ with respect to the nominal coordinate system 1, 2 as defined in Figure 10b. The properties of a reentrant lattice (solid black line) with k_{τ}/k_{τo} < 10^{−1.5} is compared to its hexagonal configuration (dashed blue line). Note that the radii of the polar plots in Figures 12a–12c are logarithmic. It is evident that the reentrant lattice is much softer for θ = 45 degrees (Figures 12a and 12b solid black line) than θ = 0 degrees, becoming over half an order of magnitude softer. As the Young’s moduli soften with changing θ the shear modulus stiffens by almost one order of magnitude (Figure 12c). Simultaneously, there are substantial changes in the Poisson’s ratios, though they remain positive. So, even when E_{1} and E_{2} are nearly equal the reentrant lattice is far from isotropic in contrast to the hexagonal lattice.
Figure 12. Variation of mechanical properties with coordinate axis orientation θ of the hexagonal (blue – –) and reentrant lattice (black –) with k_{τ}/k_{τo} = 10^{−1.5}, , α = 40. (a) log_{10}(E_{x}/E_{o}), (b) log_{10}(E_{y}/E_{o}), (c) log_{10}(G_{xy}/G_{o}), (d) ν_{xy}, (e) ν_{yx}. 
6 Conclusion
The equivalent continuum properties of magnetoelastic lattices are calculated through the extension of an existing homogenization method. Using 1D structures it is shown that many stiffnesses can be achieved by reconfiguring the structure between different stable equilibria and by changing the lattice magnetization. Using 2D structures the isotropic mechanical properties of hexagonal lattices are shown to vary with lattice physical properties and lattice magnetization. This includes the ability to change sign of the Poisson’s ratio of the structure simply by changing the lattice magnetization magnitude and the ability to stiffen or soften the lattice. The effect of reconfiguration from the hexagonal to the reentrant lattice on mechanical properties is presented, and the variation of the equivalent orthotropic properties of the reentrant lattice are shown to be tunable with changing magnetization magnitude.
Acknowledgments
This research was funded by a grant from the Air Force Office of Scientific Research Monitored by Dr. David Stargel, whose input and guidance are acknowledged and appreciated. The work of William Storrs, who performed the numerical simulations for validation, is also acknowledged and appreciated.
References
 M. Schaeffer, M. Ruzzene, Wave propagation in multistable magnetoelastic lattices, International Journal of Solids and Structures 56–57 (2014) 78–95. [Google Scholar]
 M. Schaeffer, M. Ruzzene, Wave propagation in reconfigurable magnetoelastic kagome lattice structures, Journal of Applied Physics 117 (2015) 194903. [CrossRef] [Google Scholar]
 P. Deymier, Acoustic Metamaterials and Phononic Crystals, Springer Series in SolidState Sciences, Springer, New York, 2013. [CrossRef] [Google Scholar]
 L.J. Gibson, M.F. Ashby, Cellular solids: structure and properties, Cambridge university press, Cambridge, UK, 1999. [Google Scholar]
 M.I. Hussein, M.J. Leamy, M. Ruzzene, Dynamics of phononic materials and structures: Historical origins, recent progress, and future outlook, Applied Mechanics Reviews 66 (2014) 040802. [CrossRef] [Google Scholar]
 M.H. Lu, L. Feng, Y.F. Chen, Phononic crystals and acoustic metamaterials, Materials Today 12 (2009) 34–42. [CrossRef] [Google Scholar]
 T.A. Schaedler, A.J. Jacobsen, W.B. Carter, Toward lighter, stiffer materials, Science 341 (2013) 1181–1182. [CrossRef] [Google Scholar]
 H.N. Wadley, Multifunctional periodic cellular metals, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 364 (2006) 31–68. [CrossRef] [Google Scholar]
 M. Born, K. Huang, Dynamical theory of crystal lattices, Clarendon Press, Oxford, 1954. [Google Scholar]
 S. Gonella, M. Ruzzene, Homogenization and equivalent inplane properties of twodimensional periodic lattices, International Journal of Solids and Structures 45 (2008) 2897–2915. [CrossRef] [Google Scholar]
 A.A. Maradudin, E.W. Montroll, G.H. Weiss, I. Ipatova, Theory of lattice dynamics in the harmonic approximation, Vol. 12, Academic press, New York, 1963. [Google Scholar]
 A. Suiker, A. Metrikine, R. De Borst, Comparison of wave propagation characteristics of the Cosserat continuum model and corresponding discrete lattice models, International Journal of Solids and Structures 38 (2001) 1563–1583. [CrossRef] [Google Scholar]
 J.N. Grima, R. CaruanaGauci, M.R. Dudek, K.W. Wojciechowski, R. Gatt, Smart metamaterials with tunable auxetic and other properties, Smart Materials and Structures 22 (2013) 084016. [CrossRef] [Google Scholar]
 E. du Trémolet de Lacheisserie, D. Gignoux, M. Schlenker (Eds.), Magnetism: fundamentals, Springer, New York, 2005. [Google Scholar]
 K.W. Yung, P.B. Landecker, D.D. Villani, An analytic solution for the force between two magnetic dipoles, Magnetic and Electrical Separation 9 (1998) 39–52. [CrossRef] [Google Scholar]
 H. Goldstein, Classical Mechanics, AddisonWesley Publishing Company, Inc, Cambridge, MA, 1953. [Google Scholar]
 S.S. Rao, Engineering Optimization: Theory and Practice. 4th ed., John Wiley & Sons, New York, 2009. [Google Scholar]
 R.D. Cook, D.S. Malkus, M.E. Plesha, Concepts and applications of finite element analysis. 3rd ed., John Wiley & Sons, New York, 1989. [Google Scholar]
 M. Abramowitz, I.A. Stegun, et al., Handbook of mathematical functions, Vol. 1, Dover, New York, 1972. [Google Scholar]
 M. Schaeffer, M. Ruzzene, Dynamic reconfiguration of magnetoelastic lattices, Comptes Rendus Mécanique 343 (12) (2015) 670–679. [CrossRef] [Google Scholar]
 J.N. Grima, R. CaruanaGauci, K.W. Wojciechowski, K.E. Evans, Smart hexagonal truss systems exhibiting negative compressibility through constrained angle stretching, Smart Materials and Structures 22 (2013) 084015. [CrossRef] [Google Scholar]
Cite this article as: Schaeffer M & Ruzzene M: Homogenization of 1D and 2D magnetoelastic lattices. EPJ Appl. Metamat. 2015, 2, 13.
All Figures
Figure 1. Periodic equilibrium examples for 1D structures with one (a), two (b–d), and four (e–g) particle unit cells. Red and blue circles denote magnets with dipole moments perpendicular to the lattice plane pointing towards and away from the reader respectively. Black lines denote elastic connections between magnets. Reprinted from M. Schaeffer and M. Ruzzene, Wave propagation in multistable magnetoelastic lattices, Int. J. Solids Struct. 56–57 (2015) 78–95, copyright (2015), with permission from Elsevier. 

In the text 
Figure 2. Periodic equilibrium examples for 2D structures with two particle unit cells. Red and blue circles denote magnets with dipole moments perpendicular to the lattice plane pointing towards and away from the reader respectively. Black lines denote elastic connections between magnets. Particles in the unit cell have been highlighted by placing yellow dots in their centers. Reprinted from M. Schaeffer and M. Ruzzene, Wave propagation in multistable magnetoelastic lattices, Int. J. Solids Struct. 56–57 (2015) 78–95, copyright (2015), with permission from Elsevier. 

In the text 
Figure 3. Exact analytical solution of normalized modulus of elasticity for: the structure pictured in Figure 1a (black –) and the structure pictured in Figure 1b (red – –). 

In the text 
Figure 4. Normalized modulus of elasticity of structures with an even number of magnets per unit cell; configuration in Figure 1b (black circles), configuration in Figure 1c (red squares), configuration in Figure 1e (red – –), configuration in Figure 1f (black –), configuration in Figure 1g (blue – ·). 

In the text 
Figure 5. Normalized density of structures with an even number of magnets per unit cell; configuration in Figure 1b (black circles), configuration in Figure 1c (red squares), configuration in Figure 1d (black rings), configuration in Figure 1e (red – –), configuration in Figure 1f (black –), configuration in Figure 1g (blue – ·). 

In the text 
Figure 6. Comparison between normalized modulus of elasticity calculated by homogenization and numerical simulation for the configuration in Figure 1a (blue – · and black squares respectively), the configuration in Figure 1b (red – – and black circles respectively), and the configuration in Figure 1c (black – and black diamonds respectively). 

In the text 
Figure 7. Comparison between normalized modulus of elasticity calculated by homogenization and numerical simulation for the configuration in Figure 1d (red – – and black circles respectively). 

In the text 
Figure 8. Normalized moduli and Poisson’s ratio for the hexagonal lattice with repelling magnets (a, c, and e) and with attracting nearest neighbors (b, d, and f). k _{ τ }/k _{ τo } = 10^{2} (red – –), k _{ τ }/k _{ τo } = 10^{0} (black triangles), k _{ τ }/k _{ τo } = 10^{−1} (blue – ·), and k _{ τ }/k _{ τo } = 10^{−3} (black –). 

In the text 
Figure 9. Poisson’s ratio for the hexagonal lattice with repelling magnets (a) and with attracting nearest neighbors (b) for values of k _{ τ } that allow changes in to change the sign of ν. k _{ τ }/k _{ τo } = 10^{−0.1} (blue circles), k _{ τ }/k _{ τo } = 10^{−0.2} (red squares), k _{ τ }/k _{ τo } = 10^{−0.3} (black –), k _{ τ }/k _{ τo } = 10^{−0.4} (red – –), k _{ τ }/k _{ τo } = 10^{−0.5} (blue – ·), and k _{τ}/k _{ τo } = 10^{−0.6} (black diamonds). 

In the text 
Figure 10. Reentrant unit cell with coordinate system 1, 2 pictured (a), and θ defined as the angle between coordinate axes 1, 2 and x, y (b). 

In the text 
Figure 11. Normalized equivalent orthotropic properties for the reentrant lattice with α = 40 and k _{ τ }/k _{ τo } = 10^{−1} (blue · –), and k _{ τ }/k _{ τo } = 10^{−3} (black –). 

In the text 
Figure 12. Variation of mechanical properties with coordinate axis orientation θ of the hexagonal (blue – –) and reentrant lattice (black –) with k_{τ}/k_{τo} = 10^{−1.5}, , α = 40. (a) log_{10}(E_{x}/E_{o}), (b) log_{10}(E_{y}/E_{o}), (c) log_{10}(G_{xy}/G_{o}), (d) ν_{xy}, (e) ν_{yx}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.