Computation of the input acoustic impedance of a cylindrical tube

bandrieux

Hello,

I am doing an acoustic study using Code_Aster. The purpose of this study is the computation of the input acoustic impedance of a cylindrical tube. I am using Salome-Meca with Code_Aster on Debian and my problems are linked with the use of Code_Aster.

I used the testcase AHLV100 (Waveguide with anechoic output) and I tried to adapt it. I meshed the air in my tube with Salome, I defined the input face ('ENTREE'), the output ('SORTIE'), the lateral faces ('FACES'). I print the values of pressure and intensity of a group of nodes at the input of the tube. (with the intensity I can have the particle velocity and then compute the impedance Z=P/V)

I need to have a pressure=0 at the output, a speed=0 on all other faces. I am not sure of what I am doing at this point with my script.

First I define my mesh (DEFI_MATERIAU, AFFE_MATERIAU, AFFE_MODELE). I choose the boundary conditions (AFFE_CHAR_ACOU), maybe I am not doing something right at this step?

Then I use MACRO_MATR_ASSE, CALC_VECT_ELEM, ASSE_VECTEUR, DYNA_LINE_HARM, CALC_ELEM and CALC_NO in order to build a table.

With this script the results are wrong, the pressure and the intensity are 0 over a certain frequency.

Maybe I am doing something wrong with the boundary conditions.

Can someone help me?

Thank you,

Benoit

The code :

#ETUDE TUBE



DEBUT();

MAIL=LIRE_MAILLAGE(FORMAT='MED',);

AIR=DEFI_MATERIAU(FLUIDE=_F(RHO=1.3,
                            CELE_C=('RI',343.0,0.0,),),);

CHAMPMAT=AFFE_MATERIAU(MAILLAGE=MAIL,
                       AFFE=_F(TOUT='OUI',
                               MATER=AIR,),);

GUIDE=AFFE_MODELE(MAILLAGE=MAIL,
                  VERIF='MAILLE',
                  AFFE=_F(TOUT='OUI',
                          PHENOMENE='ACOUSTIQUE',
                          MODELISATION='3D',),);

CHARACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,
                        VITE_FACE=_F(GROUP_MA='ENTREE',
                                     VNOR=('RI',0.014,0.,),),);

IMPEACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,
                        IMPE_FACE=_F(GROUP_MA='SORTIE',
                                     IMPE=('RI',445.9,0.,),),);

VITACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,
                       VITE_FACE=_F(GROUP_MA='FACES',
                                    VNOR=('RI',0.014,0.0,),),);

MACRO_MATR_ASSE(MODELE=GUIDE,
                CHAM_MATER=CHAMPMAT,
                CHARGE=IMPEACOU,
                NUME_DDL=CO('NUM'),
                MATR_ASSE=(_F(MATRICE=CO('MATASK'),
                              OPTION='RIGI_ACOU',),
                           _F(MATRICE=CO('MATASM'),
                              OPTION='MASS_ACOU',),
                           _F(MATRICE=CO('MATASI'),
                              OPTION='AMOR_ACOU',),),);

VECTELEM=CALC_VECT_ELEM(OPTION='CHAR_ACOU',
                        CHAM_MATER=CHAMPMAT,
                        CHARGE=CHARACOU,);

VECTASS=ASSE_VECTEUR(VECT_ELEM=VECTELEM,
                     NUME_DDL=NUM,);

#-----------------------CALCUL HARMONIQUE-------------------------------

LIFREQ=DEFI_LIST_REEL(DEBUT=10,
                      INTERVALLE=_F(JUSQU_A=10000,
                                    NOMBRE=999,),);

PRESHARM=DYNA_LINE_HARM(MATR_MASS=MATASM,
                        MATR_RIGI=MATASK,
                        MATR_AMOR=MATASI,
                        LIST_FREQ=LIFREQ,
                        EXCIT=_F(VECT_ASSE=VECTASS,
                                 COEF_MULT=-1.,
                                 PHAS_DEG=90.,
                                 PUIS_PULS=1,),);

