Télécharger hbm_jeffcott_contact_alfa.dgibi

Retour à la liste

Numérotation des lignes :

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

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