Accueil | Cast3M

* fichier :  rotor4.dgibi
************************************************************************
************************************************************************
*
************************************************************************      
*                                                                      *
*  Mots-clés : Vibrations, calcul modal, machines tournantes,          *
*              poutre, modes complexes                                 *
*                                                                      *
*                                                                      *
*  Test de GYROSCOPIQUE, AMOR, CAMPBELL pour les elements de poutre    *
*  Etude d'une machine tournante dans le repere inertiel (ou fixe)     *
*  Elements: TUYAU et SECTION (fibre)                                  *
*            POUT, TIMO en exemple                                     *
*                                                                      *
*   1-Influence des déformations de cisaillement                       *
*   2-Etude de la stabilité avec un amortissement corotatif            *
*                                                                      *
*   Modélisation d'un tube                                             *
*                                                                      *
*   Auteur: Didier Combescure Juillet 2003                             *
*                                                                      *    
************************************************************************ 
*GRAPH = 'Y';
GRAPH = 'N';
*
*L1 = 0.50 ;
L1 = 1.20;                                                
Re1 = 0.1 ;
epa1 = 0.0012;                                                             
E1 = 2.e+11 ;                                                           
Nu1 = 0.2 ;                                                             
Rho1 = 1600. ;
Visc1 = 0.00001*E1; 
*Visc1 = 0.0*E1; 
*
* Description de la section
*
OPTI 'DIME' 2 'ELEM' QUA4;
P0S = 0. 0.;
P1S = (Re1 - (0.5*epa1)) 0.;
P2S = (Re1 + (0.5*epa1)) 0.; 
*
P12S = D 1 P1S P2S;
SURS = P12S ROTA 50 360. P0S; 
ELIM 0.0001 SURS;
*                                                                     
OPTI 'DIME' 3 'ELEM' SEG2;                                                                     
*                                                         
*
P0 = 0. 0. 0.;
P1 = L1  0.  0. ;
GeoRot = D 100 P0 P1; 
*
* OMEG=1
* On donne une vitesse de rotation de 1 rad/s
*
* MODELE SECTION
*
kappa = 0.52;
*kappa = 1.0;
*kappB = 1.00;
*
MOSEC1 = MODE SURS MECANIQUE ELASTIQUE QUAS;
MASEC1 = MATE MOSEC1 'YOUNG' E1 'NU' Nu1 'RHO' Rho1 
               ALPY kappa ALPZ kappa 'VISQ' Visc1;
*
*
* Elément TIMO - FIBRE
*    
Mod1 = MODE GEOROT 'MECANIQUE' ELASTIQUE SECTION TIMO;                                          
Mat1 = MATE Mod1 'MODS' MOSEC1 'MATS' MASEC1 
           'OMEG' (2.*pi) VECT (0. 0. 1.) ;      
*
*
* Elément TUYAU
*    
*        
Mod2 = MODE GEOROT 'MECANIQUE' TUYA;                                          
Mat2 = MATE Mod2 'YOUNG' E1 'NU' Nu1 'RHO' Rho1
        'RAYO' Re1 'EPAI' EPA1 'OMEG' (2.*pi) 'VISQ' Visc1;
*
* Elément POUT
*    
*        
Mod3 = MODE GEOROT 'MECANIQUE' POUT;                                          
Mat3 = MATE Mod3 'YOUNG' E1 'NU' Nu1 'RHO' Rho1 'SECT' (pi * ((Re1 ** 2)-((Re1 - epa1)**2)) ) 'INRY' (pi * ((Re1 ** 4)-((Re1 - epa1)**4))/4.) 'INRZ' (pi * ((Re1 ** 4)-((Re1 - epa1)**4))/4.) 'TORS' (pi * ((Re1 ** 4)-((Re1 - epa1)**4))/2.) 'OMEG' (2.*pi) 'VISQ' Visc1;    
*
*
* Timoshenko
*
Mod4 = MODE GEOROT 'MECANIQUE' TIMO;                                          
Car4 = CARA Mod4 'SECT' (pi * ((Re1 ** 2)-((Re1 - epa1)**2)) )
                 'INRY' (pi * ((Re1 ** 4)-((Re1 - epa1)**4))/4.) 
                 'INRZ' (pi * ((Re1 ** 4)-((Re1 - epa1)**4))/4.)
                 'TORS' (pi * ((Re1 ** 4)-((Re1 - epa1)**4))/2.)
                 'SECY' (kappa*(pi * ((Re1 ** 2)-((Re1 - epa1)**2))) )
                 'SECZ' (kappa*(pi * ((Re1 ** 2)-((Re1 - epa1)**2))) )
                 'OMEG' (2.*pi);    
