Télécharger dync01.dgibi

Retour à la liste

Numérotation des lignes :

  1. *********************************************************************
  2. * *
  3. * RÉPONSE FORCEE D'UN OSCILLATEUR DE DUFFING *
  4. * OPTION EXPOSANT DE L'OPERATEUR DYNE (raideur cubique) *
  5. * *
  6. * CREATION : 25/01/2018 *
  7. * *
  8. * MANIP KOALA (CEA,DEN,DM2S,SEMT,DYN) *
  9. * *
  10. * *
  11. * _____ *
  12. * | | *
  13. * | m | *
  14. * |_____| *
  15. * | *
  16. * /| / | *
  17. * /|----[k]----[xi]----(knl)---- *
  18. * /| / ----> F(t) *
  19. * |----> x *
  20. * *
  21. *********************************************************************
  22.  
  23.  
  24. *********************************************************************
  25. * Parametres globaux
  26. *********************************************************************
  27. GRAPH = FAUX;
  28.  
  29. OPTION DIME 2 ELEM SEG2 ;
  30. OPTION 'TRAC' 'PSC' 'EPTR' 5 'POTR' 'TIMES_16';
  31. OPTION 'EPSI' 'LINEAIRE';
  32.  
  33.  
  34. *********************************************************************
  35. * PARAMETRES DU MODELE
  36. *********************************************************************
  37. * nombre de modes a calculer
  38. NMODE = 1;
  39. * nombre de harmoniques a calculer et NFFT
  40. NHBM = 5; NFFT = 256;
  41. * Amplitude chargement
  42. Fext0 = 1.;
  43. * Parametre non-lineaire fnl=knl*x^3
  44. knl = 0.25 ;
  45.  
  46.  
  47. *********************************************************************
  48. * DEFINITION DU MODELE
  49. *********************************************************************
  50.  
  51. *--------------------- Characteristiques du modele ----------------
  52. *masse frequence propre (Hz), raideur
  53. xi = 0.1/2.;
  54. m = 1.;
  55. k = 1.;
  56. *point de la masse
  57. p1 = poin 0. 0. ;
  58.  
  59. *------------------Calcul des Modes Propres --------------
  60.  
  61. p1A = point 0. 0.;
  62. phi1 = manu chpo p1 1 'UX' 1. NATURE DIFFUS;
  63.  
  64. TMOD = TABLE 'BASE_MODALE';
  65. TMOD . 'MODES' = TABL 'BASE_DE_MODES';
  66. TMOD . 'MODES' . 'MAILLAGE' = manu 'POI1' p1;
  67. TMOD . 'MODES' . 1 = TABL 'MODE';
  68. TMOD . 'MODES' . 1 . 'POINT_REPERE' = p1A ;
  69. TMOD . 'MODES' . 1 . 'NUMERO_MODE' = 1;
  70. TMOD . 'MODES' . 1 . 'FREQUENCE' = 1./(2.*PI);
  71. TMOD . 'MODES' . 1 . 'MASSE_GENERALISEE' = m;
  72. TMOD . 'MODES' . 1 . 'DEFORMEE_MODALE' = phi1;
  73.  
  74.  
  75. *********************************************************************
  76. * LIAISONS, CHARGEMENT ET AMORTISSEMENT
  77. *********************************************************************
  78.  
  79. *----------------------------------------------*
  80. * CHARGEMENT = EXCITATION HARMONIQUE
  81. *----------------------------------------------*
  82. GFC1 = MANU CHPO p1A 1 'FALF' Fext0 ;
  83.  
  84. *----------------------------------------------*
  85. * Les FORCES due a la raideur non-lineaire
  86. *----------------------------------------------*
  87.  
  88. TL3 = TABLE 'LIAISON_ELEMENTAIRE' ;
  89. TL3 . 'TYPE_LIAISON' = MOT 'COUPLAGE_DEPLACEMENT';
  90. TL3 . 'SUPPORT' = p1A ;
  91. TL3 . 'ORIGINE' = p1A ;
  92. TL3 . 'COEFFICIENT' = -1.*knl ;
  93. TL3 . 'EXPOSANT' = 3 ;
  94.  
  95. *----------------------------------------------*
  96. * Regroupement des LIAISONS
  97. *----------------------------------------------*
  98. TAB_LIAI = TABLE 'LIAISON' ;
  99. TAB_LIAI . 'LIAISON_A' = TABLE 'LIAISON_A';
  100. TAB_LIAI . 'LIAISON_A' . 1 = TL3;
  101.  
  102. *----------------------------------------------*
  103. * L'AMORTISSEMENT (attention en % pour AMOR)
  104. *----------------------------------------------*
  105. TAB_AMOR = TABLE 'AMORTISSEMENT';
  106. TAB_AMOR . 'AMORTISSEMENT' = AMOR TMOD (100.*(prog xi));
  107.  
  108. *********************************************************************
  109. * CALCUL NON LINEAIRE REPONSE FORCEE : TABLES DYNC
  110. *********************************************************************
  111. *nouvelle syntaxe pour TAB_CHAR
  112. * chargement spatial : GFC1
  113. * sur l'harmonique fondamentale 1 ==> F=Fext*cos(wt)
  114. TAB_CHAR = TABLE 'CHARGEMENT';
  115. TAB_CHAR . 1 = GFC1;
  116. * constant avec la frequence, si F=w^2*(...) ==> BALOURD = 1 (eg. JEFFCOTT)
  117. * TAB_CHAR . 'BALOURD' = 0;
  118.  
  119. * Approximation pour la solution initiale à fréquence fixe
  120. VEC_INIT = TABLE 'INITIAL';
  121. Wext = 0.01;
  122. VEC_INIT . 'FREQUENCE' = Wext ;
  123. *bp VEC_INIT . 2 = 0.5;
  124. *bp : on utilise desormais la convention j=0 +1 -1 +2 -1 ... +nhbm -nhbm
  125. *bp ainsi qu'un chpoint
  126. Q0 = MANU CHPO p1A 1 'ALFA' 0.5 ;
  127. VEC_INIT . +1 = Q0;
  128.  
  129. * Paramètres numériques pour la continuation
  130. PAR_NUM = TABLE 'PARAMETRES_NUMERIQUES';
  131. PAR_NUM . 'TYPE' = MOT 'FORC';
  132. PAR_NUM . 'VAL_FIN' = 2.5;
  133. PAR_NUM . 'DS0' = 0.01;
  134. PAR_NUM . 'DSMAX' = 0.05;
  135. PAR_NUM . 'DSMIN' = 1.E-4;
  136. PAR_NUM . 'ITERMOY' = 3.4;
  137. PAR_NUM . 'ITERMAX' = 30;
  138. PAR_NUM . 'ANGLE_MIN' = 0.;
  139. PAR_NUM . 'ANGLE_MAX' = 50.;
  140. PAR_NUM . 'ISENS' = 1.;
  141. PAR_NUM . 'NBPAS' = 1000;
  142. *----------------------------------------------------------------------*
  143. * Appel a l'operateur
  144. * opti impi 3;
  145. TCON = DYNC TMOD TAB_CHAR TAB_LIAI TAB_AMOR VEC_INIT PAR_NUM
  146. NHBM NFFT;
  147. * opti impi 0 donn 5;
  148. *----------------------------------------------------------------------*
  149.  
  150.  
  151. *********************************************************************
  152. * CALCUL NON LINEAIRE MODE NON-LINEAIRE : TABLES DYNC
  153. *********************************************************************
  154.  
  155. * qq modifs dans les donnees d'entree
  156. PAR_NUM . 'TYPE' = MOT 'NNM';
  157. * PAR_NUM . 'VAL_INI' = 0.001; ---> VEC_INIT . 'ENERGIE'
  158. PAR_NUM . 'VAL_FIN' = 25;
  159. * PAR_NUM . 'MODE_DEBUT' = 1; ---> VEC_INIT . 'MODE' = 1
  160. *
  161. *-----------------------------------------------------------------------
  162. * Appel a l'operateur
  163. * rem. : a) pour l'instant, pas de TAB_AMOR
  164. * b) VEC_INIT basée sur l'energie initiale
  165. VEC_INIT = TABLE 'INITIAL';
  166. VEC_INIT . 'MODE' = 1;
  167. VEC_INIT . 'ENERGIE' = 0.001;
  168. opti impi 3;
  169. TNNM = DYNC TMOD TAB_CHAR TAB_LIAI PAR_NUM VEC_INIT NHBM NFFT;
  170. opti impi 0;
  171. *-----------------------------------------------------------------------
  172.  
  173.  
  174. ************************************************************************
  175. * POST-TRAITEMENT
  176. ************************************************************************
  177.  
  178. * BACKBONE
  179. * --------
  180. * frequence
  181. FNNM = TNNM . 'REPONSE' . 'FREQUENCE';
  182. * MAX de la NORME2 des COEFFICIENTS de FOURIER pour chaque freq
  183. ANNM = TNNM . 'REPONSE' . 'NORME_DEPLACEMENT';
  184. EVNNM = EVOL MANU '\w' FNNM '|Q|' ANNM;
  185.  
  186. * COURBE DE REPONSE FORCE
  187. * -----------------------
  188. * frequence
  189. FREQ = TCON . 'REPONSE' . 'FREQUENCE';
  190. * MAX de la NORME2 des COEFFICIENTS de FOURIER pour chaque freq
  191. AMPS = TCON . 'REPONSE' . 'NORME_DEPLACEMENT';
  192. XY = EVOL MANU '\w' FREQ '|Q|' AMPS;
  193. * stabilite
  194. STAB = TCON . 'REPONSE' . 'STABILITE';
  195. TABSTB = TABLE;
  196. TABSTB . 'LIGNE_VARIABLE' = TABL ENTI;
  197. TABSTB . 'LIGNE_VARIABLE' . 1 = STAB;
  198. TITR 'Duffing : courbe de reponse';
  199. DESS XY TABSTB ;
  200.  
  201. TITR 'Duffing : courbe de reponse + NNM';
  202. DESS (XY ET EVNNM) TABSTB XBOR 0.0 2.5 XGRA 0.5 YBOR 0.0 5.0 YGRA 1.0; ;
  203.  
  204.  
  205. * BIFURCATIONS
  206. * ------------
  207.  
  208. * recup
  209. NBIF = DIME TCON . 'BIFURCATION' . 'TYPE';
  210. prLP_Q = PROG; prLP_w = PROG;
  211. prBP_Q = PROG; prBP_w = PROG;
  212. prPD_Q = PROG; prPD_w = PROG;
  213. prNS_Q = PROG; prNS_w = PROG;
  214. REPE BBIF NBIF;
  215. typbif = EXTR TCON . 'BIFURCATION' . 'TYPE' &BBIF;
  216. Qbif = EXTR TCON . 'BIFURCATION' . 'NORME_DEPLACEMENT' &BBIF;
  217. wbif = EXTR TCON . 'BIFURCATION' . 'FREQUENCE' &BBIF;
  218. SI (EGA typbif 'LP'); prLP_Q = prLP_Q et Qbif; prLP_w = prLP_w et wbif; FINSI;
  219. SI (EGA typbif 'BP'); prBP_Q = prBP_Q et Qbif; prBP_w = prBP_w et wbif; FINSI;
  220. SI (EGA typbif 'PD'); prPD_Q = prPD_Q et Qbif; prPD_w = prPD_w et wbif; FINSI;
  221. SI (EGA typbif 'NS'); prNS_Q = prNS_Q et Qbif; prNS_w = prNS_w et wbif; FINSI;
  222. FIN BBIF;
  223.  
  224. * creation d'1 evolution + table de dessin
  225. ibif = 1;
  226. evBIF = VIDE 'EVOLUTIO';
  227. * Limit Point
  228. si ((dime prLP_w) > 0); ibif = ibif + 1;
  229. evLP = EVOL MANU LEGE 'LP' '\w' prLP_w '|Q|' prLP_Q;
  230. evBIF = evBIF et evLP;
  231. TABSTB . ibif = MOT 'NOLI MARQ LOSA';
  232. finsi;
  233. * Branch Point
  234. si ((dime prBP_w) > 0); ibif = ibif + 1;
  235. evBP = EVOL MANU LEGE 'BP' '\w' prBP_w '|Q|' prBP_Q;
  236. evBIF = evBIF et evBP;
  237. TABSTB . ibif = MOT 'NOLI MARQ CARR';
  238. finsi;
  239. * Period Doubling Point
  240. si ((dime prPD_w) > 0); ibif = ibif + 1;
  241. evPD = EVOL MANU LEGE 'PD' '\w' prPD_w '|Q|' prPD_Q;
  242. evBIF = evBIF et evPD;
  243. TABSTB . ibif = MOT 'NOLI MARQ TRID';
  244. finsi;
  245. * Neimark Sacker Point
  246. si ((dime prNS_w) > 0); ibif = ibif + 1;
  247. evNS = EVOL MANU LEGE 'NS' '\w' prNS_w '|Q|' prNS_Q;
  248. evBIF = evBIF et evNS;
  249. TABSTB . ibif = MOT 'NOLI MARQ TRIU';
  250. finsi;
  251.  
  252. TITRE 'Duffing : courbe de reponse + bifurcations';
  253. DESS (XY ET evBIF) TABSTB ;
  254. DESS (XY ET evBIF) TABSTB XBOR 0.1 0.5 XGRA 0.1;
  255. DESS (XY ET evBIF) TABSTB XBOR 2.10 2.25 XGRA 0.05 YBOR 4.2 4.5 YGRA 0.05;
  256.  
  257. * STABILITE via les EXPOSANTS DE FLOQUET
  258. * --------------------------------------
  259. TEXP_R = TCON . 'REPONSE' . 'EXPOSANT_REEL';
  260. TEXP_I = TCON . 'REPONSE' . 'EXPOSANT_IMAGINAIRE';
  261. exp_R = VIDE 'EVOLUTIO';
  262. exp_I = VIDE 'EVOLUTIO';
  263. coco = MOTS 'ROUG' 'BLEU' 'VERT' 'AZUR' 'ORAN' 'VIOL' 'TURQ' ;
  264. repe bexp;
  265. si (non (exis TEXP_R &bexp)); quit bexp ; finsi;
  266. coco_i = EXTR coco &bexp;
  267. exp_R = exp_R ET (evol coco_i 'MANU' '\w' FREQ '\m_{R}' TEXP_R . &bexp);
  268. exp_I = exp_I ET (evol coco_i 'MANU' '\w' FREQ '\m_{I}' TEXP_I . &bexp);
  269. fin bexp;
  270.  
  271. TITR 'Duffing : exposants de Floquet';
  272. DESS exp_R;
  273. DESS exp_I;
  274.  
  275.  
  276. * DETAIL DES HARMONIQUES DU 1ER MODE
  277. * ----------------------------------
  278. * imode = 1;
  279. * le POINT_REPERE est p1A
  280. mycoco = @PALETTE (nhbm + 1);
  281. TABSTB2 = TABLE;
  282. TABSTB2 . 'LIGNE_VARIABLE' = TABL ENTI;
  283. ev1 = VIDE 'EVOLUTIO';
  284. j = 0;
  285. repe bhbm (nhbm + 1);
  286. si (EGA j 0);
  287. AMPj = TCON . REPONSE . COEFFICIENTS . j . p1A;
  288. sinon;
  289. AMPjC = TCON . REPONSE . COEFFICIENTS . j . p1A;
  290. AMPjS = TCON . REPONSE . COEFFICIENTS . (-1*j) . p1A;
  291. AMPj = ((AMPjC**2) + (AMPjS**2))**0.5;
  292. finsi;
  293. coco = EXTR mycoco &bhbm;
  294. ev1 = ev1 ET (EVOL coco 'MANU' 'LEGE' (chai 'j='j) '\w' FREQ '|Q^{j}|' AMPj);
  295. TABSTB2 . 'LIGNE_VARIABLE' . &bhbm = STAB;
  296. j = j + 1;
  297. fin bhbm;
  298.  
  299. TITR 'Duffing : Reponse du 1er mode';
  300. DESS ev1 'LEGE' TABSTB2;
  301. DESS ev1 'LEGE' XBOR 0.1 0.5 XGRA 0.1 TABSTB2;
  302.  
  303.  
  304. * TEMPOREL ET PORTRAITS DE PHASE POUR QUELQUES POINTS
  305. * ---------------------------------------------------
  306.  
  307. * base de fonctions T (et leur derivees T2) sur une periode adimensionnee
  308. T = TABL;
  309. T2 = TABL;
  310. time = PROG 0. PAS (360./NFFT) 360.;
  311. * terme constant
  312. T . 1 = PROG (NFFT + 1)*1.;
  313. T2 . 1 = PROG (NFFT + 1)*0.;
  314. * harmonique +j et -j
  315. j = 0;
  316. REPE bhbm nhbm; j = j + 2;
  317. T . j = COS (&bhbm * time) ;
  318. T . (j + 1) = SIN (&bhbm * time) ;
  319. T2 . j = -1. * &bhbm * T . (j + 1);
  320. T2 . (j + 1) = &bhbm * T . j ;
  321. FIN bhbm;
  322.  
  323. * boucle sur les solutions a tracer
  324. wprog = PROG 0.23 0.39 1.0 2.0 ;
  325. REPE BW (DIME wprog);
  326. w = EXTR wprog &BW;
  327. ipas = POSI 1. 'DANS' (MASQ FREQ 'EGSUPE' w);
  328. w = EXTR FREQ ipas;
  329. * harmonique 0
  330. Qprog = PROG (EXTR TCON . REPONSE . COEFFICIENTS . 0 . p1A ipas);
  331. * harmonique +j et -j
  332. REPE bhbm nhbm; j= &bhbm; jm = -1 * j;
  333. Qprog = Qprog et (EXTR TCON . REPONSE . COEFFICIENTS . j . p1A ipas)
  334. et (EXTR TCON . REPONSE . COEFFICIENTS . jm . p1A ipas);
  335. FIN bhbm;
  336. * q(t) = [T] * Q
  337. qtime = COLI T Qprog;
  338. dotqtime = COLI T2 Qprog;
  339. * trace du portrait de phase
  340. evphase = EVOL 'AZUR' 'MANU' 'q' qtime 'dq/dt' dotqtime;
  341. DESS evphase 'TITR' (CHAI 'Duffing : w=' FORMAT '(F6.3)' w) 'CARR';
  342. FIN BW;
  343.  
  344.  
  345. ************************************************************************
  346. * TEST DE BON FONCTIONNEMENT
  347. ************************************************************************
  348.  
  349. * listreel des erreurs
  350. prERR = PROG;
  351.  
  352. * amplitude max et w correspondant a la 1ere resonance (harmonique 1)
  353. wref1 = 2.210;
  354. Qref1 = 4.475;
  355. imax1 wmax1 Qmax1 = MAXI XY;
  356. prERR = prERR et (wmax1 - wref1) et (Qmax1 - Qref1);
  357.  
  358. * resonance harmonique 3
  359. wref3 = 0.38;
  360. Qref3 = 1.03;
  361. XY3 = EXTR XY 'COMP' 0.3 1.5;
  362. FREQ3 = EXTR XY3 'ABSC';
  363. AMPS3 = EXTR XY3 'ORDO';
  364. dXY3 = (ENLE AMPS3 1) - (ENLE AMPS3 (dime AMPS3));
  365. imax3= POSI 1. 'DANS' (dXY3 MASQ EGINFE 0.);
  366. wmax3 = EXTR FREQ3 imax3;
  367. Qmax3 = EXTR AMPS3 imax3;
  368. prERR = prERR et (wmax3 - wref3) et (Qmax3 - Qref3);
  369.  
  370. * test
  371. ERRMAX = MAXI prERR 'ABS';
  372. SI (ERRMAX > 0.05);
  373. LIST prERR;
  374. ERRE 5;
  375. FINSI;
  376.  
  377. * nombre et/ou type et/ou localisation des bifurcations
  378. SI ( (NBIF NEG 2) OU ((dime prLP_w) NEG 2) );
  379. ERRE 5;
  380. FINSI;
  381.  
  382.  
  383. ************************************************************************
  384. * PERFORMANCES
  385. ************************************************************************
  386. TEMP IMPR MAXI CPU;
  387.  
  388. * opti donn 5 trac X;
  389. FIN ;
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  

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