* fichier : coneq.dgibi
*=---------------------------------------------------------------------
*= Transport d'un cone par un champ de vitesse à rotationnel constant.
*= Comparaison de schéma en temps EF implicites
*=---------------------------------------------------------------------
*
*------------------
* Options générales
*------------------
*
GRAPH = 'N' ;
ERR1=5.e-5 ;
COMPLET = FAUX ;
SI ( COMPLET ) ;
DIAM=2. ;
SINON ;
DIAM=0.5 ;
FINSI ;
DISCR='QUAF';
'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
'OPTI' 'ISOV' 'SULI' ;
'OPTI' 'ECHO' 1 'TRACER' 'PSC' ;
*
*=========
* MAILLAGE
*=========
*
*
*- Création des points supports des DROITES
*
L = 2.D0 ; LS2 = L / 2.D0 ;
H = L ; HS2 = H / 2.D0 ;
X0 = -1.D0 * LS2 ; X1 = X0 + L ;
Y0 = -1.D0 * HS2 ; Y1 = Y0 + H ;
INUMX = 20 ;
INUMY = INUMX ;
XNUMY = 'FLOT' INUMX ;
INUX1 = INUMX - 1 ;
INUY1 = INUMY - 1 ;
X01 = X0 + X1 * 0.5D0 ;
Y01 = Y0 + Y1 * 0.5D0 ;
DX = X1 - X0 / INUMX ;
DY = Y1 - Y0 / INUMY ;
DX1 = DX / 2.D0 ;
DY1 = DY / 2.D0 ;
*
B1 = X0 Y0 ;
B3 = X1 Y0 ;
H1 = X0 Y1 ;
H3 = X1 Y1 ;
*
*
*- Création des DROITES frontieres
*
DRBAS = B1 'DROI' INUMX B3 ;
DRDRO = B3 'DROI' INUMY H3 ;
DRHAU = H3 'DROI' INUMX H1 ;
DRGAU = H1 'DROI' INUMY B1 ;
DRCOT = DRBAS 'ET' DRDRO 'ET' DRGAU 'ET' DRHAU ;
PELIM = DX1 / (5.D0 * INUMX) ;
*
*- Creation maillage GEOMETRIQUE
*
PTOT1 = 'DALL' DRBAS DRDRO DRHAU DRGAU ;
*
*- Creation maillage HYBRIDE y compris sous-objets (cond. limites)
*
PTOTC = CHAN PTOT1 QUAF ;
MDRCOT = CHAN DRCOT QUAF ;
ELIM (PTOTC et MDRCOT) 1.E-5 ;
HYTOT = 'MODE' PTOTC 'NAVIER_STOKES' DISCR;
$DRCOT= 'MODE' MDRCOT 'NAVIER_STOKES' DISCR ;
DRCOT=doma $drcot 'MAILLAGE' ;
XVOL= DOMA HYTOT 'VOLUME' ;
DOMA HYTOT 'IMPR' ;
ptot1=DOMA HYTOT 'MAILLAGE';
*
*
*- Coordonnées des points par classe de points
*
XS YS = 'COOR' ( DOMA HYTOT 'MAILLAGE') ;
*
* XCEN : Abscisse du sommet du cone
* YCEN : Ordonnée du sommet du cone
* RCEN : Rayon du cone
*
XCEN = 0.D0 ;
YCEN =-0.5D0 ;
RCEN = 3.0D0 * DX ;
*================
* INITIALISATIONS
*================
*
* ------------------------
* = Conditions initiales =
* ------------------------
* h(x,y) = f(r)
* avec f(r) = max( h(1.-r/rcen) , 0.)
* On se sert des égalités suivantes :
* m = min(f,g) = (f+g - abs(f-g)) / 2
* M = max(f,g) = (f+g + abs(f-g)) / 2
*
HAUT = 1.D0 ;
*
*- Charge aux sommets
*
CRS1 = (XS-XCEN)*(XS-XCEN) + ((YS-YCEN)*(YS-YCEN)) ** 0.5D0 ;
CRS2 = 1.D0 - (CRS1 / RCEN) * HAUT ;
CRS3 = 'ABS' CRS2 ;
CRS4 = CRS2 + CRS3 / 2.D0 ;
SI ( COMPLET ) ;
daxh=drbas plus (0. 0.5) ;
daxv=drgau plus (1. 0.) ;
elim (daxh et daxv et ptot1) 1.e-5 ;
SINON ;
daxh=drbas plus (0. 1.) ;
daxv=drgau plus (0.5 0.) ;
elim (daxh et daxv et ptot1) 1.e-5 ;
FINSI ;
*
*- Vitesses aux sommets
*
SPEEX = 'NOMC' 'UX' YS 'NATU' 'DISCRET' ;
SPEEY = 'NOMC' 'UY' XS 'NATU' 'DISCRET' ;
SPEEY = (-1.D0) * SPEEY ;
TN = KCHT HYTOT SCAL SOMMET CRS4 ;
evolth=evol 'CHPO' TN daxh ;
evoltv=evol 'CHPO' TN daxv ;
s2c= SOMT ((NOEL HYTOT TN)*XVOL) ;
mess ' Solution exacte projetée ' s2c ;
DEBPROC TEST ;
ARGU OPDFDT*MOT OPKONV*MOT BETA*FLOTTANT CFL*FLOTTANT
DIAM*FLOTTANT GRAPH*MOT TIT*MOT;
DT=CFL*0.1 ;
NBIT=enti (3.146 * DIAM / DT) + 1 ;
coef = 1.e-20 ;
RV = 'EQEX' HYTOT ITMA NBIT 'FIDT' 1 'DTI' 0.
'OPTI' 'EF' 'SEMI' BETA OPKONV 'MATCONS' 'MTCV'
'ZONE' HYTOT 'OPER' 'KONV' 1. 'UN' COEF DT 'INCO' 'TN'
'CLIM' 'TN' 'TIMP' DRCOT 0. ;
Si ( EGA OPDFDT 'BDF2');
RV = 'EQEX' RV
'OPTI' 'EF' 'CENTREE' BDF2
'ZONE' HYTOT 'OPER' 'DFDT' 1. 'TNM' 'TNMM' DT 'UN' coef 'INCO' 'TN'
;
Sinon ;
RV = 'EQEX' RV
'OPTI' 'EF' OPDFDT
'ZONE' HYTOT 'OPER' 'DFDT' 1. 'TNM' DT 'UN' coef 'INCO' 'TN'
;
FINSI ;
RV.'INCO' = TABLE 'INCO' ;
RV.'INCO'.'TN' = KCHT HYTOT SCAL SOMMET 0. ;
RV.'INCO'.'TNM' = KCHT HYTOT SCAL SOMMET CRS4 ;
RV.'INCO'.'TNMM' = KCHT HYTOT SCAL SOMMET CRS4 ;
RV.'INCO'.'T0' = KCHT HYTOT SCAL SOMMET 0. ;
RV.'INCO'.'UN' = KCHT HYTOT VECT SOMMET SPEEX SPEEY ;
SI ('EGA' graph 'O' );
u = vect RV.'INCO'.'UN' 1. ux uy ;
trac u ptot1 'NCLK';
trac RV.'INCO'.'TNM' ptot1 'NCLK';
Finsi ;
rv.'OMEGA'= 1. ;
rv.'NITER'=1 ;
exec rv ;
maxc=maxi (RV.'INCO'.'TN') ; minc=mini (RV.'INCO'.'TN');
mess 'MAXC=' maxc 'MINC=' minc ;
TN=(RV.'INCO'.'TN') ;
s2c= SOMT ((NOEL HYTOT TN)*XVOL) ;
mess ' S2C= ' s2c ;
SI ('EGA' graph 'O' );
P1 = 'PROG' -5. -5. 2. ;
MONTAGNE (RV.'INCO'.'TN') ptot1 2. P1 'TITRE' tit 'CACHE' ;
'TRACER' RV.'INCO'.'TN' ptot1 ('PROG' -0.1 pas 0.1 0.9) drcot
'TITRE' TIT 'NCLK' ;
FINSI ;
FINPROC RV ;
* SCHEMA CRANK-NICHOLSON GENERALISE GALERKIN (CNG)
CFL=0.3 ;
BETA=1. ;
OPDFDT='CENTREE' ;
OPKONV='CNG' ;
TIT= CHAI 'CNG (1/4 de TOUR)' ;
rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
evolth=evolth et evolh ;
evoltv=evoltv et evolv ;
SI(NON COMPLET ) ;
lrr=prog
0.0000 -2.55800E-04 -1.72779E-03 -2.61253E-03 -2.39863E-03
-7.29778E-04 1.29162E-03 2.72106E-02 0.33912 0.70209
0.92622 0.65468 0.33257 3.10799E-02 -7.50800E-03
-2.72013E-03 -6.27409E-03 -4.72605E-03 -1.68431E-03 4.86326E-03
0.0000;
lr=extr evolv 'ORDO';
list lr ;
ER=resu( abs (lr - lrr) ) ;
mess ' Erreur sur CNG ' er ;
si ( er > 1.e-4) ; erreur 5 ; finsi ;
FINSI ;
* SCHEMA Tenseur Visqueux (TVISQ)
CFL=0.2 ;
BETA=0. ;
OPDFDT='CENTREE' ;
OPKONV='TVISQUEU' ;
TIT= CHAI 'TVISQUEU (1/4 de TOUR)' ;
rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
evolth=evolth et evolh ;
evoltv=evoltv et evolv ;
SI(NON COMPLET ) ;
lrr=prog
0.0000 -6.50725E-04 1.84987E-04 -3.96571E-04 1.01095E-03
-8.97916E-04 -2.64323E-03 1.68876E-02 0.31908 0.66562
0.92950 0.68354 0.35128 3.29776E-02 3.66575E-03
-2.57892E-04 5.48000E-04 1.06891E-03 -4.26246E-04 -8.88270E-04
0.0000;
lr=extr evolv 'ORDO';
list lr ;
ER=resu( abs (lr - lrr) ) ;
mess ' Erreur sur Tenseur Visqueux ' er ;
si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
FINSI ;
* SCHEMA CRANK-NICHOLSON GALERKIN (CN)
CFL=1. ;
BETA=0.5 ;
OPDFDT='CENTREE' ;
OPKONV='CENTREE' ;
TIT= CHAI 'CN (1/4 de TOUR)' ;
rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
evolth=evolth et evolh ;
evoltv=evoltv et evolv ;
SI ( NON COMPLET ) ;
lrr=prog
0.0000 -1.03945E-03 -2.15865E-05 1.46772E-03 -8.37973E-04
-5.20953E-04 9.73026E-03 5.68738E-02 0.24455 0.61662
0.85196 0.80630 0.39733 7.72915E-02 -5.45626E-02
3.51786E-03 2.85817E-03 9.89801E-04 -5.25209E-03 -7.05692E-03
0.0000;
lr=extr evolv 'ORDO';
list lr ;
ER=resu( abs (lr - lrr) ) ;
mess ' Erreur sur CN ' er ;
si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
FINSI ;
* SCHEMA CRANK-NICHOLSON PETROV-GALERKIN (CN-PG)
CFL=1. ;
BETA=0.5 ;
OPDFDT='SUPG' ;
OPKONV='SUPG' ;
TIT= CHAI 'CN-PG (1/4 de TOUR)' ;
rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
evolth=evolth et evolh ;
evoltv=evoltv et evolv ;
SI (NON COMPLET ) ;
lrr=prog
0.0000 2.74117E-08 5.34498E-07 7.58818E-07 2.94414E-05
7.57305E-04 8.80454E-03 5.71758E-02 0.24644 0.61474
0.84870 0.78601 0.39482 8.58106E-02 -4.81594E-02
2.41366E-03 9.25329E-04 -5.22345E-05 -2.55730E-04 -1.64667E-04
0.0000;
lr=extr evolv 'ORDO';
list lr ;
ER=resu( abs (lr - lrr) ) ;
mess ' Erreur sur CN-PG ' er ;
si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
FINSI ;
* SCHEMA EULER IMPLICITE PETOV-GALERKIN (IMPL-PG)
CFL=1. ;
BETA=1. ;
OPDFDT='SUPG' ;
OPKONV='SUPG' ;
TIT= CHAI 'IMPL-PG (1/4 de TOUR)' ;
rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
evolth=evolth et evolh ;
evoltv=evoltv et evolv ;
SI ( NON COMPLET ) ;
lrr=prog
0.0000 -9.04250E-06 6.03244E-05 -1.72068E-04 2.18717E-03
2.76778E-02 8.96770E-02 0.19189 0.32630 0.45198
0.51491 0.48040 0.35676 0.20231 8.34051E-02
1.93513E-02 9.73832E-04 2.47997E-04 4.06084E-05 2.72083E-05
0.0000;
lr=extr evolv 'ORDO';
list lr ;
ER=resu( abs (lr - lrr) ) ;
mess ' Erreur sur IMPL-PG ' er ;
si ( er > 1.e-4) ; erreur 5 ; finsi ;
FINSI ;
* SCHEMA BDF2
CFL=1. ;
BETA=0.5 ;
OPDFDT='BDF2' ;
OPKONV='SUPG' ;
TIT= CHAI 'BDF2 (1/4 de TOUR)' ;
rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
evolth=evolth et evolh ;
evoltv=evoltv et evolv ;
SI(NON COMPLET ) ;
lrr=prog
0.0000 5.49984E-08 2.57798E-07 -9.63347E-07 1.52810E-05
2.97655E-04 4.82244E-03 3.65515E-02 0.16889 0.45454
0.75960 0.84691 0.62439 0.24678 -3.80415E-02
-0.10517 -1.80764E-03 -1.91432E-03 1.90002E-04 4.46618E-04
0.0000;
lr=extr evolv 'ORDO';
list lr ;
ER=resu( abs (lr - lrr) ) ;
mess ' Erreur sur BDF2 ' er ;
si ( er > 1.e-4) ; erreur 5 ; finsi ;
FINSI ;
si ('EGA' graph 'O' );
titre 'Coupe ox à y=1/2 ' ;
TAB1=TABLE;
TAB1.'TITRE'=TABLE ;
TAB1.1='TIRR MARQ CROI REGU';
TAB1.'TITRE' . 1 = mot 'T initial ' ;
TAB1.2='MARQ ETOI REGU';
TAB1.'TITRE' . 2 = mot 'CNG CFL 0.9';
TAB1.3='MARQ ETOI REGU';
TAB1.'TITRE' . 3 = mot 'TVISQ CFL 0.5';
TAB1.4='MARQ CROI REGU ';
TAB1.'TITRE' . 4 = mot 'CN CFL 1. ';
TAB1.5='MARQ LOSA REGU ';
TAB1.'TITRE' . 5 = mot 'CNPG CFL 1.' ;
TAB1.6='MARQ CARR REGU ';
TAB1.'TITRE' . 6 = mot 'IMPLPG CFL 1.' ;
TAB1.7='MARQ PLUS REGU ';
TAB1.'TITRE' . 7 = mot 'BDF2 CFL 1.' ;
DESS EVOLTH 'TITX' 'Ox' LEGE TAB1 ;
titre 'Coupe oy à x=0 ' ;
DESS EVOLTV 'TITX' 'Oy' LEGE TAB1 ;
FINSI ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer