* fichier : SOURCE3.dgibi
* section : thermique
*----------------------------------------------------------------------*
* SOURCE3.DGIBI *
*----------------------------------------------------------------------*
*
* Objet :
* -------
*
* Verfication / validation d'un modele de source de chaleur.
* Cas d'une SOURCE GAUSSIENNE ELLIPTIQUE.
*
* Description :
* -------------
* Comparaison des flux nodaux equivalents obtenus avec le modele a
* ceux obtenus en construisant le champ "a la main", puis en l'integrant
* avec l'operateur SOURCE.
*
* Type de calcul : Aucun
* Mode de calcul : 2D PLAN, AXIS et 3D
* Type d'element : TRI3 TRI6 QUA4 QUA8 TET4 TE10 PYR5 PY13 PRI6 PR15
* CUB8 CU20
* Objectifs : Ecart relatif entre flux integres < 1.e-12
* Ecart relatif entre la chaleur totale fournie (QTOT)
* et la resultante sur le champ integre < 1.e-12
*
*----------------------------------------------------------------------*
*
* IG1 vrai : traces actives
IG1 = faux ;
*
*------------------------ 2D ELEMENTS LINEAIRES -----------------------*
*
* Parametres du maillage :
lo1 = 50.e-3 ;
ep1 = 20.e-3 ;
* Parametres de la source :
QT1 = 1.e3 ;
OR1 = 30.e-3 ep1 ;
DR1 = 0. 2. ;
RG1 = 5.e-3 ;
ZG1 = 9.e-3 ;
* Maillage :
l1
= (0 0) droi 10 (0 ep1
) ;s1
= l1
tran 50 (lo1
0) ;s2
= l2
tran 50 (lo1 ep1
) ;s0 = s1 et s2 ;
si IG1 ;
mot1
= chai mot1 ' '
(ltyp1
extr &bm1
) ; fin bm1 ;
mot1
= chai ' Maillage test en 2D
(' mot1 '
)'
;fins;
* MODE PLAN :
* -----------
*
*
moq1
= mode s0 thermique source gaussienne elliptique
;maq1
= mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 'DIRE' DR1 'ZGAU' ZG1
;
mo1
= mode s0 thermique
;
x1 y1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) ;x1 y1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) ;
z1
= (x1
* (DR1
coor 1)) + (y1
* (DR1
coor 2)) / (norm DR1
) ;r1 = r1 - (z1 * z1) ;
xqx1 = 4. / PI / RG1 / ZG1 ;
sqref1 = exp (-2. * ((r1 / RG1 / RG1) + (z1 * z1 / ZG1 / ZG1))) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
si IG1 ;
trac sq1 S0
titr '
Flux nodaux du champ de source Gaussienne du MODELE '
; trac sqref1 s0
titr '
Flux nodaux du champ de source Gaussienne REFERENCE'
; trac (sq1
- sqref1
) s0
titr ' Difference des
2 champs'
; fins ;
err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2 = (abs (qtmoq1 - qt1)) / (abs qt1) ;
si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
fins ;
* MODE AXIS :
* -----------
*
*
moq1
= mode s0 thermique source gaussienne elliptique
;maq1
= mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 'DIRE' DR1 'ZGAU' ZG1
;
mo1
= mode s0 thermique
;
x1 y1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) ;x1 y1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) ;
z1
= (x1
* (DR1
coor 1)) + (y1
* (DR1
coor 2)) / (norm DR1
) ;r1 = r1 - (z1 * z1) ;
rc2 = 2. ** 0.5 ;
rcpi1 = pi ** 0.5 ;
xik1 = 0.25 * RG1 * RG1 * (exp (-2. * x0 * x0 / RG1 / RG1)) ;
xik2 = 0.5 * RG1 * x0 * rcpi1 / rc2 ;
xik3
= RG1
* x0
* (erf (x0
* rc2
/ RG1
)) * rcpi1
/ (rc2
* rc2
* rc2
) ;xqx1 = PI * ZG1 * rcpi1 / rc2 ;
xqx1 = xqx1 * (xik1 + xik2 + xik3) ;
sqref1 = exp (-2. * ((r1 / RG1 / RG1) + (z1 * z1 / ZG1 / ZG1))) * QT1 / xqx1 ;
sqref1
= sour mo1 sqref1
;
si IG1 ;
trac sq1 S0
titr '
Flux nodaux du champ de source Gaussienne du MODELE '
; trac sqref1 s0
titr '
Flux nodaux du champ de source Gaussienne REFERENCE'
; trac (sq1
- sqref1
) s0
titr ' Difference des
2 champs'
; fins ;
err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2 = (abs (qtmoq1 - qt1)) / (abs qt1) ;
si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
fins ;
*---------------------- 2D ELEMENTS QUADRATIQUES ----------------------*
* Passage en EF quadratiques :
si IG1 ;
mot1
= chai mot1 ' '
(ltyp1
extr &bm1
) ; fin bm1 ;
mot1
= chai ' Maillage test en 2D
(' mot1 '
)'
;fins;
* MODE PLAN :
* -----------
*
*
moq1
= mode s0 thermique source gaussienne elliptique
;maq1
= mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 'DIRE' DR1 'ZGAU' ZG1
;
mo1
= mode s0 thermique
;
x1 y1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) ;x1 y1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) ;
z1
= (x1
* (DR1
coor 1)) + (y1
* (DR1
coor 2)) / (norm DR1
) ;r1 = r1 - (z1 * z1) ;
xqx1 = 4. / PI / RG1 / ZG1 ;
sqref1 = exp (-2. * ((r1 / RG1 / RG1) + (z1 * z1 / ZG1 / ZG1))) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
si IG1 ;
trac sq1 S0
titr '
Flux nodaux du champ de source Gaussienne du MODELE '
; trac sqref1 s0
titr '
Flux nodaux du champ de source Gaussienne REFERENCE'
; trac (sq1
- sqref1
) s0
titr ' Difference des
2 champs'
; fins ;
err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2 = (abs (qtmoq1 - qt1)) / (abs qt1) ;
si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
fins ;
* MODE AXIS :
* -----------
*
*
moq1
= mode s0 thermique source gaussienne elliptique
;maq1
= mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 'DIRE' DR1 'ZGAU' ZG1
;
mo1
= mode s0 thermique
;
x1 y1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) ;x1 y1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) ;
z1
= (x1
* (DR1
coor 1)) + (y1
* (DR1
coor 2)) / (norm DR1
) ;r1 = r1 - (z1 * z1) ;
rc2 = 2. ** 0.5 ;
rcpi1 = pi ** 0.5 ;
xik1 = 0.25 * RG1 * RG1 * (exp (-2. * x0 * x0 / RG1 / RG1)) ;
xik2 = 0.5 * RG1 * x0 * rcpi1 / rc2 ;
xik3
= RG1
* x0
* (erf (x0
* rc2
/ RG1
)) * rcpi1
/ (rc2
* rc2
* rc2
) ;xqx1 = PI * ZG1 * rcpi1 / rc2 ;
xqx1 = xqx1 * (xik1 + xik2 + xik3) ;
sqref1 = exp (-2. * ((r1 / RG1 / RG1) + (z1 * z1 / ZG1 / ZG1))) * QT1 / xqx1 ;
sqref1
= sour mo1 sqref1
;
si IG1 ;
trac sq1 S0
titr '
Flux nodaux du champ de source Gaussienne du MODELE '
; trac sqref1 s0
titr '
Flux nodaux du champ de source Gaussienne REFERENCE'
; trac (sq1
- sqref1
) s0
titr ' Difference des
2 champs'
; fins ;
err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2 = (abs (qtmoq1 - qt1)) / (abs qt1) ;
si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
fins ;
*------------------------ 3D ELEMENTS LINEAIRES -----------------------*
*
* Maillage du volume :
* Ajour de tetraedres :
v0
= v0
et (v3
chan tet4
) ;
* Ajout d'un element pyramide :
py0 = py1 et py2 et py3 ;
v0
= (v0
diff el1
) et py0
;
si IG1 ;
mot1
= chai mot1 ' '
(ltyp1
extr &bm1
) ; fin bm1 ;
mot1
= chai ' Maillage test en 3D
(' mot1 '
)'
;fins;
* Parametres de la source en 3D
OR1 = 30.e-3 (0.5 * ep1) ep1 ;
DR1 = 0. 0. 1. ;
* Modeles :
moq1
= mode V0 thermique source gaussienne elliptique
;maq1
= mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 'DIRE' DR1 'ZGAU' ZG1
;
mo1
= mode V0 thermique
;
x1 y1 z1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) (z1
- (or1
coor 3)) ;x1 y1 z1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) (chan cham z1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) + (z1 * z1) ;
z1
= (x1
* (DR1
coor 1)) + (y1
* (DR1
coor 2)) + (z1
* (DR1
coor 3)) / (norm DR1
) ;r1 = r1 - (z1 * z1) ;
xqx1 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / ZG1 ;
sqref1 = exp (-2. * ((r1 / RG1 / RG1) + (z1 * z1 / ZG1 / ZG1))) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
si IG1 ;
trac sq1 V0
titr '
Flux nodaux du champ de source Gaussienne du MODELE '
coup or1
(0 1 0) (0 0 1) ; trac sqref1 V0
titr '
Flux nodaux du champ de source Gaussienne REFERENCE'
coup or1
(0 1 0) (0 0 1) ; trac (sq1
- sqref1
) V0
titr ' Difference des
2 champs'
coup or1
(0 1 0) (0 0 1) ; fins ;
err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2 = (abs (qtmoq1 - qt1)) / (abs qt1) ;
si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
fins ;
*---------------------- 3D ELEMENTS QUADRATIQUES ----------------------*
*
si IG1 ;
mot1
= chai mot1 ' '
(ltyp1
extr &bm1
) ; fin bm1 ;
mot1
= chai ' Maillage test en 3D
(' mot1 '
)'
;fins;
* Modeles :
moq1
= mode V0 thermique source gaussienne elliptique
;maq1
= mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 'DIRE' DR1 'ZGAU' ZG1
;
mo1
= mode V0 thermique
;
x1 y1 z1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) (z1
- (or1
coor 3)) ;x1 y1 z1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) (chan cham z1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) + (z1 * z1) ;
z1
= (x1
* (DR1
coor 1)) + (y1
* (DR1
coor 2)) + (z1
* (DR1
coor 3)) / (norm DR1
) ;r1 = r1 - (z1 * z1) ;
xqx1 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / ZG1 ;
sqref1 = exp (-2. * ((r1 / RG1 / RG1) + (z1 * z1 / ZG1 / ZG1))) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
si IG1 ;
trac sq1 V0
titr '
Flux nodaux du champ de source Gaussienne du MODELE '
; trac sqref1 V0
titr '
Flux nodaux du champ de source Gaussienne REFERENCE'
; trac (sq1
- sqref1
) V0
titr ' Difference des
2 champs'
; fins ;
err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2 = (abs (qtmoq1 - qt1)) / (abs qt1) ;
si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
fins ;
*---------------/!\ SOURCES SUR PLUSIEURS SOUS-ZONES /!\---------------*
*
* Maillage :
si IG1 ;
trac (v1
et v2
) titr ' Maillage a
2 sous
-zones'
; fins ;
* DEUX SOURCES : GAUSSIENNE & UNIFORME :
* --------------------------------------
*
* Sur 2 sous-zones disctinctes :
* ------------------------------
*
* Modeles / caracteristiques / sources :
moq1
= mode v1 thermique source gaussienne
;moq2
= mode v2 thermique source
;maq1
= mate moq1 qtot qt1 orig or1 rgau rg1
;maq2
= mate moq2 qvol 1.
e9 ;moq0 = moq1 et moq2 ;
maq0 = maq1 et maq2 ;
mo1
= mode v1 thermique
;x1 y1 z1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) (z1
- (or1
coor 3)) ;x1 y1 z1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) (chan cham z1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) + (z1 * z1) ;
xqx1 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / RG1 ;
sqref1 = exp (-2. * (r1 / RG1 / RG1)) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
mo2
= mode v2 thermique
;sqref2
= sour mo2 v2 1.
e9 ;
sqref0 = sqref1 et sqref2 ;
si IG1 ;
trac sq0 v0
titr '
Flux nodaux des deux modeles de source '
; trac (sq0
- sqref0
) v0
titr ' Difference des deux champs'
; fins ;
err0
= (maxi abs (sq0
- sqref0
)) / (maxi abs sqref0
) ;err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2
= (maxi abs (sq2
- sqref2
)) / (maxi abs sqref2
) ;si (err1 > 1.e-12) ;
fins ;
* Sur 2 sous-zones "a cheval" :
* -----------------------------
*
* Modeles / caracteristiques / sources :
moq1
= mode v0 thermique source gaussienne
;moq2
= mode v2 thermique source
;maq1
= mate moq1 qtot qt1 orig or1 rgau rg1
;maq2
= mate moq2 qvol 1.
e9 ;moq0 = moq1 et moq2 ;
maq0 = maq1 et maq2 ;
mo1
= mode v0 thermique
;x1 y1 z1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) (z1
- (or1
coor 3)) ;x1 y1 z1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) (chan cham z1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) + (z1 * z1) ;
xqx1 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / RG1 ;
sqref1 = exp (-2. * (r1 / RG1 / RG1)) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
mo2
= mode v2 thermique
;sqref2
= sour mo2 v2 1.
e9 ;
sqref0 = sqref1 et sqref2 ;
si IG1 ;
trac sq0 V0
titr '
Flux nodaux des deux modeles de source '
; trac (sq0
- sqref0
) V0
titr ' Difference des deux champs'
; fins ;
err0
= (maxi abs (sq0
- sqref0
)) / (maxi abs sqref0
) ;err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2
= (maxi abs (sq2
- sqref2
)) / (maxi abs sqref2
) ;si (err1 > 1.e-12) ;
fins ;
* DEUX SOURCES GAUSSIENNES : SPHERIQUE & ELLIPTIQUE :
* ---------------------------------------------------
* Sur 2 sous-zones disctinctes :
* ------------------------------
*
* Parametres de la source en 3D
qt2 = 0.5 * qt1 ;
OR2 = 60.e-3 (0.6*ep1) ep1 ;
DR2 = 0. 0. 2. ;
* Modeles / caracteristiques / sources :
moq1
= mode v1 thermique source gaussienne
;moq2
= mode v2 thermique source gaussienne elliptique
;maq1
= mate moq1 qtot qt1 orig or1 rgau rg1
;maq2
= mate moq2 qtot qt2 orig or2 rgau rg1 dire dr2 zgau zg1
;moq0 = moq1 et moq2 ;
maq0 = maq1 et maq2 ;
mo1
= mode v1 thermique
;x1 y1 z1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) (z1
- (or1
coor 3)) ;x1 y1 z1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) (chan cham z1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) + (z1 * z1) ;
xqx1 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / RG1 ;
sqref1 = exp (-2. * (r1 / RG1 / RG1)) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
mo2
= mode v2 thermique
;x2 y2 z2
= (x2
- (or2
coor 1)) (y2
- (or2
coor 2)) (z2
- (or2
coor 3)) ;x2 y2 z2
= (chan cham x2 mo2 stresses
) (chan cham y2 mo2 stresses
) (chan cham z2 mo2 stresses
) ;r2 = (x2 * x2) + (y2 * y2) + (z2 * z2) ;
z2
= (x2
* (dr2
coor 1)) + (y2
* (dr2
coor 2)) + (z2
* (dr2
coor 3)) / (norm dr2
) ;r2 = r2 - (z2 * z2) ;
xqx2 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / ZG1 ;
sqref2 = exp (-2. * ((r2 / RG1 / RG1) + (z2 * z2 / ZG1 / ZG1))) * QT2 * xqx2 ;
sqref2
= sour mo2 sqref2
;
sqref0 = sqref1 et sqref2 ;
si IG1 ;
trac sq0 V0
titr '
Flux nodaux des deux modeles de source '
; trac (sq0
- sqref0
) V0
titr ' Difference des deux champs'
; fins ;
err0
= (maxi abs (sq0
- sqref0
)) / (maxi abs sqref0
) ;err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2
= (maxi abs (sq2
- sqref2
)) / (maxi abs sqref2
) ;si (err1 > 1.e-12) ;
fins ;
* Sur 2 sous-zones "a cheval" :
* -----------------------------
*
* Modeles / caracteristiques / sources :
moq1
= mode v0 thermique source gaussienne
;moq2
= mode v0 thermique source gaussienne elliptique
;maq1
= mate moq1 qtot qt1 orig or1 rgau rg1
;maq2
= mate moq2 qtot qt2 orig or2 rgau rg1 dire dr2 zgau zg1
;moq0 = moq1 et moq2 ;
maq0 = maq1 et maq2 ;
mo1
= mode v0 thermique
;x1 y1 z1
= (x1
- (or1
coor 1)) (y1
- (or1
coor 2)) (z1
- (or1
coor 3)) ;x1 y1 z1
= (chan cham x1 mo1 stresses
) (chan cham y1 mo1 stresses
) (chan cham z1 mo1 stresses
) ;r1 = (x1 * x1) + (y1 * y1) + (z1 * z1) ;
xqx1 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / RG1 ;
sqref1 = exp (-2. * (r1 / RG1 / RG1)) * QT1 * xqx1 ;
sqref1
= sour mo1 sqref1
;
mo2
= mode v0 thermique
;x2 y2 z2
= (x2
- (or2
coor 1)) (y2
- (or2
coor 2)) (z2
- (or2
coor 3)) ;x2 y2 z2
= (chan cham x2 mo2 stresses
) (chan cham y2 mo2 stresses
) (chan cham z2 mo2 stresses
) ;r2 = (x2 * x2) + (y2 * y2) + (z2 * z2) ;
z2
= (x2
* (dr2
coor 1)) + (y2
* (dr2
coor 2)) + (z2
* (dr2
coor 3)) / (norm dr2
) ;r2 = r2 - (z2 * z2) ;
xqx2 = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / ZG1 ;
sqref2 = exp (-2. * ((r2 / RG1 / RG1) + (z2 * z2 / ZG1 / ZG1))) * QT2 * xqx2 ;
sqref2
= sour mo2 sqref2
;
sqref0 = sqref1 et sqref2 ;
si IG1 ;
trac sq0 V0
titr '
Flux nodaux des deux modeles de source '
; trac (sq0
- sqref0
) V0
titr ' Difference des deux champs'
; fins ;
err0
= (maxi abs (sq0
- sqref0
)) / (maxi abs sqref0
) ;err1
= (maxi abs (sq1
- sqref1
)) / (maxi abs sqref1
) ;err2
= (maxi abs (sq2
- sqref2
)) / (maxi abs sqref2
) ;si (err1 > 1.e-12) ;
fins ;
*-------------------------------- FIN ---------------------------------*
fin ;