Next Article in Journal
Off-Shell Quantum Fields to Connect Dressed Photons with Cosmology
Next Article in Special Issue
Steady Solitary and Periodic Waves in a Nonlinear Nonintegrable Lattice
Previous Article in Journal
In-Situ Energy Dispersive X-ray Reflectivity Applied to Polyoxometalate Films: An Approach to Morphology and Interface Stability Issues in Organic Photovoltaics
Previous Article in Special Issue
A Cosserat Model of Elastic Solids Reinforced by a Family of Curved and Twisted Fibers
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Failure of Classic Elasticity in Predicting Elastic Wave Propagation in Gyroid Lattices for Very Long Wavelengths

by
Giuseppe Rosi
1,*,
Nicolas Auffray
2 and
Christelle Combescure
2
1
CNRS, MSME, Univ Paris Est Creteil, F-94010 Creteil, France
2
CNRS, MSME, Univ Gustave Eiffel, F-77447 Marne-la-Vallée, France
*
Author to whom correspondence should be addressed.
Symmetry 2020, 12(8), 1243; https://doi.org/10.3390/sym12081243
Submission received: 28 April 2020 / Revised: 5 June 2020 / Accepted: 1 July 2020 / Published: 28 July 2020
(This article belongs to the Special Issue Recent Advances in the Study of Symmetry and Continuum Mechanics)

Abstract

:
In this work we investigate the properties of elastic waves propagating in gyroid lattices. First, we rigorously characterize the lattice from the point of view of crystallography. Second, we use Bloch–Floquet analysis to compute the dispersion relations for elastic waves. The results for very long wavelengths are then compared to those given by classic elasticity for a cubic material. A discrepancy is found in terms of the polarization of waves and it is related to the noncentrosymmetry of the gyroid. The gyroid lattice results to be acoustically active, meaning that transverse waves exhibit a circular polarization when they propagate along an axis of rotational symmetry. This phenomenon is present even for very long wavelengths and is not captured by classic elasticity.

1. Introduction

Architectured materials are those that possess an inner geometry [1]. This multi-scale spatial arrangement of the constitutive materials allows for achieving mechanical properties that are not present in bulk material itself [2]. Although this appears to be an engineering-based approach to materials design, it should be noted that this strategy is, in fact, central in nature where biomaterials must perform many functions from a small and limited set of elementary chemical elements [3,4]. Therefore, to enhance some target properties, regular patterns often emerge. The best-known example is the honeycomb, where bees need to maximize the volume of the cells while minimizing the quantity of matter (wax) used [5]. Another example is the iridescent color of the wings of some butterflies. This phenomenon is due to the non-centrosymmetric mesostructure of the material constituting the wings, which acts a photonic crystal [6,7].
The study of elastic waves propagating in architectured materials is of particular interest since unconventional effects due to the local organization of the matter can emerge on a macroscale. In order to study these phenomena adequately, two points of view can be adopted. Either to describe all the details of the architecture or to consider an effective continuum as replacement. The first option is very general since no particular modeling assumptions are involved. However, since the inner geometry of the material has to be explicitly described and meshed, the numerical cost is often prohibitive for actual applications. Moreover, the computed solution often contains many unnecessary details for practical use. The second option, which is based on elastodynamic homogenization [8,9,10] amounts to substitute the initial heterogeneous material by an equivalent homogeneous continuum. This equivalence is only valid under specific assumptions on the range of variation of some intrinsic parameters, and hence more restrictive than the first approach. However, within the validity domain of the method, the physics, up to a certain order is correctly described. Indeed, an effective theory is a reduced model obtained by filtering the actual physics so as to retain, in the continuum formulation, only the most prominent effects. Depending on the targeted applications, the effective model can be of different degrees of richness. This results in a fairly important reduction in the computational cost, which is very interesting for optimizing an architectured material, since in that case the numerical model has to be computed many times along the process.
For infinite Periodic Architectured Materials (PAM), which are the subject of this paper, the condition under which the complete wave problem can be substituted by an effective one relies on the ratio ( η ) between the size of the periodic unit cell (L) and the wavelength ( λ ) of the mechanical field. When this ratio approaches zero, the classical Long-Wavelength (LW) approximation is obtained and, provided that the frequency ω is also low (Low-Frequency (LF) approximation), the heterogeneous material can be replaced by a classical effective continuum. This situation, which has been well investigated, is completely contained in what is called LF-LW elastodynamics homogenization [8,9,10]. When the equivalent medium is a classical continuum (Cauchy continuum), the effective behavior is non-sensitive to certain features of the inner geometry such as noncentrosymmetry, chirality, or a n > 4 -fold axis of rotational anisotropy [11,12,13].
Now, when the scale separation ratio η is small, but not vanishingly small, elastic waves propagating through the matter interact with the inner architecture. In this situation, several propagation quantities, such as the phase and group velocities, become frequency dependent and the wave propagation is dispersive [14]. Non-standard dependencies on the architecture, that were left over in LF-LW approximation, may thus appear. These situations, which are outside the frame of standard elastodynamic homogenization, can nevertheless be modeled if the Cauchy equivalent continuum is replaced by a generalized continuum [15,16,17]. In this work we focus on bulk propagation, however it is important to notice that effects near boundaries, such as surface waves [18,19], are also of particular interest.
Wave propagation in non-centrosymmetric or chiral materials, the two concepts being distinct (see Appendix A for a dictionary of point groups and their associated properties), has been a subject of interest among physicists for centuries, mainly in the field of optics and electromagnetism. The first experiments showing the interaction of light with chiral molecules like sugar goes back to the beginning of the 19th century [20]. The effect which is associated to electromagnetic waves propagating in non-centrosymmetric crystals is the rotation of the plane of polarization when the wave propagates along an optical axis, i.e., an axis of rotational symmetry. The rotation is due to the decomposition of the linearly polarized transverse wave into two circularly polarized waves with opposite handedness and different phase velocities [21]. This phenomenon is known under the name of “optical activity”. The analogue of this effect can be observed for elastic waves and is known as “acoustical activity” [21]. It is interesting to remark that optically-active crystals are also found to be acoustically active and that, as it will be shown in this paper, this effect can occur also in the LF-LW regime.
Recently, the interest in investigating the properties of materials based on chiral and non-centrosymmetric architectures has grown. To this end it is important to point out that chirality and noncentrosymmetry are not equivalent, and that their impact on the physics of the problem can be different.
The set of transformations that let the unit cell of an architectured material invariant constitutes its symmetry group. The material is said chiral if its symmetry group contains only rotations and it is said to be centrosymmetric when it contains the inversion [22]. It is important to observe that in a 2D space the inversion is a rotation (preserving the material orientation) while in 3D, it is a transformation reversing the material orientation. Since the nature of the inversion depends on the dimension of the space, the implication between chirality and centrosymmetry are not the same in 2D and 3D. In 2D, the chirality and centrosymmetry are independent [15], while in 3D chiral materials are necessarily non-centrosymmetric [23].
Several works can be found in the literature that investigate 2D chiral elasticity and focus on their unusual mechanical properties such as negative Poisson ratio [24,25]. Concerning wave propagation, these architectures have been extensively studied [26,27,28] and the need for a generalized continuum theory in order to capture the onset of dispersion and anisotropy at higher frequencies has been pointed out [15,16]. In all these cases, the unit cells under investigation are chiral and centrosymmetric. Again, it is worth noting that such a combination is possible pnly in 2D. When moving to 3D, the picture becomes more complex and due to habits coming from the 2D situation, an ambiguity between the two definitions can be usually found in the literature. For instance, the well-known and studied in-plane hexachiral and tetrachiral patterns [29,30], are no more chiral once extruded in 3D.
Concerning wave propagation in non-centrosymmetric materials, if the phenomenon can be studied in 2D [15], the effects become even more interesting in 3D as the polarization of waves are then richer. As a consequence, interest in 3D non-centrosymmetric metamaterials have recently emerged for their strech-twist coupling [31] or their acoustical activity [32]. The effects related to size-dependent properties and characteristic lengths are also exploited and investigated in [33] and a micropolar generalized model is used to investigate acoustical activity in [32].
The present work is focused on the features of elastic wave propagating in non-centrosymmetric architectured materials. Among them, those based on gyroid architectures are probably the most commonly used. In electromagnetics they are widely studied as metamaterials [6], in acoustics as phononic crystals [34], and in biomechanics as bone substitutes [35,36]. In this paper, we will highlight a particular situation for which the solution predicted by classic continuum mechanics is wrong even for very long wavelengths. It is important to note that the sensitivity of the mechanical behavior to the lack of centrosymmetry can also manifest in statics [37,38,39].
The paper is organized as follows: In Section 2 the gyroid lattice is described. In Section 3 the Bloch–Floquet analysis is introduced along with some necessary definitions for polarization studies. Dispersion analysis is performed and discussed in Section 4. Section 5 compares the results from Section 4 to those obtained in the LW-LF approximation. Finally, some conclusions are drawn in Section 6.

Notations

Throughout this paper, the Euclidean space E 3 is equipped with a rectangular Cartesian coordinate system with origin O and an orthonormal basis B = { e 1 , e 2 , e 3 } . Upon the choice of a reference point O in E 3 , the Euclidean space and its underlying vector space E 3 can be considered as coincident. As a consequence, points will be designated by their vector positions with respect to O . For the sake of simplicity, E 3 will, from now on, simply be denoted E . In the following, r will designate the position vector of a point P , and, with respect to B ,
r = x e 1 + y e 2 + z e 3 .
When needed, Einstein summation convention is used, i.e., when an index appears twice in an expression, it implies summation of that term over all the values of the index. The dot operator ( · ) stands for the scalar product, the ∧ for the cross product, and δ i j is the Kronecker delta.
Moreover, the following convention is retained:
  • Blackboard fonts will denote tensor spaces: T ;
  • Tensors of order > 1 will be denoted using uppercase Roman Bold fonts: T ;
  • Vectors will be denoted by lowercase Roman Bold fonts: t .
The orthogonal group in R 3 is defined as O ( 3 ) = { Q GL ( 3 ) | Q T = Q 1 } , in which GL ( 3 ) denotes the set of invertible transformations acting on R 3 .

2. The Gyroid Lattice

The gyroid is a triply periodic minimal surface introduced for the first time in [40,41]. Since it is a minimal surface, it has a zero mean curvature, meaning that every point on the surface is a saddle point with equal and opposite principal curvatures [42]. It is periodic with respect to three orthogonal space vectors and is chiral, meaning that the surface only possesses rotation symmetry elements or, equivalently, that it does not possess any symmetry plane nor symmetry center [43].

2.1. Parametrization of the Gyroid Lattice

The gyroid’s morphology is usually described using a level surface given by the following equation:
ϕ ( x , y , z ; a , b ) = 0
with a and b real parameters and x , y , z coordinates of the position vector r .
In this paper, we will focus on a particular gyroid defined by:
ϕ ( x , y , z ; a , b ) : = sin 2 π a x cos 2 π a y + sin 2 π a y cos 2 π a z + sin 2 π a z cos 2 π a x b ,
that exists only in the range | b | < 2 [43]. Indeed, beyond this peculiar value of b, it is possible to show that the surface described by ϕ presents discontinuities located at the borders of the fundamental (or asymmetrical) unit cell. A proof of this geometric constraint is provided in Appendix C. Due to its chiral nature, the gyroid surface exists in two enantiomorphic forms: dextrogyre and levogyre. The surface described by the implicit Equation (1) will be, arbitrarily, chosen to be the dextrogyre form. The levogyre form of the implicit equation is easily obtained by applying, for instance, the transformation x x in Equation (1). From the definition of the surface, one can then obtain a volume by defining the presence of matter for points satisfying the following inequality:
ϕ ( x , y , z ; a , b ) > 0 .
The parameter a controls the spatial period while the parameter b controls the porosity p, defined as the unity minus the ratio between the volume of the gyroid lattice and that of the unit cell. Examples of unit cells of such solids obtained with different values of b are plotted in Figure 1.
The relationship between the porosity p and the parameter b is not analytical but can be estimated numerically. As can be observed in Figure 2 this relationship is almost linear and for porosities between 0.2 and 0.8, the following linearized formula can be used to estimate the porosity:
p = 0.325 b + 0.5 .

2.2. Symmetry Properties