PRESHARM=CALC_ELEM(reuse =PRESHARM,
                   RESULTAT=PRESHARM,
                   TOUT_ORDRE='OUI',
                   OPTION='PRAC_ELNO',);

PRESHARM=CALC_NO(reuse =PRESHARM,
                 RESULTAT=PRESHARM,
                 TOUT_ORDRE='OUI',
                 OPTION='PRAC_NOEU',);

PRESSAC=CREA_TABLE(RESU=_F(RESULTAT=PRESHARM,
                           NOM_CHAM='PRAC_NOEU',
                           TOUT_ORDRE='OUI',
                           TOUT_CMP='OUI',
                           GROUP_NO='ENTRN',),
                   TYPE_TABLE='TABLE',);

IMPR_TABLE(TABLE=PRESSAC,
           FORMAT='TABLEAU',
           NOM_PARA=('NOM_CHAM','NUME_ORDRE','PRES_R','PRES_I','DB',),
           SEPARATEUR=' ,',);

PRESHARM=CALC_ELEM(reuse =PRESHARM,
                   RESULTAT=PRESHARM,
                   TOUT_ORDRE='OUI',
                   OPTION='INTE_ELNO',);

PRESHARM=CALC_NO(reuse =PRESHARM,
                 RESULTAT=PRESHARM,
                 TOUT_ORDRE='OUI',
                 OPTION='INTE_NOEU',);

INTAC=CREA_TABLE(RESU=_F(RESULTAT=PRESHARM,
                         NOM_CHAM='INTE_NOEU',
                         TOUT_ORDRE='OUI',
                         TOUT_CMP='OUI',
                         GROUP_NO='ENTRN',),
                 TYPE_TABLE='TABLE',);

IMPR_TABLE(TABLE=INTAC,
           FORMAT='TABLEAU',
           NOM_PARA=('NOM_CHAM','NUME_ORDRE','INTX_R','INTY_R','INTZ_R','INTX_I','INTY_I','INTZ_I',),
           SEPARATEUR=' ,',);

FIN();

CLF

This is from memory and my cellphone...

Pressure release is not obvious in CA.

Pres=0 does not work. It produces zero pressure at the bc node, but that is it.

Try to put PHI=0 where you want pressure release. Phi zero implies zero acoustic displacement, ie max particle velocity and hence,pressure release.

The other way to go is to use the peseanteur element to impose a pressure release bc.

Claes

CLF

I notice you are a student. Congrats. CA is quite a good code for vibroacoustics

Claes

YvonGoth

Hello,

I notice that your max frequency is pretty high. Is your mesh set for such a high frequency?.

With the lambda/6 rule, you should use a mesh with a max size of 6 mm.

I don't see any boundary condition setting p=0. You should have something like

CHARACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,

PRES_IMPO=_F(GROUP_MA='SORTIE', PRES=('RI',0.,0.,)))

The comment of CLF about setting the 0 pressure condition by PHI=0 is right, but for the fluide_structure modele. For the ACOUSTIQUE phenomena, the direct setting of the pressure should work.

To set a boundary condition of fixed impedance type, your code seems OK.

Bon courage,

Yvon

bandrieux

Hello,

Thank you for your answers. I already tried to set the pressure with :

PRESACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,

PRES_IMPO=_F(GROUP_MA='SORTIE', PRES=('RI',0.,0.,)))

I think I am doing something wrong with MACRO_MATR_ASSE or CALC_VECT_ELEM. When I choose the boundary conditions with AFFE_CHAR_ACOU should I had something in MACRO_MATR_ASSE or in CALC_VECT_ELEM?

I want to have an anechoic output so I think I should keep CHARGE=IMPEACOU in MACRO_MATR_ASSE, but for CALC_VECT_ELEM should I keep only CHARGE=CHARACOU (with CHARACOU corresponding to the velocity=0 at the input?

Is there another solution to impose boundary walls for example for the tube or the velocity=0 is enough?

Here is my code and I attached my project :

#ETUDE TUBE

DEBUT();

MAIL=LIRE_MAILLAGE(FORMAT='MED',);

AIR=DEFI_MATERIAU(FLUIDE=_F(RHO=1.3,
                            CELE_C=('RI',343.0,0.0,),),);

CHAMPMAT=AFFE_MATERIAU(MAILLAGE=MAIL,
                       AFFE=_F(TOUT='OUI',
                               MATER=AIR,),);

GUIDE=AFFE_MODELE(MAILLAGE=MAIL,
                  VERIF='MAILLE',
                  AFFE=_F(TOUT='OUI',
                          PHENOMENE='ACOUSTIQUE',
                          MODELISATION='3D',),);

CHARACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,
                        VITE_FACE=_F(GROUP_MA='ENTREE',
                                     VNOR=('RI',0.014,0.,),),);

IMPEACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,
                        IMPE_FACE=_F(GROUP_MA='SORTIE',
                                     IMPE=('RI',445.9,0.,),),);

VITACOU=AFFE_CHAR_ACOU(MODELE=GUIDE,
                       VITE_FACE=_F(GROUP_MA='FACES',
                                    VNOR=('RI',0.014,0.0,),),);

PRESACOU=AFFE_CHAR_ACOU(MODELE=GUIDE, 
			     PRES_IMPO=_F(GROUP_MA='SORTIE',
					      PRES=('RI',0.,0.,),),);

MACRO_MATR_ASSE(MODELE=GUIDE,
                CHAM_MATER=CHAMPMAT,
                CHARGE=IMPEACOU,
                NUME_DDL=CO('NUM'),
                MATR_ASSE=(_F(MATRICE=CO('MATASK'),
                              OPTION='RIGI_ACOU',),
                           _F(MATRICE=CO('MATASM'),
                              OPTION='MASS_ACOU',),
                           _F(MATRICE=CO('MATASI'),
                              OPTION='AMOR_ACOU',),),);

VECTELEM=CALC_VECT_ELEM(OPTION='CHAR_ACOU',
                        CHAM_MATER=CHAMPMAT,
                        CHARGE=CHARACOU,);

VECTASS=ASSE_VECTEUR(VECT_ELEM=VECTELEM,
                     NUME_DDL=NUM,);

#-----------------------CALCUL HARMONIQUE-------------------------------

LIFREQ=DEFI_LIST_REEL(DEBUT=10,
                      INTERVALLE=_F(JUSQU_A=10000,
                                    NOMBRE=999,),);

PRESHARM=DYNA_LINE_HARM(MATR_MASS=MATASM,
                        MATR_RIGI=MATASK,
                        MATR_AMOR=MATASI,
                        LIST_FREQ=LIFREQ,
                        EXCIT=_F(VECT_ASSE=VECTASS,
                                 COEF_MULT=-1.,
                                 PHAS_DEG=90.,
                                 PUIS_PULS=1,),);

PRESHARM=CALC_ELEM(reuse =PRESHARM,
                   RESULTAT=PRESHARM,
                   TOUT_ORDRE='OUI',
                   OPTION='PRAC_ELNO',);

PRESHARM=CALC_NO(reuse =PRESHARM,
                 RESULTAT=PRESHARM,
                 TOUT_ORDRE='OUI',
                 OPTION='PRAC_NOEU',);

PRESSAC=CREA_TABLE(RESU=_F(RESULTAT=PRESHARM,
                           NOM_CHAM='PRAC_NOEU',
                           TOUT_ORDRE='OUI',
                           TOUT_CMP='OUI',
                           GROUP_NO='ENTRN',),
                   TYPE_TABLE='TABLE',);

IMPR_TABLE(TABLE=PRESSAC,
           FORMAT='TABLEAU',
           NOM_PARA=('NOM_CHAM','NUME_ORDRE','PRES_R','PRES_I','DB',),
           SEPARATEUR=' ,',);

PRESHARM=CALC_ELEM(reuse =PRESHARM,
                   RESULTAT=PRESHARM,
                   TOUT_ORDRE='OUI',
                   OPTION='INTE_ELNO',);

PRESHARM=CALC_NO(reuse =PRESHARM,
                 RESULTAT=PRESHARM,
                 TOUT_ORDRE='OUI',
                 OPTION='INTE_NOEU',);

INTAC=CREA_TABLE(RESU=_F(RESULTAT=PRESHARM,
                         NOM_CHAM='INTE_NOEU',
                         TOUT_ORDRE='OUI',
                         TOUT_CMP='OUI',
                         GROUP_NO='ENTRN',),
                 TYPE_TABLE='TABLE',);

IMPR_TABLE(TABLE=INTAC,
           FORMAT='TABLEAU',
           NOM_PARA=('NOM_CHAM','NUME_ORDRE','INTX_R','INTY_R','INTZ_R','INTX_I','INTY_I','INTZ_I',),
           SEPARATEUR=' ,',);

FIN();

I don't think it is a problem of mesh, I tried with finer mesh and the results are the same (for 40Hz for example PRES_R=10^-20!).

Benoit

Attached file:
0404TUBE.tar, 778240kb

YvonGoth

Hello,

Null velocity is the natural boundary condition for the F.E. acoustic problem, so it has not to be imposed on rigid walls.

Please see the .comm file attached for examples of use of different acoustic boundary conditions.

Regard,

Yvon

Attached file:
cylacou.comm, 3235kb

bandrieux

Hello,

Yhank you very much for your reply, it helps me to understand how Code_Aster works.

I don't understand where I make a mistake but when I try to use your examples or mine I still have the same problem. The pressure field is ok for the first frequencies (I can observe the field with the post-pro in Salome-meca) but then for the next frequencies I have nothing (fields are 0!).

I meshed my tube with the mesh module of Salome, run my study with Aster in Salome and export in a table format RESULTAT (with a table CALC_ELEM then CALC_NO). I tried with the format MED that's how I saw the fields in Salome but the problem is the same.

I hope I will find a solution, if you have an idea...

Thank you very much,

Benoit

bandrieux

Hi,

Thanks to Claes and Yvon I achieve to compute the acoustic impedance of my tube. I think my mistakes were more on the mesh and on the sizes of my geometry.

Now I need to add the losses on the walls the tube. I use the operator FORMULE to add the formula that modelize the losses:

LVOR = FORMULE(VALE='sqrt((2*1.85e-5)/((2*pi*FREQ)*RHO))',


NOM_PARA=('FREQ','RHO',),);

Ys = FORMULE(VALE='0.5*(1-i)*((2*pi*FREQ)/(RHO*CELE_C**2))*(LVOR(FREQ,RHO)*(sin(pi/2))**2+(1.4-1)*(LVOR(FREQ,RHO)/(sqrt(0.708))))',
NOM_PARA=('FREQ','LVOR','RHO','CELE_C',),);

Zparois = FORMULE(VALE='1/Ys(FREQ,LVOR,RHO,CELE_C)',

NOM_PARA='Ys',);

I want to use the impedance of the walls calculated (Zparois) in AFFE_CHAR_ACOU for my list of frequencies. I don't understand how to do this. The wall impedance depends on the frequency.

I try to have a look at the test-case zzzz100a but I can't see something like this.

If you know how to do this it will be very helpful...

Thank you very much,

Benoit

CLF

Hi Benoit

Congratulations to your success. Always good to get more CA acoustics users. :)

Yvon knows Python much better than me (I know nothing), but the way to handle frequency dependence in Code Aster is to use a Python FOR loop that loops the Impedance definition for each frequency and then executes it.This is doable as CA runs in Python.

If you take a look at the forum and search for FOR loops and Python, I am sure you will find some usable examples.

Also, watch out when using the complex wave speed. Yvon reported a bug in CA for this a while back, but I am not sure it has been fixed yet.

Sincerely

Claes

bandrieux

Hi,

