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

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