Télécharger palier_stationnaire_coq4.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : palier_stationnaire_coq4.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. *
  6. * Mots-clés : machines tournantes, palier, hydrodynamique, lubrification
  7. * equation de Reynolds
  8. * ETUDE DU CHAMP DE PRESSION D'UNE LAME FLUIDE
  9. * ENTRE 2 CYLINDRES CONCENTRIQUES par analogie avec la thermique
  10. * Creation : BP, 10/06/2013
  11. *
  12. ************************************************************************
  13.  
  14. *** OPTIONS GENERALES **************************************************
  15.  
  16. COMPLET = faux;
  17. GRAPH = faux;
  18. OPTI DIME 3 ELEM QUA4;
  19. opti trac PSC EPTR 5 POTR 'HELVETICA_16';
  20. eye1 = -10. -50. 30.;
  21. OEIL = eye1;
  22.  
  23.  
  24. *** PARAMETRES DU PALIER *****************************************
  25.  
  26. * Longueur Rayon intérieur jeu
  27. L = 10.E-3; R = 0.5*L;
  28. C = 0.5E-3;
  29. * on en deduit :
  30. D = 2.*R; LsurD = L / D;
  31.  
  32. * % comportement du fluide
  33. * (masse volumique, viscosité cinematique et dynamique)
  34. rho = 1000.;
  35. nu = 150.E-6;
  36. mu = rho*nu;
  37.  
  38. * discretisation
  39. * discretisation circonferentielle, axiale
  40. * nth = 200 ; nL = 100 ;
  41. nth = 100 ; nL = 40 ;
  42. * nth = 24 ; nL = 10 ;
  43.  
  44. * densites associees
  45. dth = (2.*pi*R) / nth; dL = L / nL;
  46. dx = (dth*dL)**0.5;
  47. delim = 1.E-4 * (mini (prog dth dL));
  48. DENS dx;
  49.  
  50.  
  51. *** MAILLAGE 2D **************************
  52.  
  53. p0 = 0. 0. 0.;
  54. p1 = (2.*pi*R) 0. 0.;
  55. l1 = p0 droi nth p1;
  56. l2 = l1 plus (0. (L/2) 0.);
  57. Sfilm = l1 regl nL l2;
  58. l0 = INVE (Sfilm COTE 4 COUL ORAN);
  59. l2pi = Sfilm COTE 2 COUL JAUN;
  60. trac (Sfilm et l0 et l2pi);
  61. ConSfilm = CONT Sfilm;
  62. Sfilm0 = Sfilm plus (0. 0. 0.);
  63. nbno_Sfilm = NBNO Sfilm;
  64.  
  65. * MODELE *********************
  66.  
  67. MOfilm = Sfilm MODE THERMIQUE ISOTROPE 'COQ4';
  68.  
  69. * CONDITONS DE CHARGEMENT **************************
  70.  
  71. * vitesse de rotation > 0 si rotor
  72. OMEGA = 100. * 2*pi ;
  73. * vitesse de rotation imposée du centre du cylindre interieur
  74. vphi = 0. ;
  75. ve = 0.;
  76. phi = 0.;
  77.  
  78. * excentrements imposés
  79. si COMPLET;
  80. presurC = prog 0.01 0.02 0.05 PAS 0.05 0.9 0.92
  81. 0.94 PAS 0.01 0.96 PAS 0.002
  82. 0.97 PAS 0.01 0.99 0.999;
  83. sino;
  84. presurC = prog 0.01 0.1 PAS 0.1 0.9 0.95 0.99;
  85. fins;
  86. ne = dime presurC;
  87.  
  88. * un trait pour l'affichage
  89. chatrait = chai '--------------------------------------------'
  90. '-------------------------------';
  91. * listreels pour le stockage des resultats
  92. prFX = prog; prFY = prog;
  93. prFXG = prog; prFYG = prog;
  94. prCXX = prog; prKXX = prog;
  95. prCYX = prog; prKYX = prog;
  96. prCXY = prog; prKXY = prog;
  97. prCYY = prog; prKYY = prog;
  98. prCXXG = prog; prKXXG = prog;
  99. prCYXG = prog; prKYXG = prog;
  100. prCXYG = prog; prKXYG = prog;
  101. prCYYG = prog; prKYYG = prog;
  102. pr_errQ = prog ne*0.;
  103.  
  104. ie = 0;
  105. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> boucle sur les valeurs d'excentrements
  106. REPE BOUE ne; ie = ie + 1 ;
  107.  
  108. esurC = extr presurC ie;
  109. * cha_e = extr mots_e ie;
  110. e = esurC * C ;
  111. chacha = chai '======================================== ' ie ' / ' ne
  112. ' ==== e/C=' FORMAT '(F6.3)' esurC;
  113. MESS chacha;
  114.  
  115.  
  116. * QUELQUES CHPOINT UTILES ****************
  117.  
  118. un0 = MANU 'CHPO' Sfilm 'SCAL' 1.;
  119.  
  120. * theta en degrés
  121. s0 = coor Sfilm 1;
  122. theta0 = s0 * (180./(pi*R)) ;
  123. cth0 = cos theta0;
  124. sth0 = sin theta0;
  125.  
  126. * phi - theta en degrés
  127. phi0 = MANU 'CHPO' Sfilm 'SCAL' phi;
  128. phimth0 = phi0 - theta0;
  129. cphimth0 = cos phimth0;
  130. sphimth0 = sin phimth0;
  131.  
  132. * hauteur du film
  133. h0 = (C*un0) - (e * cphimth0);
  134. h02 = h0**2;
  135.  
  136.  
  137. * QUELQUES CHAML UTILES ****************
  138.  
  139. un_e = MANU 'CHML' Mofilm 'SCAL' 1. 'STRESSES';
  140.  
  141. * theta et (phi - theta) en degrés
  142. theta_e = chan 'CHAM' theta0 MOfilm 'STRESSES';
  143. cth_e = cos theta_e;
  144. sth_e = sin theta_e;
  145. phimth_e = chan 'CHAM' phimth0 MOfilm 'STRESSES';
  146. c_e = cos phimth_e;
  147. s_e = sin phimth_e;
  148.  
  149. * hauteur du film
  150. h_e = (C*un_e) - (e * c_e) ;
  151. h3_e = h_e**3;
  152.  
  153. * epaisseur de la coque (arbitrairement assez petit mais pas trop)
  154. epfilm = 1.E-4*C;
  155. mess 'epaisseur du film = ' epfilm;
  156.  
  157.  
  158. * MATERIAU *********************
  159.  
  160. * definition de la conductivite = mchaml = h(theta)^3
  161. MAfilm = MATE MOfilm 'K' h3_e 'EPAI' epfilm;
  162. * trac MAfilm MOfilm 'TITR' 'conductivite K (h^{3}) et EPAI (\e_{0})';
  163.  
  164.  
  165. * CALCUL DES MATRICES ********************
  166.  
  167. * conductivity
  168. CONfilm = (1./epfilm) * (COND MOfilm MAfilm);
  169.  
  170. * CL en pression
  171. * CLp0 = BLOQ 'T' (l1 et l2);
  172. CLp0 = BLOQ 'T' l2;
  173. CLppi= RELA 'T' l0 - 'T' l2pi;
  174.  
  175. * * CL en flux pour la symetrie de la solution
  176. * MOsym = MODE l1 THERMIQUE CONVECTION;
  177. * MAsym = MAT MOsym
  178. * CLdp0 = CONV
  179.  
  180. * CL pour n'avoir qu'une seule Temperature
  181. CLsup = RELA 'TSUP' Sfilm - 'T' Sfilm;
  182. CLinf = RELA 'TINF' Sfilm - 'T' Sfilm;
  183.  
  184. * assemblage
  185. CONtot = CONfilm et CLp0 et CLppi
  186. et CLsup et CLinf;
  187.  
  188.  
  189. * 2nd MEMBRE *****************************
  190.  
  191. * 1/ via CHPOINTS
  192. * du a la rotation de l'arbre,
  193. * a la translation et a la rotation du centre de l'arbre
  194. FOMEGA = -6. * mu * OMEGA * e * sphimth0;
  195. FTRANS = -12. * mu * ve * cphimth0;
  196. FROTA = 12. * mu * e * vphi * sphimth0;
  197. * mess (maxi FOMEGA abs) (maxi FTRANS abs) (maxi FROTA abs);
  198. Qtot1 = NOMC (FOMEGA + FTRANS + FROTA) 'Q' 'NATU' 'DISCRET';
  199. Qtot1 = -1. * (FLUX MOfilm Qtot1 'SUPE');
  200.  
  201. * 2/ via MCHAML (permet de tester FLUX)
  202. FOMEGA = -6. * mu * OMEGA * e * s_e;
  203. FTRANS = -12. * mu * ve * c_e;
  204. FROTA = 12. * mu * e * vphi * s_e;
  205. * mess (maxi FOMEGA abs) (maxi FTRANS abs) (maxi FROTA abs);
  206. Qtot2 = NOMC (FOMEGA + FTRANS + FROTA) 'Q' ;
  207. Qtot2 = -1. * (FLUX MOfilm Qtot2 'SUPE');
  208.  
  209. *bp, 2020-03-02 : on utilise l'un des 2 mode de calcul du flux
  210. Qtot = Qtot2;
  211. si GRAPH; TRAC Qtot Sfilm 'TITR' '2nd membre'; finsi;
  212.  
  213. * on verifie que les 2 sont tres proches ici
  214. * ecart en norme infinie :
  215. * deltaQ = MAXI 'ABS' (Qtot1 - Qtot2);
  216. * maxQ = MAXI 'ABS' Qtot;
  217. * REMP pr_errQ ie (deltaQ / maxQ);
  218. * TRAC ((Qtot1 - Qtot2) / maxQ) Sfilm 'TITR' 'Ecart sur le 2nd membre';
  219. * on observe un ecart maxi de 3.8E-03 en theta=0+ et 0- car, pour ces elements, les
  220. * approximations 0+ et 0- ne se compensent pas lors du passage en chpoint qui somme.
  221. * --> suite a l'observation ci-dessus, on utilise plutot la resultante (norme 2) :
  222. deltaQ = (RESU ((Qtot1 - Qtot2)**2))**0.5;
  223. deltaQ = (MAXI 'ABS' deltaQ) / nbno_Sfilm;
  224. REMP pr_errQ ie deltaQ ;
  225.  
  226.  
  227.  
  228. * RESOLUTION *****************************
  229. * regime stationnaire
  230. p_0 = RESO CONtot Qtot;
  231.  
  232. * POST TRAITEMENT ************************
  233.  
  234. si GRAPH;
  235. * tracé de la pression
  236. TRAC eye1 (exco p_0 T) Sfilm ConSfilm 'DIX'
  237. 'TITR' 'Pression (T) du film - calcul Cast3m' ;
  238.  
  239. * evolution de la pression delon theta
  240. p_0theta = extr (EVOL 'ROUG' 'CHPO' p_0 'T' l1) 'COUR' 1;
  241. DESS (p_0theta) MIMA
  242. TITX 'R\q' TITY 'p' 'TITR' 'p(\q)';
  243.  
  244. * evolution de la pression delon Z
  245. * il faut trouver thmax tq pmax = max(p)
  246. ithmax thmax pmax = maxi p_0theta;
  247. lthmax = l0 plus (thmax 0. 0.);
  248. elim lthmax Sfilm (delim);
  249. p_0z = extr (EVOL 'ROUG' 'CHPO' p_0 'T' lthmax) 'COUR' 1;
  250. DESS p_0z
  251. TITX 'Z' TITY 'p' 'TITR' 'p(Z)';
  252. fins;
  253.  
  254. * CALCUL DE LA FORCE RESULTANTE **********
  255.  
  256. * hypothese de film complet
  257. p_0cham = chan 'CHAM' (p_0 exco 'T') MOfilm 'STRESSES';
  258. FX = 2. * (INTG (p_0cham * cth_e) MOfilm);
  259. FY = 2. * (INTG (p_0cham * sth_e) MOfilm);
  260. * mess FX FY;
  261.  
  262. * hypothese de Gumbel (p>0)
  263. * p_0pos = 0.5 *(p_0cham + (abs p_0cham));
  264. masq0 = MASQ p_0cham 'EGSUPE' 0.;
  265. p_0pos = masq0 * p_0cham;
  266. FXG = 2. * (INTG (p_0pos * cth_e) MOfilm);
  267. FYG = 2. * (INTG (p_0pos * sth_e) MOfilm);
  268. * mess FXG FYG;
  269.  
  270.  
  271. * AFFICHAGE DES RESULTATS ************************
  272. OPTI ECHO 0;
  273. mess chatrait;
  274. mess '-------------' ie ' : e/C=' esurC '-------------';
  275. mess 'hyp : film complet & Gumbel ' ;
  276. mess chatrait;
  277. mess 'FX = ' FX ' & ' FXG;
  278. mess 'FY = ' FY ' & ' FYG;
  279. mess chatrait;
  280. OPTI ECHO 1;
  281.  
  282. * STOCKAGE DES RESULTATS ************************
  283.  
  284. * efforts
  285. prFX = prFX et FX;
  286. prFY = prFY et FY;
  287. prFXG = prFXG et FXG;
  288. prFYG = prFYG et FYG;
  289.  
  290.  
  291. fin BOUE ;
  292. *<<<<<<<<<<<<<<<<<<<<<<< fin de boucle sur les valeurs d'excentrements
  293.  
  294.  
  295.  
  296. * opti donn 5 trac X;
  297. * TRACE DES RESULTATS STATIONNAIRES ************************
  298.  
  299. * table de dessin
  300. Tdess1 = tabl;
  301. Tdess1 . 2 = mot 'MARQ TRIU REGU';
  302. Tdess1 . 2 = mot 'MARQ PLUS REGU';
  303. Tdess1 . 3 = mot 'NOLI MARQ L LOSA';
  304. Tdess1 . 'TITRE' = tabl;
  305. Tdess1 . 'TITRE' . 1 = mot 'Film complet' ;
  306. Tdess1 . 'TITRE' . 2 = mot 'Gumbel' ;
  307. Tdess1 . 'TITRE' . 3 = chai '[Frene] L/D=' FORMAT '(F4.1)'LsurD ;
  308.  
  309. * calcul norme de la charge - angle de calage
  310. prFAMP = ((prFX**2) + (prFY**2))**0.5;
  311. prbeta = atg prFY prFX;
  312. prFAMPG = ((prFXG**2) + (prFYG**2))**0.5;
  313. prbetaG = atg prFYG prFXG;
  314.  
  315. * normalisation et comparaison avec Frene
  316. * Frene, p139 L/D = 1/2
  317. e_Frene = prog 0.1 PAS 0.1 0.9 0.95;
  318. si (LsurD ega 0.5);
  319. S_Frene = prog 4.32 2.03 1.21 0.784 0.508 0.318 0.184 0.0912 0.0309
  320. 0.0116;
  321. b_Frene = prog 82 75 68.5 61.53 55 48 41 33 23.5 17.;
  322. finsi;
  323. si (LsurD ega 1.);
  324. S_Frene = prog 1.33 0.631 0.388 0.260 0.178 0.120 0.0776 0.0443
  325. 0.0185 0.00831;
  326. b_Frene = prog 79.5 74 68 62.5 56.5 50.5 44 36 26 19;
  327. finsi;
  328. si (LsurD ega 2.);
  329. S_Frene = prog 0.559 0.271 0.173 0.122 0.0893 0.0654 0.0463 0.0297
  330. 0.0143 0.00707;
  331. b_Frene = prog 75 71 67 62.5 58 52.5 46.5 39 29 21;
  332. finsi;
  333. evAMP_F = EVOL MANU e_Frene (S_Frene**-1);
  334. evbeta_F = EVOL MANU e_Frene (-1.*b_Frene);
  335.  
  336. * evolution + dessin : norme de la charge
  337. evFAMP = EVOL 'BLEU' MANU 'e/C' presurC '|F|' prFAMP;
  338. evFAMPG = EVOL 'ROUG' MANU 'e/C' presurC '|F|' prFAMPG;
  339. DESS (evFAMP et evFAMPG)
  340. GRIL 'POIN' 'GRIS'
  341. TITX 'e/C' POSX 'CENT'
  342. TITY '|F| (N)' POSY 'CENT' LOGY
  343. 'TITR' 'Evolution de la charge admissible - Résultats Cast3M'
  344. Tdess1 LEGE 'NO';
  345. INTGr = mu * L * R * OMEGA * ((R/C)**2) / pi;
  346. DESS ((evFAMP/INTGr) et (evFAMPG/INTGr) et evAMP_F)
  347. GRIL 'POIN' 'GRIS'
  348. TITX 'e/C' POSX 'CENT'
  349. TITY '|F|/F*' POSY 'CENT' LOGY
  350. 'TITR' 'Evolution de la charge admissible - Résultats Cast3M'
  351. Tdess1 LEGE 'NO';
  352.  
  353. * evolution + dessin : angle de calage beta
  354. evbeta = EVOL 'BLEU' MANU 'e/C' presurC '\b' prbeta;
  355. evbetaG = EVOL 'ROUG' MANU 'e/C' presurC '\b' prbetaG;
  356. DESS (evbeta et evbetaG)
  357. GRIL 'POIN' 'GRIS'
  358. TITX 'e/C' POSX 'CENT'
  359. TITY '\b' POSY 'CENT' YBOR -90. 0. YGRA 10.
  360. 'TITR' 'Evolution de l angle de calage - Résultats Cast3M'
  361. Tdess1 LEGE 'NO';
  362. DESS (evbeta et evbetaG et evbeta_F)
  363. GRIL 'POIN' 'GRIS'
  364. TITX 'e/C' POSX 'CENT'
  365. TITY '\b' POSY 'CENT' YBOR -90. 0. YGRA 10.
  366. 'TITR' 'Evolution de l angle de calage - Résultats Cast3M'
  367. Tdess1 LEGE 'NO';
  368. * representation polaire
  369. evXY = EVOL 'BLEU' MANU 'eX' (esurC * (cos prbeta))
  370. 'eY' (esurC * (sin prbeta));
  371. evXYG = EVOL 'ROUG' MANU 'eX' (esurC * (cos prbetaG))
  372. 'eY' (esurC * (sin prbetaG));
  373. DESS (evXY et evXYG)
  374. GRIL 'POIN' 'GRIS' CARR
  375. TITX 'e_{X}' POSX 'CENT' XBOR 0. 1. XGRA 0.2
  376. TITY 'e_{Y}' POSY 'CENT' YBOR -1. 0. YGRA 0.2
  377. 'TITR' 'Evolution de la position stationnaire - Résultats Cast3M'
  378. Tdess1 LEGE 'NO';
  379.  
  380.  
  381. * opti donn 5 trac X;
  382.  
  383. * TEST DE BON FONCTIONNEMENT ******************************************
  384.  
  385. * test de la solution interpolee en 0.5 par comparaison avec solution de [Frene]
  386. Fref = ipol 0.5 evAMP_F;
  387. FG = ipol 0.5 (evFAMPG/INTGr) ;
  388. ecart = abs (FG / Fref - 1.);
  389. mess 'ecart=' ecart;
  390. si (ecart > 0.15); ERRE 5; fins;
  391.  
  392. * test de l'ecart entre les 2 modes de FLUX (chpoint vs mchaml)
  393. list pr_errQ;
  394. ecartQ = MAXI 'ABS' pr_errQ;
  395. mess 'ecartQ=' ecartQ;
  396. si (ecartQ > 1.E-10); ERRE 5; fins;
  397.  
  398.  
  399. FIN ;
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  

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