Download palier_stationnaire_coq4.dgibi
* fichier : palier_stationnaire_coq4.dgibi
************************************************************************
*
* 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.);
* 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;
* 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;
ne = dime presurC;
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;
* 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 *****************************
* 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);
Qtot = NOMC (FOMEGA + FTRANS + FROTA) 'Q' 'NATU' 'DISCRET';
Qtot = -1. * (FLUX MOfilm Qtot 'SUPE');
* TRAC Qtot Sfilm 'TITR' '2nd membre';
* 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';
Sommer = mu * L * R * OMEGA * ((R/C)**2) / pi;
DESS ((evFAMP/Sommer) et (evFAMPG/Sommer) 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 ******************************************
Fref = ipol 0.5 evAMP_F;
FG = ipol 0.5 (evFAMPG/Sommer) ;
ecart = abs (FG / Fref - 1.);
si (ecart > 0.15); ERRE 5; fins;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer