Télécharger dy_dev12.dgibi

Retour à la liste

* fichier :  dy_dev12.dgibi
************************************************************************      
*    
*    Cas-Test de la liaison palier de l'operateur DYNE       
*                                                                       
*                                                                       
* Masse ponctuelle sur 1 squeeze-film (NL)
*   Chargement statique = poids propre = -m*g  -->  Wcharg
* + Chargement tournant = balourd = m*u*w^2    -->  Wbal 
* avec ou sans centering spring de raideur k
*                                                                       
*
* ref : JI Inayat-Hussain et al, On the bifurcations of a rigid rotor
*       response in squeeze-film dampers, Journal of Fluids and 
*       Structures 17 (2003) 433-459
*                                                                       
*
* Auteur : BP, 2015-05-05 
*                                                                       
*
* Mots-clés : Vibrations non-lineaires, DYNE, squeeze-film,      
*             dynamique des machines tournantes                 
*                                                                       
************************************************************************      
*
OPTI 'DIME' 3 'ELEM' SEG2 ;
COMPLET = FAUX;
* COMPLET = VRAI;
* GRAPH   = FAUX;
GRAPH   = VRAI; 
 
* FLCOUR=VRAI pour utiliser le modele de palier_court (sinon rodeli)
* FLCOUR = FAUX;
FLCOUR = VRAI;
 
* devogelaere ou diff centrees ?
FLDEVO = vrai; 
FLDEVO = faux; 
 
* donnees du squeeze :
Lpalier = 0.05;
Rarbre  = 0.02;
* mu = 0.02;
rhofl   = 900.;
xjeu    = 0.50E-3;   xjeum1 = 1./xjeu;
* adimensionne
B = 0.015; valB = mot'0015';
* B = 0.09; valB = mot'009';
* B = 0.40; valB = mot'040';
 
* donnees du rotor (rigide) :
M1 = 1.;
 
*amplitude force statique
g = 10.;
* poids = M1*g;
 
* poids propre adimensionne
W = 0.05; valW = mot '005';
* W = 0.1; valW = mot '010';
* W = 0.5; valW = mot '050';
 
* on deduit la vitesse de rotation
si (EGA W 0.); 
  poids = 0.;
  OMEGHZ = 100. ;
  OMEGA = OMEGHZ * (2*pi);
sinon;
  poids = M1*g;
  OMEGA  = (g / (xjeu * W))**0.5;
  OMEGHZ = OMEGA/(2*pi);
mess '\W = ' OMEGA 'rad/s = ' OMEGHZ 'Hz' ;
 
* on deduit la viscosite
mu = B * M1 * OMEGA * (xjeu**3) / (Rarbre * (Lpalier**3));
mess '\m = ' mu ' unite SI ' ;
 
* amplitude balourd
* Fbalad = 0.05; valU = mot'005';
 Fbalad = 0.10; valU = mot'010';
 Fbalad = 0.20; valU = mot'020';
 Fbalad = 0.30; valU = mot'030';
 Fbalad = 0.355; valU = mot'0355';
*  Fbalad = 0.363; valU = mot'0363';
*  Fbalad = 0.3635; valU = mot'03635';
* *  Fbalad = 0.364; valU = mot'0364';
* *  Fbalad = 0.3641; valU = mot'03641';
* *  Fbalad = 0.370; valU = mot'0370';
*  Fbalad = 0.385; valU = mot'0385';
* * Fbalad = 0.40; valU = mot'040';
* Fbal = Fbalad * M1 * xjeu * (OMEGA**2);
Fbal = Fbalad * M1 * xjeu * (OMEGA**2);
mess 'Fbal='Fbal 'N';
 
* centering spring
* S = 0.10; valS = mot '010'; 
S = 0.30; valS = mot '030'; 
* S = 0.50; valS = mot '050'; 
k1 = M1 * ((S * OMEGA)**2);
mess 'k='k1 'N/m';
 
