Procédures | Cast3M

* PRODT     PROCEDUR  JC220346  10/11/18    21:15:01     6795           
************************************************************************
* NOM         : PRODT
* DESCRIPTION : Renvoie une mesure du tenseur gradient des vitesses.
*               Cette grandeur intervient notamment dans l'expression
*               du terme de production turbulente.
************************************************************************
* HISTORIQUE :  14/03/2000 : MAGN : création de la procédure
* HISTORIQUE :  28/10/2010 : JCARDO : généralisation de PRODT
*               1) ajout de 3 autres mesures pour le tenseur gradient,
*                  en plus du taux de déformation: - taux de rotation
*                                                  - tenseur complet
*                                                  - expression mixte
*               2) prise en compte du cas où NCO est différent de NDIM
*               3) gestion de composantes de noms quelconques dans UN
*               4) modification de la syntaxe et de la notice
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS DE COMPLÉTER LES COMMENTAIRES
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*   Syntaxe :
*   _________
*
*          CHPO1 = PRODT $MD (| 'TODEF' |) UN (GB TN) ;
*                             | 'TOROT' |
*                             | 'COMPL' |
*                             | 'MIXTE' |
*
*
************************************************************************
 
 
DEBP PRODT ;
 
 
* Dimension de l'espace physique
NDIM = VALE 'DIME' ;
 
SI (EGA NDIM 1) ;
    MESS '******************************************************' ;
    MESS '            /!\    ERREUR dans PRODT    /!\           ' ;
    MESS ' ' ;
    MESS ' Cette procédure est inutile pour les problèmes 1D.' ;
    MESS '******************************************************' ;
    QUIT PRODT ;
FINS ;
 
 
 
* +====================================================================+
* |                                                                    |
* |             L E C T U R E   D E S   A R G U M E N T S              |
* |                                                                    |
* +====================================================================+
 
* +========+
* | MODÈLE |
* +========+
 
ARGU $MD*'MMODEL' ;
 
 
* +=========+
* | MOT-CLÉ |
* +=========+
 
ARGU MCLE/'MOT' ;
SI (NON (EXIS MCLE)) ;
    MCLE = 'TODEF' ;
FINS ;
 
B_TODEF = (EGA MCLE 'TODEF') ;
B_TOROT = (EGA MCLE 'TOROT') ;
B_COMPL = (EGA MCLE 'COMPL') ;
B_MIXTE = (EGA MCLE 'MIXTE') ;
 
 
* +==================+
* | CHAMP DE VITESSE |
* +==================+
 
ARGU UN*'CHPOINT' ;
LCO = EXTR UN 'COMP' ;
NCO = DIME LCO ;
 
 
* +==============================================+
* | CHAMP DE TEMPÉRATURE ET VECTEUR FLOTTABILITÉ |
* +==============================================+
 
ARGU TN/'CHPOINT' ;
 
B_TN = FAUX ;
SI (EXIS TN) ;
 
    ARGU GB/'POINT' ;
    SI (NON (EXIS GB)) ;
        MESS '******************************************************' ;
        MESS '/!\ ERREUR dans PRODT :' ;
        MESS '    Champ de température fourni => on attendait GB' ;
        MESS '******************************************************' ;
        ERRE 5 ;
        QUIT PRODT ;
    FINS ;
 
    B_TN = VRAI ;
FINS ;
 
 
 
 
* +====================================================================+
* |                                                                    |
* |               C A L C U L   D E S   R É S U L T A T S              |
* |                                                                    |
* +====================================================================+
 
* Extraction des composantes du tenseur gradient
* ==============================================
 
UNX  = EXCO UN (EXTR LCO 1) ;
GDUX = KOPS 'GRADS' UNX $MD ;
DUDX = EXCO GDUX 'UX' ;
DUDY = EXCO GDUX 'UY' ;
 
SI ((NCO > 1) ET (NDIM > 1)) ;
    UNY  = EXCO UN (EXTR LCO 2) ;
    GDUY = KOPS 'GRADS' UNY $MD ;
    DVDX = EXCO GDUY 'UX' ;
    DVDY = EXCO GDUY 'UY' ;
 
    SI ((NCO > 2) ET (NDIM EGA 3)) ;
        UNZ  = EXCO UN (EXTR LCO 3) ;
        GDUZ = KOPS 'GRADS' UNZ $MD ;
        DWDX = EXCO GDUZ 'UX' ;
        DWDY = EXCO GDUZ 'UY' ;
        DWDZ = EXCO GDUZ 'UZ' ;
    SINON ;
        DWDX = 0. ;
        DWDY = 0. ;
        DWDZ = 0. ;
    FINS ;
 
    SI (NDIM EGA 3) ;
        DUDZ = EXCO GDUX 'UZ' ;
        DVDZ = EXCO GDUY 'UZ' ;
    SINON ;
        DUDZ = 0. ;
        DVDZ = 0. ;
    FINS ;
 
SINON ;
    SI (NDIM EGA 3) ;
        DUDZ = EXCO GDUX 'UZ' ;
    SINON ;
        DUDZ = 0. ;
    FINS ;
    DVDX = 0. ;
    DVDY = 0. ;
    DVDZ = 0. ;
    DWDX = 0. ;
    DWDY = 0. ;
    DWDZ = 0. ;
FINS ;
 
 
 
* Calcul du carré des normes des parties symétrique et anti-symétrique
* ====================================================================
 
A2 = (DUDY*DUDY) + (DUDZ*DUDZ) + (DVDX*DVDX) + (DVDZ*DVDZ)
                               + (DWDX*DWDX) + (DWDY*DWDY) ;
A3 = (DUDY*DVDX) + (DUDZ*DWDX) + (DVDZ*DWDY) ;
SI (NON B_TOROT) ;
    A1 = (DUDX*DUDX) + (DVDY*DVDY) + (DWDZ*DWDZ) ;
FINS ;
 
 
* Mot-clé 'TODEF': taux de déformation => 2.SijSij (par défaut)
SI (B_TODEF) ;
    P = (2. * (A1 + A3)) + A2 ;
FINS ;
 
* Mot-clé 'TOROT': taux de rotation => 2.RijRij
SI (B_TOROT) ;
*   /!\ Valeur absolue pour éviter les zéros "négatifs"
    P = ABS ( A2 + (-2. * A3) ) ;
FINS ;
 
* Mot-clé 'COMPL': tenseur gradient complet => UijUij
SI (B_COMPL) ;
    P = A1 + A2 ;
FINS ;
 
* Mot-clé 'MIXTE': Expression proposée par Daclès-Mariani et al.
SI (B_MIXTE) ;
    CPROD = 2. ;
    DRR = (ABS (A2 + (-2. * A3))) ** 0.5 ;
    DSS = ((2. * (A1 + A3)) + A2) ** 0.5 ;
    DMX = KOPS (DSS - DRR) '>|' 0. ;
    P = (DRR + (CPROD * DMX)) ** 2 ;
FINS ;
 
 
 
 
* Prise en compte de la contribution des forces de flottabilité
* =============================================================
* (Production turbulente moindre si stratification stable)
 
SI (B_TN) ;
*   Prandtl turbulent
    SGT = 0.7 ;
 
*   Calcul du produit scalaire gB.Grad(T)
    GDT  = KOPS 'GRADS' TN $MD ;
    DTDX = EXCO GDT 'UX' ;
    DTDY = EXCO GDT 'UY' ;
    AX   = COOR 1 GB ;
    AY   = COOR 2 GB ;
    G    = (DTDX*AX) + (DTDY*AY) ;
    SI (EGA NDIM 3) ;
        DTDZ = EXCO GDT 'UZ' ;
        AZ   = COOR 3 GB ;
        G    = G + (DTDZ*AZ) ;
    FINS ;
 
*   G > 0 => dT/dz < 0 => Stratification instable
*   G < 0 => dT/dz > 0 => Stratification stable
    IG = MASQ G 'INFERIEUR' 0. ;
    G  = G * IG * (1./SGT) ;
 
*   Modélisation selon W.Rodi
    C3  = 0.8 ;
    PG  = P + G ;
    IPG = MASQ (ABS PG) 'INFERIEUR' 1.E-30 ;
    RF  = (-1.) * G * (INVE (PG + (IPG*1.E-30))) ;
    P   = PG * (1. + (C3 * RF)) ;
FINS ;
 
 
* Empêche une éventuelle division par zéro
P = P + 1.E-30 ;
 
RESP P ;
FINP ;