Télécharger la plaquette Cast3M


fr en

Télécharger dync01.dgibi

Retour à la liste

Numérotation des lignes : 

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

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