* fichier : shearkei.dgibi
'OPTI' dime 2 elem QUA4;  
'OPTI' 'EPTR' 6 ;

 COMPLET = FAUX ;
*COMPLET = VRAI ;
 GRAPH = FAUX;
*GRAPH = VRAI;


 DISCR = 'LINE';
 DISCR = 'QUAF';
 DISCR = 'MACRO';
KPRES = 'CENTREP1';
KSUPG = 'SUPG' ;  

Ntauo=1. ;
Tadim = 100.;
Tadim = 20. ;

IKAS ='L1n';
IKAS ='M02n';

IKAS ='M002n';
IKAS ='M002n+';
IKAS ='M002m+';
IKAS ='M002n-';
IKAS ='Q01f';
IKAS ='M01f';
IKAS ='M1n';
IKAS ='Q1n';
IKAS ='Q01f';
IKAS ='M01n';
IKAS ='M05n';

*'OPTI' 'TRAC' 'PSC';
*'OPTI' 'FTRA' (chai '/test2/magn/SHEARKEI1/' ikas '.ps');

Si(EGA IKAS 'M002n');
Ntauo=0.02;Tadim = 80.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M002n+');
Ntauo=0.02;Tadim = 100.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' ;
KEA= 0.03; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M002m+');
Ntauo=0.02;Tadim = 100.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' ;
KEA= 0.015; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M002n-');
Ntauo=0.02;Tadim = 80.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' ;
KEA= 0.003; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M02n');
Ntauo=0.2;Tadim = 100.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M05n');
Ntauo=0.5;Tadim = 100.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M1n');
Ntauo=1.;Tadim = 100.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'Q1n');
Ntauo=1.;Tadim = 80.;DISCR = 'QUAF';
ALGOKEPS= mots 'Nut' ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'Q01n');
Ntauo=0.1;Tadim =20.;DISCR = 'QUAF';
ALGOKEPS= mots 'Nut'           ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
*KEA= 0.03; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'Q01f');
Ntauo=0.1;Tadim =20.;DISCR = 'QUAF';
ALGOKEPS= mots 'Nut' 'Filtre';
Echl=0.01;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
*KEA= 0.03; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M01n');
Ntauo=0.1;Tadim = 100.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut'         ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
*KEA= 0.03; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'M01f');
Ntauo=0.1;Tadim = 100.;DISCR = 'MACRO';
ALGOKEPS= mots 'Nut' 'Filtre';
Echl=0.001;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
*KEA= 0.03; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;

Si(EGA IKAS 'L1n');
KPRES = 'CENTRE';
Ntauo=0.5;Tadim = 60.;DISCR = 'LINE';
ALGOKEPS= mots 'Nut' ;
KE = 0.01; EE = 1.e0 ;
*KE = 0.03;
Finsi ;


Si(NON COMPLET);
 DISCR = 'LINE'; KPRES = 'CENTRE';
 Ntauo=0.1;Tadim = 20.;
 ALGOKEPS= mots 'Nut'  ;
KEA= 0.01; KEB=KEA/2. ; EE = 1.e0 ;
Finsi ;


NITMA =Tadim/Ntauo ENTI ;

*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*
*            CALCUL D'UNE COUCHE DE MELANGE
*            ------------------------------
*

   Ua = 42.8 ;
   Ub = 25.2 ;
   DeltaU= Ua - Ub ;
   DeltaU2= DeltaU * DeltaU ;
   Tauo   = 0.03 / DeltaU ;
   DT = Ntauo * Tauo;
mess '--------------------------------------- ' ;
mess 'Ua=' Ua '  Ub=' Ub  'DELTAU=' DELTAU
     'DELTAU2=' DELTAU2 ' Tauo=' Tauo           ;
mess ' Tauo = ' Tauo ' DT =  ' DT ;
mess ' NITMA= ' NITMA                           ;
mess '--------------------------------------- ' ;
*-------------------------------------------------------------------


LX = 0.75;
LY0= 0.18;
LY0= 0.09;
LY =LY0*2.;

p0=0. 0. ;  
p1=0. LY0;
p2=0. ((-1.)*LY0);
q0=LX 0.;
q1=LX LY0;
q2=LX ((-1.)*LY0);

dn1=0.03*LY0;   dn2=0.3*LY0;

ent1=p0 d dini dn1 dfin dn2 p1 ;  
ent2 = ent1 syme 'POINT' p0;  
m1= ent1 trans dini dn1 dfin dn2 (q1 moins p1);  
m2= m1 syme 'DROITE' p0 q0 ;
mf = m1 et m2;
elim mf 1.e-5;
*trace mf ;  
*
*
haut = cote 2 m1;  
bas = cote 2 m2 ;  
axe = cote 4 m2 ;  
elim (mf et ent2) 1.e-5 ;  
*trace (ent1 et ent2 et haut et bas et axe);  
Mmf   = chan QUAF mf ;  
Ment1 = chan QUAF ent1;  
Ment2 = chan QUAF ent2;  
Mhaut = chan QUAF haut;  
Mbas  = chan QUAF bas ;  
Maxe  = chan QUAF axe;  
elim (Mmf et Ment1 et Ment2 et Mhaut et Mbas et MAxe)1.e-5;  

$mf = mode Mmf 'NAVIER_STOKES' DISCR;  
$ent1=mode Ment1 'NAVIER_STOKES' DISCR;  
ent1 = doma $ent1 maillage ;
$ent2=mode Ment2 'NAVIER_STOKES' DISCR;  
ent2 = doma $ent2 maillage ;

$bas=mode Mbas 'NAVIER_STOKES' DISCR;  
bas = doma $bas maillage ;
$haut=mode Mhaut 'NAVIER_STOKES' DISCR;  
haut= doma $haut maillage ;
$axe=mode Maxe 'NAVIER_STOKES' DISCR;  
axe = doma $axe  maillage inve ;

doma $mf 'IMPR';  
mf = doma $mf maillage;
cmf = cont mf ;

entt=ent1 et ent2 ;
enta = diff (chan ent2 poi1) (chan entt poi1);
entb = diff (chan ent1 poi1) (chan entt poi1);

axe1 = chan POI1 axe ;
cflstab = (coor 1 (poin (elem axe1 (lect 2))1))/Ua;
CFL = DT/cflstab;
CFLI = ENTI CFL ;
mess '--------------------------------------- ' ;
mess '--------------------------------------- ' ;
mess ' CFLstab = ' cflstab ' CFL = ' CFL ;
mess '--------------------------------------- ' ;
mess '--------------------------------------- ' ;

nba = nbno axe1;
tbzone= table ;
nz=0 ;
tbaret= table ;
tbaret.1 = entt;
repeter BCL (nba-1);
pi = poin axe1 (&bcl+1);
entti=entt plus (pi moins p0);
elim (mf et entti) 1.e-5 ;
tbaret.(&bcl+1) = entti;
*trace (axe et entti);
Si ((coor 1 pi) > 0.35);
nz = nz+1;
tbzone.nz= entti;
Finsi ;
FIN BCL ;
xaxe = evol chpo (coor 1 axe) axe;
xaxe = extr xaxe ABSC;

NU =1.E-5;
cnu = 0.09 ;

  MUF  = KCHT $mf SCAL SOMMET  1.e-2;

  RV= EQEX $mf 'ITMA' NITMA 'OMEGA' 1. 'NITER' 1 'FIDT' 20
 'OPTI' 'IMPL' EF  KSUPG  KPRES
 'ZONE'  $mf    'OPER'    'KEPSILON' 1. 'UN' NU DT 'INCO' 'KN' 'EN'
 'ZONE'  $mf    'OPER'    'NS'       1. 'UN' 'MUF' 'INCO' 'UN'
 'OPTI' 'EFM1'
 'ZONE' $mf  'OPER' 'DFDT' 1. 'UN' DT 'UN' 'MUF' 'INCO' 'UN'
   ;         
  RV= EQEX RV
  CLIM 'UN' VIMP (HAUT et BAS et ent1 et ent2 ) 0.
  CLIM 'UN' UIMP ENT1  Ua 'UN' UIMP ENTB  Ub
* CLIM 'UN' VIMP (ENTA et ENT2 et HAUT et BAS) 0.
* CLIM 'KN' TIMP (ENTA et ENT2 et HAUT et BAS) KE
* CLIM 'EN' TIMP (ENTA et ENT2 et HAUT et BAS) EE
* CLIM 'UN' UIMP ENT1  Ua 'UN' UIMP ENT2  Ub
  CLIM 'UN' VIMP (ENTA et ENT2               ) 0.
  CLIM 'KN' TIMP ENT1 KEA 'KN' TIMP ENTB KEB
  CLIM 'EN' TIMP (ENTA et ENT2               ) EE
  ;

  RVP = EQEX
  'OPTI' 'EF' KPRES
  ZONE  $mf  OPER  KBBT  (-1.)     INCO 'UN' 'PRES'
