* fichier : konv_impl_murs.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ *********************************************************** *********************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la **** **** solution des **** **** Equations d'Euler pour un gaz parfait **** **** OPERATEURS PRIM, PRET, KONV **** **** Implicit: calcul du jacobien du residu **** **** **** **** Cas gaz monoespece, "calorically perfect" **** **** **** **** Methodes: VLH **** **** **** **** A. BECCANTINI DRN/DMT/SEMT/LTMF AOUT 2000 **** *********************************************************** *********************************************************** 'OPTION' 'DIME' 2 ; 'OPTION' 'ELEM' QUA4 ; 'OPTION' 'ECHO' 0 ; 'OPTION' 'TRAC' 'X' ; * *** GRAPH * GRAPH = FAUX ; * GRAPH = VRAI ; ERRTOL = 1.0D-3 ; *************************** ***** DOMAINE SPATIAL **** *************************** A1 = 0.0D0 0.0D0; A2 = 1.0D0 0.0D0; A3 = 1.0D0 1.0D0; A4 = 0.0D0 1.0D0; L12 = A1 'DROIT' 1 A2; L23 = A2 'DROIT' 1 A3; L34 = A3 'DROIT' 1 A4; L41 = A4 'DROIT' 1 A1; DOM10 = 'DALL' L12 L23 L34 L41 'PLANE'; * *** Etat * ro = 1.11 ; p = 1234.7; gam = 1.4D0; cson = (gam '*' p '/' ro) '**' 0.5 ; un = 0.2 * cson ; ut = 0.3 * cson ; * *** retgd * ecin = 0.5D0 '*' ro '*' ((un '*' un ) '+' (ut '*' ut )); ret = (p '/' (gam '-' 1.0)) '+' ecin ; * *** flux en (n,t) * **************************************************** **************************************************** ******** Boucle sur les angles ********* **************************************************** **************************************************** DANGLE = (360.0 '/' 7.15) ; ANGLE = 0.0 ; 'REPETER' BLOC 8 ; * *** Rotation * ANGLE = ANGLE '+' DANGLE ; ORIG = 0.0D0 0.0D0; 'MESSAGE' ; 'MESSAGE' ('CHAINE' 'Angle de rotation= ' ANGLE); 'MESSAGE' ; DOMTOT = DOM10 'TOURNER' ANGLE ORIG; $DOMTOT = 'MODELISER' DOMTOT 'EULER'; TDOMTOT = 'DOMA' $DOMTOT 'VF'; MDOMTOT = TDOMTOT . 'QUAF' ; 'ELIM' MDOMTOT 1.E-6 ; PCON = ('DOMA' $DOMTOT 'CENTRE') 'POIN' 1 ; 'SI' GRAPH; 'TRACER' (('DOMA' $DOMTOT 'MAILLAGE' ) 'ET' ('DOMA' $DOMTOT 'CENTRE')) 'TITRE' 'Domaine et centre' ; 'FINSI' ; *********************** **** Les CHPOINTs **** *********************** ux = (un '*' ('COS' ANGLE)) '-' (ut '*' ('SIN' ANGLE)); uy = (un '*' ('SIN' ANGLE)) '+' (ut '*' ('COS' ANGLE)); ro0 = ro ; gnx0 = ro '*' ux ; gny0 = ro '*' uy ; ret0 = ret ; GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro0 ; GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx0 gny0) ; RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret0 ; *************************** **** L'operateur PRIM**** *************************** VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ; *************************** **** L'operateur PRET**** *************************** ORDESP = 1; ORDTEM = 1; ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM $DOMTOT RN VITESSE PRES GAMMAN ; *************************** **** L'operateur KONV**** *************************** METO = 'VLH' ; LISTINCO = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RETN' ; CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO $DOMTOT ROF VITF PF GAMF LISTINCO ; IJACO = 'KONV' 'VF' 'PERFMONO' 'JACOCONS' $DOMTOT LISTINCO METO RN VITESSE PRES GAMMAN ; *********************** ***** Le residu ******* *********************** DEBRN0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ; DEBGNX0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ; DEBGNY0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ; DEBRETN0 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ; * ***** Le jacobien du residu et le residu * UNRO = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 1.0 'RUXN' 0.0 'RUYN' 0.0 'RETN' 0.0 'NATURE' 'DISCRET' ; UNGX = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 0.0 'RUXN' 1.0 'RUYN' 0.0 'RETN' 0.0 'NATURE' 'DISCRET' ; UNGY = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 0.0 'RUXN' 0.0 'RUYN' 1.0 'RETN' 0.0 'NATURE' 'DISCRET' ; UNRET = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 4 'RN' 0.0 'RUXN' 0.0 'RUYN' 0.0 'RETN' 1.0 'NATURE' 'DISCRET' ; DDRHO = 'KOPS' IJACO 'MULT' UNRO ; * * ***** DDRHO contient: dRES_RN '/' dRN ('RN') ; * dRES_GXN '/' dRN ('RUXN') ; * dRES_GYN '/' dRN ('RUYN') ; * dRES_RETN '/' dRN ('RETN') ; * DDGX contient ... * * DDGX = 'KOPS' IJACO 'MULT' UNGX ; DDGY = 'KOPS' IJACO 'MULT' UNGY ; DDRET = 'KOPS' IJACO 'MULT' UNRET ; DRR = 'EXTRAIRE' DDRHO PCON 'RN' ; DGXR = 'EXTRAIRE' DDRHO PCON 'RUXN' ; DGYR = 'EXTRAIRE' DDRHO PCON 'RUYN' ; DRETR = 'EXTRAIRE' DDRHO PCON 'RETN' ; DRGX = 'EXTRAIRE' DDGX PCON 'RN' ; DGXGX = 'EXTRAIRE' DDGX PCON 'RUXN' ; DGYGX = 'EXTRAIRE' DDGX PCON 'RUYN' ; DRETGX = 'EXTRAIRE' DDGX PCON 'RETN' ; DRGY = 'EXTRAIRE' DDGY PCON 'RN' ; DGXGY = 'EXTRAIRE' DDGY PCON 'RUXN' ; DGYGY = 'EXTRAIRE' DDGY PCON 'RUYN' ; DRETGY = 'EXTRAIRE' DDGY PCON 'RETN' ; DRRET = 'EXTRAIRE' DDRET PCON 'RN' ; DGXRET = 'EXTRAIRE' DDRET PCON 'RUXN' ; DGYRET = 'EXTRAIRE' DDRET PCON 'RUYN' ; DRETRET = 'EXTRAIRE' DDRET PCON 'RETN' ; ********************************************************************* ********************************************************************* ***** On calcule les residues pour ro1 = ro0 * (1.'+' DELTA) ******* ********************************************************************* ********************************************************************* DELTA = 1.0D-4 ; ro1 = ro0 '*' (1.0 '+' DELTA) ; gnx1 = gnx0 ; gny1 = gny0 ; ret1 = ret0 ; GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ; GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ; RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ; VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ; ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM $DOMTOT RN VITESSE PRES GAMMAN ; CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO $DOMTOT ROF VITF PF GAMF LISTINCO ; *********************** ***** Le residu ******* *********************** DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ; DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ; DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ; DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ; * **** On calcule le jacobien numeriquement * DRRN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ro0) ; DGXRN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ro0) ; DGYRN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ro0) ; DRETRN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ro0) ; ERR1 = ('ABS' (DRRN '-' DRR)) '*' (ro0 '/' (ro0 '*' cson)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 1'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGXRN '-' DGXR)) '*' (ro0 '/' (ro0 '*' cson '*' cson)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 2'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGYRN '-' DGYR)) '*' (ro0 '/' (ro0 '*' cson '*' cson)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 3'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DRETRN '-' DRETR)) '*' (ro0 '/' (ret0 '*' cson)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 4'; 'ERREUR' 5 ; 'FINSI' ; *************************************************************************** *************************************************************************** ***** On calcule les residues pour gnx1 = gnx0 '+' (DELTA ro0 cson) ****** *************************************************************************** *************************************************************************** DELTA = 1.0D-4 ; ro1 = ro0 ; gnx1 = gnx0 '+' (DELTA '*' ro0 '*' cson) ; gny1 = gny0 ; ret1 = ret0 ; GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ; GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ; RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ; VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ; ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM $DOMTOT RN VITESSE PRES GAMMAN ; CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO $DOMTOT ROF VITF PF GAMF LISTINCO ; *********************** ***** Le residu ******* *********************** DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ; DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ; DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ; DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ; * **** On calcule le jacobien numeriquement * DRGXN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ro0 '*' cson) ; DGXGXN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ro0 '*' cson) ; DGYGXN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ro0 '*' cson) ; DRETGXN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ro0 '*' cson) ; ERR1 = ('ABS' (DRGXN '-' DRGX)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 5'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGXGXN '-' DGXGX)) '/' cson ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 6'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGYGXN '-' DGYGX)) '/' cson ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 7'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DRETGXN '-' DRETGX)) '*' (ro0 / ret0) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 8'; 'ERREUR' 5 ; 'FINSI' ; *************************************************************************** *************************************************************************** ***** On calcule les residues pour gny1 = gny0 '+' (DELTA ro0 cson) ****** *************************************************************************** *************************************************************************** DELTA = 1.0D-4 ; ro1 = ro0 ; gnx1 = gnx0 ; gny1 = gny0 '+' (DELTA '*' ro0 '*' cson) ; ret1 = ret0 ; GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ; GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ; RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ; VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ; ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM $DOMTOT RN VITESSE PRES GAMMAN ; CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO $DOMTOT ROF VITF PF GAMF LISTINCO ; *********************** ***** Le residu ******* *********************** DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ; DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ; DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ; DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ; * **** On calcule le jacobien numeriquement * DRGYN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ro0 '*' cson) ; DGXGYN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ro0 '*' cson) ; DGYGYN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ro0 '*' cson) ; DRETGYN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ro0 '*' cson) ; ERR1 = ('ABS' (DRGYN '-' DRGY)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 9'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGXGYN '-' DGXGY)) '/' cson ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 10'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGYGYN '-' DGYGY)) '/' cson ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 11'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DRETGYN '-' DRETGY)) '*' (ro0 / ret0) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 12'; 'ERREUR' 5 ; 'FINSI' ; *************************************************************************** *************************************************************************** ***** On calcule les residues pour ret1 = ret0 '*' (1. '+' DELTA) ****** *************************************************************************** *************************************************************************** DELTA = 1.0D-4 ; ro1 = ro0 ; gnx1 = gnx0 ; gny1 = gny0 ; ret1 = ret0 * (1. '+' DELTA) ; GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gam ; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ro1 ; GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (gnx1 gny1) ; RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ret1 ; VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ; ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM $DOMTOT RN VITESSE PRES GAMMAN ; CHPRES DT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO $DOMTOT ROF VITF PF GAMF LISTINCO ; *********************** ***** Le residu ******* *********************** DEBRN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 1) ; DEBGNX1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 2) ; DEBGNY1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 3) ; DEBRETN1 = 'EXTRAIRE' CHPRES PCON ('EXTRAIRE' LISTINCO 4) ; * **** On calcule le jacobien numeriquement * DRRETN = (DEBRN1 '-' DEBRN0) '/' (DELTA '*' ret0) ; DGXRETN = (DEBGNX1 '-' DEBGNX0) '/' (DELTA '*' ret0) ; DGYRETN = (DEBGNY1 '-' DEBGNY0) '/' (DELTA '*' ret0) ; DRETRETN = (DEBRETN1 '-' DEBRETN0) '/' (DELTA '*' ret0) ; ERR1 = ('ABS' (DRRETN '-' DRRET)) '*' (ret0 '/' (ro0 '*' cson)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 13'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGXRETN '-' DGXRET)) '*' (ret0 '/' (ro0 '*' cson * cson)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 14'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DGYRETN '-' DGYRET)) '*' (ret0 '/' (ro0 '*' cson * cson)) ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 15'; 'ERREUR' 5 ; 'FINSI' ; ERR1 = ('ABS' (DRETRETN '-' DRETRET)) '/' cson ; 'SI' (ERR1 > ERRTOL) ; 'MESSAGE' 'Problem 16'; 'ERREUR' 5 ; 'FINSI' ; 'FIN' BLOC ; 'FIN' ;