* fichier : rayo-2D-4-bis.dgibi ************************************************************* * * * Calcul d'une plaque infinie (de largeur 2L) soumise à * * de la convection et du rayonnement. * * * * Modélisation plane. * * * * Auteurs : Michel Bulik & Nadia Coulon * * * * Date : Octobre 1995 * * * * Références : * * ------------ * * [1] Klaus-Jürgen Bathe & Mohammad R. Khoshgoftaar, * * Finite element formulation and solution of * * non-linear heat transfer, Nuclear Engineering and * * Design, v. 51 (1979), pp. 389-401 * * * * [2] J. Joly, Cas tests non linéaires de validation pour * * DELFINE, Note technique EMT.SMTS.TTMF 84/29 * * * * [3] Haji-Sheikh, A. & Sparrow, A. M., The solution of * * Heat Conduction Problems by Probability Methods, * * Transactions of ASME, Journal of Heat Transfer, 89 * * (1967), pp. 121-130 * * * ************************************************************* * * La solution de référence a été obtenue par une méthode de * Monte-Carlo décrite dans [3] et que réalise le programme * suivant en C. Après sa compilation il convient de demander * l'éditions de liens avec les subroutines quarti, cubic et * quadra de Castem2000. * * #include* * #define NB_TOTAL 100000 * #define NB_INTERVAL 20 * #define NB_PAS 8000 * * double actuel[NB_INTERVAL+1], suivant[NB_INTERVAL+1], phi[NB_PAS+1] ; * int les_nk[NB_PAS+1] ; * * #define Bi 0.2 * #define Gamma 1.0 * * double petit_l, delta_t ; * double A_0, A_1, A_2, A_3, A_4 ; * * char nomfichier[128] ; * FILE *f_out ; * * void main() * { * int i, j, k, approx ; * * /*** Initialisation ***/ * for(i=0;i<=NB_INTERVAL;i++) actuel[i] = 0 ; * actuel[0] = NB_TOTAL ; * les_nk[0] = actuel[0] ; * * /*** Boucle sur les pas - calcul des n_k ***/ * for(i=1;i<=NB_PAS;i++) { * * for(j=0;j<=NB_INTERVAL;j++) suivant[j] = 0 ; * * for(j=0;j<=NB_INTERVAL;j++) { * if(j==0) suivant[j+1] += actuel[j * else if(j==NB_INTERVAL) suivant[j-1] += actuel[j * else { * suivant[j-1] += actuel[j * suivant[j+1] += actuel[j * } * } * * for(j=0;j<=NB_INTERVAL;j++) actuel[j] = suivant[j] ; * * les_nk[i] = actuel[0] ; * } * * petit_l = 1.0 / NB_INTERVAL ; * * A_1 = 2.0 + (Bi * petit_l) ; * A_2 = 0.0 ; * A_3 = 0.0 ; * A_4 = Gamma * petit_l ; * * phi[0] = 1.0 ; * * /*** Boucle sur les pas - calcul des phi ***/ * for(i=1;i<=NB_PAS;i++) { * * double phiact ; * double q1, q2, q3, q4 ; * int nroot ; * * A_0 = -1.0 - phi[i-1] ; * if(i > 1) { * double lephi, tototo ; * * for(j=1;j<=i-1;j++) { * lephi = phi[i - j] ; * tototo = les_nk[j] * lephi * petit_l * * (Bi + (Gamma * lephi * lephi * l * tototo /= NB_TOTAL ; * A_0 += tototo ; * } * } * * if( Gamma > 0.0) { * quarti(&A_4,&A_3,&A_2,&A_1,&A_0,&q1,&q2,&q3,&q4, * phiact = q1 ; * } * else { * phiact = -A_0 / A_1 ; * } * * phi[i] = phiact ; * } * * sprintf(nomfichier,"tw.Bi=%3.1f.Gamma=%3.1f.xy\0",Bi,Gamma) ; * f_out = fopen(nomfichier,"w") ; * if(f_out == NULL) exit(5) ; * delta_t = petit_l * petit_l / 2 ; * for(i=1;i<=NB_PAS;i++) fprintf(f_out,"%15g %15g\n",i*delta_t,phi * fclose(f_out); * } * ************************************************************************* 'DEBP' pipo PRECED*'TABLE' tt*'FLOTTANT' IAPPEL*'ENTIER'; etab=preced.'WTABLE'; CTE_SB = ETAB.'CTE_STEFAN_BOLTZMANN' ; ICALC=VRAI; 'SI' ( 'EXIS' etab 'DEJA_CALC') ; 'SI' ( 'EGA' etab.'DEJA_CALC' TT); ICALC=FAUX; 'FINSI'; 'FINSI'; 'SI' ICALC; etab.'DEJA_CALC'=TT; 'FINSI'; TREF=ETAB.'TREF'; TEMM=TT; CTE_SB = ETAB.'CTE_STEFAN_BOLTZMANN' ; U0=ETAB.'THER_COURANT' + TREF; MAT_CHPO=(U0 exco 'T' 'Q' )*0.; MAT_CHPO= CHANG 'ATTRIBUT' MAT_CHPO 'NATURE' 'DISCRET'; mo_ray = preced.'MON_MOD_RAY' ; ma_ray = preced.'MON_MAT_RAY'; CH_EMI = ma_ray; CH_EMI2 = CHAN NOEUD moray CH_EMI ; MAIL_RAD = 'EXTR' CH_EMI 'MAIL' ; CH_RAD = REDU U0 MAIL_RAD ; CH_TER = 'TIRE' (ETAB.'CHARGEMENT') 'TERA' TEMM ; CH_TER = CH_TER + TREF ; *--------- Calcul des différents CHPOINT sur le maillage ---------- *------------- Cas du rayonnement à l'infini ---------------------- mcv = ('EXTR' mo_ray 'MAIL') 'MODE' thermique CONVECTION ; ma_inf = mate mo_ray 'EMIS' 1.; ; CH_TE1 = 'REDU' CH_TER MAIL_RAD; MAT_RAD= HRAYO MCV mo_ray ma_ray CH_RAD mo_ray ma_inf CH_TE1 CTE_SB ; * CHAL_TES = 'CONVECTION' MCV MAT_RAD CH_TE1 ; MAT_CHPO = MAT_CHPO + CHAL_TES ; * RIG_RAD = 'CONDUCTIVITE' MCV MAT_RAD ; TAA=TABLE; TAA.'ADDI_SECOND'= MAT_CHPO; TAA.'ADDI_MATRICE'=RIG_rad ; TAA.'RAYO_VALEUR' = RAYO_VAL; 'FINPROC' TAA; *** Options ... opti dime 2 elem seg2 echo 1 ; *** Paramètres ... L = 1. ; ep = L / 10 ; graph = faux ; Tempini = 273. ; Tempext = 0. * Tempini ; cte_sb = 5.673e-8 ; Gamma = 1. ; si ( neg Gamma 0. ) ; lambda = (cte_sb * Tempini * Tempini * Tempini * L) / Gamma ; mess 'Lambda = ' lambda ; valemis = 1. ; sinon ; lambda = 1. ; mess 'Lambda choisie arbitrairement = ' lambda ; valemis = 0. ; finsi ; Bi = 0.2 ; *** Points ... dens 0.1 ; p1 = 0 0 ; p2 = 0 ep ; vechoriz = L 0 ; *** Lignes ... li1 = p1 d 1 p2 ; *** Surface ... opti elem qua4 ; su1 = li1 tran vechoriz dini 0.1 dfin 0.01 ; li2 = cote 3 su1 ; p3 = li2 poin proc vechoriz ; li3 = cote 2 su1 ; si(graph) ; titr 'Le maillage de la plaque' ; trac su1 ; finsi ; *** Modèles ... mocnd = mode su1 thermique ; macnd = mate mocnd 'K' lambda 'RHO' 1. 'C' lambda ; mocnv = mode li2 thermique convection CONS 'CCONV'; macnv = mate mocnv 'H' (Bi*lambda/L) ; moray = mode li2 thermique rayonnement INFINI CONS 'CRAYO' ; maray = mate moray 'EMIS' valemis ; *** Préparation de la table pour PASAPAS ... tabnl = table ; tabnl . MODELE = mocnd et mocnv et moray; tabnl . CARACTERISTIQUES = macnd et macnv et maray ; lr1 = prog 0 100 ; lr2 = prog 1 1 ; ev1 = evol manu 't' lr1 'g(t)' lr2 ; chtext = manu chpo li2 1 'T' Tempext ; tabnl . CHARGEMENT =( char 'TECO' ev1 chtext) et ( char 'TERA' ev1 chtext) ; tabnl . TEMPS_CALCULES = prog 0. pas 2.e-4 0.001 pas 5.e-4 0.01 pas 0.005 0.1 pas 0.05 1. pas 0.2 10. ; tabnl . TEMPERATURES = table ; tabnl . TEMPERATURES . 0 = manu chpo su1 1 'T' Tempini ; * tabnl .'MON_MOD_RAY'= moray ; * tabnl .'MON_MAT_RAY' = maray; * tabnl.'PROCEDURE_CHAR_THER'=vrai; tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ; tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ; tabnl . 'RELAXATION_THETA' = 1. ; *** Appel à PASAPAS ... pasapas tabnl ; *** Petit post-traitement ... listt = prog ; listtp1 = prog ; listtp1r = prog ; listtp3 = prog ; listtp3r = prog ; nbpas = dime (tabnl . TEMPS) ; repeter surpas nbpas ; lindice = &surpas - 1 ; listt = listt et (prog (tabnl . TEMPS . lindice)) ; valtp1 = extr (tabnl . TEMPERATURES . lindice) T p1 ; listtp1 = listtp1 et (prog valtp1) ; valtp1r = (valtp1 - Tempext) / (Tempini - Tempext) ; listtp1r = listtp1r et (prog valtp1r) ; valtp3 = extr (tabnl . TEMPERATURES . lindice) T p3 ; listtp3 = listtp3 et (prog valtp3) ; valtp3r = (valtp3 - Tempext) / (Tempini - Tempext) ; listtp3r = listtp3r et (prog valtp3r) ; * ... Vérification s'il n'y a pas d'oscillations spatiales ... si(faux) ; chtit = chai 'Profil de temperature au temps ' (tabnl . TEMPS . lindice) ; titr chtit ; profil = evol chpo (tabnl . TEMPERATURES . lindice) 'T' li3 ; dess profil ; finsi ; fin surpas ; titr 'Evolution de la temperature au centre de la plaque' ; evtp1 = evol manu 't' listt 'T' listtp1 ; titr 'Evolution de la temperature relative au centre de la plaque' ; evtp1r = evol manu 't' listt 'T' listtp1r ; titr 'Evolution de la temperature au bord de la plaque' ; evtp3 = evol manu 't' listt 'T' listtp3 ; titr 'Evolution de la temperature relative au bord de la plaque' ; evtp3r = evol manu 't' listt 'T' listtp3r ; si(graph) ; dess evtp1 ; dess evtp1r ; dess evtp3 ; dess evtp3r ; finsi ; *** Vérification si c'est OK ... listtref = prog 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. ; listfref = prog 0.566607 0.437195 0.348971 0.283422 0.232317 0.19139 0.158117 0.130835 0.108357 0.0897874 ; titr 'Evolution de la temperature relative au bord de la plaque' ' (M-C)' ; evref = evol manu 't' listtref 'T' listfref ; si(graph) ; dess (evref et evtp3r) ; finsi ; nbtests = dime listtref ; maxerr = 0. ; repeter surtst nbtests ; instant = extr listtref &surtst ; valref = extr listfref &surtst ; tcalc = peche tabnl TEMPERATURES instant ; valcalc = extr tcalc T p3 ; valcalcr = (valcalc - Tempext) / (Tempini - Tempext) ; errrel = valcalcr - valref ; mess 'Erreur relative à t =' instant 'est de' errrel ; errrela = abs errrel ; si(errrela > maxerr) ; maxerr = errrela ; finsi ; fin surtst ; si(maxerr > 0.01) ; erre 5 ; finsi ; *** Bye ... fin ;
© Cast3M 2003 - All rights reserved.
Disclaimer