Accueil | Cast3M

Télécharger 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.);
  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 ;