* fichier : cl_D_3.dgibi
*C***********************************************************************
*C* *
*C* PROJET : Opérateur CLMI *
*C* NOM : cl_D_3 .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 sans gradient de *
*C* pression *
*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 : 02/05/2000 *
*C* CREATION : 06/04/2000 *
*C* *
*C***********************************************************************
'OPTION' dime 2 elem seg2 ;
'OPTION' trace x ;
GRAPH=FAUX ;
COMPLET=FAUX ;
* ****************************
* *** Maillage et Domaine ***
* ****************************
*Maillage
P1 = 0.3 0.;
p2 = 10. 0. ;
nb1 = 10;
p1p2 = 'DROIT' p1 p2 'DINI' 0.01 'DFIN' 1. ;
*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' ;
* ****************************
* * Propriétés phy & Données *
* ****************************
*Données
Ue0 = 1.;
NU = 1.d-6 ;
*Valeurs initiales données par la loi de puissance
D20=8.67d-4;
Q10=6.42*D20;
*Pas de temps et nombre de pas
DT = 1.;
NBIT =25 ;
Si complet;
NBIT =100;
Finsi ;
*Champ de vitesse à l'extérieur
*vitesse en puissance
K=1.;
m=(0.);
m1= 'DOMA' $mt 'SOMMET';
X = 'COORDONNEE' 1 m1;
X1 = (X**m);
UEX = K*X1;
UE = 'KCHT' $mt scal 'SOMMET' UEX;
X2 = (X**(m-1));
DUEX = K*m*X2 ;
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' 4 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' 4 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>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
'OPTION' trace psc;
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;
Si graph ;
*autre dessin
a = 'EXTRAIRE' ev1s 'SCAL';
log = 'LOG' a;
l10= 'LOG' 10;
loga = (1./l10)*log;
b = 'EXTRAIRE' ev5s 'SCAL';
evH = 'EVOL' 'MANU' 'Rtheta' a 'H' b;
'DESSIN' evH
'LOGX'
'TITRE' 'facteur de forme en couche limite turbulente'
'TITX' 'Log(Rtheta)'
'TITY' 'H'
'GRIL';
Finsi ;
'FIN';
© Cast3M 2003 - All rights reserved.
Disclaimer