Application Exchange

Phononic Band-Gap Structure Eigenfrequency Analysis
Filename Size
BandGapUnitCell1.mph 19.9 MB
Download all files (Zip-archive) ~ 15.9 MB

Phononic Band-Gap Structure Eigenfrequency Analysis

Nagi Elabbasi, Veryst Engineering

Phononic crystals are artificially manufactured structures, or materials, with periodic constitutive or geometric properties designed to influence the characteristics of mechanical wave propagation. They can be engineered to isolate vibration in a certain frequency range. Vibration in that frequency range, called a band gap, is attenuated by a mechanism of wave interferences within the periodic system.

To illustrate, we created this model involving a 2D periodic structure with a unit cell composed of a stiff inner core and a softer outer matrix material, designed to have a band gap around 60-70 kHz. We applied Bloch boundary conditions to constrain the displacements of the unit cell, and set up a complex Eigenfrequency Study with a Parametric Sweep spanning the wave vectors that represent the boundaries of the irreducible Brillouin zone. When we plot the wave propagation frequencies for all wave numbers, a band gap appears as a region where no wave propagation branches exist.

User Comments

Ahmed Nagaty
Feb 11, 2016 at 11:37am UTC

please i am using comsol 4.4 and this module can't open with ....if you could help me in that by a pdf or any thing please ,because really i need it ..
thank you

Nagi Elabbasi
Entry submitter
Feb 11, 2016 at 11:36pm UTC

