* fichier : soudage1.dgibi
***********************************************************************
* NOM : soudage1.dgibi
* ___ pb. d'advection-diffusion avec rayonnement par une
* méthode d'éléments finis
* Equation de l'energie + modèle de plaque plane
*
* Biblio :
* Rapport DM2S SFME/LTMF/RT/02-052/A
* S. Gounand
* Deux cas-tests en vue d'une modélisation numérique d'un
* processus de soudage.
*
* LANGAGE : GIBIANE-CAST3M
* AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
* mél : gounand@semt2.smts.cea.fr
*
******************************************************************
*
* Paramètres : * interact = VRAI si on travaille en interactif
* FAUX sinon
* * graph = VRAI si on veut des graphiques
* FAUX sinon
* * typdisc = 'LINE' : éléments finis linéaires
* 'QUAF' : éléments finis quadratiques
* * grossier = VRAI : maillage grossier
* FAUX : maillage fin
*
*interact= VRAI ;
*graph = VRAI ;
interact= FAUX ;
graph = FAUX ;
*typdisc = 'LINE';
typdisc = 'QUAF';
grossier = FAUX ;
*
******************************************************************
'OPTION' 'DIME' 2 ;
'SI' ('EGA' typdisc 'LINE') ;
'OPTION' 'ELEM' QUA4 ;
'SINON' ;
'OPTION' 'ELEM' QUA8 ;
'FINSI' ;
'OPTION' 'ISOV' 'SULI' ;
nbisov = 15 ;
'SI' ('NON' interact) ;
'OPTION' 'TRAC' 'PSC' ;
'OPTION' 'ECHO' 0 ;
'FINSI' ;
*
formflot ='(1PE9.2)' ;
*
*=======================================================================
* DEBUT DE LA PROCEDURE TRACVIT
*=======================================================================
'DEBPROC' TRACVIT ;
'ARGUMENT' vit*'CHPOINT ' ;
'ARGUMENT' tit/'MOT ' ;
'ARGUMENT' echvit/'FLOTTANT' ;
ladens = (** ('/' ('MESURE' mt) ('NBEL' mt))
('/' 1.D0 ('VALEUR' 'DIME'))) ;
'SI' ('NON' ('EXISTE' echvit)) ;
echvit = 'MAXIMUM' vit 'ABS' ;
'SINON' ;
echvit = 'ABS' echvit ;
'FINSI' ;
'SI' ('<' echvit 1.D-30) ;
echvit = 1.D0 ;
'FINSI' ;
uref = '/' ('*' ladens 2.0D0) echvit ;
vecvit = 'VECTEUR' vit uref 'UX' 'UY' 'JAUN' ;
'SI' ('EXISTE' tit) ;
titvit = 'CHAINE' 'FORMAT' formflot tit ' echvit=' echvit titglob ;
'SINON' ;
titvit = 'CHAINE' 'FORMAT' formflot 'VIT echvit=' echvit titglob ;
'FINSI' ;
'TRACER' vecvit mt 'TITR' titvit ;
'FINPROC' ;
*=======================================================================
* FIN DE LA PROCEDURE TRACVIT
*=======================================================================
*=======================================================================
* DEBUT DE LA PROCEDURE TRACCHML
*=======================================================================
'DEBPROC' TRACCHML ;
'ARGUMENT' chml*'CHPOINT ' ;
'ARGUMENT' tit/'MOT ' ;
rescal = 'KCHA' $mt chml 'CHAM' ;
tituap = 'CHAINE' tit ' ' titglob ;
'TRACER' rescal modbid 'TITR' tituap ;
'FINPROC' ;
*=======================================================================
* FIN DE LA PROCEDURE TRACCHML
*=======================================================================
*=======================================================================
* DEBUT DE LA PROCEDURE CALCRAYO
*=======================================================================
'DEBPROC' CALCRAYO ;
'ARGUMENT' rx*'TABLE ' ;
rv = rx . 'EQEX' ;
rvirayo = rx . 'ARG1' ;
tn = 'NOMC' 'SCAL' (rv . 'INCO' . 'T') 'NATURE' 'DISCRET' ;
tn4 = '-' (tn '*' tn '*' tn '*' tn) (t0 '*' t0 '*' t0 '*' t0) ;
rs = '*' tn4 ('*' -2.D0 sigma) ;
srs = 'NOMC' 'SOUR' rs 'NATURE' 'DISCRET' ;
srs = 'KCHT' $mt 'SCAL' 'SOMMET' 'COMP' 'SOUR' srs ;
srs2 = 'NOEL' $mt srs ;
isr = 'MAXIMUM' ('RESULT' ('*' srs2 ('DOMA' $mt 'VOLUME'))) ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'integrale source rayonnement='
isr ' W' ) ;
*rescal = srs2 ;
*tituap = 'CHAINE' 'FORMAT' formflot 'Src. rayo. (en W.m-2) ; int='
* isr ' W' ;
*TRACCHML rescal tituap ;
*rv . 'INCO' . 'PRAYO' = '*' srs2 copuis ;
rv . 'INCO' . rvirayo = '*' srs2 copuis ;
as2 ama1 = 'KOPS' 'MATRIK' ;
'RESPRO' as2 ama1 ;
'FINPROC' ;
*=======================================================================
* FIN DE LA PROCEDURE CALCRAYO
*=======================================================================
*
*
* Donnes materielles (systeme SI)
* Geometrie
* long : longeur de la plaque
* larg : largeur de la plaque
* epai : epaisseur de la plaque
* xs, ys : coordonnées de la source gaussienne
* ras : rayon de la source gaussienne
* dg : densité grosse (pour le maillage)
* dp : densité petite (pour le maillage)
long = 35.D-2 ;
larg = 13.5D-2 ;
*epai = 2.5D-2 ;
epai = 2.5D-3 ;
xs = 7.D-2 ;
ys = 0.D-2 ;
*ras = 3.D-3 ;
ras = 1.5D-3 ;
'SI' (grossier) ;
dg = 5.D-2 ;
dp = 2.D-2 ;
'SINON' ;
dg = 2.5D-2 ;
dp = '/' ras 50.D0 ;
'FINSI' ;
* Vitesse défilement de la plaque
ux = 1.D-2 ;
*ux = 0.23333D-2 ;
* Physique
* rho : masse volumique
* cp : chaleur spécifique à pression constante
* mu : viscosite (dynamique)
* kappa : conductivité thermique
* puis : puissance de la source gaussienne
* sigma : constante de Stefan
rho = 7800.D0 ;
cp = 803.D0 ;
kappa = 25.D0 ;
puis = 900.D0 ;
sigma = 5.6710D-8 ;
* Conditions aux limites
t0 = 300.D0 ;
tg = t0 ;
*
* titre global pour les dessins
*
titglob = 'CHAINE' 'FORMAT' formflot
' puis=' puis ' W' ;
'MESSAGE' '**************' ;
'MESSAGE' 'TP1 :' ;
'MESSAGE' titglob ;
'MESSAGE' '**************' ;
*
* Géométrie discrétisée (maillage régulier)
*
pA = 0.D0 0.D0 ;
pB = xs 0.D0 ;
pC = long 0.D0 ;
pD = long larg ;
pE = xs larg ;
pF = 0.D0 larg ;
bas1 = 'DROIT' pA pB 'DINI' dg 'DFIN' dp ;
bas2 = 'DROIT' pB pC 'DINI' dp 'DFIN' dg ;
dro = 'DROIT' pC pD 'DINI' dp 'DFIN' dg ;
hau2 = 'DROIT' pD pE 'DINI' dg 'DFIN' dp ;
hau1 = 'DROIT' pE pF 'DINI' dp 'DFIN' dg ;
gau = 'DROIT' pF pA 'DINI' dg 'DFIN' dp ;
bas = 'ET' bas1 bas2 ;
hau = 'ET' hau2 hau1 ;
mt = 'DALLER' bas dro hau gau ;
'SI' graph ;
titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
' NBELEM=' ('NBEL' mt) titglob ;
'TRACER' mt 'TITRE' titgeo ;
'FINSI' ;
*
* Modèle
*
_mt = 'CHANGER' mt 'QUAF' ;
_bas = 'CHANGER' bas 'QUAF' ;
_dro = 'CHANGER' dro 'QUAF' ;
_hau = 'CHANGER' hau 'QUAF' ;
_gau = 'CHANGER' gau 'QUAF' ;
'ELIMINATION' (_mt 'ET' _bas 'ET' _dro 'ET' _hau 'ET' _gau)
('/' dp 100.D0) ;
modbid = 'MODELISER' mt 'THERMIQUE' ;
$mt = 'MODELISER' _mt 'NAVIER_STOKES' typdisc ;
$bas = 'MODELISER' _bas 'NAVIER_STOKES' typdisc ;
$dro = 'MODELISER' _dro 'NAVIER_STOKES' typdisc ;
$hau = 'MODELISER' _hau 'NAVIER_STOKES' typdisc ;
$gau = 'MODELISER' _gau 'NAVIER_STOKES' typdisc ;
mtc = 'DOMA' $mt 'CENTRE' ;
*
* Construction de la source gaussienne
*
xmtc ymtc = 'COORDONNEE' mtc ;
xxs = '-' xmtc xs ;
yys = '-' ymtc ys ;
fes = '*' ('+' ('*' xxs xxs) ('*' yys yys))
('/' -3.D0 ('*' ras ras)) ;
gs = '*' ('EXP' fes)
('/' 3.D0 (PI '*' ras '*' ras)) ;
sgs = 'NOMC' 'SOUR' ('*' gs puis) 'NATURE' 'DISCRET' ;
sgs = 'KCHT' $mt 'SCAL' 'CENTRE' 'COMP' 'SOUR' sgs ;
*
* Rescaling de la source (attention, on ne modélise qu'une moitié
* de la source : symetrie par rapport à Ox)
*
isgavc = 'MAXIMUM' ('RESULT' ('*' sgs ('DOMA' $mt 'VOLUME'))) ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. source gauss avant correction='
isgavc ' W' ) ;
fcor = '/' ('/' puis 2.D0) isgavc ;
sgs = '*' sgs fcor ;
isgavc = 'MAXIMUM' ('RESULT' ('*' sgs ('DOMA' $mt 'VOLUME'))) ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. source gauss apres correction='
isgavc ' W' ) ;
'SI' graph ;
rescal = sgs ;
tituap = 'CHAINE' 'Src. gauss. (en W.m-2)' ;
TRACCHML rescal tituap ;
'FINSI' ;
stot = sgs ;
*
* conditions initiales
*
* temperature
tini = 'KCHT' $mt 'SCAL' 'SOMMET' 'COMP' 'T' t0 ;
* vitesse
uini = 'KCHT' $mt 'VECT' 'SOMMET' (ux 0.D0) ;
*
* conditions aux limites
*
* temperature
mailt = 'DOMA' $gau 'SOMMET' ;
climt = 'MANUEL' 'CHPO' mailt 1 'T' tg 'NATURE' 'DISCRET' ;
*
climtot = climt ;
*
* Mise en place du calcul numérique
*
codiff = '/' kappa (rho '*' cp) ;
copuis = '/' 1.D0 (epai '*' rho '*' cp) ;
* options + équation de l'énergie
rv = 'EQEX' 'NITER' 7 'OMEGA' 1.D0
'OPTI' 'EF' 'IMPL' 'CENTREE'
'ZONE' $mt 'OPER' 'KONV' 1.D0 uini codiff 'INCO' 'T'
'OPTI' 'EF' 'IMPL' 'CENTREE'
'ZONE' $mt 'OPER' 'LAPN' codiff 'INCO' 'T' 'T'
'OPTI' 'EF' 'IMPL'
'ZONE' $mt 'OPER' 'FIMP' ('*' stot copuis) 'INCO' 'T'
;
* rayonement
rv = 'EQEX' rv
'OPTI' 'EF' 'IMPL'
'ZONE' $mt 'OPER' 'CALCRAYO' 'PRAYO'
'OPTI' 'EF' 'IMPL'
'ZONE' $mt 'OPER' 'FIMP' 'PRAYO' 'INCO' 'T'
;
rv . 'CLIM' = climtot ;
rv . 'INCO' = 'TABLE' 'INCO' ;
rv . 'INCO' . 'T' = 'KCHT' $mt 'SCAL' 'SOMMET' 'COMP' 'T' tini ;
*
* Calcul proprement dit
*
'TEMPS' 'ZERO' ;
EXEC rv ;
tcpu = 'TEMPS' 'NOEC' ;
tcpus = '/' ('FLOTTANT' tcpu) 100.D0 ;
'MESSAGE' ('CHAINE' 'tcpus=' tcpus) ;
*'TEMPS' 'IMPR' ;
*
* Post-traitement des variables
*
vart = rv . 'INCO' . 'T' ;
covarfd = -1.D0 '*' epai '*' kappa ;
varfd = '*' ('KOPS' vart 'GRADS' $mt) covarfd ;
covarfc = epai '*' rho '*' cp ;
varfc = '*' uini ('NOMC' 'SCAL' vart 'NATURE' 'DISCRET') ;
varfc = '*' varfc covarfc ;
'SI' graph ;
rescal = vart ;
tituap = 'CHAINE' 'Temperature (en K)' titglob ;
'TRACER' rescal mt mt 'TITR' tituap nbisov ;
rescal = varfd ;
tituap = 'CHAINE' 'Flux diffusif (en W.m-1)' ;
TRACVIT rescal tituap ;
rescal = varfc ;
tituap = 'CHAINE' 'Flux convectif (en W.m-1)' ;
TRACVIT rescal tituap ;
'FINSI' ;
*
* Calcul des intégrales des flux diffusifs
*
ifdbas = 'DBIT' $bas varfd ;
ifddro = 'DBIT' $dro varfd ;
ifdhau = 'DBIT' $hau varfd ;
ifdgau = 'DBIT' $gau varfd ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux diff. bas =' ifdbas ' W') ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux diff. dro =' ifddro ' W') ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux diff. hau =' ifdhau ' W') ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux diff. gau =' ifdgau ' W') ;
*
* Calcul des intégrales des flux convectifs
*
ifcbas = 'DBIT' $bas varfc ;
ifcdro = 'DBIT' $dro varfc ;
ifchau = 'DBIT' $hau varfc ;
ifcgau = 'DBIT' $gau varfc ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux conv. bas =' ifcbas ' W') ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux conv. dro =' ifcdro ' W') ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux conv. hau =' ifchau ' W') ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. flux conv. gau =' ifcgau ' W') ;
*
* Calcul de l'énergie apportée par la source
*
isg = 'MAXIMUM' ('RESULT' ('*' stot ('DOMA' $mt 'VOLUME'))) ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. source gauss =' isg ' W') ;
*
* Calcul de l'énergie rayonnée
*
srtot = '/' (rv . 'INCO' . 'PRAYO') copuis ;
isr = 'MAXIMUM' ('RESULT' ('*' srtot ('DOMA' $mt 'VOLUME'))) ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'int. source rayo. =' isr ' W') ;
*
* Bilan d'energie (gain/perte en %)
*
ifd = ifdbas '+' ifddro '+' ifdhau '+' ifdgau ;
ifc = ifcbas '+' ifcdro '+' ifchau '+' ifcgau ;
be = ifd '+' ifc '+' isg '+' isr ;
berp = '*' ('/' be ('/' puis 2.D0)) 100.D0 ;
'MESSAGE' ('CHAINE' 'FORMAT' formflot
'Verif bilan energie =' berp ' %') ;
'SI' interact ;
'OPTION' 'DONN' 5 ;
'SINON' ;
'SI' ('>' ('ABS' berp) 1.D-2) ;
'ERREUR' 5 ;
'SINON' ;
'ERREUR' 0 ;
'FINSI' ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer