* fichier : cl_E_4.dgibi
*C***********************************************************************
*C* *
*C* PROJET : Opérateur CLMI *
*C* NOM : cl_E_4 .dgibi *
*C* DESCRIPTION : Jeu de données pour le calcul de la couche limite *
*C* turbulente en utilisant les relations de fermeture de *
*C* Cousteix - Cas de la plaque plane avec gradient de *
*C* pression positif constant -Méthode de Head *
*C* LANGAGE : Gibiane *
*C* AUTEUR : Guillaume VENCO - DRN/DMT/SEMT/LTMF *
*C* *
*C***********************************************************************
*C* *
*C* APPELES : Opérateur CLMI *
*C* *
*C***********************************************************************
*C* *
*C* *
*C***********************************************************************
*C* *
*C* VERSION : 19/05/2000 *
*C* CREATION : 19/04/2000 *
*C* *
*C***********************************************************************
* 19/04/2000 Source: TTMF2 Nom: Cl_E_3.dgibi
*Couche limite plaque plane turbulente
* avec gradient de pression négatif constant
* Fermetures de HEAD
'OPTION' dime 2 elem seg2 ;
GRAPH=FAUX ;
COMPLET=FAUX;
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>>Maillage et domaine>>>>>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
P1 = 0.15 0.;
p2 = 1. 0. ;
nb1 = 10;
p1p2 = 'DROIT' p1 p2 'DINI' 0.0001 'DFIN' 0.01 ;
*p1p2 = 'DROIT' p1 p2 20;
mt = p1p2 ;
* 'TRACER' mt ;
*Définition du modèle
Mmt = 'CHANGER' mt 'QUAF' ;
Mail = 'CHANGER' Mmt 'POI1';
Pt0 = Mail 'ELEM' 1;
$mt = 'MODELISER' Mmt 'NAVIER_STOKES' 'LINE' ;
*Données
Ue0 = 1.;
NU = 1.d-6 ;
*Valeurs initiales données par la loi de puissance
D20=4.33d-4;
Q10=6.37*D20;
*Pas de temps et nombre de pas
DT = 0.1;
NBIT =20;
Si COMPLET;
NBIT =200;
Finsi;
*Champ de vitesse à l'extérieur
*vitesse en puissance
K=1.;
m=(1.);
m1= 'DOMA' $mt 'SOMMET';
X = 'COORDONNEE' 1 m1;
*X1 = (X**m);
*UEX = K*X1;
X1 = ((X-0.15)**m);
UEX = 2.-(K*X1);
UE = 'KCHT' $mt scal 'SOMMET' UEX;
X2 = ((X-0.15)**(m-1));
DUEX = (-1.)*K*m*X2 ;
DUEX=-1.;
DUE = 'KCHT' $mt scal 'SOMMET' DUEX;
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>Création de la table RV et de la table RK>>>>>>
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*Dans CLMI le 1er coef correspond au type de fermeture
*le 2ème coef correspond au n°de l'équation à traiter
*1=QDM 2=Energie cinétique
rv = 'EQEX' 'OMEGA' 1. 'NITER' 1 'ITMA' 0 'FIDT' 1
'OPTI' 'EF' 'IMPL' 'CENTREE' 'BDF2'
'ZONE' $mt
'OPER' 'DFDT' 1. 'D2NM' 'D2N2' DT 'INCO' 'D2' ;
* 'OPER' 'DFDT' 1. 'D2NM' DT 'INCO' 'D2' ;
rv = 'EQEX' rv
'OPTI' 'EF' 'IMPL' 'SUPG'
'ZONE' $mt
'OPER' 'CLMI' 5 1 UE DUE 'D2NM' 'Q1NM'
'INCO' 'D2' ;
rv = 'EQEX' rv
'OPTI' 'EF' 'IMPL' 'CENTREE' 'BDF2'
'ZONE' $mt
'OPER' 'DFDT' 1. 'Q1NM' 'Q1N2' DT 'INCO' 'Q1' ;
* 'OPER' 'DFDT' 1. 'D1NM' DT 'INCO' 'D1' ;
rv = 'EQEX' rv
'OPTI' 'EF' 'IMPL' 'SUPG'
'ZONE' $mt
'OPER' 'CLMI' 5 3 UE DUE 'Q1NM' 'D2NM'
'INCO' 'Q1' ;
rv = 'EQEX' rv
'CLIM' 'D2' 'TIMP' pt0 D20
'CLIM' 'Q1' 'TIMP' pt0 Q10;
rv.inco = 'TABLE' inco ;
rv.inco.'D2'= 'KCHT' $mt scal sommet D20;
rv.inco.'D2NM'= 'KCHT' $mt scal sommet D20;
rv.inco.'D2N2'= 'KCHT' $mt scal sommet D20;
rv.inco.'Q1'= 'KCHT' $mt scal sommet Q10;
rv.inco.'Q1NM'= 'KCHT' $mt scal sommet Q10;
rv.inco.'Q1N2'= 'KCHT' $mt scal sommet Q10;
ltps = 'PROG' 0. ;
lteta = 'PROG' 0. ;
TPS = 0. ;
nupt = 0 ;
mts1 = 'DOMA' $mt 'SOMMET' ;
pt1 = 'POIN' mts1 2 ;
'LISTE' pt1 ;
*'OPTION' donn 5 ;
'REPETER' BCLT nbit ;
nupt=nupt+1;
tps=tps+dt ;
mess ' Pas de temps DT=' dt ' Temps = ' tps ;
ltps=ltps et (prog tps) ;
EXEC rv ;
teta='EXTRAIRE' (rv.inco.'D2') 'SCAL' pt1 ;
mess ' Teta=' teta ;
lteta= lteta et ('PROG' teta);
YNM = rv.'INCO'.'D2NM';
ZNM = rv.'INCO'.'Q1NM';
rv.'INCO'.'D2N2' = rv.'INCO'.'D2NM';
rv.'INCO'.'D2NM' = rv.'INCO'.'D2';
rv.'INCO'.'Q1N2' = rv.'INCO'.'Q1NM';
rv.'INCO'.'Q1NM' = rv.'INCO'.'Q1';
'FIN' bclt ;
li4 = 'PROG' 'EXPO' 'A' -1. 'B' 0. ltps;
evt2 = 'EVOL' 'ROUGE' manu ltps li4 ;
evt= evol manu ltps lteta;
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>>>>>>>>>>>> résultats>>>>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*évolution de la vitesse extérieure et du gradient
vit = 'EVOL' 'CHPO' UE mt;
gvit = 'EVOL' 'ROUG' 'CHPO' DUE mt;
*quantité de mouvement
D2S = (rv.inco.'D2');
RD2S = UE*D2S/NU;
ev1s = 'EVOL' 'ROUGE' 'CHPO' RD2S mt;
*différence entre (d-d1)
Q1S = (rv.inco.'Q1');
RQ1S = UE*Q1S/NU;
ev2s = 'EVOL' 'ROUGE' 'CHPO' RQ1S mt;
*Epaisseur de déplacement
D1S = (rv.'INCO'.'D1');
RD1S = UE*D1S/NU;
ev3s = 'EVOL' 'ROUGE' 'CHPO' RD1S mt;
*évolution du frottement à la paroi
CfS = rv.'INCO'.'CF';
ev4s = 'EVOL' 'ROUGE' 'CHPO' CfS mt;
*évolution du facteur de forme H final
HS = rv.'INCO'.'H';
ev5s = 'EVOL' 'ROUGE''CHPO' HS mt;
*évolution du facteur H*
HHS = Q1S*(D2S**(-1));
ev6s = 'EVOL' 'ROUGE' 'CHPO' HHS mt;
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*solution théorique cf:SCHLICHTING >>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
RX = X*UE*(1./NU);
*Quantité de mvt
*----Rteta donné par Nikuradze
RD2N = (0.01738/1.304)*(RX**(0.861));
ev1n = 'EVOL' 'BLEU' 'CHPO' RD2N mt;
*----Rteta par la loi de puissance
RD2P = (1./1.e-6)*(0.036*X*(RX**(-1./5)))*UE;
ev1p = 'EVOL' 'VERT' 'CHPO' RD2P mt;
*différence (d-d1) théorique
Q1P=(0.37-(1.2857*0.036))*(X*(RX**(-1./5)));
RQ1P=UE*Q1P/NU;
ev2p='EVOL' 'VERT' 'CHPO' RQ1P mt;
*Epaisseur de déplacement
*----Rdelta1 par la Nikuradze
RD1N = 1.304*RD2N;
ev3n = 'EVOL' 'BLEU' 'CHPO' RD1N mt;
*----Rdelta1 par la loi de puissance
RD1P = (9./7.)*RD2P;
ev3p = 'EVOL' 'VERT' 'CHPO' RD1P mt;
*Coefficient de frottement
*----Formule par Nikuradze
CfN = 0.02296*(RX**(-0.139));
ev4n = 'EVOL' 'BLEU' 'CHPO' CfN mt;
*----Formule par loi de puissance
CfP = 0.0576*(RX**(-1./5.));
ev4p = 'EVOL' 'VERT' 'CHPO' CfP mt;
*facteur de forme
*----H par Nikuradze
HN = 'KCHT' $mt scal sommet 1.304;
ev5n = 'EVOL' 'BLEU' 'CHPO' HN mt;
*----H par loi de puissance
HP = 'KCHT' $mt scal sommet 1.2857;
ev5p = 'EVOL' 'VERT' 'CHPO' HP mt;
*coef H* théorique
HHP='KCHT' $mt scal sommet 7.89;
ev6p= 'EVOL' 'BLEU' 'CHPO' HHP mt;
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*Courbes>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Si GRAPH ;
'DESSIN' (evt) 'TITR' 'Convergence de la valeur au dernier noeud'
'TITY' 'Valeurs de Y' 'TITX' 'Temps';
'DESSIN' (vit 'ET' gvit)
'TITRE' 'Vitesse et gradient de vitesse extérieure';
TAB1 = 'TABLE' ;
TAB1.'TITRE' = 'TABLE' ;
TAB1.'TITRE' . 1 = MOT ' Simulation';
TAB1.'TITRE' . 2 = MOT ' Nikuradze';
TAB1.'TITRE' . 3 = MOT ' Loi de puissance';
TAB2 = 'TABLE' ;
TAB2.'TITRE' = 'TABLE' ;
TAB2.'TITRE' . 1 = MOT ' Simulation';
TAB2.'TITRE' . 2 = MOT ' Loi de puissance';
'DESSIN' (ev1s 'ET' ev1n 'ET' ev1p)
'TITRE' 'Epaisseur de quantité de mouvement'
'TITY' 'Rteta'
'TITX' 'Abscisses'
'GRIL'
'LEGE' TAB1;
'DESSIN' (ev2s 'ET' ev2p) 'TITRE' 'différence entre delta et delta1'
'TITY' 'Rdelta-Rdelta1'
'TITX' 'Abscisses'
'GRIL'
'LEGE' TAB2;
'DESSIN' (ev3s 'ET' ev3n 'ET' ev3p)
'TITRE' 'Epaisseur de déplacement'
'TITY' 'Rdelta1'
'TITX' 'Abscisses'
'GRIL'
'LEGE' TAB1;
'DESSIN' (ev4s 'ET' ev4n 'ET' ev4p)
'TITRE' 'Coefficient de frottement à la paroi'
'TITX' 'Abscisses'
'TITY' 'Cf'
'GRIL'
'LEGE' TAB1;
'DESSIN' (ev5s 'ET' ev5n 'ET' ev5p)
'TITRE' 'Evolution du facteur de forme'
'TITX' 'Abscisses'
'TITY' 'H=delta1/theta'
'LEGE' TAB1;
'DESSIN' (ev6s 'ET' ev6p) 'TITRE' 'Evolution du facteur H*'
'TITX' 'Abscisses'
'TITY' 'H'
'LEGE' TAB1;
Finsi;
*calcul de l'erreur au sens de la norme L2
diag = 'DOMA' $mt 'XXDIAGSI';
r = (RD2S-RD2N)**2;
r2= (RD2S-RD2P)**2;
EL2N = 'SOMT'(diag*r);
EL2P = 'SOMT'(diag*r2);
'MESSAGE' 'Erreur sur teta au sens de la norme L2=';
'LISTE' EL2N;
'LISTE' EL2P;
r3=(RQ1S-RQ1P)**2;
EL3P= 'SOMT'(diag*r3);
'MESSAGE' 'Erreur sur (delta-delta1) au sens de la norme L2=';
'LISTE' EL3P;
*SAUVEGARDE
*'OPTION' 'SAUVER' '/u2/ttmf2/guillaume/c20/turb4H_sauv';
*SAUVER;
**************************************************************************
**************************************************************************
*Facteur de forme
x = 'EXTRAIRE' ev5s 'ABSC';
hs = 'EXTRAIRE' ev5s 'ORDO';
hn = 'EXTRAIRE' ev5n 'ORDO';
hp = 'EXTRAIRE' ev5p 'ORDO';
nlist = 'DIME' x;
*'OPTION' echo 0;
*'OPTION' impr 'turb4H_H.data';
'REPETER' BLIST nlist;
xi = 'EXTRAIRE' x &BLIST;
hsi = 'EXTRAIRE' hs &BLIST;
hni = 'EXTRAIRE' hn &BLIST;
hpi = 'EXTRAIRE' hp &BLIST;
* 'MESSAGE' ('CHAINE' xi ' ' hsi ' ' hni ' ' hpi);
'FIN' BLIST;
*'OPTION' impr 'Bidon';
*Epaisseur de quantité de mouvement
x = 'EXTRAIRE' ev1s 'ABSC';
d2s = 'EXTRAIRE' ev1s 'ORDO';
d2n = 'EXTRAIRE' ev1n 'ORDO';
d2p = 'EXTRAIRE' ev1p 'ORDO';
nlist = 'DIME' x;
*'OPTION' echo 0;
*'OPTION' impr 'turb4H_D2.data';
'REPETER' BLIST nlist;
xi = 'EXTRAIRE' x &BLIST;
d2si = 'EXTRAIRE' d2s &BLIST;
d2ni = 'EXTRAIRE' d2n &BLIST;
d2pi = 'EXTRAIRE' d2p &BLIST;
* 'MESSAGE' ('CHAINE' xi ' ' d2si ' ' d2ni ' ' d2pi);
'FIN' BLIST;
*'OPTION' impr 'Bidon';
*Epaisseur de déplacement
x = 'EXTRAIRE' ev3s 'ABSC';
d1s = 'EXTRAIRE' ev3s 'ORDO';
d1n = 'EXTRAIRE' ev3n 'ORDO';
d1p = 'EXTRAIRE' ev3p 'ORDO';
nlist = 'DIME' x;
*'OPTION' echo 0;
*'OPTION' impr 'turb4H_D1.data';
'REPETER' BLIST nlist;
xi = 'EXTRAIRE' x &BLIST;
d1si = 'EXTRAIRE' d1s &BLIST;
d1ni = 'EXTRAIRE' d1n &BLIST;
d1pi = 'EXTRAIRE' d1p &BLIST;
* 'MESSAGE' ('CHAINE' xi ' ' d1si ' ' d1ni ' ' d1pi);
'FIN' BLIST;
*'OPTION' impr 'Bidon';
*Coeff de frottement
x = 'EXTRAIRE' ev4s 'ABSC';
cfs = 'EXTRAIRE' ev4s 'ORDO';
cfn = 'EXTRAIRE' ev4n 'ORDO';
cfp = 'EXTRAIRE' ev4p 'ORDO';
nlist = 'DIME' x;
*'OPTION' echo 0;
*'OPTION' impr 'turb4H_CF.data';
'REPETER' BLIST nlist;
xi = 'EXTRAIRE' x &BLIST;
cfsi = 'EXTRAIRE' cfs &BLIST;
cfni = 'EXTRAIRE' cfn &BLIST;
cfpi = 'EXTRAIRE' cfp &BLIST;
* 'MESSAGE' ('CHAINE' xi ' ' cfsi ' ' cfni ' ' cfpi);
'FIN' BLIST;
*'OPTION' impr 'Bidon';
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer