Accueil | Cast3M

* fichier :  rayo-2D-5.dgibi
************************************************************************
************************************************************************
 
*************************************************************
*                                                           *
*   Calcul d'une plaque infinie (de largeur 2L) avec tempé- *
*   rature imposée au milieu et soumise au rayonnement.     *
*   Il s'agit de trouver la température au bord en régime   *
*   permanent.                                              *
*                                                           *
*   Néanmoins nous cherchons la solution en faisant un      * 
*   calcul transitoire avec PASAPAS pour tester celle-ci    *
*   en configuration température imposée + rayonnement vers *
*   l'infini.                                               *
*                                                           *
*   Modélisation plane.                                     *
*                                                           *
*   Auteurs : Michel Bulik & Nadia Coulon                   *
*                                                           *
*   Date : Avril 1997                                       *
*                                                           *
*   on introduit une emissivite a l infini differente de 1  *    
*      et une temperature a l'infini differente de 0K       *    
*************************************************************
*
*** Options ...
 
    opti dime 2 elem seg2 echo 1 ;
 
*** Paramètres ...
 
    L = 1. ;
    ep = L / 10 ;
 
    graph = faux ;
 
   Tempini = 1000. ;
   Tempext = 2. * Tempini ;
 
    cte_sb = 5.673e-8 ;
    Gamma = 1. ;
    si ( neg Gamma 0. ) ;
       lambda = (cte_sb * Tempini * Tempini * Tempini * L) / Gamma ;
       mess 'Lambda = ' lambda ;
       valemis = 0.8 ;
        e0  = 0.5 ; 
 
    sinon ;
       lambda = 1. ;
       mess 'Lambda choisie arbitrairement = ' lambda ;
       valemis = 0. ;
       e0  = 1. ;
    finsi ;
 
*** 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 ;
 
    moray = mode li2 thermique  rayonnement  INFINI ;
    maray = mate moray 'EMIS' valemis 'E_IN' e0;
 
*** Température imposée ...
 
    blt = bloq 'T' li1 ;
    ti  = depi blt Tempini ;
    lr1 = prog   0 100 ;
    lr2 = prog   1   1 ;
    ev1 = evol manu 't' lr1 'T' lr2 ;
    chti = char 'TIMP' ev1 ti ;
 
*** Température extérieure ...
 
    chtext = manu chpo li2 1 'T' Tempext nature diffus ;
    chte = char 'TERA' ev1 chtext ;
 
*** Préparation de la table pour PASAPAS ...
 
    tabnl = table ;
 
    tabnl . MODELE = mocnd et moray ;
    tabnl . CARACTERISTIQUES = macnd et maray ;
    tabnl . BLOCAGES_THERMIQUES = blt ;
    tabnl . CHARGEMENT = chti et chte ;
 
    tabnl . TEMPS_CALCULES = prog 0. pas  0.1  1. 
                                     pas  0.2  2. ;
 
    tabnl . TEMPERATURES = table ;
    tabnl . TEMPERATURES . 0 = manu chpo su1 1 'T' Tempini ;
 
*    tabnl . RAYONNEMENT = table ;
*    tabnl . RAYONNEMENT . 1 = table ;
*    tabnl . RAYONNEMENT . 1 . TYPE = 'INFINI' ;
*    tabnl . RAYONNEMENT . 1 . MODELE = moray ;
 
    tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ;
    tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ;
 
    tabnl . 'RELAXATION_THETA' = 1. ;
 
*** Appel à PASAPAS ...
 
    pasapas tabnl ;
 
 
*** Petit post-traitement ...
 
    listt = prog ;
    listtp3 = prog ;
    listtp3r = prog ;
    nbpas = dime (tabnl . TEMPS) ;
    repeter surpas nbpas ;
       lindice = &surpas - 1 ;
       listt = listt et (prog (tabnl . TEMPS . lindice)) ;
       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 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 evtp3  ;
       dess evtp3r ;
    finsi ;
 
*** Vérification si c'est OK ...
 
************************************************************************
*                                                                      *
*   On va résoudre analytiquement le problème en écrivant l'égalité    *
*   des flux (T = température au bord de la plaque) :                  *
*                                                                      *
*            T_imp - T                                                 *
*   lambda ------------- = epsilon * sigma * (T^4 - T_ext^4)           *
*                L                                                     *
*   avec epsilon(e,e0)                                                 *
*                                                                      *
************************************************************************
 
    ee = ((1./valemis) + (1./e0) - 1.); ee = 1./ee ;
    mess ' ee: ' ee ;
 
    A_0 = -1 * ((lambda * Tempini / L) + 
                             (ee * cte_sb * (Tempext**4))) ;
    A_1 = lambda / L ;
    A_2 = 0.0 ;
    A_3 = 0.0 ;
    A_4 = ee * cte_sb ;
 
 
    r1  = racp A_0 A_1 A_2 A_3 A_4 ;
 
*    mess 'Zéros du polynôme :' r1 ' ' r2 ' ' r3 ' ' r4 ;
 
    sol_anlt = r1 extr 2 ;
 
    sol_calc = extr listtp3 (dime listtp3) ;
    mess ' anal ' sol_anlt ; mess ' calc ' sol_calc ;
    erreur = 100 * (sol_calc - sol_anlt) / sol_anlt ;
    mess 'Erreur relative = ' erreur '%' ;
 
    si(abs(erreur) > 0.01) ;
       erre 5 ;
    sinon ;
       erre 0 ;
    finsi ;
 
*** Bye ...
 
    fin ;