Télécharger la plaquette Cast3M


fr en

Télécharger dync02.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : dync02.dgibi
  2. *
  3. ************************************************************************
  4. *
  5. * Calcul d'un rotor de type Jeffcott avec contact frottant
  6. * avec la methode HBM (DYNC)
  7. *
  8. *                 ___                           ^ Z
  9. *                 _|_ kc,mu                     |__
  10. *               _______  jeu                    |_ \_ 
  11. *   |   k,c    |       |                        | \  \
  12. *   |----------|   m   |----------|             +--|--|-----> Y
  13. *   |          |_______|                           R  R+jeu
  14. *                 ___    jeu                       
  15. *                 _|_ kc,mu
  16. *                
  17. *  ^ Z
  18. *  |
  19. *  +----> X
  20. *  
  21. *
  22. * Ref : [Xie et al., MSSP, 2015]
  23. * Auteur : PBZ CEA-ENSTA, 2020-07-07
  24. *
  25. ************************************************************************
  26.  
  27.  
  28. *********************************************************************
  29. * OPTIONS GENERALES
  30. *********************************************************************
  31.  
  32.   OPTI DIME 3 ELEM SEG2 ;
  33.   OPTION 'TRAC' 'PSC' 'EPTR' 8 'POTR' 'HELVETICA_16';
  34.  
  35.  
  36. *********************************************************************
  37. *     PARAMETRES 
  38. *********************************************************************
  39.  
  40. * nombre de modes a calculer
  41.   NMODE = 2;
  42. * nombre de harmoniques a calculer et NFFT
  43.   NHBM = 7;
  44.   NFFT = 256;
  45.  
  46. * mu de frottement  
  47. *  mu = 0.05;
  48.   mu = 0.2;
  49.  
  50.  
  51. *********************************************************************
  52. *     DEFINITION DU MODELE 
  53. *********************************************************************
  54.  
  55. * rotor data
  56.   m = 1.  ;
  57.   k = 100.;
  58.   c = 5.  ;
  59.   jeu  = 0.105;
  60.   R = 20.*jeu ;
  61.   ebal = 0.1;
  62. *   Fbal = ebal * m * (Omega**2) ; avec Omega = vitesse de rotation en rad/s
  63.  
  64. * contact parameters
  65.   kchoc= 25. * k;
  66.   mu_s = mu;
  67.   mu_d = mu;
  68.   si (mu_s < mu_d); erreur 21; finsi;
  69.   Kt = 10.*kchoc;
  70.   Ct = 0.5*(m*Kt)**0.5;
  71.  
  72. * deduced parameter
  73.   w   = (k/m)**0.5;                 
  74.   wHz = w / (2.*pi);
  75.   wchoc = (kchoc/m)**0.5;           
  76.   xi  = 0.5 * c / ((kchoc*m)**0.5); 
  77.   mess 'w=' w '(' wHz ')Hz < wchoc=' wchoc; 
  78.   mess 'xi=' xi;
  79.  
  80. * critical timestep (Diff Centrees, no damping)
  81.   wadhe = ((k+Kt)/m)**0.5;
  82.   dtc = 2. / wadhe;
  83.  
  84. * !!! on adimensionne apres !!! 
  85.  
  86.  
  87. ************************************************************************
  88. *                       BASE MODALE 
  89. ************************************************************************
  90.  
  91. * point physique
  92.   p1 = 0. 0. 0.;
  93.  
  94. * base modale = 2 modes
  95.   palfa1 = 0. 0. 0.;
  96.   phi1 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 1. 'UZ' 0. 'NATURE' 'DIFFUS';
  97.   palfa2 = 0. 0. 0.;
  98.   phi2 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 0. 'UZ' 1. 'NATURE' 'DIFFUS';
  99.  
  100.   TBAS = TABL 'BASE_MODALE';
  101.   TBAS . 'MODES' = TABL 'BASE_DE_MODES'; 
  102.   TBAS . 'MODES' . 1 = TABL 'MODES';
  103.   TBAS . 'MODES' . 1 . 'POINT_REPERE'      = palfa1;
  104.   TBAS . 'MODES' . 1 . 'FREQUENCE'         = wHz/wchoc;
  105.   TBAS . 'MODES' . 1 . 'MASSE_GENERALISEE' = 1.;
  106.   TBAS . 'MODES' . 1 . 'DEFORMEE_MODALE'   = phi1;
  107.   TBAS . 'MODES' . 2 = TABL 'MODES';
  108.   TBAS . 'MODES' . 2 . 'POINT_REPERE'      = palfa2;
  109.   TBAS . 'MODES' . 2 . 'FREQUENCE'         = wHz/wchoc;
  110.   TBAS . 'MODES' . 2 . 'MASSE_GENERALISEE' = 1.;
  111.   TBAS . 'MODES' . 2 . 'DEFORMEE_MODALE'   = phi2;
  112.  
  113.  
  114. ************************************************************************
  115. *                       AMORTISSEMENT 
  116. ************************************************************************
  117.  
  118.   TAMOR = TABL 'AMORTISSEMENT';
  119.   Ctot = MANU 'RIGIDITE' 'TYPE' 'AMORTISSEMENT' 
  120.         (palfa1 et palfa2) (MOTS 'ALFA') (prog (c / (m*wchoc)));
  121.   TAMOR . 'AMORTISSEMENT' = Ctot  ;
  122.  
  123.  
  124. ************************************************************************
  125. *                           LIAISONS 
  126. ************************************************************************
  127.  
  128. * contact-frottant point_cercle_frottement
  129.   TL1 = TABL 'LIAISON_ELEMENTAIRE';
  130.   TL1 . 'TYPE_LIAISON' = mot 'POINT_CERCLE_FROTTEMENT';
  131.   TL1 . 'SUPPORT'      = p1;
  132.   TL1 . 'NORMALE'      = (1. 0. 0.) ;
  133.   TL1 . 'EXCENTRATION' = (0. 0. 0.) ;
  134.   TL1 . 'RAIDEUR'      = 1.;
  135.   TL1 . 'RAYON'        = 1. ;
  136.   TL1 . 'COEFFICIENT_GLISSEMENT' = mu_d;
  137.   TL1 . 'COEFFICIENT_ADHERENCE'  = mu_s;
  138.   TL1 . 'RAIDEUR_TANGENTIELLE'    = -1.*Kt;
  139. *  TL1 . 'AMORTISSEMENT_TANGENTIEL'= Ct;
  140.   TL1 . 'AMORTISSEMENT_TANGENTIEL'= 1.E-4;
  141.   TL1 . 'VITESSE_ENTRAINEMENT'    = 0. ;
  142. *  TL1 . 'VITESSE_ENTRAINEMENT'    = -1.*R*Omega;
  143.  
  144. * stockage  
  145.   TLB  = TABL 'LIAISON_B';
  146.   TLB . 1 = TL1;
  147.   TLIA = TABL 'LIAISON';
  148.   TLIA . 'LIAISON_B' = TLB; 
  149.  
  150.  
  151. ************************************************************************
  152. *                           CHARGEMENT
  153. ************************************************************************
  154.  
  155.   Fbal = ebal*m/(m*jeu) ;
  156.   Fchpo_Y  = MANU 'CHPO' 1 palfa1 'FALF' Fbal;
  157.   Fchpo_Z  = MANU 'CHPO' 1 palfa2 'FALF' Fbal;
  158. *+++ BALOURD == 1 pour ajouter w**2 a Fbal  
  159.  
  160.  
  161. *********************************************************************
  162. *     CALCUL NON LINEAIRE
  163. *********************************************************************
  164.  
  165. * chargement spatial : GFC1
  166.       TAB_CHAR = TABLE 'CHARGEMENT';
  167.       TAB_CHAR .  1 =  Fchpo_Y;
  168.       TAB_CHAR . -1 =  Fchpo_Z;
  169. *Pour ajouter w**2 dans la definition du chargement      
  170.       TAB_CHAR . 'BALOURD' = 1;
  171.  
  172. * Approximation pour la solution initiale à fréquence fixe       
  173.       VEC_INIT = TABLE 'INITIAL';
  174.       VEC_INIT . 'FREQUENCE' = 0.5/wchoc ;
  175. *      VEC_INIT . 3 = 0.1;
  176.  
  177. * Paramètres numériques pour la continuation
  178.       PAR_NUM = TABLE 'PARAMETRES_NUMERIQUES';
  179.       PAR_NUM . 'TYPE'    = MOT 'FORC';
  180.       PAR_NUM . 'VAL_FIN' = 60./wchoc;
  181.       PAR_NUM . 'DS0'     =   0.01;
  182. *      PAR_NUM . 'DSMAX'   =   0.1;
  183.       PAR_NUM . 'DSMAX'   =   0.025;
  184.       PAR_NUM . 'DSMIN'   =  1.E-4;
  185.       PAR_NUM . 'ITERMOY' =    3.4;  
  186.       PAR_NUM . 'ITERMAX'    =     20;
  187.       PAR_NUM . 'ANGLE_MIN'  =     0.;
  188.       PAR_NUM . 'ANGLE_MAX'  =    40.;
  189.       PAR_NUM . 'ISENS'     =   1.;
  190.        PAR_NUM . 'NBPAS'     =    2000;
  191. *      PAR_NUM . 'NBPAS'     =    200;
  192.       PAR_NUM . 'TOLERANCE' = 1.E-7;
  193. *      PAR_NUM . 'CALC_JAC'  =  VRAI;       
  194.       PAR_NUM . 'CALC_JAC'  =  FAUX;       
  195.  
  196. *----------------------------------------------------------------------*
  197. * Appel a l'operateur 
  198. * opti impi 2;
  199.  opti impi 1;
  200.       TCON = DYNC TBAS TAB_CHAR TLIA TAMOR VEC_INIT PAR_NUM
  201.                         NHBM NFFT;
  202.  opti impi 0;
  203. *----------------------------------------------------------------------*   
  204.  
  205.  
  206. ************************************************************************
  207. *  POST-TRAITEMENT
  208. ************************************************************************
  209.  
  210. * COURBE DE REPONSE
  211. * -----------------
  212.  
  213. * frequence
  214.   FREQ = TCON . 'REPONSE' . 'FREQUENCE';
  215. * MAX de la NORME2 des COEFFICIENTS de FOURIER pour chaque freq      
  216.   AMPS = TCON . 'REPONSE' . 'NORME_DEPLACEMENT';
  217. * evolution
  218.   XY = EVOL MANU '\w' FREQ  '|Q|' AMPS; 
  219.  
  220. * stabilite      
  221.   STAB = TCON . 'REPONSE' . 'STABILITE';
  222.   TABSTB = TABLE;
  223.   TABSTB . 'LIGNE_VARIABLE' = TABL ENTI;
  224.   TABSTB . 'LIGNE_VARIABLE' . 1 = STAB;  
  225.  
  226.   TITR 'Jeffcott : courbe de reponse';
  227.   DESS XY TABSTB ;
  228.   DESS XY TABSTB XBOR 0.14 0.16 ;
  229.  
  230.  
  231. * BIFURCATIONS
  232. * ------------
  233.  
  234. * recup
  235.   NBIF = DIME TCON . 'BIFURCATION' . 'TYPE';
  236.   prLP_Q = PROG;  prLP_w = PROG;
  237.   prBP_Q = PROG;  prBP_w = PROG;
  238.   prPD_Q = PROG;  prPD_w = PROG;
  239.   prNS_Q = PROG;  prNS_w = PROG;
  240.   REPE BBIF NBIF;
  241.     typbif = EXTR TCON . 'BIFURCATION' . 'TYPE'               &BBIF;
  242.     Qbif   = EXTR TCON . 'BIFURCATION' . 'NORME_DEPLACEMENT'  &BBIF;
  243.     wbif   = EXTR TCON . 'BIFURCATION' . 'FREQUENCE'          &BBIF;
  244.     SI (EGA typbif 'LP'); prLP_Q = prLP_Q et Qbif; prLP_w = prLP_w et wbif; FINSI;
  245.     SI (EGA typbif 'BP'); prBP_Q = prBP_Q et Qbif; prBP_w = prBP_w et wbif; FINSI;
  246.     SI (EGA typbif 'PD'); prPD_Q = prPD_Q et Qbif; prPD_w = prPD_w et wbif; FINSI;
  247.     SI (EGA typbif 'NS'); prNS_Q = prNS_Q et Qbif; prNS_w = prNS_w et wbif; FINSI;
  248.   FIN  BBIF;
  249.  
  250. * creation d'1 evolution + table de dessin
  251.   ibif  = 1;
  252.   evBIF = VIDE 'EVOLUTIO';
  253. * Limit Point
  254.   si ((dime prLP_w) > 0);  ibif = ibif + 1;
  255.     evLP = EVOL MANU LEGE 'LP' '\w' prLP_w '|Q|' prLP_Q;
  256.     evBIF = evBIF et evLP;
  257.     TABSTB . ibif = MOT 'NOLI MARQ LOSA';
  258.   finsi;
  259. * Branch Point
  260.   si ((dime prBP_w) > 0);  ibif = ibif + 1;
  261.     evBP = EVOL MANU LEGE 'BP' '\w' prBP_w '|Q|' prBP_Q;
  262.     evBIF = evBIF et evBP;
  263.     TABSTB . ibif = MOT 'NOLI MARQ CARR';
  264.   finsi;
  265. * Period Doubling Point
  266.   si ((dime prPD_w) > 0);  ibif = ibif + 1;
  267.     evPD = EVOL MANU LEGE 'PD' '\w' prPD_w '|Q|' prPD_Q;
  268.     evBIF = evBIF et evPD;
  269.     TABSTB . ibif = MOT 'NOLI MARQ TRID';
  270.   finsi;
  271. * Neimark Sacker Point
  272.   si ((dime prNS_w) > 0);  ibif = ibif + 1;
  273.     evNS = EVOL MANU LEGE 'NS' '\w' prNS_w '|Q|' prNS_Q;
  274.     evBIF = evBIF et evNS;
  275.     TABSTB . ibif = MOT 'NOLI MARQ TRIU';
  276.   finsi;
  277.  
  278.   DESS (XY ET evBIF )TABSTB 'TITR' 'Jeffcott : courbe de reponse + bifurcations';
  279.  
  280.  
  281. * STABILITE via les EXPOSANTS DE FLOQUET
  282. * --------------------------------------
  283.  
  284. * On assigne une couleur aux exposants de FLOQUET et on leur assigne une couleur
  285.   TEXP_R = TCON . 'REPONSE' . 'EXPOSANT_REEL';
  286.   TEXP_I = TCON . 'REPONSE' . 'EXPOSANT_IMAGINAIRE';
  287.   exp_R = VIDE 'EVOLUTIO';
  288.   exp_I = VIDE 'EVOLUTIO';
  289.    coco = MOTS 'ROUG' 'BLEU' 'VERT' 'AZUR' 'ORAN' 'VIOL' 'TURQ' ;
  290.   repe bexp;
  291.     si (non (exis TEXP_R &bexp)); quit bexp ; finsi;
  292.     coco_i = EXTR coco &bexp;
  293.     exp_R = exp_R ET (evol coco_i 'MANU' '\w' FREQ '\m_{R}' TEXP_R . &bexp);
  294.     exp_I = exp_I ET (evol coco_i 'MANU' '\w' FREQ '\m_{I}' TEXP_I . &bexp);
  295.   fin bexp;
  296.  
  297.   TITRE 'Jeffcott : exposants de Floquet';
  298.   DESS exp_R ;
  299.   DESS exp_I ;
  300.  
  301.  
  302. * REPONSE DE L'HARMONIQUE 1 
  303. * -------------------------
  304. * evolution de l'harmonique 1 des modes en y et en z
  305.   j = 1;
  306.   y1 = TCON . REPONSE . COEFFICIENTS . j . palfa1;
  307.   z1 = TCON . REPONSE . COEFFICIENTS . j . palfa2;
  308.   ev1 = (EVOL 'AZUR' 'MANU' 'LEGE' 'mode 1 (y)' '\w'  FREQ 'Q^{j=1}' y1)
  309.     et (EVOL 'ORAN' 'MANU' 'LEGE' 'mode 2 (Z)' '\w'  FREQ 'Q^{j=1}' Z1);
  310.   DESS ev1 'TITR' 'Jeffcott : Reponse de l harmonique 1 des 2 modes' LEGE SO;
  311.  
  312.  
  313. ************************************************************************
  314. * TEST DE BON FONCTIONNEMENT
  315. ************************************************************************
  316.  
  317. * listreel des erreurs
  318.   prERR = PROG;
  319.  
  320. * amplitude max et w correspondant a la 1ere resonance (harmonique 1)
  321.   wref1 = 0.8860;
  322.   Qref1 = 4.5324;
  323.   imax1 wmax1 Qmax1 = MAXI XY;
  324.   prERR = prERR et (wmax1 - wref1) et (Qmax1 - Qref1);
  325.  
  326. * nombre et/ou type et/ou localisation des bifurcations
  327.   wrefNS = 0.28929;
  328.   QrefNS = 1.5662 ;
  329.   wNS = EXTR prNS_w 1;
  330.   QNS = EXTR prNS_Q 1;
  331.   prERR = prERR et (wNS - wrefNS) et (QNS - QrefNS);
  332.  
  333. * test
  334.   ERRMAX = MAXI prERR 'ABS';
  335.   SI (ERRMAX > 0.05);
  336.     LIST prERR;opti donn 5 trac X;
  337.     ERRE 5;
  338.   FINSI;
  339.  
  340.  
  341.  
  342. ************************************************************************
  343. * PERFORMANCES
  344. ************************************************************************
  345.   TEMP IMPR MAXI CPU;
  346.  
  347. * opti donn 5 trac X;
  348. FIN ; 
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  

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