Thanks to its triply periodic nature, the gyroid structure can be considered as a crystal [34] with Body Centred Cubic (BCC) Bravais lattice and point group O , using group theoretic notation [44], or 432 using the Hermann–Maugin notation. It should be noted that this cubic group only contains rotations geometry of the gyroid is hence chiral. To be more specific there are 3 different cubic point groups: O , O , and O Z 2 c . The first one just contains rotations, the group is hence chiral and non-centrosymetric. The second, O , possesses symmetry planes but not the inversion, the group is achiral and non-centrosymetric. The last group O Z 2 c is centrosymmetric hence achiral. Some details are provided in the Appendix A, and more information can be found in [23]. The symmetry of the spatial structure is described by the space group, which details how transformations from the Bravais lattice and the point group are combined in the actual crystal. The space group ( S G ) of the gyroid crystal is, using Hermann–Maugin notation, I 4 1 32 (space group #214 in the International Tables of Crystallography [45]), where the I stands for Body-Centered (BC), meaning that the conventional unit cell defined in crystallography is not primitive, but body-centered (more details provided in Section 2.3). This space group contains screw axes and, as such, is not symmorphic. A space group is called symmorphic if, apart from the lattice translations, all generating symmetry operations leave one common point fixed. Permitted as generators are thus only the point-group operations: Rotations, reflections, inversions, and rotoinversions. The symmorphic space groups may be easily identified because their Hermann–Mauguin symbol does not indicate a glide or screw operation.
If C stands for the “crystal” structure, and ⭑ for the group action as defined in Appendix B:
r C , g S G , r = g r C .
From the generating transformations defined in Appendix B and using the equation of the gyroid surface (c.f. Equation (1)) it is straightforward to verify that:
g S G , ϕ ( g r ) = ϕ ( r ) .

2.3. Unit Cell

Due to the periodicity of the geometry, the study of the gyroid structure can be restricted to a unit cell:
r E , ( r 0 , t ) T , R , r = r 0 + t
where T designates a unit cell and R a periodicity lattice. Note that for a given lattice R , the choice of T is not unique. It can be convenient to chose a reference unit cell from which other unit cells will be defined using lattice vectors a i and a triplet ( n 1 , n 2 , n 3 ) combined as follow: t = n 1 a 1 + n 2 a 2 + n 3 a 3 , n i Z . The triplet ( 0 , 0 , 0 ) is then associated to the reference unit cell. Once a lattice basis chosen, the considered unit cell is defined as follows,
T C = { r 0 E | r 0 = r a 1 + s a 2 + t a 3 , ( r , s , t ) [ 0 , 1 [ 3 }
and the associated periodicity lattice is given by:
R C = { t E | t = n 1 a 1 + n 2 a 2 + n 3 a 3 , n i Z } .
These geometrical sets can be described using lattice vectors a i , gathered into a basis as B = { a 1 , a 2 , a 3 } . Note that they can as well be described with respect to the basis B = { e 1 , e 2 , e 3 } of E .
Among all the possible unit cells, some are special and have been given a standard name in the crystallography community: The Conventional Unit Cell (CUC) and Primitive Unit Cell (PUC).
The CUC of the a BCC lattice is depicted in Figure 3. It is defined as the smallest cell having its edges along the symmetry directions of the Bravais lattice. Notice that, for BC lattice, this unit cell is not minimal and a so-called PUC can be considered instead. For a continuous structure tilling of the space, the primitive unit cell is defined as the smallest tile that generates the whole tiling using only translations. As such the primitive unit cell is a fundamental domain with respect to translational symmetries only.

2.3.1. BCC Conventional Unit Cell

For a BCC lattice, the conventional unit cell is defined as depicted in Figure 3. As its faces are perpendicular to Bravais lattice directions, despite its non minimality, this unit cell is easy to use for numerical computations.
In this case, the conventional lattice vectors a i , are chosen such that a i e i = 0 .

2.3.2. BCC Primitive Unit Cell

For a BCC lattice, two possible PUC are represented in Figure 4.
The primitive lattice vectors a i are not unique and the ones for the PUC depicted in Figure 4a are defined as:
b 1 = a 2 ( e 1 + e 2 e 3 ) , b 2 = a 2 ( e 1 + e 2 + e 3 ) , b 3 = a 2 ( e 1 + e 2 + e 3 ) ,
while the ones presented in Figure 4b are defined as:
b 1 = a e 1 , b 2 = a e 2 , b 3 = a 2 ( e 1 + e 2 + e 3 ) .
They form P = { b i } 1 i 3 and P = { b i } 1 i 3 two other bases of E 3 , which metric tensors are given by:
g ( P ) = ( b i · b j ) = a 2 4 3 1 1 1 3 1 1 1 3 , g ( P ) = ( b i · b j ) = a 2 4 4 0 2 0 4 2 2 2 3 .
Being defined by a more symmetrical set of vectors P , only the first primitive unit cell will be considered here after.

2.4. Reciprocal Basis and Brillouin Zone

The vector space dual to E is symbolized by E , and is formally defined as the space of linear forms on E ,
l E , u E , l ( u ) = α R .
Upon the choice of a scalar product the two spaces can be identified:
( l E , v E ) , u E , l ( u ) = v · u = α R
and from a basis of E a basis of E can be constructed. In the field of physics, E corresponds to the space of wavevectors, and a generic element of E is denoted by k . For our applications, it is fundamental to introduce the reciprocal lattice R of R :
R = { ξ E | ξ = ξ 1 a 1 + ξ 2 a 2 + ξ 3 a 3 , ξ i Z } .
The vectors ( a 1 , a 2 , a 3 ) constitute the lattice basis B of R and verify:
a i · a j = δ i j , where δ i j = 1 if i = j , 0 if i j .
B can be computed from any lattice basis ( a 1 , a 2 , a 3 ) of R according to:
a 1 = a 2 a 3 a 1 · ( a 2 a 3 ) , a 2 = a 3 a 1 a 2 · ( a 3 a 1 ) , a 3 = a 1 a 2 a 3 · ( a 1 a 2 ) .
Due to the following property,
( t , ξ ) R × R , e 2 π i ξ · t = 1
vectors of the reciprocal lattice are the supports of R -periodic functions on E since,
( r 0 , t ) T × R , f ( r 0 + t ) = ξ R λ ξ e 2 π i ξ · ( r 0 + t ) = ξ R λ ξ e 2 π i ξ · r 0 e 2 π i ξ · t = f ( r 0 ) .
In addition, use will be made of the First Brillouin Zone (FBZ) T of the reciprocal lattice R defined as:
T : = { k E | ξ R 0 , k < k ξ }
Using the reciprocal lattice R and the FBZ T , any wavevector k can be expressed as:
k E , ( k 0 , ξ ) T , R , k = k 0 + ξ .
We can geometrically interpret T as the set of wavevectors k which are closer to the null wavevector than to any other wavevector ξ of the reciprocal lattice R . It is the Wigner–Seitz cell of the reciprocal lattice, this cell is uniquely defined and independent of the choice of T . Similarly to the primitive unit cell in the direct space, the FBZ is a fundamental domain with respect to translations.
Physically, the wavelength λ is defined as the inverse of the wavenumber, which is the norm of the wavevector: λ = 1 / k . Then, wavevectors belonging to T have wavelengths that are greater than the wavelength of the periodicity lattice. When k 0 the wavelength becomes infinite, solicitations varying with almost null wavenumber are said to be scale separated with respect to the periodicity lattice. This is usually the regime in which the LW approximation of elastodynamics homogenization holds.
The FBZ can be further reduced if we consider also the symmetry operations of the point group. The result is an Irreducible Brillouin Zone (IBZ), that is delimited by points of high symmetry, summarized for the considered gyroid lattice in Table 1. In this table, the high symmetry points are given in the non-orthogonal reciprocal lattice basis P , dual to the primitive lattice basis P as well as in the orthonormal lattice basis B which coincides with its dual in the reciprocal space. The path obtained connecting these high symmetry points along the edges of the IBZ is often used to characterize the photonic and phononic properties of the lattice [6,26]. However, it has been pointed out that this choice is not always reliable as some relevant information, e.g., about band gaps, could be missing [46].
In this paper we consider the basis P , that corresponds to the one depicted in Figure 4a. The reciprocal lattice is itself a Bravais lattice, and in the case of BCC lattice, it is a Face Centered Cubic (FCC) lattice. Using Equation (6), the reciprocal basis P is equal to:
a 1 = 1 a ( e 1 + e 2 ) , a 2 = 1 a ( e 2 + e 3 ) , a 3 = 1 a ( e 1 + e 3 ) .
The metric tensor of P is the inverse of the one of P :
g ( P ) = g ( P ) 1 = 1 a 2 2 1 1 1 2 1 1 1 2 .

3. Analysis Tools

In the previous section we characterized the lattice from the point of view of crystallography. In this next section we will use these results to compute the elastodynamic response of the lattice. The objective of the present section is to provide the analysis tools to be used to perform the computation and to interpret the results.

3.1. Bloch–Floquet Analysis

Since the material is periodic, the dispersion diagram will be computed using Bloch–Floquet analysis [47]. The elastodynamics equation for the periodic continuum reads:
div C ( r ) : u ( r ) ̲ = ρ ( r ) u ¨ ( r )
where ρ ( r ) is the R -periodic mass density and C ( r ) is the R -periodic fourth-order elasticity tensor. As we saw in Section 2.3, each cell of the assembly can be identified by the triplet ( n 1 , n 2 , n 3 ) , where the triplet ( 0 , 0 , 0 ) is conventionally assigned to the reference unit cell. The position of a point r of the ( n 1 , n 2 , n 3 ) -cell is obtained from the position of a point in the reference unit cell r 0 by Equation (4) where t = n p a p . Being R -periodic, ρ ( r ) and C ( r ) verify:
( r 0 , t ) T × R , ρ ( r 0 + t ) = ρ ( r 0 ) , C ( r 0 + t ) = C ( r 0 )
Thanks to the Bloch–Floquet theorem [47], elementary solutions to the Equation (9) over C can be searched for in the form of Bloch-waves:
u k ( r 0 ) = U k ( r 0 ) e 2 π i f t k · r 0 , U k C 3 ,
where U k is the complex polarisation vector which is R -periodic in space and constant in time, f is the frequency of the Bloch-wave, and k its wavevector U k describes the movement of matter as the wave propagates. It is important to remark that wavevectors k follow the so-called “crystallographer’s definition” which consists in dropping the often seen 2 π coefficient. This implies, for instance, that the norm k , i.e., the wavenumber, is directly the inverse of the wavelength λ , which is more convenient for physical interpretation of the results. In the case of an homogeneous material the polarization vector becomes constant in space and the classical plane wave solution is retrieved. Since the displacement field in Equation (10) is complex valued, its real part should be computed in order to retrieve the physical solution.
From its definition as a Bloch-wave, the displacement at a point r image of the r 0 T by a translation t R has the following expression:
u k ( r ) = u k ( r 0 + t ) = U k ( r 0 + t ) e 2 π i f t k · ( r 0 + t ) = u k ( r 0 ) e 2 π i k · t .
The physical meaning is that the displacement vector at two homologous points (two points r 1 and r 2 are said homologous if r 1 r 2 R ) only differs by a phase factor.
Additionally, the Bloch-wave expression in Equation (10) has the interesting property to be also R -periodic. Indeed, due to its R -periodicity, U k ( r 0 ) can be decomposed as a Fourier series, leading to the equivalent expression for u k ( r 0 ) :
u k ( r 0 ) = ξ R U ˜ k + ξ e 2 π i f t ( k + ξ ) · r 0
where U ˜ k + ξ stands for the Fourier coefficients of the series expansion. Using this particular form, the R -periodicity of the Bloch-waves is easily proven using the change of variable ξ ˜ = ξ + ξ :
u k + ξ ( r 0 ) = ξ R U ˜ k + ξ + ξ e 2 π i f t ( k + ξ + ξ ) · r 0 = ξ ˜ R U ˜ k + ξ ˜ e 2 π i f t ( k + ξ ˜ ) · r 0 = u k ( r 0 ) .
The main consequence of this property is that the characterization of the elastodynamics behavior of a periodic material does not require to investigate the mechanical response to all the k E but can be restricted to the study of k T via R -periodicity of the wavevector, where T corresponds to the First Brillouin Zone (FBZ), as introduced in Section 2.3.

3.2. Polarization of Waves in Homogeneous Materials

Before presenting the results, it is useful to recall some definitions concerning the polarization of elastic waves in homogeneous materials. Let us go back to the Bloch-wave ansatz introduced in Equation (10), since the material is now considered homogeneous, the polarization vector U k is constant both in space and time. In the most general case, the complex polarization vector U k , that will be denoted U from now on for the sake of simplicity, can be decomposed in its real and imaginary parts as follows:
U = U R + i U C .
An interpretation of this decomposition, and of its consequences on wave propagation, can be obtained by considering the real part of Equation (10):
u ^ = Re ( u ) = U R cos 2 π ( f t k · r ) U C sin 2 π ( f t k · r ) .
Since the vectors U R and U C are independent, the polarization of the displacement can be very rich. Its precise nature is directly related to conditions on U R and U C , as summarized in Table 2. It is important to remark that different conventions are used to define the handedness of circularly polarized waves. In this paper, we will consider that a wave is right handed if it follows the curl of the fingers of a right hand whose thumb is directed towards the wave propagation, away from the source. In the table, the unit normal vector defining the direction of propagation is defined by n = k k .
A complex polarization vector can lead to a phase shift between the components of the displacement vector u ^ , and thus to a polarization that is other than linear, see Figure 5 for illustration of a linear and two circular polarizations with opposite handedness.

4. Dispersion Analysis Using Finite Elements Analysis (FEA)

In order to investigate the ultrasonic properties of the gyroid lattice, and given the periodicity of the architecture as described in the previous sections, an approach based on Bloch–Floquet analysis will be followed. For the sake of simplicity, the conventional unit cell depicted in Figure 6 is used to define the numerical model. For the investigation of the elastodynamic properties of the gyroid crystal, the wavevector will be restricted to the boundaries of the Irreducible Brillouin Zone (IBZ), as depicted in red in Figure 7a. The high symmetry points of this IBZ are defined in Table 1. The model has been implemented using the commercial software Comsol Multiphysics and considering titanium as constitutive material, the parameters of which are displayed in Table 3. The mesh of the unit cell is presented in Figure 6 and consists of 66,232 tetrahedral elements. Lagrange quadratic elements are used, for a total of 329,277 degrees of freedom. Periodic Bloch–Floquet conditions are implemented by imposing them as displacement conditions at the boundaries, following Equation (11). Then, the wavenumber in k is imposed and the corresponding frequencies are retrieved by solving the eigenvalue problem. The computation of each wavenumber took an average of 109 seconds on a workstation equipped with an Intel(R) Xeon(R) CPU E5-1650 v2 at 3.50 GHz using six cores.
The results of the dispersion analysis are depicted in Figure 7a. It can be observed, qualitatively, that these results are similar to those obtained for electromagnetic waves in [6] (see Figure 8). In particular, the behavior of the acoustic branches, i.e., those branches starting from the origin Γ , corresponding to transverse waves (gray lines in Figure 7a) is remarkably similar.
Since the objective of the paper is to investigate the behavior of the lattice within the LW-LF approximation given by classic continuum mechanics, the phase velocities and polarization of waves have been computed for large values of the wavelength with respect to the size of the unit cell. In this case, k = 16 . 7 m 1 , which corresponds to a wavelength close to 60 times the size of the unit cell. For each mode, the polarization vector has been estimated by computing the average of the complex displacement of the eigen-mode over the unit cell, and the results are listed in Table 4. We will now analyze the results along the following directions of propagation, also depicted in Figure 7b:
  • [ 001 ] : This direction is going from the center of the fundamental cell to the middle of a face. It corresponds to an axis of rotation of order 4 (rotations of π / 2 rad);
  • [ 011 ] : This direction is going from the center of the fundamental cell to the middle of an edge. It corresponds to an axis of rotation of order 2 (rotations of π rad);
  • [ 111 ] : This direction is going from the center of the fundamental cell to a vertex. It corresponds to an axis of rotation of order 3 (rotations of 2 π / 3 rad).
Using the conditions listed in Table 2, we have characterized the polarization for each of the above propagation directions. The results are summarized in Table 4.
Along the direction [ 001 ] a longitudinal wave propagating at 3018.5 m/s can be observed. The two transverse waves have eigenmodes with complex amplitude and propagate with close phase velocity that tends to a common value of of 1932.2 m/s for infinite wavelengths, as it can be also observed in Figure 9. These complex amplitudes correspond to two circularly polarized transverse waves with opposite handedness. In the direction [ 011 ] , a longitudinal wave propagating at 3249.8 m/s can be observed. One transverse wave is linearly polarized in direction [ 100 ] , and propagates with a phase velocity of 1931.7 m/s. The last solution corresponds again to a transverse wave, linearly polarized along ( 0 , 1 , 1 ) , with velocity 1510.9 m/s. Finally, direction [ 111 ] has a linearly polarized longitudinal wave at 3322.8 m/s, and two circularly polarized waves with opposite handedness and propagating with close velocity, converging to 1663.4 m/s.
In summary, circularly polarized waves exist only if the direction of propagation is along a rotation axis of symmetry of order greater than 2. Moreover, as can be seen in Figure 9, for both [ 001 ] and [ 111 ] directions, the circularly polarized waves with opposite handedness propagate with the same phase velocity only in the infinite wavelength limit, and they start to diverge as the wavenumber increases. In particular, for direction [ 001 ] the right handed wave becomes slower than the left handed one, while the opposite phenomenon can be observed for direction [ 111 ] . This is due to the chirality of the unit cell. As one can notice, phase velocities are different even for a very large wavelength compared to the size of the unit cell, i.e., λ / a 10 . Since only the phase velocity is affected, and not the amplitude, this effect can be interpreted as the elastic equivalent of circular birefringence in optics. This means that if a linearly polarized wave passes through a gyroid lattice, the polarization plane of the incident wave will be rotated. This is due to the phase difference (retardance) between the two circular components, which produces a rotation of the polarization plane. The concept is illustrated in Figure 10. Moreover, since phase velocity is involved in reflection of waves at boundaries via the Snell–Descartes law, and in particular in the definition of the Brewster angle of total reflection, gyroid lattices show the potential for being used as elastic circular polarizing filters.
Furthermore, the overall dispersion is normal for direction [ 001 ] , i.e., phase velocity decreases when increasing the frequency, and anomalous for direction [ 001 ] (see [48] for the definition). The anisotropy of the material and the dispersive properties could also have consequences on surface and guided waves propagating in presence of boundaries [49,50], as well as in reflection/transmission problems [51]. These effects will be investigated in further works.

5. Long-Wavelength and Low-Frequency Approximation and Classic Elasticity

In this last section we will introduce and identify the equivalent homogenized model in the framework of classic linear elasticity. This equivalent homogenized model is characterized by a couple of effective tensors ρ and C H in such a way that the displacement field v is solution of the following problem:
div C H : v ( r ) ̲ = ρ v ¨ ( r )
where v verifies < u > = v , < . > denotes the spatial average operator over T and u is the displacement field solution to the heterogeneous problem Equation (9), as done for instance in [10].
Since the effective continuum is homogeneous, we consider a plane wave solution with k = f / c n where c is the phase velocity of the wave and n a unitary vector. The substitution of this wave solution and of a linear elastic constitutive law into Equation (12) leads to following equation:
Γ · U = ρ c 2 U
where Γ = n · C H · n is the Christoffel, or acoustic, tensor. The solution of the eigenvalue problem stated in Equation (13) for a given direction n gives the phase velocities and polarizations of waves in the effective continuum.
In classical elasticity, a material with cubic symmetry is defined by three independent material constants. Using Mandel notation [52], the corresponding elastic tensor for a material having its symmetry axis parallel to B reads:
[ C ˜ H ] = c 11 c 12 c 12 0 0 0 c 12 c 11 c 12 0 0 0 c 12 c 12 c 11 0 0 0 0 0 0 c 44 0 0 0 0 0 0 c 44 0 0 0 0 0 0 c 44 B .
It is worth noting that classical elasticity is insensitive to the lack of centrosymmetry [53]. The symmetry class of the elasticity tensor in the cubic system is O Z 2 c meaning that even if the material symmetry group of the unit cell does not contain the inversion, the symmetry group of elasticity tensor will inherit it.
In the case of cubic materials, the solutions of Equation (13) listed in Table 5 are directly obtained.
Using the phase velocities computed from the Bloch–Floquet analysis, parameters c 11 , c 12 , c 44 , and ρ can be identified and the homogeneous equivalent properties listed in Table 6 are then deduced. It can be noticed that, as presented in Table 5, since we considered the propagation along the rotational axes of symmetry, for each direction we observe a purely longitudinal wave and two purely transverse waves.
We now move on to comparing phase velocities and polarizations obtained from the Bloch–Floquet analysis with the ones forecast by the Long-Wavelength and Low-Frequency approximation using classical elasticity. We start with direction [ 001 ] . As already mentioned, this direction corresponds to a rotational axis of symmetry of order 4. In this case, as the elasticity tensor is non sensitive to chirality, the symmetry group of the physical phenomenon (the symmetry group of the physical phenomenon is the intersection of the symmetry group of the constitutive equations and the symmetry group of the mechanical solicitation) is conjugate to D 4 Z 2 c . Indeed, as the acoustic tensor defined in Equation (13) is a second-order tensor, the Hermann theorem of Crystal physics [11] predicts its behavior as transversely isotropic, i.e., O ( 2 ) Z 2 c . As a consequence Γ must have an eigenspace of dimension 2. All the directions of wave propagation for which this is verified are called acoustic axis of the material. Moving to the results presented in Table 5, we can see that the classic theory indeed predicts one faster longitudinal wave and two slower transverse waves propagating with the same phase velocity.
In the previous section we saw that Bloch–Floquet analysis identifies these waves as circularly polarized with opposite handedness, which is of course equivalent. Indeed, even for very large wavelengths, the numeric evaluation of the polarization provided by Table 4 corresponds to the eigensystem (eigenvalues plus eigenvectors) of a Hermitian acoustic tensor. In such case, the space corresponding to the double eigenvalue (which is real due to Hermitian nature of the matrix) is two dimensional over the complex field C . However, in the case of a classic Cauchy continuum the acoustic tensor is symmetric real, and the eigen-space corresponding to the double eigenvalue is two dimensional over the real field R . Since the 2D vector space over C can be considered as a four dimensional vector space over R the span is not equivalent. Moreover, as presented in Section 4, when the ratio between the wavelength and the size of the unit cell becomes finite, this 2D eigenspace splits into two 1D eigenspaces with different phase velocities. It is important to notice again that this effect occurs in the LF-LW regime, where the classical elastodynamic homogenization is supposed to hold, or give at least approximated results while preserving the physics of the problem. Similar results are obtained for propagation along [ 111 ] , which corresponds to the rotational axis of symmetry of order 3. In this case the symmetry group of the physical phenomenon is conjugate to D 3 Z 2 c , and thus again transverse isotropy is imposed to the acoustic tensor. Finally, the direction of propagation [ 011 ] is along to a rotation axis of symmetry of order 2, the physical point group is thus conjugate to D 2 Z 2 c . Here, the symmetry class of the acoustic tensor is D 2 Z 2 c , and all the eigenspaces are unidimensional. In this last case, as this kind of symmetry can be seen by second order tensors, the results from FEA on the heterogeneous material and classic elasticity are in agreement in the LF-LW regime.
In this section we have shown that some discrepancies can be observed when using an overall homogeneous continuum of Cauchy type. Classical elasticity (as opposed to generalized elasticity) is not rich enough to capture certain specific physical phenomena related to the symmetries of the material. In particular, if phase velocities are correctly estimated the polarizations are incorrectly predicted. Moreover, as it is well known, the onset of dispersion when frequency or wavenumber increase cannot be described in the classic Cauchy model.

6. Conclusions

In this work it was shown that a classical continuum model could not capture the correct behavior of elastic waves propagating in gyroid lattices. This is due to the fact that the classic continuum mechanics was insensitive to the lack of centrosymmetry of the architectured material. However, it is a well-established belief that the effects of non-centrosymmetry are only related to waves having a wavelength which compares to the size of the microstructure. Here we demonstrate that the solution given by the classical theory failed to predict the correct response, even in the Long Wavelength-Low Frequency domain.
In order to capture the onset of this unconventional behavior, called acoustical activity, the elastic continuum model needs to be enriched. Different strategies of enrichment are possible. In particular, the use of strain-gradient elasticity will be investigate in a forthcoming study.
The main practical consequence of the results presented in this work was in the evidence that circularly polarized waves could be observed in gyroid lattices, and that classic models failed to describe such effect. This could have a practical interest, since devices based on manipulation of circular polarization are frequently used in optics and electromagnetism. However, if one wants to exploit the same effects in mechanics, it appears important not to rely on classical theories of elasticity. Finally, it should be noted that, in this work, we addressed bulk wave propagation in an infinite medium. The interaction of these waves with boundaries, in the case of reflection/transmission problems or in the case of guided propagation, also deserves to be investigated.

Author Contributions

Conceptualization, G.R., N.A. and C.C.; Formal analysis, G.R., N.A. and C.C.; Methodology, G.R., N.A. and C.C.; Validation, G.R.; Software, G.R.; Writing—original draft, G.R., N.A. and C.C.; Writing—review & editing, G.R., N.A. and C.C. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-19-CE08-0005 (project MaxOasis). This work was partially funded by CNRS/IRP Coss&Vita between Fédération Francilienne de Mécanique (F2M, CNRS FR2609) and M&MoCS.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
BCCBody Centered Cubic
BZBrillouin Zone
IBZIrreducible Brillouin Zone
FEAFinite Elements Analysis
FCCFace Centered Cubic
LFLow Frequency
LWLong Wavelength

Appendix A. Dictionary

To obtain the normal forms for the different classes the generators provided in the following table have been used:
Table A1. The set of group generators used to construct matrix representation for each symmetry class.
Table A1. The set of group generators used to construct matrix representation for each symmetry class.
GroupGenerators
Z 2 P e 3
Z n R e 3 ; 2 π n
D n R e 3 ; 2 π n , R ( e 1 ; π )
Z 2 n , n 2 R e 3 ; π n
D 2 n h n 2 R e 3 ; π n , R ( e 1 , π )
D n v R e 3 ; 2 π n , P e 1
T R ( e 3 ; π ) , R ( e 1 ; π ) , R ( e 1 + e 2 + e 3 ; 2 π 3 )
O R ( e 3 ; π 2 ) , R ( e 1 ; π ) , R ( e 1 + e 2 + e 3 ; 2 π 3 )
O R ( e 3 ; π 2 ) , P e 2 e 3
Where the following elements of O ( 3 ) will be used in this study:
  • R ( v ; θ ) SO ( 3 ) the rotation about v R 3 through an angle θ [ 0 ; 2 π ] ;
  • P n O ( 3 ) SO ( 3 ) the reflection through the plane normal to n ( P n = 1 2 n n ).
Type I Subgroups
Table A2. Dictionary between different group notations for Type I subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.
Table A2. Dictionary between different group notations for Type I subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.
SystemHermann−Maugin Schonflies Group Nature
Triclinic 1 Z 1 1 I ¯ CP
Monoclinic 2 C 2 Z 2 I ¯ CP
Orthotropic 222 D 2 D 2 I ¯ C P ¯
Trigonal 3 C 3 Z 3 I ¯ CP
Trigonal 32 D 3 D 3 I ¯ C P ¯
Tetragonal 4 C 4 Z 4 I ¯ CP
Tetragonal 422 D 4 D 4 I ¯ C P ¯
Hexagonal 6 C 6 Z 6 I ¯ CP
Hexagonal 622 D 6 D 6 I ¯ C P ¯
C SO ( 2 ) I ¯ CP
2 D O ( 2 ) I ¯ C P ¯
Cubic 23 T T I ¯ C P ¯
Cubic 432 O O I ¯ C P ¯
532 I I I ¯ C P ¯
SO ( 3 ) I ¯ C P ¯
Type II Subgroups
Table A3. Dictionary between different group notations for Type II subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.
Table A3. Dictionary between different group notations for Type II subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.
SystemHermann−Maugin Schonflies Group Nature
Triclinic 1 ¯ C i Z 2 c I C ¯ P ¯
Monoclinic 2 / m C 2 h Z 2 Z 2 c I C ¯ P ¯
Orthotropic m m m D 2 h D 2 Z 2 c I C ¯ P ¯
Trigona 3 ¯ S 6 , Z 3 i Z 3 Z 2 c I C ¯ P ¯
Trigonal 3 ¯ m D 3 d D 3 Z 2 c I C ¯ P ¯
Tetragonal 4 / m C 4 h Z 4 Z 2 c I C ¯ P ¯
Tetragonal 4 / m m m D 4 h D 4 Z 2 c I C ¯ P ¯
Hexagonal 6 / m C 6 h Z 6 Z 2 c I C ¯ P ¯
Hexagonal 6 / m m m D 6 h D 6 Z 2 c I C ¯ P ¯
/ m C h SO ( 2 ) Z 2 c I C ¯ P ¯
/ m m D h O ( 2 ) Z 2 c I C ¯ P ¯
Cubic m 3 ¯ T h T Z 2 c I C ¯ P ¯
Cubic m 3 ¯ m O h O Z 2 c I C ¯ P ¯
5 ¯ 3 ¯ m I h I Z 2 c I C ¯ P ¯
/ m / m O ( 3 )
Type III Subgroups
Table A4. Dictionary between different group notations for Type III subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.
Table A4. Dictionary between different group notations for Type III subgroups. The last column indicates the nature of the group: C = Chiral, P = Polar, I = Centrosymetric, and overline indicates that the property is missing.
SystemHermann−Maugin Schonflies Group Nature
Monocinic m C s Z 2 I ¯ C ¯ P
Orthotropic 2 m m C 2 v D 2 v I ¯ C ¯ P
Trigonal 3 m C 3 v D 3 v I ¯ C ¯ P
Tetragonal 4 ¯ S 4 Z 4 I ¯ C ¯ P ¯
Tetragonal 4 m m C 4 v D 4 v I ¯ C ¯ P
Tetragonal 4 ¯ 2 m D 2 d D 4 h I ¯ C ¯ P ¯
Hexagonal 6 ¯ C 3 h Z 6 I ¯ C ¯ P ¯
Hexagonal 6 m m C 6 v D 6 v I ¯ C ¯ P
Hexagonal 6 ¯ 2 m D 3 h D 6 h I ¯ C ¯ P ¯
Cubic 4 ¯ 3 m T d O I ¯ C ¯ P ¯
m C v O ( 2 ) I ¯ C ¯ P

Appendix B. Generators of Space Group #214

Consider the affine space E 3 , the vector space R 3 acts on E 3 by translations. The affine group Aff ( E 3 ) of E 3 , which is the set of all affine invertible transformations is constructed as the semidirect product of R 3 by GL ( 3 ) , the general linear group of R 3 :
Aff ( E 3 ) = GL ( R 3 ) R 3
as such, an affine transformation is given by a pair ( Q , v ) GL ( R 3 ) × R 3 . Composition of transformations follows from the construction of Aff ( E 3 ) as a semi-direct product, to be explicit:
( Q 2 , v 2 ) ( Q 1 , v 1 ) = ( Q 2 Q 1 , Q 2 v 1 + v 2 ) .
Elements of Aff ( E 3 ) can be nicely represented by ( 4 × 4 ) block matrices:
Q v 0 1
the internal law in Aff ( E 3 ) following the matrix product in M 4 , 4 .
For our needs, we are interested not in the full affine group but in the group of isometries of E 3 , this group Euc ( E 3 ) is a subgroup of Aff ( E 3 ) and defined as the semi direct product of the orthogonal group and the spatial translation of R 3 :
Euc ( E 3 ) = O ( R 3 ) R 3 .
Space groups can be considered as discrete subgroups of Euc ( E 3 ) .
The generators of the space group I 4 1 32 (No. 214) are given in the following table in various notations [54]:
Table A5. Generators of Group I 4 1 32 (No. 214).
Table A5. Generators of Group I 4 1 32 (No. 214).
SeitzMathMatrices in Conventional Basis B
{ 2 001 | 1 2 0 1 2 } R ( π , e 3 ) ; 1 2 ( e 1 + e 3 ) 1 0 0 1 2 0 1 0 0 0 0 1 1 2 0 0 0 1 B
{ 2 010 | 0 1 2 1 2 } R ( π , e 2 ) ; 1 2 ( e 2 + e 3 ) 1 0 0 0 0 1 0 1 2 0 0 1 1 2 0 0 0 1 B
{ 3 111 + | 0 } R ( 2 π 3 ; e 1 + e 2 + e 3 ) , 0 ̲ 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 1 B
{ 2 110 | 3 4 1 4 1 4 } R ( π , e 1 + e 2 ) ; 1 4 ( 3 e 1 + e 2 + e 3 ) 0 1 0 3 4 1 0 0 1 4 0 0 1 1 4 0 0 0 1 B
{ 1 | 1 2 1 2 1 2 } Id ; 1 2 ( e 1 + e 2 + e 3 ) 1 0 0 1 2 0 1 0 1 2 0 0 1 1 2 0 0 0 1 B

Appendix C. Proof

The gyroid lattice is defined from an implicit equation (Equation (1)) that creates a periodic surface. For a given value of parameter b (= 2 ), this surface is found to become singular, thus creating an unrealistic discontinuous solid. This section presents an explanation for the admissible variation range of gyroid parameter: | b | < 2 .
Let us first restrict the variation range of variables x , y , z in Equation (2) to [ 0 , 1 2 ] in order to work in the fundamental domain of function ϕ . The gyroid lattice restricted to this domain is presented in Figure A1a).
Figure A1. (a) The gyroid restricted to its fundamental domain along with symmetry axes C 3 (plain) and C 2 (dashed) and (b) evolution of parameter b as a function of x position along e 1 axis, dashed line corresponds to b = 2 .
Figure A1. (a) The gyroid restricted to its fundamental domain along with symmetry axes C 3 (plain) and C 2 (dashed) and (b) evolution of parameter b as a function of x position along e 1 axis, dashed line corresponds to b = 2 .
Symmetry 12 01243 g0a1
The fundamental domain of the gyroid is invariant with respect to the following symmetry operations:
-
Rotation of angle 2 π 3 along the axis defined by equations y = z = x , plotted in plain line in Figure A1 and corresponding to the transformation ( x , y , z ) ( y , z , x ) . The directing vector of this axis is ( 1 , 1 , 1 ) in orthonormal basis B and passes through point ( 0 , 0 , 0 ) ;
-
Three rotations of angle π about the three axes defined by equations { y = 1 4 x , z = 1 8 } , { z = 1 4 x , y = 1 8 } , and { z = 1 4 y , x = 1 8 } and plotted in dashed lines in Figure A1. These axes correspond to transformations ( x , y , z ) ( 1 4 y , 1 4 x , 1 4 z ) , ( x , y , z ) ( 1 4 z , 1 4 y , 1 4 x ) , ( x , y , z ) ( 1 4 x , 1 4 z , 1 4 y ) , respectively. The directing vector of these axes are, in orthonormal basis B , ( 1 , 1 , 0 ) , ( 1 , 0 , 1 ) and ( 0 , 1 , 1 ) and they pass through points ( 0 , 1 4 , 1 8 ) , ( 0 , 1 8 , 1 4 ) and ( 1 8 , 0 , 1 4 ) , respectively.
It is trivial to see that these transformations leave Equation (2) unchanged thus defining symmetry operations. As a consequence, the point symmetry group of the fundamental domain of the gyroid is conjugated to D 3 . For the sake of simplicity, we will only consider generating operations of 2 π 3 rotation about ( x , x , x ) axis and π rotation about ( x , 1 / 4 x , 1 / 8 ) axis, denoted C 3 and C 2 in the following, respectively; the two other π rotations being generated by combination of these two generators.
If the gyroid surface intersects one of the rotation axes non-orthogonally, then the surface automatically becomes degenerate. The expression of the normal director to the gyroid surface at its intersection point with generating symmetry axes and the equation defining this intersection point are summarized in the following Table A6:
Table A6. The expression of the normal director to the gyroid surface at the intersection point with its symmetry axes and expression of this intersection point.
Table A6. The expression of the normal director to the gyroid surface at the intersection point with its symmetry axes and expression of this intersection point.
Sym. axis & Directing VectorNormal to the Gyroid SurfaceIntersection Point
C 3 ( 1 , 1 , 1 ) B cos 2 2 π x sin 2 2 π x ( 1 , 1 , 1 ) B ( x , x , x ) with 3 cos 2 π x sin 2 π x = b
C 2 ( 1 , 1 , 0 ) B sin 2 π x 2 2 cos 2 π x ( 1 , 1 , 0 ) B ( x , x , 0 ) with 2 cos 2 π x + sin 2 2 π x = b
Note that the normal director to the gyroid surface depends on variable x which, itself, is determined by parameter b through the non-linear equation defining the intersection point at which the normal director is computed.
One can easily check that the normal directors are generically colinear with directing vectors of C 3 and C 2 operations. However, for given values of variable x (or equivalently of parameter b), the normal to the gyroid surface is null, thus leading to singularity of the gyroid surface. These values are x = 0 –and thus b = 2 –leading to singularity of the gyroid surface at its intersection with axis C 2 and x = 1 8 –and thus b = 3 2 –leading to singularity of the gyroid surface at its intersection with both axes.
Finally, the equation defining the intersection point between the gyroid surface and the C 2 axis (see Figure A1b) depends on x and parameter b: 2 cos x + sin x 2 b = 0 . By plotting this equation considering b as a function of x, we can see that there are two intersection points between the gyroid surface and the C 2 axis for values of b over 2 thus showing that the gyroid surface forms a closed domain in these directions leading to an unrealistic discontinuous solid.

References

  1. Schaedler, T.A.; Carter, W.B. Architected Cellular Materials. Annu. Rev. Mater. Res. 2016, 46, 187–210. [Google Scholar] [CrossRef]
  2. Ashby, M.; Bréchet, Y. Designing hybrid materials. Acta Mater. 2003, 51, 5801–5821. [Google Scholar] [CrossRef]
  3. Fratzl, P.; Weinkamer, R. Nature’s hierarchical materials. Prog. Mater. Sci. 2007, 52, 1263–1334. [Google Scholar]
  4. Estrin, Y.; Bréchet, Y.; Dunlop, J.; Fratzl, P. Architectured Materials in Nature and Engineering; Springer: Berlin/Heidelberg, Germany, 2019. [Google Scholar]
  5. Hales, T.C. The Honeycomb Conjecture. Discret. Comput. Geom. 2001, 25, 1–22. [Google Scholar] [CrossRef] [Green Version]
  6. Dolan, J.A.; Wilts, B.D.; Vignolini, S.; Baumberg, J.J.; Steiner, U.; Wilkinson, T.D. Optical Properties of Gyroid Structured Materials: From Photonic Crystals to Metamaterials. Adv. Opt. Mater. 2015, 3, 12–32. [Google Scholar] [CrossRef] [Green Version]
  7. Wilts, B.D.; Michielsen, K.; Raedt, H.D.; Stavenga, D.G. Iridescence and spectral filtering of the gyroid-type photonic crystals in Parides sesostris wing scales. Interface Focus 2012, 2, 681–687. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Boutin, C.; Auriault, J. Rayleigh scattering in elastic composite materials. Int. J. Eng. Sci. 1993, 31, 1669–1689. [Google Scholar]
  9. Parnell, W.; Abrahams, I. Homogenization for wave propagation in periodic fibre-reinforced media with complex microstructure. i—theory. J. Mech. Phys. Solids 2008, 56, 2521–2540. [Google Scholar] [CrossRef]
  10. Nassar, H.; He, Q.C.; Auffray, N. Willis elastodynamic homogenization theory revisited for periodic media. J. Mech. Phys. Solids 2015, 77, 158–178. [Google Scholar] [CrossRef] [Green Version]
  11. Hermann, C. Tensoren und Kristallsymmetrie. Zeitschrift Kristallogr. 1934, 32–48. [Google Scholar] [CrossRef]
  12. Olive, M.; Auffray, N. Symmetry classes for even-order tensors. Math. Mech. Complex Syst. 2013, 1, 177–210. [Google Scholar] [CrossRef]
  13. Olive, M.; Auffray, N. Symmetry classes for odd-order tensors. ZAMM J. Appl. Math. Mech. Z. Für Angew. Math. Und Mech. 2014, 94, 421–447. [Google Scholar] [CrossRef] [Green Version]
  14. DiVincenzo, D.P. Dispersive corrections to continuum elastic theory in cubic crystals. Phys. Rev. B 1986, 34, 5450–5465. [Google Scholar] [CrossRef] [PubMed]
  15. Auffray, N.; Dirrenberger, J.; Rosi, G. A complete description of bi-dimensional anisotropic strain-gradient elasticity. Int. J. Solids Struct. 2015, 69–70, 195–206. [Google Scholar] [CrossRef]
  16. Rosi, G.; Auffray, N. Anisotropic and dispersive wave propagation within strain-gradient framework. Wave Motion 2016, 63, 120–134. [Google Scholar] [CrossRef] [Green Version]
  17. Eremeyev, V.A. On the material symmetry group for micromorphic media with applications to granular materials. Mech. Res. Commun. 2018, 94, 8–12. [Google Scholar] [CrossRef]
  18. Eremeyev, V.A.; Rosi, G.; Naili, S. Comparison of anti-plane surface waves in strain-gradient materials and materials with surface stresses. Math. Mech. Solids 2018, 24, 2526–2535. [Google Scholar] [CrossRef]
  19. Eremeyev, V.A.; Rosi, G.; Naili, S. Transverse surface waves on a cylindrical surface with coating. Int. J. Eng. Sci. 2019, 103188. [Google Scholar] [CrossRef]
  20. Arago, F. Mémoire sur une Modification Remarquable Qu’éProuvent Les Rayons Lumineux Dans Leur Passage à Travers Certains Corps Diaphanes et sur Quelques Autres Nouveaux Phénomènes D’Optique; Institut National de France: Paris, France, 1811; Volume 1, pp. 93–134. [Google Scholar]
  21. Portigal, D.L.; Burstein, E. Acoustical Activity and Other First-Order Spatial Dispersion Effects in Crystals. Phys. Rev. 1968, 170, 673–678. [Google Scholar] [CrossRef]
  22. Sivardière, J. Description de la Symétrie; EDP: Ulysse, France, 2004. [Google Scholar]
  23. Auffray, N.; He, Q.C.; Le Quang, H. Complete symmetry classification and compact matrix representations for 3D strain gradient elasticity. Int. J. Solids Struct. 2019, 159, 197–210. [Google Scholar] [CrossRef] [Green Version]
  24. Prall, D.; Lakes, R.S. Properties of a chiral honeycomb with a Poisson’s ratio of—1. Int. J. Mech. Sci. 1997, 39, 305–307, 309–314. [Google Scholar] [CrossRef]
  25. Lakes, R. Elastic and viscoelastic behavior of chiral materials. Int. J. Mech. Sci. 2001, 43, 1579–1589. [Google Scholar] [CrossRef] [Green Version]
  26. Spadoni, A.; Ruzzene, M.; Gonella, S.; Scarpa, F. Phononic properties of hexagonal chiral lattices. Wave Motion 2009, 46, 435–450. [Google Scholar] [CrossRef] [Green Version]
  27. Liu, X.N.; Hu, G.K.; Sun, C.T.; Huang, G.L. Wave propagation characterization and design of two-dimensional elastic chiral metacomposite. J. Sound Vib. 2011, 330, 2536–2553. [Google Scholar] [CrossRef]
  28. Liu, X.N.; Huang, G.L.; Hu, G.K. Chiral effect in plane isotropic micropolar elasticity and its application to chiral lattices. J. Mech. Phys. Solids 2012, 60, 1907–1921. [Google Scholar] [CrossRef] [Green Version]
  29. Dirrenberger, J.; Forest, S.; Jeulin, D. Effective elastic properties of auxetic microstructures: Anisotropy and structural applications. Int. J. Mech. Mater. Des. 2013, 9, 21–33. [Google Scholar] [CrossRef]
  30. Bacigalupo, A.; Gambarotta, L. Homogenization of periodic hexa-and tetrachiral cellular solids. Compos. Struct. 2014, 116, 461–476. [Google Scholar] [CrossRef] [Green Version]
  31. Fernandez-Corbaton, I.; Rockstuhl, C.; Ziemke, P.; Gumbsch, P.; Albiez, A.; Schwaiger, R.; Frenzel, T.; Kadic, M.; Wegener, M. New Twists of 3D Chiral Metamaterials. Adv. Mater. (Deerfield Beach, Fla.) 2019, 31, e1807742. [Google Scholar] [CrossRef] [Green Version]
  32. Chen, Y.; Frenzel, T.; Guenneau, S.; Kadic, M.; Wegener, M. Mapping acoustical activity in 3D chiral mechanical metamaterials onto micropolar continuum elasticity. J. Mech. Phys. Solids 2020, 103877. [Google Scholar] [CrossRef]
  33. Ziemke, P.; Frenzel, T.; Wegener, M.; Gumbsch, P. Tailoring the characteristic length scale of 3D chiral mechanical metamaterials. Extrem. Mech. Lett. 2019, 32, 100553. [Google Scholar] [CrossRef]
  34. Chen, Y.; Yao, H.; Wang, L. Acoustic band gaps of three-dimensional periodic polymer cellular solids with cubic symmetry. J. Appl. Phys. 2013, 114, 043521. [Google Scholar] [CrossRef]
  35. Rammohan, A.V.; Lee, T.; Tan, V.B.C. A Novel Morphological Model of Trabecular Bone Based on the Gyroid. Int. J. Appl. Mech. 2015, 07, 1550048. [Google Scholar] [CrossRef]
  36. Ma, S.; Tang, Q.; Feng, Q.; Song, J.; Han, X.; Guo, F. Mechanical behaviours and mass transport properties of bone-mimicking scaffolds consisted of gyroid structures manufactured using selective laser melting. J. Mech. Behav. Biomed. Mater. 2019, 93, 158–169. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Poncelet, M.; Somera, A.; Morel, C.; Jailin, C.; Auffray, N. An experimental evidence of the failure of Cauchy elasticity for the overall modeling of a non-centro-symmetric lattice under static loading. Int. J. Solids Struct. 2018, 147, 223–237. [Google Scholar] [CrossRef] [Green Version]
  38. Abdoul-Anziz, H.; Seppecher, P.; Bellis, C. Homogenization of frame lattices leading to second gradient models coupling classical strain and strain-gradient terms. Math. Mech. Solids 2019, 24, 3976–3999. [Google Scholar] [CrossRef] [Green Version]
  39. Yvonnet, J.; Auffray, N.; Monchiet, V. Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior. Int. J. Solids Struct. 2020, 191, 434–448. [Google Scholar] [CrossRef]
  40. Schoen, A.H. Infinite Periodic Minimal Surfaces without Self-Intersections; Nasa Technical Notes TN D-5541; NTRS: Washington, DC, USA, 1970; p. 100. [Google Scholar]
  41. Schoen, A.H. Reflections concerning triply-periodic minimal surfaces. Interface Focus 2012, 2, 658–668. [Google Scholar] [CrossRef] [Green Version]
  42. Dacorogna, B. Introduction to the Calculus of Variations; World Scientific Publishing Company: Singapore, 2014. [Google Scholar]
  43. Wohlgemuth, M.; Yufa, N.; Hoffman, J.; Thomas, E.L. Triply periodic bicontinuous cubic microdomain morphologies by symmetries. Macromolecules 2001, 34, 6083–6089. [Google Scholar] [CrossRef]
  44. Golubitsky, M.; Stewart, I.; Schaeffer, D.G. Singularities and Groups in Bifurcation Theory; Applied Mathematical Sciences; Springer: New York, NY, USA, 1988; Volume 69. [Google Scholar] [CrossRef]
  45. Hahn, T.; Shmueli, U.; Wilson, A.J.C.; Prince, E. International Tables for Crystallography; D. Reidel Publishing Company: Dordrecht, The Netherlands, 2005. [Google Scholar]
  46. Craster, R.V.; Antonakakis, T.; Makwana, M.; Guenneau, S. Dangers of using the edges of the Brillouin zone. Phys. Rev. B 2012, 86, 115130. [Google Scholar] [CrossRef] [Green Version]
  47. Gazalet, J.; Dupont, S.; Kastelik, J.C.; Rolland, Q.; Djafari-Rouhani, B. A tutorial survey on waves propagating in periodic media: Electronic, photonic and phononic crystals. Perception of the Bloch theorem in both real and Fourier domains. Wave Motion 2013, 50, 619–654. [Google Scholar] [CrossRef]
  48. Achenbach, J. Wave Propagation in Elastic Solids; Elsevier: Amsterdam, The Netherlands, 1984. [Google Scholar]
  49. Zakharenko, A.A. On cubic crystal anisotropy for waves with Rayleigh-wave polarization. Nondestruct. Test. Eval. 2006, 21, 61–77. [Google Scholar] [CrossRef]
  50. Rosi, G.; Nguyen, V.H.; Naili, S. Surface waves at the interface between an inviscid fluid and a dipolar gradient solid. Wave Motion 2015, 53, 51–65. [Google Scholar] [CrossRef]
  51. Gourgiotis, P.A.; Georgiadis, H.G.; Neocleous, I. On the reflection of waves in half-spaces of microstructured materials governed by dipolar gradient elasticity. Wave Motion 2013, 50, 437–455. [Google Scholar] [CrossRef]
  52. Mandel, J. Généralisation de la théorie de plasticité de WT Koiter. Int. J. Solids Struct. 1965, 1, 273–295. [Google Scholar] [CrossRef]
  53. Forte, S.; Vianello, M. Symmetry classes for elasticity tensors. J. Elast. 1996, 43, 81–108. [Google Scholar] [CrossRef]
  54. Aroyo, M.I.; Perez-Mato, J.M.; Capillas, C.; Kroumova, E.; Ivantchev, S.; Madariaga, G.; Kirov, A.; Wondratschek, H. Bilbao Crystallographic Server: I. Databases and crystallographic computing programs. Z. Für Krist. Cryst. Mater. 2006, 221, 15–27. [Google Scholar]
Figure 1. The unit cells of gyroid lattices obtained for a = 1 mm and: (a) b = 0 , (b) b = 1 , (c) b = 1.3 . Despite what the angle of view may suggest, all these structures are simply connected.
Figure 1. The unit cells of gyroid lattices obtained for a = 1 mm and: (a) b = 0 , (b) b = 1 , (c) b = 1.3 . Despite what the angle of view may suggest, all these structures are simply connected.
Symmetry 12 01243 g001
Figure 2. The relationship between the parameter b and the porosity, dashed line represents a linear approximation of real porosity plotted in plain line. The evaluation of the porosity is obtained numerically.
Figure 2. The relationship between the parameter b and the porosity, dashed line represents a linear approximation of real porosity plotted in plain line. The evaluation of the porosity is obtained numerically.
Symmetry 12 01243 g002
Figure 3. The conventional unit cell. The lattice vectors are indicated in blue.
Figure 3. The conventional unit cell. The lattice vectors are indicated in blue.
Symmetry 12 01243 g003
Figure 4. Two examples of Primitive Unit Cells (PUC). In (a,b) the Conventional Unit Cell are in blue and the PUC are in red. The PUC in (a) is considered in the present paper, as their lattice vectors (indicated in red) are more symmetrical than those of the one in (b).
Figure 4. Two examples of Primitive Unit Cells (PUC). In (a,b) the Conventional Unit Cell are in blue and the PUC are in red. The PUC in (a) is considered in the present paper, as their lattice vectors (indicated in red) are more symmetrical than those of the one in (b).
Symmetry 12 01243 g004
Figure 5. Some examples of polarizations: (a) Linear, (b) circular right handed, and (c) circular left handed.
Figure 5. Some examples of polarizations: (a) Linear, (b) circular right handed, and (c) circular left handed.
Symmetry 12 01243 g005
Figure 6. The meshed unit cell used in simulations.
Figure 6. The meshed unit cell used in simulations.
Symmetry 12 01243 g006
Figure 7. The dispersion diagram of the Gyroid lattice computed along the boundaries of the Irreducible Brillouin Zone (IBZ) and directions of propagation with respect to the unit cell. (a) The dispersion relation of the Gyroid lattice computed along the boundaries of the IBZ. (b) The direction of propagation with respect to the unit cell.
Figure 7. The dispersion diagram of the Gyroid lattice computed along the boundaries of the Irreducible Brillouin Zone (IBZ) and directions of propagation with respect to the unit cell. (a) The dispersion relation of the Gyroid lattice computed along the boundaries of the IBZ. (b) The direction of propagation with respect to the unit cell.
Symmetry 12 01243 g007
Figure 8. The photonic band diagram of a single gyroid photonic crystal. Reproduced with permission from [6].
Figure 8. The photonic band diagram of a single gyroid photonic crystal. Reproduced with permission from [6].
Symmetry 12 01243 g008
Figure 9. The phase velocity of the circularly polarized waves in function of the wavenumber or of the ratio between the wavelength λ and the size of the unit cell a for propagation direction [ 001 ] (a) and [ 111 ] (b). Right handed waves are in black, left handed waves are in gray. The dashed horizontal lines correspond to the phase velocity to which they converge to for infinite wavelength.
Figure 9. The phase velocity of the circularly polarized waves in function of the wavenumber or of the ratio between the wavelength λ and the size of the unit cell a for propagation direction [ 001 ] (a) and [ 111 ] (b). Right handed waves are in black, left handed waves are in gray. The dashed horizontal lines correspond to the phase velocity to which they converge to for infinite wavelength.
Symmetry 12 01243 g009
Figure 10. An illustration of circular birefringence observed along direction [100] at f c = 200 kHz, corresponding to a wavelength around 10 times the size of the unit cell. A linearly polarized wave entering the material is subjected to a rotation of 1.3 degrees/wavelength, then after 20 wavelengths the rotation is θ t = 26.0 deg. This illustration does not account for the changes in amplitude due to the reflections at boundaries.
Figure 10. An illustration of circular birefringence observed along direction [100] at f c = 200 kHz, corresponding to a wavelength around 10 times the size of the unit cell. A linearly polarized wave entering the material is subjected to a rotation of 1.3 degrees/wavelength, then after 20 wavelengths the rotation is θ t = 26.0 deg. This illustration does not account for the changes in amplitude due to the reflections at boundaries.
Symmetry 12 01243 g010
Table 1. The high symmetry points of the gyroid lattice. The group notations are detailed in Appendix A.
Table 1. The high symmetry points of the gyroid lattice. The group notations are detailed in Appendix A.
SymmetryCoordinatesCoordinatesPoint GroupPoint GroupIllustration of the
Pointw.r.t. P w.r.t. B (Math.)(H-M)First Brillouin Zone
( k 1 , k 2 , k 3 ) ( x 1 , x 2 , x 3 )
Γ ( 0 , 0 , 0 ) ( 0 , 0 , 0 ) O 432 Symmetry 12 01243 i001
H ( 1 2 , 1 2 , 1 2 ) ( 0 , 0 , 1 a ) O 432
P ( 1 4 , 1 4 , 1 4 ) ( 1 2 a , 1 2 a , 1 2 a ) D 3 32
N ( 0 , 1 2 , 0 ) ( 0 , 1 2 a , 1 2 a ) D 2 222
Table 2. The polarizations of plane waves and conditions on the complex amplitude.
Table 2. The polarizations of plane waves and conditions on the complex amplitude.
PolarizationCondition
Longitudinal polarization U = α n with α C
Transverse polarization U R and U C belong to the plane orthogonal to n
Linear polarization U U * = 0 , with U * complex conjugate of U
Circular polarization U · U = 0
      ↳ Right handedness n · ( U R U C ) < 0
      ↳ Left handedness n · ( U R U C ) > 0
Elliptic polarization U U * 0 and U · U 0
Table 3. The parameters used in the numerical simulations, corresponding to bulk titanium.
Table 3. The parameters used in the numerical simulations, corresponding to bulk titanium.
Mass Density [kg/m 2 ]Young Modulus [GPa]Poisson Ratio [1]Porosity [1]Unit Cell Size [mm]
ρ b E b ν b p a
4506115.70.3210.71
Table 4. The phase velocity and polarization in the very long wavelength approximation.
Table 4. The phase velocity and polarization in the very long wavelength approximation.
DirectionPhase VelocityPolarizationType of Wave
[m/s] ( x 1 , x 2 , x 3 )
3018.5 ( 0 , 0 , 1 ) Longitudinal
[ 001 ] 1933.4 ( 1 , i , 0 ) Circular L ⥁
1931.0 ( 1 , i , 0 ) Circular R ⥀
3249.8 ( 0 , 1 , 1 ) Longitudinal
[ 011 ] 1931.7 ( 1 , 0 , 0 ) Transverse
1510.9 ( 0 , 1 , 1 ) Transverse
3322.8 ( 1 , 1 , 1 ) Longitudinal
[ 111 ] 1664.6 ( 1 , 0.54 i 0.88 , 0.46 + i 0.88 ) Circular R ⥀
1662.3 ( 1 , 0.45 + i 0.85 , 0.55 i 0.85 ) Circular L ⥁
Table 5. The phase velocity and polarization in classic elasticity.
Table 5. The phase velocity and polarization in classic elasticity.
DirectionPhase Velocities [m/s]PolarizationType of Wave
c 11 ρ 3018.5 ( 0 , 0 , 1 ) Longitudinal
[ 001 ] c 44 2 ρ 1932.2 ( 0 , 1 , 0 ) Transverse
c 44 2 ρ 1932.2 ( 1 , 0 , 0 ) Transverse
c 11 + c 12 + c 44 2 ρ 3249.8 ( 0 , 1 , 1 ) Longitudinal
[ 011 ] c 44 2 ρ 1932.2 ( 1 , 0 , 0 ) Transverse
c 11 c 12 2 ρ 1510.9 ( 0 , 1 , 1 ) Transverse
c 11 + 2 c 12 + 2 c 44 3 ρ 3322.8 ( 1 , 1 , 1 ) Longitudinal
[ 111 ] 2 c 11 2 c 12 + c 44 6 ρ 1663.4 ( 1 , 1 , 0 ) Transverse
2 c 11 2 c 12 + c 44 6 ρ 1663.4 ( 1 , 1 , 2 ) Transverse
Table 6. The material properties in the very long wavelength approximation.
Table 6. The material properties in the very long wavelength approximation.
Elastic Coefficients Mass Density
c 11 [GPa] c 12 [GPa] c 44 = 2 c 2323 [GPa] ρ [kg/m 3 ]
12.326.14510.0921351.8

Share and Cite

MDPI and ACS Style

Rosi, G.; Auffray, N.; Combescure, C. On the Failure of Classic Elasticity in Predicting Elastic Wave Propagation in Gyroid Lattices for Very Long Wavelengths. Symmetry 2020, 12, 1243. https://doi.org/10.3390/sym12081243

AMA Style

Rosi G, Auffray N, Combescure C. On the Failure of Classic Elasticity in Predicting Elastic Wave Propagation in Gyroid Lattices for Very Long Wavelengths. Symmetry. 2020; 12(8):1243. https://doi.org/10.3390/sym12081243

Chicago/Turabian Style

Rosi, Giuseppe, Nicolas Auffray, and Christelle Combescure. 2020. "On the Failure of Classic Elasticity in Predicting Elastic Wave Propagation in Gyroid Lattices for Very Long Wavelengths" Symmetry 12, no. 8: 1243. https://doi.org/10.3390/sym12081243

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop