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.  
  226. GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI;
  227. * GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX;
  228.  
  229. *
  230. *-------- Definition des options du calcul :
  231. * complet = faux;
  232. * CALCUL PAR DFT ? (Transformation de Fourier Discret
  233. CALC_DFT = FAUX;
  234. *si analyser la stabilite?
  235. OP_STAB = FAUX;
  236. * *si adimensionnement du temps?
  237. * Adim_t = VRAI;
  238. *
  239. *-------- Definition des parametres du calcul --------
  240. *nombre d'harmoniques
  241. nhbm = 3 ;
  242. *nhbm = 7 ;
  243. * nhbm = 9 ;
  244. * nhbm = 11 ;
  245. *nhbm = 15 ;
  246. *plage d'etude de frequence en Hz
  247. t0 = 0.0;
  248. t1 = 10.;
  249. dt = 0.10;
  250. *des iterations
  251. itmoy = 6;
  252. ds = 1.;
  253. dsmax = 1.;
  254. dsmin = 1.E-3 ;
  255. * si AFT, OTFR: nb de points 2**OTFR pour TFR
  256. * OTFR = 7;
  257. OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2));
  258. mess 'OTFR = ' OTFR;
  259.  
  260. *-------- Definition des parametres du contact ---------------
  261. *coefficient de frottement
  262. pmu = 0.11 ;
  263. si (pmu ega 0.00 1.E-4); mopmu = mot '000'; finsi;
  264. si (pmu ega 0.05 1.E-4); mopmu = mot '005'; finsi;
  265. si (pmu ega 0.11 1.E-4); mopmu = mot '011'; finsi;
  266. si (pmu ega 0.15 1.E-4); mopmu = mot '015'; finsi;
  267. si (pmu ega 0.20 1.E-4); mopmu = mot '020'; finsi;
  268. si (pmu ega 0.25 1.E-4); mopmu = mot '025'; finsi;
  269. *rigidite du contact
  270. pkc = 2500. ;
  271. * pkc = 0.;
  272. *jeu initial
  273. ph0 = 0.105 ;
  274. *
  275. *parametre de l'adimensionnement = 20% du jeu
  276. u1max = 0.20*ph0;
  277. *
  278. *---------Definition des parametres du rotor --------
  279. *masse constant
  280. pmass = 1. ;
  281. *amortissement
  282. pamor = 5. ;
  283. *rigidite
  284. prigi = 100. ;
  285. *force du balourd
  286. pbal = 0.1;
  287.  
  288. * frequences caracteristiques du pb (en Hz)
  289. w0 = (prigi/pmass)**0.5 / (2.*pi);
  290. w0c = ((prigi+pkc)/pmass)**0.5 / (2.*pi);
  291. wc = ((pkc)/pmass)**0.5 / (2.*pi);
  292. mess w0 w0c wc;
  293.  
  294. * prefix :
  295. prefix = chai 'hbm_jeffcott_contact_' mopmu '_n' nhbm ;
  296. TITRE prefix;
  297. si GRAPH;
  298. ficps = chai prefix '.ps';
  299. opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps;
  300. finsi;
  301.  
  302. *------------- geometrie --------
  303. *
  304. P1 = 0. 0. ; P2 = 1. 0. ;
  305.  
  306. *----------- construction des matrices caracteristiques ----------
  307. *
  308. *OPTI DONN 5;
  309. LIM1 = MOTS UX UY;
  310. MASS1 = (MASS 'UX' pmass P2 ) ET (MASS 'UY' pmass P2 );
  311. RI1 = MANU 'RIGI' P2 LIM1 (prog prigi 0. 0. prigi);
  312. AMOR1 = MANU 'RIGI' P2 LIM1 (prog pamor 0. 0. pamor);
  313.  
  314. *------- on définit une unique table pour tout (HBM et CONTINU) -------
  315. TAB1 = TABLE ;
  316.  
  317. *--------- remplissage des matrices "temporelles" ---------
  318. * TAB1 . 'MASSE_CONSTANTE' = MASS1 ;
  319. * TAB1 . 'AMORTISSEMENT_CONSTANT' = AMOR1 ;
  320. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici
  321. TAB1 . 'MASSE_CONSTANTE' = (2.*pi **2) * MASS1 ;
  322. TAB1 . 'AMORTISSEMENT_CONSTANT' = (2.*pi) * AMOR1 ;
  323. TAB1 . 'RIGIDITE_CONSTANTE' = RI1 ;
  324. TAB1 . 'N_HARMONIQUE' = nhbm;
  325.  
  326. *------- remplissage des resultats attendus sur ddls "temporels" -------
  327. mycoul6 = MOTS 'VIOL' 'BLEU' 'BLEU' 'TURQ' 'TURQ' 'VERT' 'VERT'
  328. 'OLIV' 'OLIV' 'JAUN' 'JAUN' 'ORAN' 'ORAN'
  329. 'ROUG' 'ROUG' 'ROSE' 'ROSE' ;
  330. mycoul6 = mycoul6 et (MOTS 24*'GRIS');
  331. TAB1 . 'RESULTATS' = tabl;
  332. TAB1 . 'RESULTATS' . 1 = tabl;
  333. TAB1 . 'RESULTATS' . 1 . 'POINT_MESURE' = P2;
  334. TAB1 . 'RESULTATS' . 1 . 'COMPOSANTE' = mot 'UX';
  335. TAB1 . 'RESULTATS' . 1 . 'TITRE' = mot 'UX';
  336. TAB1 . 'RESULTATS' . 1 . 'COULEUR' = mot 'BLEU';
  337. TAB1 . 'RESULTATS' . 1 . 'COULEUR_HBM' = mycoul6;
  338. TAB1 . 'RESULTATS' . 2 = tabl;
  339. TAB1 . 'RESULTATS' . 2 . 'POINT_MESURE' = P2;
  340. TAB1 . 'RESULTATS' . 2 . 'COMPOSANTE' = mot 'UY';
  341. TAB1 . 'RESULTATS' . 2 . 'TITRE' = mot 'UY';
  342. TAB1 . 'RESULTATS' . 2 . 'COULEUR' = mot 'ROSE';
  343. TAB1 . 'RESULTATS' . 2 . 'COULEUR_HBM' = mycoul6;
  344.  
  345.  
  346. *--------- passage en frequentiel ---------
  347. HBM TAB1;
  348.  
  349.  
  350. MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP TEMPORELLES D UN NOEUD';
  351. LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT';
  352. MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP FREQUENTIELLES D UN NOEUD';
  353. LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT_HBM';
  354.  
  355.  
  356. *----------- definition du chargement -----------------------------------
  357. *
  358. * ici, on definit directement le chargement en frequentiel
  359. FP11 = MANU 'CHPO' p2 2 'F3' pbal 'G4' (pbal) 'NATURE' 'DISCRET';
  360. LIX1 = PROG 0. pas 0.1 (t1+30.) ;
  361. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici
  362. LIY1 = (2.*pi*LIX1)**2 ;
  363. EV1 = EVOL MANU 't' LIX1 'F(t)' LIY1 ;
  364. CHA1 = CHAR 'MECA' FP11 EV1 ;
  365. si GRAPH;
  366. dess EV1 POSX CENT POSY CENT ;
  367. finsi;
  368.  
  369.  
  370. *----------- resolution par la procedure CONTINU -----------------------
  371.  
  372. * OPTIONS DE CALCULS
  373.  
  374. TAB1 . 'CHARGEMENT' = CHA1;
  375. * calcul des forces NL
  376. TAB1 . 'PROCEDURE_CHARMECA'= VRAI ;
  377. TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT';
  378. TAB1 . 'N_PT_TFR' = OTFR;
  379. TAB1 . 'CALC_CHPO' = VRAI;
  380. * TAB1 . 'CALC_CHPO' = faux; -> charmeca pas prevu !
  381. *
  382. *plage de frequence en Hz
  383. LIS1T = prog t0 PAS dt t1;
  384. TAB1 . 'TEMPS_CALCULES' = LIS1T;
  385. TAB1 . 'MAXIPAS' = 500;
  386. TAB1 . 'MAXIPAS' = 300;
  387. TAB1 . 'PRECISION' = 1.E-5;
  388. TAB1 . 'ACCELERATION' = 4;
  389.  
  390. *------parametres du contact
  391. TAB1 . 'COEFF_FNL' = PROG pmu pkc ph0;
  392. TAB1 . 'POINT_FNL' = P2;
  393. TAB1 . 'COMP_FNL' = mots 'UX' 'UY';
  394.  
  395. *------parametres du continuation
  396. TAB1 . 'NB_ITERATION' = itmoy;
  397. TAB1 . 'MAXITERATION' = 30;
  398. TAB1 . 'WTABLE' = tabl;
  399. TAB1 . 'WTABLE' . 'DS' = ds;
  400. TAB1 . 'WTABLE' . 'DSMAX' = dsmax;
  401. TAB1 . 'WTABLE' . 'DSMIN' = dsmin;
  402. TAB1 . 'MAXI_DEPLACEMENT' = u1max;
  403. TAB1 . 'MAXSOUSPAS' = 3;
  404. * *resultats a sortir
  405. * TAB1 . 'RESULTATS' = TRES1;
  406.  
  407. * stabilité
  408. TAB1 . 'STABILITE' = vrai;
  409. *
  410.  
  411. * opti debug 1 ;
  412. *calcul avec continuation
  413. CONTINU TAB1;
  414.  
  415. TEMP 'IMPR' 'MAXI' 'CPU';
  416.  
  417. * *---------------------- sauv avant reprise ----------------------------
  418. * ficxdr = chai prefix '-avantREPRISE.xdr';
  419. * mess ficxdr;
  420. * opti sauv ficxdr;
  421. * sauv ;
  422.  
  423. *------------------------- reprise -----------------------------
  424.  
  425. * on impose l'absence de contact
  426. TAB1 . 'COEFF_FNL' = PROG (0.*pmu) (0.*pkc) (1E6*ph0);
  427.  
  428. * on fait croire qu'on vient dans la direction opposee de maniere
  429. * a etre capable de traiter le point de rebroussement
  430. WTAB1 = TAB1 . 'WTABLE';
  431. * OUBL TAB1 'WTABLE';
  432. list WTAB1;
  433. WTAB2 = TABL ;
  434. WTAB2 . 'DTIME0' = -1. * WTAB1 . 'DTIME0';
  435. WTAB2 . 'DDEP0' = -1. * WTAB1 . 'DDEP0';
  436. TAB1 . 'WTABLE' = WTAB2;
  437.  
  438. CONTINU TAB1;
  439.  
  440.  
  441. *------------------------- post-traitement -----------------------------
  442.  
  443.  
  444. ********** tracé de toutes les harmoniques individuellement **********
  445.  
  446. wprog = TAB1 . 'TEMPS_PROG';
  447. * wprog = (extr evtot cour 1) extr ABSC;
  448. evtot = TAB1 . 'RESULTATS_HBM' . 'RESULTATS_EVOL';
  449.  
  450. TDESS1 = TABL;
  451. TDESS1 . 'TITRE' = TAB1 . 'RESULTATS_HBM' . 'TITRE';
  452.  
  453. si GRAPH;
  454. dess evtot
  455. TITX '\w(Hz)' POSX CENT
  456. TITY 'U_{k}' POSY CENT TDESS1 'LEGE';
  457. dess evtot YBOR -0.5 0.3 XBOR 0. 12
  458. TITX '\w(Hz)' POSX CENT
  459. TITY 'U_{k}' POSY CENT TDESS1 'LEGE';
  460. finsi;
  461.  
  462.  
  463. ********** traduction des resultats frequentiels en temporels **********
  464.  
  465. HBM_POST TAB1;
  466.  
  467.  
  468. ************* tracé du max d'amplitude max_t(u(t)) *********************
  469.  
  470. evuyamp = TAB1 . 'RESULTATS' . 'RESULTATS_EVOL';
  471. evuyamp1 = extr evuyamp COUR 1;
  472.  
  473. TDESS2 = TABL;
  474. TDESS2 . 2 = mot 'TIRR';
  475. TDESS2 . 'TITRE' = TAB1 . 'RESULTATS' . 'TITRE';
  476.  
  477. si GRAPH;
  478. dess evuyamp1
  479. TITX '\w(Hz)' POSX CENT
  480. TITY 'max|U(t)|' POSY CENT TDESS2 LEGE 'NO';
  481. dess evuyamp1
  482. TITX '\w(Hz)' POSX CENT
  483. TITY 'max|U(t)|' POSY CENT TDESS2 LEGE ;
  484. finsi;
  485.  
  486. ************* tracé du max d'amplitude avec la stabilite ! *************
  487.  
  488. istab = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'STABILITE';
  489. Tdess2 = TABL;
  490. Tdess2 . 'LIGNE_VARIABLE' = TABL;
  491. Tdess2 . 'LIGNE_VARIABLE' . 1 = istab;
  492. si GRAPH;
  493. dess evuyamp1
  494. TITX '\w(Hz)' POSX CENT
  495. TITY 'max|U(t)|' POSY CENT Tdess2;
  496. Tdess2 . 1 = mot 'MARQ S PLUS';
  497. dess evuyamp1
  498. TITX '\w(Hz)' POSX CENT
  499. TITY 'max|U(t)|' POSY CENT Tdess2;
  500. finsi;
  501.  
  502. * autres tracé de la stabilite : lambda_R en fonction de W
  503. mycoul = mots 'VIOL' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'OLIV'
  504. 'JAUN' 'ORAN' 'ROUG' 'ROSE';
  505. ncoul = dime mycoul;
  506. wprog2 = enle wprog (dime wprog);
  507. evlitot = VIDE 'EVOLUTIO' ;
  508. evlrtot = VIDE 'EVOLUTIO' ;
  509. nl = dime TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL';
  510. il = 0;
  511. icoul = 0;
  512. repe bl nl; il = il + 1; icoul = icoul + 1;
  513. si(icoul > ncoul); icoul = 1; finsi; coul1 = extr mycoul icoul;
  514. lr1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL' . il;
  515. evlrtot = evlrtot
  516. et (evol coul1 MANU '\w (Hz)' wprog2 '\l_{R}' lr1);
  517. li1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_IMAG' . il;
  518. evlitot = evlitot
  519. et (evol coul1 MANU '\w (Hz)' wprog2 '\l_{I}' li1);
  520. fin bl;
  521. wprog3 = prog (mini wprog2) (maxi wprog2);
  522. evzero = evol manu '\w (Hz)' wprog3 '\l_{R}' (0.*wprog3);
  523. Tdess3 = tabl;
  524. Tdess3 . 1 = mot 'TIRR';
  525.  
  526. si GRAPH;
  527. dess (evzero et evlrtot) 'TITY' '\l_{R} (s^{-1})' Tdess3;
  528. dess (evlitot) 'TITY' '\l_{I} (Hz)' ;
  529. finsi;
  530.  
  531. * stabilite : trace dynamique pour faire une animation
  532. si GRAPH;
  533. ficps = chai prefix '-ORBITE.ps';
  534. opti 'FTRA' ficps;
  535. TABORB = tabl;
  536. TABORB . 'PAS' = 3;
  537. TABORB . 'QUEUE' = MOT 'INFINIE';
  538. TABORB . 'EVOL_FIXE' = evzero ;
  539. ORBITE evlrtot TABORB;
  540. finsi;
  541.  
  542.  
  543.  
  544. *------------------------- sauvegarde -----------------------------
  545.  
  546. si SAVE ;
  547. ficxdr = chai prefix '.xdr';
  548. mess ficxdr;
  549. opti sauv ficxdr;
  550. sauv ;
  551. finsi;
  552.  
  553.  
  554.  
  555. *--------------------- test de non regression ---------------------
  556.  
  557. si TNR ;
  558. * reference : these de Lihan Xie fig3.10 page 84
  559. * sauf qu'ici, on ne prend que 3 harmoniques
  560.  
  561. * recup valeurs de la courbes de reponse
  562. wamp = extr evuyamp 'ABSC' 1;
  563. wadim = wamp / wc;
  564. uyamp = extr evuyamp 'ORDO' 1;
  565. ramp = uyamp;
  566. * R=UY puisque l'on a des orbites circulaires,
  567. * sinon il faudrait calculer r(t)=uy(t)**2+uz(t)**2 , puis R=max(r(t))
  568. namp = DIME wamp;
  569.  
  570. * recherche de la mise en contact
  571. idebut = POSI 1 'DANS' (ENTI (MASQ ramp EGSUPE ph0));
  572. wdeb0 = extr wadim (idebut - 1);
  573. wdeb1 = extr wadim idebut;
  574. si ((wdeb0 > 0.16) ou (wdeb1 < 0.16));
  575. mess 'Mise en contact devrait se produire vers w/w0~0.16 !';
  576. mess 'PB car : ' wdeb0 ' < 0.16 < 'wdeb1 ' non verifie !';
  577. ERRE 5;
  578. finsi;
  579.  
  580. * verification que le calcul a bien attenit les w=10Hz
  581. wmax = maxi wamp;
  582. si (wmax < 10.);
  583. mess 'On devrait calculer jusqu a 10Hz au moins !';
  584. mess wmax;
  585. ERRE 5;
  586. finsi;
  587.  
  588. * verification de l'amplitude max atteinte
  589. * rem : rmax = 4.4807 lors de la creation du cas test,
  590. * mais on arrondit pour rref
  591. rmax = (maxi uyamp) / ph0;
  592. rref = 4.5 ;
  593. si ((ABS (rmax - rref)) > 0.05);
  594. mess 'Le max|U| devrait etre de ' rref ' !';
  595. mess rmax;
  596. ERRE 5;
  597. finsi;
  598.  
  599. * recherche des changements de stabilite
  600. istab1 = (istab enle 1) - (istab enle (DIME istab));
  601. istab1 = (posi 1 'DANS' istab1 'TOUS')
  602. et (posi -1 'DANS' istab1 'TOUS');
  603. si (neg (dime istab1) 4);
  604. mess 'Il devrait y avoir 4 points de changement de stabilite !';
  605. list istab1;
  606. ERRE 5;
  607. finsi;
  608.  
  609. finsi;
  610.  
  611.  
  612.  
  613. fin ;
  614.  
  615.  
  616.  
  617.  
  618.  
  619.  
  620.  

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