Y0 = W / (S**2);
mess 'static eccentricity ratio Y0='Y0;
 
 
* valeurs adimensionnées du squeeze :
cha1 = chai 'Squeeze-film : B=' FORMAT '(F5.3)' B 
                        ' - W=' FORMAT '(F5.3)' W  
                        ' - S=' FORMAT '(F5.3)' S 
                        ' - U=' FORMAT '(F5.3)' Fbalad ;
mess cha1;
TITR cha1;
 
* prefixe
si FLCOUR; mopal = mot 'court';
sinon;     mopal = mot 'rodeli';
finsi;
si FLDEVO; moalgo = mot 'DV';
sinon;     moalgo = mot 'DC';
finsi;
 
prefix = chai 'dy_dev12-'mopal'-'moalgo
              '-B'valB'-W'valW'-S'valS'-U'valU;
iCI = 0;  prefix = chai prefix '-CI0';
* iCI = 1;  prefix = chai prefix '-CI1';
* iCI = 2;  prefix = chai prefix '-CI2';
* iCI = 3;  prefix = chai prefix '-CI3';
mess prefix ;
 
ficxdr = chai prefix'.xdr';
ficps = chai prefix'.ps';
OPTI 'TRAC' 'PSC' 'POTR' 'HELVETICA_16' 'EPTR' 6 'FTRA' ficps;
 
 
************************************************
* Nombre de pas de temps, pas de temps
************************************************
 
si (COMPLET);  
  si (((Fbalad > 0.36 ) et (Fbalad <eg 0.37)) 
       ou (Fbalad ega 0.385)); Nperiod = 1000 ; 
  sinon;                       Nperiod =  400 ;
  finsi;
 Nperiod =  400 ;
sinon;          Nperiod =  80 ; 
finsi;
* nombre de points par periode, nombre de period
nepp1 = 8;
npp1 = 2**nepp1;
NPas1 = npp1 * Nperiod;
dt1 = 1. / (OMEGHZ * npp1);
Tfinal = NPas1 * dt1;
* Sortie tous les NSor1 pas de temps
* NSor1 = 2**3 ;
NSor1 = 1 ;
NbMod1 = 2;
 
 
************************************************
* Creation des MODES
************************************************
 
P1 = 0.  0. 0. ;
Mail1 = MANU 'POI1' P1;
*
moprim = mots 'UY' 'UZ';
MMass1= MANU 'RIGIDITE' 'TYPE' 'MASSE' Mail1 moprim (prog M1 0 0 M1);
MRigi1= MANU 'RIGIDITE' 'TYPE' 'RIGIDITE' Mail1 moprim (prog k1 0 0 k1);
*
chpuy = manu 'CHPO' Mail1 2 'UY' 1. 'UZ' 0. 'NATURE' 'DIFFUS';
chpuz = manu 'CHPO' Mail1 2 'UY' 0. 'UZ' 1. 'NATURE' 'DIFFUS';
lchp = SUIT chpuy chpuz;
 
TBasR1 = VIBR 'PROC' (prog -0.1)  (lect NbMod1) 
         MRigi1 MMass1  lchp;
 
mess 'mode 1 : f=' TBasR1 . 'MODES' . 1 . 'FREQUENCE'
              'M=' TBasR1 . 'MODES' . 1 . 'MASSE_GENERALISEE';
mess 'mode 2 : f=' TBasR1 . 'MODES' . 2 . 'FREQUENCE'
              'M*=' TBasR1 . 'MODES' . 2 . 'MASSE_GENERALISEE';
 
*
************************************************
* Creation de la table des CONDITIONS INITIALES
************************************************
 
TInit1 = TABLE 'INITIAL' ;
*
si (ega iCI 0); UY0 = 0.       ; UZ0 =  0.       ; finsi;
* * 1er jeu de CI :
si (ega iCI 1); UY0 = 0.1*xjeu ; UZ0 = -0.1*xjeu ; finsi;
* 2eme jeu de CI :
si (ega iCI 2); UY0 = 0.0*xjeu ; UZ0 =  0.5*xjeu ; finsi;
* 3eme jeu de CI :
si (ega iCI 3); UY0 = 0.0*xjeu ; UZ0 = -0.8*xjeu ; finsi;
* si (ega iCI 3); UY0 = 0.1*xjeu ; UZ0 = -0.95*xjeu ; finsi;
 
Dep0 = MANU 'CHPO' Mail1 2 'UY' UY0 'UZ' UZ0 NATURE DIFFUS ;
*
REPE bb1 NbMod1 ;
 i1 = &bb1 ;
 Def1 = TBasR1 . 'MODES' . i1 . 'DEFORMEE_MODALE' ;
 M_G1 = TBasR1 . 'MODES' . i1 . 'MASSE_GENERALISEE' ;
 P_R1 = TBasR1 . 'MODES' . i1 . 'POINT_REPERE' ;
 qi1 = (YTMX Def1 MMass1 Dep0) / M_G1 ;
 f_r1 = TBasR1 . 'MODES' . i1 . 'FREQUENCE';
 mess 'Mode' i1 'f=' f_r1 ' M*=' qi1;
 SI (EGA i1 1) ;
  TInit1.'DEPLACEMENT' = MANU 'CHPO' P_R1 1 'ALFA' qi1 NATURE DIFFUS ;
  TInit1.'VITESSE' = MANU 'CHPO' P_R1 1 'ALFA' 0. NATURE DIFFUS ;
 SINON ;
  TInit1.'DEPLACEMENT' = TInit1.'DEPLACEMENT' ET 
                      (MANU 'CHPO' P_R1 1 'ALFA' qi1 NATURE DIFFUS) ;
  TInit1.'VITESSE' = TInit1.'VITESSE' ET
                  (MANU 'CHPO' P_R1 1 'ALFA' 0. NATURE DIFFUS) ;
 FINSI ;
FIN bb1 ;
 
*
**********************************************
* Creation de la table des LIAISONS (paliers)
**********************************************
*
TLiai1 = TABLE 'LIAISON' ;
TLiai1 . LIAISON_B = TABLE 'LIAISON_B' ;
PAL1 = TABLE 'LIAISON_ELEMENTAIRE' ;
TLiai1 . LIAISON_B . 1 = PAL1 ;
PAL1 . TYPE_LIAISON     = 'PALIER_FLUIDE' ;
PAL1 . POINT_SUPPORT    = P1 ;
PAL1 . LONGUEUR_PALIER  = Lpalier ;
PAL1 . RAYON_ARBRE      = Rarbre ;
PAL1 . VISCOSITE_FLUIDE = mu ;
PAL1 . RHO_FLUIDE       = rhofl ;
PAL1 . PRESSION_ADMISSION = 0. ;
PAL1 . VITESSE_ROTATION = 0.;
si (FLCOUR);
  PAL1 . MODELE_PALIER = 'PALIER_COURT' ;
  PAL1 . JEU_USINAGE = xjeu ;
sinon;
  PAL1 . MODELE_PALIER = 'RODELI' ;
  PAL1 . AFFICHAGE = 0 ;
  PAL1 . GEOMETRIE_PALIER = TABLE ;
  PAL1 . GEOMETRIE_PALIER . NOMBRE_LOBES = 1 ;
  PAL1 . GEOMETRIE_PALIER . CRITERE_ARRET = 1e-6 ;
  PAL1 . GEOMETRIE_PALIER . 1 = TABLE ;
  PAL1 . GEOMETRIE_PALIER . 1 . ASYMETRIE = 0. ;
  PAL1 . GEOMETRIE_PALIER . 1 . JEU_USINAGE = xjeu ;
  PAL1 . GEOMETRIE_PALIER . 1 . PRECHARGE = 0. ;
  PAL1 . GEOMETRIE_PALIER . 1 . ANGLE_DEBUT = 0. ;
  PAL1 . GEOMETRIE_PALIER . 1 . AMPL_ANGLE = 2 * pi ;
  PAL1 . GEOMETRIE_PALIER . 1 . COEF_SUR = 1.715 ;
  PAL1 . GEOMETRIE_PALIER . 1 . NB_MAILLES = 120 ;
*   PAL1 . GEOMETRIE_PALIER . 1 . NB_MAILLES = 360 ;
finsi;
 
 
 
***********************************************
* Creation de la table des VARIABLES DE SORTIE
***********************************************
TSort1 = TABLE 'SORTIE' ;
TSort1 . 'VARIABLE' = TABLE 'VARIABLE' ;
TSort1 . 'LIAISON_B' = TABLE 'LIAISON_B' ;
TSort1 . 'LIAISON_B' . PAL1 = VRAI ;
 
 
*************************************************************
* Creation de la table du CHARGEMENT
*************************************************************
 
LTps1 = PROG 0.0 'PAS' dt1 'NPAS' (NPas1+1) ;
 
*** evol temporelle ***
 
*gravite
LFW1 = PROG  (NPas1+2)*poids;
evW1 = evol ROSE MANU LTps1 LFW1;
Chp_FZ1 = MANU CHPO 2 Mail1 'FY' 0. 'FZ' -1. ;
 
 
*balourd 
LFdY1 = PROG 'SINU' OMEGHZ 'AMPL' Fbal LTps1 ;
LFdZ1 = PROG 'SINU' OMEGHZ 'PHAS' 90. 'AMPL' Fbal LTps1 ;
Ev_FdY1 = EVOL 'BLEU' MANU 'LEGE' 'FY' 't(s)' LTps1 'FY' LFdY1 ;
Ev_FdZ1 = EVOL 'OCEA' MANU 'LEGE' 'FZ' 't(s)' LTps1 'FZ' LFdZ1 ;
 
dess ((-1.*evW1) et Ev_FdY1 et (-1.*Ev_FdZ1))  
   'POSX' 'CENT' 'TITX' 't(s)' 'XBOR' 0. (1./OMEGHZ)
   'POSY' 'CENT' 'TITY' 'F(N)' 'LEGE' 'NO';
dess ((-1.*evW1) et Ev_FdY1 et (-1.*Ev_FdZ1))  
   'POSX' 'CENT' 'TITX' 't(s)' 'XBOR' 0. Tfinal
   'POSY' 'CENT' 'TITY' 'F(N)' 'LEGE' 'NO';
 
*** chpoints ***
Chp_FdY1 = MANU CHPO 2 Mail1 'FY' 1. 'FZ' 0. ;
Chp_FdZ1 = MANU CHPO 2 Mail1 'FY' 0. 'FZ' -1. ;
 
*** chargement *** 
Chatot = (CHAR Chp_FdY1 Ev_FdY1) 
      ET (CHAR Chp_FdZ1 Ev_FdZ1) 
      ET (CHAR Chp_FZ1  evW1) ;
 
TChar1 = TABLE 'CHARGEMENT' ;
TChar1 . 'BASE_A' = PJBA TBasR1 Chatot ;
 
 
****************************************************************
* Cacul temporel 1 : Appel a DYNE
****************************************************************
si FLDEVO;  algodyne = mot 'DE_VOGELAERE';
sinon;      algodyne = mot 'DIFFERENCES_CENTREES';
finsi; 
 
TResu1 = DYNE algodyne TBasR1 TChar1 TInit1 TSort1
              TLiai1 dt1 NPas1 NSor1 ;
 
