Proceedings of the Institute of Acoustics

 

 

Towards bridging nanoscale and macroscale acoustics of porous solids

 

Alan Sam¹, Université Grenoble Alpes – Laboratoire Interdisciplinaire de Physique (LIPhy), Grenoble, France

Benoit Coasne², Université Grenoble Alpes – Laboratoire Interdisciplinaire de Physique (LIPhy), Grenoble, France

Rodolfo Venegas³, Universidad Austral de Chile (UACh) – Instituto de Acústica, Valdivia, Chile

 

ABSTRACT

 

This work presents results towards the bridging of nanoscale and macroscale acoustics of porous solids. We use molecular dynamics (MD) simulations to estimate the mechanical parameters as well as the acoustic properties of a nanoporous solid. These simulations consist in applying tensile and shear deformation to a prototypical nanoporous material (i.e. pure silica zeolite) and calculating the mechanical properties from the linear regime of the stress-strain curve. The specific outcomes of the simulations are phase velocities and mechanical parameters that can be used as inputs in meso/macroscopic models. To further exemplify the proposed nano-macro modelling strategy, wave propagation in a tube with an array of zeolite-made plates is studied, thereby providing a first approach to the acoustic modelling of multiscale metamaterials.

 

1. INTRODUCTION

 

Over the past few decades, the advancements in nanoscience and nanotechnologies have fostered overwhelming research on nanoporous materials – especially to decipher mechanisms to adsorb and transport fluid with a plethora of novel phenomena specific to the nanometer scale [1–3]. From a practical viewpoint, such research efforts are motivated by the ultralow density, high surface area and high strength-to-weight ratio of these nanoporous materials which make them suitable candidates for catalytic, gas-sensing, filtration and separation applications [4–6]. Effective utilization of such nanoscale confinement for various applications, however, is largely influenced by their mechanical performance. Despite numerous unique properties exhibited by nanoporous materials, their low-dense structure may cause mechanical deformation that could become highly significant when guest molecules are adsorbed within their porosity at the solid/fluid interface. In recent years, experimental and theoretical studies have shown significant progress toward explaining the mechanics of deformation in nanoporous media [7–10]. On the other hand, despite the intrinsic link between acoustic and mechanical properties, understanding acoustic wave propagation in nanoporous materials remains to be fully explored. Considering the intrinsic molecular, i.e. nanoscopic lengthscale, inherent to physical phenomena in these nanostructures, many questions remain unanswered regarding the acoustic properties of this class of materials. In this context, a fundamental understanding of acoustic wave propagation in these nanoporous media is highly interesting with many potential applications in the field of acoustics and nanofluidics.

 

A major challenge in conducting experiments for determining mechanical and acoustic properties of nanoporous materials is the difficulty in synthesizing predetermined architectures with well controlled microstructures. Moreover, when compared to mechanical tests for bulk i.e. dense materials, conducting precise measurements with solids involving an intrinsic nanoscopic lengthscale is rather complex. In this regard, molecular dynamics (MD) simulations have proved to be a powerful and efficient tool for understanding the physical mechanisms at the nanoscale. In the present study, we demonstrate how the mechanical properties of a nanoporous material with a prescribed atomic structure can be calculated by using molecular dynamics simulation. Then, we extend this study by showing how all information obtained at this vanishing lengthscale (i.e. molecular scale) can be upscaled to predict the acoustic behavior of mesoscopic architectures built from nanoporous material.

 

As a prototypical nanoporous material, we use a pure zeolite solid which consists of a crystalline phase of silica (SiO₂) having a very small porosity at the nanometer scale (pore size of the order of a few Å to 1-2 nm) [4, 11]. In more detail, we used a crystalline cubic zeolite material, i.e. RHO zeolite, with a 3-dimensional tetrahedral framework in which each oxygen atom is shared by two tetrahedra consisting of silicon (Si) atoms connected to 4 oxygen atoms. RHO zeolite possesses very small apertures with pore sizes of 3-4 Å enabling them to selectively adsorb molecules and exclude the large-sized molecules based on size exclusion and shape-selective basis. Besides its suitability for many applications involving adsorption, the cubic structure of RHO zeolite considerably reduces the complexity of the present study – especially from an acoustic point of view. For a material like RHO zeolite with a cubic structure, the total number of independent elastic constants reduce from 21 to just 3. Here, by applying uniaxial tensile and shear strain, the three independent elastic stiffnesses of RHO zeolite are predicted from molecular simulation. With the knowledge of mechanical constants obtained from MD, we calculated the longitudinal and transverse speeds of sound in different directions of acoustic wave propagation by applying Kelvin-Christoffel’s equation. As illustrated in the last part of this paper, the outcomes of the molecular simulation study can be used as inputs for an upscaled model of wave propagation in a metamaterial, thereby exemplifying our nano-macro modelling strategy.

 

