Télécharger hbm_duffing_mu.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : hbm_duffing_mu.dgibi
  2. ******************************************************************
  3. *
  4. * Systeme a un degre de liberte
  5. *
  6. * Oscillateur de Duffing calcule en k harmoniques :
  7. *
  8. * m*x'' + c*x' + k*x + a*x^3 = f*cos(w*t)
  9. *
  10. * Continuation + HBM + AFT
  11. *
  12. ******************************************************************
  13.  
  14. *-----------------------------------------------------------------------
  15. * CHARMECA.PROCEDUR
  16. * Calcul de la force non lineaire fnl(x,t) et de sa derivee dfnl(x,t)/dx
  17. * dans le domaine temporel
  18. *-----------------------------------------------------------------------
  19. *
  20. * EN ENTREE :
  21. *
  22. * TAB1
  23. * >> parametres génériques :
  24. * . 'DEP_NL' : LISTCHPO DE DEPLACEMENTS TEMPORELLES
  25. * ou (pas encore programmé)
  26. * . i . 'CHPOINT' = CHPOINT unitaire des deplacements NL
  27. * . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
  28. * >> parametres propres a ce cas test :
  29. * . 'COEFF_FNL' : coefficient de non linearite
  30. * . 'POINT_FNL' : point de non linearite
  31. * . 'COMP_FNL' : composante de non linearite
  32. *
  33. *
  34. * EN SORTIE :
  35. *
  36. * TAB1
  37. * >> parametres génériques :
  38. * . 'FNL_T' : LISTCHPO DES FORCES NON LINEAIRES TEMPORELLES
  39. * ou TABLE d'indice : (pas encore programmé)
  40. * . i . 'CHPOINT' = CHPOINT unitaire des forces non lineaires
  41. * . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
  42. * . 'KNL_T' : DERIVEE DE FORCES NON LINEAIRES sous forme de TABLE :
  43. * . i . j = 'LISTREEL' evolution temporelle de dFi/dxj
  44. *
  45. *-----------------------------------------------------------------------
  46. DEBPROC CHARMECA TAB1*'TABLE';
  47.  
  48. ******** parametres d'entree propres a ce cas test **********
  49.  
  50. * parametre de continuation
  51. time = TAB1 . 'TIME';
  52.  
  53. *coefficient de la force nl
  54. SI (exis TAB1 . 'COEFF_FNL');
  55. ALev = TAB1 . 'COEFF_FNL';
  56. al = IPOL ALev time;
  57. SINON;
  58. AL = 0.;
  59. MESS 'CHARMECA: PAS D EFFET NON LINEAIRE !';
  60. FINSI;
  61.  
  62. SI (EXIS TAB1 . 'POINT_FNL');
  63. pNL = TAB1 . 'POINT_FNL';
  64. SINON;
  65. MESS 'CHARMECA: IL MANQUE LA GEOMETRIE NL dans POINT_FNL'; ERRE 21;
  66. FINSI;
  67.  
  68. SI (EXIS TAB1 . 'COMP_FNL');
  69. COMP_FNL = TAB1 . 'COMP_FNL';
  70. n_CNL = DIME COMP_FNL;
  71. SINON;
  72. MESS 'CHARMECA: IL MANQUE LA COMPOSANTE NL dans COMP_FNL'; ERRE 21;
  73. FINSI;
  74.  
  75. ******** parametres d'entree génériques **********
  76.  
  77. *DEPNL : list de champ par points, deplacements temporels du point nl
  78. *
  79. SI (EXIS TAB1 'DEP_NL');
  80. DEPNL = TAB1 . 'DEP_NL';
  81. IFCHPO = ega (type DEPNL) 'LISTCHPO';
  82. si (non IFCHPO);
  83. si (neg (type DEPNL) 'TABLE');
  84. mess 'CHARMECA: DEP_NL doit etre de type LISTCHPO ou TABLE';
  85. ERRE 21;
  86. sinon;
  87. si (neg (type DEPNL . 1) 'LISTREEL');
  88. mess 'CHARMECA: DEP_NL . i doit etre un LISTREEL';
  89. ERRE 21;
  90. finsi;
  91. finsi;
  92. finsi;
  93. SINON;
  94. MESS 'CHARMECA: IL MANQUE LA LISTE DES DEPLACEMENTS dans DEP_NL';
  95. ERRE 21;
  96. FINSI;
  97.  
  98. * calcul de fnl ET dfnl/dx ? dans le doute on calcule tout...
  99. SI (NEG (TYPE IFFNL ) 'LOGIQUE'); IFFNL = VRAI; FINSI;
  100. SI (NEG (TYPE IFKNL) 'LOGIQUE'); IFKNL = VRAI; FINSI;
  101.  
  102.  
  103.  
  104. ******** Calcul de FNL ********
  105.  
  106. *---------------- Cas DEPNL LISTCHPO
  107. SI IFCHPO;
  108.  
  109. * preparation des donnees de sortie
  110. FNLTEMP = VIDE 'LISTCHPO';
  111. FNLm = prog;
  112. KNLp = prog;
  113. CNL = EXTR COMP_FNL 1;
  114.  
  115. * nbt : nombre de pas de temps
  116. nbt = DIME DEPNL;
  117. ib = 0;
  118. REPE BFOR nbt; ib = ib + 1;
  119.  
  120. DEP_ib = extr DEPNL ib;
  121. DEP_i = DEP_ib ;
  122. FNL_i = -1.*al* (DEP_i **3) ;
  123. FNLTEMP = FNLTEMP et FNL_i;
  124. FNLm = FNLm et (MAXI FNL_i 'ABS');
  125.  
  126. * pour le Calcul de DFNL :
  127. SI (IFKNL);
  128. xDEP_i = extr DEP_i pNL CNL;
  129. KNLp = KNLp et (3.*al* (xDEP_i **2));
  130. FINSI;
  131.  
  132. FIN BFOR;
  133. FNLmax = maxi (ABS FNLm);
  134.  
  135. *---------------- Cas DEPNL TABLE de LISTREEL
  136. SINON;
  137.  
  138. * DEPNL : table de listreel, indicée par les ddls NL
  139.  
  140. * nbt : nombre de pas de temps
  141. nbt = DIME DEPNL . 1;
  142. ib = 0;
  143. FNLTEMP = TABLE; FNLTEMP . 1 = PROG;
  144. KNLp = prog;
  145. REPE bfor nbt; ib = ib + 1;
  146. DEP_i0 = EXTR DEPNL . 1 ib;
  147. DEP_i = DEP_i0 ;
  148. FNLTEMP . 1 = FNLTEMP . 1 et (-1.*al* (DEP_i **3));
  149. SI (IFKNL);
  150. KNLp = KNLp et ( 3.*al* (DEP_i **2));
  151. FINSI;
  152. FIN bfor;
  153. FNLmax = maxi (ABS FNLTEMP . 1);
  154.  
  155. FINSI;
  156.  
  157.  
  158. ******** Calcul de DFNL ********
  159.  
  160. SI (IFKNL);
  161. * uniquement sous forme de LISTREELs
  162. KNLT1A = TABL;
  163. KNLT1A . 1 = tabl;
  164. KNLT1A . 1 . 1 = KNLp;
  165. FINSI;
  166.  
  167.  
  168. TAB1 . 'FNL_T' = FNLTEMP;
  169. TAB1 . 'KNL_T' = KNLT1A;
  170. TAB1 . 'FNLMAX' = FNLmax;
  171.  
  172. FINPROC;
  173. *-----------------------------------------------------------------------
  174.  
  175.  
  176.  
  177. ************************************************************************
  178. *
  179. * dynamique non lineaire
  180. * reponse forcee d'un oscillateur de type Duffing
  181. * reponse par HBM
  182. *
  183. ************************************************************************
  184.  
  185. *opti echo 0 ;
  186. opti debug 1;
  187. graph = faux ;
  188. OPTI DIME 2 ELEM SEG2 MODE PLAN ;
  189. OPTI EPSI LINEAIRE;
  190.  
  191. GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI;
  192. * GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX;
  193.  
  194.  
  195. *-------- Definition des parametres du calcul :
  196.  
  197. *nombre d'harmoniques utilises
  198. * nhbm = 1;
  199. nhbm = 3;
  200. * nhbm = 7;
  201. * nhbm = 11;
  202.  
  203. *parametre de l'adimensionnement
  204. u1max = 1.;
  205.  
  206. *des iterations
  207. itmoy = 6;
  208. *ds = 0.1;
  209. ds = 1.;
  210. dsmax = 1.;
  211. dsmin = 1.E-3 ;
  212.  
  213. * si AFT, OTFR tel que le nb de points pour TFR = 2**OTFR
  214. * OTFR = 7;
  215. OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2));
  216. mess 'OTFR = ' OTFR;
  217.  
  218.  
  219. *-------- Definition des donnees d'oscillateur ---------------
  220.  
  221. * M d2u/dt2 + C du/dt + K u + pfnl*u^3 = F(t)
  222. *mass
  223. pmass = 1.;
  224. *amortissement
  225. pamor = 0.1 ;
  226. *rigidite
  227. prigi = 1.;
  228. *force d'exterieur
  229. pfext = 0.5;
  230. *coefficient de la force non lineaire
  231. pfnl = 0.25;
  232. *pfnl = 1.;
  233. * pfnl = 5.;
  234. fprec = 1.E-3;
  235. si (ega pfnl 0.25 fprec) ; mofnl = mot '025'; finsi;
  236. si (ega pfnl 1.00 fprec) ; mofnl = mot '1'; finsi;
  237. si (ega pfnl 5.00 fprec) ; mofnl = mot '5'; finsi;
  238.  
  239. * prefix :
  240. prefix = chai 'hbm_duffing_mu_' mofnl '_n' nhbm ;
  241. si GRAPH ;
  242. ficps = chai prefix '.ps';
  243. opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps;
  244. finsi;
  245.  
  246. *plage d'etude de mu
  247. tcalc = prog 0.0 pas 0.1 5;
  248. evmu = evol manu 't' tcalc '\l' tcalc;
  249.  
  250.  
  251. *------------- geometrie ---------------
  252. *
  253. P1 = 0. 0. ;
  254. P2 = 1. 0. ;
  255. ST = P1 D 1 P2 ;
  256.  
  257. *----------- construit des matrices caracteristiques ----------
  258. *
  259. LIM1 = MOTS UY ;
  260. MASS1 = MASS 'UY' pmass P2 ;
  261. RI1 = MANU 'RIGI' P2 LIM1 (prog prigi);
  262. AMOR1 = MANU 'RIGI' P2 LIM1 (prog pamor);
  263.  
  264.  
  265.  
  266. *------- on définit une unique table pour tout (HBM et CONTINU) -------
  267. TAB1 = TABLE ;
  268.  
  269. *--------- remplissage des matrices "temporelles" ---------
  270.  
  271. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN rad/s ==> pas de 2pi
  272. TAB1 . 'MASSE_CONSTANTE' = MASS1 ;
  273. TAB1 . 'AMORTISSEMENT_CONSTANT' = AMOR1 ;
  274. TAB1 . 'RIGIDITE_CONSTANTE' = RI1 ;
  275. TAB1 . 'N_HARMONIQUE' = nhbm;
  276.  
  277. *------- remplissage des resultats attendus sur ddls "temporels" -------
  278. mycoul12 = mots VIOL BLEU BLEU TURQ TURQ OCEA OCEA VERT VERT OLIV OLIV
  279. JAUN JAUN ORAN ORAN ROUG ROUG ROSE ROSE AZUR AZUR GRIS GRIS DEFA DEFA ;
  280.  
  281. TAB1 . 'RESULTATS' = tabl;
  282. TAB1 . 'RESULTATS' . 1 = tabl;
  283. TAB1 . 'RESULTATS' . 1 . 'POINT_MESURE' = P2;
  284. TAB1 . 'RESULTATS' . 1 . 'COMPOSANTE' = mot 'UY';
  285. TAB1 . 'RESULTATS' . 1 . 'TITRE' = mot 'U';
  286. TAB1 . 'RESULTATS' . 1 . 'COULEUR' = mot 'BLEU';
  287. TAB1 . 'RESULTATS' . 1 . 'COULEUR_HBM' = mycoul12;
  288.  
  289.  
  290. *--------- passage en frequentiel ---------
  291.  
  292. HBM TAB1;
  293.  
  294.  
  295. MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP TEMPORELLES';
  296. LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT';
  297. MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP FREQUENTIELLES';
  298. LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT_HBM';
  299.  
  300. MESS 'RIGIDITE_HBM'; LIST TAB1 . 'RIGIDITE_HBM' ;
  301. MESS 'BLOCAGES_HBM'; LIST TAB1 . 'BLOCAGES_HBM' ;
  302. MESS 'AMORTISSEMENT_HBM'; LIST TAB1 . 'AMORTISSEMENT_HBM' ;
  303. MESS 'MASSE_HBM'; LIST TAB1 . 'MASSE_HBM' ;
  304. MESS 'MASSE_HBM_1'; LIST TAB1 . 'MASSE_HBM_1' ;
  305. MESS 'MASSE_HBM_2'; LIST TAB1 . 'MASSE_HBM_2' ;
  306. MESS 'AMORTISSEMENT_HBM_1'; LIST TAB1 . 'AMORTISSEMENT_HBM_1';
  307.  
  308.  
  309. *----------- definition du chargement ----------------------------------
  310.  
  311. * ici, on definit directement le chargement en frequentiel
  312. * F2 : composante en cos 1wt
  313. FP11 = MANU 'CHPO' p2 1 'F2' pfext 'NATURE' 'DISCRET' ;
  314. *
  315. LIX1 = PROG 0. pas 0.1 (10.) ;
  316. LIY1 = PROG (dime LIX1)*1.;
  317. EV1 = EVOL MANU 'w' LIX1 'F(w)' LIY1 ;
  318. CHA1 = CHAR MECA FP11 EV1 ;
  319.  
  320. * frequence associee (en rad/s) = constante = 1.6*w0
  321. Wev = evol manu 't' LIX1 'w' (1.6*LIY1);
  322. TAB1 . 'FREQUENCE' = Wev;
  323.  
  324.  
  325. *----------- OPTIONS DE CALCUL ----------------------
  326.  
  327. TAB1 . 'GRANDS_DEPLACEMENTS' = FAUX ;
  328. TAB1 . 'GEOMETRIE' = P2;
  329. TAB1 . 'CHARGEMENT' = CHA1;
  330. TAB1 . 'PRECISION' = 1.E-5;
  331. TAB1 . 'TEMPS_CALCULES' = tcalc;
  332. TAB1 . 'MAXIPAS' = 500;
  333. TAB1 . 'ACCELERATION' = 4;
  334. TAB1 . 'NB_ITERATION' = itmoy;
  335. TAB1 . 'MAXITERATION' = 30;
  336. TAB1 . 'WTABLE' = tabl;
  337. TAB1 . 'WTABLE' . 'DS' = ds;
  338. TAB1 . 'WTABLE' . 'DSMAX' = dsmax;
  339. TAB1 . 'WTABLE' . 'DSMIN' = dsmin;
  340. TAB1 . 'MAXI_DEPLACEMENT' = u1max;
  341. TAB1 . 'PAS_SAUVES' = 2;
  342.  
  343. * calcul des forces NL
  344. TAB1 . 'PROCEDURE_CHAR_MECA'= vrai ;
  345. TAB1 . 'CALC_CHPO' = faux;
  346. TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT';
  347. TAB1 . 'N_PT_TFR' = OTFR;
  348.  
  349. * parametres du cas test lié a AFT et CHARMECA
  350. TAB1 . 'COEFF_FNL' = evmu;
  351. TAB1 . 'POINT_FNL' = P2;
  352. TAB1 . 'COMP_FNL' = mots 'UY' ;
  353.  
  354. * stabilité
  355. TAB1 . 'STABILITE' = MOTS 'DIAG' 'HILL';
  356.  
  357.  
  358. *----------- resolution par la procedure CONTINU -----------------------
  359.  
  360. *calcul avec continuation
  361. CONTINU TAB1;
  362.  
  363. TEMP 'IMPR' 'MAXI';
  364.  
  365. *------------------------- post-traitement -----------------------------
  366.  
  367.  
  368. ********** tracé de toutes les harmoniques individuellement **********
  369.  
  370. wprog = TAB1 . 'TEMPS_PROG';
  371. evtot = TAB1 . 'RESULTATS_HBM' . 'RESULTATS_EVOL';
  372.  
  373. TDESS1 = TABL;
  374. TDESS1 . 'TITRE' = TAB1 . 'RESULTATS_HBM' . 'TITRE';
  375.  
  376. si GRAPH;
  377. dess evtot
  378. TITX '\m' POSX CENT
  379. TITY 'U_{k}' POSY CENT TDESS1 LEGE NO;
  380. finsi;
  381.  
  382. ******************* uy2p : amplitude en norme 2 ******************
  383.  
  384. ihbm = 0;
  385. uy2p = TAB1 . 'RESULTATS_HBM' . 1 . 'RESULTATS' ** 2;
  386. repe BUYHBM nhbm; ihbm = ihbm + 1;
  387. ikcos = 2*ihbm;
  388. iksin = 2*ihbm + 1;
  389. uy2p = uy2p
  390. + (0.5* ( (TAB1 . 'RESULTATS_HBM' . ikcos . 'RESULTATS' ** 2)
  391. + (TAB1 . 'RESULTATS_HBM' . iksin . 'RESULTATS' ** 2)) );
  392. fin BUYHBM;
  393.  
  394. *** evolution + tracé ***
  395. evuy2 = evol VIOL MANU wprog uy2p;
  396. si GRAPH;
  397. dess evuy2
  398. TITX '\m' POSX CENT
  399. TITY '|U|_{2}' POSY CENT ;
  400. finsi;
  401.  
  402. ********** traduction des resultats frequentiels en temporels **********
  403.  
  404. HBM_POST TAB1;
  405.  
  406. ************* tracé du max d'amplitude max_t(u(t)) *********************
  407.  
  408. evuyamp = TAB1 . 'RESULTATS' . 'RESULTATS_EVOL';
  409. si GRAPH;
  410. dess evuyamp
  411. TITX '\m' POSX CENT
  412. TITY 'max|U(t)|' POSY CENT ;
  413. finsi;
  414.  
  415. ************* tracé du max d'amplitude avec la stabilite ! *************
  416.  
  417. istab = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'STABILITE';
  418. Tdess2 = TABL;
  419. Tdess2 . 'LIGNE_VARIABLE' = TABL;
  420. Tdess2 . 'LIGNE_VARIABLE' . 1 = istab;
  421. si GRAPH;
  422. dess evuyamp
  423. TITX '\m' POSX CENT
  424. TITY 'max|U(t)|' POSY CENT Tdess2;
  425. * zoom
  426. Tdess2 . 1 = mot 'MARQ S PLUS';
  427. dess evuyamp
  428. TITX '\m' POSX CENT XBOR 1.55 1.70 YBOR 2.8 3.1
  429. TITY 'max|U(t)|' POSY CENT Tdess2;
  430. dess evuyamp
  431. TITX '\m' POSX CENT XBOR 1.28 1.31 YBOR 0.8 1.5
  432. TITY 'max|U(t)|' POSY CENT Tdess2;
  433. finsi;
  434.  
  435. * autres tracé de la stabilite : lambda_R en fonction de W
  436. mycoul = mots 'VIOL' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'OLIV'
  437. 'JAUN' 'ORAN' 'ROUG' 'ROSE';
  438. ncoul = dime mycoul;
  439. wprog2 = enle wprog (dime wprog);
  440. evlrtot = VIDE 'EVOLUTIO' ;
  441. evlitot = VIDE 'EVOLUTIO' ;
  442. nl = dime TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL';
  443. il = 0;
  444. icoul = 0;
  445. repe bl nl; il = il + 1; icoul = icoul + 1;
  446. si(icoul > ncoul); icoul = 1; finsi; coul1 = extr mycoul icoul;
  447. lr1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL' . il;
  448. evlrtot = evlrtot
  449. et (evol coul1 MANU '\m' wprog2 '\l_{R}' lr1);
  450. li1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_IMAG' . il;
  451. evlitot = evlitot
  452. et (evol coul1 MANU '\m' wprog2 '\l_{I}' li1);
  453. fin bl;
  454. wprog3 = prog (mini wprog2) (maxi wprog2);
  455. evzero = evol manu '\m' wprog3 '\l_{R}' (0.*wprog3);
  456. Tdess3 = tabl;
  457. Tdess3 . 1 = mot 'TIRR';
  458. si GRAPH;
  459. dess (evzero et evlrtot)
  460. 'TITY' '\l_{R} (s^{-1})' Tdess3;
  461. dess (evlitot) 'TITY' '\l_{I} (Hz)' ;
  462. finsi;
  463.  
  464.  
  465. *-------------------------- sauvegarde -----------------------------
  466.  
  467. si SAVE ;
  468. ficxdr = chai prefix '.xdr';
  469. mess 'sauvegarde dans ' ficxdr;
  470. opti sauv ficxdr;
  471. finsi;
  472.  
  473.  
  474. *--------------------- test de non regression ---------------------
  475.  
  476. si TNR ;
  477.  
  478. * recup valeurs de la courbes de reponse
  479. wamp = extr evuyamp 'ABSC';
  480. uyamp = extr evuyamp 'ORDO';
  481. namp = DIME uyamp;
  482.  
  483. * on teste la presence de 2 points limites
  484. dw1 = ( wamp enle 1) - ( wamp enle namp);
  485. dw1dw0 = (dw1 enle 1) * (dw1 enle (namp - 1));
  486. mptlim = (lect 0) et (enti 'PROCH' (masq dw1dw0 'EGINFE' 0.)) et 0;
  487. nptlim = somm mptlim;
  488. si (nptlim neg 2);
  489. mess 'Il devrait y avoir 2 points limites !';
  490. list nptlim;
  491. ERRE 5;
  492. finsi;
  493.  
  494. * on verifie que pour ces 2 points il y a chgt de stabilité
  495. iptlim = posi 1 'DANS' mptlim 'TOUS' ;
  496. istab1 = (istab enle 1) - (istab enle (DIME istab));
  497. istab1 = (posi 1 'DANS' istab1 'TOUS')
  498. et (posi -1 'DANS' istab1 'TOUS');
  499. * sumstab = (extr istab (iptlim - 1)) + (extr istab (iptlim));
  500. * * sumstab doit valoir 1 (ni 0 ni 2)
  501. * si ( ((mini sumstab) neg 1) ou ((maxi sumstab) neg 1) );
  502. * ==> on ne fait pas le test car precis a 1 pas pres !!!
  503. dstab = maxi (iptlim - istab1) 'ABS';
  504. si (dstab >eg 2);
  505. mess 'Les 2 points limites doivent correspondre a un changement '
  506. 'de stabilite!';
  507. list sumstab;
  508. ERRE 5;
  509. finsi;
  510.  
  511. * on termine par verifier l'amplitude en ces 2 points + l'amplitude max
  512. Utest = (extr uyamp iptlim) et (maxi uyamp) ;
  513. * on compare a une precedente execution (lors de la creation du cas test)
  514. Uref = prog 0.48057 3.1837 3.1847;
  515. si ( (MAXI (Utest - Uref) 'ABS') > 0.01 );
  516. mess 'Erreur dans les amplitudes predites !';
  517. list Utest;
  518. ERRE 5;
  519. finsi;
  520.  
  521. finsi;
  522.  
  523.  
  524. fin ;
  525.  
  526.  
  527.  
  528.  

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