Accueil | Cast3M

* fichier :  rayo-axi-3.dgibi
 
*******************************************
*******************************************
 
*************************************************************
*                                                           *
*   Calcul d'un cylindre infini (de rayon L) soumis à de    *
*   la convection et du rayonnement.                        *
*                                                           *
*   Modélisation axisymétrique.                             *
*                                                           *
*   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] James SUCEC & Ashok KUMAR, Transient Cooling of a   *
*       Solid Cylinder by Combined Convection ans Radiation *
*       at its Surface, International Journal for Numerical *
*       Methods in Engineering, vol. 6 (1973), pp. 297-312  *
*                                                           *
*************************************************************
*
* La solution de référence a été obtenue par une méthode de
* différences finies 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 la subroutine gauscl de Castem2000.
*
*       
*       #include<stdio.h>
*       #include<math.h>
*       
*       #define N       26      /* Nombre de points */
*       
*       double  A[N*N], phi_n[N], phi_n1[N] ;
*       
*       void main()
*       {
*               double  L = 1. ;
*               double  eta = 0.55 ;
*               double  Bi = 1.0 ;
*               double  Gamma = 1.0 ;
*               double  alpha = 1.0 ;
*       
*               double  delta_r, delta_R, delta_t, M ;
*               int     i, j, k, dimension, iprint ;
*       
*               delta_t = 0.001 ;
*               delta_r = L / (N - 1) ;
*               delta_R = delta_r / L ;
*               M = alpha * delta_t / (delta_r*delta_r) ;
*               dimension = N ;
*               iprint = 0 ;
*       
*               /* État initial */
*               for(i=0;i<N;i++) phi_n1[i]=1 ;
*               /* La matrice */
*               /* La première ligne */
*               A[0] = 1.0 + 4*M ;
*               A[N] = -4*M ;
*               /* Les lignes 2,...,N-1 */
*               for(j=2;j<=N-1;j++) {
*                       k = (j-1)*N + j - 1 ;
*                       A[k-N] = -1*M*(1.0-(1.0/(2*(j-1)))) ;
*                       A[k  ] = 1.0 + 2*M ;
*                       A[k+N] = -1*M*(1.0+(1.0/(2*(j-1)))) ;
*               }
*               /* La dernière ligne */
*               A[N*N-N-1] = -2*M*(1-delta_R/2) ;
*               A[N*N  -1] = 1.0 + 2*M*(1-delta_R/2) + 2*delta_R*M*Bi ;
*       
*               /* Boucle sur les pas de temps */
*               for(j=0;j<10000;j++) {
*       
*                       /* On met le résultat dans le second membre */
*                       for(i=0;i<N;i++) phi_n[i] = phi_n1[i] ;
*                       /* puis on modifie le dernier terme */
*                       phi_n[N-1] -= (2*M*delta_R*Gamma) * (pow(phi_n[N-1],4)*pow(1-eta,3) + 4*pow(phi_n[N-1],3)*pow(1-eta,2)*eta
*                                                               + 6*pow(phi_n[N-1],2)*(1-eta)*eta*eta + 4*phi_n[N-1]*pow(eta,3)  ) ;
*       
*                       gauscl(A,phi_n1,phi_n,&dimension,&dimension,&iprint) ;
*       
*                       /* Impression des résultats */
*                       printf("%13g %13g %13g\n",(j+1)*delta_t,phi_n1[0],phi_n1[N-1]) ;
*               }
*       }
*       
**************************************************************
*
*** Options ...
 
    opti dime 2 mode axis elem seg2 echo 1 ;
 
*** Paramètres ...
 
    L = 1. ;
    ep = L / 10 ;
    eta = 0.55 ;
 
    graph = faux ;
 
    Tempini = 273. ;
    Tempext = eta * Tempini ;
    cte_sb = 5.673e-8 ;
    Gamma = 1.0 ;
    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 = 1.0 ;
 
*** 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 du cylindre' ;
       trac su1 ;  
    finsi ;
 
*** Modèles ...
 
    mocnd = mode su1 thermique ;
    macnd = mate mocnd 'K' lambda 'RHO' 1. 'C' lambda ;
 
    mocnv = mode li2 thermique convection ;
    macnv = mate mocnv 'H' (Bi*lambda/L) ;
 
    moray = mode li2  thermique rayonnement  infini ;
    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 'Text' 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 . 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 ;
    listtp1 = prog ;
    listtp1r = prog ;
    listtp3 = prog ;
    listtp3r = prog ;
    nbpas = dime (tabnl . TEMPS) ;
    repeter surpas (nbpas-1) ;
       lindice = &surpas ;
       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 sur l axe du cylindre' ;
    evtp1 = evol manu 't' listt 'T' listtp1 ;
    titr 'Evolution de la temperature relative sur l axe du cylindre' ;
    evtp1r = evol manu 't' listt 'T' listtp1r ;
 
    titr 'Evolution de la temperature au bord du cylindre' ;
    evtp3 = evol manu 't' listt 'T' listtp3 ;
    titr 'Evolution de la temperature relative au bord du cylindre' ;
    evtp3r = evol manu 't' listt 'T' listtp3r ;
 
    si(graph) ;
       dess evtp1  LOGX ;
       dess evtp1r LOGX ;
       dess evtp3  LOGX ;
       dess evtp3r LOGX ;
    finsi ;
 
*** Vérification si c'est OK ...
 
    listtref = prog 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ;
    listbord = prog 0.465573 0.352181 0.276436 0.218669 0.173342
                    0.137545 0.109216 0.0867689 0.0689675 0.0548392 ;
    listcent = prog 0.95046  0.769248 0.598471 0.46495  0.36275 
                    0.284227 0.223483 0.176205 0.139228 0.110195 ;
 
    titr 'Evolution de la T relative au bord du cylindre (D-F)' ;
    evbord = (evol manu 't' listtref 'T' listbord) coul ROUG ;
 
    titr 'Evolution de la T relative au centre du cylindre (D-F)' ;
    evcent = (evol manu 't' listtref 'T' listcent) coul JAUN ;
 
    si(graph) ;
       dess (evtp3r et evbord) LOGX ;
       dess (evtp1r et evcent) LOGX ;
    finsi ;
 
    nbtests = dime listtref ;
    maxerr = 0. ;
    repeter surtst nbtests ;
       instant = extr listtref &surtst ;
       tcalc = peche tabnl TEMPERATURES instant ;
 
       valref  = extr listbord &surtst ;
       valcalc = extr tcalc T p3 ;
       valcalcr = (valcalc - Tempext) / (Tempini - Tempext) ;
       errrel = valcalcr - valref ;
       mess 'Erreur au bord   à t =' instant 'est de' errrel ;
       errrela = abs errrel ;
       si(errrela > maxerr) ;
          maxerr = errrela ;
       finsi ;
 
       valref  = extr listcent &surtst ;
       valcalc = extr tcalc T p1 ;
       valcalcr = (valcalc - Tempext) / (Tempini - Tempext) ;
       errrel = valcalcr - valref ;
       mess 'Erreur au centre à t =' instant 'est de' errrel ;
       errrela = abs errrel ;
       si(errrela > maxerr) ;
          maxerr = errrela ;
       finsi ;
    fin surtst ;
    si(maxerr > 0.02) ;
       erre 5 ;
    finsi ;
 
*** Bye ...
 
   fin ;