Télécharger vibr12.dgibi

Retour à la liste

* fichier :  vibr12.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*                                                                      *
*  Mots-cles : Vibrations, calcul modal, modes complexes,              *
*              interaction fluide-structure, instabilite               *
*                                                                      *
*  TEST : VIBR12                                                       *
*                                                                      *
*  CALCUL D INSTABILITE FLUIDE-ELASTIQUE SOUS ECOULEMENT AXIAL CONFINE * 
*  2 CONFIGURATIONS : ENCASTRE-LIBRE ou AA                             *
*                                                                      *
*  Ref : GIBERT Page 409-410                                           *
*        GREGORY, PAIDOUSSIS, 1966, "Unstable oscillation of tubular   *
*        cantilevers conveying fluid "                                 *
*                                                                      *
*                               ^ Z                                    *
*  ->                           |                                      *
*  -> V    +--------------------+--------------------+  --> X          *
*  ->      P0                                       P2                 *
*                                                                      *
*                                                                      *
************************************************************************  
 
************************************
*            OPTIONS               *
************************************
OPTI DIME 3 ELEM SEG2 ;
 
* Conditions limites : 'EL' ou 'AA'
CL = mot 'AA'   ;   
CL = mot 'EL'   ;  
prefix = CHAI 'vibr12-' CL;
 
* traces
* GRAPH = VRAI ;
* OPTI TRAC PSC EPTR 10 POTR 'HELVETICA_16' 'FTRA' (CHAI prefix'.ps');
GRAPH = FAUX ;
 OPTI TRAC PSC EPTR 10 POTR 'HELVETICA_16' 'FTRA' (CHAI prefix'.ps');
 
 
 
************************************
*   Caracteristiques du systeme    *
************************************
 
*--- Materiau Acier
Rho_A = 8100.  ;        
You_A = 2.0e11 ;     
Nu_A  = 0.3    ;   
 
*--- Longueurs : Longueur Tube = 2L
L  = .5             ;
 
*--- Diametre int. Tube
Di   = 10.0e-3    ;
 
*--- Rapport Beta = Mf / (Mf+ms)
Beta = 0.2        ;
 
*--- choix d'une methode de calcul pour :
*    | Rhof : Masse volumique fluide
*    | De : diametre externe 
*    sachant que : |  ms         = Mf * (1-Beta) / Beta
*                  | (De^2-Di^2) = Rhof*Di^2/Rhos * (1-Beta)/Beta
Meth = 1 ;
SI (EGA Meth 1) ;
*   On calcul le Diametre externe De pour respecter le rapport Beta
*   a masse volumique fluide fixe
    Rhof  = 1000. ;
    De2  = ( (((Rhof/Rho_A) * ( (1-Beta) / Beta ) )) + 1 ) * (Di*Di)    ;
    De   = De2**0.5   ;
SINON ;
*   On calcul Rhof pour respecter Beta ;
*   a diametre exterieur fixe
    De   = 12e-3      ;
    Rhof = (((( (De*De) - (Di*Di)) / ((1-Beta)/Beta) ) ) / (Di*Di) ) *Rho_A ;
FINSI ;
 
*--- Epaisseur du tube
Ep = (De - Di) / 2 ;
 
Mess 'Di=' Di ' De=' De  'Rhof=' Rhof ;
 
*--- Section Tube
S_T = (pi/4)*((De**2)-(Di**2)) ;      
 
*--- Moments quadratiques 
Ie_T= (pi/64)*((De**4)-(Di**4)) ;  
 
*--- Masse ajoutee lineique interne
S_I    = (pi/4.0)*(Di**2) ; 
Me     = Rhof*S_I    ;
Rho_eq = Me/S_T         ;
 
* Nbre de modes reels et amortissement modaux en % (cf operateur AMOR)
Nbmod  = 4 ;
L_AMOR = PROG Nbmod * 0.0 ;
 
* Facteurs pour adimensionner [Gibert,p.408] :
*  V : VR = RVR*V  
*  w : OMEG = Rw*w 
RVR =  (2.0*L) * (( ((Rho_A*S_T) + Me) / (You_A*Ie_T)  )**0.5) ;
Rw  =  (2.0*L)* RVR ;
* Autre adimensionnement [Axisa,Paidoussis]
si (ega CL 'EL');
RVR =  (2.0*L) * (( Me / (You_A*Ie_T)  )**0.5) ;
finsi;
 
* Juste pour verif
Beta_c= Rhof*S_I / ((Rhof*S_I) + (Rho_A*S_T)) ;
 
*--- Vitesse axial du fluide
pas_VR = .1 ;
VRmax  = 12.  ;
VRmin  = 0.  ;
LVR     = PROG VRmin PAS pas_VR VRmax ;
LVa     = (1.0 * LVR )  /  RVR ;
 
mess '  Rho A  Rho_eq   Beta Beta_c  Rw  RVr De'  ;
mess '  ' Rho_A '  ' Rho_eq '  ' Beta '  ' Beta_c '  ' Rw '  ' RVR '  ' De;
 
 
**************************************
*         Maillage du systeme        *
**************************************
 
P0 = (-1.0*L) 0. 0. ;
P2 = P0 PLUS ((2.0*L)  0. 0.) ;
 
* Lfin : longueur objectif d'1 element
* L1   : longueur reelle
*Lfin = 0.001 ;
Lfin = 0.0025 ;
nL   = ENTI 'PROCHE' (2.0*L / Lfin);
Tube = P0 DROI nL P2  ;
* Longueur reelle des elements
* L1 = (MESURE Tube) / nL ; 
L1 = 2.0*L / nL;
 
**************************************
*        Modele et Materiau          *
**************************************
 
Modl1 = MODE Tube MECANIQUE ELASTIQUE ISOTROPE 'TUYA' ;
Matr1 = MATE Modl1 'YOUN' You_A 'NU' Nu_A 'RHO' (Rho_A + Rho_eq) 'EPAI' Ep 'RAYO' (De/2);
 
 
**************************************
*   Masses, Rigidites ... du systeme *
**************************************
 
Mstru = MASS Modl1 Matr1 ;
Kstru = RIGI Modl1 Matr1 ;
 
*--- Conditions aux limites : Probleme plan EL ou EA 
Cl1 = BLOQ 'UX' Tube ;
Cl2 = BLOQ 'UY' Tube ;
Cl3 = BLOQ 'RX' Tube ;
Cl4 = BLOQ 'RZ' Tube ;
* EL : Encastre-Libre
SI (EGA CL 'EL') ;
    Cl5 = BLOQ 'RY' P0   ;
SINON ; 
* AA : Appuye-Appuye
  SI (EGA CL 'AA') ;
    Cl5 = BLOQ 'UZ' P2   ;
  SINON ;
    MESS 'CL doit etre AA ou EL !';
    ERRE 21;
  FINSI ;
FINSI ;
Cl6 = BLOQ 'UZ' P0   ;
Kstru = Kstru ET Cl1 ET Cl2 ET Cl3 ET Cl4 ET CL5 ET CL6 ;
 
 
*****************************************
*             Analyse modale            *
*****************************************
 
*Basr = VIBR 'PROC' (PROG 1.) (LECT Nbmod) Kstru Mstru ;
Basr = VIBR 'IRAM' 1. Nbmod Kstru Mstru ;
Basr.'MODES' = NNOR (Basr.'MODES') (MOTS 'UZ') ;
Basrm = Basr.'MODES' ;
 
OPTI OEIL (0.0 -1000.0 0.0) ;
POSTVIBR Basr ;
 
 
***********************************
*    Rigidites suplementaires    *
***********************************
 
PAfu1 = (PROG    -.5        0.          0.       (L1/10.)         .5        0.            0.        (-1*L1/10.)) ;
PAfu2 = (PROG    0.         -.5      (-1*L1/10.)    0.            0.         .5        (L1/10.)          0.)  ;
PAfu3 = (PROG    0.      (L1/10.)       0.          0.            0.     (-1*L1/10.) (-1*L1*L1/60.)      0.)  ;
PAfu4 = (PROG (-1*L1/10.)   0.          0.          0.       (L1/10.)       0.            0.        (-1*L1*L1/60.)) ;
PAfu5 = (PROG    -.5        0.          0.       (-1*L1/10.)       .5       0.            0.        (L1/10.)) ;
PAfu6 = (PROG    0.         -.5      (L1/10.)       0.            0.         .5      (-1*L1/10.)         0.)  ;
PAfu7 = (PROG    0.      (-1*L1/10.) (L1*L1/60.)    0.            0.     (L1/10.)         0.             0.)  ;
PAfu8 = (PROG (L1/10.)      0.          0.       (L1*L1/60.) (-1*L1/10.)    0.            0.             0.)  ;
PAfu    = PAfu1 ET PAfu2 ET PAfu3 ET PAfu4 ET PAfu5 ET PAfu6 ET PAfu7 ET PAfu8 ;
Afu     = MANU 'RIGI' 'TYPE' 'AMORTISSEMENT' 
          Tube (MOTS 'UY' 'UZ' 'RY' 'RZ') 'QUEL' PAfu ;
 
PKfu1 = (PROG (-6./(5.*L1))     0.           0.           -1.1      (6./(5.*L1))     0.            0.           -.1) ;
PKfu2 = (PROG    0.         (-6./(5.*L1))    1.1           0.          0.         (6./(5.*L1))      .1          0.)  ;
PKfu3 = (PROG    0.              .1       (-2.*L1/15.)     0.          0.            -.1        (L1/30.)        0.)  ;
PKfu4 = (PROG    -.1            0.           0.        (-2.*L1/15.)     .1           0.            0.         (L1/30.)) ;
PKfu5 = (PROG (6./(5.*L1))      0.           0.             .1      (-6./(5.*L1))    0.            0.           1.1) ;
PKfu6 = (PROG    0.         (6./(5.*L1))     -.1           0.          0.         (-6./(5.*L1))   -1.1          0.)  ;
PKfu7 = (PROG    0.              .1       (L1/30.)         0.          0.            -.1        (-2.*L1/15.)    0.)  ;
PKfu8 = (PROG    -.1            0.           0.        (L1/30.)         .1           0.            0.        (-2.*L1/15.)) ;
PKfu   = PKfu1 ET PKfu2 ET PKfu3 ET PKfu4 ET PKfu5 ET PKfu6 ET PKfu7 ET PKfu8 ;
Kfu    = MANU 'RIGI' 'TYPE' 'RIGIDITE' 
         Tube (MOTS 'UY' 'UZ' 'RY' 'RZ') 'QUEL' PKfu ;
 
 
list 'RESUME' Afu;
list 'RESUME' Kfu;
 
 
*****************************************
* Mise en forme du tableau de resultats *
*****************************************
 
Mod = TABLE ;
Mod.Reel = TABLE ;
Mod.Imag = TABLE ;
Repeter Bloc6 (2*Nbmod) ;
  i = &Bloc6 ;
  Mod.Reel.i = PROG ;
  Mod.Imag.i = PROG ;
Fin Bloc6 ;
 
*****************************************
*        Projection des matrices        *
*****************************************
 
*--- Struture Mt,Ksp
Mt  = PJBA MStru Basr ;
Ksp = PJBA Kstru Basr ;
*--- Amortissement modale Ct
Ct = AMOR Basr L_AMOR ;
 
*--- Matrice amortissement due a l ecoulement 
*    a multiplier par 2MV pour Afp et MV² pour Kfp
Afp = PJBA Afu Basr ;
Kfp = PJBA Kfu Basr  ;
 
NbModC = 2*Nbmod;
mocomp = MOTS 'ALFA';
numodC = LECT 1 'PAS' 1 NbModC;
 
* ===================================== *
*      Boucle sur la vitesse axiale     *
* ===================================== *
REPE b_Vax (DIME LVa) ;
  k = &b_Vax ;
  Va= Extr LVa k ;
  Vr= Extr LVR k ;
 
  Mess ' ' ;
  Mess '#' k  ': Vitesse axiale =' Va ' Vitesse reduite =' Vr ; 
  Mess '-------------------------------------------------------------------------------' ;
  Mess '  Mode   fn_d       Amort.          P. Reeln1        P. Imag.  Stabilite  ' ;
 
* Matrices fluides totale : Amortissement et Raideur
  At = Ct  ET ((2.*Me*Va)*Afp) ;
  Kt = Ksp ET (Me*Va*Va*Kfp);
 
* Calcul des modes complexes 
  Basc = VIBC Mt Kt At ;
  Bascm = Basc . 'MODES' ;
 
*--- Extraction des frequences
  Prf = PROG ;
  Pif = PROG ;
  REPE Bloc2 NbModC ;
    i = &Bloc2 ;
    Modn = Bascm.i ;
    Prf = Prf ET (PROG Modn.'FREQUENCE_REELLE') ;
    Pif = Pif ET (PROG Modn.'FREQUENCE_IMAGINAIRE') ;
  FIN bloc2 ;
 
*--- Tri croissant des frequences
*   NbModC  = DIME Prf ; 
*   Prfc Pifc numodC = ORDO Prf Pif numodC0;
* ou plutot suivi des modes (auteur : Benoit Prabel)
  si (ega k 1);
*   1er pas : on ordonne par fr croissant  
    Prfc Pifc numodCc = ORDO Prf Pif numodC;
    Bascmc  = COPI Bascm;
    repe b0 NbModC; 
        Bascmc . &b0 = Bascm . (extr numodCc &b0);
    fin  b0;
  sinon;
    xcou01 = prog;
*   boucle sur les modes du pas precedent --------------------
    repe b0 NbModC; 
       DEFR0 = Bascm0 . &b0 . 'DEFORMEE_MODALE_REELLE';
       DEFI0 = Bascm0 . &b0 . 'DEFORMEE_MODALE_IMAGINAIRE';
       frqR0 = Bascm0 . &b0 . 'FREQUENCE_REELLE';
       frqI0 = Bascm0 . &b0 . 'FREQUENCE_IMAGINAIRE';
*      boucle sur les modes du pas actuel ----------------    
       repe b1 NbModC; 
         DEFR1 = Bascm . &b1 . 'DEFORMEE_MODALE_REELLE';
         DEFI1 = Bascm . &b1 . 'DEFORMEE_MODALE_IMAGINAIRE';
         frqR1 = Bascm . &b1 . 'FREQUENCE_REELLE';
         frqI1 = Bascm . &b1 . 'FREQUENCE_IMAGINAIRE';
*        ecart en frequence^2
         dfrq1 = ((frqR1 - frqR0)**2) + ((frqI1 - frqI0)**2);
         ref1 =  ((frqR1**2) + (frqI1**2)) + ((frqR0**2) + (frqI0**2));
*        produit scalaire des vecteurs propres
         cos1 = (XTY DEFR1 DEFR0 mocomp mocomp)
              + (XTY DEFI1 DEFI0 mocomp mocomp); 
         sin1 = (XTY DEFI1 DEFR0 mocomp mocomp)
              - (XTY DEFR1 DEFI0 mocomp mocomp); 
         modul1 = (cos1**2) + (sin1**2);
         xsuiv1 = (abs (1. - modul1)) + (0.5*dfrq1 / ref1);
         xcou01 = xcou01 et xsuiv1;
       fin  b1;
*      fin de boucle sur les modes du pas actuel ---------
    fin  b0;
*   fin de boucle sur les modes du pas precedent -------------    
*   SUIVI = apprairage 
*         = minimisation du cout (cad des distances) par permutation
    couttot numodCc Prfc Pifc 
    = ORDO numodC Prf Pif 'COUT' xcou01 ;
*   on modifie la base modale
    Bascmc  = COPI Bascm;
    repe b0 NbModC; 
        Bascmc . &b0 = Bascm . (extr numodCc &b0);
    fin  b0;
  finsi;
* pour le prochain pas  
  Bascm0 = Bascmc; 
 
*--- Sauvegarde
  REPE Bloc5 NbModC ;  i = &Bloc5 ;
    f  = EXTR Prfc i ;
    if = EXTR Pifc i ;
 
    Msg = 'STAB.' ;
    SI ((< if -1.D-10) ET (>EG f 0.) ET (f < 1e-5 )) ;
      Msg = 'INST.' ;
      SI (EGA VSauv 'N') ;
        Vcrit = Va_CO - (p_fin/2.) ;
        VSauv = 'O' ;
MESS 'bipbip' ;
      FINSI ;
    FINSI ;
 
    SI (> f 0.) ;
*--- Valeurs propres relatives a lambda = i*omega = i*2*PI*f
      MESS  i  f  (if/f)  (-2*PI*if)  (2*PI*f)  Msg  (((f**2) + (if**2))**0.5)  (if/(((f**2) + (if**2))**0.5));
    SINON ;
      SI (EGA f 0.) ;
        MESS i f '     ------         ' '       ' (-2*PI*if) (2*PI*f)  Msg (((f**2) + (if**2))**0.5)  (if/(((f**2) + (if**2))**0.5)) ;
      FINSI ;
    FINSI ;
    Mod.Reel.i = Mod.Reel.i ET  f ;
    Mod.Imag.i = Mod.Imag.i ET  if ;
 
  FIN Bloc5 ;
 
  SI GRAPH ;
    MESS '-------------------------------------------------------------------------------' ;
    MESS ' ' ;
  FINSI ;
