Télécharger la plaquette Cast3M


fr en

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.  
  231. GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI;
  232. * GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX;
  233. *
  234. *-------- Definition des options du calcul :
  235. * complet = faux;
  236. * CALCUL PAR DFT ?  (Transformation de Fourier Discret
  237. CALC_DFT = FAUX;
  238. *si analyser la stabilite?
  239. OP_STAB = FAUX;
  240. * *si adimensionnement du temps?
  241. * Adim_t = VRAI;
  242. *
  243. *-------- Definition des parametres du calcul --------
  244. *nombre d'harmoniques
  245.  nhbm = 3 ;
  246. *nhbm = 7 ;
  247. * nhbm = 9 ;
  248. * nhbm = 11 ;
  249. *nhbm = 15 ;
  250. *plage d'etude de frequence en Hz
  251. t0 = 0.0;
  252. t1 = 10.; 
  253. dt = 0.10;
  254. *des iterations
  255. itmoy = 6;
  256. ds = 1.;
  257. dsmax = 1.;
  258. dsmin = 1.E-3 ;
  259. * si AFT,  OTFR: nb de points 2**OTFR pour TFR
  260. * OTFR = 7;
  261. OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2)); 
  262. mess 'OTFR = ' OTFR;
  263.  
  264. *-------- Definition des parametres du contact ---------------
  265. *coefficient de frottement
  266. pmu = 0.11 ; 
  267. si (pmu ega 0.00 1.E-4); mopmu = mot '000'; finsi;
  268. si (pmu ega 0.05 1.E-4); mopmu = mot '005'; finsi;
  269. si (pmu ega 0.11 1.E-4); mopmu = mot '011'; finsi;
  270. si (pmu ega 0.15 1.E-4); mopmu = mot '015'; finsi;
  271. si (pmu ega 0.20 1.E-4); mopmu = mot '020'; finsi;
  272. si (pmu ega 0.25 1.E-4); mopmu = mot '025'; finsi;
  273. *rigidite du contact
  274. pkc = 2500. ;
  275. * pkc = 0.;
  276. *jeu initial
  277. ph0 = 0.105 ;
  278. *
  279. *parametre de l'adimensionnement = 20% du jeu
  280. u1max = 0.20*ph0;
  281. *
  282. *---------Definition des parametres du rotor --------
  283. *masse constant
  284. pmass = 1. ;
  285. *amortissement
  286. pamor = 5. ;
  287. *rigidite
  288. prigi = 100. ;
  289. *force du balourd
  290. pbal = 0.1;
  291. * frequences caracteristiques du pb (en Hz)
  292. w0  = (prigi/pmass)**0.5  /  (2.*pi);
  293. w0c = ((prigi+pkc)/pmass)**0.5  /  (2.*pi); 
  294. wc  = ((pkc)/pmass)**0.5  /  (2.*pi); 
  295. mess w0 w0c wc;
  296.  
  297. * prefix :
  298. prefix = chai 'hbm_jeffcott_contact_alfa_' mopmu '_n' nhbm ;
  299. TITRE prefix;
  300. si GRAPH;
  301.   ficps = chai prefix '.ps';
  302.   opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps;
  303. finsi;
  304.  
  305. *------------- geometrie --------
  306. *
  307. P1 = 0.    0.    ; P2 = 1. 0.  ; 
  308.  
  309. *----------- construction des matrices caracteristiques ----------       
  310. *                                                                               
  311. *OPTI DONN 5;
  312. LIM1 = MOTS 'ALFA';
  313. MASS1 = MANU 'RIGI' 'TYPE' 'MASSE' (P1 et P2) LIM1 (prog pmass) ; 
  314. RI1   = MANU 'RIGI' (P1 et P2) LIM1 (prog prigi) ; 
  315. AMOR1 = MANU 'RIGI' 'TYPE' 'AMORTISSEMENT' (P1 et P2) LIM1 (prog pamor); 
  316.  
  317. *------- on définit une unique table pour tout (HBM et CONTINU) -------
  318. TAB1 = TABLE ;
  319.  
  320. *--------- remplissage des matrices "temporelles" ---------
  321. * TAB1 . 'MASSE_CONSTANTE'        =  MASS1 ;
  322. * TAB1 . 'AMORTISSEMENT_CONSTANT' =  AMOR1 ;
  323. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici
  324. TAB1 . 'MASSE_CONSTANTE'        = (2.*pi **2) * MASS1 ;
  325. TAB1 . 'AMORTISSEMENT_CONSTANT' = (2.*pi) * AMOR1 ;
  326. TAB1 . 'RIGIDITE_CONSTANTE'     = RI1 ; 
  327. TAB1 . 'N_HARMONIQUE'           = nhbm;
  328.  
  329. *------- remplissage des resultats attendus sur ddls "temporels" -------
  330. mycoul6 = MOTS 'VIOL' 'BLEU' 'BLEU' 'TURQ' 'TURQ' 'VERT' 'VERT' 
  331.                       'OLIV' 'OLIV' 'JAUN' 'JAUN' 'ORAN' 'ORAN'
  332.                       'ROUG' 'ROUG' 'ROSE' 'ROSE' ;
  333. mycoul6 =  mycoul6 et (MOTS 24*'GRIS');                     
  334. TAB1 . 'RESULTATS' = tabl; 
  335. TAB1 . 'RESULTATS' . 1 = tabl; 
  336. TAB1 . 'RESULTATS' . 1 . 'POINT_MESURE' = P1; 
  337. TAB1 . 'RESULTATS' . 1 . 'COMPOSANTE'  = mot 'ALFA'; 
  338. TAB1 . 'RESULTATS' . 1 . 'TITRE'       = mot 'UX'; 
  339. TAB1 . 'RESULTATS' . 1 . 'COULEUR'     = mot 'BLEU'; 
  340. TAB1 . 'RESULTATS' . 1 . 'COULEUR_HBM' = mycoul6; 
  341. TAB1 . 'RESULTATS' . 2 = tabl; 
  342. TAB1 . 'RESULTATS' . 2 . 'POINT_MESURE' = P2; 
  343. TAB1 . 'RESULTATS' . 2 . 'COMPOSANTE'  = mot 'ALFA'; 
  344. TAB1 . 'RESULTATS' . 2 . 'TITRE'       = mot 'UY'; 
  345. TAB1 . 'RESULTATS' . 2 . 'COULEUR'     = mot 'ROSE'; 
  346. TAB1 . 'RESULTATS' . 2 . 'COULEUR_HBM' = mycoul6; 
  347.  
  348.  
  349. *--------- passage en frequentiel ---------
  350.   HBM TAB1;
  351.  
  352.  
  353. MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP TEMPORELLES D UN NOEUD';
  354. LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT';
  355. MESS '>>>>>>>>>>>>>>>COMPOSANTES DEP FREQUENTIELLES D UN NOEUD';
  356. LIST TAB1 . 'COMPOSANTES' . 'DEPLACEMENT_HBM'; 
  357.  
  358.  
  359. *----------- definition du chargement -----------------------------------       
  360. *
  361. * ici, on definit directement le chargement en frequentiel                                                       
  362. FP11 = (MANU 'CHPO' p1 1 'F2' pbal 'NATURE' 'DISCRET')
  363.     et (MANU 'CHPO' p2 1 'G2' pbal 'NATURE' 'DISCRET');
  364. LIX1 = PROG 0. pas 0.1 (t1+30.) ;  
  365. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN Hz => mettre le 2pi ici
  366. LIY1 = (2.*pi*LIX1)**2 ;
  367. EV1 = EVOL MANU 't' LIX1 'F(t)' LIY1 ;  
  368. CHA1 = CHAR 'MECA' FP11 EV1 ; 
  369. si GRAPH;
  370.   dess EV1 POSX CENT POSY CENT ;
  371. finsi;
  372.  
  373.  
  374. *----------- resolution par la procedure CONTINU ----------------------- 
  375.  
  376. * OPTIONS DE CALCULS
  377.  
  378. TAB1 . 'CHARGEMENT' = CHA1;
  379. * calcul des forces NL
  380. TAB1 . 'PROCEDURE_CHARMECA'= VRAI ;
  381. TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT';
  382. TAB1 . 'N_PT_TFR' = OTFR;
  383. TAB1 . 'CALC_CHPO' = VRAI;
  384. * TAB1 . 'CALC_CHPO' = faux; -> charmeca pas prevu !
  385. *
  386. *plage de frequence en Hz
  387. LIS1T = prog t0 PAS dt t1; 
  388. TAB1 . 'TEMPS_CALCULES' = LIS1T;
  389. TAB1 . 'MAXIPAS' = 300;
  390. TAB1 . 'PRECISION' = 1.E-5;
  391. TAB1 . 'ACCELERATION' = 4;
  392.  
  393. *------parametres du contact
  394. TAB1 . 'COEFF_FNL' = PROG pmu pkc ph0;
  395. TAB1 . 'POINT_FNL' = (P1 et P2);
  396. TAB1 . 'COMP_FNL' = mots 'ALFA';
  397.  
  398. *------parametres du continuation
  399. TAB1 . 'NB_ITERATION' = itmoy;
  400. TAB1 . 'MAXITERATION' = 30;
  401. TAB1 . 'WTABLE' = tabl;
  402. TAB1 . 'WTABLE' . 'DS' = ds;
  403. TAB1 . 'WTABLE' . 'DSMAX' = dsmax;
  404. TAB1 . 'WTABLE' . 'DSMIN' = dsmin;
  405. TAB1 . 'MAXI_DEPLACEMENT' = u1max;
  406. TAB1 . 'MAXSOUSPAS' = 3;
  407.  
  408. * stabilité 
  409. TAB1 . 'STABILITE'  = vrai;
  410.  
  411.  
  412. * 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.  
  614.  
  615.  
  616.  
  617.  

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