A Generalized Strain Energy-Based Homogenization Method for 2-D and 3-D Cellular Materials with and without Periodicity Constraints

: A generalized strain energy-based homogenization method for 2-D and 3-D cellular materials with and without periodicity constraints is proposed using Hill’s Lemma and the matrix method for spatial frames. In this new approach, the equilibrium equations are enforced at all boundary and interior nodes and each interior node is allowed to translate and rotate freely, which differ from existing methods where the equilibrium conditions are imposed only at the boundary nodes. The newly formulated homogenization method can be applied to cellular materials with or without symmetry. To illustrate the new method, four examples are studied: two for a 2-D cellular material and two for a 3-D pentamode metamaterial, with and without periodic constraints in each group. For the 2-D cellular material, an asymmetric microstructure with or without periodicity constraints is analyzed, and closed-form expressions of the effective stiffness components are obtained in both cases. For the 3-D pentamode metamaterial, a primitive diamond-shaped unit cell with or without periodicity constraints is considered. In each of these 3-D cases, two different representative cells in two orientations are examined. The homogenization analysis reveals that the pentamode metamaterial exhibits the cubic symmetry based on one representative cell, with the effective Poisson’s ratio 𝑣̅ being nearly 0.5. Moreover, it is revealed that the pentamode metamaterial with the cubic symmetry can be tailored to be a rubber-like material (with 𝑣̅ ≅ 0.5 ) or an auxetic material (with 𝑣̅ < 0).

Various homogenization methods have been developed using classical elasticity. For example, Warren and Kraynik [12] proposed an analytical method to homogenize lowdensity open-cell foams based on solving force and moment equilibrium equations at joints. Tollenaere and Caillerie [13] used an asymptotic expansion method in homogenizing 2-D lattice truss structures. Li et al. [14] developed a micromechanics model and obtained closed-form formulas for predicting effective elastic properties of 3-D open-cell foams based on Castigliano's second theorem. Demiray et al. [15] homogenized 2-D and 3-D hyperelastic foams undergoing large deformations by employing a strain energybased scheme and periodic conditions. Martinsson and Babuška [16] provided a homogenization method for materials with periodic truss and frame microstructures utilizing an asymptotic analysis built upon Fourier transforms. Freund et al. [17] presented a two-scale computational homogenization technique for 2-D cellular structures based on the virtual work principle. Norris [18] proposed a homogenization method for periodic lattice structures and derived analytical formulas for effective properties of elastic networks, including pentamode metamaterials. Ongaro et al. [19] employed a strain energy-based approach to homogenize 2-D cellular structures of honeycomb cells filled with an elastic material. Ai and Gao [20] provided an analytical model for predicting effective elastic properties of 2-D periodic star-shaped re-entrant lattice structures with the orthotropic symmetry and negative Poisson's ratios by using Castigliano's second theorem and the Timoshenko beam theory. Materials with negative Poisson ratios are also known as auxetic materials, which have emerged as a class of metamaterials that can find important engineering applications (e.g., [6,[20][21][22][23][24][25][26]). Recently, Czarnecki and Łukasiak [27] applied the asymptotic homogenization method for periodic media to estimate effective moduli of 2-D auxetic cellular materials, which were also compared with those obtained through optimization.
In the current study, a generalized strain energy-based homogenization method for 2-D and 3-D cellular materials with and without periodicity constraints is developed, which is built on classical elasticity and has no restriction on shape, symmetry or number of struts in a unit cell. In this new approach, the nodal equilibrium equations are enforced at all boundary and interior nodes, unlike in existing classical elasticity-based methods where the nodal equilibrium is imposed only at the boundary nodes and, as a result, the equilibrium equations are often not satisfied at the interior nodes by the approximate solutions obtained (e.g., [35]).
The rest of the paper is organized as follows. In Section 2, the generalized strain energy-based homogenization method is formulated for cellular materials with and without periodicity constraints. In Section 3, the newly proposed method is applied to homogenize 2-D and 3-D cellular materials in four example problems, which leads to closed-form formulas for effective elastic stiffness and compliance components. In Section 4, the paper concludes with a summary.

