* fichier : dyne06.dgibi * ************************************************************************ * * DYNE06 : * * Calcul d'un rotor de type Jeffcott avec contact frottant * * ___ ^ Z * _|_ kc,mu |__ * _______ jeu |_ \_ * | k,c | | | \ \ * |----------| m |----------| +--|--|-----> Y * | |_______| R R+jeu * ___ jeu * _|_ kc,mu * * ^ Z * | * +----> X * * Test de l'option VITESSE_ENTRAINEMENT pour la liaison DYNE * POINT_CERCLE_FROTTEMENT qui ajoute une sur-vitesse tangentielle * * Ref : [Xie et al., MSSP, 2015] * * Auteur : BP, 2020-03-25 * ************************************************************************ ************************************************************************ * OPTIONS et DONNEES ************************************************************************ * options OPTI DIME 3 ELEM SEG2; GRAPH = FAUX; COMPLET = FAUX; OPTI TRAC PSC EPTR 8 POTR HELVETICA_16; prefix = chai 'DYNE06'; TITR prefix; * rotor data m = 1. ; k = 100.; c = 5. ; jeu = 0.105; R = 20.*jeu ; ebal = 0.1; * Fbal = ebal * m * (Omega**2) ; avec Omega = vitesse de rotation en rad/s * contact parameters kchoc= 25. * k; mu_s = 0.2; mu_d = 0.2; si (mu_s < mu_d); erreur 21; finsi; Kt = 10.*kchoc; Ct = 0.5*(m*Kt)**0.5; * deduced parameter w = (k/m)**0.5; wHz = w / (2.*pi); wchoc = (kchoc/m)**0.5; xi = 0.5 * c / ((kchoc*m)**0.5); mess 'w=' w '(' wHz ')Hz < wchoc=' wchoc; mess 'xi=' xi; * critical timestep (Diff Centrees, no damping) wadhe = ((k+Kt)/m)**0.5; dtc = 2. / wadhe; ************************************************************************ * CALCUL ANALYTIQUE LINEAIRE ************************************************************************ * si pas de contact, alors pas de couplage entre y et z * on resout pour y par exemple : * * m y" + c y' + ky = meW^2 cos Wt * * on pose y(t) = A cos Wt + B sin Wt * y' = -Aw sin Wt + Bw cos Wt * y" = -Aw^2 cos Wt - Bw^2 sin Wt * * [ k-W^2m Wc ] (A) = (meW^2) termes en cos * [ -Wc k-W^2m ] (B) = ( 0 ) termes en sin * * 2eme ligne : B = Wc/k-W^2m * A * 1ere ligne : A = meW^2 / [ k-W^2m + W^2c^2/k-W^2m] 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; Wanal = wchoc * Wadim; Wanal2 = Wanal**2; * m=1 -> ca simplifie aux1 = k - Wanal2; aux2 = (c**2) * Wanal2; denom1 = aux1 + (aux2/aux1); A_anal = (ebal * Wanal2) / denom1; B_anal = (c*Wanal) / aux1 * A_anal; * evA = evol VIOL manu Wanal A_anal; * evB = evol ROSE manu Wanal B_anal; * dess (evA et evB); Amp_anal = ( (A_anal**2) + (B_anal**2) )**0.5; evAmp = evol VIOL manu Wanal Amp_anal; evjeu = evol manu Wanal (prog (dime Wanal)*jeu); * dess (evjeu et evAmp); * en adimensionne evjeu2 = evol manu Wadim (prog (dime Wanal)*1); evAmp2 = evol VIOL manu 'w/w_{c}' Wadim 'r/j' (Amp_anal/jeu); si GRAPH; dess (evjeu2 et evAmp2) TITRE 'solution analytique lineaire'; finsi; * opti donn 5 trac X; ************************************************************************ * BASE MODALE ************************************************************************ * point physique p1 = 0. 0. 0.; * base modale = 2 modes palfa1 = 0. 0. 0.; phi1 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 1. 'UZ' 0. 'NATURE' 'DIFFUS'; palfa2 = 0. 0. 0.; phi2 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 0. 'UZ' 1. 'NATURE' 'DIFFUS'; TBAS = TABL 'BASE_MODALE'; TBAS . 'MODES' = TABL 'BASE_DE_MODES'; TBAS . 'MODES' . 1 = TABL 'MODES'; TBAS . 'MODES' . 1 . 'POINT_REPERE' = palfa1; TBAS . 'MODES' . 1 . 'FREQUENCE' = wHz; TBAS . 'MODES' . 1 . 'MASSE_GENERALISEE' = m; TBAS . 'MODES' . 1 . 'DEFORMEE_MODALE' = phi1; TBAS . 'MODES' . 2 = TABL 'MODES'; TBAS . 'MODES' . 2 . 'POINT_REPERE' = palfa2; TBAS . 'MODES' . 2 . 'FREQUENCE' = wHz; TBAS . 'MODES' . 2 . 'MASSE_GENERALISEE' = m; TBAS . 'MODES' . 2 . 'DEFORMEE_MODALE' = phi2; ************************************************************************ * AMORTISSEMENT ************************************************************************ TAMOR = TABL 'AMORTISSEMENT'; Ctot = MANU 'RIGIDITE' 'TYPE' 'AMORTISSEMENT' (palfa1 et palfa2) (MOTS 'ALFA') (prog c); TAMOR . 'AMORTISSEMENT' = Ctot ; ************************************************************************ * LIAISONS ************************************************************************ * contact-frottant point_cercle_frottement TL1 = TABL 'LIAISON_ELEMENTAIRE'; TL1 . 'TYPE_LIAISON' = mot 'POINT_CERCLE_FROTTEMENT'; TL1 . 'SUPPORT' = p1; TL1 . 'NORMALE' = (1. 0. 0.) ; TL1 . 'EXCENTRATION' = (0. 0. 0.) ; TL1 . 'RAIDEUR' = kchoc; TL1 . 'RAYON' = jeu ; TL1 . 'COEFFICIENT_GLISSEMENT' = mu_d; TL1 . 'COEFFICIENT_ADHERENCE' = mu_s; TL1 . 'RAIDEUR_TANGENTIELLE' = Kt; TL1 . 'AMORTISSEMENT_TANGENTIEL'= Ct; * stockage TLB = TABL 'LIAISON_B'; TLB . 1 = TL1; TLIA = TABL 'LIAISON'; TLIA . 'LIAISON_B' = TLB; ************************************************************************ * CONDITIONS INITIALES ************************************************************************ * TINIT = TABL 'INITIAL'; * TINIT . 'DEPLACEMENT' * = (MANU 'CHPO' palfa1 1 'ALFA' ... 'NATURE' 'DISCRET') * et (MANU 'CHPO' palfa2 1 'ALFA' ... 'NATURE' 'DISCRET'); * TINIT . 'VITESSE' * = (MANU 'CHPO' palfa1 1 'ALFA' ... 'NATURE' 'DIFFUS') * et (MANU 'CHPO' palfa2 1 'ALFA' ... 'NATURE' 'DIFFUS'); *u=v=0 ici --> inutile ************************************************************************ * SORTIES ************************************************************************ TSORT = TABLE 'SORTIE'; * +on stocke les deplacements modaux TSORV = TABLE 'VARIABLE'; TSORV . 'TYPE_SORTIE' = MOT 'LISTREEL'; TSORT . 'VARIABLE' = TSORV; * + on stocke les variables du contact TSORLB = TABLE 'LIAISON_B'; TSORLB . TL1 = VRAI; TSORT . 'LIAISON_B' = TSORLB; * pour le DESSin des evolutions Tdess1 = TABL ; Tdess1 . 2 = MOT 'TIRR'; TDSP1 = TABL; TDSP1 . 'TITRE' = TABL; ************************************************************************ * PARAMETRES TEMPORELS ************************************************************************ * NOMBRE DE PERIODES SIMULEES, NOMBRE DE PAS PAR PERIODE * NPERIOD = 8; *debug NPERIOD = 14; NPERIOD = 80; NPPP = 2**10; * NOMBRE DE PAS TOTAL, SORTIE TOUS LES NSORT NPAS = NPERIOD*NPPP; NSORT = 8; * NUMERO DES POINTS SORTIS APPARTENANT AUX nlast DERNIERES PERIODES * nlast = 2; nlast = 20; llast = LECT (NPPP / NSORT * (NPERIOD - nlast)) 'PAS' 1 (NPPP / NSORT * NPERIOD + 1); *----------------------------------------------------------------------- *--------------> BOUCLE SUR LA VITESSE DE ROTATION TESTEE *----------------------------------------------------------------------- * vitesse de rotation (prOmega en rad/s) si COMPLET; prWadim = PROG 0.1 0.15 0.154 0.2 0.290 0.3 0.33 0.4 0.5 1. 1.2 ; coco_dyne = MOTS 'INDI' 'MARI' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'BOUT' 'OR' 'ORAN' 'ROUG' 'BRUN' ; sino; prWadim = PROG 0.15 0.2 0.33 0.4 1. ; coco_dyne = MOTS 'MARI' 'AZUR' 'BOUT' 'OR' 'ROUG' ; finsi; prOmega = wchoc * prWadim; * RESULTATS * ymax = prog; * zmax = prog; rmax = prog; rmin = prog; ev_dsp = VIDE 'EVOLUTIO'; REPE BOMEGA (DIME prOmega); *debug &BOMEGA=1; *debug &BOMEGA=2; *debug &BOMEGA=4; coco = EXTR coco_dyne &BOMEGA; Omega = EXTR prOmega &BOMEGA; OMEGHZ = Omega / (2.*PI); MESS '>>> Calcul #' &BOMEGA ' pour OMEGA='Omega 'rad/s (adim ='(Omega/wchoc) ')'; FMT = MOT '(F6.3)'; cha1 = chai '\W=' (enti Omega) ' rad/s - \w/\w_{c}=' FORMAT FMT (Omega/wchoc); TITRE cha1; TDSP1 . 'TITRE' . &BOMEGA = chai '\w/\w_{c}='FORMAT FMT (Omega/wchoc); * ON COMPLETE : *-LIAISON *** TL1 . 'VITESSE_ENTRAINEMENT' = Omega; TL1 . 'VITESSE_ENTRAINEMENT' = -1.*R*Omega; * *-PERIODE, PAS DE TEMPS, TEMPS FINAL TPERIOD = 1. / OMEGHZ; DT = TPERIOD / NPPP; TFINAL = NPAS*DT; si (DT < dtc); mess 'Pas de temps =' DT ' < critique=' dtc; sinon; erre 21; finsi; ************************************************************************ * CHARGEMENT ************************************************************************ * nombre de periodes pour la rampe de chargement NPrampe = 10; * evolution temporelle du chargement t_char = PROG (-1.*DT) 'PAS' DT (TFINAL + DT); rampe = BORN (t_char/(NPrampe*TPERIOD)) 'COMPRIS' 0. 1.; F_Y = COS (360.*OMEGHZ * t_char); F_Z = SIN (360.*OMEGHZ * t_char); *debug F_Y = 1.; *debug F_Z = 0.; evF_Y = EVOL 'ROUG' MANU t_char (rampe * F_Y); evF_Z = EVOL 'BRIQ' MANU t_char (rampe * F_Z); * DESS (evF_Y et evF_Z) TITRE 'Chargement' XBOR 0. TFINAL XGRA TPERIOD; * description "spatiale" du chargement Fbal = ebal*m * (Omega**2); Fchpo_Y = MANU 'CHPO' 1 palfa1 'FALF' Fbal; Fchpo_Z = MANU 'CHPO' 1 palfa2 'FALF' Fbal; CHAR_Y = CHAR evF_Y Fchpo_Y ; CHAR_Z = CHAR evF_Z Fchpo_Z ; * table pour DYNE TCHARG = TABLE 'CHARGEMENT' ; TCHARG . 'BASE_A' = CHAR_Y et CHAR_Z ; ************************************************************************ * CALCUL DYNE : INTEGRATION TEMPORELLE ************************************************************************ * TDYNE = DYNE 'DE_VOGELAERE' TBAS TAMOR TLIA TSORT TCHARG TDYNE = DYNE 'DIFFERENCES_CENTREES' TBAS TAMOR TLIA TSORT TCHARG NPAS DT NSORT ; ************************************************************************ * POST TRAITEMENT ************************************************************************ * temps t pr_time = TDYNE . 'TEMPS_DE_SORTIE'; * deplacement modaux pr_y = TDYNE . 'DEPLACEMENT' . palfa1; pr_z = TDYNE . 'DEPLACEMENT' . palfa2; * cercle r=R+jeu pr_360deg = prog 0. PAS 2.5 360.; ev_butee = EVOL 'DEFA' MANU 'u_{Y}' (jeu*(COS pr_360deg)) 'u_{Z}' (jeu*(SIN pr_360deg)); t_bidon = prog 0. TFINAL; ev_butee_t = EVOL 'DEFA' MANU 't' t_bidon 'r' (prog 2*jeu); * deplacement u(t) et r(t) * ev_y = EVOL coco MANU 'LEGE' 'u_{Y}' 't' pr_time 'u' pr_y; * ev_z = EVOL coco MANU 'LEGE' 'u_{Z}' 't' pr_time 'u' pr_z; * DESS (ev_y et ev_z) Tdess1; pr_r = ((pr_y**2)+(pr_z**2))**0.5; ev_r = EVOL coco MANU 'LEGE' 'r' 't' pr_time 'r' pr_r; si GRAPH; DESS (ev_butee_t et ev_r); finsi; * orbite complete ev_xy = EVOL coco MANU 'u_{Y}' pr_y 'u_{Z}' pr_z; * DESS (ev_xy) 'CARR'; si GRAPH; DESS (ev_butee et ev_xy) 'CARR'; finsi; * restriction aux nlast dernieres periodes pr_tim2 = extr pr_time llast; pr_y2 = extr pr_y llast; pr_z2 = extr pr_z llast; pr_r2 = ((pr_y2**2)+(pr_z2**2))**0.5; ev_xy2 = EVOL coco MANU 'u_{Y}' pr_y2 'u_{Z}' pr_z2; si GRAPH; DESS (ev_butee et ev_xy2) 'CARR'; finsi; ev_y2 = EVOL coco MANU 't' pr_tim2 'u_{Y}' pr_y2; * ev_r2 = EVOL coco MANU 'LEGE' 'r' 't' pr_tim2 'r' pr_r2; * DSP ev_dsp = ev_dsp et (DSPR ev_y2 11 'FMAX' 40. coco); * stockage * ymax = ymax et (MAXI 'ABS' pr_y2); * zmax = zmax et (MAXI 'ABS' pr_z2); rmax = rmax et ( MAXI 'ABS' pr_r2); rmin = rmin et ( MINI 'ABS' pr_r2); * vitesse tangentielle pr_vy = TDYNE . 'VITESSE' . palfa1; pr_vz = TDYNE . 'VITESSE' . palfa2; ev_v2 = EVOL DEFA MANU LEGE '|v|' 't' pr_time 'v' (((pr_vy**2)+(pr_vz**2))**0.5); pr_vt = TDYNE . TL1 . 'VITESSE_TANGENTIELLE'; ev_vt = EVOL coco MANU LEGE 'v_{t}' 't' pr_time 'v' pr_vt; si GRAPH; DESS (ev_vt et ev_v2); finsi; FIN BOMEGA; *----------------------------------------------------------------------- *<-------------- FIN DE BOUCLE SUR LA VITESSE DE ROTATION TESTEE *----------------------------------------------------------------------- * evolution de r avec la vitesse de rotation ev_rmax = EVOL MANU 'w/w_{c}' prWadim 'max(r)/j' (rmax/jeu); Tdessr = tabl ; Tdessr . 1 = mot 'NOLI MARQ PLUS'; si GRAPH; dess ev_rmax Tdessr 'TITRE' 'DYNE06'; finsi; * DSP si GRAPH; DESS ev_dsp 'LOGY' 'LEGE' TDSP1 'TITX' 'f(Hz)' 'TITY' 'DSP(u_{Y})' 'TITRE' 'DYNE06'; finsi; ************************************************************************ * TEST DE VALIDATION ************************************************************************ * ON VA CHECKER LE REGIME QUASI-PERIODIQUE A W/WC=0.33 i033 = POSI 0.33 'DANS' prWadim; * recherche du nombre de pics ev_dsp033 = EXTR ev_dsp 'COUR' i033; freq033 = EXTR ev_dsp033 'ABSC'; dsp033 = EXTR ev_dsp033 'ORDO'; * recherche des maxi N033 = dime dsp033; vari033 = (enle dsp033 N033) - (enle dsp033 1); chg033 = (enle vari033 (N033 - 1)) * (enle vari033 1); * imax033 = (POSI 1. 'DANS' (MASQ chg033 INFERIEUR 0.) 'TOUS') + 1; * on filtre sur les valeurs > 10E-6 masq6 = MASQ (enle (enle dsp033 N033) 1) 'EGSUPE' 1.E-6; chg033 = masq6 * chg033; imax033 = (POSI 1. 'DANS' (MASQ chg033 INFERIEUR 0.) 'TOUS') + 1; list imax033; w033 = EXTR freq033 imax033; list w033; nw033 = DIME w033; * TNR SI (nw033 > 2); MESS 'OK: Presence de plusieurs pics a w/wc=0.33 car quasi-periodique'; sinon; ERREUR 5; finsi; * recherche des bornes du mvt radial rmin033 = (extr rmin i033)/jeu; rmax033 = (extr rmax i033)/jeu; mess rmin033 rmax033 ; * valeurs de [Peletan et al., 2014, Nonlinear Dynamics] rminPel = 0.70; rmaxPel = 1.45; errmin = ABS (rmin033 - rminPel); errmax = ABS (rmax033 - rmaxPel); * TNR SI ((errmin < 0.1) et (errmax < 0.1)); MESS 'OK: Amplitude du regime quasi-periodique a w/wc=0.33 correctes'; sinon; ERREUR 5; finsi; FIN ;