Télécharger dyne06.dgibi

Retour à la liste

Numérotation des lignes :

  1. ************************************************************************
  2. *
  3. * DYNE06 :
  4. *
  5. * Calcul d'un rotor de type Jeffcott avec contact frottant
  6. *
  7. * ___ ^ Z
  8. * _|_ kc,mu |__
  9. * _______ jeu |_ \_
  10. * | k,c | | | \ \
  11. * |----------| m |----------| +--|--|-----> Y
  12. * | |_______| R R+jeu
  13. * ___ jeu
  14. * _|_ kc,mu
  15. *
  16. * ^ Z
  17. * |
  18. * +----> X
  19. *
  20. * Test de l'option VITESSE_ENTRAINEMENT pour la liaison DYNE
  21. * POINT_CERCLE_FROTTEMENT qui ajoute une sur-vitesse tangentielle
  22. *
  23. * Ref : [Xie et al., MSSP, 2015]
  24. *
  25. * Auteur : BP, 2020-03-25
  26. *
  27. ************************************************************************
  28.  
  29.  
  30. ************************************************************************
  31. * OPTIONS et DONNEES
  32. ************************************************************************
  33.  
  34. * options
  35. OPTI DIME 3 ELEM SEG2 EPSI LINEAIRE;
  36. GRAPH = FAUX;
  37. COMPLET = FAUX;
  38. OPTI TRAC PSC EPTR 8 POTR HELVETICA_16;
  39. prefix = chai 'DYNE06';
  40. TITR prefix;
  41.  
  42. * rotor data
  43. m = 1. ;
  44. k = 100.;
  45. c = 5. ;
  46. jeu = 0.105;
  47. R = 20.*jeu ;
  48. ebal = 0.1;
  49. * Fbal = ebal * m * (Omega**2) ; avec Omega = vitesse de rotation en rad/s
  50.  
  51. * contact parameters
  52. kchoc= 25. * k;
  53. mu_s = 0.2;
  54. mu_d = 0.2;
  55. si (mu_s < mu_d); erreur 21; finsi;
  56. Kt = 10.*kchoc;
  57. Ct = 0.5*(m*Kt)**0.5;
  58.  
  59. * deduced parameter
  60. w = (k/m)**0.5;
  61. wHz = w / (2.*pi);
  62. wchoc = (kchoc/m)**0.5;
  63. xi = 0.5 * c / ((kchoc*m)**0.5);
  64. mess 'w=' w '(' wHz ')Hz < wchoc=' wchoc;
  65. mess 'xi=' xi;
  66.  
  67. * critical timestep (Diff Centrees, no damping)
  68. wadhe = ((k+Kt)/m)**0.5;
  69. dtc = 2. / wadhe;
  70.  
  71.  
  72. ************************************************************************
  73. * CALCUL ANALYTIQUE LINEAIRE
  74. ************************************************************************
  75. * si pas de contact, alors pas de couplage entre y et z
  76. * on resout pour y par exemple :
  77. *
  78. * m y" + c y' + ky = meW^2 cos Wt
  79. *
  80. * on pose y(t) = A cos Wt + B sin Wt
  81. * y' = -Aw sin Wt + Bw cos Wt
  82. * y" = -Aw^2 cos Wt - Bw^2 sin Wt
  83. *
  84. * [ k-W^2m Wc ] (A) = (meW^2) termes en cos
  85. * [ -Wc k-W^2m ] (B) = ( 0 ) termes en sin
  86. *
  87. * 2eme ligne : B = Wc/k-W^2m * A
  88. * 1ere ligne : A = meW^2 / [ k-W^2m + W^2c^2/k-W^2m]
  89.  
  90. Wadim = prog 0.02 pas 0.01 0.19 0.198 0.202 0.21 pas 0.01 0.5 pas 0.1 1.2;
  91. Wanal = wchoc * Wadim;
  92. Wanal2 = Wanal**2;
  93. * m=1 -> ca simplifie
  94. aux1 = k - Wanal2;
  95. aux2 = (c**2) * Wanal2;
  96. denom1 = aux1 + (aux2/aux1);
  97. A_anal = (ebal * Wanal2) / denom1;
  98. B_anal = (c*Wanal) / aux1 * A_anal;
  99. * evA = evol VIOL manu Wanal A_anal;
  100. * evB = evol ROSE manu Wanal B_anal;
  101. * dess (evA et evB);
  102. Amp_anal = ( (A_anal**2) + (B_anal**2) )**0.5;
  103. evAmp = evol VIOL manu Wanal Amp_anal;
  104. evjeu = evol manu Wanal (prog (dime Wanal)*jeu);
  105. * dess (evjeu et evAmp);
  106. * en adimensionne
  107. evjeu2 = evol manu Wadim (prog (dime Wanal)*1);
  108. evAmp2 = evol VIOL manu 'w/w_{c}' Wadim 'r/j' (Amp_anal/jeu);
  109. si GRAPH; dess (evjeu2 et evAmp2) TITRE 'solution analytique lineaire'; finsi;
  110. * opti donn 5 trac X;
  111.  
  112.  
  113. ************************************************************************
  114. * BASE MODALE
  115. ************************************************************************
  116.  
  117. * point physique
  118. p1 = 0. 0. 0.;
  119.  
  120. * base modale = 2 modes
  121. palfa1 = 0. 0. 0.;
  122. phi1 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 1. 'UZ' 0. 'NATURE' 'DIFFUS';
  123. palfa2 = 0. 0. 0.;
  124. phi2 = MANU 'CHPO' p1 3 'UX' 0. 'UY' 0. 'UZ' 1. 'NATURE' 'DIFFUS';
  125.  
  126. TBAS = TABL 'BASE_MODALE';
  127. TBAS . 'MODES' = TABL 'BASE_DE_MODES';
  128. TBAS . 'MODES' . 1 = TABL 'MODES';
  129. TBAS . 'MODES' . 1 . 'POINT_REPERE' = palfa1;
  130. TBAS . 'MODES' . 1 . 'FREQUENCE' = wHz;
  131. TBAS . 'MODES' . 1 . 'MASSE_GENERALISEE' = m;
  132. TBAS . 'MODES' . 1 . 'DEFORMEE_MODALE' = phi1;
  133. TBAS . 'MODES' . 2 = TABL 'MODES';
  134. TBAS . 'MODES' . 2 . 'POINT_REPERE' = palfa2;
  135. TBAS . 'MODES' . 2 . 'FREQUENCE' = wHz;
  136. TBAS . 'MODES' . 2 . 'MASSE_GENERALISEE' = m;
  137. TBAS . 'MODES' . 2 . 'DEFORMEE_MODALE' = phi2;
  138.  
  139.  
  140. ************************************************************************
  141. * AMORTISSEMENT
  142. ************************************************************************
  143.  
  144. TAMOR = TABL 'AMORTISSEMENT';
  145. Ctot = MANU 'RIGIDITE' 'TYPE' 'AMORTISSEMENT'
  146. (palfa1 et palfa2) (MOTS 'ALFA') (prog c);
  147. TAMOR . 'AMORTISSEMENT' = Ctot ;
  148.  
  149.  
  150. ************************************************************************
  151. * LIAISONS
  152. ************************************************************************
  153.  
  154. * contact-frottant point_cercle_frottement
  155. TL1 = TABL 'LIAISON_ELEMENTAIRE';
  156. TL1 . 'TYPE_LIAISON' = mot 'POINT_CERCLE_FROTTEMENT';
  157. TL1 . 'SUPPORT' = p1;
  158. TL1 . 'NORMALE' = (1. 0. 0.) ;
  159. TL1 . 'EXCENTRATION' = (0. 0. 0.) ;
  160. TL1 . 'RAIDEUR' = kchoc;
  161. TL1 . 'RAYON' = jeu ;
  162. TL1 . 'COEFFICIENT_GLISSEMENT' = mu_d;
  163. TL1 . 'COEFFICIENT_ADHERENCE' = mu_s;
  164. TL1 . 'RAIDEUR_TANGENTIELLE' = Kt;
  165. TL1 . 'AMORTISSEMENT_TANGENTIEL'= Ct;
  166. * stockage
  167. TLB = TABL 'LIAISON_B';
  168. TLB . 1 = TL1;
  169. TLIA = TABL 'LIAISON';
  170. TLIA . 'LIAISON_B' = TLB;
  171.  
  172. ************************************************************************
  173. * CONDITIONS INITIALES
  174. ************************************************************************
  175.  
  176. * TINIT = TABL 'INITIAL';
  177. * TINIT . 'DEPLACEMENT'
  178. * = (MANU 'CHPO' palfa1 1 'ALFA' ... 'NATURE' 'DISCRET')
  179. * et (MANU 'CHPO' palfa2 1 'ALFA' ... 'NATURE' 'DISCRET');
  180. * TINIT . 'VITESSE'
  181. * = (MANU 'CHPO' palfa1 1 'ALFA' ... 'NATURE' 'DIFFUS')
  182. * et (MANU 'CHPO' palfa2 1 'ALFA' ... 'NATURE' 'DIFFUS');
  183. *u=v=0 ici --> inutile
  184.  
  185. ************************************************************************
  186. * SORTIES
  187. ************************************************************************
  188.  
  189. TSORT = TABLE 'SORTIE';
  190. * +on stocke les deplacements modaux
  191. TSORV = TABLE 'VARIABLE';
  192. TSORV . 'TYPE_SORTIE' = MOT 'LISTREEL';
  193. TSORT . 'VARIABLE' = TSORV;
  194. * + on stocke les variables du contact
  195. TSORLB = TABLE 'LIAISON_B';
  196. TSORLB . TL1 = VRAI;
  197. TSORT . 'LIAISON_B' = TSORLB;
  198.  
  199. * pour le DESSin des evolutions
  200. Tdess1 = TABL ;
  201. Tdess1 . 2 = MOT 'TIRR';
  202. TDSP1 = TABL; TDSP1 . 'TITRE' = TABL;
  203.  
  204.  
  205. ************************************************************************
  206. * PARAMETRES TEMPORELS
  207. ************************************************************************
  208.  
  209. * NOMBRE DE PERIODES SIMULEES, NOMBRE DE PAS PAR PERIODE
  210. * NPERIOD = 8;
  211. *debug NPERIOD = 14;
  212. NPERIOD = 80;
  213. NPPP = 2**10;
  214.  
  215. * NOMBRE DE PAS TOTAL, SORTIE TOUS LES NSORT
  216. NPAS = NPERIOD*NPPP;
  217. NSORT = 8;
  218.  
  219. * NUMERO DES POINTS SORTIS APPARTENANT AUX nlast DERNIERES PERIODES
  220. * nlast = 2;
  221. nlast = 20;
  222. llast = LECT (NPPP / NSORT * (NPERIOD - nlast)) 'PAS' 1 (NPPP / NSORT * NPERIOD + 1);
  223.  
  224.  
  225. *-----------------------------------------------------------------------
  226. *--------------> BOUCLE SUR LA VITESSE DE ROTATION TESTEE
  227. *-----------------------------------------------------------------------
  228.  
  229. * vitesse de rotation (prOmega en rad/s)
  230. si COMPLET;
  231. prWadim = PROG 0.1 0.15 0.154 0.2 0.290 0.3
  232. 0.33 0.4 0.5 1. 1.2 ;
  233. coco_dyne = MOTS 'INDI' 'MARI' 'BLEU' 'AZUR' 'TURQ' 'VERT'
  234. 'BOUT' 'OR' 'ORAN' 'ROUG' 'BRUN' ;
  235. sino;
  236. prWadim = PROG 0.15 0.2 0.33 0.4 1. ;
  237. coco_dyne = MOTS 'MARI' 'AZUR' 'BOUT' 'OR' 'ROUG' ;
  238. finsi;
  239. prOmega = wchoc * prWadim;
  240.  
  241. * RESULTATS
  242. * ymax = prog;
  243. * zmax = prog;
  244. rmax = prog;
  245. rmin = prog;
  246. ev_dsp = VIDE 'EVOLUTIO';
  247.  
  248. REPE BOMEGA (DIME prOmega);
  249. *debug &BOMEGA=1;
  250. *debug &BOMEGA=2;
  251. *debug &BOMEGA=4;
  252.  
  253. coco = EXTR coco_dyne &BOMEGA;
  254. Omega = EXTR prOmega &BOMEGA;
  255. OMEGHZ = Omega / (2.*PI);
  256. MESS '>>> Calcul #' &BOMEGA ' pour OMEGA='Omega 'rad/s (adim ='(Omega/wchoc) ')';
  257. FMT = MOT '(F6.3)';
  258. cha1 = chai '\W=' (enti Omega) ' rad/s - \w/\w_{c}=' FORMAT FMT (Omega/wchoc);
  259. TITRE cha1;
  260. TDSP1 . 'TITRE' . &BOMEGA = chai '\w/\w_{c}='FORMAT FMT (Omega/wchoc);
  261.  
  262. * ON COMPLETE :
  263. *-LIAISON
  264. *** TL1 . 'VITESSE_ENTRAINEMENT' = Omega;
  265. TL1 . 'VITESSE_ENTRAINEMENT' = -1.*R*Omega;
  266. *
  267. *-PERIODE, PAS DE TEMPS, TEMPS FINAL
  268. TPERIOD = 1. / OMEGHZ;
  269. DT = TPERIOD / NPPP;
  270. TFINAL = NPAS*DT;
  271. si (DT < dtc); mess 'Pas de temps =' DT ' < critique=' dtc;
  272. sinon; erre 21;
  273. finsi;
  274.  
  275.  
  276. ************************************************************************
  277. * CHARGEMENT
  278. ************************************************************************
  279.  
  280. * nombre de periodes pour la rampe de chargement
  281. NPrampe = 10;
  282.  
  283. * evolution temporelle du chargement
  284. t_char = PROG (-1.*DT) 'PAS' DT (TFINAL + DT);
  285. rampe = BORN (t_char/(NPrampe*TPERIOD)) 'COMPRIS' 0. 1.;
  286. F_Y = COS (360.*OMEGHZ * t_char);
  287. F_Z = SIN (360.*OMEGHZ * t_char);
  288. *debug F_Y = 1.;
  289. *debug F_Z = 0.;
  290. evF_Y = EVOL 'ROUG' MANU t_char (rampe * F_Y);
  291. evF_Z = EVOL 'BRIQ' MANU t_char (rampe * F_Z);
  292. * DESS (evF_Y et evF_Z) TITRE 'Chargement' XBOR 0. TFINAL XGRA TPERIOD;
  293.  
  294. * description "spatiale" du chargement
  295. Fbal = ebal*m * (Omega**2);
  296. Fchpo_Y = MANU 'CHPO' 1 palfa1 'FALF' Fbal;
  297. Fchpo_Z = MANU 'CHPO' 1 palfa2 'FALF' Fbal;
  298.  
  299. CHAR_Y = CHAR evF_Y Fchpo_Y ;
  300. CHAR_Z = CHAR evF_Z Fchpo_Z ;
  301.  
  302. * table pour DYNE
  303. TCHARG = TABLE 'CHARGEMENT' ;
  304. TCHARG . 'BASE_A' = CHAR_Y et CHAR_Z ;
  305.  
  306.  
  307.  
  308. ************************************************************************
  309. * CALCUL DYNE : INTEGRATION TEMPORELLE
  310. ************************************************************************
  311.  
  312. * TDYNE = DYNE 'DE_VOGELAERE' TBAS TAMOR TLIA TSORT TCHARG
  313. TDYNE = DYNE 'DIFFERENCES_CENTREES' TBAS TAMOR TLIA TSORT TCHARG
  314. NPAS DT NSORT ;
  315.  
  316. ************************************************************************
  317. * POST TRAITEMENT
  318. ************************************************************************
  319.  
  320. * temps t
  321. pr_time = TDYNE . 'TEMPS_DE_SORTIE';
  322. * deplacement modaux
  323. pr_y = TDYNE . 'DEPLACEMENT' . palfa1;
  324. pr_z = TDYNE . 'DEPLACEMENT' . palfa2;
  325.  
  326. * cercle r=R+jeu
  327. pr_360deg = prog 0. PAS 2.5 360.;
  328. ev_butee = EVOL 'DEFA' MANU 'u_{Y}' (jeu*(COS pr_360deg))
  329. 'u_{Z}' (jeu*(SIN pr_360deg));
  330. t_bidon = prog 0. TFINAL;
  331. ev_butee_t = EVOL 'DEFA' MANU 't' t_bidon 'r' (prog 2*jeu);
  332.  
  333. * deplacement u(t) et r(t)
  334. * ev_y = EVOL coco MANU 'LEGE' 'u_{Y}' 't' pr_time 'u' pr_y;
  335. * ev_z = EVOL coco MANU 'LEGE' 'u_{Z}' 't' pr_time 'u' pr_z;
  336. * DESS (ev_y et ev_z) Tdess1;
  337. pr_r = ((pr_y**2)+(pr_z**2))**0.5;
  338. ev_r = EVOL coco MANU 'LEGE' 'r' 't' pr_time 'r' pr_r;
  339. si GRAPH; DESS (ev_butee_t et ev_r); finsi;
  340.  
  341. * orbite complete
  342. ev_xy = EVOL coco MANU 'u_{Y}' pr_y 'u_{Z}' pr_z;
  343. * DESS (ev_xy) 'CARR';
  344. si GRAPH; DESS (ev_butee et ev_xy) 'CARR'; finsi;
  345.  
  346. * restriction aux nlast dernieres periodes
  347. pr_tim2 = extr pr_time llast;
  348. pr_y2 = extr pr_y llast;
  349. pr_z2 = extr pr_z llast;
  350. pr_r2 = ((pr_y2**2)+(pr_z2**2))**0.5;
  351. ev_xy2 = EVOL coco MANU 'u_{Y}' pr_y2 'u_{Z}' pr_z2;
  352. si GRAPH; DESS (ev_butee et ev_xy2) 'CARR'; finsi;
  353. ev_y2 = EVOL coco MANU 't' pr_tim2 'u_{Y}' pr_y2;
  354. * ev_r2 = EVOL coco MANU 'LEGE' 'r' 't' pr_tim2 'r' pr_r2;
  355.  
  356. * DSP
  357. ev_dsp = ev_dsp et (DSPR ev_y2 11 'FMAX' 40. coco);
  358.  
  359. * stockage
  360. * ymax = ymax et (MAXI 'ABS' pr_y2);
  361. * zmax = zmax et (MAXI 'ABS' pr_z2);
  362. rmax = rmax et ( MAXI 'ABS' pr_r2);
  363. rmin = rmin et ( MINI 'ABS' pr_r2);
  364.  
  365. * vitesse tangentielle
  366. pr_vy = TDYNE . 'VITESSE' . palfa1;
  367. pr_vz = TDYNE . 'VITESSE' . palfa2;
  368. ev_v2 = EVOL DEFA MANU LEGE '|v|' 't' pr_time 'v' (((pr_vy**2)+(pr_vz**2))**0.5);
  369. pr_vt = TDYNE . TL1 . 'VITESSE_TANGENTIELLE';
  370. ev_vt = EVOL coco MANU LEGE 'v_{t}' 't' pr_time 'v' pr_vt;
  371. si GRAPH; DESS (ev_vt et ev_v2); finsi;
  372.  
  373. FIN BOMEGA;
  374. *-----------------------------------------------------------------------
  375. *<-------------- FIN DE BOUCLE SUR LA VITESSE DE ROTATION TESTEE
  376. *-----------------------------------------------------------------------
  377.  
  378. * evolution de r avec la vitesse de rotation
  379. ev_rmax = EVOL MANU 'w/w_{c}' prWadim 'max(r)/j' (rmax/jeu);
  380. Tdessr = tabl ; Tdessr . 1 = mot 'NOLI MARQ PLUS';
  381. si GRAPH; dess ev_rmax Tdessr 'TITRE' 'DYNE06'; finsi;
  382.  
  383. * DSP
  384. si GRAPH;
  385. DESS ev_dsp 'LOGY' 'LEGE' TDSP1
  386. 'TITX' 'f(Hz)' 'TITY' 'DSP(u_{Y})' 'TITRE' 'DYNE06';
  387. finsi;
  388.  
  389. ************************************************************************
  390. * TEST DE VALIDATION
  391. ************************************************************************
  392.  
  393. * ON VA CHECKER LE REGIME QUASI-PERIODIQUE A W/WC=0.33
  394. i033 = POSI 0.33 'DANS' prWadim;
  395.  
  396. * recherche du nombre de pics
  397. ev_dsp033 = EXTR ev_dsp 'COUR' i033;
  398. freq033 = EXTR ev_dsp033 'ABSC';
  399. dsp033 = EXTR ev_dsp033 'ORDO';
  400. * recherche des maxi
  401. N033 = dime dsp033;
  402. vari033 = (enle dsp033 N033) - (enle dsp033 1);
  403. chg033 = (enle vari033 (N033 - 1)) * (enle vari033 1);
  404. * imax033 = (POSI 1. 'DANS' (MASQ chg033 INFERIEUR 0.) 'TOUS') + 1;
  405. * on filtre sur les valeurs > 10E-6
  406. masq6 = MASQ (enle (enle dsp033 N033) 1) 'EGSUPE' 1.E-6;
  407. chg033 = masq6 * chg033;
  408. imax033 = (POSI 1. 'DANS' (MASQ chg033 INFERIEUR 0.) 'TOUS') + 1;
  409. list imax033;
  410. w033 = EXTR freq033 imax033;
  411. list w033;
  412. nw033 = DIME w033;
  413.  
  414. * TNR
  415. SI (nw033 > 2);
  416. MESS 'OK: Presence de plusieurs pics a w/wc=0.33 car quasi-periodique';
  417. sinon;
  418. ERREUR 5;
  419. finsi;
  420.  
  421. * recherche des bornes du mvt radial
  422. rmin033 = (extr rmin i033)/jeu;
  423. rmax033 = (extr rmax i033)/jeu;
  424. mess rmin033 rmax033 ;
  425. * valeurs de [Peletan et al., 2014, Nonlinear Dynamics]
  426. rminPel = 0.70;
  427. rmaxPel = 1.45;
  428. errmin = ABS (rmin033 - rminPel);
  429. errmax = ABS (rmax033 - rmaxPel);
  430.  
  431. * TNR
  432. SI ((errmin < 0.1) et (errmax < 0.1));
  433. MESS 'OK: Amplitude du regime quasi-periodique a w/wc=0.33 correctes';
  434. sinon;
  435. ERREUR 5;
  436. finsi;
  437.  
  438.  
  439.  
  440. FIN ;
  441.  
  442.  
  443.  
  444.  
  445.  

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