Accueil | Cast3M

* fichier :  couplage_thermique2.dgibi
************************************************************************
************************************************************************
*
*------------ Cas test couplage_thermique2.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
* à l'état stationnaire les solutions 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 = 'QUAF' ;
*
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 = -3. ; b3 = 2. ;
T0 = T1 + T2 / 2. ;
*
*--------------------------------------- Début de la procédure FRON1
*
* Première étape de l'algorithme de QUARTERONI : Imposition de la
* condition aux limites de Dirichlet à la frontière commune sur le
* probleme 1
*
DEBP FRON1 ;
ARGU RVX*TABLE ;
*
RV   = RVX . 'EQEX' ;
T2F  = 'KCHT' RVX.'DOMZ' 'SCAL' 'SOMMET' (RV . 'PROB2' . 'INCO' . 'T2');
T1F  = 'KCHT' RVX.'DOMZ' 'SCAL' 'SOMMET' (b3 - (b2*t2f) / b1) ;
*
MAIL1 = 'DOMA' RVX.'DOMZ' 'SOMMET' ;
CLIM1 = 'EXCO' RV.'CLIM' 'T1' 'T1';
CLIM2 = 'REDU' CLIM1 MAIL1 ;
CLIM3 = 'NOMC' T1F 'T1' ;
RV . 'CLIM' = RV . 'CLIM' - CLIM2 + CLIM3 ;
*
s1 ma1 = 'KOPS' 'MATRIK' ;
'RESPROC' s1 ma1 ;
'FINPROC' ;
*--------------------------------------- Fin de la procédure FRON1
*
*--------------------------------------- Début de la procédure FRON2
*
* Deuxième étape de l'algorithme de QUARTERONI : Imposition du flux
* (condition aux limites de Neumann) à la frontière commune sur le
* probleme 2
*
DEBP FRON2 ;
ARGU RVX*TABLE ;
*
RV = RVX . 'EQEX' ;
S2 MA2 = 'LAPN' RV . 'PROB1' . '2LAPN' ;
T1F    = 'NOMC' 'T1' (RV . 'PROB1' . 'INCO' . 'T1') ;
FLU1   = 'KOPS' MA2 'MULT' T1F ;
FLU1   = 'NOMC' 'T2' (FLU1 * -1.) ;
MAIL1  = 'DOMA' (RVX . 'DOMZ') 'SOMMET' ;
s1     = 'REDU' FLU1 MAIL1 ;
*
s0 ma1 = 'KOPS' 'MATRIK' ;
'RESPROC' s1 ma1 ;
'FINPROC' ;
*--------------------------------------- Fin de la procédure FRON2
*
*= MAILLAGE
*
XMIN = 0. ; X1 = XMIN + L1 ; X2 = X1 + L2 ;
YMIN = 0. ; DY = 1. ; Y1 = YMIN + DY ; Y2 = Y1 + DY ;
NX = 20 ; 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 tant que régime asymptotique d'un
* transitoire. La résolution s'appuie sur l'algorithme de Quarteroni
* (découplage des deux domaines et condition mixte Dirichlet/Neumann à
* l'interface entre les deux problèmes)
* La solution transitoire a un sens, à condition de réaliser un point
* fixe à chaque pas de temps. Sinon le transitoire est faux.
*-------------------------------------------------------------------
*
*= Description du problème (table EQEX)
*
* Description de chaque sous-problème
RV1 = 'EQEX' $DOM1 'ITMA' 1 'NITER' 1 'ALFA' 1.
             'OPTI' 'EF' 'IMPL'
             'ZONE' $FRON 'OPER' 'FRON1'
             'ZONE' $DOM1 'OPER' 'LAPN' D1 'INCO' 'T1' ;
RV2 = 'EQEX' $DOM2 'ITMA' 1 'NITER' 1 'ALFA' 1.
             'OPTI' 'EF' 'IMPL'
             'ZONE' $FRON 'OPER' 'FRON2'
             'ZONE' $DOM2 'OPER' 'LAPN' D2 'INCO' 'T2' ;
*
* Transitoire
RV1 = 'EQEX' RV1
             'OPTI' 'EF' 'IMPL' 'CENTREE'
             'ZONE' $DOM1 'OPER' 'DFDT' 1. 'T1' 1. 'INCO' 'T1' ;
RV2 = 'EQEX' RV2
             'OPTI' 'EF' 'IMPL' 'CENTREE'
             'ZONE' $DOM2 'OPER' 'DFDT' 1. 'T2' 1. 'INCO' 'T2' ;
*
* Conditions aux limites
RV1 = 'EQEX' RV1
            'CLIM' 'T1' 'TIMP' DGAU  T1
                   'T1' 'TIMP' FRON  -99. ;
RV2 = 'EQEX' RV2
            'CLIM' 'T2' 'TIMP' DDRO  T2 ;
*
*= Initialisation de la table INCO et résolution
*
RV1 . 'INCO' . 'T1'  = 'KCHT' $DOM1 SCAL SOMMET T0 ;
RV2 . 'INCO' . 'T2'  = 'KCHT' $DOM2 SCAL SOMMET T0 ;
RV1 . 'INCO' . 'T1N' = 'KCHT' $DOM1 SCAL SOMMET T0 ;
RV2 . 'INCO' . 'T2N' = 'KCHT' $DOM2 SCAL SOMMET T0 ;
RV1 . 'PROB2' = RV2 ;
RV2 . 'PROB1' = RV1 ;
*
'REPE' kartroni 25 ;
     EXEC RV1 ;
     EXEC RV2 ;
'FIN' kartroni ;
*
*= Post-traitement de la solution calculée
*
EVC1 = 'EVOL' 'CHPO' DBA1 RV1.'INCO'.'T1' ;
EVC2 = 'EVOL' 'CHPO' DBA2 RV2.'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' RV1.'INCO'.'T1' DOM1 ;
CHAM2 = 'CHAN' 'CHAM' RV2.'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' ;