* fichier : hbm_jeffcott_contact_alfa.dgibi ****************************************************************** * * Systeme a deux ddl * * ROTOR JEFFCOTT AVEC CONTACT ROTOR-STATOR * .. . * m x + c x + k x = Fbal*cos wt + Fchoc_x * .. . * m y + c y + k y = Fbal*sin wt + Fchoc_y * * modele de contact : * Fn = -{r-jeu}+ kchoc * Ft = µ|Fn|sign(vrel) * * Continuation + HBM + AFT * idem hbm_jeffcott_contact.dgibi mais 1 composante (ALFA) * sur 2 noeuds au lieu de UY UZ sur le meme noeud. * pas d'autre interet => a supprimer un jour ? * ****************************************************************** *----------------------------------------------------------------------- * CHARMECA.PROCEDUR * Calcul de la force non lineaire fnl(x,t) et de sa derivee dfnl(x,t)/dx * dans le domaine temporel *----------------------------------------------------------------------- * * EN ENTREE : * * TAB1 * >> parametres génériques : * . 'DEP_NL' : LISTCHPO DE DEPLACEMENTS TEMPORELLES * ou (pas encore programmé) * . i . 'CHPOINT' = CHPOINT unitaire des deplacements NL * . i . 'LISTREEL' = evolution temporelle coefficient du chpoint * >> parametres propres a ce cas test : * . 'COEFF_FNL' : coefficient de non linearite * . 'POINT_FNL' : point de non linearite * . 'COMP_FNL' : composante de non linearite * * * EN SORTIE : * * TAB1 * >> parametres génériques : * . 'FNL_T' : LISTCHPO DES FORCES NON LINEAIRES TEMPORELLES * ou TABLE d'indice : (pas encore programmé) * . i . 'CHPOINT' = CHPOINT unitaire des forces non lineaires * . i . 'LISTREEL' = evolution temporelle coefficient du chpoint * . 'KNL_T' : DERIVEE DE FORCES NON LINEAIRES sous forme de TABLE : * . i . j = 'LISTREEL' evolution temporelle de dFi/dxj * *----------------------------------------------------------------------- DEBPROC CHARMECA TAB1*'TABLE'; ******** parametres d'entree propres a ce cas test ********** *coefficient de la force nl SI (EXIS TAB1 . 'COEFF_FNL'); pmu = EXTR TAB1 . 'COEFF_FNL' 1; pkc = EXTR TAB1 . 'COEFF_FNL' 2; ph0 = EXTR TAB1 . 'COEFF_FNL' 3; SINON; MESS 'CHARMECA: IL MANQUE Les COEFF_FNL'; ERRE 21; FINSI; SI (EXIS TAB1 . 'POINT_FNL'); pNL = TAB1 . 'POINT_FNL'; SINON; MESS 'CHARMECA: IL MANQUE LA GEOMETRIE NL dans POINT_FNL'; ERRE 21; FINSI; SI (EXIS TAB1 . 'COMP_FNL'); COMP_FNL = TAB1 . 'COMP_FNL'; SINON; MESS 'CHARMECA: IL MANQUE LA COMPOSANTE NL dans COMP_FNL'; ERRE 21; FINSI; * ******** parametres d'entree génériques ********** *DEPNL : list de champ par points, deplacements temporels du point nl * SI (EXIS TAB1 'DEP_NL'); DEPNL = TAB1 . 'DEP_NL'; IFCHPO = ega (type DEPNL) 'LISTCHPO'; si (non IFCHPO); si (neg (type DEPNL) 'TABLE'); mess 'CHARMECA: DEP_NL doit etre de type LISTCHPO ou TABLE'; ERRE 21; sinon; si (neg (type DEPNL . 1) 'LISTREEL'); mess 'CHARMECA: DEP_NL . i doit etre un LISTREEL'; ERRE 21; finsi; finsi; finsi; SINON; MESS 'CHARMECA: IL MANQUE LA LISTE DES DEPLACEMENTS dans DEP_NL'; ERRE 21; FINSI; * calcul de fnl ET dfnl/dx ? dans le doute on calcule tout... SI (NEG (TYPE IFFNL) 'LOGIQUE'); IFFNL = VRAI; FINSI; SI (NEG (TYPE IFKNL) 'LOGIQUE'); IFKNL = VRAI; FINSI; ******** Calcul de FNL ******** *nb: nombre de pas de temps *NCNL : nombre de composantes pour la force nl nb = DIME DEPNL; FNLTEMP = VIDE 'LISTCHPO'; fxrp = prog ; fyrp =prog ; dfxdxp = prog ; dfxdyp = prog ; dfydxp = prog ; dfydyp = prog ; LST_R0 = PROG; LST_X0 = PROG; LST_Y0 = PROG; NCNL = NBNO pNL; pNL1 = POIN pNL 1; pNL2 = POIN pNL 2; COMPALFA = EXTR COMP_FNL 1; * BOUCLE POUR CONSTRUIRE le LISTREEL EXCENTRICITE(t) ib = 0; REPE bfor nb; ib = ib + 1; DEP_i0 = EXTR DEPNL ib; DEP_X0 = EXTR DEP_i0 pNL1 COMPALFA; DEP_Y0 = EXTR DEP_i0 pNL2 COMPALFA; * mess 'x,y = ' DEP_X0 DEP_Y0; LST_X0 = LST_X0 ET DEP_X0; LST_Y0 = LST_Y0 ET DEP_Y0; * EXCENTRICITE R0 = ((DEP_X0**2) + (DEP_Y0**2))**0.5; LST_R0 = LST_R0 ET R0; FIN bfor; *---------- BOUCLE POUR CALCULER LA FORCE NL ----------- ib = 0; REPE Blst nb; ib = ib + 1; DEP_X0 = EXTR LST_X0 ib; DEP_Y0 = EXTR LST_Y0 ib; R0 = EXTR LST_R0 ib; SI (R0 0 * (hyp : vitesse de roation >> deplacement du centre du disque) costet1 = DEP_X0 / R0; sintet1 = DEP_Y0 / R0; fxr = (fn * costet1) - (ft * sintet1); fyr = (fn * sintet1) + (ft * costet1); si IFKNL; cstet1 = costet1*sintet1; costet2 = costet1**2; sintet2 = sintet1**2; unmjeu1 = 1. - (ph0/R0); A1 = 1. - (ph0/R0 * sintet2); B1 = ph0/R0 * cstet1; C1 = B1; D1 = 1. - (ph0/R0 * costet2); dfxdx = pkc * ( A1 - (pmu * B1)); dfxdy = pkc * ( C1 - (pmu * D1)); dfydx = pkc * ( B1 + (pmu * A1)); dfydy = pkc * ( D1 + (pmu * C1)); * On fournit en realite les termes de Knl = - dF/dX finsi; FINSI; fxrp = fxrp et fxr; fyrp = fyrp et fyr; FNLTEMP = FNLTEMP et ((MANU 'CHPO' pNL1 1 COMPALFA fxr 'NATURE' 'DISCRET') et (MANU 'CHPO' pNL2 1 COMPALFA fyr 'NATURE' 'DISCRET')); si IFKNL; dfxdxp = dfxdxp et dfxdx; dfxdyp = dfxdyp et dfxdy; dfydxp = dfydxp et dfydx; dfydyp = dfydyp et dfydy; finsi; FIN Blst; * mess 'f^nl = ' fxr fyr; * mess 'dfnl/du=' dfxdx dfxdy dfydx dfydy; TAB1 . 'FNL_T' = FNLTEMP; TAB1 . 'KNL_T' = table; si IFKNL; TAB1 . 'KNL_T' . 1 = tabl; TAB1 . 'KNL_T' . 2 = tabl; TAB1 . 'KNL_T' . 1 . 1 = dfxdxp; TAB1 . 'KNL_T' . 1 . 2 = dfxdyp; TAB1 . 'KNL_T' . 2 . 1 = dfydxp; TAB1 . 'KNL_T' . 2 . 2 = dfydyp; finsi; fprog = (fxrp**2 + fyrp**2)**0.5; fnlmax = MAXI fprog 'ABS'; si (EXIS TAB1 'FNLMAX'); fnlmax0 = TAB1 . 'FNLMAX'; sinon; fnlmax0 = 0.; finsi; * si chgt brusque (mise en contact ou perte de contact), on actualise K si ( ((fnlmax > 1.E-50) et (fnlmax0 < 1.E-50)) OU ((fnlmax < 1.E-50) et (fnlmax0 > 1.E-50)) ); * mess 'contact detecté !'; WTAB . 'RECA_K' = VRAI; finsi; TAB1 . 'FNLMAX' = fnlmax; TAB1 . 'EXCENTRICITE' = (MAXI (ABS LST_R0)); FINPROC; *----------------------------------------------------------------------- ************************************************************************ * * DEBUT DU CALCUL * * ************************************************************************ *opti echo 0 ; opti debug 1; graph = faux ; OPTI DIME 2 ELEM SEG2 MODE PLAN CONT ; GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI; * GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX; * *-------- Definition des options du calcul : * complet = faux; * CALCUL PAR DFT ? (Transformation de Fourier Discret CALC_DFT = FAUX; *si analyser la stabilite? OP_STAB = FAUX; * *si adimensionnement du temps? * Adim_t = VRAI; * *-------- Definition des parametres du calcul -------- *nombre d'harmoniques nhbm = 3 ; *nhbm = 7 ; * nhbm = 9 ; * nhbm = 11 ; *nhbm = 15 ; *plage d'etude de frequence en Hz t0 = 0.0; t1 = 10.; dt = 0.10; *des iterations itmoy = 6; ds = 1.; dsmax = 1.; dsmin = 1.E-3 ; * si AFT, OTFR: nb de points 2**OTFR pour TFR * OTFR = 7; OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2)); mess 'OTFR = ' OTFR; *-------- Definition des parametres du contact --------------- *coefficient de frottement pmu = 0.11 ; si (pmu ega 0.00 1.E-4); mopmu = mot '000'; finsi; si (pmu ega 0.05 1.E-4); mopmu = mot '005'; finsi; si (pmu ega 0.11 1.E-4); mopmu = mot '011'; finsi; si (pmu ega 0.15 1.E-4); mopmu = mot '015'; finsi; si (pmu ega 0.20 1.E-4); mopmu = mot '020'; finsi; si (pmu ega 0.25 1.E-4); mopmu = mot '025'; finsi; *rigidite du contact pkc = 2500. ; * pkc = 0.; *jeu initial ph0 = 0.105 ; * *parametre de l'adimensionnement = 20% du jeu u1max = 0.20*ph0; * *---------Definition des parametres du rotor -------- *masse constant pmass = 1. ; *amortissement pamor = 5. ; *rigidite prigi = 100. ; *force du balourd pbal = 0.1; * frequences caracteristiques du pb (en Hz) w0 = (prigi/pmass)**0.5 / (2.*pi); w0c = ((prigi+pkc)/pmass)**0.5 / (2.*pi); wc = ((pkc)/pmass)**0.5 / (2.*pi); mess w0 w0c wc; * prefix : prefix = chai 'hbm_jeffcott_contact_alfa_' mopmu '_n' nhbm ; TITRE prefix; si GRAPH; ficps = chai prefix '.ps'; opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps; finsi; *------------- geometrie -------- * P1 = 0. 0. ; P2 = 1. 0. ; *----------- construction des matrices caracteristiques ---------- * *OPTI DONN 5; LIM1 = MOTS 'ALFA'; MASS1 = MANU 'RIGI' 'TYPE' 'MASSE' (P1 et P2) LIM1 (prog pmass) ; RI1 = MANU 'RIGI' (P1 et P2) LIM1 (prog prigi) ; AMOR1 = MANU 'RIGI' 'TYPE' 'AMORTISSEMENT' (P1 et P2) LIM1 (prog pamor); *------- on définit une unique table pour tout (HBM et CONTINU) ------- TAB1 = TABLE ; *--------- remplissage des matrices "temporelles" --------- * TAB1 . 'MASSE_CONSTANTE' = MASS1 ; * TAB1 . 'AMORTISSEMENT_CONSTANT' = AMOR1 ; * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici TAB1 . 'MASSE_CONSTANTE' = (2.*pi **2) * MASS1 ; TAB1 . 'AMORTISSEMENT_CONSTANT' = (2.*pi) * AMOR1 ; TAB1 . 'RIGIDITE_CONSTANTE' = RI1 ; TAB1 . 'N_HARMONIQUE' = nhbm; *------- remplissage des resultats attendus sur ddls "temporels" ------- mycoul6 = MOTS 'VIOL' 'BLEU' 'BLEU' 'TURQ' 'TURQ' 'VERT' 'VERT' 'OLIV' 'OLIV' 'JAUN' 'JAUN' 'ORAN' 'ORAN' 'ROUG' 'ROUG' 'ROSE' 'ROSE' ; mycoul6 = mycoul6 et (MOTS 24*'GRIS'); TAB1 . 'RESULTATS' = tabl; TAB1 . 'RESULTATS' . 1 = tabl; TAB1 . 'RESULTATS' . 1 . 'POINT_MESURE' = P1; TAB1 . 'RESULTATS' . 1 . 'COMPOSANTE' = mot 'ALFA'; TAB1 . 'RESULTATS' . 1 . 'TITRE' = mot 'UX'; TAB1 . 'RESULTATS' . 1 . 'COULEUR' = mot 'BLEU'; TAB1 . 'RESULTATS' . 1 . 'COULEUR_HBM' = mycoul6; TAB1 . 'RESULTATS' . 2 = tabl; TAB1 . 'RESULTATS' . 2 . 'POINT_MESURE' = P2; TAB1 . 'RESULTATS' . 2 . 'COMPOSANTE' = mot 'ALFA'; TAB1 . 'RESULTATS' . 2 . 'TITRE' = mot 'UY'; TAB1 . 'RESULTATS' . 2 . 'COULEUR' = mot 'ROSE'; TAB1 . 'RESULTATS' . 2 . 'COULEUR_HBM' = mycoul6; *--------- passage en frequentiel --------- HBM TAB1; MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP TEMPORELLES D UN NOEUD'; LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT'; MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP FREQUENTIELLES D UN NOEUD'; LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT_HBM'; *----------- definition du chargement ----------------------------------- * * ici, on definit directement le chargement en frequentiel FP11 = (MANU 'CHPO' p1 1 'F2' pbal 'NATURE' 'DISCRET') et (MANU 'CHPO' p2 1 'G2' pbal 'NATURE' 'DISCRET'); LIX1 = PROG 0. pas 0.1 (t1+30.) ; * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici LIY1 = (2.*pi*LIX1)**2 ; EV1 = EVOL MANU 't' LIX1 'F(t)' LIY1 ; CHA1 = CHAR 'MECA' FP11 EV1 ; si GRAPH; dess EV1 POSX CENT POSY CENT ; finsi; *----------- resolution par la procedure CONTINU ----------------------- * OPTIONS DE CALCULS TAB1 . 'CHARGEMENT' = CHA1; * calcul des forces NL TAB1 . 'PROCEDURE_CHARMECA'= VRAI ; TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT'; TAB1 . 'N_PT_TFR' = OTFR; TAB1 . 'CALC_CHPO' = VRAI; * TAB1 . 'CALC_CHPO' = faux; -> charmeca pas prevu ! * *plage de frequence en Hz LIS1T = prog t0 PAS dt t1; TAB1 . 'TEMPS_CALCULES' = LIS1T; TAB1 . 'MAXIPAS' = 300; TAB1 . 'PRECISION' = 1.E-5; TAB1 . 'ACCELERATION' = 4; *------parametres du contact TAB1 . 'COEFF_FNL' = PROG pmu pkc ph0; TAB1 . 'POINT_FNL' = (P1 et P2); TAB1 . 'COMP_FNL' = mots 'ALFA'; *------parametres du continuation TAB1 . 'NB_ITERATION' = itmoy; TAB1 . 'MAXITERATION' = 30; TAB1 . 'WTABLE' = tabl; TAB1 . 'WTABLE' . 'DS' = ds; TAB1 . 'WTABLE' . 'DSMAX' = dsmax; TAB1 . 'WTABLE' . 'DSMIN' = dsmin; TAB1 . 'MAXI_DEPLACEMENT' = u1max; TAB1 . 'MAXSOUSPAS' = 3; * stabilité TAB1 . 'STABILITE' = vrai; * opti debug 1 ; TEMP ZERO; *calcul avec continuation CONTINU TAB1; TEMP 'IMPR' 'MAXI' 'CPU'; * *---------------------- sauv avant reprise ---------------------------- * ficxdr = chai prefix '-avantREPRISE.xdr'; * mess ficxdr; * opti sauv ficxdr; * sauv ; *------------------------- reprise ----------------------------- * on impose l'absence de contact TAB1 . 'COEFF_FNL' = PROG (0.*pmu) (0.*pkc) (1E6*ph0); * on fait croire qu'on vient dans la direction opposee de maniere * a etre capable de traiter le point de rebroussement WTAB1 = TAB1 . 'WTABLE'; * OUBL TAB1 'WTABLE'; list WTAB1; WTAB2 = TABL ; WTAB2 . 'DTIME0' = -1. * WTAB1 . 'DTIME0'; WTAB2 . 'DDEP0' = -1. * WTAB1 . 'DDEP0'; TAB1 . 'WTABLE' = WTAB2; CONTINU TAB1; *------------------------- post-traitement ----------------------------- ********** tracé de toutes les harmoniques individuellement ********** wprog = TAB1 . 'TEMPS_PROG'; * wprog = (extr evtot cour 1) extr ABSC; evtot = TAB1 . 'RESULTATS_HBM' . 'RESULTATS_EVOL'; TDESS1 = TABL; TDESS1 . 'TITRE' = TAB1 . 'RESULTATS_HBM' . 'TITRE'; si GRAPH; dess evtot TITX '\w(Hz)' POSX CENT TITY 'U_{k}' POSY CENT TDESS1 'LEGE'; dess evtot YBOR -0.5 0.3 XBOR 0. 12 TITX '\w(Hz)' POSX CENT TITY 'U_{k}' POSY CENT TDESS1 'LEGE'; finsi; ********** traduction des resultats frequentiels en temporels ********** HBM_POST TAB1; ************* tracé du max d'amplitude max_t(u(t)) ********************* evuyamp = TAB1 . 'RESULTATS' . 'RESULTATS_EVOL'; evuyamp1 = extr evuyamp COUR 1; TDESS2 = TABL; TDESS2 . 2 = mot 'TIRR'; TDESS2 . 'TITRE' = TAB1 . 'RESULTATS' . 'TITRE'; si GRAPH; dess evuyamp1 TITX '\w(Hz)' POSX CENT TITY 'max|U(t)|' POSY CENT TDESS2 LEGE 'NO'; dess evuyamp1 TITX '\w(Hz)' POSX CENT TITY 'max|U(t)|' POSY CENT TDESS2 LEGE ; finsi; ************* tracé du max d'amplitude avec la stabilite ! ************* istab = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'STABILITE'; Tdess2 = TABL; Tdess2 . 'LIGNE_VARIABLE' = TABL; Tdess2 . 'LIGNE_VARIABLE' . 1 = istab; si GRAPH; dess evuyamp1 TITX '\w(Hz)' POSX CENT TITY 'max|U(t)|' POSY CENT Tdess2; Tdess2 . 1 = mot 'MARQ S PLUS'; dess evuyamp1 TITX '\w(Hz)' POSX CENT TITY 'max|U(t)|' POSY CENT Tdess2; finsi; * autres tracé de la stabilite : lambda_R en fonction de W mycoul = mots 'VIOL' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'OLIV' 'JAUN' 'ORAN' 'ROUG' 'ROSE'; ncoul = dime mycoul; wprog2 = enle wprog (dime wprog); evlitot = VIDE 'EVOLUTIO' ; evlrtot = VIDE 'EVOLUTIO' ; nl = dime TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL'; il = 0; icoul = 0; repe bl nl; il = il + 1; icoul = icoul + 1; si(icoul > ncoul); icoul = 1; finsi; coul1 = extr mycoul icoul; lr1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL' . il; evlrtot = evlrtot et (evol coul1 MANU '\w (Hz)' wprog2 '\l_{R}' lr1); li1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_IMAG' . il; evlitot = evlitot et (evol coul1 MANU '\w (Hz)' wprog2 '\l_{I}' li1); fin bl; wprog3 = prog (mini wprog2) (maxi wprog2); evzero = evol manu '\w (Hz)' wprog3 '\l_{R}' (0.*wprog3); Tdess3 = tabl; Tdess3 . 1 = mot 'TIRR'; si GRAPH; dess (evzero et evlrtot) 'TITY' '\l_{R} (s^{-1})' Tdess3; dess (evlitot) 'TITY' '\l_{I} (Hz)' ; finsi; * stabilite : trace dynamique pour faire une animation si GRAPH; ficps = chai prefix '-ORBITE.ps'; opti 'FTRA' ficps; TABORB = tabl; TABORB . 'PAS' = 3; TABORB . 'QUEUE' = MOT 'INFINIE'; TABORB . 'EVOL_FIXE' = evzero ; ORBITE evlrtot TABORB; finsi; *------------------------- sauvegarde ----------------------------- si SAVE ; ficxdr = chai prefix '.xdr'; mess ficxdr; opti sauv ficxdr; sauv evtot evuyamp istab; finsi; *--------------------- test de non regression --------------------- si TNR ; * reference : these de Lihan Xie fig3.10 page 84 * sauf qu'ici, on ne prend que 3 harmoniques * recup valeurs de la courbes de reponse wamp = extr evuyamp 'ABSC' 1; wadim = wamp / wc; uyamp = extr evuyamp 'ORDO' 1; ramp = uyamp; * R=UY puisque l'on a des orbites circulaires, * sinon il faudrait calculer r(t)=uy(t)**2+uz(t)**2 , puis R=max(r(t)) namp = DIME wamp; * recherche de la mise en contact idebut = POSI 1 'DANS' (ENTI (MASQ ramp EGSUPE ph0)); wdeb0 = extr wadim (idebut - 1); wdeb1 = extr wadim idebut; si ((wdeb0 > 0.16) ou (wdeb1 < 0.16)); mess 'Mise en contact devrait se produire vers w/w0~0.16 !'; mess 'PB car : ' wdeb0 ' < 0.16 < 'wdeb1 ' non verifie !'; ERRE 5; finsi; * verification que le calcul a bien attenit les w=10Hz wmax = maxi wamp; si (wmax < 10.); mess 'On devrait calculer jusqu a 10Hz au moins !'; mess wmax; ERRE 5; finsi; * verification de l'amplitude max atteinte * rem : rmax = 4.4807 lors de la creation du cas test, * mais on arrondit pour rref rmax = (maxi uyamp) / ph0; rref = 4.5 ; si ((ABS (rmax - rref)) > 0.05); mess 'Le max|U| devrait etre de ' rref ' !'; mess rmax; ERRE 5; finsi; * recherche des changements de stabilite istab1 = (istab enle 1) - (istab enle (DIME istab)); istab1 = (posi 1 'DANS' istab1 'TOUS') et (posi -1 'DANS' istab1 'TOUS'); si (neg (dime istab1) 4); mess 'Il devrait y avoir 4 points de changement de stabilite !'; list istab1; ERRE 5; finsi; finsi; fin ;