A Refined Finite Element Formulation for the Microstructure-Dependent Analysis of Two-Dimensional (2D) Lattice Materials

A finite element approximation is proposed for the dynamic analysis of two-dimensional (2D) lattice materials. The unit cell is modeled by means of a defined number of shear deformable micro-beams. The main innovative feature concerns the presence of a microstructure-dependent scale length, which allows the consideration of the so called size-effect that can be highly relevant, due to the characteristics of the lattice at the local scale. Some numerical results show the influence of the microstructure parameter on the dynamic behavior of two-dimensional lattice materials.


Introduction
At the local scale, many periodic materials can be seen as made of interconnected beams. Such materials, also called lattice materials [1,2], have received a wide interest in the scientific and industrial community over the last years, mainly due to their high strength to weight and stiffness to weight ratios, which, of course, are an important advantage from a static point of view. Lattice materials, in fact, can be designed to be stretching-dominated, thus providing higher structural performances. Many scientific papers dealing with the characterization of their effective elastic properties are available in the literature [3][4][5][6][7]. In comparison, only a few papers focus on their dynamic behavior, which is of an identical practical interest in the engineering field. Innovative applications are, in fact, related to many aspects concerning the elastic wave propagation mechanism OPEN ACCESS within them. A small number of recently published works [8][9][10][11][12] has pioneered the study of the wave propagation phenomena in two-dimensional (2D) periodic lattice materials, with the main scope of detecting the existence of phononic band gaps, i.e. frequency ranges, within which the propagation of elastic/acoustic waves is prevented regardless of the incident wave direction. These studies look at a new generation of devices for energy absorption, noise and vibration control. Moreover, very interesting results have also been established within the field of civil engineering, where a new seismic mitigation strategy has been proposed that exhibits a close analogy with the idea of a lattice material [13].
In the context of lattice materials, the periodic topology plays a relevant role, considering the strict relation with the bending mechanism, which can be significantly dominant at the micro scale, where both the displacements and rotations can be similarly important [14]. This, in addition to experimental evidence reported by many authors [15][16][17][18][19], suggests an overcoming of the classical continuum theory, which assumes only displacements as kinematic quantities.
According to the classical continuum theory, stress-transfer locally occurs only through a force vector per unit area. As a consequence, stresses and strains can be represented by symmetric second-order real tensors. Nevertheless, for the accuracy of any prediction given via the classical theory, the condition that changes in stresses and strains (with wavelength λ) can be considered as uniform at the local scale is pivotal: D > λ l, with D denoting the global (structural) size, while l is a characteristic length representative of the microstructure. On the contrary, when dealing with a number of situations [20][21][22][23][24] (thin films, adhesive interfaces, notches, crack tips, localized deformations, boundary layers), it can be observed that the scale lengths of the problem agree with a new different condition (D > λ > l), rather than the previous one (D > λ l). More precisely, the lower the ratio λ/l > 1.0, the greater the importance of considering non-uniform stresses and strains at the local scale.
Unfortunately, due to the lack of an internal characteristic length, classical models are unable to capture the microstructure-dependent size-effect and, therefore, need to be extended by using higher order, non-local continuum theories.
In line with these last developments, the present work proposes the analysis of the microstructure-dependent size-effect on the dynamic behavior of 2D lattice materials.
Due to the very important conclusions outlined in [33,34], this influence is expected to be relevant also on the local resonance phenomena, which determine the existence of frequency band gaps within the low frequency region.
To this scope, a refined 1D finite element approximation, accounting for a micro-scale length parameter, as well as the shear deformability, has been proposed.

Frequency Gaps
Due to the complexity of the equations governing the wave propagation within a periodic material, there is no general rule that can be used to predict whether phononic band gaps exist before a full analysis has been performed.
The existence of phononic band gaps is commonly related to the destructive interference occurring in the multiple scattering and reflection of elastic waves, which propagate throughout a periodic material. It is well known that the position and width of the phononic band gaps depend on the lattice topology, as well as the elastic moduli of the material and the mass density. Nevertheless, the influence of the micro-scale constitutive parameters has yet to be analyzed in detail.
Recent experimental and numerical studies [36][37][38], however, show the existence of phononic band gaps in frequency ranges lower than those evaluable from the above cited destructive interference theory, revealing the link between these low frequency band gaps and the local resonance phenomena of the inner microstructures. Moreover, in a recent work [39], a modified 2D lattice material with square topology, obtained by connecting auxiliary cantilever beams to the primary micro-structure, was investigated. The frequency band diagrams, numerically identified, show as a result the relevant role played by the physical mechanism of local resonance.
Based on the previous considerations, it is the opinion of the authors that the study of the dynamic response of 2D lattice materials should require a more refined continuum approach, which accounts for the stiffness modification related to the size-effect. To this aim, a micro-scale length parameter has been proposed according to the hypotheses of the modified couple stress theory developed in [35].

Free Waves Propagation in 2D Periodic Lattice Materials
The study of the wave propagation within a periodic material generally involves the Bloch's theorem, which allows one to reduce the analysis of an infinite lattice to that of a single unit cell. An example of a simple two-dimensional square lattice material and the corresponding unit cell is depicted in Figure 1, where e 1 and e 2 are the unit vectors along the directions of spatial periodicity. The square topology implies in this case that e 1 and e 2 are normal to each other. In the example above (Figure 1), the microstructure is composed of interconnected beams arranged along the directions e 1 and e 2 only. In general, however, more complex microstructures can be investigated. For instance, in [39], a modified 2D lattice material has been considered, which exhibits a more elaborated unit cell.
According to the Bloch's theorem, if a plane elastic wave with angular frequency ω propagates in the lattice material (Figure 1), the displacement of an arbitrary point P is given by: In Equation (1), the symbol r indicates the position of the point P, k denotes the Bloch wave vector and u k (r) is the amplitude, which is characterized by the same spatial periodicity as the point lattice. As a consequence, it is possible to express the position of any point P as a function of the position of the corresponding point O, which is located in a reference unit cell. Once the reference cell has been selected (see Figure 1), it results: r = r 0 + n 1 e 1 + n 2 e 2 , where (n 1 , n 2 ) is an integer pair and r 0 is the position vector of O. It is then possible to update Equation (1) as follows: In brief, changes in wave amplitude from a generic cell to the adjacent one do not depend on the cell location within the point lattice. By virtue of this, it is possible to restrict the study of an infinite lattice system to the analysis of a unit cell.
It is easy to verify that the reciprocal unit vectors 1 e and 2 e assume the following form: where e 3 is the unit vector normal to the plane of the problem ( Figure 1). If the Bloch wave vector is expressed in the reciprocal space as follows: Then, from Equation (2), the following equation is obtained, which provides a periodic boundary condition for the dynamic analysis of the unit cell:

Finite Element Approximation
The analysis of the free wave propagation in 2D lattice materials can be performed by modeling a unit cell only. Rigid connections are usually considered between the micro-beams of which the unit cell is composed. Examples of the finite element (FE) mesh over the unit cell are shown in Figure 2. In this figure, both a simple 2D square lattice materials (I) and a more elaborated one (II) are sketched, the last one being composed of auxiliary cantilever beams interconnected to the primary microstructure. It is worth noting that this second example is substantially the same as reported in [39].
Regardless of how many micro-beams are interconnected to form the unit cell, the mesh here proposed is composed of finite elements characterized by six degrees of freedom Accounting for a general orientation of the finite element, it is useful to express the numeric vector e u , which refers to the local co-ordinates system ( 1 i , 2 i , 3 i ) as a function of the nodal displacements referred, instead, to the global co-ordinates system ( 1 e , 2 e , 3 e = 3 i ):  being the angle from 1 e to 1 i .
By virtue of the next position: where η indicates the following dimensionless quantity accounting for the influence of the shear deformations 12 the displacements field   T , , u v  , expressed in the local co-ordinates, can be interpolated over the generic finite element as follows: In Equation (10) with  indicating the normalized axial coordinate ( 1   , Figure 3) and pq h being four appropriate interpolating cubic functions (p = 1,2), (q = 0,1): The finite element proposed is able to account for shear strains and is locking-free. These features, as it is well-known, are essential for short transient and wave propagation analysis. A Bernoulli-Euler element, in fact, exhibiting an infinite phase velocity, because the parabolic order of the equation of motion, is useless for simulating wave propagation.
According to the non-local beam model presented in [35] (see Appendix), the generalized stresses within the generic finite element can be expressed as a function of the nodal displacements: where N, M, Y and V denote the generalized stresses of the beam model, C denotes the matrix given in Equation (A30) and B is a numeric matrix with the following entries: with reference to the global coordinates, the element stiffness matrix e K assumes the following form: On the other hand, the kinetic energy can be expressed as follows: where the consistent mass matrix e M is computed over the generic finite element by means of: By standard procedures, the equations of motion can be assembled in the following matrix form: where g M and g K are the overall mass and stiffness matrices of the unit cell and T  g  1  2  3  4  5 , , , , , , N  U U U U U U U  and F denote the nodal displacements and external forces vectors in the global reference system (  , 1 e , 2 e ), respectively, N being the number of nodes over the unit cell. If a plane elastic wave with angular frequency  is considered, Equation (27) assumes the new form:

 
The last form of the equations of motion (Equation (28)) must be coupled with the periodic boundary conditions given by Equation (6): where 1 U , 2 U , 3 U and 4 U are the nodal displacements relative to node 1 to 4, indicated in Figure 2 (connections with the four adjacent cells). By means of Equations (29) and (30) (31) where the matrix H works as a transfer operator: with: Finally, the equations of motion of the unit cell can be expressed in the following form: The condition  F 0 (free wave motion) implies that the study of the wave transmission within the periodic square lattice material reduces to the eigenvalue problem of Equation (35). For fixed values of 1 k and 2 k , to be expressed in the reciprocal space ( 1 e , 2 e ) according to Equation (5), the frequencies of the wave propagation coincide with the eigenvalues of the problem formulated in Equation (35). Varying 1 k and 2 k along the boundary of the irreducible part of the first Brillouin zone, the band structure of the lattice material can be detected. In Figure 4, the filled region (triangle OAB) indicates the irreducible part of the first Brillouin zone concerning the square topology under consideration.

Numerical Results
In order to elucidate the investigation approach, the finite element model presented in the previous section has been applied for studying the dynamic behavior of a 2D square lattice material. More in detail, a parametric analysis has been carried out by varying the role of the microstructure length l. The example under consideration concerns the scheme (I) shown in Figure 2, the geometry and the mechanical parameters being fixed, as indicated in Table 1.   ) implies that the size-effect is not present.
The eigenvalue problem given by Equation (35) has been solved via a call to a dedicated routine belonging to the IMSL Math Library. Thus, the frequency values have been identified and represented in a dimensionless form   by means of: where 1  denotes the first bending resonance frequency of a pinned-pinned beam with the same properties given in Table 1, its length being equal to the lattice constant a: In Figure   As it is easy to argue, the size-effect is found to be strongly relevant within the high frequency region when dealing with λ > 0.01.
For what concerns the dynamic behavior within the low frequency region, Table 2 indicates also in this case a considerable influence of the parameter  .
With respect to the corresponding predictions obtained by assuming λ = 0.00, relevant variations emerge, being the maximum percentage increase (+93%) found at point A when λ = 0.10.

Conclusions
As highlighted in a number of papers, the dynamic behavior of a square lattice material is significantly affected by the underlying microstructure. Due to the spatial periodicity, the hypothesis of infinite lattice points allows us to reduce the analysis to the unit cell and to investigate the dynamic behavior by means of the Bloch theorem.
It has been found by many authors that the natural frequencies predicted via a non-local beam model are always higher than those evaluated by classical beam models, due to the increased bending stiffness related to the so-called size-effect.
A more refined beam model is therefore required. For the scope of a reliable evaluation of the existence, position and width of frequency band gaps, in the opinion of the authors, the study of a 2D lattice materials should account for possible size-effects.
The paper provides a simple 1D finite element, developed in accordance with a simplified couple stress theory, which is able to simulate the size-effect via a unique microstructural parameter. Shear deformations are also accounted for.
The numerical results confirm the influence of the size-effect on the dynamic behavior of 2D periodic materials.

Appendix
In accordance with higher order theories, the local equilibrium conditions require the presence of couple stresses  m Mn (moment per unit area) in addition to the classical Cauchy stresses  t Tn , while with respect to a global continuum, equilibrium equations can be expressed as follows: In Equations (A1) and (A2), V indicates an arbitrary volume of a deformable body bounded by V  , the symbol n denotes the outer normal to V  , b and c denote, respectively, the body force and the body couple per unit mass, ρ is the mass density, x − x o denotes the distance vector of a generic point from an arbitrary pole O. Furthermore, n t and n c denote, respectively, the traction and couple (per unit surface) acting on the boundary of the body. The divergence theorem can be invoked in order to transform the 2D integrals in Equations (A1) and (A2): In Equation (A4), the symbol T w represents the axial vector related to the skew part of the stress tensor T. Since the volume V is arbitrary, the volume dependence can be discarded, leading to the following equilibrium equations: Due to the presence of the couple stresses, Equation (A6) implies that the Cauchy stress tensor T is not symmetric. Then, it can be decomposed as follows: Furthermore, it results: On the other hand, at the local scale, the deformation measures are represented by the following two tensors: The symbols Γ and Κ denote the (small) strain and curvature tensor, respectively. It is easy to verify that the following relationships exist: where i u and i θ denote the generic components of the displacement field and the micro-rotation field (i = 1,2,3). The symmetric part of Γ in Equation (A15) denotes the classical infinitesimal strain tensor, while the skew part in Equation (A16) accounts for the difference between the macro rotations j k i j k i u 1 ω = e 2 x  and the micro rotations k θ (k = 1,2,3). As a consequence, the virtual work done by the internal stresses assumes the following final expression: It is worth noting that if no difference between the micro-rotations and macro-rotations occurs, the strain tensor A vanishes and the continuum model reduces to couple stress theory. In this case, the anti-symmetric tensor L does not contribute to the work done by the internal stresses. More precisely, the derivation of the couple stress theory [29] from the general Cosserat theory [25] has been investigated in [32], where a dimensionless number drives the transition. Unfortunately, the definition of both the Young modulus and the shear modulus fall in defect if the couple stress theory is interpreted as a special case of the micro-polar (Cosserat) theory.
Whereas, the link from the classical couple stress theory to the simplified form proposed in [35] is founded upon the following additional equilibrium equation for the moment of couples: In Equation (A29), N, M, Y and V denote the generalized stresses of the beam model (axial force, bending moment, additional bending moment due to the size-effect, shear force), while u x