Télécharger dync02.dgibi

Retour à la liste

Numérotation des lignes :

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

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