Télécharger la plaquette Cast3M


fr en

Télécharger rotor_laval_poutre.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : rotor_laval_poutre.dgibi
  2. *******************************************************************************
  3. *                                                                             *
  4. *  Rotor de Laval                                                             *
  5. *  Etude dans le repere inertiel (ou fixe) avec elements poutre de TIMO       *
  6. *                                                                             *
  7. *              p2   kpal                                                      *
  8. *  z=L          +--/\/\/\--|                                                  *
  9. *               |                                                             *
  10. *               |                                                             *
  11. *               |                                                             *
  12. *            p1a|                                                             *
  13. *        =======+=======   Mdisc, Ixyz                                        *
  14. *            p1b|                                                             *
  15. *               |                                                             *
  16. *               |                                                             *
  17. *               |                                                             *
  18. *  z=0          +--/\/\/\--|                                                  *
  19. *              p0                                                             *
  20. *                                                                             *
  21. *                                                                             *
  22. *  Réf :                                                                      *
  23. *  [1] Vollan, Arne, and Louis Komzsik. "Computational techniques of rotor    *
  24. *      dynamics with the finite element method". CRC Press, 2012.]            *              
  25. *                                                                             *
  26. *  Mots-clés : Vibrations, calcul modal, machines tournantes,                 *
  27. *              poutre, modes complexes, reponse frequentielle                 *
  28. *                                                                             *
  29. *  Auteur: Benoit Prabel, Mars 2020                                           *
  30. *                                                                             *
  31. *******************************************************************************
  32.  
  33. *******************************************************************************
  34. * OPTIONS
  35. *******************************************************************************
  36.  
  37. * dimension, type d'elements geometriques, ...
  38.   OPTI 'DIME' 3 'ELEM' 'SEG2';
  39. * visualisation (sortie postscript)
  40.   GRAPH = FAUX;
  41.   OPTI 'TRAC' PSC 'EPTR' 12 'POTR' 'HELVETICA_16';
  42.  
  43.  
  44. *******************************************************************************
  45. * DONNEES
  46. *******************************************************************************
  47.  
  48. * arbre (shaft)
  49.   L      = 1.0 ;
  50.   Dshaft = 64.E-3;
  51. * Disque (disc)  
  52.   zdisc  = 0.5*L;
  53.   Mdisc  = 40.;
  54.   Izdisc = 5. ;
  55. * Materiau
  56.   E1   = 2.1E11 ;
  57.   nu1  = 0.3 ;
  58.   rho1 = 7800. ;
  59.   Visc1= 0.00001*E1;
  60. * Palier
  61.   kpal_x = 3.9E6;
  62.   kpal_y = kpal_x;
  63.   cpal_x = 1000.;
  64.   cpal_y = cpal_x;
  65.  
  66. * Calculs preliminaires pour completer les donnees
  67.  
  68. * formule des caracteristiques des poutres a section circulaire creuse :
  69. *   Section  = pi * ((Rext**2) - (Rint**2));
  70. *   Itorsion = pi * ((Rext**4) - (Rint**4)) / 2.; 
  71. *   Iflexion = pi * ((Rext**4) - (Rint**4)) / 4.; 
  72. * on a un disque plein ==> Rint = 0
  73. * on a les donnees (Mdisc et Idisc) integrees sur h et multipliee par rho1 :
  74. *   rho1 * h * pi    * (R**2) = Mdisc
  75. *   rho1 * h * pi/2. * (R**4) = Izdisc
  76. * (2)/(1)  ==>  0.5*(R**2) = Izdisc/Mdisc  ==>  Rdisc = (2*Izdisc/Mdisc)**0.5
  77.   Rdisc = (2*Izdisc/Mdisc)**0.5;
  78.   hdisc = Mdisc / (rho1 * pi * (Rdisc**2));
  79.   MESS 'Rdisc=' Rdisc ' hdisc=' hdisc;
  80.  
  81. * ==> caracteristiques geometriques de la poutre definissant le disque :
  82.   Sdisc = pi * (Rdisc**2);
  83.   Iz    = pi/2. * (Rdisc**4);
  84.   Ix    = pi/4. * (Rdisc**4);
  85.  
  86. * caracteristique de l'arbre
  87.   Rshaft = Dshaft /2.;
  88.   Sshaft = pi * (Rshaft**2);
  89.   Izshaft = pi/2. * (Rshaft**4);
  90.   Ixshaft = pi/4. * (Rshaft**4);
  91.  
  92.  
  93. *******************************************************************************
  94. * MAILLAGE
  95. *******************************************************************************
  96.  
  97. * PARAMETRES
  98.   OPTI 'ELEM' SEG2;
  99.   nL    = 80;
  100.   nVect = nL/16;
  101.  
  102. * POINTS  
  103.   p0  = 0. 0. 0.;
  104.   p1a = 0. 0. (zdisc - hdisc);
  105.   p1b = 0. 0. (zdisc + hdisc);
  106.   p2  = 0. 0. L;
  107. * points du palier  
  108.   ppal = p0 et p2;
  109.  
  110. * DROITES
  111.   d1a = DROI (nL/2) p0  p1a;
  112.   d2  = (DROI  1     p1a p1b) COUL 'BLEU';
  113.   d1b = DROI (nL/2) p1b p2 ;
  114.   d1  = d1a   et d1b;
  115.   dtot= d1 et d2;
  116.  
  117. * TRACE
  118. si GRAPH;
  119.   TRAC dtot 'TITRE' 'maillage';
  120. finsi;
  121.  
  122. * ON VEUT QUELQUES POINTS POUR LE TRACE DE VECTEUR DANS POSTVIBR
  123.   dVect = (DROI nVect p0  p1a) ET d2 ET (DROI nVect p1b p2);
  124.   pVect = CHAN dVect 'POI1';
  125.   ELIM dtot pVect (1.E-8*L);
  126.  
  127.  
  128. *******************************************************************************
  129. * MODELE, MATERIAU, MATRICES
  130. *******************************************************************************
  131. *
  132. * 1 : arbre 
  133. mod1 = MODE d1 'MECANIQUE' TIMO;
  134. mat1 = MATE mod1 'YOUNG' E1 'NU' Nu1 'RHO' Rho1
  135.         'SECT' Sshaft 'INRY' Ixshaft 'INRZ' Ixshaft 'TORS' Izshaft
  136.         'OMEG' 1. 'VISQ' Visc1;
  137. *
  138. * 2 : disque
  139. mod2 = MODE d2 'MECANIQUE' TIMO;
  140. mat2 = MATE mod2 'YOUNG' E1 'NU' Nu1 'RHO' Rho1
  141.         'SECT' Sdisc 'INRY' Ix 'INRZ' Ix 'TORS' Iz
  142.         'OMEG' 1. 'VISQ' Visc1;
  143. *
  144. mod12 = mod1 et mod2;
  145. mat12 = mat1 et mat2;
  146.  
  147. * raideur, masse, couplage gyroscopique
  148. K12  = RIGI mod12 mat12;
  149. Mtot = MASS mod12 mat12;
  150. Gtot = GYRO mod12 mat12;
  151. * amortissement + amortissement corotatif
  152. C12 KROT12 = AMOR  mod12 mat12 'COROTATIF';
  153.  
  154. * 3 : paliers 
  155. * appuis = raideurs discretes
  156. Kx3 = APPU 'UX' kpal_x ppal;
  157. Ky3 = APPU 'UY' kpal_x ppal;
  158. K3 = Kx3 et Ky3;
  159. Cx3 = APPU 'UX' cpal_x ppal;
  160. Cy3 = APPU 'UY' cpal_x ppal;
  161. C3 = Cx3 et Cy3;
  162.  
  163. * autres conditions aux limites ?
  164. bloZ = BLOQ 'UZ' 'RZ' ppal;
  165.  
  166. * assemblage
  167. Ktot = K12 et K3 et bloZ;
  168. Ctot = C12 et C3;
  169.  
  170.  
  171. *******************************************************************************
  172. * MODES REELS
  173. *******************************************************************************
  174.  
  175. * PARAMETRES : on souhaite NbModR modes
  176.   NbModR = 4;
  177.  
  178. * CALCUL
  179.   TBasR1 = VIBR 'SIMUL' 1. NbModR Ktot Mtot;
  180.  
  181. * POST-TRAITEMENT 
  182. si GRAPH;
  183.  
  184. * tableau + deformees modale automatise via la procedure postvibr
  185. * on peux customiser avec table d'options
  186.   Toptions = TABL;
  187.   Toptions .'MAILLAGE_VECTEUR' = pVect;
  188.   POSTVIBR TBasR1 Toptions;
  189.  
  190. * quel est ce mode 3 ?
  191.   phi3 = TBasR1 . 'MODES' . 3 . 'DEFORMEE_MODALE';
  192.   TITRE 'Mode 3 : mode de torsion';
  193.   TRAC (EXCO phi3 'RZ') dtot;
  194. * c'est de la torsion !
  195. * si on souhaite le supprimer, il faudrait bloquer plus de ddls RZ
  196. sinon;
  197.   POSTVIBR TBasR1 (MOTS 'TABL');
  198.  
  199. finsi;
  200.  
  201. * VIBR 'SIMUL' a detecte des modes doubles et il a donc ajoute un mode en plus
  202. * retrouvons le vrai nombre avec une mini-boucle...
  203.   NbModR = NbModR - 1 ;
  204.   REPE bb; NbModR = NbModR + 1;
  205.     SI (EXIS TBasR1 . 'MODES' (NbModR + 1));
  206.       ITER bb; 
  207.     SINON; 
  208.       QUIT bb; 
  209.     FINSI;
  210.   FIN bb;
  211.   MESS 'nombre de modes reels fournis in fine par VIBR SIMUL='NbModR;
  212.  
  213.  
  214. *******************************************************************************
  215. * CALCUL DU DIAGRAMME DE CAMPBELL 
  216. * (= EVOLUTION DES FREQUENCES COMPLEXES AVEC LA VITESSE DE ROTATION)
  217. *******************************************************************************
  218.  
  219. * OMEGA en RoundPerMinute
  220.   PR_RPM = prog 0. 'PAS' 0.1E3 10.E3;
  221.  
  222. * on choisit l'unite pour OMEGA : RoundPerMinute, rad/s ou Hz(=tr/s) 
  223. * G est calcule pour 1 rad/s --> multiplier par FAC_G
  224.   UNIT_OMEG = mot 'RPM'  ; FAC_G = (2.*pi/60.); PROMEG = PR_RPM; 
  225. * UNIT_OMEG = mot 'rad/s'; FAC_G = 1.;          PROMEG = (2.*pi/60.) * PR_RPM ;
  226. * UNIT_OMEG = mot 'Hz'   ; FAC_G = 2.*pi;       PROMEG = PR_RPM / 60.;
  227.   cha_x = chai '\W ('UNIT_OMEG')';
  228.   NOMEG = DIME PROMEG;
  229.  
  230. * PROJECTION DES MATRICES ASSEMBLEES SUR LA BASE REELLE
  231.   M1P = PJBA TBasR1 Mtot ;
  232. * G est calcule pour 1 rad/s -> on mutliplie par FAC_G
  233.   G1P = PJBA TBasR1 Gtot ;
  234.   K1P = PJBA TBasR1 (K12 et K3) ;
  235.   C1P = PJBA TBasR1 Ctot;
  236.   KROT1P = PJBA TBasR1 KROT12;
  237.  
  238. * REM : il serait possible d'utiliser directement la procedure campbell,
  239. * mais il semble plus pedagogique de reecrire ici la boucle et le calcul
  240.  
  241. * CREATION DE 2*NbModC LISTREELS DE TAILLE NOMEG
  242.   NbModC = 2*NbModR;
  243.   TfreqR = TABL;
  244.   TfreqI = TABL;
  245.   REPE BmodC NbModC;
  246.     TfreqR . &BmodC = PROG NOMEG*0.;
  247.     TfreqI . &BmodC = PROG NOMEG*0.;
  248.   FIN  BmodC;
  249. * + 2 LISTREELS TEMPORAIRES DE TRAVAIL
  250.   prWorkR = PROG NbModC*0.;
  251.   prWorkI = PROG NbModC*0.;
  252.  
  253. * ON BOUCLE SUR LES FREQUENCES DE ROTATION Omega_j ---------------------
  254.   REPE BOMEG NOMEG;  
  255.     Omega_j = EXTR PROMEG &BOMEG;
  256.     MESS 'Calcul des modes complexes pour \W=' Omega_j ' ' UNIT_OMEG;
  257.  
  258. *   on va resoudre : [ [K + W*C'] + iw [C + W*G] - w^2 [M] ] * \psi = 0
  259. *   G est calcule pour 1 rad/s -> on mutliplie par FAC_G
  260.     K_j = K1P ET (Omega_j * KROT1P);
  261.     C_j = C1P ET (Omega_j * FAC_G * G1P); 
  262.     M_j = M1P;    
  263.     TbasC_j = VIBC M_j K_j C_j ;
  264.  
  265. *   extraction des frequences et tri par ordre croissant
  266.     SI (&BOMEG EGA 1);
  267.       ORDOVIBC TbasC_j;
  268. *   puis en minimisant la distance au resultat precedent
  269.     SINON;
  270.       ORDOVIBC TbasC_j TbasC_jm1;
  271.     FINSI;
  272.     prwR_j = TbasC_j . 'LISTE_FREQUENCES_REELLES'     ;
  273.     prwI_j = TbasC_j . 'LISTE_FREQUENCES_IMAGINAIRES' ;
  274. *   pour le prochain pas 
  275.     TbasC_jm1 = TbasC_j;    
  276.  
  277. *   on stocke dans les listreels finaux
  278.     REPE BmodC NbModC;
  279.       REMP (TfreqR . &BmodC) &BOMEG (EXTR prwR_j &BmodC);
  280.       REMP (TfreqI . &BmodC) &BOMEG (EXTR prwI_j &BmodC);
  281.     FIN  BmodC;
  282.  
  283.   FIN  BOMEG ;
  284. * FIN DE LA BOUCLE SUR LES FREQUENCES DE ROTATION ----------------------
  285.  
  286.  
  287. * POST-TRAITEMENT GRAPHIQUE 
  288. * on ne trace que les modes tq wR>0 -> i ={NbModC/2+1 ... NbModC}
  289.   IFhalf = VRAI;
  290.   si IFhalf;  NC = NbModC/2; ideb = NC;
  291.   sinon;      NC = NbModC;   ideb = 0;
  292.   finsi;
  293.  
  294.   colors = @PALETTE NC;
  295.   evfreqR = VIDE 'EVOLUTIO';
  296.   evfreqI = VIDE 'EVOLUTIO';
  297.   REPE BmodC NC;
  298.     coco = EXTR colors &BmodC;
  299.     i = ideb + &BmodC;
  300.     evfreqR = evfreqR 
  301.     et (EVOL coco 'MANU' cha_x PROMEG 'w_{R} (Hz)' (TfreqR . i));
  302.     evfreqI = evfreqI 
  303.     et (EVOL coco 'MANU' cha_x PROMEG 'w_{I} (/s)' (TfreqI . i));
  304.   FIN  BmodC;
  305.  
  306. si GRAPH;
  307.   TITRE 'Campbell diagram';
  308.   DESS evfreqR ;
  309.   DESS evfreqI ;
  310. finsi;
  311.  
  312.  
  313. * POST-TRAITEMENT GRAPHIQUE 
  314. * On recombine les deformees Complexes issues de VIBC
  315.   RECOVIBC TbasC_j TBasR1;
  316. * Listing + trace des deformees Complexes issues de VIBC
  317.   Topt = TABL;
  318.   Topt . 'MAILLAGE_VECTEUR' = ppal;
  319.   POSTVIBR TbasC_j Topt;
  320.  
  321.  
  322. * opti donn 5 trac X;
  323. *******************************************************************************
  324. * TEST DE NON REGRESSION
  325. *******************************************************************************
  326.  
  327. * ON TESTE JUSTE LE SUIVI
  328.  
  329. * preparation au calcul de la derivee 
  330. * abscisse
  331.   dOMEG  = ((ENLE PROMEG NOMEG) - (ENLE PROMEG 1));
  332.   NdOMEG = NOMEG - 1;
  333. * message
  334.   opti echo 0;
  335.   MESS (CHAI 'Mode'*4 'E[w]'*18 'E[dwR/dW]'*33 'max[dwR/dW]'*48 'E[dwI/dW]'*63 'max[dwI/dW]'*78);
  336. * les courbes etant asse lineaire xtol = 10 est largement suffisant
  337.   xtol = 10.;
  338.  
  339. * boucle sur les modes Complexes (i.e. sur les courbes)
  340.   REPE BmodC NbModC;
  341. *   --- PARTIE REELLE ---
  342.     wR = TfreqR . &BmodC;
  343.     wRmoy = (SOMM wR) / NOMEG;
  344. *   calcul de la derivee 
  345.     dwRdW = ((ENLE wR NOMEG) - (ENLE wR 1)) / dOMEG;
  346. *   moyenne et valeur max
  347.     dwRdWmoy = (SOMM dwRdW) / NdOMEG;
  348.     dwRdWmax = MAXI 'ABS' dwRdW;
  349. *   --- PARTIE IMAGINAIRE ---
  350.     wI = TfreqR . &BmodC;
  351.     wImoy = (SOMM wI) / NOMEG;
  352. *   calcul de la derivee 
  353.     dwIdW = ((ENLE wI NOMEG) - (ENLE wI 1)) / dOMEG;
  354. *   moyenne et valeur max
  355.     dwIdWmoy = (SOMM dwIdW) / NdOMEG;
  356.     dwIdWmax = MAXI 'ABS' dwIdW;
  357. *   --- MESSAGE + TEST ---
  358. *   message
  359.     MESS (CHAI &BmodC*4 wRmoy*18 dwRdWmoy*33 dwRdWmax*48 dwIdWmoy*63 dwIdWmax*78);    
  360. *   test sur la derivee pour verifier la continuite
  361.     dwdWref = MAXI 'ABS' (PROG dwRdWmoy dwIdWmoy 1.E-6);
  362.     SI (dwRdWmax > (xtol*dwdWref)); ERRE 5; FINSI;
  363.     SI (dwIdWmax > (xtol*dwdWref)); ERRE 5; FINSI;
  364.   FIN  BmodC;
  365.   opti echo 1;
  366.  
  367.  
  368. FIN ;
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  

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