Unfortunately I do not have an older version of the model. If you read this related recent COMSOL blog post ( you should hopefully be able to recreate the model in version 4.4.

Nagi Elabbasi
Entry submitter
Feb 12, 2016 at 2:26pm UTC

If it helps here is a copy of the parameter definitions:

k 0.5 0-3 spans irreducible Brillouin zone (0-1 horizontal, 1-2 vertical, 2-3 diagonal)
L1 0.01[m] Periodic cell length
L2 L1/2.5 Size of inner material
kx if(k<1,pi/L1*k,if(k<2,pi/L1,(3-k)*pi/L1)) X component of wave number
ky if(k<1,0,if(k<2,(k-1)*pi/L1,(3-k)*pi/L1)) Y component of wave number
E1 2e9[Pa] Modulus of Material 1
E2 200e9[Pa] Modulus of Material 2
Rho1 1e3[kg/m^3] Density of Material 1
Rho2 8e3[kg/m^3] Density of Material 2
Poisson1 0.45 Poisson ratio of Material 1
Poisson2 0.34 Poisson ratio of Material 2
ShearMod1 E1/2/(1+Poisson1) Shear modulus of Material 1
Factor1 L1/sqrt(ShearMod1/Rho1) Frequency normalization factor

Ahmed Nagaty
Feb 13, 2016 at 7:43pm UTC

thank you

Saeid Jamilan
Feb 14, 2016 at 12:23am UTC

Thanks for your help in advance. How it is possible to get band gap diagram for TM or TE ploarization? I am using COMSOL 5.1

Nagi Elabbasi
Entry submitter
Feb 16, 2016 at 8:51pm UTC

You should be able to use a very similar approach for photonic band gap calculation. Check out COMSOL Application ID 798 called Bandgap Analysis of a Photonic Crystal.

Saeid Jamilan
Feb 17, 2016 at 1:17am UTC

Many thanks for your response and help. I used that instructions and could simulate. The problem is that it is simulating for TE Polarization. I do not know how to adjust COMSOL settings to extract band gap data for TM polarization.

Ahmed Nagaty
Feb 19, 2016 at 2:10pm UTC

how can i control the filling fraction in the phononic crystal to know its effect on the band gap ??

said said
Feb 20, 2016 at 2:14am UTC

please what is the version which are you worked this model?

Ahmed Nagaty
Feb 22, 2016 at 7:03pm UTC

version 5.2

Yuning Guo
Mar 4, 2016 at 9:35pm UTC

Hello Sir, does it possible to obtain equifrequency contour of Phononic crystal in Comsol? If possible, could you tell me how to achieve it? Thanks a lot.

Ahmed Nagaty
Mar 6, 2016 at 8:56pm UTC

hi , i have aproblem in calculation of phononic band gap for 2D phononic crystal , i have calculated it but i have a oroblem in the wave vector axis , i want to be MX (coordinates of the first brillion zone not numbers)...

can any one help me in that

Thibaut Jacqmin
Mar 14, 2016 at 9:11pm UTC

Thanks for the nice model. Maybe it would be nice if you could export it as a *.m (matlab) file so that even people who don't have the right version of Comsol can see how to implement it. I thank you very much in advance.

Zhaoxuan Zhang
Mar 17, 2016 at 8:51am UTC

Hi, Mr. Elabbasi,
I'm using COMSOL 5.1 version and I cannot open your mpm file, could you upload another mph file?

Edoardo Belloni
Mar 24, 2016 at 9:03am UTC

Dear Mr. Elabbasi,
thanks for the nice model. I have a question: in the Floquet Periodicity Boundary Conditions, you set (kx,kx) and (ky,ky) in direction X and Y. Why this choice?
I think it is (kx,ky) and (kx,ky) in both directions. Also comparing with Comsol Application ID 798 (Bandgap Analysis of a Photonic Crystal), my choice for the k parameter is applied.
Thanks for the help!

Clinton Potts
Mar 26, 2016 at 3:24am UTC

Dr. Elabbasi

This example was extremely helpful. I was curious how you applied the excitation frequency in your GIF example. It is an extremely interesting demonstration but I cannot seem to replicate it.

said said
Apr 1, 2016 at 9:54am UTC

Comment removed by said said.

said said
Apr 1, 2016 at 9:55am UTC

Mr. Nagi hello please i want to calculate the band structure of a hexagonal mesh how exchange parameters in your loop?

Nagi Elabbasi
Entry submitter
Apr 1, 2016 at 1:54pm UTC

Thank you Edoardo for pointing out the mistake. Yes it should be (kx, ky) for both boundary conditions instead of (kx, kx) for one and (ky, ky) for the other. It does not affect the results though in any way because the dot product of k_F and (r_dst – r_src) will zero out the wrong terms. Thank you for pointing it out though, I will upload a corrected version soon.

Nagi Elabbasi
Entry submitter
Apr 1, 2016 at 1:55pm UTC

Dear Clinton, you are referring to the animations I showed in this Blog ( I applied the excitation as a horizontal displacement at the left boundary in a Frequency Domain analysis, and I ran that analysis at 27, 67.5 and 88 kHz.

Nagi Elabbasi
Entry submitter
Apr 1, 2016 at 1:58pm UTC

Dear Said, you follow the same approach basically. You model the unit cell with the proper Floquet boundary conditions, and solve over the wave numbers that span the Irreducible Brillouin Zone of the hexagonal structure. The IBZ still has three “sides” as in this example so the changes to the model should be relatively minor.

lei xiao
May 10, 2016 at 8:12am UTC

Comment removed by lei xiao.

lei xiao
May 10, 2016 at 8:14am UTC

Dr. Elabbasi
This example is really helpful . Thanks a lot.
I want to repeat your work in this example; but I do not know how to apply the horizontal displacement excitation and I did not find similar examples in the documentation. Would you mind share your whole model?

Zhuhua Tan
May 16, 2016 at 7:17am UTC

Hi Nagi,

It is a greatful example for beginners!

Could you tell me why define kx ky as following:
kx if(k<1,pi/L1*k,if(k<2,pi/L1,(3-k)*pi/L1)) X component of wave number
ky if(k<1,0,if(k<2,(k-1)*pi/L1,(3-k)*pi/L1)) Y component of wave number


Nagi Elabbasi
Entry submitter
May 22, 2016 at 11:49pm UTC

Thanks Zhuhua! If you follow the logic of the IF statements you will find that for 0<k<1 kx and ky define the “horizontal” edge of the IBR, for 1<k<2 they define the vertical edge, and for 2<k<3 they define the diagonal edge.

工 王
Jun 12, 2016 at 8:25am UTC


Omar Najera
Nov 7, 2016 at 8:49pm UTC

this model is only for the diagonal direction or for all directions in the crystal?

Henrik Asendorf
Dec 6, 2016 at 9:22am UTC


I am trying to use your model to calculate the photonic bandgap of a pillar in air. In your model there are parameters defining Youngs modulus and Poissons ratio of the materials. Since air does not have a Possoin ratio how can I modify this model to make it work for my system?

Xukun Su
Dec 25, 2016 at 9:50am UTC

Dear Mr. Elabbasi,
Thanks a lot for this model and all of your patience! It helped a lot for a beginner!

I got one question. What is 'Frequency normalization factor' you defined in Parameters? And why is it defined by the expression 'L1/sqrt(ShearMod1/Rho1)'?

Thanks a lot!!

And Merry Christmas!

Kiumars Aryana
Mar 9, 2017 at 3:01pm UTC

Comment removed by Kiumars Aryana.

Kiumars Aryana
Mar 10, 2017 at 10:08pm UTC

Dear Nagi,

I tried to replicate the animation but I don't get smooth shapes as yours. I was wondering if it has sth to do with my mesh or it's something else? by the way, is that stress field that you animated?

Thank you,

汶菊 徐
Jun 3, 2017 at 9:31am UTC

Thank you very much!

Han He
Aug 14, 2017 at 9:12pm UTC

Hi Nagi,
It is a great material. I should have read it earlier.
However, for the kx:
kx if(k<1,pi/L1*k,if(k<2,pi/L1,(3-k)*pi/L1)) X component of wave number
if 1<k<2 is along the verical edge, why the kx is pi/L1 instead of 0?

Login to comment on this entry.