* fichier : dynasp.dgibi
* Test DYNASP K-epsilon
GRAPH = 'N' ;
complet = FAUX ;
titolo_g = chaine 'Dynasp-35 ' ;
titolo = chaine titolo_g ;
titr titolo ;
opti dime 2 elem qua4 ;
si ( complet ) ;
n1 = 4 ; n2 = 4*n1 ; n3 = -10*n1 ;
sinon ;
n1 = 2 ; n2 = 4*n1 ; n3 = -5*n1 ;
finsi ;
n4 = -1*n1 ; n5 = -1*n2 ;
dmed = 2.000/n3 ; dmin = dmed/1.5 ; dmax = dmed*1.5 ;
pa0 = 0.000 0.000 ; pb0 = 0.0435 0.000 ; pc0 = 0.200 0.000 ;
pa1 = 0.000 1.000 ; pb1 = 0.0435 1.000 ; pc1 = 0.200 1.000 ;
pa2 = 0.000 2.000 ; pb2 = 0.0435 2.000 ; pc2 = 0.200 2.000 ;
pa3 = 0.000 0.500 ; pb3 = 0.0435 0.500 ; pc3 = 0.200 0.500 ;
pa4 = 0.000 1.500 ; pb4 = 0.0435 1.500 ; pc4 = 0.200 1.500 ;
axei = pa0 droi n3 pa1 dini dmax dfin dmin ;
verti = pb0 droi n3 pb1 dini dmax dfin dmin ;
sidei = pc0 droi n3 pc1 dini dmax dfin dmin ;
axes = pa1 droi n3 pa2 dini dmin dfin dmax ;
verts = pb1 droi n3 pb2 dini dmin dfin dmax ;
sides = pc1 droi n3 pc2 dini dmin dfin dmax ;
weti = axei regl n1 verti ;
dryi = verti regl n2 sidei ;
wets = axes regl n1 verts ;
drys = verts regl n2 sides ;
bot = (cote 4 weti) et (cote 4 dryi) ;
hori = (cote 4 wets) et (cote 4 drys) ;
in = (cote 2 wets) ;
dtop = (cote 2 drys) ;
top = in et dtop ;
wet = weti et wets ;
dry = dryi et drys ;
axe = axei et axes ;
side = sidei et sides ;
mt = dry et wet ;
elim mt 0.00001 ;
mt = orie mt ;
cnt = cont mt ;
wall = bot et side et dtop ;
ppa3 = mt poin proc pa3 ;
ppb3 = mt poin proc pb3 ;
ppc3 = mt poin proc pc3 ;
ppa4 = mt poin proc pa4 ;
ppb4 = mt poin proc pb4 ;
ppc4 = mt poin proc pc4 ;
sez3 = ppa3 droi n1 ppb3 droi n2 ppc3 ;
sez4 = ppa4 droi n1 ppb4 droi n2 ppc4 ;
elim (mt et sez3 et sez4) 0.00001 ;
tass mt ;
C = 3.50e-2 ;
Vi = 0.16 ;
Dp = 4.70e-4 ;
S = 0.400*0.087 ;
C_s = C/S ;
* Options
opti isov lign ;
* Reading mesh
yp = 0.010 ;
rapp = (0.400 - yp)/0.400 ;
depl mt affi rapp (0.0 0.0) (1.0 0.0) ;
* trac mt ;
* Dominions
$mt = doma mt impr ;
$bot = doma bot incl $mt 1.0e-4 ;
$axe = doma axe incl $mt 1.0e-4 ;
$side = doma side incl $mt 1.0e-4 ;
$top = doma top incl $mt 1.0e-4 ;
$in = doma in incl $mt 1.0e-4 ;
$cnt = doma cnt incl $mt 1.0e-4 ;
$wet = doma wet incl $mt 1.0e-4 ;
$dry = doma dry incl $mt 1.0e-4 ;
$wall = doma wall incl $mt 1.0e-4 ;
$dtop = doma dtop incl $mt 1.0e-4 ;
volmt = doma $mt volume ;
* Physical properties
condf = 2.62e-2 ;
nuf = 1.50e-5 ;
rof = 1.20e-0 ;
cpf = 1.00e+3 ;
diff = condf/(rof*cpf) ;
nup = 1.00e-20 ;
rop = 2.50e+3 ;
cpp = 8.40e+2 ;
difp = nup ;
rrho = rof/rop ;
qvol = 0.00 ;
gx = 0.00e+0 ;
gy = -9.81e+0 ;
gdr = 0.0 (-9.81*(rop-rof)/rop) ;
zero = 0.0 0.0 ;
beta = 1.0 1.0 ;
* alpha = (3.14/6.0)*n_dot*(Dp**3.0)/Vi ;
alpha = C_s/rop/Vi ;
n_dot = (6.0/3.14)*C_s/(rop*(Dp**3.0)) ;
ustella = 30.0*nuf/yp ;
L0 = 0.30 ;
K0 = 1.0e-4 ;
Rec = 50.0 ;
Uref = 4. ;
* mess 'Cs' C_s 'Dp' Dp 'n_dot' n_dot ; opti donn 5 ;
tabbil = table ;
tabbil.'TF' = mot 'TF' ;
tabbil.'TP' = mot 'TP' ;
tabbil.'VP' = mot 'VN' ;
tabbil.'ALPHA' = mot 'A' ;
tabbil.'ROCPF' = (rof*cpf) ;
tabbil.'ROCPP' = (rop*cpp) ;
tabbil.'VOL' = volmt ;
tabbil.'CNT' = $cnt ;
tabbil.'TOP' = $top ;
tabbil.'BOT' = $bot ;
* Table for bif
tabbif = table ;
tabbif.'RHOF' = rof ;
tabbif.'RHOP' = rop ;
tabbif.'DPART' = Dp ;
tabbif.'NUF' = nuf ;
tabbif.'ALPHA' = mot 'A' ;
tabbif.'UFLUID' = mot 'UN' ;
tabbif.'UPART' = mot 'VN' ;
tabbif.'KFLUID' = mot 'KFL' ;
tabbif.'KPART' = mot 'KPR' ;
* tabbif.'HPART' = mot 'HP' ;
* tabbif.'HFLUID' = mot 'HF' ;
* tabbif.'LAMBDAF' = condf ;
* tabbif.'ROCPP' = (rop*cpp) ;
* tabbif.'ROCPF' = (rof*cpf) ;
* tabbif.'TGASN' = mot 'TF' ;
* tabbif.'TGASE' = mot 'TFE' ;
* tabbif.'TPARTN' = mot 'TP' ;
* tabbif.'TPARTE' = mot 'TPE' ;
* Points for time histories
nc = 10 ;
npti = nbno mt ;
ps = npti / nc ;
pf = ps * nc ;
lh = elem ($mt.sommet) 'POI1' (lect ps pas ps pf) ;
tran = khis
'UN' 1 lh 'UN' 2 lh
'VN' 1 lh 'VN' 2 lh
'KN' lh 'EN' lh
* 'TP' lh 'TF' lh
'A' lh
;
si ( complet ) ;
nbit = 10000 ;
sinon ;
nbit = 800 ;
finsi ;
* Equations + boudary conditions
* Navier-Stokes + K-epsilon + particle mass conservation
rv = EQEX $mt ITMA nbit ALFA 0.70
* ZONE $mt OPER NS nuf zero INCO UN
OPTI 'NOCONS' 'SUPG' 'MMPG'
ZONE $mt OPER NSKE nuf 'NT' zero INCO UN KN EN
ZONE $wall OPER FPU nuf 'USTAR' yp INCO UN KN EN
ZONE $mt OPER NS nup gdr INCO VN
OPTI 'CONS' 'SUPG'
ZONE $mt OPER TSCA nup 'VN' qvol 'TN' INCO A
;
* Coupling terms on q.d.m equations
rv = EQEX rv
OPTI 'NOCONS'
ZONE $mt OPER FROT 'KFL' beta 'VN' INCO UN
* OPTI 'IMPL'
ZONE $mt OPER FROT 'KPR' beta 'UN' INCO VN
;
rv = EQEX rv
ZONE $mt OPER FILTREKE Uref L0 INCO KN EN
ZONE $mt OPER PPRE 'VN' rrho
ZONE $mt OPER BIF tabbif
;
* Coupling terms on energy equations
* rv = EQEX rv
* ZONE $mt OPER TSCA difp 'VN' qvol INCO TP
* ZONE $mt OPER TSCA diff 'UN' qvol INCO TF
* ZONE $mt OPER ECHI 'HP' 'TFE' INCO TP
* ZONE $mt OPER ECHI 'HF' 'TPE' INCO TF
* ZONE $mt OPER BILA tabbil
;
* Boundary coditions
rv = EQEX rv
* CLIM UN UIMP cnt 0.0 UN VIMP side 0.0
CLIM UN VIMP top 0.0 UN VIMP bot 0.0
CLIM UN UIMP axe 0.0 UN UIMP side 0.0
CLIM VN UIMP axe 0.0 VN UIMP side 0.0
CLIM VN VIMP in (-1*Vi) VN VIMP dtop 0.0
CLIM A TIMP in alpha A TIMP dtop 0.0
* CLIM TP TIMP in 20.0
;
* Pressure equation
rp = EQPR $mt KTYPI 1
ZONE $mt OPER PRESSION 0.0 ;
* Initial conditions
alphaw = kcht $wet scal sommet (alpha ) ;
tpwet = kcht $wet scal sommet (20.0 ) ;
rv.titolo = titolo ;
rv.INCO = table INCO ;
rv.INCO.'UN' = kcht $mt vect sommet (1.0e-5 1.0e-5) ;
rv.INCO.'VN' = kcht $mt vect sommet (1.0e-5 (-1*Vi)) ;
rv.INCO.'KN' = kcht $mt scal sommet (1.0e-5 ) ;
rv.INCO.'EN' = kcht $mt scal sommet (1.0e-5 ) ;
rv.INCO.'NT' = kcht $mt scal centre nuf ;
rv.INCO.'A' = kcht $mt scal sommet (0.0 ) alphaw ;
rv.INCO.'TN' = kcht $mt scal sommet (0.0 ) ;
* rv.INCO.'TP' = kcht $mt scal sommet (0.0 ) tpwet ;
* rv.INCO.'TF' = kcht $mt scal sommet (0.0 ) ;
* rv.INCO.'TPE' = kcht $mt scal centre (0.0 ) ;
* rv.INCO.'TFE' = kcht $mt scal centre (0.0 ) ;
rv.INCO.'KFL' = kcht $mt vect centre (0.0 0.0 ) ;
rv.INCO.'KPR' = kcht $mt vect centre (0.0 0.0 ) ;
* rv.INCO.'HF' = kcht $mt scal centre (0.0 ) ;
* rv.INCO.'HP' = kcht $mt scal centre (0.0 ) ;
rv.INCO.'USTAR' = kcht $wall scal centre ustella ;
rv.PRESSION = rp ;
rv.'HIST' = tran ;
rv.'YP' = yp ;
rv.'NUF' = nuf ;
* Executing
exec rv ;
* Saving results
* repix rv ;
* sauv rv ;
* Procedura test
debproc test ;
argu rv*table complet*logique ;
pt1 = mt poin proc (0.0 1.260) ;
pt2 = mt poin proc (0.0 0.740) ;
u2 = exco uy (rv.inco.'UN') uy ;
v2 = exco uy (rv.inco.'VN') uy ;
upt1 = abs (extr (rv.inco.'UN') uy pt1) ;
upt2 = abs (extr (rv.inco.'UN') uy pt2) ;
vpt1 = abs (extr (rv.inco.'VN') uy pt1) ;
vpt2 = abs (extr (rv.inco.'VN') uy pt2) ;
apt1 = abs (extr (rv.inco.'A') scal pt1) ;
apt2 = abs (extr (rv.inco.'A') scal pt2) ;
si ( complet ) ;
upt1_l = 0.53580 ;
upt2_l = 0.59501 ;
vpt1_l = 3.0245 ;
vpt2_l = 3.5610 ;
apt1_l = 1.18535E-03 ;
apt2_l = 1.00603E-03 ;
sinon ;
upt1_l = 0.43961 ;
upt2_l = 0.38635 ;
vpt1_l = 3.0859 ;
vpt2_l = 3.4686 ;
apt1_l = 6.78202E-04 ;
apt2_l = 5.94036E-04 ;
finsi ;
upt1_e = abs ((upt1 - upt1_l)/upt1_l) ;
upt2_e = abs ((upt2 - upt2_l)/upt2_l) ;
vpt1_e = abs ((vpt1 - vpt1_l)/vpt1_l) ;
vpt2_e = abs ((vpt2 - vpt2_l)/vpt2_l) ;
apt1_e = abs ((apt1 - apt1_l)/apt1_l) ;
apt2_e = abs ((apt2 - apt2_l)/apt2_l) ;
erlim = 5.00e-2 ;
si (upt1_e > erlim) ;
mess 'Error on Ugas : err_Upt1= ' upt1_e ' > ' erlim ;
erreur 5 ;
finsi ;
si (upt2_e > erlim) ;
mess 'Error on Ugas : err_Upt2= ' upt2_e ' > ' erlim ;
erreur 5 ;
finsi ;
si (vpt1_e > erlim) ;
mess 'Error on Vpart: err_Vpt1= ' vpt1_e ' > ' erlim ;
erreur 5 ;
finsi ;
si (vpt2_e > erlim) ;
mess 'Error on Vpart: err_Vpt2= ' vpt2_e ' > ' erlim ;
erreur 5 ;
finsi ;
si (apt1_e > erlim) ;
mess 'Error on alpha: err_Apt1= ' apt1_e ' > ' erlim ;
erreur 5 ;
finsi ;
si (apt2_e > erlim) ;
mess 'Error on alpha: err_Apt2= ' apt2_e ' > ' erlim ;
erreur 5 ;
finsi ;
mess 'upt1' upt1 'upt2' upt2 'vpt1' vpt1 'vpt2' vpt2 ;
mess 'apt1' apt1 'apt2' apt2 ;
finproc ;
* Procedura visu
debproc visu rv*table ;
titr rv.titolo ;
tdom = rv.domaine ;
mesh = rv.domaine.maillage ;
cmesh = cont mesh ;
pt1 = poin mesh proc (0.0 0.740) ;
pt2 = poin mesh proc (0.4 0.740) ;
pt3 = poin mesh proc (0.0 1.260) ;
pt4 = poin mesh proc (0.4 1.260) ;
pt12 = poin mesh droit pt2 pt1 0.0001 ;
pt34 = poin mesh droit pt3 pt4 0.0001 ;
meshl = chan mesh lign ;
line1 = elem meshl appuye pt12 ;
line2 = elem meshl appuye pt34 ;
line1 = inve line1 ;
line2 = inve line2 ;
f1 = kcht $wet scal sommet 1.0 ;
filtro = kcht $mt scal sommet 0.0 f1 ;
un = rv.inco.'UN' ;
vn = rv.inco.'VN' ;
vn = vn * filtro ;
a = rv.inco.'A' ;
* tp = rv.inco.'TP' ;
* tf = rv.inco.'TF' ;
u1 = exco ux un ux ;
u2 = exco uy un uy ;
v1 = exco ux vn ux ;
v2 = exco uy vn uy ;
unvt = vect un 5.0e-2 ux uy jaun ;
vnvt = vect vn 5.0e-3 ux uy roug ;
titr (chaine titolo_g 'Uair') ; trac unvt mesh cmesh ;
titr (chaine titolo_g 'Upart') ; trac vnvt mesh cmesh ;
txa = mot 'axe (m)' ;
txo = mot 'X (m)' ;
* txl = mot 'section 0.74 (m)' ;
tyu = mot 'velocity (m/s)' ;
tyua = mot 'air velocity (m/s)' ;
tyup = mot 'u_I (m/s)' ;
tynt = mot 'nu_t (m^2/s)' ;
tya = mot 'alpha' ;
tyt = mot 'temp. (C)' ;
t1 = chaine titolo_g 'U;V /axe ' ;
t2 = chaine titolo_g 'Z=0.74 (m)' ;
t3 = chaine titolo_g 'Z=1.26 (m)' ;
t4 = chaine titolo_g 'Air veloc.' ;
t5 = chaine titolo_g 'Alpha/axe ' ;
t6 = chaine titolo_g 'Alpha ' ;
t7 = chaine titolo_g 'Y+' ;
tab1 = table ;
tab1.1 = 'MARQ CROI REGU ' ;
tab1.2 = 'MARQ LOSA REGU ' ;
tab1.'TITRE' = table ;
tab1.'TITRE'. 1 = mot 'AIR' ;
tab1.'TITRE'. 2 = mot 'PARTICLE' ;
tab2 = table ;
tab2.1 = 'MARQ LOSA REGU ' ;
tab2.2 = 'MARQ TRIA REGU ' ;
tab2.3 = 'MARQ CARR REGU ' ;
tab2.4 = 'MARQ ETOI REGU ' ;
tab2.'TITRE' = table ;
tab2.'TITRE'. 1 = mot 'H=1.50M' ;
tab2.'TITRE'. 2 = mot 'H=1.00M' ;
tab2.'TITRE'. 3 = mot 'H=0.50M' ;
tab2.'TITRE'. 4 = mot 'H=0.74M' ;
tab3 = table ;
tab3.1 = 'MARQ LOSA REGU ' ;
tab3.2 = 'MARQ TRIA REGU ' ;
tab3.'TITRE' = table ;
tab3.'TITRE'. 1 = mot 'H=0.74M' ;
tab3.'TITRE'. 2 = mot 'H=1.26M' ;
u2axe = evol chpo u2 uy axe ;
v2axe = evol chpo v2 uy axe ;
u2top = evol chpo u2 uy top ;
v2top = evol chpo v2 uy top ;
u2int = evol chpo u2 uy (inve hori) ;
v2int = evol chpo v2 uy (inve hori) ;
u2ln1 = evol chpo u2 uy line1 ;
v2ln1 = evol chpo v2 uy line1 ;
u2ln2 = evol chpo u2 uy line2 ;
v2ln2 = evol chpo v2 uy line2 ;
* tpaxe = evol chpo tp scal axe ;
* tfaxe = evol chpo tf scal axe ;
* tpln1 = evol chpo tp scal line1 ;
* tfln1 = evol chpo tf scal line1 ;
aaxe = evol chpo a scal axe ;
atop = evol chpo a scal (top ) ;
aln1 = evol chpo a scal line1 ;
aln2 = evol chpo a scal line2 ;
abot = evol chpo a scal (inve bot ) ;
dess (u2axe et v2axe) mima titx txa tity tyu titr t1 gril lege tab1 ;
* @excel1 u2axe 'pp' ;
dess (u2ln1 et v2ln1) mima titx txo tity tyu titr t2 gril lege tab1 ;
dess (u2ln2 et v2ln2) mima titx txo tity tyu titr t3 gril lege tab1 ;
dess (u2ln1 et u2ln2) mima titx txo tity tyua gril titr t4 lege tab3;
dess (aaxe ) mima titx txa tity tya gril titr t5 ;
dess (aln1 et aln2 ) mima titx txo tity tya gril titr t6 lege tab3;
* dess (tfaxe et tpaxe) mima titx txa tity tyt gril lege tab1 ;
* dess (tfln et tpln ) mima titx txl tity tyt gril lege tab1 ;
si (exist rv 'INCO' 'USTAR') ;
yp = rv.'YP' ;
nuf = rv.'NUF' ;
ue = rv.'INCO'.'USTAR' ;
us = elno $wall ue ;
Yplus = us * yp / nuf ;
Ype = evol chpo Yplus scal wall ;
dess Ype mima tity 'Y+' gril titr t7 ;
finsi ;
si (exist rv 'INCO' 'HP' ) ;
titr (chaine rv.titolo ' Hp' ) ;
trac (elno tdom (rv.'INCO'.'HP' )) mesh cmesh ;
finsi ;
si (exist rv 'INCO' 'HF' ) ;
titr (chaine rv.titolo ' Hf' ) ;
trac (elno tdom (rv.'INCO'.'HF' )) mesh cmesh ;
finsi ;
si (exist rv 'INCO' 'H' ) ;
titr (chaine rv.titolo ' H ' ) ;
trac (rv.'INCO'.'H') mesh cmesh ;
finsi ;
si (exist rv 'INCO' 'TF') ;
titr (chaine rv.titolo ' Tf' ) ;
trac (rv.inco.'TF') mesh cmesh ;
finsi ;
si (exist rv 'INCO' 'TP') ;
titr (chaine rv.titolo ' Tp' ) ;
trac (rv.inco.'TP') mesh cmesh ;
finsi ;
si (exist rv 'INCO' 'A') ;
titr (chaine rv.titolo ' alpha') ;
trac (rv.inco.'A' ) mesh cmesh ;
finsi ;
si (exist rv 'INCO' 'KN') ;
kn = rv.inco.'KN' ;
en = rv.inco.'EN' ;
L0 = (kops (kops kn '**' 1.5) '*' (kops en '**' (-1.0))) ;
titr (chaine rv.titolo ' kn') ;
trac kn mesh cmesh ;
titr (chaine rv.titolo ' en') ;
trac en mesh cmesh ;
mess 'min L0' (mini L0) ; mess 'max L0' (maxi L0) ;
titr (chaine rv.titolo ' L0') ;
trac L0 mesh cmesh ;
* up = kn * 2.0 / 3.0 ;
* up = up ** 0.5 ;
* unor = kops un 'PSCA' un ;
** up = up / unor ;
* upln = evol chpo up scal line1 ;
* dess (upln ) mima titx txl tity tyup gril ;
finsi ;
si (exist rv 'INCO' 'NT') ;
nt = elno tdom (rv.inco.'NT') ;
titr (chaine rv.titolo ' nt') ;
nte1 = evol chpo nt scal line1 ;
nte2 = evol chpo nt scal line2 ;
nte = nte1 et nte2 ;
dess (nte ) titx txo tity tynt gril lege tab3 ;
nr = nt / nuf ;
trac nt mesh cmesh ;
titr (chaine rv.titolo ' nr') ;
trac nr mesh cmesh ;
finsi ;
* evoluzioni nel tempo
si (exist rv 'HIST' '1UN') ;
dess rv.hist.'1UN' ;
dess rv.hist.'2UN' ;
finsi ;
si (exist rv 'HIST' '1VN') ;
dess rv.hist.'1VN' ;
dess rv.hist.'2VN' ;
finsi ;
si (exist rv 'HIST' 'KN') ;
dess rv.hist.'KN' ;
dess rv.hist.'EN' ;
finsi ;
si (exist rv 'HIST' 'A' ) ;
dess rv.hist.'A' ;
finsi ;
si (exist rv 'HIST' 'TF') ;
dess rv.hist.'TF' ;
finsi ;
si (exist rv 'HIST' 'TP') ;
dess rv.hist.'TP' ;
finsi ;
finproc ;
test rv complet ;
si ( 'EGA' graph 'O' );
visu rv ;
finsi ;
fin ;
© Cast3M 2003 - All rights reserved.
Disclaimer