C YLA1BT SOURCE CB215821 20/11/25 13:43:59 10792 C YLAP1B SOURCE LEPOTIER 02/12/17 21:19:20 4510 SUBROUTINE YLA1BT(ICOGRT,MPTEMC, $ MPVOLU,MPNORM,MPSURF,MELEFL, $ KRFACE,KRCENT,LCLIMT,KRTIMP,LCLIMQ,KRQIMP, $ LMIXT,KRMIXT,IJACO, $ IMPR,IRET) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C NOM : YLA1BT C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien C VF 2D. C Ici, on ne calcule que les contributions à la matrice C jacobienne faisant intervenir les coefficients pour le C calcul des gradients de température (ICOGRT). C (contributions à d Res_{\rho e_t} / d var C var prenant successivement les valeurs : C \rho, \rho u, \rho v, \rho e_t) C C C C C LANGAGE : ESOPE C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) C mél : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES (UTIL) : AJMTK : ajoute un objet de type MATSIM (non C standard) à un objet de type MATRIK. C APPELE PAR : YLAP1A : Calcul de la matrice jacobienne du C résidu du laplacien VF 2D. C*********************************************************************** C ENTREES : ICOGRT (type MCHELM) : coefficients pour le C calcul du gradient de la température aux C interfaces. C MPROC (type MPOVAL) : masse volumique par C élément. C MPVITC (type MPOVAL) : vitesse par élément. C MPTEMC (type MPOVAL) : température par élément. C MPVOLU (type MPOVAL) : volume des éléments. C MPNORM (type MPOVAL) : normale aux faces. C MPSURF (type MPOVAL) : surface des faces. C MELEFL (type MELEME) : connectivités face-(centre C gauche, centre droit). C KRFACE (type MLENTI) : tableau de repérage dans C le maillage des faces des éléments. C KRCENT (type MLENTI) : tableau de repérage dans C le maillage des centres des éléments. C LCLIMT (type logique) : .TRUE. => CL de Dirichlet C sur la température. C KRTIMP (type MLENTI) : tableau de repérage dans C maillage des CL de Dirichlet sur la C température. C LCLIMQ (type logique) : .TRUE. => CL de Dirichlet C sur le flux de chaleur. C KRTIMP (type MLENTI) : tableau de repérage dans C maillage des CL de Dirichlet sur la C température. C KRQIMP (type MLENTI) : tableau de repérage dans C maillage des CL de Dirichlet sur le flux de C chaleur. C LMIXT (type logique) : .TRUE. => CL mixtes C KRMIXT (type MLENTI) : tableau de repérage dans C maillage des CL mixtes. C NOMINC (type MLMOTS) : noms des inconnues. C ICLAU : option pour ne calculer C que la thermique C KAPPA (type réel) : conductivité thermique (SI) C CV (type réel) : chaleur spécifique à volume C constant (SI). C ENTREES/SORTIES : IJACO (type MATRIK) : matrice jacobienne du C résidu du laplacien VF 2D. C SORTIES : - C CODE RETOUR (IRET) : = 0 si tout s'est bien passé C*********************************************************************** C VERSION : v1, 08/08/2001, version initiale C HISTORIQUE : v1, 08/08/2001, création C HISTORIQUE : 11/02/2003, ajout de l'option 'MIXT' C HISTORIQUE : C*********************************************************************** C Prière de PRENDRE LE TEMPS de compléter les commentaires C en cas de modification de ce sous-programme afin de faciliter C la maintenance ! C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC SMCHPOI POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL,MPTEMC.MPOVAL POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL -INC SMCHAML POINTEUR ICOGRT.MCHELM,JCOGRT.MCHAML POINTEUR KCDTDX.MELVAL,KCDTDY.MELVAL -INC SMELEME POINTEUR MELEFL.MELEME POINTEUR MCOGRT.MELEME -INC SMLENTI POINTEUR KRTIMP.MLENTI,KRQIMP.MLENTI,KRMIXT.MLENTI POINTEUR KRCENT.MLENTI,KRFACE.MLENTI -INC SMLMOTS POINTEUR NOMINC.MLMOTS -INC SMMATRIK POINTEUR IJACO.MATRIK * * Objet matrice élémentaire simplifié * SEGMENT GMATSI INTEGER POIPR1(NPP1,NEL1) INTEGER POIDU1(1,NEL1) INTEGER POIPR2(NPP2,NEL2) INTEGER POIDU2(2,NEL2) POINTEUR LMATSI(0).MATSIM ENDSEGMENT * Contributions de la part du gradient de température (CTGRT) POINTEUR CTGRT.GMATSI SEGMENT MATSIM CHARACTER*8 NOMPRI,NOMDUA REAL*8 VALMA1(1,NPP1,NEL1) REAL*8 VALMA2(2,NPP2,NEL2) ENDSEGMENT POINTEUR RETRET.MATSIM * REAL*8 KAPPA,CV * INTEGER IMPR,IRET,ICLAU * LOGICAL LCLIMT,LCLIMQ,LMIXT LOGICAL LMUR LOGICAL LCTRB1,LCTRB2 * INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2 INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2 INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL INTEGER NLCENP,NLCEND,NLFACE,NLCLQ,NLCLT INTEGER NPTEL * REAL*8 BETAX,BETAY,CNX,CNY REAL*8 SIGNOR,SURFFA,VOLUEL REAL*8 RHOP,UP,VP,TP REAL*8 FACTOR REAL*8 DTDRHO,DTDROU,DTDROV,DTDRET C * * Executable statements * IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap1b.eso' * On calcule la partie de d Res_{\rho e_t} / d var * (var prend successivement les valeurs : * \rho, \rho u, \rho v, \rho e_t) * faisant intervenir les coefficients pour le calcul des gradients de * température (ICOGRT). * C'est la partie contenant le terme : - \vect{q} \pscal \vect{n} * Les noms de matrices élémentaires (type MATSIM) associées sont : * RETRHO, RETROU, RETROV, RETRET IF (LCLIMT) THEN SEGACT KRTIMP ENDIF IF (LCLIMQ) THEN SEGACT KRQIMP ENDIF IF (LMIXT) THEN SEGACT KRMIXT ENDIF SEGACT KRCENT SEGACT KRFACE SEGACT MELEFL SEGACT MPSURF SEGACT MPNORM SEGACT MPVOLU SEGACT MPTEMC SEGACT ICOGRT NSOUCH=ICOGRT.IMACHE(/1) DO 1 ISOUCH=1,NSOUCH MCOGRT=ICOGRT.IMACHE(ISOUCH) JCOGRT=ICOGRT.ICHAML(ISOUCH) SEGACT JCOGRT KCDTDX=JCOGRT.IELVAL(1) SEGDES JCOGRT SEGACT KCDTDX SEGACT MCOGRT NELEM=MCOGRT.NUM(/2) NPTEL=MCOGRT.NUM(/1) NPP1=NPTEL-1 NPP2=NPTEL-1 NEL1=NELEM NEL2=NELEM IEL1=1 IEL2=1 c BOUCLE POUR EVALUER LE NOMBRE DE FACES FRONTIERES c POUR CHAQUE SOUS DOMAINE NFRON = 1 DO 1200 IELEM=1,NELEM NGFACE=MCOGRT.NUM(1,IELEM) NLFACE=KRFACE.LECT(NGFACE) IF (NLFACE.EQ.0) THEN WRITE(IOIMP,*) 'Erreur de programmation n1' GOTO 9999 ENDIF NGCGAU=MELEFL.NUM(1,NLFACE) NGCDRO=MELEFL.NUM(3,NLFACE) IF ( NGCGAU.EQ.NGCDRO) THEN NFRON = NFRON + 1 ENDIF 1200 CONTINUE c WRITE(6,*) 'NFRON= ',NFRON NEL1 = NFRON c WRITE(6,*) 'NPP1= ',NPP1,'NPP2= ',NPP2 c WRITE(6,*) 'NEL1= ',NEL1,'NEL2= ',NEL2 SEGINI RETRET SEGINI CTGRT RETRET.NOMPRI(1:4)='RETN' RETRET.NOMPRI(5:8)=' ' RETRET.NOMDUA(1:4)='RETN' RETRET.NOMDUA(5:8)=' ' DO 12 IELEM=1,NELEM * Le premier point du support de ICOGRT est un point FACE NGFACE=MCOGRT.NUM(1,IELEM) NLFACE=KRFACE.LECT(NGFACE) IF (NLFACE.EQ.0) THEN WRITE(IOIMP,*) 'Erreur de programmation n1' GOTO 9999 ENDIF * On calcule la contribution à la matrice jacobienne IJACO de la face * NGFAC (points duaux : centres à gauche et à droite de la face) * (points primaux : une partie (bicoz conditions aux limites) * de ceux du stencil pour le calcul du gradient * à la face, ils doivent être des points centres) * Si le flux de chaleur sur la face est imposé par les conditions * aux limites, la contribution de la face à IJACO est nulle. LCTRB1=.TRUE. c IF (LCLIMQ) THEN c NLCLQ=KRQIMP.LECT(NGFACE) c IF (NLCLQ.NE.0) THEN c LCTRB1=.FALSE. c ENDIF c ENDIF * WRITE(6,*) 'NGFACE= ',NGFACE,'LCTRB1= ', LCTRB1 IF (LCTRB1) THEN NGCGAU=MELEFL.NUM(1,NLFACE) NGCDRO=MELEFL.NUM(3,NLFACE) LMUR=(NGCGAU.EQ.NGCDRO) * On distingue le cas où la face est un bord du maillage (mur) * du cas où la face est interne au maillage IF (.NOT.LMUR) THEN NPD=2 ELSE NPD=1 ENDIF NPP=NPTEL-1 * IPD=1 : point à gauche du point NGFACE * IPD=2 : point à droite du point NGFACE DO 122 IPD=1,NPD NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE) IF (.NOT.LMUR) THEN CTGRT.POIDU2(IPD,IEL2)=NPDUAL ELSE CTGRT.POIDU1(IPD,IEL1)=NPDUAL ENDIF NLCEND=KRCENT.LECT(NPDUAL) IF (NLCEND.EQ.0) THEN WRITE(IOIMP,*) 'Erreur grave n1' GOTO 9999 ENDIF DO 124 IPP=1,NPP NPPRIM=MCOGRT.NUM(IPP+1,IELEM) LCTRB2=.TRUE. IF (LCLIMT) THEN NLCLT=KRTIMP.LECT(NPPRIM) IF (NLCLT.NE.0) THEN LCTRB2=.FALSE. ENDIF ENDIF IF (LCLIMQ) THEN NLCLQ=KRQIMP.LECT(NPPRIM) IF (NLCLQ.NE.0) THEN LCTRB2=.FALSE. ENDIF ENDIF IF (LMIXT) THEN NLCLM=KRMIXT.LECT(NPPRIM) IF (NLCLM.NE.0) THEN LCTRB2=.FALSE. ENDIF ENDIF IF (.NOT.LCTRB2) THEN * Lorsque une contribution est nulle, on fixe artificiellement le * point primal égal au point dual. IF (.NOT.LMUR) THEN CTGRT.POIPR2(IPP,IEL2)=NPDUAL RETRET.VALMA2(IPD,IPP,IEL2)=0.D0 ELSE CTGRT.POIPR1(IPP,IEL1)=NPDUAL RETRET.VALMA1(IPD,IPP,IEL1)=0.D0 ENDIF ELSE * Les contributions valent : * (d Res_{\rho e_t})_d / (d var)_p = * +/-1 (normale sortante, rentrante) (1/V_d) * (S_f) * \kappa * * [ ((n_x * \beta_x) +(n_y * \beta_y)) * * ((dT)_p / (d var)_p)] * avec : * (dT)_p / (d \rho)_p = (1 / \rho_p) * * ( (((u_p)^2 + (v_p)^2) / (2 * c_v)) - T ) * (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v) * (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v) * (dT)_p / (d \rho e_t)_p = 1 / (\rho_p * c_v) * \beta_x : coefficients pour le calcul de dT/dx * \beta_y : coefficients pour le calcul de dT/dy * NLCENP=KRCENT.LECT(NPPRIM) IF (NLCENP.EQ.0) THEN WRITE(IOIMP,*) 'Erreur grave n2' GOTO 9999 ENDIF * normale sortante pour IPD=1, rentrante pour IPD=2 SIGNOR=(-1.D0)**(IPD+1) VOLUEL=MPVOLU.VPOCHA(NLCEND,1) SURFFA=MPSURF.VPOCHA(NLFACE,1) CNX =MPNORM.VPOCHA(NLFACE,1) CNY =MPNORM.VPOCHA(NLFACE,2) BETAX =KCDTDX.VELCHE(IPP+1,IELEM) FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA $ *(BETAX) DTDRET=1.D0 IF (.NOT.LMUR) THEN CTGRT.POIPR2(IPP,IEL2)=NPPRIM RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET c WRITE(6,*) 'IPD=',IPD,'IPP=',IPP,'IEL2=',IEL2 c WRITE(6,*) 'FACT=',FACTOR*DTDRET c WRITE(6,*) 'SIGNOR=',SIGNOR,'VOLUEL=',VOLUEL, c & 'SURFFA=',SURFFA,'BETAX= ',BETAX ELSE CTGRT.POIPR1(IPP,IEL1)=NPPRIM RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET c WRITE(6,*) 'IPD=',IPD,'IPP=',IPP,'IEL1=',IEL1 c WRITE(6,*) 'FACT=',FACTOR*DTDRET c WRITE(6,*) 'SIGNOR=',SIGNOR,'VOLUEL=',VOLUEL, c & 'SURFFA=',SURFFA,'BETAX= ',BETAX ENDIF ENDIF 124 CONTINUE 122 CONTINUE IF (.NOT.LMUR) THEN IEL2=IEL2+1 ELSE IEL1=IEL1+1 ENDIF ENDIF 12 CONTINUE NPP1=NPTEL-1 NPP2=NPTEL-1 NEL1=IEL1-1 NEL2=IEL2-1 c WRITE(6,*) 'APRES SEGADJ' c WRITE(6,*) 'NPP1= ',NPP1,'NPP2= ',NPP2 c WRITE(6,*) 'NEL1= ',NEL1,'NEL2= ',NEL2 SEGDES MCOGRT SEGDES KCDTDX SEGADJ RETRET SEGADJ CTGRT CTGRT.LMATSI(**)=RETRET * On accumule les matrices résultantes dans IJACO CALL AJMTK(CTGRT,IJACO,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGSUP RETRET SEGSUP CTGRT * 1 CONTINUE SEGDES ICOGRT SEGDES MPTEMC SEGDES MPVOLU SEGDES MPNORM SEGDES MPSURF SEGDES MELEFL SEGDES KRFACE SEGDES KRCENT IF (LCLIMQ) THEN SEGDES KRQIMP ENDIF IF (LCLIMT) THEN SEGDES KRTIMP ENDIF IF (LMIXT) THEN SEGDES KRMIXT ENDIF * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine ylap1b' RETURN * * End of subroutine YLAP1B * END