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

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