* 
*
Mat4 = MATE Mod4 'YOUNG' E1 'NU' Nu1 'RHO' Rho1 'VISQ' Visc1;  
Mat4 = Mat4 et Car4;
*
MESH3D = POUT2MAS Mod1 Mat1 'GAUSS';
SI (EGA GRAPH 'Y');
  OEIL = 10. 10. 10.;
  ELIM 0.0001 MESH3D;
  Titre 'Maillage poutre avec section fibre';
  trac OEIL MESH3D CACH;
  trac OEIL (MESH3D et GeoRot);
FINSI;
*
*                                                                      
* Cond. limites : 2 appuis 
*
*
Blo2 = (BLOQ 'UZ' P0) ET
       (BLOQ 'UY' P0) ET
       (BLOQ 'UX' P0);   
 
Blo3 = (BLOQ 'UZ' P1) ET
       (BLOQ 'UY' P1) ET
       (BLOQ 'UX' P1);                             
Blo1 = BLOQ 'RX' Georot;                           
*                                                                      
MMass1 = MASS Mod1 (Mat1) ;                                   
MRigi1 = RIGI Mod1 (Mat1) ;   
*
MMass2 = MASS Mod2 (Mat2) ;                                   
MRigi2 = RIGI Mod2 (Mat2) ;   
*
MMass3 = MASS Mod3 (Mat3) ;                                   
MRigi3 = RIGI Mod3 (Mat3) ;   
*
MMass4 = MASS Mod4 (Mat4) ;                                   
MRigi4 = RIGI Mod4 (Mat4) ; 
*  
* Matrice de couplage gyroscopique
*
 CG1 = GYROSCOPIQUE Mod1 Mat1;
 CG2 = GYROSCOPIQUE Mod2 Mat2;
 CG3 = GYROSCOPIQUE Mod3 Mat3;
 CG4 = GYROSCOPIQUE Mod4 Mat4;
* 
CDAM1 CDAH1 = 'AMOR' Mod1 Mat1 'COROTATIF';
CDAM2 CDAH2 = 'AMOR' Mod2 Mat2 'COROTATIF';
CDAM3 CDAH3 = 'AMOR' Mod3 Mat3 'COROTATIF';
CDAM4 CDAH4 = 'AMOR' Mod4 Mat4 'COROTATIF';                                 
*                                                                       
TBasR1 = VIBR 'PROCHE' (prog 100.) (lect 2)
           (MRigi1 ET Blo1 ET Blo2 et Blo3) 
             MMass1  ;
TBasR2 = VIBR 'PROCHE' (prog 100.) (lect 2)
           (MRigi2 ET Blo1 ET Blo2 et Blo3) 
             MMass2  ;
TBasR3 = VIBR 'PROCHE' (prog 100.) (lect 2)
           (MRigi3 ET Blo1 ET Blo2 et Blo3) 
             MMass3  ;
TBasR4 = VIBR 'PROCHE' (prog 100.) (lect 2)
           (MRigi4 ET Blo1 ET Blo2 et Blo3) 
             MMass4  ;
*  
TMod1 = TBasR1 . 'MODES' ;
TMod2 = TBasR2 . 'MODES' ;               
TMod3 = TBasR3 . 'MODES' ;
TMod4 = TBasR4 . 'MODES' ;
NbMod3 = (DIME TMod3) - 2 ;
TInd3 = INDEX (TMod3 . 1) ;
*
SI (EGA GRAPH 'Y');
 oeil = 100. 1. 0. ;
 REPETER bbb1 NbMod3 ;
  i = &bbb1 ;
  Frq1 = TMod1 . i . 'FREQUENCE' ;
  Def1 = TMod1 . i . 'DEFORMEE_MODALE' ;
  Frq2 = TMod2 . i . 'FREQUENCE' ;
  Def2 = TMod2 . i . 'DEFORMEE_MODALE' ;
  Frq3 = TMod3 . i . 'FREQUENCE' ;
  Def3 = TMod3 . i . 'DEFORMEE_MODALE' ;
  Frq4 = TMod4 . i . 'FREQUENCE' ;
  Def4 = TMod4 . i . 'DEFORMEE_MODALE' ;
  TITRE Frq1 ;
  trac oeil (defo georot def1);
  TITRE Frq2 ;
  trac oeil (defo georot def2);
  TITRE Frq3 ;
  trac oeil (defo georot def3);
  TITRE Frq4 ;
  trac oeil (defo georot def4);
 FIN bbb1 ;
FINSI;
*
* On teste pour le modele a fibre
*
  MAmor1 = CDAM1;
  MAmor2 = CDAM2;
  MAmor3 = CDAM3;
  MAmor4 = CDAM4;
*
*-Projection des matrices assemblées sur la base réelle
*
  MMasP1 = PJBA TBasR1 MMass1 ;
  MRigP1 = PJBA TBasR1 MRigi1 ;
  MAmoP1 = PJBA TBasR1 MAmor1 ;
  MAmHP1 = PJBA TBasR1 CDAH1;
  Cg1P = PJBA TBasR1 Cg1 ;
*
  MMasP2 = PJBA TBasR2 MMass2 ;
  MRigP2 = PJBA TBasR2 MRigi2 ;
  MAmoP2 = PJBA TBasR2 MAmor2 ;
  MAmHP2 = PJBA TBasR2 CDAH2;
  Cg2P = PJBA TBasR2 Cg2 ;
*
  MMasP3 = PJBA TBasR3 MMass3 ;
  MRigP3 = PJBA TBasR3 MRigi3 ;
  MAmoP3 = PJBA TBasR3 MAmor3 ;
  MAmHP3 = PJBA TBasR3 CDAH3;
  Cg3P = PJBA TBasR3 Cg3 ;
*
  MMasP4 = PJBA TBasR4 MMass4 ;
  MRigP4 = PJBA TBasR4 MRigi4 ;
  MAmoP4 = PJBA TBasR4 MAmor4 ;
  MAmHP4 = PJBA TBasR4 CDAH4;
  Cg4P = PJBA TBasR4 Cg4 ;
*
* Vitesse de rotation (0 Hz)
* 
Omeg1 =0.;
*
* Calcul des modes complexes avec VIBC - Bernoulli
*
   BasCi1 = VIBC MMasP4 (MRigP4 et (Omeg1*MAmHP4)) 
                        (MAmoP4 et (Omeg1*Cg4P)) TBasR4 ;
*
  ModBasC = BasCi1.'MODES';
  NbMod1 = (DIME ModBasC) - 2;
*
FRQR = TABLE;
FRQI = TABLE;
FRQ = TABLE;
AMOR= TABLE;
*
mess 'Frequences et amortissement au repos';
* 
  REPETER bbb2 NbMod1 ;
   i2 = &bbb2 ;
   FrqR1 = BasCi1.'MODES'.i2.'FREQUENCE_REELLE' ;
   FrqI1 = BasCi1.'MODES'.i2.'FREQUENCE_IMAGINAIRE' ;
   DefR1 = BasCi1.'MODES'.i2.'DEFORMEE_MODALE_REELLE';
   DefI1 = BasCi1.'MODES'.i2.'DEFORMEE_MODALE_IMAGINAIRE';
   FRQR.i2 = FrqR1;
   FRQI.i2 = FrqI1;
   FRQ.i2 = ((FrqR1 ** 2) + (FrqI1 ** 2)) ** 0.5 ;
   AMOR.i2   = FrqI1 /(FRQ.i2 );
   mess 'Valeur propre réelle';
   list (FrqR.i2); 
   mess 'Valeur propre imaginaire';
   list (FrqI.i2);
   mess 'Fréquence';
   list (FRQ.i2); 
   mess 'Amortissement';
   list (AMOR.i2);
  FIN bbb2;
*
*
* Vitesse de rotation (1000 Hz)
* 
Omeg1 =1000.;
*
* Calcul des modes complexes avec VIBC - Timoshenko
*
   BasCi1 = VIBC MMasP4 (MRigP4 et (Omeg1*MAmHP4)) 
                        (MAmoP4 et (Omeg1*Cg4P)) TBasR4;
*
  ModBasC = BasCi1.'MODES';
  NbMod1 = (DIME ModBasC) - 2;
*
FRQR = TABLE;
FRQI = TABLE;
FRQ = TABLE;
AMOR= TABLE;
* 
mess 'Frequences et amortissement a 1000 Hz';
* 
  REPETER bbb2 NbMod1 ;
   i2 = &bbb2 ;
   FrqR1 = BasCi1.'MODES'.i2.'FREQUENCE_REELLE' ;
   FrqI1 = BasCi1.'MODES'.i2.'FREQUENCE_IMAGINAIRE' ;
   DefR1 = BasCi1.'MODES'.i2.'DEFORMEE_MODALE_REELLE';
   DefI1 = BasCi1.'MODES'.i2.'DEFORMEE_MODALE_IMAGINAIRE';
   FRQR.i2 = FrqR1;
   FRQI.i2 = FrqI1;
   FRQ.i2 = ((FrqR1 ** 2) + (FrqI1 ** 2)) ** 0.5 ;
   AMOR.i2   = FrqI1 /(FRQ.i2);
   mess 'Valeur propre réelle';
   list (FrqR.i2); 
   mess 'Valeur propre imaginaire';
   list (FrqI.i2);
   mess 'Fréquence';
   list (FRQ.i2); 
   mess 'Amortissement';
   list (AMOR.i2);
  FIN bbb2;
*
* Calcul du diagramme de Campbell - Bernoulli
*
TAB1 = TABLE;
TAB1.'BASE_MODALE' = TBasR3;
TAB1.'MASS_PROJ' = MMasP3;
TAB1.'RIGI_PROJ' = MRigP3;
TAB1.'GYRO_PROJ' = Cg3P;
TAB1.'AMOR_PROJ' = MAmoP3;
TAB1.'KROT_PROJ' = MAmHP3; 
TAB1.'AFFICHAGE' = VRAI;
*PRFREQ = (prog 0. pas 2.5 100.)*100/(2.*pi);
PRFREQ = (prog 0. pas 25 100.)*100/(2.*pi);
CAMPBELL TAB1 PRFREQ;
*
EVFRb1 = (TAB1. 1). 'FREQUENCE_MODULE' ;
EVFRb2 = (TAB1. 4). 'FREQUENCE_MODULE' ;
EVAMb1 = (TAB1. 1). 'AMORTISSEMENT' ;
EVAMb2 = (TAB1. 4). 'AMORTISSEMENT' ;
*
SI (EGA GRAPH 'Y');
  dess (EVFRb1  et EVFRb2 );
  dess (EVAMb1 et EVAMb2);
FINSI;
*
* Calcul du diagramme de Campbell - Timoshenko
*
TAB1 = TABLE;
TAB1.'BASE_MODALE' = TBasR4 ;
TAB1.'MASS_PROJ' = MMasP4;
TAB1.'RIGI_PROJ' = MRigP4;
TAB1.'GYRO_PROJ' = Cg4P;
TAB1.'AMOR_PROJ' = MAmoP4;
TAB1.'KROT_PROJ' = MAmHP4;
TAB1.'AFFICHAGE' = VRAI;
*PRFREQ = (prog 0. pas 2.5 100.)*100/(2.*pi);
PRFREQ = (prog 0. pas 25. 100.)*100/(2.*pi);
CAMPBELL TAB1 PRFREQ;
*
EVFRT1 = ((TAB1. 1). 'FREQUENCE_MODULE') COUL rouge ;
EVFRT2 = ((TAB1. 4). 'FREQUENCE_MODULE') COUL rouge ;
EVAMT1 = ((TAB1. 1). 'AMORTISSEMENT') COUL rouge ;
EVAMT2 = ((TAB1. 4). 'AMORTISSEMENT') COUL rouge ;
*
EvFreOme = evol manu (prog 0. 1600.) (prog 0. (1600.));
EvFreZer = evol manu (prog 0. 1600.) (prog 0. 0.);
*
SI (EGA GRAPH 'Y');
 'DESS' (EVFRT1  et EVFRT2 et EvFreOme)
 'TITR' 'Frequences' 'TITX' 'Freq rot (Tour/s)'
 'TITY' 'Frequences (Hz)';
 'DESS' (EVFRT1  et EVFRT2)
 'TITR' 'Frequences' 'TITX' 'Freq rot (Tour/s)'
 'TITY' 'Frequences (Hz)';
 'DESS' (EVAMT1 et EVAMT2 et EvFreZer)
 'TITR' 'Amortissement' 'TITX' 'Freq rot (Tour/s)'
 'TITY' 'Amortissement (%)'; 