2. METHODS

 

A series of full atomistic calculations of mechanical test cases are implemented via classical molecular dynamics (MD). This allows us to derive a simplified set of parameters to describe the RHO zeolite behavior from mechanical and, hence, acoustic viewpoints. In this section, we introduce the initial structural configuration of the supercell used in our molecular simulation strategy. Then, we present the details of the mechanical tests performed in MD to deform the zeolite supercell as well as the force-field used to model the bonded and non-bonded interactions of the silicon and oxygen atoms in the zeolite. It is stressed that, due to the preliminary nature of this study, no fluid saturating the material is considered in the MD simulations.

 

2.1. Structure

 

The positions of O and Si atoms in RHO zeolite were generated from the database of International Zeolite Association [12] (shown in Figure 1 a)). The experimental lattice parameter for this cubic zeolite is a = b = c = 14.9190 Å. To reduce finite-size effects while maintaining the computational time reasonable, the zeolite unit cell has tripled in each space direction. The resulting supercell is shown in Figure 1 b. A physical cut-off of 1.75 Å is used to define the Si-O bonds in the supercell. In Figure 1 c), a schematic of SiO₄ tetrahedron with Si-O bond length and O-Si-O angle are represented. In all MD simulations, periodic boundary conditions were implemented along all 3 directions in space.

 

 

Figure 1: a) Unit cell for RHO zeolite having a lattice parameter a = 14.9190 Å. The orange and blue atoms show the silicon and oxygen atoms, respectively. b) Super cell built by multiplying unit cell 3 times in x, y and z directions. The length of the super cell in each direction is 44.757 Å. c) Schematic of the SiO₄ tetrahedron with Si-O bond length of 1.75 Å and O-Si-O bond angle of 109.47^◦.

 

2.2. Molecular dynamics

 

The zeolite supercell was first equilibrated at constant external stress (1 atm in each direction) and temperature (300 K). The temperature was maintained constant using a Nosé-Hoover thermostat [13] with a relaxation time of 0.005 ps. The stress was imposed in each direction of space (uncoupled in each direction) using a Parrinello and Rahman barostat [14] with a relaxation time of 0.05 ps. After equilibration in this Nσ T, different mechanical tests on RHO zeolite were performed by keeping the system volume constant. In the tensile test, we deform the supercell by stretching along the x-direction at a uniform rate. To apply pure shear, the bottom of the zeolite was fixed along the x-axis, while the top surface was sheared with a constant stretching along the x-direction at a uniform rate. In both cases, the strain rate applied was 1.0 × 10⁻⁶ s⁻¹. We conducted test simulations on five independent structural configurations to estimate statistical errors associated with the results obtained from the mechanical tests. All molecular simulations carried out for the present study were implemented using the LAMMPS package [15] with a time step of 0.0005 ps.

 

In MD, the macroscopic stress developed while inducing tensile or shear strain to a solid material is calculated by taking into account the energy contribution from particle motion as well as virial stress.

 

The components of the macroscopic stress tensor (σ_ij) in a volume V, therefore, can be evaluated using Eq. 1.

 

 

where  and are the i th and j th component velocities of particle a having mass m. The virial contribution (W_ij) due to intra and intermolecular interactions can be expressed as in Eq.2;

 

 

where is the total force exerted by particle b on a along the jth vector component when the particles a and b are separated by a distance along the ith vector component.

 

2.3. Force-field

 

We employed the force-field proposed by Ghysels et al. [16] to describe the atomic interactions in zeolites with reasonably good accuracy. With this force-field, all Si-O, O-O and Si-Si interactions are described using Lennard-Jones (LJ) interaction potentials. These Lennard-Jones potentials employed to model such interactions were truncated at a distance of 13 Å . The angle between the bonds, Si-O-Si and O-Si-O, were described through harmonic functions while the contribution of torsion is neglected. The long-range electrostatic interactions in the supercell were calculated using the particle-particle particle-mesh (PPPM) solver [17] with an accuracy of 10⁻⁵. The partial charges on Si and O atoms in the zeolite are +2.1 e and -1.05 e, respectively. We used a cut-off distance of 13 Å for computing the electrostatic interaction. The parameters used to model the harmonic bonded and non-bonded LJ interactions are given in Table 1.

 

Table 1: The parameters used to model nonbonded (LJ) and bonded (harmonic) interactions in zeolite RHO super cell are given.

 

 

3. RESULTS AND DISCUSSION

 

3.1. Mechanical properties

 

The coefficients of the elastic stiffness fourth-order tensor C are fundamental parameters that describe the mechanical properties of materials. Hooke’s law is a linear response model in which the stress (σ) is assumed to be proportional to the strain (ϵ): σ_ij = C_ijkl ϵ_kl where Einstein summation over repeated indices is implicitly used. For a 3D material, C has 81 components. Applying the symmetry of stress (σ_ij = σ_ji) and strain (ϵ_kl = ϵ_lk) tensors lead to the symmetry in stiffness tensor (C_ijkl = C_klij). As a result, the number of independent elastic constants reduces from 81 to 21 for any linear elastic anisotropic material. Using Voigt’s notation, the elastic stiffness matrix of an anisotropic material can be represented as C_{i jkl} → C_mn, with m = i if i = j and m = 9 − (i + j) if i , j, and n = k if k = l and n = 9 − (k + jl) if k , l. Additional symmetry constraints apply depending on the crystal structure of the material, thereby further reducing the number of independent elastic constants. For a cubic material, like RHO zeolite considered in our study, only 3 elastic constants, C₁₁, C₁₂ and C₄₄, are required to describe its mechanical behavior. The stiffness matrix corresponding to a material having cubic symmetry is represented as:

 

 

MD simulations are used to predict the mechanical constants C₁₁, C₁₂ and C₄₄ for RHO zeolite through uniaxial tensile and shear tests. A schematic of these mechanical tests is shown in Figure 2. A constant strain rate of 1.0 × 10⁻⁶s⁻¹is applied to the material for inducing tensile and shear deformation. In tensile test (Figure 2 a)), RHO zeolite is deformed along the x-direction with the tensile strain applied to the supercell defined as ϵ_x = ∆L/L₀. For inducing shear deformation to the supercell, the top surface of the zeolite is displaced along the x-direction while keeping the bottom surface fixed in position (Figure 2 b). The total shear deformation applied to the zeolite structure can be evaluated in terms of the engineering shear strain γ_xy = ∆L / L₀.

 

In Fig. 3, the components of the stress tensor calculated using Eq.1 are plotted as a function of the strain applied to deform the zeolite. A linear fit to the stress-strain (σ vs ϵ) curve in the elastic regime yields the 3 independent elastic constants for RHO zeolite. To confirm the fitting of σ vs ϵ curve is in the linear elastic regime, the strain limit is chosen at approximately 10% of the strain applied for the mechanical failure of the material. The material constants C₁₁ and C₁₂ are predicted from the tensile test by taking the ratio of the stress components σ_x and σ_y to the tensile strain applied along x-direction (ϵ_x). The values of C₁₁ and C₁₂ are obtained as 106.81 and 71.96 GPa, respectively. For predicting C₄₄, i.e. the shear modulus G of a material, the stress developed (τ_xy) upon shearing the upper surface of the supercell in x-direction is calculated while keeping the lower surface of zeolite fixed in position. The value of C₄₄ = G for RHO zeolite obtained from the shear test using MD is 38.86 GPa. To verify the validity of our results, considering the cubic symmetry of the material under study, it is checked that tensile tests in y and z-directions and shear tests on xy and zx planes yield the same values. From the tensile tests we found C₁₁ = C₂₂ = C₃₃ and C₁₂ = C₁₃ = C₂₃ and the elastic constants C₅₅ and C₆₆ calculated from the shear tests to be the same as C₄₄. Further, the elastic stiffness matrix of RHO zeolite is found to obey the necessary and sufficient elastic stability conditions (C₁₁ - C₁₂ > 0; C₁₁ + 2C₁₂ > 0; C₄₄ > 0) for the cubic system [18].

 

Using the elastic parameters C_ij and noting that the Zener ratio A = 2C₄₄ / (C₁₁−C₁₂) = 2.23, which reveals weak anisotropy; we determined a single, representative bulk modulus K, shear modulus G, Young’s modulus E and Poisson’s ratio ν for RHO zeolite as follows [19].

 

 

Figure 2: Schematic of the uniaxial tensile (a) and shear (b) tests performed on zeolite RHO super cell using MD simulations. The undeformed configuration of the zeolite is represented as opaque while the deformed configuration after tensile and shear tests is shown transparent. The orange and blue spheres correspond to the O and Si atoms in zeolite, respectively. L₀ is the length of the undeformed system. a) For the tensile test, the 3D zeolite structure is deformed along x-direction. The tensile strain ϵ_x is defined as the ratio of total change in length due to tensile deformation ∆L to the initial length of the super cell (L₀). b) To apply shear, the top of the zeolite (yz plane) is deformed along x-direction while keeping the bottom fixed in position. The shear strain γ_xy is defined as the ratio between the amount of deformation and the length perpendicular to the surface in which the deformation is applied.

 

 

Figure 3: Stress-strain curves (σ_{i j} versus ϵ_{i j}) for zeolite RHO obtained from the mechanical tests using MD simulations. The elastic constants C₁₁ and C₁₂ are calculated from the uniaxial tensile test while C₄₄ is determined by performing shear test on the zeolite structure. The red dashed line shows a linear fit to stress-strain curve in the small strain regime to yield the respective elastic constant. C₁₁, C₁₂ and C₄₄ are represented in a), b) and c) respectively.

 

The elastic constants and mechanical properties for RHO zeolite are shown in Table 2. The bulk modulus K calculated from our study is found to be in good agreement with the value obtained from a previous DFT study [20].

 

Table 2: Elastic stiffness and mechanical parameters of RHO zeolite as obtained using MD simulations. All units are in GPa except for Poisson’s ratio which is dimensionless. The bulk modulus K from DFT (ab initio) calculations [20] is shown in parentheses.

 

 

3.2. Acoustic properties

 

The longitudinal and transverse modes of sound propagation in a material with respect to any symmetry axes is correlated to the stiffness tensor through Kelvin-Christoffel’s equation.

 

 

where p_k represents the unit displacement vector (or polarization vector), v the phase velocity and ρ the mass density of the material. As indicated by Equation 5, the determination of wave motion is an eigenvalue problem where the Christoffel’s tensor Γ_ik = C_{i jkl}n_jn_l is related to the stiffness tensor C and the direction cosines vector n. Due to the stiffness tensor symmetry, Γ_ik − ρv²δ_ik is also symmetric so that the resulting three eigenvalues on solving Equation 5 are real. We estimated the longitudinal [v_L] and transverse [v_T1 and v_T2] speeds of sound along axes [100], [110] and [111] for RHO zeolite by solving Equation 5. The results of the speeds of sound averaged over 5 simulations with different initial configurations and the standard error of the average are shown in Table 3.

 

Table 3: Longitudinal [v_L] and transverse [v_T1 and v_T2] speeds of sound in RHO zeolite for different directions of wave propagation. The formulas for propagation velocities along [100], [110] and [111] axes of symmetry derived by solving Eq.5 are given. The velocity and density (ρ) are in m/s and kg/m³, respectively. The standard error on the speeds of sound calculated from 5 simulations with different initial configurations are shown in brackets.

 

 

The longitudinal or transverse modes of acoustic velocities are different along the three different axes of symmetry. This is consistent with the results of elastic constants obtained from mechanical tests where RHO zeolite is found to be not perfectly isotropic. For a completely isotropic material, the Zener ratio must be equal to one [21], and for RHO zeolite this is not satisfied. The weak anisotropy in RHO zeolite thus attributes to the orientation-dependent acoustic velocity.

 

3.3. Example of nano-macro modelling strategy

 

This section exemplifies our nano-macro strategy for the acoustic modelling of a multiscale metamaterial. This strategy consists in using some of the previously calculated mechanical parameters as inputs of an upscaled model of wave propagation in an air-filled tube in which a periodic array of zeolite-made plates clamped onto a perfectly rigid solid frame is placed. Thus, two upscaling processes, i.e. one from the molecular to the mesoscopic scale, another one from the mesoscopic to the macroscopic scale; are considered. The second upscaling process has been reported in Ref. [22] and revealed that wave propagation in the material is governed by a macroscopic mass balance equation and a fluid flow constitutive law (see equations 2.43 and 2.44 in Ref. [22], respectively). The effective parameters of the material are the frequency-dependent effective compressibility C(ω) and visco-elasto-inertial permeability K(ω). The former is here calculated by using the JCAL model [23], with numerically calculated static thermal permeability and thermal characteristic length. On the other hand, K(ω) is calculated from the finite-element solution of a local fluid-structure interaction problem for which the divergence-free fluid velocity is governed by the Stokes equation coupled with the equation that governs the scalar out-of-plane displacement of the clamped zeolite-made plates, as shown by equations 2.24–2.31 in Ref. [22]. It is emphasised that i) the plates are modelled as Love-Kirchoff plates with plate modulus E / (1 − ν ²) and superficial density ρt, where t is the plate thickness, and ii) because of the low porosity and high apparent flow resistivity of the zeolite-made plates, only acoustic wave propagation through the solid frame of the zeolite is considered. It is also noted that normal pressure and temperature condition is considered in what follows.

 

Figure 4a) illustrates the geometry of the metamaterial, while Figure 4b) shows the magnitude and phase of K(ω) / K₀, with K₀ = K(ω → 0). At low frequencies, the plates can be considered as quasi-rigid and the effects of fluid viscosity dominate the acoustic behaviour. As the frequency increases, elasto-inertial effects become prominent. Noteworthy is that in between the frequencies f_a = 530 Hz and f_g = 600 Hz, which correspond to an anti-resonance and a resonance frequency, respectively; an atypical band emerges. In this band, the real part of the effective density [defined as = η / jωK where η is the dynamic viscosity of air] is negative while the real part of the effective compressibility is positive. This leads to the emergence of a sub-wavelength band gap (i.e. a frequency band where acoustic wave propagation is forbidden), as can be seen in Figures 4c) and 4d) where the speed of sound [defined as c(ω) = ω / k_c(ω)] and attenuation coefficient (defined as −ℑ(k_c (ω)), where is the effective wave number) exhibit slow and large values for frequencies in between f_a and f_g, respectively. The results show the potential and prospects of designing, from the molecular to the macroscopic scale, multiscale metamaterials for noise reduction applications.

 

4. CONCLUSIONS

 

This paper presented results towards the bridging of nanoscale and macroscale acoustics of porous solids. We have estimated the mechanical parameters and the acoustic properties of a prototypical zeolite material with cubic symmetry using molecular dynamics simulation. For the estimation of mechanical constants using MD, the zeolite structure was subjected to different mechanical tests, i.e. tensile and shear tests. From the tensile test, the mechanical constants C₁₁ and C₁₂ were determined by deforming the zeolite along the axial direction. The elastic constant C₄₄ was determined by applying a shear deformation to the upper surface of the zeolite while keeping the lower surface fixed to its position. Using the mechanical constants obtained from MD, we calculated the direction dependent longitudinal and transverse speeds of sound by solving Kelvin-Christoffel’s equation. In addition, the outcomes of the molecular simulation study were used as inputs of an upscaled model of wave propagation in a metamaterial comprising an air-filled tube and a periodic array of zeolite made plates clamped on to a rigid frame. We, thus, considered two upscaling processes, i.e. i) from molecular to mesoscopic scale and ii) from mesoscopic to macroscopic scale. This provided an example for our nano-macro modelling strategy. Finally, the results of this work showed the potential for designing, from the molecular to the macroscopic scale, multiscale metamaterials for noise reduction applications.

 

 