Thank you very much for your reply. I think I have a problem with MACRO_MATR_ASSE, I made a for loop including from AFFE_CHAR_ACOU until the end of the script (necessary in order to have a wall impedance dependent on the frequency).

How is it possible in

'MACRO_MATR_ASSE(MODELE=ModelDef,

CHAM_MATER=Champmat,

CHARGE=Cac[n],

NUME_DDL=CO('NUM[n]'),

MATR_ASSE=(_F(MATRICE=CO('AMK[n]'), OPTION='RIGI_ACOU',),

_F(MATRICE=CO('AMM[n]'), OPTION='MASS_ACOU',),),);'

to update the number of iterations n (in AMK[n], AMM[n] and NUM[n], they are between quotes)?

I have the error :

CR phase d'initialisation

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! erreur de syntaxe, Erreur de nom: name 'NUM' is not defined ligne 54 !

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

fin CR phase d'initialisation

Thank you very much,

Benoit

The complete code

DEBUT();



Cac=[None]*500;
Ldv=[None]*500;
Vass=[None]*500;
Pr=[None]*500;
Prac=[None]*500;
Inac=[None]*500;
pas=2;

Case1=LIRE_MAILLAGE(FORMAT='MED',);

Air=DEFI_MATERIAU(FLUIDE=_F(RHO=1.3,
                            CELE_C=('RI',343.0,0.0,),),);

Champmat=AFFE_MATERIAU(MAILLAGE=Case1,
                       AFFE=_F(TOUT='OUI',
                               MATER=Air,),);

ModelDef=AFFE_MODELE(MAILLAGE=Case1,
                     VERIF='MAILLE',
                     AFFE=_F(TOUT='OUI',
                             PHENOMENE='ACOUSTIQUE',
                             MODELISATION='3D',),);

for n in range(1,50):
	frequence=n*pas;
	
#	Zparois=n*440;

	Cac[n]=AFFE_CHAR_ACOU(MODELE=ModelDef,
                        PRES_IMPO=_F(GROUP_MA='End2',
                                     PRES=('RI',0.0,0.0,),),
#			IMPE_FACE=_F(GROUP_MA='FACES',
#                                     IMPE=('RI',Zparois,0.0,),),	     
                        VITE_FACE=_F(GROUP_MA='End1',
                                     VNOR=('RI',0.014,0.0,),),);

	MACRO_MATR_ASSE(MODELE=ModelDef,
                CHAM_MATER=Champmat,
                CHARGE=Cac[n],
                NUME_DDL=CO('NUM[n]'),
                MATR_ASSE=(_F(MATRICE=CO('AMK[n]'),
                              OPTION='RIGI_ACOU',),
                           _F(MATRICE=CO('AMM[n]'),
                              OPTION='MASS_ACOU',),),);

	Ldv[n]=CALC_VECT_ELEM(OPTION='CHAR_ACOU',
                        CHAM_MATER=Champmat,
                        CHARGE=Cac[n],);

	Vass[n]=ASSE_VECTEUR(VECT_ELEM=Ldv[n],
                     NUME_DDL=NUM[n],);

	ARK[n]=COMB_MATR_ASSE(COMB_R=_F(PARTIE='REEL',
                                  MATR_ASSE=AMK[n],
                                  COEF_R=1.,),);

	ARM[n]=COMB_MATR_ASSE(COMB_R=_F(PARTIE='REEL',
                                  MATR_ASSE=AMM[n],
                                  COEF_R=1.,),);


	Pr[n]=DYNA_LINE_HARM(MATR_MASS=AMM,
			MATR_RIGI=AMK,
                       	FREQ=frequence,
                       	EXCIT=_F(VECT_ASSE=Vass[n],
				COEF_MULT=1.0,
                               	PHAS_DEG=90.0,
                               	PUIS_PULS=1,),);

	Pr[n]=CALC_ELEM(reuse =Pr[n],
                   RESULTAT=Pr[n],
                   TOUT_ORDRE='OUI',
                   OPTION='PRAC_ELNO',);

	Pr[n]=CALC_NO(reuse =Pr[n],
                 RESULTAT=Pr[n],
                 TOUT_ORDRE='OUI',
                 OPTION='PRAC_NOEU',);
		 
	Prac[n]=CREA_TABLE(RESU=_F(RESULTAT=Pr[n],
                           NOM_CHAM='PRAC_NOEU',
                           TOUT_CMP='OUI',
                           GROUP_NO='ENTRN',),
                   TYPE_TABLE='TABLE',);
		   
	IMPR_TABLE(TABLE=Prac[n],
           FORMAT='TABLEAU',
           FORMAT_C='REEL_IMAG',
           NOM_PARA=('RESULTAT','PRES_R','PRES_I'),
           SEPARATEUR=' ,',); 


	Pr[n]=CALC_ELEM(reuse =Pr[n],
                   RESULTAT=Pr[n],
                   TOUT_ORDRE='OUI',
                   OPTION='INTE_ELNO',);

	Pr[n]=CALC_NO(reuse =Pr[n],
                 RESULTAT=Pr[n],
                 TOUT_ORDRE='OUI',
                 OPTION='INTE_NOEU',);
		 
	Inac[n]=CREA_TABLE(RESU=_F(RESULTAT=Pr[n],
                          NOM_CHAM='INTE_NOEU',
                          TOUT_CMP='OUI',
                          GROUP_NO='ENTRN',),
                  TYPE_TABLE='TABLE',);

	IMPR_TABLE(TABLE=Inac[n],
           FORMAT='TABLEAU',
           FORMAT_C='REEL_IMAG',
           NOM_PARA=('RESULTAT','INTX_R','INTY_R','INTZ_R','INTX_I','INTY_I','INTZ_I',),
           SEPARATEUR=' ,',);
	   
	n=n+1;