*    CLIM 'PRES' TIMP Sortiep 0.
  ; 
    rvp.'METHINV'.TYPINV=2 ; 
    rvp.'METHINV'.IMPINV=0 ; 
    rvp.'METHINV'.NITMAX=250; 
    rvp.'METHINV'.PRECOND=3 ; 
    rvp.'METHINV'.RESID  =1.e-10 ; 
    rvp.'METHINV' . 'FCPRECT'=300 ; 
    rvp.'METHINV' . 'FCPRECI'=300 ;  

    RV.'PROJ'= RVP ;  

  RV.INCO=TABLE 'INCO' ;  
  RV.INCO.'PRES'=kcht $mf SCAL KPRES 0.;  
  RV.INCO.'UN'=kcht $mf VECT SOMMET (1. 0.);
  RV.INCO.'KN'=kcht $mf SCAL SOMMET   KEB  ;
  RV.INCO.'EN'=kcht $mf SCAL SOMMET   EE   ;  
  RV.INCO.'MUF'=kcht $mf SCAL SOMMET NU   ;
* RV.INCO.'Ro' =kcht $mf SCAL CENTRE  1.   ;
  RV.INCO.'Ro' = 1.   ;
  RV.INCO.'Mu' =kcht $mf SCAL SOMMET   NU   ;
  RV.INCO.'Echl'=Echl;

  RV.'NOMVI' = 'UN' ;  
  RV.'ALGO_KEPSILON'=ALGOKEPS;

 na = 10;
 nna=nba/na ;
 lh= elem axe1 (lect 1 pas na (nna*na+1)) ;
       his = khis 'UN' 1 lh
                  'UN' 2 lh
                  'KN'   lh
                  'EN'   lh  ;  
       his.'KFIH'=1;  
   rv.'HIST'=his ;  

* temp place ;  
*option donn 5 ;
      temp place ;
  exec rv ;
      temp place ;

Si GRAPH;
************** Controle sur les historiques *****************************
 his=rv.'HIST' ;
 dessin his.'TABD' his.'1UN' ;
 dessin his.'TABD' his.'2UN' ;
 dessin his.'TABD' his.'KN' ;
 dessin his.'TABD' his.'EN' ;
************** Controle sur les historiques *****************************
Finsi;



*************************************************************************
un = rv.inco.'UN';  

************** Evolution de l'épaisseur de la couche de mélange **********************
uux = EXCO un 'UX';
dudy = exco (kops 'GRADS' uux $mf) 'UY';

deltaw = prog ;
repeter BCL nba;
ai = tbaret.&bcl;
u1 = redu uux ai;
du1dy = redu dudy ai;
delti = ((maxi u1) - (mini u1))/ (maxi du1dy);
deltaw = deltaw et (prog delti);
FIN BCL ;

evdelta= evol manu xaxe deltaw ;

TAB1 = TABLE ;
tab1.1 = 'MARQ LOSA NOLI' ;
tab1.2 = 'MARQ TRIA REGU' ;
tab1.3 = 'MARQ CROI REGU' ;
tab1.4 = 'MARQ TRIB REGU' ;
tab1.5 = 'MARQ PLUS REGU' ;
tab1.'TITRE'    = TABLE         ;
tab1.'TITRE'. 1 = 'MOT' 'Expérience CEAT' ;
tab1.'TITRE'. 2 = 'MOT' (chai 'CAST3M K-eps: CFL=' CFLI);

   dltexpx= prog 0.55 0.7 1. 1.4 1.8 2.35 2.8 3.5 4.15 4.85 5.5 6.25
  6.9 7.6 8.3 9. 10.1 10.4;
   dltexpx= dltexpx* (0.8/11.1);
   dltexp = prog 1.4 1.5 1.65 1.85 2.05 2.4 2.35 2.7 3.2 3.6 4. 4.35
   4.7 5.3 5.75 6.2 6.95 7.05;
   dltexp = dltexp * (0.05/8.65);

dlexp = prog 1.4 1.5 1.65 1.85 2.05 2.4 2.35 2.7 3.2 3.6 4. 4.35 4.7
             5.3 5.75 6.2 6.95 7.05;
dlexp = dlexp * (0.05/8.65);
xexp  = prog 0.55 0.7 1. 1.4 1.8 2.35 2.8 3.5 4.15 4.85 5.5 6.25 6.9
             7.6 8.3 9. 10.1 10.4 ;
xexp = xexp * (0.8/11.1);
 edltexp = evol manu xexp dlexp;
evd = edltexp et evdelta;
******** FIN * Evolution de l'épaisseur de la couche de mélange **********************

* Graphiques
Si GRAPH;

dess evd TITR (chai ' Delta_w (m) Tadim=' Tadim ' CFL=' CFLI)
 TITX ' X (m) ' TITY 'Delta_w (m)' LEGE tab1;

kn = rv.inco.'KN';
ung = vect un 0.00055 ux uy jaune ;
trace ung mf ;  

*trace rv.inco.'Lw' mf cmf (prog 1.e-2 pas 0.5 10.) TITR 'Lw' ;

trace (exco un 'UX') mf TITR ' Composante UX de la vitesse';

trace (exco un 'UY') mf cmf TITR ' Composante UY de la vitesse';


************** Autosimilarité du profil transversal de vitesse ***********************
ysec = prog ;
delt = prog ;
nbz = dime tbzone;
nbt = dime tbaret;

teta= rv.inco.'TKTE';
Ksi = rv.inco.'Ksi';
unxa=(1./DeltaU)*((exco un 'UX') - Ub);
knxa=(1./DeltaU/DeltaU)*kn;
tetaxa=teta;
ksia=ksi;


repeter BCLZ (dime tbzone);
ar = tbzone.&bclz inve ;
deltawsi=extr (nbt - nbz + &bclz) deltaw   ;
taui    =(extr (nbt - nbz + &bclz) deltaw)/DeltaU   ;

ev1 = evol chpo unxa ar;
evk = evol chpo knxa ar;
evt = evol chpo tetaxa ar;
evksi = evol chpo ksia  ar;
ysec = (extr ev1 'ABSC');
ny= dime ysec ;
ysec0 = prog ny*LY0;
yseci = (ysec - ysec0)/deltawsi;
uxa=extr ev1 'ORDO';
kxa=extr evk 'ORDO';
txa=(extr evt 'ORDO')/(prog ny*taui);
ksixa=extr evksi 'ORDO';
eva = evol manu yseci uxa;
evak= evol manu yseci kxa;
evat= evol manu yseci txa;
evaksi= evol manu yseci ksixa;

Si (EGA &bclz 1);
evux = eva;
evkx = evak;
evtx = evat;
evksix = evaksi;
Sinon ;
evux = evux et eva;
evkx = evkx et evak;
evtx = evtx et evat;
evksix = evksix et evaksi;
Finsi ;
FIN BCLZ;

 kex = prog 5. 13.5 52. 116. 166. 166. 115. 50. 2.;
xkex = prog -37. -29. -21. -12. -3.5 3.5 12. 21. 45.;
 kex = kex * (0.5e-2/183.*7.);
xkex = xkex* (3./114.);
evkex = evol manu xkex kex ;
xuex1= prog 65. 56. 50. 44. 38. 35. 31.5 28.5 25. 22. 19. 16. 13. 10. 6.5
 3. 0. ;
xuex1=xuex1*(-3./115.);
 uex1=prog 181. 181. 181.5 182.5 182. 182. 181.5 181. 175. 172. 166. 156.5
 145. 133.5 119. 108.5 95.;
 uex1 = ((prog 17*183.) - uex1)*(1./183.);
xuex2=  prog 3. 5.5 9.5 12.5 15. 18. 22. 25. 28. 31. 34.5 37. 41. 44.
 47. 50. 53. 56. ;
xuex2=xuex2*(3./115.);
 uex2= prog 80. 70. 54. 41. 32.5 24. 18.5 12.5 10. 6.5 6.5 3. 3. 3.
 1. 1. 1. 1.;
 uex2 = ((prog 18*183.) - uex2)*(1./183.);
uex = uex1 et uex2;
xuex = xuex1 et xuex2;
evuex = evol manu xuex  uex;

TAB1='TABLE' ;
TAB1 . 1 = 'MOT' 'MARQ LOSA NOLI';
*TAB1 . 2 = 'MOT' 'TIRR ';
TAB1 . 'TITRE'     = 'TABLE' ;
TAB1 . 'TITRE' . 1 = 'MOT' 'Expérience CEAT/LEAT';

dess (evuex et evux)
 TITR (chai '(U-Ub)/(Ua-Ub) : Tadim=' Tadim ' CFL=' CFLI)
 TITX '(Y-Y0)/delta_w' TITY '(U-Ub)/(Ua-Ub)' 'LEGE' TAB1 ;

dess (evkex et evkx)
 TITR (chai 'k/(Ua-Ub)**2. : Tadim=' Tadim ' CFL=' CFLI)
 TITX '(Y-Y0)/delta_w' TITY 'k/(Ua-Ub)**2.' 'LEGE' TAB1 ;

 dess evTx  TITR (chai 'TETA : Tadim=' Tadim ' CFL=' CFLI)
 TITX '(Y-Y0)/delta_w' TITY 'TETA' ;

 dess evksix  TITR (chai 'Ksi : Tadim=' Tadim ' CFL=' CFLI)
 TITX '(Y-Y0)/delta_w' TITY 'Ksi' ;

********** FIN Autosimilarité du profil transversal de vitesse ***********************

pe = rv.inco.'PRESSION' ;
pn = elno pe $mf KPRES;
trace pn mf cmf TITR ' Pression' ;


trace kn mf cmf TITR 'KN';
en = rv.inco.'EN';
trace en mf cmf TITR 'EN';
nut = cnu * kn * kn * (inve en) ;
nutsnu= nut* (1./NU);
trace nutsnu mf cmf TITR ' Nut / Nu';

ksi = rv.inco.'Ksi';
tkti= rv.inco.'TKTI' ;
To  = rv.inco.'TKTE' ;
trace tkti mf cmf TITR ' TKTI';
trace To   mf cmf TITR ' To  ';
trace ksi  mf cmf TITR ' ksi ';
Finsi;
* Fin Graphiques

Si (NON COMPLET);
xd= extr edltexp 'ABSC';
yd= extr edltexp 'ORDO';
xc= extr evdelta 'ABSC';
yc= extr evdelta 'ORDO';
list yc;

ycr=prog
    5.41126E-03  5.41313E-03  5.41535E-03  5.41784E-03  5.42038E-03
    5.42327E-03  5.42664E-03  5.43045E-03  5.43462E-03  5.43922E-03
    5.44424E-03  5.44978E-03  5.45584E-03  5.46250E-03  5.46982E-03
    5.47788E-03  5.48679E-03  5.49666E-03  5.50763E-03  5.51989E-03
    5.53364E-03  5.54914E-03  5.56668E-03  5.58681E-03  5.60970E-03
    5.63609E-03  5.66645E-03  5.70168E-03  5.74247E-03  5.79062E-03;
ycr=ycr et (prog
    5.84654E-03  5.91233E-03  5.98995E-03  6.08212E-03  6.18944E-03
    6.31623E-03  6.46651E-03  6.64144E-03  6.84154E-03  7.06852E-03
    7.32907E-03  7.61250E-03  7.90744E-03  8.19769E-03  8.48321E-03
    8.68052E-03  8.83170E-03  8.95027E-03  9.34207E-03  9.83521E-03
    1.03896E-02  1.10250E-02  1.18099E-02  1.26858E-02  1.36029E-02);
ycr=ycr et (prog
    1.44548E-02  1.53243E-02  1.61430E-02  1.68317E-02  1.77894E-02
    1.89391E-02  1.98979E-02  2.07947E-02  2.15067E-02  2.25427E-02
    2.38273E-02  2.52953E-02  2.62061E-02  2.72421E-02  3.12650E-02
    3.40889E-02  3.30795E-02 );

xd= (prog (mini xc)) et xd et (prog (maxi xc));
yd= (prog (mini yd)) et yd et (prog (maxi yd));
yd=yd - (prog 20*0.0106);

*edltexpr= evol manu xd yd;
 edltexpr= evol manu xc ycr;


*evd = edltexpr et evdelta ; dess evd;
m=(somme edltexpr 'ABSO') resu ;
ddeltaw= resu (somme (edltexpr - evdelta) 'ABSO')/m;
mess ' Deltaw : m=' m ' ddeltaw=' ddeltaw;
 Si (ddeltaw > 1.e-3 );ERREUR 5;Finsi;
FINSI;
 

 FIN;