Télécharger la plaquette Cast3M


fr en

Télécharger dyne06.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : dyne06.dgibi
  2. *
  3. ************************************************************************
  4. *
  5. * DYNE06 : 
  6. *
  7. * Calcul d'un rotor de type Jeffcott avec contact frottant
  8. *
  9. *                 ___                           ^ Z
  10. *                 _|_ kc,mu                     |__
  11. *               _______  jeu                    |_ \_ 
  12. *   |   k,c    |       |                        | \  \
  13. *   |----------|   m   |----------|             +--|--|-----> Y
  14. *   |          |_______|                           R  R+jeu
  15. *                 ___    jeu                       
  16. *                 _|_ kc,mu
  17. *                
  18. *  ^ Z
  19. *  |
  20. *  +----> X
  21. *  
  22. * Test de l'option VITESSE_ENTRAINEMENT pour la liaison DYNE
  23. * POINT_CERCLE_FROTTEMENT qui ajoute une sur-vitesse tangentielle
  24. *
  25. * Ref : [Xie et al., MSSP, 2015]
  26. *
  27. * Auteur : BP, 2020-03-25
  28. *
  29. ************************************************************************
  30.  
  31.  
  32. ************************************************************************
  33. *                      OPTIONS et DONNEES
  34. ************************************************************************
  35.  
  36. * options 
  37.   OPTI DIME 3 ELEM SEG2;
  38.   GRAPH   = FAUX;
  39.   COMPLET = FAUX;
  40.   OPTI TRAC PSC EPTR 8 POTR HELVETICA_16;
  41.   prefix = chai 'DYNE06';
  42.   TITR prefix;
  43.  
  44. * rotor data
  45.   m = 1.  ;
  46.   k = 100.;
  47.   c = 5.  ;
  48.   jeu  = 0.105;
  49.   R = 20.*jeu ;
  50.   ebal = 0.1;
  51. *   Fbal = ebal * m * (Omega**2) ; avec Omega = vitesse de rotation en rad/s
  52.  
  53. * contact parameters
  54.   kchoc= 25. * k;
  55.   mu_s = 0.2;
  56.   mu_d = 0.2;
  57.   si (mu_s < mu_d); erreur 21; finsi;
  58.   Kt = 10.*kchoc;
  59.   Ct = 0.5*(m*Kt)**0.5;
  60.  
  61. * deduced parameter
  62.   w   = (k/m)**0.5;                 
  63.   wHz = w / (2.*pi);
  64.   wchoc = (kchoc/m)**0.5;           
  65.   xi  = 0.5 * c / ((kchoc*m)**0.5); 
  66.   mess 'w=' w '(' wHz ')Hz < wchoc=' wchoc; 
  67.   mess 'xi=' xi;
  68.  
  69. * critical timestep (Diff Centrees, no damping)
  70.   wadhe = ((k+Kt)/m)**0.5;
  71.   dtc = 2. / wadhe;
  72.  
  73.  
  74. ************************************************************************
  75. *                       CALCUL ANALYTIQUE LINEAIRE 
  76. ************************************************************************
  77. * si pas de contact, alors pas de couplage entre y et z
  78. * on resout pour y par exemple :
  79. *
  80. *     m y" + c y' + ky = meW^2 cos Wt
  81. *
  82. * on pose y(t) =  A  cos Wt + B  sin Wt
  83. *         y'   = -Aw sin Wt + Bw cos Wt
  84. *         y"   = -Aw^2 cos Wt - Bw^2 sin Wt
  85. *
  86. * [ k-W^2m    Wc   ] (A) = (meW^2)      termes en cos
  87. * [  -Wc    k-W^2m ] (B) = (  0  )      termes en sin
  88. *
  89. * 2eme ligne : B = Wc/k-W^2m * A
  90. * 1ere ligne : A = meW^2 / [ k-W^2m + W^2c^2/k-W^2m]
  91.  
  92. Wadim = prog 0.02 pas 0.01 0.19 0.198 0.202 0.21 pas 0.01 0.5 pas 0.1 1.2;
  93. Wanal = wchoc * Wadim;
  94. Wanal2 = Wanal**2;
  95. * m=1 -> ca simplifie
  96. aux1  = k - Wanal2;
  97. aux2  = (c**2) * Wanal2;
  98. denom1 = aux1 + (aux2/aux1);
  99. A_anal = (ebal * Wanal2) / denom1;
  100. B_anal = (c*Wanal) / aux1 * A_anal;
  101. * evA = evol VIOL manu Wanal A_anal;
  102. * evB = evol ROSE manu Wanal B_anal;
  103. * dess (evA et evB);
  104. Amp_anal = ( (A_anal**2) + (B_anal**2) )**0.5;
  105. evAmp = evol VIOL manu  Wanal Amp_anal;
  106. evjeu = evol manu Wanal (prog (dime Wanal)*jeu);
  107. * dess (evjeu et evAmp);
  108. * en adimensionne
  109. evjeu2 = evol manu Wadim (prog (dime Wanal)*1);
  110. evAmp2 = evol VIOL manu 'w/w_{c}' Wadim 'r/j' (Amp_anal/jeu);
  111. si GRAPH; dess (evjeu2 et evAmp2) TITRE 'solution analytique lineaire'; finsi;
  112. * opti donn 5 trac X;
  113.  
  114.  
  115. ************************************************************************
  116. *                       BASE MODALE 
  117. ************************************************************************
  118.  
  119. * point physique
  120.   p1 = 0. 0. 0.;
  121.  
  122. * base modale = 2 modes
  123.   palfa1 = 0. 0. 0.;
  124.   phi1 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 1. 'UZ' 0. 'NATURE' 'DIFFUS';
  125.   palfa2 = 0. 0. 0.;
  126.   phi2 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 0. 'UZ' 1. 'NATURE' 'DIFFUS';
  127.  
  128.   TBAS = TABL 'BASE_MODALE';
  129.   TBAS . 'MODES' = TABL 'BASE_DE_MODES'; 
  130.   TBAS . 'MODES' . 1 = TABL 'MODES';
  131.   TBAS . 'MODES' . 1 . 'POINT_REPERE'      = palfa1;
  132.   TBAS . 'MODES' . 1 . 'FREQUENCE'         = wHz;
  133.   TBAS . 'MODES' . 1 . 'MASSE_GENERALISEE' = m;
  134.   TBAS . 'MODES' . 1 . 'DEFORMEE_MODALE'   = phi1;
  135.   TBAS . 'MODES' . 2 = TABL 'MODES';
  136.   TBAS . 'MODES' . 2 . 'POINT_REPERE'      = palfa2;
  137.   TBAS . 'MODES' . 2 . 'FREQUENCE'         = wHz;
  138.   TBAS . 'MODES' . 2 . 'MASSE_GENERALISEE' = m;
  139.   TBAS . 'MODES' . 2 . 'DEFORMEE_MODALE'   = phi2;
  140.  
  141.  
  142. ************************************************************************
  143. *                       AMORTISSEMENT 
  144. ************************************************************************
  145.  
  146.   TAMOR = TABL 'AMORTISSEMENT';
  147.   Ctot = MANU 'RIGIDITE' 'TYPE' 'AMORTISSEMENT' 
  148.         (palfa1 et palfa2) (MOTS 'ALFA') (prog c);
  149.   TAMOR . 'AMORTISSEMENT' = Ctot  ;
  150.  
  151.  
  152. ************************************************************************
  153. *                           LIAISONS 
  154. ************************************************************************
  155.  
  156. * contact-frottant point_cercle_frottement
  157.   TL1 = TABL 'LIAISON_ELEMENTAIRE';
  158.   TL1 . 'TYPE_LIAISON' = mot 'POINT_CERCLE_FROTTEMENT';
  159.   TL1 . 'SUPPORT'      = p1;
  160.   TL1 . 'NORMALE'      = (1. 0. 0.) ;
  161.   TL1 . 'EXCENTRATION' = (0. 0. 0.) ;
  162.   TL1 . 'RAIDEUR'      = kchoc;
  163.   TL1 . 'RAYON'        = jeu ;
  164.   TL1 . 'COEFFICIENT_GLISSEMENT' = mu_d;
  165.   TL1 . 'COEFFICIENT_ADHERENCE'  = mu_s;
  166.   TL1 . 'RAIDEUR_TANGENTIELLE'    = Kt;
  167.   TL1 . 'AMORTISSEMENT_TANGENTIEL'= Ct;
  168. * stockage  
  169.   TLB  = TABL 'LIAISON_B';
  170.   TLB . 1 = TL1;
  171.   TLIA = TABL 'LIAISON';
  172.   TLIA . 'LIAISON_B' = TLB; 
  173.  
  174. ************************************************************************
  175. *                      CONDITIONS INITIALES
  176. ************************************************************************
  177.  
  178. *   TINIT = TABL 'INITIAL';
  179. *   TINIT . 'DEPLACEMENT' 
  180. *     = (MANU 'CHPO' palfa1 1 'ALFA' ... 'NATURE' 'DISCRET')
  181. *    et (MANU 'CHPO' palfa2 1 'ALFA' ... 'NATURE' 'DISCRET');
  182. *   TINIT . 'VITESSE'     
  183. *     = (MANU 'CHPO' palfa1 1 'ALFA' ... 'NATURE' 'DIFFUS')
  184. *    et (MANU 'CHPO' palfa2 1 'ALFA' ... 'NATURE' 'DIFFUS');
  185. *u=v=0 ici --> inutile  
  186.  
  187. ************************************************************************
  188. *                           SORTIES
  189. ************************************************************************
  190.  
  191.   TSORT = TABLE 'SORTIE'; 
  192. * +on stocke les deplacements modaux 
  193.   TSORV = TABLE 'VARIABLE';
  194.   TSORV . 'TYPE_SORTIE' = MOT 'LISTREEL';
  195.   TSORT . 'VARIABLE' = TSORV;
  196. * + on stocke les variables du contact 
  197.   TSORLB = TABLE 'LIAISON_B';
  198.   TSORLB . TL1 = VRAI;
  199.   TSORT . 'LIAISON_B' = TSORLB;
  200.  
  201. * pour le DESSin des evolutions
  202.   Tdess1 = TABL ;
  203.   Tdess1 . 2 = MOT 'TIRR';
  204.   TDSP1 = TABL;  TDSP1 . 'TITRE' = TABL;
  205.  
  206.  
  207. ************************************************************************
  208. *                      PARAMETRES TEMPORELS
  209. ************************************************************************
  210.  
  211. * NOMBRE DE PERIODES SIMULEES, NOMBRE DE PAS PAR PERIODE
  212. *  NPERIOD = 8;
  213. *debug  NPERIOD = 14;
  214.   NPERIOD = 80;
  215.   NPPP    = 2**10;
  216.  
  217. * NOMBRE DE PAS TOTAL, SORTIE TOUS LES NSORT
  218.   NPAS  = NPERIOD*NPPP;
  219.   NSORT = 8;
  220.  
  221. * NUMERO DES POINTS SORTIS APPARTENANT AUX nlast DERNIERES PERIODES
  222. *  nlast = 2;
  223.   nlast = 20;
  224.   llast = LECT (NPPP / NSORT * (NPERIOD - nlast)) 'PAS' 1 (NPPP / NSORT * NPERIOD + 1);
  225.  
  226.  
  227. *-----------------------------------------------------------------------
  228. *--------------> BOUCLE SUR LA VITESSE DE ROTATION TESTEE  
  229. *-----------------------------------------------------------------------
  230.  
  231. * vitesse de rotation (prOmega en rad/s)
  232. si COMPLET;
  233.   prWadim =   PROG   0.1   0.15   0.154   0.2    0.290   0.3  
  234.                      0.33  0.4   0.5    1.      1.2   ;
  235.   coco_dyne = MOTS 'INDI' 'MARI' 'BLEU' 'AZUR'  'TURQ'  'VERT'  
  236.                     'BOUT' 'OR' 'ORAN' 'ROUG'    'BRUN' ;
  237. sino;
  238.   prWadim =   PROG   0.15    0.2    0.33   0.4    1.     ;
  239.   coco_dyne = MOTS  'MARI'  'AZUR' 'BOUT' 'OR'  'ROUG'   ;
  240. finsi;
  241.   prOmega   = wchoc  * prWadim;
  242.  
  243. * RESULTATS
  244. *   ymax = prog;
  245. *   zmax = prog;
  246.   rmax = prog;
  247.   rmin = prog;
  248.   ev_dsp = VIDE 'EVOLUTIO';
  249.  
  250.   REPE BOMEGA (DIME prOmega);
  251. *debug  &BOMEGA=1;
  252. *debug  &BOMEGA=2;
  253. *debug  &BOMEGA=4;
  254.  
  255.   coco  = EXTR coco_dyne &BOMEGA;
  256.   Omega = EXTR prOmega &BOMEGA;  
  257.   OMEGHZ = Omega / (2.*PI);
  258.   MESS '>>> Calcul #' &BOMEGA ' pour OMEGA='Omega 'rad/s (adim ='(Omega/wchoc) ')';
  259.   FMT = MOT '(F6.3)';
  260.   cha1 = chai '\W=' (enti Omega) ' rad/s - \w/\w_{c}=' FORMAT FMT (Omega/wchoc);
  261.   TITRE cha1;
  262.   TDSP1 . 'TITRE' . &BOMEGA = chai '\w/\w_{c}='FORMAT FMT (Omega/wchoc);
  263.  
  264. * ON COMPLETE :
  265. *-LIAISON
  266. ***  TL1 . 'VITESSE_ENTRAINEMENT'    = Omega;
  267.   TL1 . 'VITESSE_ENTRAINEMENT'    = -1.*R*Omega;
  268. *   
  269. *-PERIODE, PAS DE TEMPS, TEMPS FINAL
  270.   TPERIOD = 1. / OMEGHZ;
  271.   DT      = TPERIOD / NPPP;
  272.   TFINAL  = NPAS*DT;
  273.   si (DT < dtc);  mess 'Pas de temps =' DT '  <  critique=' dtc; 
  274.   sinon;          erre 21;
  275.   finsi;
  276.  
  277.  
  278. ************************************************************************
  279. *                           CHARGEMENT
  280. ************************************************************************
  281.  
  282. * nombre de periodes pour la rampe de chargement
  283. NPrampe = 10;
  284.  
  285. * evolution temporelle du chargement
  286. t_char = PROG (-1.*DT) 'PAS' DT (TFINAL + DT);
  287. rampe  =  BORN (t_char/(NPrampe*TPERIOD)) 'COMPRIS' 0. 1.;
  288. F_Y = COS (360.*OMEGHZ * t_char);
  289. F_Z = SIN (360.*OMEGHZ * t_char);
  290. *debug F_Y = 1.;
  291. *debug F_Z = 0.;
  292. evF_Y = EVOL 'ROUG' MANU t_char (rampe * F_Y);
  293. evF_Z = EVOL 'BRIQ' MANU t_char (rampe * F_Z);
  294. * DESS (evF_Y et evF_Z) TITRE 'Chargement' XBOR 0. TFINAL XGRA TPERIOD;
  295.  
  296. * description "spatiale" du chargement
  297. Fbal = ebal*m * (Omega**2);
  298. Fchpo_Y  = MANU 'CHPO' 1 palfa1 'FALF' Fbal;
  299. Fchpo_Z  = MANU 'CHPO' 1 palfa2 'FALF' Fbal;
  300.  
  301. CHAR_Y = CHAR evF_Y Fchpo_Y ;
  302. CHAR_Z = CHAR evF_Z Fchpo_Z ;
  303.  
  304. * table pour DYNE
  305.   TCHARG = TABLE 'CHARGEMENT' ;
  306.   TCHARG . 'BASE_A' = CHAR_Y et CHAR_Z ;
  307.  
  308.  
  309.  
  310. ************************************************************************
  311. * CALCUL DYNE : INTEGRATION TEMPORELLE
  312. ************************************************************************
  313.  
  314. * TDYNE = DYNE 'DE_VOGELAERE'         TBAS TAMOR TLIA TSORT TCHARG
  315.   TDYNE = DYNE 'DIFFERENCES_CENTREES' TBAS TAMOR TLIA TSORT TCHARG  
  316.                               NPAS DT   NSORT ;
  317.  
  318. ************************************************************************
  319. * POST TRAITEMENT 
  320. ************************************************************************
  321.  
  322. * temps t
  323.   pr_time = TDYNE . 'TEMPS_DE_SORTIE'; 
  324. * deplacement modaux
  325.   pr_y = TDYNE . 'DEPLACEMENT' . palfa1;
  326.   pr_z = TDYNE . 'DEPLACEMENT' . palfa2;
  327.  
  328. * cercle r=R+jeu
  329.   pr_360deg = prog 0. PAS 2.5 360.;
  330.   ev_butee  = EVOL 'DEFA' MANU 'u_{Y}' (jeu*(COS pr_360deg)) 
  331.                                'u_{Z}' (jeu*(SIN pr_360deg));
  332.   t_bidon = prog 0. TFINAL;
  333.   ev_butee_t = EVOL 'DEFA' MANU 't' t_bidon 'r'  (prog 2*jeu);                          
  334.  
  335. * deplacement u(t) et r(t)
  336. *   ev_y = EVOL coco MANU  'LEGE' 'u_{Y}'  't' pr_time 'u' pr_y; 
  337. *   ev_z = EVOL coco MANU  'LEGE' 'u_{Z}'  't' pr_time 'u' pr_z;
  338. *   DESS (ev_y et ev_z) Tdess1;
  339.   pr_r = ((pr_y**2)+(pr_z**2))**0.5;
  340.   ev_r = EVOL coco MANU  'LEGE' 'r'  't' pr_time 'r' pr_r;
  341. si GRAPH;   DESS (ev_butee_t et ev_r); finsi;
  342.  
  343. * orbite complete
  344.   ev_xy = EVOL coco MANU    'u_{Y}' pr_y 'u_{Z}' pr_z; 
  345. *   DESS (ev_xy) 'CARR';
  346. si GRAPH;   DESS (ev_butee et ev_xy) 'CARR'; finsi;
  347.  
  348. * restriction aux nlast dernieres periodes 
  349.   pr_tim2 = extr pr_time llast;
  350.   pr_y2 = extr pr_y llast;
  351.   pr_z2 = extr pr_z llast;
  352.   pr_r2 = ((pr_y2**2)+(pr_z2**2))**0.5;
  353.   ev_xy2 = EVOL coco MANU    'u_{Y}' pr_y2 'u_{Z}' pr_z2; 
  354. si GRAPH;   DESS (ev_butee et ev_xy2) 'CARR'; finsi;
  355.   ev_y2 = EVOL coco MANU   't' pr_tim2 'u_{Y}' pr_y2;
  356. *   ev_r2 = EVOL coco MANU  'LEGE' 'r'  't' pr_tim2 'r' pr_r2;
  357.  
  358. * DSP 
  359.   ev_dsp = ev_dsp et (DSPR ev_y2 11 'FMAX' 40. coco);
  360.  
  361. * stockage  
  362. *   ymax = ymax et (MAXI 'ABS' pr_y2);
  363. *   zmax = zmax et (MAXI 'ABS' pr_z2);
  364.   rmax = rmax et ( MAXI 'ABS' pr_r2);
  365.   rmin = rmin et ( MINI 'ABS' pr_r2);
  366.  
  367. * vitesse tangentielle
  368.   pr_vy = TDYNE . 'VITESSE' . palfa1;
  369.   pr_vz = TDYNE . 'VITESSE' . palfa2;
  370.   ev_v2 = EVOL DEFA MANU  LEGE '|v|' 't' pr_time 'v' (((pr_vy**2)+(pr_vz**2))**0.5);
  371.   pr_vt = TDYNE . TL1 . 'VITESSE_TANGENTIELLE';
  372.   ev_vt = EVOL coco MANU  LEGE 'v_{t}' 't' pr_time 'v' pr_vt;
  373. si GRAPH;   DESS (ev_vt et ev_v2); finsi;
  374.  
  375.   FIN BOMEGA;
  376. *-----------------------------------------------------------------------
  377. *<-------------- FIN DE BOUCLE SUR LA VITESSE DE ROTATION TESTEE  
  378. *-----------------------------------------------------------------------
  379.  
  380. * evolution de r avec la vitesse de rotation
  381.   ev_rmax = EVOL MANU 'w/w_{c}' prWadim 'max(r)/j' (rmax/jeu);
  382.   Tdessr = tabl ; Tdessr . 1 = mot 'NOLI MARQ PLUS';
  383. si GRAPH;  dess ev_rmax Tdessr 'TITRE'  'DYNE06'; finsi;
  384.  
  385. * DSP 
  386. si GRAPH; 
  387.   DESS ev_dsp 'LOGY' 'LEGE' TDSP1 
  388.   'TITX' 'f(Hz)' 'TITY' 'DSP(u_{Y})' 'TITRE' 'DYNE06';
  389. finsi;
  390.  
  391. ************************************************************************
  392. * TEST DE VALIDATION 
  393. ************************************************************************
  394.  
  395. * ON VA CHECKER LE REGIME QUASI-PERIODIQUE A W/WC=0.33
  396.   i033 = POSI 0.33 'DANS' prWadim;
  397.  
  398. * recherche du nombre de pics
  399.   ev_dsp033 = EXTR ev_dsp 'COUR' i033;
  400.   freq033 = EXTR ev_dsp033 'ABSC';
  401.   dsp033 = EXTR ev_dsp033 'ORDO';
  402. * recherche des maxi
  403.   N033 = dime dsp033;
  404.   vari033 = (enle dsp033 N033) - (enle dsp033 1);
  405.   chg033  = (enle vari033 (N033 - 1)) * (enle vari033 1); 
  406. *   imax033 = (POSI 1. 'DANS' (MASQ chg033 INFERIEUR 0.) 'TOUS') + 1;
  407. * on filtre sur les valeurs >  10E-6
  408.   masq6 = MASQ (enle (enle dsp033 N033) 1) 'EGSUPE' 1.E-6;
  409.   chg033 = masq6 * chg033;
  410.   imax033 = (POSI 1. 'DANS' (MASQ chg033 INFERIEUR 0.) 'TOUS') + 1;
  411. list imax033;
  412.   w033 = EXTR  freq033 imax033;
  413. list  w033;
  414.   nw033 = DIME w033;
  415.  
  416. * TNR   
  417.   SI (nw033 > 2);
  418.     MESS 'OK: Presence de plusieurs pics a w/wc=0.33 car quasi-periodique';
  419.   sinon;
  420.     ERREUR 5;
  421.   finsi;
  422.  
  423. * recherche des bornes du mvt radial
  424.   rmin033 = (extr rmin i033)/jeu;
  425.   rmax033 = (extr rmax i033)/jeu;
  426.   mess rmin033 rmax033 ;
  427. * valeurs de [Peletan et al., 2014, Nonlinear Dynamics]
  428.   rminPel = 0.70;
  429.   rmaxPel = 1.45;
  430.   errmin = ABS (rmin033 - rminPel);
  431.   errmax = ABS (rmax033 - rmaxPel);
  432.  
  433. * TNR   
  434.   SI ((errmin < 0.1) et (errmax < 0.1));
  435.     MESS 'OK: Amplitude du regime quasi-periodique a w/wc=0.33 correctes';
  436.   sinon;
  437.     ERREUR 5;
  438.   finsi;
  439.  
  440.  
  441.  
  442. FIN ;
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  

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