FIN();

delmas

Before MACRO_MATR_ASSE

try this :

NUM[n]=CO('NUM_%d' % (n))

and write in MACRO_MATR_ASSE

NUME_DDL=NUM[n],

bandrieux

Hi,

thank you for the help I have a program working with the python for loop but I can't see any difference between the model with losses and without losses (the losses should reduce the height of the peaks for the high frequencies).

Claes what do you mean with the use of the complex acoustic speed? I can't find anything on the forum about that. If I don't use the complex wave speed I have the same results.

In my computation I always have the real part of the pressure equals to zero as well as the real parts of the acoustic intensity and I don't know if it is normal (when I compute the acoustic impedance I have a real part equals to zero).

What is the influence of the amplitude of the velocity source I use? For example if I want to have an equivalent of a plane wave of 1Pa of amplitude, what would it be with a velocity source?

Thank you very much,

Benoit

CLF

Hi Benoit

True, there should be a difference between cases with complex and zero imaginary wave speed, but you are correct, in CA there is no difference. This bug was reported a while back by Yvon.

In this bug report Yvon also shows a workaround to the problem.

About velocity source - you only find a simple relation for free wave conditions, where pressue/particle velocity = rho*c (= 407 for air at 20c).

Sincerely

Claes

CLF

Hi Again

This is Yvon's post I refer to:

http://www.code-aster.org/forum2/viewtopic.php?id=13394

Claes

ThomasDESOZA

Hi,

I'm not familiar with acoustics but are you suggesting that it is possible to input a complex celerity (with non zero imaginary part) and yet this won't affect the result ?

If so this needs to be fixed.

TdS

bandrieux

Hi,

and thank you Claes for the help.

Thomas, I only try with a zero imaginary part but Yvon in the post (http://www.code-aster.org/forum2/viewtopic.php?id=13394) said that the sound celerity was declared as REAL so I understand that a complex speed won't work in CA.

Benoit

delmas

Merci pour le feedback. Ce problème a été remonté dans notre bugtracker et devrait être corrigé d'ici la fin du mois de juin.

delmas

Ce problème sera corrigé cette semaine dans la version 12.0.14.