uy1 = TResu1 . PAL1 . 'UY_ARBRE';
uz1 = TResu1 . PAL1 . 'UZ_ARBRE';
vy1 = TResu1 . PAL1 . 'VY_ARBRE';
vz1 = TResu1 . PAL1 . 'VZ_ARBRE';
fy1 = TResu1 . PAL1 . 'FY_PALIER';
fz1 = TResu1 . PAL1 . 'FZ_PALIER';
tprog = TResu1 . 'TEMPS_DE_SORTIE' ;
nt = dime tprog;
 
 
********************************
* graphiques
********************************
si (GRAPH);
 
  evuy1 = EVOL BLEU 'MANU' 'LEGE' 'UY/e' 't' tprog 'UY/e' (uy1*xjeum1);
  evuz1 = EVOL TURQ 'MANU' 'LEGE' 'UZ/e' 't' tprog 'UZ/e' (uz1*xjeum1);
  evvy1 = EVOL BLEU 'MANU' 'LEGE' 'VY' 't' tprog 'VY' (vy1);
  evvz1 = EVOL TURQ 'MANU' 'LEGE' 'VZ' 't' tprog 'VZ' (vz1);
  evfy1 = EVOL BLEU 'MANU' 't' tprog 'FY' fy1;
  evfz1 = EVOL TURQ 'MANU' 't' tprog 'FZ' fz1;
  evyz1 = EVOL BLEU 'MANU' 'UY/e' (uy1*xjeum1) 'UZ/e' (uz1*xjeum1);
  t1 = prog -180. PAS 1. 180.;
  y1 = sin t1; z1 = -1.* (cos t1);
  evjeu = evol 'MANU' y1 z1;
  Tjeu = tabl; Tjeu . 2 = mot 'TIRR';
 
  si (non complet);
    dess (evuy1 et evuz1) 
    'XBOR' 0. Tfinal TITX 'U/e' POSX CENT POSY CENT LEGE NO;
    dess (evvy1 et evvz1) 
    'XBOR' 0. Tfinal POSX CENT POSY CENT LEGE NO;
    dess (evfy1 et evfz1) 
    'XBOR' 0. Tfinal POSX CENT POSY CENT LEGE NO;
  finsi;
    dess (evyz1 et evjeu)
      POSX CENT POSY CENT 'CARR' Tjeu;
 
* derniere periode
   si complet;  neplast = 7 ; 
   sinon;       neplast = 4 ; 
   finsi;
  nplast = 2**neplast ; 
  nlast = (Nperiod - nplast) * npp1 / NSor1;
  llast = lect nlast pas 1 (NPas1 / NSor1);
  tprogf = extr tprog llast;
  tprogf = tprogf - (extr tprogf 1);
  uy1f = extr uy1 llast;
  uz1f = extr uz1 llast;
  evuy1f = EVOL BLEU 'MANU' 'LEGE' 'UY/e' 
           't' tprogf 'UY/e' (uy1f*xjeum1);
  evuz1f = EVOL TURQ 'MANU' 'LEGE' 'UZ/e'
          't' tprogf 'UZ/e' (uz1f*xjeum1);
  evyz1f = EVOL BLEU 'MANU' 'UY/e' (uy1f*xjeum1) 'UZ/e' (uz1f*xjeum1);
  dess (evyz1f et evjeu)
    POSX CENT POSY CENT 'CARR' Tjeu;
  dess (evyz1f)
    POSX CENT POSY CENT 'CARR' Tjeu;
 
* DSP   
  ndsp = nepp1 + neplast ;
  dspy1 = DSPR ndsp evuy1f 'BLEU';
  dspz1 = DSPR ndsp evuz1f 'TURQ';
* adimensionnement des frequences (+ pratique)  
  dspy1 = EVOL 'BLEU' 'MANU' 'LEGE' 'UY/e' 
               '\w/\W' ((EXTR dspy1 'ABSC')/OMEGHZ) 
               'UY/e' (EXTR dspy1 'ORDO');
  dspz1 = EVOL 'TURQ' 'MANU' 'LEGE' 'UZ/e' 
               '\w/\W' ((EXTR dspz1 'ABSC')/OMEGHZ) 
               'UZ/e' (EXTR dspz1 'ORDO');
  si (non complet);
    dess (dspy1 et dspz1) 
    POSX CENT POSY CENT              LOGY YBOR 1.E-15 1. LEGE NE ;
  finsi;
  dess (dspy1 et dspz1) 
  POSX CENT POSY CENT XBOR 0. 20.  LOGY YBOR 1.E-15 1. LEGE NE ;
  dess (dspy1 et dspz1) 
  POSX CENT POSY CENT XBOR 0. 4. XGRA 0.5 LOGY YBOR 1.E-15 1. LEGE NE;
  ymax  = maxi (EXTR dspy1 'ORDO'); 
*   nymax = enti ((log ymax) / (log 10.)) + 1;
  nymax = enti 'SUPERIEUR' ((log ymax) / (log 10.));
  nymin = nymax - 6; 
  ymin = 10.**nymin;
  ymax = 10.**nymax;
  dess (dspy1 et dspz1) 
  POSX CENT POSY CENT XBOR 0. 4. XGRA 0.5 LOGY YBOR ymin ymax LEGE NE;
  dess (dspy1 et dspz1) 
  POSX CENT POSY CENT XBOR 0. 2. XGRA 0.25 LOGY YBOR ymin ymax LEGE NE;
 
finsi;
 
 
* *excentricite et angle de calage
* prex1 = ((uy1**2) + (uz1**2))**0.5;
* mas00 =  (masq (abs uy1) 'EGINFE' 1.E-30 ) 
*       *  (masq (abs uz1) 'EGINFE' 1.E-30 );
* uy100 = (((prog nt*1.) - mas00) * uy1)  + mas00;
* prang1 = ATG uy100 (-1.*uz1);
* mess 'excentricite ='; list prex1;
* mess 'angle de calage ='; list prang1;
* tfin = extr tprog nt;
* an1fin = extr prang1 nt;
* ex1fin = extr prex1 nt;
* mess 'excentricite relative (t=' tfin ')=' (ex1fin*xjeum1);
*  
 
 
****************************************************************
* SAUVEGARDE DES COURBES
****************************************************************
si (FAUX);
 
TEVOL1 = INDEX '*EVOLUTIO';
list  TEVOL1;
opti sauv ficxdr;
sauv TEVOL1;
 
finsi ;
 
 
****************************************************************
* Tests de Non Regression
****************************************************************
 
* on verifie que l'on obtient un regime 4T-Periodique
 
* recup des 4 dernieres periodes
Tuy = tabl;
Tuz = tabl;
nfin = NPas1 / NSor1;
repe b4 5;
  nlast = (Nperiod - &b4) * npp1 / NSor1;
  llast = lect nlast pas 1 nfin;
  Tuy .  &b4 = extr uy1 llast;
  Tuz .  &b4 = extr uz1 llast;
  mess 'nlast - nfin' nlast nfin (dime Tuy .  &b4);
  nfin = nlast;
fin  b4;
 
* ecart uy(t)-uy(t-1*T)
yecart = prog;
zecart = prog;
repe b4 4 ;
  yecart = yecart et (MAXI (Tuy . 1 - Tuy . (&b4 + 1)) 'ABS');
  zecart = zecart et (MAXI (Tuz . 1 - Tuz . (&b4 + 1)) 'ABS');
fin  b4;
yecart = yecart*xjeum1;
zecart = zecart*xjeum1;
* resultats obtenus par BP en 2016-07-07 (avec COMPLET=vrai):
* 0.84326      0.45235      0.71956      5.85832E-11
* 1.3414      0.24801       1.2888      5.62050E-11
* resultats obtenus par BP en 2016-07-07 (avec COMPLET=faux):
* 0.84326      0.45235      0.71956      3.73157E-08
* 1.3414       0.24801       1.2888      2.20570E-08
list yecart;
list zecart;
opti echo 0;
saut lign ;
mess '-------------------------------------------------------------';
repe b4 4 ;
  mess ' Max_t{ y(t) - y(t-'&b4'*T)} = ' (extr yecart &b4);
fin  b4;
mess '-------------------------------------------------------------';
opti echo 1;
 
SI (    ((extr yecart 1) < 0.1) ou ((extr yecart 2) < 0.1)
     ou ((extr yecart 3) < 0.1) ou ((extr yecart 4) > 1.E-6)
     ou ((extr zecart 1) < 0.1) ou ((extr zecart 2) < 0.1)
     ou ((extr zecart 3) < 0.1) ou ((extr zecart 4) > 1.E-6) );
 ERRE 5;
SINON;
 ERRE 0;
FINSI;
 
 
fin ; 
 
 
 
 
 
 
 
 
 

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