* fichier : palier_stationnaire_coq4.dgibi ************************************************************************ * Section : Thermique Hydraulique * Section : Fluides Thermique * Section : Mécanique Dynamique ************************************************************************ ************************************************************************ * * Mots-clés : machines tournantes, palier, hydrodynamique, lubrification * equation de Reynolds * ETUDE DU CHAMP DE PRESSION D'UNE LAME FLUIDE * ENTRE 2 CYLINDRES CONCENTRIQUES par analogie avec la thermique * Creation : BP, 10/06/2013 * ************************************************************************ *** OPTIONS GENERALES ************************************************** COMPLET = faux; GRAPH = faux; OPTI DIME 3 ELEM QUA4; opti trac PSC EPTR 5 POTR 'HELVETICA_16'; eye1 = -10. -50. 30.; OEIL = eye1; *** PARAMETRES DU PALIER ***************************************** * Longueur Rayon intérieur jeu L = 10.E-3; R = 0.5*L; C = 0.5E-3; * on en deduit : D = 2.*R; LsurD = L / D; * % comportement du fluide * (masse volumique, viscosité cinematique et dynamique) rho = 1000.; nu = 150.E-6; mu = rho*nu; * discretisation * discretisation circonferentielle, axiale * nth = 200 ; nL = 100 ; nth = 100 ; nL = 40 ; * nth = 24 ; nL = 10 ; * densites associees dth = (2.*pi*R) / nth; dL = L / nL; dx = (dth*dL)**0.5; delim = 1.E-4 * (mini (prog dth dL)); DENS dx; *** MAILLAGE 2D ************************** p0 = 0. 0. 0.; p1 = (2.*pi*R) 0. 0.; l1 = p0 droi nth p1; l2 = l1 plus (0. (L/2) 0.); Sfilm = l1 regl nL l2; l0 = INVE (Sfilm COTE 4 COUL ORAN); l2pi = Sfilm COTE 2 COUL JAUN; trac (Sfilm et l0 et l2pi); ConSfilm = CONT Sfilm; Sfilm0 = Sfilm plus (0. 0. 0.); nbno_Sfilm = NBNO Sfilm; * MODELE ********************* MOfilm = Sfilm MODE THERMIQUE ISOTROPE 'COQ4'; * CONDITONS DE CHARGEMENT ************************** * vitesse de rotation > 0 si rotor OMEGA = 100. * 2*pi ; * vitesse de rotation imposée du centre du cylindre interieur vphi = 0. ; ve = 0.; phi = 0.; * excentrements imposés si COMPLET; presurC = prog 0.01 0.02 0.05 PAS 0.05 0.9 0.92 0.94 PAS 0.01 0.96 PAS 0.002 0.97 PAS 0.01 0.99 0.999; sino; presurC = prog 0.01 0.1 PAS 0.1 0.9 0.95 0.99; fins; ne = dime presurC; * un trait pour l'affichage chatrait = chai '--------------------------------------------' '-------------------------------'; * listreels pour le stockage des resultats prFX = prog; prFY = prog; prFXG = prog; prFYG = prog; prCXX = prog; prKXX = prog; prCYX = prog; prKYX = prog; prCXY = prog; prKXY = prog; prCYY = prog; prKYY = prog; prCXXG = prog; prKXXG = prog; prCYXG = prog; prKYXG = prog; prCXYG = prog; prKXYG = prog; prCYYG = prog; prKYYG = prog; pr_errQ = prog ne*0.; ie = 0; *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> boucle sur les valeurs d'excentrements REPE BOUE ne; ie = ie + 1 ; esurC = extr presurC ie; * cha_e = extr mots_e ie; e = esurC * C ; SAUT LIGN; chacha = chai '======================================== ' ie ' / ' ne ' ==== e/C=' FORMAT '(F6.3)' esurC; MESS chacha; * QUELQUES CHPOINT UTILES **************** un0 = MANU 'CHPO' Sfilm 'SCAL' 1.; * theta en degrés s0 = coor Sfilm 1; theta0 = s0 * (180./(pi*R)) ; cth0 = cos theta0; sth0 = sin theta0; * phi - theta en degrés phi0 = MANU 'CHPO' Sfilm 'SCAL' phi; phimth0 = phi0 - theta0; cphimth0 = cos phimth0; sphimth0 = sin phimth0; * hauteur du film h0 = (C*un0) - (e * cphimth0); h02 = h0**2; * QUELQUES CHAML UTILES **************** un_e = MANU 'CHML' Mofilm 'SCAL' 1. 'STRESSES'; * theta et (phi - theta) en degrés theta_e = chan 'CHAM' theta0 MOfilm 'STRESSES'; cth_e = cos theta_e; sth_e = sin theta_e; phimth_e = chan 'CHAM' phimth0 MOfilm 'STRESSES'; c_e = cos phimth_e; s_e = sin phimth_e; * hauteur du film h_e = (C*un_e) - (e * c_e) ; h3_e = h_e**3; * epaisseur de la coque (arbitrairement assez petit mais pas trop) epfilm = 1.E-4*C; mess 'epaisseur du film = ' epfilm; * MATERIAU ********************* * definition de la conductivite = mchaml = h(theta)^3 MAfilm = MATE MOfilm 'K' h3_e 'EPAI' epfilm; * trac MAfilm MOfilm 'TITR' 'conductivite K (h^{3}) et EPAI (\e_{0})'; * CALCUL DES MATRICES ******************** * conductivity CONfilm = (1./epfilm) * (COND MOfilm MAfilm); * CL en pression * CLp0 = BLOQ 'T' (l1 et l2); CLp0 = BLOQ 'T' l2; CLppi= RELA 'T' l0 - 'T' l2pi; * * CL en flux pour la symetrie de la solution * MOsym = MODE l1 THERMIQUE CONVECTION; * MAsym = MAT MOsym * CLdp0 = CONV * CL pour n'avoir qu'une seule Temperature CLsup = RELA 'TSUP' Sfilm - 'T' Sfilm; CLinf = RELA 'TINF' Sfilm - 'T' Sfilm; * assemblage CONtot = CONfilm et CLp0 et CLppi et CLsup et CLinf; * 2nd MEMBRE ***************************** * 1/ via CHPOINTS * du a la rotation de l'arbre, * a la translation et a la rotation du centre de l'arbre FOMEGA = -6. * mu * OMEGA * e * sphimth0; FTRANS = -12. * mu * ve * cphimth0; FROTA = 12. * mu * e * vphi * sphimth0; * mess (maxi FOMEGA abs) (maxi FTRANS abs) (maxi FROTA abs); Qtot1 = NOMC (FOMEGA + FTRANS + FROTA) 'Q' 'NATU' 'DISCRET'; Qtot1 = -1. * (FLUX MOfilm Qtot1 'SUPE'); * 2/ via MCHAML (permet de tester FLUX) FOMEGA = -6. * mu * OMEGA * e * s_e; FTRANS = -12. * mu * ve * c_e; FROTA = 12. * mu * e * vphi * s_e; * mess (maxi FOMEGA abs) (maxi FTRANS abs) (maxi FROTA abs); Qtot2 = NOMC (FOMEGA + FTRANS + FROTA) 'Q' ; Qtot2 = -1. * (FLUX MOfilm Qtot2 'SUPE'); *bp, 2020-03-02 : on utilise l'un des 2 mode de calcul du flux Qtot = Qtot2; si GRAPH; TRAC Qtot Sfilm 'TITR' '2nd membre'; finsi; * on verifie que les 2 sont tres proches ici * ecart en norme infinie : * deltaQ = MAXI 'ABS' (Qtot1 - Qtot2); * maxQ = MAXI 'ABS' Qtot; * REMP pr_errQ ie (deltaQ / maxQ); * TRAC ((Qtot1 - Qtot2) / maxQ) Sfilm 'TITR' 'Ecart sur le 2nd membre'; * on observe un ecart maxi de 3.8E-03 en theta=0+ et 0- car, pour ces elements, les * approximations 0+ et 0- ne se compensent pas lors du passage en chpoint qui somme. * --> suite a l'observation ci-dessus, on utilise plutot la resultante (norme 2) : deltaQ = (RESU ((Qtot1 - Qtot2)**2))**0.5; deltaQ = (MAXI 'ABS' deltaQ) / nbno_Sfilm; REMP pr_errQ ie deltaQ ; * RESOLUTION ***************************** * regime stationnaire p_0 = RESO CONtot Qtot; * POST TRAITEMENT ************************ si GRAPH; * tracé de la pression TRAC eye1 (exco p_0 T) Sfilm ConSfilm 'DIX' 'TITR' 'Pression (T) du film - calcul Cast3m' ; * evolution de la pression delon theta p_0theta = extr (EVOL 'ROUG' 'CHPO' p_0 'T' l1) 'COUR' 1; DESS (p_0theta) MIMA TITX 'R\q' TITY 'p' 'TITR' 'p(\q)'; * evolution de la pression delon Z * il faut trouver thmax tq pmax = max(p) ithmax thmax pmax = maxi p_0theta; lthmax = l0 plus (thmax 0. 0.); elim lthmax Sfilm (delim); p_0z = extr (EVOL 'ROUG' 'CHPO' p_0 'T' lthmax) 'COUR' 1; DESS p_0z TITX 'Z' TITY 'p' 'TITR' 'p(Z)'; fins; * CALCUL DE LA FORCE RESULTANTE ********** * hypothese de film complet p_0cham = chan 'CHAM' (p_0 exco 'T') MOfilm 'STRESSES'; FX = 2. * (INTG (p_0cham * cth_e) MOfilm); FY = 2. * (INTG (p_0cham * sth_e) MOfilm); * mess FX FY; * hypothese de Gumbel (p>0) * p_0pos = 0.5 *(p_0cham + (abs p_0cham)); masq0 = MASQ p_0cham 'EGSUPE' 0.; p_0pos = masq0 * p_0cham; FXG = 2. * (INTG (p_0pos * cth_e) MOfilm); FYG = 2. * (INTG (p_0pos * sth_e) MOfilm); * mess FXG FYG; * AFFICHAGE DES RESULTATS ************************ OPTI ECHO 0; saut lign; mess chatrait; mess '-------------' ie ' : e/C=' esurC '-------------'; mess 'hyp : film complet & Gumbel ' ; mess chatrait; mess 'FX = ' FX ' & ' FXG; mess 'FY = ' FY ' & ' FYG; mess chatrait; OPTI ECHO 1; * STOCKAGE DES RESULTATS ************************ * efforts prFX = prFX et FX; prFY = prFY et FY; prFXG = prFXG et FXG; prFYG = prFYG et FYG; fin BOUE ; *<<<<<<<<<<<<<<<<<<<<<<< fin de boucle sur les valeurs d'excentrements * opti donn 5 trac X; * TRACE DES RESULTATS STATIONNAIRES ************************ * table de dessin Tdess1 = tabl; Tdess1 . 2 = mot 'MARQ TRIU REGU'; Tdess1 . 2 = mot 'MARQ PLUS REGU'; Tdess1 . 3 = mot 'NOLI MARQ L LOSA'; Tdess1 . 'TITRE' = tabl; Tdess1 . 'TITRE' . 1 = mot 'Film complet' ; Tdess1 . 'TITRE' . 2 = mot 'Gumbel' ; Tdess1 . 'TITRE' . 3 = chai '[Frene] L/D=' FORMAT '(F4.1)'LsurD ; * calcul norme de la charge - angle de calage prFAMP = ((prFX**2) + (prFY**2))**0.5; prbeta = atg prFY prFX; prFAMPG = ((prFXG**2) + (prFYG**2))**0.5; prbetaG = atg prFYG prFXG; * normalisation et comparaison avec Frene * Frene, p139 L/D = 1/2 e_Frene = prog 0.1 PAS 0.1 0.9 0.95; si (LsurD ega 0.5); S_Frene = prog 4.32 2.03 1.21 0.784 0.508 0.318 0.184 0.0912 0.0309 0.0116; b_Frene = prog 82 75 68.5 61.53 55 48 41 33 23.5 17.; finsi; si (LsurD ega 1.); S_Frene = prog 1.33 0.631 0.388 0.260 0.178 0.120 0.0776 0.0443 0.0185 0.00831; b_Frene = prog 79.5 74 68 62.5 56.5 50.5 44 36 26 19; finsi; si (LsurD ega 2.); S_Frene = prog 0.559 0.271 0.173 0.122 0.0893 0.0654 0.0463 0.0297 0.0143 0.00707; b_Frene = prog 75 71 67 62.5 58 52.5 46.5 39 29 21; finsi; evAMP_F = EVOL MANU e_Frene (S_Frene**-1); evbeta_F = EVOL MANU e_Frene (-1.*b_Frene); * evolution + dessin : norme de la charge evFAMP = EVOL 'BLEU' MANU 'e/C' presurC '|F|' prFAMP; evFAMPG = EVOL 'ROUG' MANU 'e/C' presurC '|F|' prFAMPG; DESS (evFAMP et evFAMPG) GRIL 'POIN' 'GRIS' TITX 'e/C' POSX 'CENT' TITY '|F| (N)' POSY 'CENT' LOGY 'TITR' 'Evolution de la charge admissible - Résultats Cast3M' Tdess1 LEGE 'NO'; INTGr = mu * L * R * OMEGA * ((R/C)**2) / pi; DESS ((evFAMP/INTGr) et (evFAMPG/INTGr) et evAMP_F) GRIL 'POIN' 'GRIS' TITX 'e/C' POSX 'CENT' TITY '|F|/F*' POSY 'CENT' LOGY 'TITR' 'Evolution de la charge admissible - Résultats Cast3M' Tdess1 LEGE 'NO'; * evolution + dessin : angle de calage beta evbeta = EVOL 'BLEU' MANU 'e/C' presurC '\b' prbeta; evbetaG = EVOL 'ROUG' MANU 'e/C' presurC '\b' prbetaG; DESS (evbeta et evbetaG) GRIL 'POIN' 'GRIS' TITX 'e/C' POSX 'CENT' TITY '\b' POSY 'CENT' YBOR -90. 0. YGRA 10. 'TITR' 'Evolution de l angle de calage - Résultats Cast3M' Tdess1 LEGE 'NO'; DESS (evbeta et evbetaG et evbeta_F) GRIL 'POIN' 'GRIS' TITX 'e/C' POSX 'CENT' TITY '\b' POSY 'CENT' YBOR -90. 0. YGRA 10. 'TITR' 'Evolution de l angle de calage - Résultats Cast3M' Tdess1 LEGE 'NO'; * representation polaire evXY = EVOL 'BLEU' MANU 'eX' (esurC * (cos prbeta)) 'eY' (esurC * (sin prbeta)); evXYG = EVOL 'ROUG' MANU 'eX' (esurC * (cos prbetaG)) 'eY' (esurC * (sin prbetaG)); DESS (evXY et evXYG) GRIL 'POIN' 'GRIS' CARR TITX 'e_{X}' POSX 'CENT' XBOR 0. 1. XGRA 0.2 TITY 'e_{Y}' POSY 'CENT' YBOR -1. 0. YGRA 0.2 'TITR' 'Evolution de la position stationnaire - Résultats Cast3M' Tdess1 LEGE 'NO'; * opti donn 5 trac X; * TEST DE BON FONCTIONNEMENT ****************************************** * test de la solution interpolee en 0.5 par comparaison avec solution de [Frene] Fref = ipol 0.5 evAMP_F; FG = ipol 0.5 (evFAMPG/INTGr) ; ecart = abs (FG / Fref - 1.); mess 'ecart=' ecart; si (ecart > 0.15); ERRE 5; fins; * test de l'ecart entre les 2 modes de FLUX (chpoint vs mchaml) list pr_errQ; ecartQ = MAXI 'ABS' pr_errQ; mess 'ecartQ=' ecartQ; si (ecartQ > 1.E-10); ERRE 5; fins; FIN ;