Figure 4: (a) Geometry of the metamaterial. (b) Normalised magnitude of the visco-elasto-inertial permeability. (c) Normalised real part of the effective speed of sound. (d) Attenuation coefficient. The inset plots show the normalised phase of the respective parameter.

 

ACKNOWLEDGEMENTS

 

This work, which was selected as an IRGA project from Univ. Grenoble Alpes and Universidad Austral de Chile, is supported by the French National Research Agency in the framework of the “Investissements d’Avenir” program (ANR-15-IDEX-02). Calculations were performed using the Froggy platform of the GRICAD infrastructure (https://gricad.univ-grenoblealpes.fr), which is supported by the Rhône-Alpes region (GRANT CPER07-13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the French Research Agency. Support from the Chilean National Agency for Research and Development (ANID) through FONDECYT Grant 1211310 is also acknowledged.

 

REFERENCES

 

1. L. Bocquet and E. Charlaix. Nanofluidics, from bulk to interfaces. Chem. Soc. Rev., 39:1073–1095, 2010.

2. B. Coasne, A. Galarneau, R. J.-M. Pellenq, and F. Di Renzo. Adsorption, intrusion and freezing in porous silica: the view from the nanoscale. Chem. Soc. Rev., 42:4141–4171, 2013.

3. B. Coasne. Multiscale adsorption and transport in hierarchical porous materials. New J. Chem., 40:4078–4094, 2016.

4. P. Van Der Voort, K. Leus, and E. De Canck. Introduction to Porous Materials. John Wiley & Sons, 2019.

5. M. L. Terranova, S. Orlanducci, and M. Rossi. Carbon nanomaterials for gas adsorption. CRC Press, 2012.

6. S. Valencia and F. Rey. New Developments in Adsorption/Separation of Small Molecules by Zeolites. Springer, 2020.

7. O. Coussy. Mechanics and physics of porous solids. John Wiley & Sons, 2011.

8. L. Brochard, M. Vandamme, and R. J.-M. Pellenq. Poromechanics of microporous media. Journal of the Mechanics and Physics of Solids, 60(4):606–622, 2012.

9. F.-X. Coudert. Responsive metal–organic frameworks and framework materials: under pressure, taking the heat, in the spotlight, with friends. Chemistry of Materials, 27(6):1905–1916, 2015.

10. B. Coasne et al. Poroelastic theory applied to the adsorption-induced deformation of vitreous silica. J. Phys. Chem. B, 118(49):14519–14525, 2014.

11. V. Van Speybroeck et al. Advances in theory and their application within the field of zeolite chemistry. Chem. Soc. Rev., 44(20):7044–7111, 2015.

12. Ch. Baerlocher and L. B. McCusker. Database of zeolite structures.

13. D. J. Evans and B. L. Holian. The Nose–Hoover thermostat. J. Chem. Phys., 83(8):4069–4074, 1985.

14. M. Parrinello and A. Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys., 52(12):7182–7190, 1981.

15. A. P. Thompson et al. LAMMPS – a flexible simulation tool… Comp. Phys. Comm., 271:108171, 2022.

16. A. Ghysels et al. Shape-selective diffusion of olefins in 8-ring solid acid microporous zeolites. J. Phys. Chem. C, 119(41):23721–23734, 2015.

17. R. W. Hockney and J. W. Eastwood. Computer simulation using particles. CRC Press, 2021.

18. F. Mouhat and F.-X. Coudert. Necessary and sufficient elastic stability conditions… Phys. Rev. B, 90(22):224104, 2014.

19. Allan F. Bower. Applied mechanics of solids. CRC Press, 2009.

20. R. Poloni and J. Kim. Predicting low-k zeolite materials. J. Mater. Chem. C, 2(13):2298–2300, 2014.

21. C. M. Zener and S. Siegel. Elasticity and anelasticity of metals. J. Phys. Chem., 53(9):1468–1468, 1949.

22. R. Venegas and C. Boutin. Acoustics of permeo-elastic materials. J. Fluid Mech., 828:135–174, 2017.

23. J.-F. Allard. Propagation of sound in porous media: Modeling sound absorbing materials, 1994.

 

---

¹ alan.sam@univ-grenoble-alpes.fr

² benoit.coasne@univ-grenoble-alpes.fr

³ rodolfo.venegas@uach.cl