* fichier : rotor_laval_poutre.dgibi ******************************************************************************* * * * Rotor de Laval * * Etude dans le repere inertiel (ou fixe) avec elements poutre de TIMO * * * * p2 kpal * * z=L +--/\/\/\--| * * | * * | * * | * * p1a| * * =======+======= Mdisc, Ixyz * * p1b| * * | * * | * * | * * z=0 +--/\/\/\--| * * p0 * * * * * * Réf : * * [1] Vollan, Arne, and Louis Komzsik. "Computational techniques of rotor * * dynamics with the finite element method". CRC Press, 2012.] * * * * Mots-clés : Vibrations, calcul modal, machines tournantes, * * poutre, modes complexes, reponse frequentielle * * * * Auteur: Benoit Prabel, Mars 2020 * * * ******************************************************************************* ******************************************************************************* * OPTIONS ******************************************************************************* * dimension, type d'elements geometriques, ... OPTI 'DIME' 3 'ELEM' 'SEG2'; * visualisation (sortie postscript) GRAPH = FAUX; OPTI 'TRAC' PSC 'EPTR' 12 'POTR' 'HELVETICA_16'; ******************************************************************************* * DONNEES ******************************************************************************* * arbre (shaft) L = 1.0 ; Dshaft = 64.E-3; * Disque (disc) zdisc = 0.5*L; Mdisc = 40.; Izdisc = 5. ; * Materiau E1 = 2.1E11 ; nu1 = 0.3 ; rho1 = 7800. ; Visc1= 0.00001*E1; * Palier kpal_x = 3.9E6; kpal_y = kpal_x; cpal_x = 1000.; cpal_y = cpal_x; * Calculs preliminaires pour completer les donnees * formule des caracteristiques des poutres a section circulaire creuse : * Section = pi * ((Rext**2) - (Rint**2)); * Itorsion = pi * ((Rext**4) - (Rint**4)) / 2.; * Iflexion = pi * ((Rext**4) - (Rint**4)) / 4.; * on a un disque plein ==> Rint = 0 * on a les donnees (Mdisc et Idisc) integrees sur h et multipliee par rho1 : * rho1 * h * pi * (R**2) = Mdisc * rho1 * h * pi/2. * (R**4) = Izdisc * (2)/(1) ==> 0.5*(R**2) = Izdisc/Mdisc ==> Rdisc = (2*Izdisc/Mdisc)**0.5 Rdisc = (2*Izdisc/Mdisc)**0.5; hdisc = Mdisc / (rho1 * pi * (Rdisc**2)); MESS 'Rdisc=' Rdisc ' hdisc=' hdisc; * ==> caracteristiques geometriques de la poutre definissant le disque : Sdisc = pi * (Rdisc**2); Iz = pi/2. * (Rdisc**4); Ix = pi/4. * (Rdisc**4); * caracteristique de l'arbre Rshaft = Dshaft /2.; Sshaft = pi * (Rshaft**2); Izshaft = pi/2. * (Rshaft**4); Ixshaft = pi/4. * (Rshaft**4); ******************************************************************************* * MAILLAGE ******************************************************************************* * PARAMETRES OPTI 'ELEM' SEG2; nL = 80; nVect = nL/16; * POINTS p0 = 0. 0. 0.; p1a = 0. 0. (zdisc - hdisc); p1b = 0. 0. (zdisc + hdisc); p2 = 0. 0. L; * points du palier ppal = p0 et p2; * DROITES d1a = DROI (nL/2) p0 p1a; d2 = (DROI 1 p1a p1b) COUL 'BLEU'; d1b = DROI (nL/2) p1b p2 ; d1 = d1a et d1b; dtot= d1 et d2; * TRACE si GRAPH; TRAC dtot 'TITRE' 'maillage'; finsi; * ON VEUT QUELQUES POINTS POUR LE TRACE DE VECTEUR DANS POSTVIBR dVect = (DROI nVect p0 p1a) ET d2 ET (DROI nVect p1b p2); pVect = CHAN dVect 'POI1'; ELIM dtot pVect (1.E-8*L); ******************************************************************************* * MODELE, MATERIAU, MATRICES ******************************************************************************* * * 1 : arbre mod1 = MODE d1 'MECANIQUE' TIMO; mat1 = MATE mod1 'YOUNG' E1 'NU' Nu1 'RHO' Rho1 'SECT' Sshaft 'INRY' Ixshaft 'INRZ' Ixshaft 'TORS' Izshaft 'OMEG' 1. 'VISQ' Visc1; * * 2 : disque mod2 = MODE d2 'MECANIQUE' TIMO; mat2 = MATE mod2 'YOUNG' E1 'NU' Nu1 'RHO' Rho1 'SECT' Sdisc 'INRY' Ix 'INRZ' Ix 'TORS' Iz 'OMEG' 1. 'VISQ' Visc1; * mod12 = mod1 et mod2; mat12 = mat1 et mat2; * raideur, masse, couplage gyroscopique K12 = RIGI mod12 mat12; Mtot = MASS mod12 mat12; Gtot = GYRO mod12 mat12; * amortissement + amortissement corotatif C12 KROT12 = AMOR mod12 mat12 'COROTATIF'; * 3 : paliers * appuis = raideurs discretes Kx3 = APPU 'UX' kpal_x ppal; Ky3 = APPU 'UY' kpal_x ppal; K3 = Kx3 et Ky3; Cx3 = APPU 'UX' cpal_x ppal; Cy3 = APPU 'UY' cpal_x ppal; C3 = Cx3 et Cy3; * autres conditions aux limites ? bloZ = BLOQ 'UZ' 'RZ' ppal; * assemblage Ktot = K12 et K3 et bloZ; Ctot = C12 et C3; ******************************************************************************* * MODES REELS ******************************************************************************* * PARAMETRES : on souhaite NbModR modes NbModR = 4; * CALCUL TBasR1 = VIBR 'SIMUL' 1. NbModR Ktot Mtot; * POST-TRAITEMENT si GRAPH; * tableau + deformees modale automatise via la procedure postvibr * on peux customiser avec table d'options Toptions = TABL; Toptions .'MAILLAGE_VECTEUR' = pVect; POSTVIBR TBasR1 Toptions; * quel est ce mode 3 ? phi3 = TBasR1 . 'MODES' . 3 . 'DEFORMEE_MODALE'; TITRE 'Mode 3 : mode de torsion'; TRAC (EXCO phi3 'RZ') dtot; * c'est de la torsion ! * si on souhaite le supprimer, il faudrait bloquer plus de ddls RZ sinon; POSTVIBR TBasR1 (MOTS 'TABL'); finsi; * VIBR 'SIMUL' a detecte des modes doubles et il a donc ajoute un mode en plus * retrouvons le vrai nombre avec une mini-boucle... NbModR = NbModR - 1 ; REPE bb; NbModR = NbModR + 1; SI (EXIS TBasR1 . 'MODES' (NbModR + 1)); ITER bb; SINON; QUIT bb; FINSI; FIN bb; MESS 'nombre de modes reels fournis in fine par VIBR SIMUL='NbModR; ******************************************************************************* * CALCUL DU DIAGRAMME DE CAMPBELL * (= EVOLUTION DES FREQUENCES COMPLEXES AVEC LA VITESSE DE ROTATION) ******************************************************************************* * OMEGA en RoundPerMinute PR_RPM = prog 0. 'PAS' 0.1E3 10.E3; * on choisit l'unite pour OMEGA : RoundPerMinute, rad/s ou Hz(=tr/s) * G est calcule pour 1 rad/s --> multiplier par FAC_G UNIT_OMEG = mot 'RPM' ; FAC_G = (2.*pi/60.); PROMEG = PR_RPM; * UNIT_OMEG = mot 'rad/s'; FAC_G = 1.; PROMEG = (2.*pi/60.) * PR_RPM ; * UNIT_OMEG = mot 'Hz' ; FAC_G = 2.*pi; PROMEG = PR_RPM / 60.; cha_x = chai '\W ('UNIT_OMEG')'; NOMEG = DIME PROMEG; * PROJECTION DES MATRICES ASSEMBLEES SUR LA BASE REELLE M1P = PJBA TBasR1 Mtot ; * G est calcule pour 1 rad/s -> on mutliplie par FAC_G G1P = PJBA TBasR1 Gtot ; K1P = PJBA TBasR1 (K12 et K3) ; C1P = PJBA TBasR1 Ctot; KROT1P = PJBA TBasR1 KROT12; * REM : il serait possible d'utiliser directement la procedure campbell, * mais il semble plus pedagogique de reecrire ici la boucle et le calcul * CREATION DE 2*NbModC LISTREELS DE TAILLE NOMEG NbModC = 2*NbModR; TfreqR = TABL; TfreqI = TABL; REPE BmodC NbModC; TfreqR . &BmodC = PROG NOMEG*0.; TfreqI . &BmodC = PROG NOMEG*0.; FIN BmodC; * + 2 LISTREELS TEMPORAIRES DE TRAVAIL prWorkR = PROG NbModC*0.; prWorkI = PROG NbModC*0.; * ON BOUCLE SUR LES FREQUENCES DE ROTATION Omega_j --------------------- REPE BOMEG NOMEG; Omega_j = EXTR PROMEG &BOMEG; MESS 'Calcul des modes complexes pour \W=' Omega_j ' ' UNIT_OMEG; * on va resoudre : [ [K + W*C'] + iw [C + W*G] - w^2 [M] ] * \psi = 0 * G est calcule pour 1 rad/s -> on mutliplie par FAC_G K_j = K1P ET (Omega_j * KROT1P); C_j = C1P ET (Omega_j * FAC_G * G1P); M_j = M1P; TbasC_j = VIBC M_j K_j C_j ; * extraction des frequences et tri par ordre croissant SI (&BOMEG EGA 1); ORDOVIBC TbasC_j; * puis en minimisant la distance au resultat precedent SINON; ORDOVIBC TbasC_j TbasC_jm1; FINSI; prwR_j = TbasC_j . 'LISTE_FREQUENCES_REELLES' ; prwI_j = TbasC_j . 'LISTE_FREQUENCES_IMAGINAIRES' ; * pour le prochain pas TbasC_jm1 = TbasC_j; * on stocke dans les listreels finaux REPE BmodC NbModC; REMP (TfreqR . &BmodC) &BOMEG (EXTR prwR_j &BmodC); REMP (TfreqI . &BmodC) &BOMEG (EXTR prwI_j &BmodC); FIN BmodC; FIN BOMEG ; * FIN DE LA BOUCLE SUR LES FREQUENCES DE ROTATION ---------------------- * POST-TRAITEMENT GRAPHIQUE * on ne trace que les modes tq wR>0 -> i ={NbModC/2+1 ... NbModC} IFhalf = VRAI; si IFhalf; NC = NbModC/2; ideb = NC; sinon; NC = NbModC; ideb = 0; finsi; colors = @PALETTE NC; evfreqR = VIDE 'EVOLUTIO'; evfreqI = VIDE 'EVOLUTIO'; REPE BmodC NC; coco = EXTR colors &BmodC; i = ideb + &BmodC; evfreqR = evfreqR et (EVOL coco 'MANU' cha_x PROMEG 'w_{R} (Hz)' (TfreqR . i)); evfreqI = evfreqI et (EVOL coco 'MANU' cha_x PROMEG 'w_{I} (/s)' (TfreqI . i)); FIN BmodC; si GRAPH; TITRE 'Campbell diagram'; DESS evfreqR ; DESS evfreqI ; finsi; * POST-TRAITEMENT GRAPHIQUE * On recombine les deformees Complexes issues de VIBC RECOVIBC TbasC_j TBasR1; * Listing + trace des deformees Complexes issues de VIBC Topt = TABL; Topt . 'MAILLAGE_VECTEUR' = ppal; POSTVIBR TbasC_j Topt; * opti donn 5 trac X; ******************************************************************************* * TEST DE NON REGRESSION ******************************************************************************* * ON TESTE JUSTE LE SUIVI * preparation au calcul de la derivee * abscisse dOMEG = ((ENLE PROMEG NOMEG) - (ENLE PROMEG 1)); NdOMEG = NOMEG - 1; * message opti echo 0; MESS (CHAI 'Mode'*4 'E[w]'*18 'E[dwR/dW]'*33 'max[dwR/dW]'*48 'E[dwI/dW]'*63 'max[dwI/dW]'*78); * les courbes etant asse lineaire xtol = 10 est largement suffisant xtol = 10.; * boucle sur les modes Complexes (i.e. sur les courbes) REPE BmodC NbModC; * --- PARTIE REELLE --- wR = TfreqR . &BmodC; wRmoy = (SOMM wR) / NOMEG; * calcul de la derivee dwRdW = ((ENLE wR NOMEG) - (ENLE wR 1)) / dOMEG; * moyenne et valeur max dwRdWmoy = (SOMM dwRdW) / NdOMEG; dwRdWmax = MAXI 'ABS' dwRdW; * --- PARTIE IMAGINAIRE --- wI = TfreqR . &BmodC; wImoy = (SOMM wI) / NOMEG; * calcul de la derivee dwIdW = ((ENLE wI NOMEG) - (ENLE wI 1)) / dOMEG; * moyenne et valeur max dwIdWmoy = (SOMM dwIdW) / NdOMEG; dwIdWmax = MAXI 'ABS' dwIdW; * --- MESSAGE + TEST --- * message MESS (CHAI &BmodC*4 wRmoy*18 dwRdWmoy*33 dwRdWmax*48 dwIdWmoy*63 dwIdWmax*78); * test sur la derivee pour verifier la continuite dwdWref = MAXI 'ABS' (PROG dwRdWmoy dwIdWmoy 1.E-6); SI (dwRdWmax > (xtol*dwdWref)); ERRE 5; FINSI; SI (dwIdWmax > (xtol*dwdWref)); ERRE 5; FINSI; FIN BmodC; opti echo 1; FIN ;