Accueil | Cast3M

* fichier :  couplage_thermique3.dgibi
************************************************************************
************************************************************************
*
*------------ Cas test couplage_thermique3.dgibi -------------------
*
*
* Tests des opérateurs de Castem-fluide
*
*
*- Options générales
*
OPTI DIME 2 ELEM QUA8 ECHO 0 ;
GRAPH = 'N' ;
EPS0  = 1.D-6 ;
*
*-------------------------------------------------------------------
* Couplage thermique entre deux domaines solides l'échange entre les
* deux domaines étant réalisés en assurant la continuité du flux et
* une relation entre les champs inconnues à la frontière : on compare
* à les solutions stationnaires calculée et analytique.
*-------------------------------------------------------------------
* La géométrie est un barreau constitué de deux parties de longueur
* respective L1 et L2 et de diffusivité thermique D1 et D2. On a la
* relation b1*Tg + b2*Tg = b3 à l'interface, bi étant connus et Tg
* (resp. Td) représentant la valeur à la frontière coté gauche (resp.
* coté droit). On note T1(x) (resp. T2(x)) la température dans le
* premier barreau (resp. dans le deuxième).
*
* Chargement : Flux nul sur les frontières haute et basse (1D) et
* température imposées aux deux extrémités (notée T1 et T2).
*
* Solution stationnaire : Soit Tg (resp. Td) la température au niveau
* de l'interface entre les deux barreaux vue du premier barreau (resp.
* du deuxième). A l'état stationnaire, l'égalité des gradients et la
* condition de saut à l'interface donne :
*    Tg*det = (a1*T1 + a2*T2)*b2 - b3*a2
*    Td*det = b3*a1 - (a1*T1 + a2*T2)*b1
* avec a1=D1/L1, a2=D2/L2 et det=a1b2 - a2b1
* On a alors
*    T1(x) = T1 + (Tg-T1) x/L1 pour x in [0,L1]
*    T2(x) = T2 + (Td-T2) (L1+L2-x)/L2 pour x in [L1,L1+L2]
*
* Application numérique : D1=2 ; L1=1 ; T1=0 ; D2=3 ; L2=5 ; T2=1 ;
*                         Relation de saut : 5*Tg - 3*Td = 2
* On obtient alors a1=2 et a2=3/5; Tg=1/3 et Td=-1/9 ;
*    T1(x) = x/3
*    T2(x) = (2*x-3)/9
*
* ATTENTION : La convergence de l'algorithme de Quarteroni n'est pas
* garantie pour d'autres choix de data ! 
*
*-------------------------------------------------------------------
* Auteur : F.Dabbene                                           05/06
*-------------------------------------------------------------------
*
*= Données physico-numériques
*
* DISCR=LINE/MACRO/QUAF
DISCR = 'LINE' ;
*
L1 = 1. ; D1 = 2. ; T1 = 0. ; a1 = D1/L1 ;
L2 = 5. ; D2 = 3. ; T2 = 1. ; a2 = D2/L2 ;
*b1 = 7. ; b2 = 11. ; b3 = 2.  ;
*b1 =-5. ; b2 = 5.  ; b3 = 10. ;
*b1 = 1. ; b2 = -1. ; b3 = -2. ;
b1 = 5. ; b2 = -3. ; b3 = 2.  ;
T0 = T1 + T2 / 2. ;
*
*--------------------------------------- Début de la procédure KEXEC
*
'DEBP' kexec ;
*
************************  K E X E C ************************************
*
* Résolution implicite d'un problème de mécanique des fluides décrit
* par l'intermédiaire d'une table structurée par l'opérateur EQEX.
*
* E/  rv : Table décrivant le problème à résoudre
*  /S rv : A l'indice INCO, table contenant la solution au temps final
*
* L'utilisateur peut par l'intermédiaire de procédures personnelles 
* stocker les résultats aux temps intermédiaires, faire des tracés, ...
*
************************  K E X E C ************************************
*
'ARGUMENT' rv*'TABLE   ' ;
*
*========================
* TESTS et INITIALISATION
*========================
*
*
* On n'accepte que les formulations en modèle NAVIER_STOKES 
'SI' ('EGA'  (rv . 'NAVISTOK') 0) ;
   'MESS' 'Pas de modèle NAVIER_STOKES' ;
   'ERRE' 5 ;
'FINSI'  ;
*
* On n'accepte qu'un schéma en temps d'Euler
* ischt=0 si le schéma en temps n'est ni BDF2, ni BDF4 
ischt  = 0 ;
dimope = 'DIME' (rv . 'LISTOPER') ;
'REPETER' bloc0 dimope ;
      nomper  = 'EXTRAIRE' &bloc0 (rv . 'LISTOPER') ;
      notable = 'MOT' ('TEXTE' ('CHAINE' &bloc0 nomper)) ;
      'SI' ('EGA' nomper 'DFDT    ') ;
         ischt   = ischt + rv . notable . 'KOPT' . 'ISCHT' ;
      'FINSI' ;
'FIN' bloc0 ;
'SI' ('NEG' ischt 0) ;
   'MESS' 'Pas de schéma en temps autre que Euler' ;
   'ERRE' 5 ;
'FINSI'  ;
*
* Coefficient de relaxation mis à 1 si non initialisé
'SI' ('NON' ('EXISTE' rv 'OMEGA')) ;
   'MESS' '** WARNING ** OMEGA NON défini --> OMEGA=1' ;
   omeg = 1.D0 ;
'SINON' ;
   omeg = rv . 'OMEGA' ;  
'FINSI' ;
*
* Création de la table pour historique
'SI' ('NON' ('EXISTE' rv 'HIST')) ;
   rv . 'HIST' = 'TABLE' ;
'FINSI' ;
*
* IMPKRES : Niveau d'impression pour KRES
* IMPTCRR : Fréquence d'impression pour TCRR (affichage résidu)
IMPKRES = 0 ;
IMPTCRR = RV . 'IMPR' ;
*
*========================
* RESOLUTION suivant EQEX
*========================
*
*
* 1
* Boucle en temps
*----------------
ITMA = rv . 'ITMA' ;
'SI' ('&lt;EG' ITMA 1) ; ITMA = 1 ; 'FINSI' ;
'REPETER' bloc1 ITMA ;
*
* 2
* Boucle de point fixe interne à un pas de temps
*-----------------------------------------------
'REPETER' bloci (rv . 'NITER') ;
* st mat : Second membre et matrice des opérateurs DFDT
* sf mau : Second membre er matrice des opérateurs sauf DFDT
st mat  = 'KOPS' 'MATRIK' ;
sf mau  = 'KOPS' 'MATRIK' ;
*
* 3
* Boucle sur les opérateurs de EQEX
*----------------------------------
'REPETER' bloc2 dimope ;
      nomper  = 'EXTRAIRE' &bloc2 (rv . 'LISTOPER') ;
      notable = 'MOT' ('TEXTE' ('CHAINE' &bloc2 nomper)) ;
      'SI' ('EGA' nomper 'DFDT    ') ;
         msi mai = ('TEXTE' nomper) (rv . notable) ;
         mat     = mat 'ET' mai ;
         st      = st  'ET' msi ;
      'SINO' ;
         msi mai = ('TEXTE' nomper) (rv . notable) ;
         mau     = mau 'ET' mai ;
         sf      = sf  'ET' msi ;
      'FINSI' ;
'FIN' bloc2 ;
*--------------------------------------
* Fin Boucle sur les opérateurs de EQEX
* 3
*
s2  = sf 'ET' st ;
ma1 = mau 'ET' mat ;
'SI' ('EXISTE' rv 'CLIM') ;
   s1 = rv . 'CLIM' ;
'SINON' ;
   s1 = ' ' ;
'FINSI'  ;
rv . 'S2' = s2 ;
rv . 'METHINV' . 'MATASS' = ma1 ;
rv . 'METHINV' . 'MAPREC' = ma1 ;
res = 'KRES'  ma1 'TYPI' (rv . 'METHINV')
                  'CLIM' s1
                  'SMBR' s2
                  'IMPR' IMPKRES ;
*'LIST' res ;
'SI' ('EXIS' (rv . 'INCO') 'LX') ;
   rv . 'INCO' . 'LXOK' = 1 ;
   rv . 'INCO' . 'LX' = 'EXCO' 'LX' res ;
'FINSI' ;
eps = 'TCRR' res omeg (rv . 'INCO') 'IMPR' IMPTCRR ;
'MENAGE' ;
'FIN' bloci ;
*---------------------------------------------------
* Fin Boucle de point fixe interne à un pas de temps
* 2
*
irt = 0 ;
'SI' ('EGA' (rv . 'ITMA') 0) ;
   irt = 'TCNM' rv 'NOUP';
'SINON'                                         ;
   irt = 'TCNM' rv ;
'FINSI' ;
'MENAGE' ;
'SI' ('EGA' irt 1) ;
   'MESSAGE' ' Temps final atteint : ' (rv . 'PASDETPS' . 'TPS') ;
   'QUITTER' bloc1 ;
'FINSI' ;
'FIN' bloc1 ;
*--------------------
* Fin Boucle en temps
* 1
*
************************  K E X E C ************************************
'FINPROC' ;
*--------------------------------------- Fin de la procédure KEXEC
*
*
*--------------------------------------- Début de la procédure KRELA
*
'DEBPROC' krela rvx*table ;
*
RV   = RVX . 'EQEX' ;
fron = 'DOMA' RVX.'DOMZ' 'MAILLAGE' ;
chp1 mat1 = 'KOPS' 'MATRIK' ;
*
* Imposition du flux pour chaque sous-problème
chp1 mat1 = 'LAPN' rv . '1LAPN' ;
tf1  = 'NOMC' 'T1' rv . 'INCO' . 'T1' ;
val1 = 'KOPS' mat1 'MULT' tf1 ;
val1 = 'NOMC' 'SCAL' val1 ;
chp2 mat2 = 'LAPN' rv . '2LAPN' ;
tf2  = 'NOMC' 'T2' rv . 'INCO' . 'T2' ;
val2 = 'KOPS' mat2 'MULT' tf2 ;
val2 = 'NOMC' 'SCAL' val2 ;
v1 = val1 - val2 / 2. ;
v2 = val2 - val1 / 2. ;
rv . 'INCO' . 'FLU1' = 'KCHT' RVX.'DOMZ' 'SCAL' 'SOMMET' v1 ;
rv . 'INCO' . 'FLU2' = 'KCHT' RVX.'DOMZ' 'SCAL' 'SOMMET' v2 ;
*
rela1 = 'RELA' b1 'UX' fron + b2 'UY' fron ;
smr1  = 'DEPI' rela1 b3 ;
rela2 = 'KOPS' 'RIMA' rela1 ;
rela3 = 'KOPS' 'CHANINCO' rela2
               ('MOTS' 'LX'  'UX' 'UY')
               ('MOTS' 'LX'  'T1' 'T2')
               ('MOTS' 'FLX' 'FX' 'FY')
               ('MOTS' 'LX'  'T1' 'T2')
;
smr3 = 'EXCO' smr1 'FLX' 'LX' ;
*
'FINP' smr3 rela3 ;
*--------------------------------------- Fin de la procédure KRELA
*
*= MAILLAGE
*
XMIN = 0. ; X1 = XMIN + L1 ; X2 = X1 + L2 ;
YMIN = 0. ; DY = 1. ; Y1 = YMIN + DY ; Y2 = Y1 + DY ;
NX = 5 ; NY = 1 ;
*
*------------------------------------------ Points du premier domaine
P1 = XMIN YMIN ;
P2 = X1   YMIN ;
P3 = X1   Y2   ;
P4 = XMIN Y2   ;
*------------------------------------------ Points du deuxième domaine
PD1 = P2 ;
PD2 = X2  YMIN ;
PD3 = X2  Y2   ;
PD4 = P3       ;
*------------------------------------------ Lignes du premier domaine
P1P2 = P1 'DROIT' NX P2 ;
P2P3 = P2 'DROIT' NY P3 ;
P3P4 = P3 'DROIT' NX P4 ;
P4P1 = P4 'DROIT' NY P1 ;
*------------------------------------------ Lignes du deuxième domaine
PD1PD2 = PD1 'DROIT' NX PD2 ;
PD2PD3 = PD2 'DROIT' NY PD3 ;
PD3PD4 = PD3 'DROIT' NX PD4 ;
PD4PD1 = 'INVERSE' P2P3 ;
*------------------------------------------ Maillages
*
* Afin d'éviter toute confusion de maillages, ces objets sont à écraser
* après la création des modèles (à remplacer par le maillage qui
* soustend le modèle)
DOM1 = 'DALLER' P1P2 P2P3 P3P4 P4P1 ;
DOM2 = 'DALLER' PD1PD2 PD2PD3 PD3PD4 PD4PD1 ;
DGAU = P4P1 ;
DDRO = PD2PD3 ;
DBA1 = P1P2 ;
DBA2 = PD1PD2 ;
*
*= Création des MODELES
*
DDOM1 = 'CHANGER' DOM1 'QUAF' ;
DDOM2 = 'CHANGER' DOM2 'QUAF' ;
DDOMT = DDOM1 'ET' DDOM2 ;
DDBA1 = 'CHANGER' DBA1 'QUAF' ;
DDBA2 = 'CHANGER' DBA2 'QUAF' ;
DDGAU = 'CHANGER' DGAU 'QUAF' ;
DDDRO = 'CHANGER' DDRO 'QUAF' ;
DFRON = 'CHANGER' P2P3 'QUAF' ;
*
$DOM1 = 'MODELISER' DDOM1 'NAVIER_STOKES' DISCR ;
$DOM2 = 'MODELISER' DDOM2 'NAVIER_STOKES' DISCR ;
$DOMT = 'MODELISER' DDOMT 'NAVIER_STOKES' DISCR ;
$DBA1 = 'MODELISER' DDBA1 'NAVIER_STOKES' DISCR ;
$DBA2 = 'MODELISER' DDBA2 'NAVIER_STOKES' DISCR ;
$DGAU = 'MODELISER' DDGAU 'NAVIER_STOKES' DISCR ;
$DDRO = 'MODELISER' DDDRO 'NAVIER_STOKES' DISCR ;
$FRON = 'MODELISER' DFRON 'NAVIER_STOKES' DISCR ;
*
DOM1 = 'DOMA' $DOM1 'MAILLAGE' ;
DOM2 = 'DOMA' $DOM2 'MAILLAGE' ;
DOMT = 'DOMA' $DOMT 'MAILLAGE' ;
DGAU = 'DOMA' $DGAU 'MAILLAGE' ;
DDRO = 'DOMA' $DDRO 'MAILLAGE' ;
FRON = 'DOMA' $FRON 'MAILLAGE' ;
DBA1 = 'DOMA' $DBA1 'MAILLAGE' ;
DBA2 = 'DOMA' $DBA2 'MAILLAGE' ;
*
*= Solution analytique
*
DET0 = a1*b2 - (a2*b1) ; 
DET1 = a1*T1 + (a2*T2) * b2 - (a2*b3) ;
DET2 = a1*T1 + (a2*T2) * b1 * -1. + (a1*b3) ;
TG   = DET1 / DET0 ;
TD   = DET2 / DET0 ;
D1X = 'EXTR' ('EVOL' 'CHPO' DBA1 ('COOR' 1 DBA1)) 'ORDO' 1 ;
D2X = 'EXTR' ('EVOL' 'CHPO' DBA2 ('COOR' 1 DBA2)) 'ORDO' 1 ;
X1X = D1X / L1 ;
NBPT = 'NBEL' ('CHAN' 'POI1' DBA1) ;
X2X = D2X - ('PROG' NBPT*L1) / L2 ;
T1X = (TG-T1)*X1X + ('PROG' NBPT*T1) ;
T2X = (T2-TD)*X2X + ('PROG' NBPT*TD) ;
EV1 = 'EVOL' 'ROUG' 'MANU' D1X T1X ;
EV2 = 'EVOL' 'ROUG' 'MANU' D2X T2X ;
EV3 = EV1 ET EV2 ;
*
*
*
*
*                  MODELISATION DU PROBLEME
*
*-------------------------------------------------------------------
* L'état stationnaire est obtenu en résolvant le problème contraint par
* la relation sur les températures d'interface, la condition de flux
* continu étant imposée à chaque pas d'un processus itératif (sorte
* d'algorithme de Quarteroni sans découplage mais avec condition de
* Neumann/Neumann à l'interface entre les deux problèmes) 
*-------------------------------------------------------------------
*
*= Description du problème (table EQEX)
*
* Description de chaque sous-problème
RV = 'EQEX' $DOMT 'ITMA' 1 'NITER' 200 'ALFA' 1.
            'OPTI' 'EF' 'IMPL'
            'ZONE' $DOM1 'OPER' 'LAPN' D1 'INCO' 'T1'
            'ZONE' $DOM2 'OPER' 'LAPN' D2 'INCO' 'T2' ;
*
* Traitement de l'interface
RV = 'EQEX' RV
            'OPTI' 'EF' 'IMPL' 'CENTREE'
            'ZONE' $FRON 'OPER' 'FIMP' 'FLU1' 'INCO' 'T1'
            'ZONE' $FRON 'OPER' 'FIMP' 'FLU2' 'INCO' 'T2'
            'ZONE' $FRON 'OPER' KRELA
            ;
*
* Conditions aux limites (extrémités du domaine)
RV = 'EQEX' RV
            'CLIM' 'T1' 'TIMP' DGAU  T1
            'CLIM' 'T2' 'TIMP' DDRO  T2 ;
*
*= Initialisation de la table INCO et résolution
*
RV . 'INCO' . 'T1'   = 'KCHT' $DOM1 SCAL SOMMET T0 ;
RV . 'INCO' . 'T2'   = 'KCHT' $DOM2 SCAL SOMMET T0 ;
RV . 'INCO' . 'FLU1' = 'KCHT' $FRON SCAL SOMMET 0. ;
RV . 'INCO' . 'FLU2' = 'KCHT' $FRON SCAL SOMMET 0. ;
*
chp1 mat1 = 'KOPS' 'MATRIK' ;
RV . 'INCO' . 'LX'   = chp1 ;
*
*
* Méthode d'inversion du problème : méthode directe
rv . 'METHINV' . 'TYPINV'  = 1      ;
*-------------------- Utilisé en itératif seulement
rv . 'METHINV' . 'IMPINV'  = 0      ;
rv . 'METHINV' . 'NITMAX'  = 100    ;
rv . 'METHINV' . 'PRECOND' = 3      ;
rv . 'METHINV' . 'RESID'   = 1.e-6  ;
rv . 'METHINV' . 'FCPRECT' = 1      ;
rv . 'METHINV' . 'FCPRECI' = 1      ;
rv . 'METHINV' . 'PCMLAG'  = 'APR2' ;
*
KEXEC RV ;
*
*= Post-traitement de la solution calculée
*
EVC1 = 'EVOL' 'CHPO' DBA1 RV.'INCO'.'T1' ;
EVC2 = 'EVOL' 'CHPO' DBA2 RV.'INCO'.'T2' ;
CT1  = 'EXTR' EVC1 'ORDO' 1 ;
CT2  = 'EXTR' EVC2 'ORDO' 1 ;
EVC3 = ('EVOL' 'VERT' 'MANU' D1X CT1)
  'ET' ('EVOL' 'VERT' 'MANU' D2X CT2) ;
*
*= Tracés
*
'SI' ('EGA' GRAPH 'O') ;
CHAM1 = 'CHAN' 'CHAM' RV.'INCO'.'T1' DOM1 ;
CHAM2 = 'CHAN' 'CHAM' RV.'INCO'.'T2' DOM2 ;
'TRAC' $DOMT (CHAM1 'ET' CHAM2) ;
*
TAB1 = 'TABLE' ;
TAB1 . 'TITRE' = 'TABLE' ;
TAB1 . 'TITRE' . 2 = 'Solution castem' ;
TAB1 . 'TITRE' . 1 = '---------------' ;
TAB1 . 'TITRE' . 4 = 'Solution exacte' ;
TAB1 . 'TITRE' . 3 = '---------------' ;
TAB1 . 1 = 'MARQ CROI NOLI' ;
TAB1 . 2 = 'MARQ CROI NOLI' ;
'DESS' (EVC3 'ET' EV3) 'TITR' 'Temperature en fonction de x'
                       'TITX' 'x' 'TITY' 'T' 'MIMA'
                       'LEGE' TAB1 ;
'FINSI' ;
*
*= Controle erreur
*
SOM1 = 'ABS' (('MAXI' ('INTG' EV3)) - ('MAXI' ('INTG' EVC3))) ;
TEST = SOM1 '<' EPS0 ;
'SI' TEST ;
     'ERRE' 0 ;
'SINO' ;
     'MESS'  'Analytic Integration -  Computed = Difference '  ;
     'MESS'  ('MAXI' ('INTG' EV3))'-' ('MAXI' ('INTG' EVC3)) '=' SOM1 ;
     'ERRE' 5 ;
'FINSI' ;
'FIN' ;