Matrix Method for Spatial Frames
According to the matrix method for spatial frames (e.g., [36,37]), each frame member is regarded as weightless and loaded only at its two end points (nodes). It is rigidly connected to other members at the two end nodes, each having six degrees of freedom-three translational and three rotational displacements-if unconstrained.
For a 3-D frame member with two end nodes I and J, denoted as "I, J " and shown in Figure 1, which is made from an isotropic linear elastic material and has a uniform circular cross-section, the stiffness matrix is given by (1) where the 6 × 6 sub-stiffness matrix IJ k represents the force at node I due to a unit displacement at node J. In the local coordinate system {x, y, z} with the base vectors { 1 , 2 , 3 }, as shown in Figure 1, the sub-stiffness matrices can be written as (e.g., [37] (2) where the superscript "T" denotes the transpose of the matrix, L is the length of the frame member, and , , and a s b t k k k k are, respectively, the axial stiffness along the x-axis, transverse shear stiffness about the y-or z-axis, bending stiffness about the y-or z-axis, and torsional stiffness about the x-axis. When the Timoshenko beam theory is used [37], 43 12 , , , , 12 12 which can be reduced to those based on the Bernoulli-Euler beam theory given by   (7) which is related to Δ (see Equation (6) ij  , (9) in which Q is the orthogonal coordinate transformation matrix from the local coordinate system {x, y, z} to the global coordinate system {X, Y, Z} (i.e., i = Q ), and the summation on i and j (with i, j ∈{1, 2, 3}) is implied. Note that is also proper orthogonal with 1 T − = and det 1 = , where the superscript "−" denotes the inverse of the matrix.
The nodal force vector of the member "I, J " that satisfies the equilibrium is given by (e.g., [37]) == F M F M Γ kΔ (10) or in terms of the components in the local coordinate system,  (12) where , T = Kk (13) is the stiffness matrix of the member "I, J " in the global coordinate system {X, Y, Z}.
For a cellular material containing W frame members, the total strain energy is given by 1 1 , 2 W T n n n n= =  D K D (14) where n ({1, 2, …, W}) represents the nth frame member, and . T n n n n = Kk (15)

Hill's Lemma
Hill's lemma [38] enables the prediction of effective properties of a heterogeneous material through constructing a homogeneous comparison solid based on the strain energy equivalence (e.g., [39][40][41]).
For the Cauchy continuum, Hill's lemma reads (e.g., [38,39,42,43]) where Ω is the region occupied by a heterogeneous material, Ω is the closed surface of Ω, V is the volume of Ω, dS is an area element on Ω, σ and ε are, respectively, the Cauchy stress and infinitesimal strain tensors,    denotes the volume-averaged quantity, u is the displacement vector, E and  are, respectively, the volume-averaged (effective) stress and strain tensors, x is the position vector, and n is the unit outward normal to Ω.  (18) in which the overhead bar represents the prescribed quantity. From Equations (17) and (18), the displacement u at each boundary node can be obtained. However, the rotation θ can vary independently at the boundary nodes while satisfying the moment equilibrium equations there (e.g., [44]). Accordingly, the displacement vector for the nth frame member "I, J " (with n{1, 2, …, W}) with the end nodes I and J both lying on the boundary of the cellular material can be rewritten as If node I lies on the boundary and node J is located in the interior, then the displacement vector for the nth frame member "I, J " has the form: (20) Finally, if both the end nodes I and J lie in the interior of the cellular material, then = Du θ u θ (21)

Generalized Homogenization Method
Consider a general (asymmetric) 3-D cellular material with no periodicity, as shown in Figure 2. The cellular material is composed of P boundary frame members (struts), each of which has at least one node on the boundary, and Q interior struts, each of which has its two end nodes located inside the boundary. The cellular material also features N nodes that lie on the boundary and M nodes that are located inside the boundary. The cellular material can be completely characterized by N, M, P, Q and the coordinates of each node. To homogenize this cellular material with the equilibrium satisfied at all nodes (on the boundary and in the interior), the displacement boundary conditions in Equation (17) are prescribed for the boundary nodes. The remaining nodes inside the boundary are allowed to translate and rotate without any constraint. Thus, for the N boundary nodes, the following holds: where n ({1, 2, …, N}) denotes the nth boundary node, and N is the total number of struts emanating from the nth boundary node. This leaves the rotation vector at each boundary node, n θ (n{1, 2, …, N}), unspecified, while satisfying the nodal moment equilibrium. Then, the displacement vector for the pth boundary strut "I, J ", with node I lying on the boundary and node J located in the interior, can be obtained from Equations (20) and (22) (28) Then, the volume-averaged strain energy, called the strain energy density function, in the cellular material with P boundary struts and Q interior struts can be obtained as , PQ u VV + == (29) where V is the volume of the region enclosed by the bounding surface of the cellular material.
The strain energy density function u in Equation (29) contains a number of unknown displacement and rotation components. For a general 3-D case, these unknowns include 3N rotation components at the N boundary nodes ( n θ ; n{1, 2, …, N}), 3M displacement components at the M interior nodes ( m u ; m{1, 2, …, M}), and 3M rotation components at the M interior nodes ( m θ ; m{1, 2, …, M}). As a result, there are totally 6M + 3N unknown displacement and rotation components in u. Hence, an equal number of equations are required to solve for these unknowns so that the strain energy density function will not contain any undetermined displacement or rotation component.
Since each boundary node is allowed to rotate freely, enforcing the moment equilibrium at all the N boundary nodes provides 3N equations. These moment balance equations at the boundary nodes read Furthermore, since each interior node is allowed to translate and rotate without constraints, enforcing the force and moment equilibrium at all the M interior nodes gives 6M equations, which are where M is the total number of struts emanating from the mth interior node.
Solving Equation (30)  Substituting these determined kinematic variables into Equation (29) will give the strain energy density function u in its final form without containing any unknown. The effective stiffness tensor can then be obtained from u as 2 2 ,, uu  ==  C   (32) where E and  are the constant strain and stress tensors defined in Equation (18), and C is the effective stiffness (elasticity) tensor.

Extension to Periodic Materials
For periodic cellular materials, the homogenization can be performed using an extended version of the approach based on Hill's lemma discussed in Section 2.3.
For a periodic material, Hill's lemma given in Equation (16) can be extended to (e.g., [40]) where * u is the periodic part of the displacement field, which is the same for each periodic pair of nodes on the unit cell boundary. Note that for each periodic pair of two boundary nodes the traction vector = t σn is anti-periodic with +− += t t 0 (e.g., [40,45,46]). Accordingly, the surface integral of the product ( ) * − u σ n  vanishes (e.g., [40,43]), and hence Equation (33) is equivalent to Equation (16). Applying the Hill-Mandel condition to Equation (33) then gives or on x u σn n  (34) as the periodic BCs.
Consider a general 3-D periodic cellular material, with a unit cell shown in Figure 3a or Figure 3b. Each unit cell contains P boundary struts, each of which has one node on the boundary, and Q interior struts, each of which has its two nodes located inside the boundary. Each unit cell also includes B periodic pairs of nodes that lie on the boundary and M nodes that are located inside the boundary.  A periodic pair can contain two or more boundary nodes, which depends on the unit cell structure. For example, the unit cell shown in Figure 3a possesses three periodic pairs (B = 3), each having two distinct nodes, while the unit cell in Figure 3b contains one periodic pair (B = 1), with all the boundary nodes included in this pair. In general, for a periodic cellular material, a periodic pair can be defined as a collection of boundary nodes that satisfy the following conditions (e.g., [40,45]):  (36) where H represents the total number of nodes contained in the periodic pair. The kinematic boundary conditions shown in Equation (34) can be extended to represent the periodic boundary conditions for a cellular material in the following form: Then, it follows from Equations (12) and (38) that the force vector for the pth strut "I, J " has the form: and from Equations (14) and (38) that the total strain energy stored in the P boundary struts can be obtained as The displacement vector and force vector for the qth interior strut can be computed from Equations (26) and (27), respectively, and the total strain energy stored in the Q interior struts can be determined using Equation (28).
Then, the strain energy density function u in the unit cell with P boundary struts and Q interior struts can be obtained from Equations (28), (29) and (40). There are The 6B unknown periodic parts of the displacement and rotation components at the boundary nodes can be obtained by enforcing the anti-periodicity conditions of forces and moments for each periodic pair given in Equation (36). In addition, the 6M displacement and rotation components at the M interior nodes can be identified by imposing the force and moment equilibrium at those nodes according to Equation (31). Using these determined displacement and rotation components in Equation (29) will give the final expression of u, which can be used in Equation (32) to find the effective stiffness tensor C .

Case Studies
In this section, 2-D and 3-D cellular materials with and without periodicity constraints are homogenized by applying the generalized strain energy-based method proposed in Section 2. The mathematical formulation here is facilitated by using symbolic operations in MATLAB. For simplicity, in all cases considered here, struts are taken to be Bernoulli-Euler beams with uniform circular cross sections whose stiffness constants are given in Equation (4).

2-D Homogenization without Periodicity Constraints
Consider a 2-D cellular material shown in Figure 4, which is asymmetric. The unit cell is composed of three struts with equal length L (i.e., (1) (1) , , , a a a a s s s where ka, ks and kb are listed in Equation (4). The area of the unit cell is   The coordinate transformation tensor for each of the three struts with respect to the global coordinate system {X, Y, Z} is given by In addition, the position vector of each boundary node with respect to the origin o (the center of the unit cell; an interior node) reads Substituting Equations (18) and (44) into Equation (22) gives the boundary conditions at the three boundary nodes as 3 12 3 12 ,, Substituting Equations (47)- (49) into Equation (46) or (30), which gives the moment equilibrium equations at the three boundary nodes in Figure 4b, leads to the rotation components at the three boundary nodes as Using Equation (47)-(49) and (50) in Equation (31), the force and moment equilibrium equations at the interior node o, results in Upon solving the system of three equations in Equation (51), the displacement and rotation components at the interior node can be obtained as Substituting Equations (1), (2), (9), (15), (25), (41), (43), (45), (50) and (53) into Equation (29) gives the stain energy density function as Clearly, Equation (54) does not contain any unknown kinematic variable. (32) and (54) that the effective stiffness matrix C for the 2-D cellular material can be obtained from the coefficient matrix of the following constitutive equations:

It then follows from Equations
Note that only one specific configuration shown in Figure 4 is considered in the 2-D case studies here. However, the inclination angle of struts in the 2-D asymmetric cellular material displayed in Figure 4a could be set as an adjustable variable, which would lead to the orientation dependence of the effective elastic properties upon using the current homogenization method. Such orientation dependence exhibited by 2-D cellular materials has been extensively studied (e.g., [20,25,47]).

2-D homogenization with Periodicity Constraints
In this sub-section, periodicity constraints are imposed on the unit cell shown in Figure 4b, which is further illustrated in Figure 5. It can be seen from Figure 5 that the three nodes on the boundary of the unit cell form one periodic pair (with B = 1 and H = 3). Based on the conditions listed in Equations (35) and (36), the periodicity constraints for the unit cell shown in Figure 5 take the following form: Substituting Equations (45) and (57) (the first set) into Equation (37) yields the periodic boundary conditions at the three boundary nodes as ** 12 ** 12 ** 12 Using Equations (1), (2), (7), (9), (13), (43) and (58) in Equation (12) gives the nodal force vector for each of the three struts as Substituting Equations (59)-(61) into Equation (57) (the second set) results in Upon solving Equation (62), the periodic parts of the displacement and rotation components at the three boundary nodes are obtained as A comparison of Equations (67) and (68) with Equations (55) and (56) shows that the effective stiffness components obtained from the unit cell with the periodicity constraints (see Equation (68)) do not depend on the strut bending stiffness kb, while those determined from the same unit cell without the periodicity constraints (see Equation (56)) are dependent on kb as well as ka and ks.

Pentamode Metamaterial
Pentamode metamaterials were first proposed by Milton and Cherkaev [48] using a diamond-shaped primitive unit cell shown in Figure 6b. This unit cell can be characterized by three independent primitive lattice vectors (red arrows in Figure 6b). Struts extending from each of the four vertices meet at the center of the cube C, which is located inside the primitive cell. This basic structure is repeated periodically in the primitive lattice vector directions to generate the periodic lattice material shown in Figure 6a. The unit cell shown in Figure 6b represents the ideal (perfect) pentamode metamaterial having a zero-shear resistance, where each strut (composed of two truncated cones) overlaps with the other struts at a single point. In the current study, each strut is taken to be a circular cylinder with a uniform cross-section. This allows for a simplified model, while retaining the essential features of the pentamode metamaterial.
The volume and volume fraction of the primitive diamond-shaped unit cell shown in Figure 6b are given by, for struts with a uniform circular cross-section, where Vd and fd are, respectively, the volume and volume fraction of the diamond unit cell, and L and d are, respectively, the length and diameter of each strut. Note that ) (e.g., [20]).
Pentamode metamaterials have been shown to display the isotropic, transversely isotropic or orthotropic symmetry, depending on whether the eigenvalues are distinct or repeated [18]. In addition, pentamode metamaterials can be characterized by three different unit cells: primitive diamond, cubic, and parallelepiped [11]. The parallelepiped unit cell can be obtained from the cubic unit cell through a simple rotation of the coordinate system.
Homogenization in two different coordinate systems shown in Figure 7 is considered in this subsection to understand how the effective engineering constants transform under the coordinate system change.

Homogenization without Periodicity Constraints
The generalized strain energy-based homogenization method proposed in Section 2 is applied here to predict the effective elastic properties of the pentamode metamaterial using the diamond-shaped unit cell viewed in two different orientations (coordinate systems) shown in Figure 8a,b. Each of the two representative cells associated with the two orientations in Figure 7 contains four struts of equal lengths L and identical stiffness constants ka, ks and kb, which are the same as those included in the diamond-shaped unit cell shown in Figure 7. In each case, the volume fraction is that of the pentamode metamaterial represented by the primitive diamond-shaped unit cell.  The difference between the two representative cells is that they are 45 apart on the X-Y plane, while the Z-axis is in the same vertical direction in both cases (see Figure 7). That is, the base vectors of the two coordinate systems shown in Figure 8a,b are related through To construct the pentamode metamaterial from the representative cell in orientation # 1 (the blue cube) in Figures 7 and 8a or 9a, the 4-strut structure in the blue cube is first repeated once along each of the three primitive lattice vector directions (the red arrows in Figure 9a), which results in a new unit cell containing 16 struts enclosed in the pink cube, as displayed in Figure 9b. Then, repeating the 16-strut unit cell in the pink cube along the X1, Y1 and Z1 directions will generate the 3-D periodic pentamode lattice structure, as shown in Figure 9c. For the representative cell shown in Figure 8a, which will be referred to as orientation # 1, the coordinate transformation matrices for the four struts with respect to the global coordinate system {X1, Y1, Z1} can be obtained as 1  1  1  2  2  2  2  21  22  23  21  22  23   1  1  1  2  2  2  31  32  33  31  32  33   3  3  3  3  4  4  4  4  21  22  23  21  22  23   3  3  1  4  31  32  33  31  32   1  1  1  1  1  1  3  3  3  3  3  3   ,  ,   1  1  1  1  1  1  3  3  3  3  3  3 , This is due to the fact that the cross-section of each strut is circular, so that the local coordinate axes y and z can be arbitrarily oriented on a cross-section. Based on these constraints, the following transformation matrices with respect to the coordinate system {X1, Y1, Z1} have been chosen: (1) (3) , .
In addition, the position vector for each of the five nodes in the unit cell shown in  , , which does not involve any unknown kinematic variable.
Finally, from Equations (32) and (82) the effective stiffness matrix (1) C for the pentamode metamaterial based on the representative cell in orientation # 1 shown in Figure  8a can be obtained as the coefficient matrix of the following constitutive equations:    , which indicates that the pentamode metamaterial possesses the cubic symmetry in the current case with the representative cell in orientation # 1. These stiffness constants are the same as those obtained in [18] for the pentamode material with the diamond unit cell using an equilibrium equation-based approach, thereby providing a validation of the current model. The effective compliance matrix (1) based on the representative cell in orientation # 1 can be readily obtained from Equations (83) and (84) as For a cubic material, the directional dependence of the engineering constants can be described by (e.g., [49,50] (n1, n2, n3) and m = (m1, m2, m3).  ,, ,, which establish the functional relations between the effective engineering constants and the angle θ for the pentamode metamaterial using the stiffness matrix based on the representative cell in orientation # 1 shown in Figure 8a.

C derived in Equations (83) and
(84) based on the representative cell in orientation # 1. The effective compliance matrix (2) (2) 1 [] − = SC based on the representative cell in orientation # 2 can then be readily obtained from Equations (102) and (103) as (2) (2) 12 (2) (2) 12 (2) (2)  Figure 8b can be reproduced from the effective stiffness and compliance matrices derived in Equations (83) and (84) and Equations (85) and (86) based on the representative cell in orientation # 1 displayed in Figure 8a by using a coordinate transformation, as shown next. When the coordinate system # 1 is rotated by an angle of θ about the common Z-axis to the coordinate system # 2, the effective stiffness and compliance matrices satisfy the following transformation relations [51] (2) (1) (2) (1) 1 ,, These verify and support the newly proposed homogenization approach.

Homogenization with Periodicity Constraints
Periodicity constraints are considered in homogenizing the pentamode metamaterial shown in Figure 10 using the method proposed in Section 2.4. As discussed earlier, the unit cell is repetitive along three directions which are not mutually perpendicular to each other. From Figure 10, it is seen that the pentamode material can be generated by repeating the diamond-shaped unit cell along the directions BA , BD and BE . This indicates that and for the periodic pair of the four boundary nodes based on the representative cell in orientation # 2,   .
Using Equations (1), (2), (7), (9), (12), (13), (73), (110) and (111) in the anti-periodicity conditions listed in Equation (108) yields the periodic parts of the displacement and rotation at the four boundary nodes as Next, applying the equilibrium conditions at the central node o (an interior node) listed in Equation (31) which no longer depends on the strut bending stiffness kb, unlike that obtained in Equation (82) without the periodicity constraints. This is consistent with what is observed from comparing the 2-D cases with and without the periodic constraints in Section 3.1. Then, the effective stiffness matrix (1) P C for the pentamode metamaterial based on the representative cell in orientation # 1 shown in Figure 8a with the periodic boundary conditions listed in Equations (110)-(112) can be obtained from Equations (32) and (114)   Clearly, Equations (115) and (116) show that the pentamode metamaterial with the periodicity constraints based on the representative cell in orientation # 1 exhibits the cubic symmetry, which is the same as what is observed from Equations (83) and (84) for the same representative cell and orientation but without the periodicity constraints. In addition, it is seen from Equations (115) and (116) that the effective stiffness components do not depend on the strut bending stiffness kb, which differs from those given in Equations (83) and (84) for the pentamode metamaterial without imposing the periodicity constraints.

Homogenization Based on the Second Representative Cell with Periodicity Constraints
For the representative cell in orientation # 2 shown in Figure 8b, the displacement and rotation vectors at the four boundary nodes can be obtained from Equations (18) Clearly, this strain energy density function does not depend on the strut bending stiffness b k , which is different from that obtained in Equation (101) based on the representative cell shown in Figure 8b without using the periodicity constraints. This is similar to what is observed from comparing Equations (114) and (82), which are based on the representative cell shown in Figure 8a with and without the periodicity constraints.
Then, it follows from Equations (32) and (123)  (2)  (2) (2) 1 2 3  The compliance matrix (2) based on the representative cell in orientation # 2 with the periodicity constraints can then be readily obtained as (2) (2) 12 (2) (2) 12 (2) (2) It can be readily shown that the effective stiffness and compliance matrices obtained in Equations (124) and (126) based on the representative cell in orientation # 2 can be reproduced from those in Equations (115) and (117) based on the representative cell in orientation # 1 through a coordinate transformation given in Equations (106) and (107) Figure  8a), the effective Young's modulus is the smallest, the effective shear modulus is the largest, and the effective Poisson's ratio is nearly 0.5. However, at 45  =− (i.e., the representative cell in orientation # 2 shown in Figure 8b), the effective Young's modulus is the largest, the effective shear modulus is the smallest, and the effective Poisson's ratio reaches zero. Figure 11 also reveals that both the effective Young's modulus and shear modulus of the pentamode metamaterial with the periodicity constraints are higher than those without the constraints for any orientation θ ∈[−90, 0], while the effective Poisson's ratio is almost the same in the two cases for all orientations. The former reaffirms what has been qualitatively observed earlier from analyzing Equation (132). In addition, the ratios given in Equation (132) are plotted in Figure 12. It is seen from Figure 12 that the smaller r (= d/L) is, the closer the results are to the limiting values listed in Equation (133), as expected. Moreover, it is observed from Figure 12 that the ratio of the effective shear moduli

Conclusions
A generalized strain energy-based homogenization method is provided for 2-D and 3-D cellular materials with and without periodicity constraints by using Hill's lemma and the matrix method for spatial frames. In the newly proposed homogenization method, the equilibrium equations are imposed at all boundary and interior nodes, and the interior nodes are allowed to translate and rotate without constraints, unlike in the existing strain energy-based methods that do not enforce the equilibrium conditions at the interior nodes. The new method can be applied to all types of 2-D and 3-D cellular structures with no geometrical constraints, which differs from the existing methods that can only be used to accurately predict effective elastic properties of cellular materials with symmetric or simple microstructures.
An asymmetric 2-D cellular material and a 3-D pentamode metamaterial, with and without periodicity constraints in each group, are homogenized by directly using the new method. In all four cases, closed-form expressions are obtained for the components of the effective stiffness matrix. For the 3-D pentamode metamaterial, the effective stiffness and compliance components are derived using two representative cells in two different orientations for both cases with or without periodicity constraints. The results based on the representative cell in one orientation show that the pentamode metamaterial displays the cubic symmetry and can be tailor-designed to be a rubber-like incompressible material or an auxetic material with a negative Poisson's ratio.
Finally, it should be mentioned that the effective stiffness matrix obtained without the periodic constraints has been found to be different from that acquired with the periodic constraints in each of the 2-D and 3-D homogenization examples considered. These analytical results attained using the newly proposed homogenization method need to be compared with experimental data, when they become available, to further verify the new method and to determine which one, with or without the periodic constraints, is more accurate.