FreeFem++ scripts for phononic crystals
This distribution is a companion to the book Phononic Crystals, Second Edition, by Vincent Laude (De Gruyter, 2020).
Please send requests to me (using the contact form). Requests for corrections or clarifications are gratefully expected and will be gladly answered; junk email will be treated as such.
This is version 2.0.0 dated July 23, 2020.
Download
ff++_phononic_laude_v2.zip (16.6 MBytes, download from Dropbox)
https://cloud.femto-st.fr/nextcloud/index.php/s/j7bmdZ8t7ae3PwF (alternative cloud repository)
Changes from version 1.1.0, dated April 10, 2018.
All scripts have been revised, improved, and checked to work with FreeFem++ 4.6 (July 2020). The default sparse solver now executes automatically on all the cores of the computing machine.
Changes from version 1.0.0
All scripts have been rechecked for consistency. The main correction is the definition of the perfectly matched layer (PML) for pressure wave problems. The main novelty is the use of a driver of the eigenvalue solver ARPACK that was not available previously: when the matrix B in the generalized eigenvalue problem Ax=λBx is not symmetrical (it should work with any version of FreeFem++ dated after February 2017). As a result, all cases are now covered by the use of the standard ARPACK eigenvalue solver and there is no direct need to use the SLEPC solver. Another progress is the use of the sequential MUMPS solver instead of UMFPACK for many of the scripts. MUMPS can make use of all available CPU cores without explicit MPI programming, while UMFPACK essentially executes on a single core even if you have a multicore machine.
How to use these files
The files in this archive are provided as companions to the book “Phononic crystals”. Though they are freely accessible, it is implied that they can be best understood if you actually are in possession of the book. Mathematical details and variational formulations, in particular, are not described here but referred to. All you will need is free software to be able to run the examples provided and reproduce many of the figures in the book. All script files can be edited with a simple text editor.
FreeFem++ (abbreviated ff++) is a finite element language developed by Frédéric Hecht and co-workers at Pierre et Marie Curie University, Paris, France. It is available for free and runs on most computers, including Windows, Linux, Unix, and Mac. You are strongly advised to install a working version first and to read the manual that comes with it before running the scripts provided as edp files. ff++ is a complete multiphysics partial differential equation solver with a syntax close to c++. Follow the explanations for installation on the ff++ install page. ff++ scripts are typically run inside a terminal, e.g.:
> FreeFem++ script.edp
Examples of gnuplot scripts are provided to plot band structures. Gnuplot is a powerful program to produce professional looking graphs in 2D and 3D. It was used throughout the book. Look for plt files. Output is generally to the vector eps (encapsulated postscript) format; images converted to pixelated png format are provided for the convenience of Windows users. gnuplot scripts are typically run inside a terminal, e.g.:
> gnuplot script.plt
Gmsh was used specifically for 3D mesh generation. Gmsh is software developed by Christophe Geuzaine and Jean-François Remacle. Some example geo files are provided to illustrate mesh creation by gmsh in the case of artificial crystals in 3D. Look for msh files for 3D meshes exported to the gmsh format that can be imported within ff++.
Note: 2D meshes are generated directly within the ff++ scripts. See also appendix B.
Basic scripts for acoustic waves (chapters 2 and 3)
These scripts are located in folder scattered_field and are devoted to the model of pressure waves: the variable is the spatial distribution of pressure, obeying a scalar Helmholtz equation (section 3.3).
Radiation from a line with prescribed pressure, modeling a simple source of pressure waves. Scripts src2d_source.edp and src2d_source_radiation.edp illustrate the use of a perfectly matched layer (PML) or of a radiation boundary condition, respectively (figures 3.10 and 3.11). The computation is performed at a given frequency (time-harmonic computation).
Finite sonic crystal of steel rods in water. Scripts src2d_grating.edp and scr2d_pc.edp consider the cases of one row and 4 rows of rods, illuminated by a line source (figure 2.18). Radiation at infinity is implemented with a PML.
Scattered field model (section 3.3.6). The same problems are reconsidered within the scattered field model, with an incident infinite plane wave at a given frequency. Scripts scf2d_disk.edp, scf2d_grating.edp, and scf2d_pc.edp consider the case of a single rod, of a single row, and of 4 rows of rods (figures 2.17 and 3.12).
Generating band structures for 2D sonic crystals (chapter 4)
These scripts are located in folder 2D_SC and combine the model of pressure waves with a formulation of the Bloch-Floquet theorem (section 4.1.5). Eigenvalues (for frequencies) and eigenvectors (for Bloch waves) are obtained with the Arpack solver.
Sonic crystal of rigid rods in air. Script band_air_rigid.edp solves for the band structure for the square lattice (figure 4.8), the hexagonal lattice (figure 4.10), and the honeycomb lattice (figure 4.11). The rods are not meshed and their boundary receives a “wall” condition (pressure is free, particle velocity is zero in a weak sense).
Sonic crystal of steel rods in water. Script band_water_steel.edp solves for the band structure for the square lattice (figure 4.13), the hexagonal lattice (figure 4.15), and the honeycomb lattice (figure 4.17). The rods are meshed. Steel is represented as an equivalent “fluid” medium with shear elastic waves neglected.
Other example scripts are given for the nylon rods in water case (illustrating the uploading of material parameters from a file; see appendix A) and for the case the inclusion must be described using a continuous variation of material constants (i.e., when defining materials by domains of the mesh is not practical).
Please note that if you wish to change the lattice type, you also need to edit by hand some lines.
A sample gnuplot bs.plt script file is provided. It generates the band structures in the figures cited previously.
Generating band structures for 3D sonic crystals (chapter 4)
These scripts are located in folder 3D_SC. The main difference with the 2D case is mesh generation using gmsh. Sample geometry (geo) files are given for the SC, BCC, and FCC lattices. Different cases without and with a filled or empty spherical inclusion are given. Meshes (msh files) can be produced by running gmsh on the geo files provided. Note that you will need to edit the geo files and generate meshes again if you change geometrical parameters (such as the lattice constant or the inclusion diameter).
Band structure calculation is illustrated on 4 different cases. First, the rigid sphere in air model is considered in script band_air_rigid.edp. Next, the air bubble in water (band_water_air.edp, figures 4.24 and 4.25) and the tungsten carbide bead in water (band_water_tungstencarbide.edp) are considered. Finally, the close-packing case is treated with continuous variation of the material constants (band_water_tungstencarbide_continuous.edp, figure 4.26).
A sample gnuplot bs.plt script file is given. It generates the band structures in the figures cited previously.
Bulk elastic waves and plate modes (chapter 5)
Scripts for generating slowness curves and wave surfaces in anisotropic elastic media, including piezoelectric media, are contained in folders bulk and plate_waves, respectively. Variational formulations are given in Appendix 5B.
Bulk elastic waves. Both the purely elastic and piezoelectric cases are illustrated. In either case, a 2D mesh composed of a single triangle and P0 Lagrange vector elements are used. The finite element matrices that result are very small, but full (i.e., not sparse); the Lapack solver is then used to obtain eigenvalues and eigenvectors. The example ff++ script bulk_slowness.edp treats the case of silicon; just edit line 25 to change the material file and run the script again (also see appendix A below). Slowness curves are generated by scanning the wavevector direction in space and solving for the phase velocity. Wave surfaces are further obtained by evaluating the energy velocity vector at each point of the slowness curves.
Bulk elastic waves in piezoelectric media. Exactly the same procedure as above is followed, except for the slightly generalized variational formulation. The example ff++ script bulk_piezo_slowness.edp treats the case of lithium niobate with cut (XY): crystallographic axis X is orthogonal to the figure, crystallographic axis Y is the horizontal axis.
Gnuplot script bulk.plt reproduces figures 5.4 (a and b), and 5.5.
Plate waves. Elastic waves guided along a plate with 2 plane parallel surfaces (including Lamb waves) are obtained from a 1D-like 2D-mesh. The reason is that there are no 1D-meshes in ff++. Script band_plate.edp considers the case of silicon and reproduces the band structure in figure 5.14a. As before, it is not difficult to generate the similar band structure for another elastic material, or to extend the script to piezoelectric media.
Gnuplot script bs_plate.plt reproduces figure 5.14a.
Generating band structures for 2D phononic crystals (chapter 6)
These scripts are located in folder 2D_PC and combine the model of (vector) elastic waves with a formulation of the Bloch-Floquet theorem (section 6.1.4). Eigenvalues (for frequencies) and eigenvectors (for Bloch waves) are obtained with the Arpack solver.
The ff++ scripts provided are similar to those for 2D sonic crystals (chapter 4). Remember that a few lines have to be edited by hand in case you wish to change the lattice type. Most band structure examples in chapter 6 can be reproduced with these scripts. They include the following list.
|
Epoxy matrix – steel rods |
Figures 6.4a, 6.6, 6.7 |
band_epoxy_steel.edp |
|
Steel matrix – epoxy rods |
Figures 6.8, 6.9, 6.10 |
band_steel_epoxy.edp |
|
Steel matrix –cylindrical holes |
Figures 6.11, 6.12, 6.13 |
band_steel_void.edp |
|
Silicon matrix – cylindrical holes |
Figure 6.14 |
band_si_void.edp |
|
TiO2 matrix – cylindrical holes |
Figure 6.15 |
band_teo2_void.edp |
|
Lithium niobate – cylindrical holes; cuts (ZX), (YX), and (XY) |
Figures 6.16, 6.17, 6.18 |
band_lno_void.edp and band_lnonp_void.edp |
A sample gnuplot bs.plt script file is provided. It generates the band structures in the figures cited previously.
Generating band structures for 3D phononic crystals (chapter 6)
These scripts are located in folder 3D_PC and combine the model of (vector) elastic waves with a formulation of the Bloch-Floquet theorem (section 6.1.4).
Band structure calculation is illustrated on 2 different cases for the FCC lattice and steel spheres in a solid matrix, in script band_steel_polyester.edp. First, a moderate filling fraction is considered and the parallelepiped primitive unit cell with 6 faces is used (figure 6.18). Then the close-packed condition is considered for an epoxy matrix, requiring the use of the 12 faces Wigner-Seitz cell (figure 6.19). Remember that a few lines have to be edited by hand in case you wish to change the lattice type.
A sample gnuplot bs.plt script file is given. It generates the band structures in the figures cited previously.
Generating band structures for phononic crystal slabs (chapter 7)
These scripts are located in folder PCslab and combine the model of (vector) elastic waves with a formulation of the Bloch-Floquet theorem (section 6.1.4). Eigenvalues (for frequencies) and eigenvectors (for Bloch waves) are obtained with the Arpack solver.
Band structure calculation is illustrated for different lattices and for either a holey membrane (script bs_slab.edp) or a solid-solid composition (script bs_slab_fill.edp). Figures 6.4a, 6.5, 6.6, 6.8, and 6.9 can be reproduced with these scripts. Remember that a few lines have to be edited by hand in case you wish to change the lattice type. A selection of geo files for use with gmsh and the meshes generated from them are given.
A sample gnuplot bs.plt script file is given. It generates the band structures in the figures cited previously.
Coupling of acoustic and elastic waves (chapter 8)
These scripts are located in folder 2D_ae and implement the model of acoustic-elastic coupling introduced in section 8.1.
Band structure calculation is illustrated for 2D crystals of steel rods in water with various lattices (script band_water_steel.edp). Figures 8.4, 8.5, and 8.6 can be reproduced. Remember that a few lines have to be edited by hand in case you wish to change the lattice type.
A sample gnuplot bs.plt script file is given. It generates the band structures in the figures cited previously.
Complex band structure (chapter 9)
The computation of complex band structures with the finite element method is illustrated for both sonic and phononic crystals.
In the sonic crystal case, in folder 2D_SC_CBS, the variational formulation of section 9.2.2 is used in script bc_fluid.edp generating figure 9.3a and 9.3b. The corresponding two gnuplot files are provided.
In the phononic crystal case, in folder 2D_PC_CBS, the variational formulation of section 9.3.2 is used in script cbs.edp generating figures 9.7a, 9.7b, 9.8a and 9.8b. The script includes the evaluation of the polarization contents of evanescent Bloch waves. The corresponding four gnuplot files are provided.
Wave propagation in tube with grafted resonator
These scripts are located in folder 1D_LR. They illustrate in the 2D case the physics of resonators grafted along a waveguide. There is no correspondance to figure in the book (though figure 10.7 is a 3D version of the same problem).
The scripts illustrate the following basic problems:
- Obtaining the acoustic modes of a closed cavity (cavity.edp) and for a single cavity grafted on a waveguide (cavity+wg.edp). Eigenvalues (for frequencies) and eigenvectors are obtained with the Arpack solver.
- Computing the transmission through a tube by using radiation boundary conditions (source.edp). This is a time-harmonic computation with frequency considered as a parameter.
- Computing the complex band structure with the variational formulation of section 9.2.2 (cbs.edp). A sample gnuplot cbs.plt script file is given.
Super-cells for waveguide and cavity problems
The waveguide problem is illustrated in the case of a square-lattice sonic crystal of steel rods in water. Script files in folder WG_SC reproduce figures 9.14 and 9.15 (for a total of 5 sub-figures). The complex band structures for the Gamma-X and the Gamma-M directions are specifically computed (script bc_fluid.edp) and are superposed with the classical band structure (script band_fluid_2d.edp) in figures 9.14a and 9.14b. Next, the 6*1 super-cell includes (or not) a central homogeneous water region to define the waveguide. The complex band structures for the Gamma-X direction with/without the central core (script bc_fluid_wg.edp) are shown in figures 9.15a and 9.15b. Finally, the classical band structure (script band_fluid_wg.edp) is shown for comparison in figures 9.15c.
The cavity problem is illustrated in the 2D case of a square-lattice phononic crystal of tungsten rods in silicon. Script file square_lattice_cavity_modes_W_Si.edp in folder cavity considers the case of a 7*7 super-cell with a central missing rod. The script computes eigenfrequencies and outputs the modal shape of a cavity mode whose frequency falls inside a complete phononic band gap.
Appendix A: units and material files
Material files are located in folder constants. Different file extensions are used depending on the material type.
- Fully anisotropic data are formatted in text files with extension const. Material constants include (among others) the elastic, phonon viscosity, dielectric, and piezoelectric tensors. This format should be used for all anisotropic (including crystalline) materials.
- Isotropic elastic materials (without viscosity) can be presented via text files with extension isolid. The elastic tensor is computed from the longitudinal and shear velocities in this case.
- Fluids (treated as linear acoustic materials) can be presented via text files with extension fluid. Only the mass density and the elastic modulus are necessary in this case.
Material constants have physical units and all scripts comply with the S.I. system of units. In addition, numerical values for material constants are internally rescaled to avoid numerical under- and overflows. With very large matrices, even though they are sparse, the condition number is highly dependent of the choice for this internal representation. The general rule that is followed is that mesh dimensions are normalized to the lattice constant so that all lengths are close to unity. Multiplicative factors are used for material constants in a consistent manner. For instance, elastic constants are divided by 1011 since most stiff solid materials have c11 rigidities of the order of a few 100 GPa. Some of the numerical factors used in the scripts are listed below.
|
Quantity |
LMTI |
Units |
Numerical factor |
Remark |
|
Elastic constants, pressure, stiffness |
(-1,1,-2,0) |
Pa |
1011 |
Arbitrary choice |
|
Mass density |
(-3,1,0,0) |
kg/m3 |
103 |
Arbitrary choice |
|
Dielectric constants, permittivity |
(-3,-1,4,2) |
F/m |
10-11 |
Arbitrary choice |
|
Piezoelectric constants, electric displacement |
(-2,0,1,1) |
C/m2 |
1 |
Implied by constitutive relation for piezoelectricity |
|
Slowness |
(-1,0,1,0) |
s/m |
10-4 |
Implied by Christofell's equation |
|
Velocity |
(1,0,-1,0) |
m/s |
104 |
Idem |
|
Strain |
(0,0,0,0) |
- |
1 |
Adimensional |
|
Electric field |
(1,1,-3,-1) |
V/m |
1011 |
Implied by constitutive relation for electric displacement |
|
Length, displacement |
(1,0,0,0) |
m |
10-6 |
Arbitrary choice |
|
Time |
(0,0,1,0) |
s |
10-10 |
Definition of slowness |
|
Angular frequency |
(0,0,-1,0) |
Hz |
1010 |
Idem |
LMTI: Length, Mass, Time, Intensity
Appendix B: utility scripts
Utility ff++ scripts are provided in folder functions. The extension for source files included within edp files is *.idp.
- LoadMaterial.idp contains a collection of functions returning material constants read in specially formatted text files. The material constants provided are in no way guaranteed to be accurate but were instead collected by the author from various sources.
- 2DCrystalMesh.idp contains functions to build 2D meshes of the square, hexagonal, and honeycomb lattices with/without inclusion(s). The case of dual meshes for acoustic-elastic coupling is also included.
- kvectors.idp contains functions to sample the usual symmetry directions of Brillouin zones, for many lattices in both 2D and 3D.``