Télécharger la plaquette Cast3M


fr en

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 ;
  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.  
  618.  
  619.  
  620.  
  621.  

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