FIN b_Vax ;
*
SAUTER 2 LIGNES ;
 
**********************************
*           Traces               *
**********************************
* 4*2 couleurs
  numero = LECT Nbmod 'PAS' -1 1 1 'PAS' 1 Nbmod;
  coco4 = MOTS 'ROUG' 'ORAN' 'TURQ' 'BLEU' ;  
  Nbcour = 2*Nbmod ;
  Tab1 = TABLE ;
  Tab2 = TABLE ;
  REPE Ident1 Nbcour ;
*     Tab1 . &Ident1 = 'NOLI MARQ S PLEIN TRIU'  ;
    Tab1 . &Ident1 = 'NOLI MARQ S PLUS'  ;
    Tab2 . &Ident1 = 'MARQ S PLUS'  ;
  FIN Ident1 ;
 
*--- Traces relatifs a la frequence
  Evol1 = VIDE 'EVOLUTIO';
  Evol2 = VIDE 'EVOLUTIO';
*--- Traces relatifs a OMEGA 
  Evol3 = VIDE 'EVOLUTIO';
  Evol4 = VIDE 'EVOLUTIO';
  Evol5 = VIDE 'EVOLUTIO';
* --- Boucle sur les modes
  REPE Trac1 Nbcour ;
    k = &Trac1 ;
    numk = EXTR numero k; 
    si (k &lt;EG Nbmod);  legek = CHAI 'Mode' ' ' numk;
    sinon;             legek = MOT 'PAS DE LEGENDE';
    finsi;
    cocok = EXTR coco4 numk;
    SI (> (DIME Mod.Reel.k) 0) ;
      Evol1 = Evol1 ET (EVOL cocok MANU 'LEGE' legek 
      'V_{r}' LVR 'Re(f) (Hz)' (Mod.Reel.k)) ;
      Evol2 = Evol2 ET (EVOL cocok MANU  'LEGE' legek 
      'V_{r}' LVR 'Im(f) (/s)' (Mod.Imag.k)) ;
      Evol3 = Evol3 ET (EVOL cocok MANU  'LEGE' legek 
      'V_{r}' LVR 'Re(w)' (2.*PI*(Rw*Mod.Reel.k))) ;
      Evol4 = Evol4 ET (EVOL cocok MANU  'LEGE' legek 
      'V_{r}' LVR 'Im(w)' (2.*PI*Rw*(Mod.Imag.k))) ;
      Evol5 = Evol5 ET (EVOL cocok MANU  'LEGE' legek 
      'Re(w)' (2.*PI*Rw*(Mod.Reel.k)) 'Im(w)' (2.*PI*Rw*(Mod.Imag.k))) ;
    SINON;
      MESS 'Mode' k 'vide !';
    FINSI ;
  FIN Trac1 ;
 
SI GRAPH;
 
* avec legende
  DESS Evol1 'GRILL' 'GRIS'  Tab2 'LEGE';
* sans/avec zoom
  TITR 'Vitesse reduite - Frequences (partie reelle)' ;
  DESS Evol1 'GRILL' 'GRIS'  Tab1 ;
  DESS Evol1 'XBOR'  0. 12. XGRA 1. 'YBOR'   0.  120. YGRA 10 'GRILL' 'GRIS' Tab1 ;
 
  TITR 'Vitesse reduite - Frequences (partie imaginaire)' ;
  DESS Evol2 'GRILL' 'GRIS'  Tab1 ;
  DESS Evol2 'XBOR'  0. 12. XGRA 1. 'YBOR'  -0.   10. YGRA 1. 'GRILL' 'GRIS'  Tab1 ;
 
  TITR 'Frequences (partie reelle - partie imaginaire)' ;
  DESS Evol5 'GRILL' 'GRIS'  Tab1 ;
  DESS Evol5 'XBOR'  0. 70. XGRA 10. 'YBOR'   0.   20. YGRA 2. 'GRILL' 'GRIS'  Tab1 ;
FINSI ;
 
 
 
*** TEST A PREVOIR ***
 
 
 
FIN ; 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales