* fichier : couplage_thermique2.dgibi ************************************************************************ * Section : Thermique Transitoire ************************************************************************ * *------------ 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' ;