Télécharger hbm_jeffcott_contact.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : hbm_jeffcott_contact.dgibi
  2. ******************************************************************
  3. *
  4. * Systeme a deux ddl
  5. *
  6. * ROTOR JEFFCOTT AVEC CONTACT ROTOR-STATOR
  7. * .. .
  8. * m x + c x + k x = Fbal*cos wt + Fchoc_x
  9. * .. .
  10. * m y + c y + k y = Fbal*sin wt + Fchoc_y
  11. *
  12. * modele de contact :
  13. * Fn = -{r-jeu}+ kchoc
  14. * Ft = µ|Fn|sign(vrel)
  15. *
  16. * Continuation + HBM + AFT
  17. *
  18. ******************************************************************
  19.  
  20. *-----------------------------------------------------------------------
  21. * CHARMECA.PROCEDUR
  22. * Calcul de la force non lineaire fnl(x,t) et de sa derivee dfnl(x,t)/dx
  23. * dans le domaine temporel
  24. *-----------------------------------------------------------------------
  25. *
  26. * EN ENTREE :
  27. *
  28. * TAB1
  29. * >> parametres génériques :
  30. * . 'DEP_NL' : LISTCHPO DE DEPLACEMENTS TEMPORELLES
  31. * ou (pas encore programmé)
  32. * . i . 'CHPOINT' = CHPOINT unitaire des deplacements NL
  33. * . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
  34. * >> parametres propres a ce cas test :
  35. * . 'COEFF_FNL' : coefficient de non linearite
  36. * . 'POINT_FNL' : point de non linearite
  37. * . 'COMP_FNL' : composante de non linearite
  38. *
  39. *
  40. * EN SORTIE :
  41. *
  42. * TAB1
  43. * >> parametres génériques :
  44. * . 'FNL_T' : LISTCHPO DES FORCES NON LINEAIRES TEMPORELLES
  45. * ou TABLE d'indice : (pas encore programmé)
  46. * . i . 'CHPOINT' = CHPOINT unitaire des forces non lineaires
  47. * . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
  48. * . 'KNL_T' : DERIVEE DE FORCES NON LINEAIRES sous forme de TABLE :
  49. * . i . j = 'LISTREEL' evolution temporelle de dFi/dxj
  50. *
  51. *-----------------------------------------------------------------------
  52. DEBPROC CHARMECA TAB1*'TABLE';
  53.  
  54. ******** parametres d'entree propres a ce cas test **********
  55.  
  56. *coefficient de la force nl
  57. SI (EXIS TAB1 . 'COEFF_FNL');
  58. pmu = EXTR TAB1 . 'COEFF_FNL' 1;
  59. pkc = EXTR TAB1 . 'COEFF_FNL' 2;
  60. ph0 = EXTR TAB1 . 'COEFF_FNL' 3;
  61. SINON;
  62. MESS 'CHARMECA: IL MANQUE Les COEFF_FNL'; ERRE 21;
  63. FINSI;
  64.  
  65. SI (EXIS TAB1 . 'POINT_FNL');
  66. pNL = TAB1 . 'POINT_FNL';
  67. SINON;
  68. MESS 'CHARMECA: IL MANQUE LA GEOMETRIE NL dans POINT_FNL'; ERRE 21;
  69. FINSI;
  70.  
  71. SI (EXIS TAB1 . 'COMP_FNL');
  72. COMP_FNL = TAB1 . 'COMP_FNL';
  73. n_CNL = DIME COMP_FNL;
  74. SINON;
  75. MESS 'CHARMECA: IL MANQUE LA COMPOSANTE NL dans COMP_FNL'; ERRE 21;
  76. FINSI;
  77. *
  78. ******** parametres d'entree génériques **********
  79.  
  80. *DEPNL : list de champ par points, deplacements temporels du point nl
  81. *
  82. SI (EXIS TAB1 'DEP_NL');
  83. DEPNL = TAB1 . 'DEP_NL';
  84. IFCHPO = ega (type DEPNL) 'LISTCHPO';
  85. si (non IFCHPO);
  86. si (neg (type DEPNL) 'TABLE');
  87. mess 'CHARMECA: DEP_NL doit etre de type LISTCHPO ou TABLE';
  88. ERRE 21;
  89. sinon;
  90. si (neg (type DEPNL . 1) 'LISTREEL');
  91. mess 'CHARMECA: DEP_NL . i doit etre un LISTREEL';
  92. ERRE 21;
  93. finsi;
  94. finsi;
  95. finsi;
  96. SINON;
  97. MESS 'CHARMECA: IL MANQUE LA LISTE DES DEPLACEMENTS dans DEP_NL';
  98. ERRE 21;
  99. FINSI;
  100.  
  101. * calcul de fnl ET dfnl/dx ? dans le doute on calcule tout...
  102. SI (NEG (TYPE IFFNL) 'LOGIQUE'); IFFNL = VRAI; FINSI;
  103. SI (NEG (TYPE IFKNL) 'LOGIQUE'); IFKNL = VRAI; FINSI;
  104.  
  105. ******** Calcul de FNL ********
  106.  
  107. *nb: nombre de pas de temps
  108. *NCNL : nombre de composantes pour la force nl
  109. nb = DIME DEPNL;
  110. FNLTEMP = VIDE 'LISTCHPO';
  111. fxrp = prog ; fyrp =prog ;
  112. dfxdxp = prog ; dfxdyp = prog ; dfydxp = prog ; dfydyp = prog ;
  113. LST_R0 = PROG; LST_X0 = PROG; LST_Y0 = PROG;
  114.  
  115. NCNL = DIME COMP_FNL;
  116. CNL1 = mot (EXTR COMP_FNL 1);
  117. CNL2 = mot (EXTR COMP_FNL 2);
  118. * BOUCLE POUR CONSTRUIRE le LISTREEL EXCENTRICITE(t)
  119. ib = 0;
  120. REPE bfor nb; ib = ib + 1;
  121. DEP_i0 = EXTR DEPNL ib;
  122. DEP_X0 = EXTR DEP_i0 pNL CNL1;
  123. DEP_Y0 = EXTR DEP_i0 pNL CNL2;
  124. * mess 'x,y = ' DEP_X0 DEP_Y0;
  125. LST_X0 = LST_X0 ET DEP_X0;
  126. LST_Y0 = LST_Y0 ET DEP_Y0;
  127. * EXCENTRICITE
  128. R0 = ((DEP_X0**2) + (DEP_Y0**2))**0.5;
  129. LST_R0 = LST_R0 ET R0;
  130. FIN bfor;
  131. *---------- BOUCLE POUR CALCULER LA FORCE NL -----------
  132. ib = 0;
  133. REPE Blst nb; ib = ib + 1;
  134. DEP_X0 = EXTR LST_X0 ib;
  135. DEP_Y0 = EXTR LST_Y0 ib;
  136. R0 = EXTR LST_R0 ib;
  137. SI (R0 <EG ph0);
  138. fxr = 0.;
  139. fyr = 0.;
  140. si IFKNL;
  141. dfxdx = 0. ;
  142. dfxdy = 0. ;
  143. dfydx = 0. ;
  144. dfydy = 0. ;
  145. finsi;
  146. SINON;
  147. fn = -1. * pkc * (R0 - ph0);
  148. ft = pmu * fn;
  149. * car on suppose que v_tangentielle_relative >0
  150. * (hyp : vitesse de roation >> deplacement du centre du disque)
  151. costet1 = DEP_X0 / R0;
  152. sintet1 = DEP_Y0 / R0;
  153. fxr = (fn * costet1) - (ft * sintet1);
  154. fyr = (fn * sintet1) + (ft * costet1);
  155. si IFKNL;
  156. cstet1 = costet1*sintet1;
  157. costet2 = costet1**2;
  158. sintet2 = sintet1**2;
  159. unmjeu1 = 1. - (ph0/R0);
  160. A1 = 1. - (ph0/R0 * sintet2);
  161. B1 = ph0/R0 * cstet1;
  162. C1 = B1;
  163. D1 = 1. - (ph0/R0 * costet2);
  164. dfxdx = pkc * ( A1 - (pmu * B1));
  165. dfxdy = pkc * ( C1 - (pmu * D1));
  166. dfydx = pkc * ( B1 + (pmu * A1));
  167. dfydy = pkc * ( D1 + (pmu * C1));
  168. * On fournit en realite les termes de Knl = - dF/dX
  169. finsi;
  170. FINSI;
  171. fxrp = fxrp et fxr;
  172. fyrp = fyrp et fyr;
  173. FNLTEMP = FNLTEMP ET
  174. (MANU 'CHPO' pNL 2 CNL1 fxr CNL2 fyr 'NATURE' 'DISCRET');
  175. si IFKNL;
  176. dfxdxp = dfxdxp et dfxdx;
  177. dfxdyp = dfxdyp et dfxdy;
  178. dfydxp = dfydxp et dfydx;
  179. dfydyp = dfydyp et dfydy;
  180. finsi;
  181. FIN Blst;
  182. * mess 'f^nl = ' fxr fyr;
  183. * mess 'dfnl/du=' dfxdx dfxdy dfydx dfydy;
  184.  
  185. TAB1 . 'FNL_T' = FNLTEMP;
  186. TAB1 . 'KNL_T' = table;
  187. si IFKNL;
  188. TAB1 . 'KNL_T' . 1 = tabl;
  189. TAB1 . 'KNL_T' . 2 = tabl;
  190. TAB1 . 'KNL_T' . 1 . 1 = dfxdxp;
  191. TAB1 . 'KNL_T' . 1 . 2 = dfxdyp;
  192. TAB1 . 'KNL_T' . 2 . 1 = dfydxp;
  193. TAB1 . 'KNL_T' . 2 . 2 = dfydyp;
  194. finsi;
  195.  
  196. fprog = (fxrp**2 + fyrp**2)**0.5;
  197. fnlmax = MAXI fprog 'ABS';
  198. si (EXIS TAB1 'FNLMAX'); fnlmax0 = TAB1 . 'FNLMAX';
  199. sinon; fnlmax0 = 0.;
  200. finsi;
  201. * si chgt brusque (mise en contact ou perte de contact), on actualise K
  202. si ( ((fnlmax > 1.E-50) et (fnlmax0 < 1.E-50))
  203. OU ((fnlmax < 1.E-50) et (fnlmax0 > 1.E-50)) );
  204. * mess 'contact detecté !';
  205. WTAB . 'RECA_K' = VRAI;
  206. finsi;
  207.  
  208. TAB1 . 'FNLMAX' = fnlmax;
  209. TAB1 . 'EXCENTRICITE' = (MAXI (ABS LST_R0));
  210.  
  211. FINPROC;
  212. *-----------------------------------------------------------------------
  213.  
  214.  
  215.  
  216. ************************************************************************
  217. *
  218. * DEBUT DU CALCUL *
  219. *
  220. ************************************************************************
  221. *opti echo 0 ;
  222. opti debug 1;
  223. graph = faux ;
  224. OPTI DIME 2 ELEM SEG2 MODE PLAN CONT ;
  225. OPTI EPSI LINEAIRE;
  226.  
  227. GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI;
  228. * GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX;
  229.  
  230. *
  231. *-------- Definition des options du calcul :
  232. * complet = faux;
  233. * CALCUL PAR DFT ? (Transformation de Fourier Discret
  234. CALC_DFT = FAUX;
  235. *si analyser la stabilite?
  236. OP_STAB = FAUX;
  237. * *si adimensionnement du temps?
  238. * Adim_t = VRAI;
  239. *
  240. *-------- Definition des parametres du calcul --------
  241. *nombre d'harmoniques
  242. nhbm = 3 ;
  243. *nhbm = 7 ;
  244. * nhbm = 9 ;
  245. * nhbm = 11 ;
  246. *nhbm = 15 ;
  247. *plage d'etude de frequence en Hz
  248. t0 = 0.0;
  249. t1 = 10.;
  250. dt = 0.10;
  251. *des iterations
  252. itmoy = 6;
  253. ds = 1.;
  254. dsmax = 1.;
  255. dsmin = 1.E-3 ;
  256. * si AFT, OTFR: nb de points 2**OTFR pour TFR
  257. * OTFR = 7;
  258. OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2));
  259. mess 'OTFR = ' OTFR;
  260.  
  261. *-------- Definition des parametres du contact ---------------
  262. *coefficient de frottement
  263. pmu = 0.11 ;
  264. si (pmu ega 0.00 1.E-4); mopmu = mot '000'; finsi;
  265. si (pmu ega 0.05 1.E-4); mopmu = mot '005'; finsi;
  266. si (pmu ega 0.11 1.E-4); mopmu = mot '011'; finsi;
  267. si (pmu ega 0.15 1.E-4); mopmu = mot '015'; finsi;
  268. si (pmu ega 0.20 1.E-4); mopmu = mot '020'; finsi;
  269. si (pmu ega 0.25 1.E-4); mopmu = mot '025'; finsi;
  270. *rigidite du contact
  271. pkc = 2500. ;
  272. * pkc = 0.;
  273. *jeu initial
  274. ph0 = 0.105 ;
  275. *
  276. *parametre de l'adimensionnement = 20% du jeu
  277. u1max = 0.20*ph0;
  278. *
  279. *---------Definition des parametres du rotor --------
  280. *masse constant
  281. pmass = 1. ;
  282. *amortissement
  283. pamor = 5. ;
  284. *rigidite
  285. prigi = 100. ;
  286. *force du balourd
  287. pbal = 0.1;
  288.  
  289. * frequences caracteristiques du pb (en Hz)
  290. w0 = (prigi/pmass)**0.5 / (2.*pi);
  291. w0c = ((prigi+pkc)/pmass)**0.5 / (2.*pi);
  292. wc = ((pkc)/pmass)**0.5 / (2.*pi);
  293. mess w0 w0c wc;
  294.  
  295. * prefix :
  296. prefix = chai 'hbm_jeffcott_contact_' mopmu '_n' nhbm ;
  297. TITRE prefix;
  298. si GRAPH;
  299. ficps = chai prefix '.ps';
  300. opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps;
  301. finsi;
  302.  
  303. *------------- geometrie --------
  304. *
  305. P1 = 0. 0. ; P2 = 1. 0. ;
  306.  
  307. *----------- construction des matrices caracteristiques ----------
  308. *
  309. *OPTI DONN 5;
  310. LIM1 = MOTS UX UY;
  311. MASS1 = (MASS 'UX' pmass P2 ) ET (MASS 'UY' pmass P2 );
  312. RI1 = MANU 'RIGI' P2 LIM1 (prog prigi 0. 0. prigi);
  313. AMOR1 = MANU 'RIGI' P2 LIM1 (prog pamor 0. 0. pamor);
  314.  
  315. *------- on définit une unique table pour tout (HBM et CONTINU) -------
  316. TAB1 = TABLE ;
  317.  
  318. *--------- remplissage des matrices "temporelles" ---------
  319. * TAB1 . 'MASSE_CONSTANTE' = MASS1 ;
  320. * TAB1 . 'AMORTISSEMENT_CONSTANT' = AMOR1 ;
  321. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici
  322. TAB1 . 'MASSE_CONSTANTE' = (2.*pi **2) * MASS1 ;
  323. TAB1 . 'AMORTISSEMENT_CONSTANT' = (2.*pi) * AMOR1 ;
  324. TAB1 . 'RIGIDITE_CONSTANTE' = RI1 ;
  325. TAB1 . 'N_HARMONIQUE' = nhbm;
  326.  
  327. *------- remplissage des resultats attendus sur ddls "temporels" -------
  328. mycoul6 = MOTS 'VIOL' 'BLEU' 'BLEU' 'TURQ' 'TURQ' 'VERT' 'VERT'
  329. 'OLIV' 'OLIV' 'JAUN' 'JAUN' 'ORAN' 'ORAN'
  330. 'ROUG' 'ROUG' 'ROSE' 'ROSE' ;
  331. mycoul6 = mycoul6 et (MOTS 24*'GRIS');
  332. TAB1 . 'RESULTATS' = tabl;
  333. TAB1 . 'RESULTATS' . 1 = tabl;
  334. TAB1 . 'RESULTATS' . 1 . 'POINT_MESURE' = P2;
  335. TAB1 . 'RESULTATS' . 1 . 'COMPOSANTE' = mot 'UX';
  336. TAB1 . 'RESULTATS' . 1 . 'TITRE' = mot 'UX';
  337. TAB1 . 'RESULTATS' . 1 . 'COULEUR' = mot 'BLEU';
  338. TAB1 . 'RESULTATS' . 1 . 'COULEUR_HBM' = mycoul6;
  339. TAB1 . 'RESULTATS' . 2 = tabl;
  340. TAB1 . 'RESULTATS' . 2 . 'POINT_MESURE' = P2;
  341. TAB1 . 'RESULTATS' . 2 . 'COMPOSANTE' = mot 'UY';
  342. TAB1 . 'RESULTATS' . 2 . 'TITRE' = mot 'UY';
  343. TAB1 . 'RESULTATS' . 2 . 'COULEUR' = mot 'ROSE';
  344. TAB1 . 'RESULTATS' . 2 . 'COULEUR_HBM' = mycoul6;
  345.  
  346.  
  347. *--------- passage en frequentiel ---------
  348. HBM TAB1;
  349.  
  350.  
  351. MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP TEMPORELLES D UN NOEUD';
  352. LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT';
  353. MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP FREQUENTIELLES D UN NOEUD';
  354. LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT_HBM';
  355.  
  356.  
  357. *----------- definition du chargement -----------------------------------
  358. *
  359. * ici, on definit directement le chargement en frequentiel
  360. FP11 = MANU 'CHPO' p2 2 'F3' pbal 'G4' (pbal) 'NATURE' 'DISCRET';
  361. LIX1 = PROG 0. pas 0.1 (t1+30.) ;
  362. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici
  363. LIY1 = (2.*pi*LIX1)**2 ;
  364. EV1 = EVOL MANU 't' LIX1 'F(t)' LIY1 ;
  365. CHA1 = CHAR 'MECA' FP11 EV1 ;
  366. si GRAPH;
  367. dess EV1 POSX CENT POSY CENT ;
  368. finsi;
  369.  
  370.  
  371. *----------- resolution par la procedure CONTINU -----------------------
  372.  
  373. * OPTIONS DE CALCULS
  374.  
  375. TAB1 . 'CHARGEMENT' = CHA1;
  376. * calcul des forces NL
  377. TAB1 . 'PROCEDURE_CHAR_MECA'= VRAI ;
  378. TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT';
  379. TAB1 . 'N_PT_TFR' = OTFR;
  380. TAB1 . 'CALC_CHPO' = VRAI;
  381. * TAB1 . 'CALC_CHPO' = faux; -> charmeca pas prevu !
  382. *
  383. *plage de frequence en Hz
  384. LIS1T = prog t0 PAS dt t1;
  385. TAB1 . 'TEMPS_CALCULES' = LIS1T;
  386. TAB1 . 'MAXIPAS' = 500;
  387. TAB1 . 'MAXIPAS' = 300;
  388. TAB1 . 'PRECISION' = 1.E-5;
  389. TAB1 . 'ACCELERATION' = 4;
  390.  
  391. *------parametres du contact
  392. TAB1 . 'COEFF_FNL' = PROG pmu pkc ph0;
  393. TAB1 . 'POINT_FNL' = P2;
  394. TAB1 . 'COMP_FNL' = mots 'UX' 'UY';
  395.  
  396. *------parametres du continuation
  397. TAB1 . 'NB_ITERATION' = itmoy;
  398. TAB1 . 'MAXITERATION' = 30;
  399. TAB1 . 'WTABLE' = tabl;
  400. TAB1 . 'WTABLE' . 'DS' = ds;
  401. TAB1 . 'WTABLE' . 'DSMAX' = dsmax;
  402. TAB1 . 'WTABLE' . 'DSMIN' = dsmin;
  403. TAB1 . 'MAXI_DEPLACEMENT' = u1max;
  404. TAB1 . 'MAXSOUSPAS' = 3;
  405. * *resultats a sortir
  406. * TAB1 . 'RESULTATS' = TRES1;
  407.  
  408. * stabilité
  409. TAB1 . 'STABILITE' = vrai;
  410. *
  411.  
  412. * opti debug 1 ;
  413. *calcul avec continuation
  414. CONTINU TAB1;
  415.  
  416. TEMP 'IMPR' 'MAXI' 'CPU';
  417.  
  418. * *---------------------- sauv avant reprise ----------------------------
  419. * ficxdr = chai prefix '-avantREPRISE.xdr';
  420. * mess ficxdr;
  421. * opti sauv ficxdr;
  422. * sauv ;
  423.  
  424. *------------------------- reprise -----------------------------
  425.  
  426. * on impose l'absence de contact
  427. TAB1 . 'COEFF_FNL' = PROG (0.*pmu) (0.*pkc) (1E6*ph0);
  428.  
  429. * on fait croire qu'on vient dans la direction opposee de maniere
  430. * a etre capable de traiter le point de rebroussement
  431. WTAB1 = TAB1 . 'WTABLE';
  432. * OUBL TAB1 'WTABLE';
  433. list WTAB1;
  434. WTAB2 = TABL ;
  435. WTAB2 . 'DTIME0' = -1. * WTAB1 . 'DTIME0';
  436. WTAB2 . 'DDEP0' = -1. * WTAB1 . 'DDEP0';
  437. TAB1 . 'WTABLE' = WTAB2;
  438.  
  439. CONTINU TAB1;
  440.  
  441.  
  442. *------------------------- post-traitement -----------------------------
  443.  
  444.  
  445. ********** tracé de toutes les harmoniques individuellement **********
  446.  
  447. wprog = TAB1 . 'TEMPS_PROG';
  448. * wprog = (extr evtot cour 1) extr ABSC;
  449. evtot = TAB1 . 'RESULTATS_HBM' . 'RESULTATS_EVOL';
  450.  
  451. TDESS1 = TABL;
  452. TDESS1 . 'TITRE' = TAB1 . 'RESULTATS_HBM' . 'TITRE';
  453.  
  454. si GRAPH;
  455. dess evtot
  456. TITX '\w(Hz)' POSX CENT
  457. TITY 'U_{k}' POSY CENT TDESS1 'LEGE';
  458. dess evtot YBOR -0.5 0.3 XBOR 0. 12
  459. TITX '\w(Hz)' POSX CENT
  460. TITY 'U_{k}' POSY CENT TDESS1 'LEGE';
  461. finsi;
  462.  
  463.  
  464. ********** traduction des resultats frequentiels en temporels **********
  465.  
  466. HBM_POST TAB1;
  467.  
  468.  
  469. ************* tracé du max d'amplitude max_t(u(t)) *********************
  470.  
  471. evuyamp = TAB1 . 'RESULTATS' . 'RESULTATS_EVOL';
  472. evuyamp1 = extr evuyamp COUR 1;
  473.  
  474. TDESS2 = TABL;
  475. TDESS2 . 2 = mot 'TIRR';
  476. TDESS2 . 'TITRE' = TAB1 . 'RESULTATS' . 'TITRE';
  477.  
  478. si GRAPH;
  479. dess evuyamp1
  480. TITX '\w(Hz)' POSX CENT
  481. TITY 'max|U(t)|' POSY CENT TDESS2 LEGE 'NO';
  482. dess evuyamp1
  483. TITX '\w(Hz)' POSX CENT
  484. TITY 'max|U(t)|' POSY CENT TDESS2 LEGE ;
  485. finsi;
  486.  
  487. ************* tracé du max d'amplitude avec la stabilite ! *************
  488.  
  489. istab = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'STABILITE';
  490. Tdess2 = TABL;
  491. Tdess2 . 'LIGNE_VARIABLE' = TABL;
  492. Tdess2 . 'LIGNE_VARIABLE' . 1 = istab;
  493. si GRAPH;
  494. dess evuyamp1
  495. TITX '\w(Hz)' POSX CENT
  496. TITY 'max|U(t)|' POSY CENT Tdess2;
  497. Tdess2 . 1 = mot 'MARQ S PLUS';
  498. dess evuyamp1
  499. TITX '\w(Hz)' POSX CENT
  500. TITY 'max|U(t)|' POSY CENT Tdess2;
  501. finsi;
  502.  
  503. * autres tracé de la stabilite : lambda_R en fonction de W
  504. mycoul = mots 'VIOL' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'OLIV'
  505. 'JAUN' 'ORAN' 'ROUG' 'ROSE';
  506. ncoul = dime mycoul;
  507. wprog2 = enle wprog (dime wprog);
  508. evlitot = VIDE 'EVOLUTIO' ;
  509. evlrtot = VIDE 'EVOLUTIO' ;
  510. nl = dime TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL';
  511. il = 0;
  512. icoul = 0;
  513. repe bl nl; il = il + 1; icoul = icoul + 1;
  514. si(icoul > ncoul); icoul = 1; finsi; coul1 = extr mycoul icoul;
  515. lr1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL' . il;
  516. evlrtot = evlrtot
  517. et (evol coul1 MANU '\w (Hz)' wprog2 '\l_{R}' lr1);
  518. li1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_IMAG' . il;
  519. evlitot = evlitot
  520. et (evol coul1 MANU '\w (Hz)' wprog2 '\l_{I}' li1);
  521. fin bl;
  522. wprog3 = prog (mini wprog2) (maxi wprog2);
  523. evzero = evol manu '\w (Hz)' wprog3 '\l_{R}' (0.*wprog3);
  524. Tdess3 = tabl;
  525. Tdess3 . 1 = mot 'TIRR';
  526.  
  527. si GRAPH;
  528. dess (evzero et evlrtot) 'TITY' '\l_{R} (s^{-1})' Tdess3;
  529. dess (evlitot) 'TITY' '\l_{I} (Hz)' ;
  530. finsi;
  531.  
  532. * stabilite : trace dynamique pour faire une animation
  533. si GRAPH;
  534. ficps = chai prefix '-ORBITE.ps';
  535. opti 'FTRA' ficps;
  536. TABORB = tabl;
  537. TABORB . 'PAS' = 3;
  538. TABORB . 'QUEUE' = MOT 'INFINIE';
  539. TABORB . 'EVOL_FIXE' = evzero ;
  540. ORBITE evlrtot TABORB;
  541. finsi;
  542.  
  543.  
  544.  
  545. *------------------------- sauvegarde -----------------------------
  546.  
  547. si SAVE ;
  548. ficxdr = chai prefix '.xdr';
  549. mess ficxdr;
  550. opti sauv ficxdr;
  551. sauv ;
  552. finsi;
  553.  
  554.  
  555.  
  556. *--------------------- test de non regression ---------------------
  557.  
  558. si TNR ;
  559. * reference : these de Lihan Xie fig3.10 page 84
  560. * sauf qu'ici, on ne prend que 3 harmoniques
  561.  
  562. * recup valeurs de la courbes de reponse
  563. wamp = extr evuyamp 'ABSC' 1;
  564. wadim = wamp / wc;
  565. uyamp = extr evuyamp 'ORDO' 1;
  566. ramp = uyamp;
  567. * R=UY puisque l'on a des orbites circulaires,
  568. * sinon il faudrait calculer r(t)=uy(t)**2+uz(t)**2 , puis R=max(r(t))
  569. namp = DIME wamp;
  570.  
  571. * recherche de la mise en contact
  572. idebut = POSI 1 'DANS' (ENTI (MASQ ramp EGSUPE ph0));
  573. wdeb0 = extr wadim (idebut - 1);
  574. wdeb1 = extr wadim idebut;
  575. si ((wdeb0 > 0.16) ou (wdeb1 < 0.16));
  576. mess 'Mise en contact devrait se produire vers w/w0~0.16 !';
  577. mess 'PB car : ' wdeb0 ' < 0.16 < 'wdeb1 ' non verifie !';
  578. ERRE 5;
  579. finsi;
  580.  
  581. * verification que le calcul a bien attenit les w=10Hz
  582. wmax = maxi wamp;
  583. si (wmax < 10.);
  584. mess 'On devrait calculer jusqu a 10Hz au moins !';
  585. mess wmax;
  586. ERRE 5;
  587. finsi;
  588.  
  589. * verification de l'amplitude max atteinte
  590. * rem : rmax = 4.4807 lors de la creation du cas test,
  591. * mais on arrondit pour rref
  592. rmax = (maxi uyamp) / ph0;
  593. rref = 4.5 ;
  594. si ((ABS (rmax - rref)) > 0.05);
  595. mess 'Le max|U| devrait etre de ' rref ' !';
  596. mess rmax;
  597. ERRE 5;
  598. finsi;
  599.  
  600. * recherche des changements de stabilite
  601. istab1 = (istab enle 1) - (istab enle (DIME istab));
  602. istab1 = (posi 1 'DANS' istab1 'TOUS')
  603. et (posi -1 'DANS' istab1 'TOUS');
  604. si (neg (dime istab1) 4);
  605. mess 'Il devrait y avoir 4 points de changement de stabilite !';
  606. list istab1;
  607. ERRE 5;
  608. finsi;
  609.  
  610. finsi;
  611.  
  612.  
  613.  
  614. fin ;
  615.  
  616.  
  617.  

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