*
 dess (EVFRT1  et EVFRT2 et EVFRB1  et EVFRB2 et EvFreOme);
 dess (EVAMT1 et EVAMT2 et EVAMB1 et EVAMB2); 
*
finsi;
*
mess 'Frequences à 0 Hz et 1000 Hz';
mess 'Poutre de Bernoulli';
FR1b = extr (IPOL  (prog 0.01 1000.) (extr EVFRB1  absc)
                                (extr EVFRB2  ordo)) 1; 
AM1b = extr (IPOL  (prog 0.01 1000.) (extr EVAMB1   absc)
                                (extr EVAMB1   ordo)) 1; 
FR2b = extr (IPOL  (prog 0.01 1000.) (extr  EVFRB1 absc)
                                (extr EVFRB2  ordo)) 2; 
AM2b = extr (IPOL  (prog 0.01 1000.) (extr EVAMB1   absc)
                                (extr EVAMB1  ordo)) 2; 
mess 'Poutre de Timoshenko';
FR1T = extr (IPOL  (prog 0.01 1000.) (extr EVFRT1  absc)
                                (extr EVFRT2  ordo)) 1; 
AM1T = extr (IPOL  (prog 0.01 1000.) (extr EVAMT1   absc)
                                (extr EVAMT1   ordo)) 1; 
FR2T = extr (IPOL  (prog 0.01 1000.) (extr  EVFRT1 absc)
                                (extr EVFRT2  ordo)) 2; 
AM2T = extr (IPOL  (prog 0.01 1000.) (extr EVAMT1   absc)
                                (extr EVAMT1  ordo)) 2; 
 
*
* On norme 
*
EVFRb1N = EVOL MANU  (extr EVFRb1 absc)
                     ((extr EVFRb1 ordo)/FR1b );
EVFRb2N = EVOL MANU  (extr EVFRb1 absc)
                     ((extr EVFRb2 ordo)/FR1b );
EVFRT1N = EVOL MANU  (extr EVFRT1 absc) 
                     ((extr EVFRT1 ordo)/FR1T );
EVFRT2N = EVOL MANU  (extr EVFRT1 absc)
                    ((extr EVFRT2 ordo)/FR1T );
*
SI (EGA GRAPH 'Y');
  TA= table;
  ta.1 = 'TIRR';
  ta.2 = 'TIRR';
  titre 'Diagramme de Campbell (Frequence = f(Fréquence rotation en Hz)';
  DESS (EVFRb1 et EVFRb2 et EVFRT1 et EVFRT2) ta ybord 0. 4000.;
  DESS (EVFRb1 et EVFRb2 et EVFRT1 et EVFRT2) ta ;
  titre 'Diagramme de Campbell (F/F0 = f(Fréquence rotation en Hz)';
  DESS (EVFRb1N et EVFRb2N et EVFRT1N et EVFRT2N) ta ybord 0. 2.;
  DESS (EVFRb1N et EVFRb2N et EVFRT1N et EVFRT2N) ta;
FINSI;
*
*
FR1Th = 787.3;
FR2Th = 812.7;
AM1Th = 2.48E-02;
AM2Th = -0.559E-02;
Errfrq = 0.01*FR1Th;
Erramo = 0.01*AM1Th;
*
SI ((ABS(FR1T - FR1Th) > Errfrq ) OU
    (ABS(FR2T - FR2Th) > Errfrq ) OU
    (ABS(AM1T - AM1Th) > Erramo ) OU
    (ABS(AM2T - AM2Th) > Erramo ));
   ERRE 5;
SINON;
   ERRE 0; 
FINSI;
*
FIN;
*