Télécharger la plaquette Cast3M


fr en

Télécharger ylap1b.eso

Retour à la liste

Numérotation des lignes : 

ylap1b

  1. C YLAP1B    SOURCE    CB215821  20/11/25    13:44:05     10792          
  2.       SUBROUTINE YLAP1B(ICOGRT,MPROC,MPVITC,MPTEMC,
  3.      $     MPVOLU,MPNORM,MPSURF,MELEFL,
  4.      $     KRFACE,KRCENT,LCLIMT,KRTIMP,LCLIMQ,KRQIMP,
  5.      $     LMIXT,KRMIXT,NOMINC,ICLAU,
  6.      $     KAPPA,CV,
  7.      $     IJACO,
  8.      $     IMPR,IRET)
  9.       IMPLICIT INTEGER(I-N)
  10.       IMPLICIT REAL*8 (A-H,O-Z)
  11. C***********************************************************************
  12. C NOM         : YLAP1B
  13. C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
  14. C               VF 2D.
  15. C               Ici, on ne calcule que les contributions à la matrice
  16. C               jacobienne faisant intervenir les coefficients pour le
  17. C               calcul des gradients de température (ICOGRT).
  18. C               (contributions à d Res_{\rho e_t} / d var
  19. C                var prenant successivement les valeurs :
  20. C                \rho,   \rho u,    \rho v,     \rho e_t)
  21. C
  22. C
  23. C
  24. C
  25. C LANGAGE     : ESOPE
  26. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  27. C               mél : gounand@semt2.smts.cea.fr
  28. C***********************************************************************
  29. C APPELES (UTIL)   : AJMTK : ajoute un objet de type MATSIM (non
  30. C                            standard) à un objet de type MATRIK.
  31. C APPELE PAR       : YLAP1A : Calcul de la matrice jacobienne du
  32. C                    résidu du laplacien VF 2D.
  33. C***********************************************************************
  34. C ENTREES            : ICOGRT (type MCHELM) : coefficients pour le
  35. C                        calcul du gradient de la température aux
  36. C                        interfaces.
  37. C                      MPROC (type MPOVAL) : masse volumique par
  38. C                                            élément.
  39. C                      MPVITC (type MPOVAL) : vitesse par élément.
  40. C                      MPTEMC (type MPOVAL) : température par élément.
  41. C                      MPVOLU (type MPOVAL) : volume des éléments.
  42. C                      MPNORM (type MPOVAL) : normale aux faces.
  43. C                      MPSURF (type MPOVAL) : surface des faces.
  44. C                      MELEFL (type MELEME) : connectivités face-(centre
  45. C                                             gauche, centre droit).
  46. C                      KRFACE (type MLENTI) : tableau de repérage dans
  47. C                        le maillage des faces des éléments.
  48. C                      KRCENT (type MLENTI) : tableau de repérage dans
  49. C                        le maillage des centres des éléments.
  50. C                      LCLIMT (type logique) : .TRUE. => CL de Dirichlet
  51. C                        sur la température.
  52. C                      KRTIMP (type MLENTI) : tableau de repérage dans
  53. C                        maillage des CL de Dirichlet sur la
  54. C                        température.
  55. C                      LCLIMQ (type logique) : .TRUE. => CL de Dirichlet
  56. C                        sur le flux de chaleur.
  57. C                      KRTIMP (type MLENTI) : tableau de repérage dans
  58. C                        maillage des CL de Dirichlet sur la
  59. C                        température.
  60. C                      KRQIMP (type MLENTI) : tableau de repérage dans
  61. C                        maillage des CL de Dirichlet sur le flux de
  62. C                        chaleur.
  63. C                      LMIXT (type logique) : .TRUE. => CL mixtes
  64. C                      KRMIXT (type MLENTI) : tableau de repérage dans
  65. C                        maillage des CL mixtes.
  66. C                      NOMINC (type MLMOTS) : noms des inconnues.
  67. C                      ICLAU               : option pour ne calculer
  68. C                                            que la thermique
  69. C                      KAPPA (type réel) : conductivité thermique (SI)
  70. C                      CV (type réel) : chaleur spécifique à volume
  71. C                                       constant (SI).
  72. C ENTREES/SORTIES    : IJACO (type MATRIK) : matrice jacobienne du
  73. C                        résidu du laplacien VF 2D.
  74. C SORTIES            : -
  75. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  76. C***********************************************************************
  77. C VERSION    : v1, 08/08/2001, version initiale
  78. C HISTORIQUE : v1, 08/08/2001, création
  79. C HISTORIQUE : 11/02/2003, ajout de l'option 'MIXT'
  80. C HISTORIQUE :
  81. C***********************************************************************
  82. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  83. C en cas de modification de ce sous-programme afin de faciliter
  84. C la maintenance !
  85. C***********************************************************************
  86.  
  87. -INC PPARAM
  88. -INC CCOPTIO
  89. -INC SMCHPOI
  90.       POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL,MPTEMC.MPOVAL
  91.       POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
  92. -INC SMCHAML
  93.       POINTEUR ICOGRT.MCHELM,JCOGRT.MCHAML
  94.       POINTEUR KCDTDX.MELVAL,KCDTDY.MELVAL
  95. -INC SMELEME
  96.       POINTEUR MELEFL.MELEME
  97.       POINTEUR MCOGRT.MELEME
  98. -INC SMLENTI
  99.       POINTEUR KRTIMP.MLENTI,KRQIMP.MLENTI,KRMIXT.MLENTI
  100.       POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
  101. -INC SMLMOTS
  102.       POINTEUR NOMINC.MLMOTS
  104.       POINTEUR IJACO.MATRIK
  105. *
  106. * Objet matrice élémentaire simplifié
  107. *
  108.       SEGMENT GMATSI
  109.       INTEGER POIPR1(NPP1,NEL1)
  110.       INTEGER POIDU1(1,NEL1)
  111.       INTEGER POIPR2(NPP2,NEL2)
  112.       INTEGER POIDU2(2,NEL2)
  113.       POINTEUR LMATSI(0).MATSIM
  114.       ENDSEGMENT
  115. *  Contributions de la part du gradient de température (CTGRT)
  116.       POINTEUR CTGRT.GMATSI
  117.       SEGMENT MATSIM
  118.       CHARACTER*8 NOMPRI,NOMDUA
  119.       REAL*8 VALMA1(1,NPP1,NEL1)
  120.       REAL*8 VALMA2(2,NPP2,NEL2)
  121.       ENDSEGMENT
  122.       POINTEUR RETRHO.MATSIM
  123.       POINTEUR RETROU.MATSIM
  124.       POINTEUR RETROV.MATSIM
  125.       POINTEUR RETRET.MATSIM
  126. *
  127.       REAL*8 KAPPA,CV
  128. *
  129.       INTEGER IMPR,IRET,ICLAU
  130. *
  131.       LOGICAL LCLIMT,LCLIMQ,LMIXT
  132.       LOGICAL LMUR
  133.       LOGICAL LCTRB1,LCTRB2
  134. *
  135.       INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
  136.       INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
  137.       INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
  138.       INTEGER NLCENP,NLCEND,NLFACE,NLCLQ,NLCLT
  139.       INTEGER NPTEL
  140. *
  141.       REAL*8 BETAX,BETAY,CNX,CNY
  142.       REAL*8 SIGNOR,SURFFA,VOLUEL
  143.       REAL*8 RHOP,UP,VP,TP
  144.       REAL*8 FACTOR
  145.       REAL*8 DTDRHO,DTDROU,DTDROV,DTDRET
  146. C
  147. *
  148. * Executable statements
  149. *
  150.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap1b.eso'
  151. *   On calcule la partie de d Res_{\rho e_t} / d var
  152. *   (var prend successivement les valeurs :
  153. *      \rho,   \rho u,    \rho v,     \rho e_t)
  154. *   faisant intervenir les coefficients pour le calcul des gradients de
  155. *   température (ICOGRT).
  156. *   C'est la partie contenant le terme : - \vect{q} \pscal \vect{n}
  157. *   Les noms de matrices élémentaires (type MATSIM) associées sont :
  158. *   RETRHO, RETROU, RETROV, RETRET
  159.       IF (LCLIMT) THEN
  160.          SEGACT KRTIMP
  161.       ENDIF
  162.       IF (LCLIMQ) THEN
  163.          SEGACT KRQIMP
  164.       ENDIF
  165.       IF (LMIXT) THEN
  166.          SEGACT KRMIXT
  167.       ENDIF
  168.       SEGACT NOMINC
  169.       SEGACT KRCENT
  170.       SEGACT KRFACE
  171.       SEGACT MELEFL
  172.       SEGACT MPSURF
  173.       SEGACT MPNORM
  174.       SEGACT MPVOLU
  175.       SEGACT MPROC
  176.       SEGACT MPVITC
  177.       SEGACT MPTEMC
  178.       SEGACT ICOGRT
  179.       NSOUCH=ICOGRT.IMACHE(/1)
  180.       DO 1 ISOUCH=1,NSOUCH
  181.          MCOGRT=ICOGRT.IMACHE(ISOUCH)
  182.          JCOGRT=ICOGRT.ICHAML(ISOUCH)
  183.          SEGACT JCOGRT
  184.          KCDTDX=JCOGRT.IELVAL(1)
  185.          KCDTDY=JCOGRT.IELVAL(2)
  186.          SEGDES JCOGRT
  187.          SEGACT KCDTDX
  188.          SEGACT KCDTDY
  189.          SEGACT MCOGRT
  190.          NELEM=MCOGRT.NUM(/2)
  191.          NPTEL=MCOGRT.NUM(/1)
  192.          NPP1=NPTEL-1
  193.          NPP2=NPTEL-1
  194.          NEL1=NELEM
  195.          NEL2=NELEM
  196.          IEL1=1
  197.          IEL2=1
  198.          SEGINI RETRHO
  199.          SEGINI RETROU
  200.          SEGINI RETROV
  201.          SEGINI RETRET
  202.          SEGINI CTGRT
  203.          RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  204.          RETRHO.NOMPRI(5:8)='    '
  205.          RETRHO.NOMDUA(1:4)=NOMINC.MOTS(4)
  206.          RETRHO.NOMDUA(5:8)='    '
  207.          RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
  208.          RETROU.NOMPRI(5:8)='    '
  209.          RETROU.NOMDUA(1:4)=NOMINC.MOTS(4)
  210.          RETROU.NOMDUA(5:8)='    '
  211.          RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
  212.          RETROV.NOMPRI(5:8)='    '
  213.          RETROV.NOMDUA(1:4)=NOMINC.MOTS(4)
  214.          RETROV.NOMDUA(5:8)='    '
  215.          RETRET.NOMPRI(1:4)=NOMINC.MOTS(4)
  216.          RETRET.NOMPRI(5:8)='    '
  217.          RETRET.NOMDUA(1:4)=NOMINC.MOTS(4)
  218.          RETRET.NOMDUA(5:8)='    '
  219.          DO 12 IELEM=1,NELEM
  220. *     Le premier point du support de ICOGRT est un point FACE
  221.             NGFACE=MCOGRT.NUM(1,IELEM)
  222.             NLFACE=KRFACE.LECT(NGFACE)
  223.             IF (NLFACE.EQ.0) THEN
  224.                WRITE(IOIMP,*) 'Erreur de programmation n°1'
  225.                GOTO 9999
  226.             ENDIF
  227. *     On calcule la contribution à la matrice jacobienne IJACO de la face
  228. *     NGFAC (points duaux     : centres à gauche et à droite de la face)
  229. *           (points primaux   : une partie (bicoz conditions aux limites)
  230. *            de ceux du stencil pour le calcul du gradient
  231. *            à la face, ils doivent être des points centres)
  232. *       Si le flux de chaleur sur la face est imposé par les conditions
  233. *       aux limites, la contribution de la face à IJACO est nulle.
  234.             LCTRB1=.TRUE.
  235.             IF (LCLIMQ) THEN
  236.                NLCLQ=KRQIMP.LECT(NGFACE)
  237.                IF (NLCLQ.NE.0) THEN
  238.                   LCTRB1=.FALSE.
  239.                ENDIF
  240.             ENDIF
  241.             IF (LCTRB1) THEN
  242.                NGCGAU=MELEFL.NUM(1,NLFACE)
  243.                NGCDRO=MELEFL.NUM(3,NLFACE)
  244.                LMUR=(NGCGAU.EQ.NGCDRO)
  245. *     On distingue le cas où la face est un bord du maillage (mur)
  246. *                  du cas où la face est interne au maillage
  247.                IF (.NOT.LMUR) THEN
  248.                   NPD=2
  249.                ELSE
  250.                   NPD=1
  251.                ENDIF
  252.                NPP=NPTEL-1
  253. *     IPD=1 : point à gauche du point NGFACE
  254. *     IPD=2 : point à droite du point NGFACE
  255.                DO 122 IPD=1,NPD
  256.                   NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
  257.                   IF (.NOT.LMUR) THEN
  258.                      CTGRT.POIDU2(IPD,IEL2)=NPDUAL
  259.                   ELSE
  260.                      CTGRT.POIDU1(IPD,IEL1)=NPDUAL
  261.                   ENDIF
  262.                   NLCEND=KRCENT.LECT(NPDUAL)
  263.                   IF (NLCEND.EQ.0) THEN
  264.                      WRITE(IOIMP,*) 'Erreur grave n°1'
  265.                      GOTO 9999
  266.                   ENDIF
  267.                   DO 124 IPP=1,NPP
  268.                      NPPRIM=MCOGRT.NUM(IPP+1,IELEM)
  269.                      LCTRB2=.TRUE.
  270.                      IF (LCLIMT) THEN
  271.                         NLCLT=KRTIMP.LECT(NPPRIM)
  272.                         IF (NLCLT.NE.0) THEN
  273.                            LCTRB2=.FALSE.
  274.                         ENDIF
  275.                      ENDIF
  276.                      IF (LCLIMQ) THEN
  277.                         NLCLQ=KRQIMP.LECT(NPPRIM)
  278.                         IF (NLCLQ.NE.0) THEN
  279.                            LCTRB2=.FALSE.
  280.                         ENDIF
  281.                      ENDIF
  282.                      IF (LMIXT) THEN
  283.                         NLCLM=KRMIXT.LECT(NPPRIM)
  284.                         IF (NLCLM.NE.0) THEN
  285.                            LCTRB2=.FALSE.
  286.                         ENDIF
  287.                      ENDIF
  288.                      IF (.NOT.LCTRB2) THEN
  289. *     Lorsque une contribution est nulle, on fixe artificiellement le
  290. *     point primal égal au point dual.
  291.                         IF (.NOT.LMUR) THEN
  292.                            CTGRT.POIPR2(IPP,IEL2)=NPDUAL
  293.                            RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
  294.                            RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
  295.                            RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
  296.                            RETRET.VALMA2(IPD,IPP,IEL2)=0.D0
  297.                         ELSE
  298.                            CTGRT.POIPR1(IPP,IEL1)=NPDUAL
  299.                            RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
  300.                            RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
  301.                            RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
  302.                            RETRET.VALMA1(IPD,IPP,IEL1)=0.D0
  303.                         ENDIF
  304.                      ELSE
  305. *     Les contributions valent :
  306. *     (d Res_{\rho e_t})_d / (d var)_p =
  307. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f) * \kappa
  308. *     * [ ((n_x * \beta_x) +(n_y * \beta_y)) *
  309. *         ((dT)_p / (d var)_p)]
  310. *     avec :
  311. *          (dT)_p / (d \rho)_p = (1 / \rho_p) *
  312. *                ( (((u_p)^2 + (v_p)^2) / (2 * c_v)) - T )
  313. *          (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v)
  314. *          (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v)
  315. *          (dT)_p / (d \rho e_t)_p =   1 / (\rho_p * c_v)
  316. *     \beta_x : coefficients pour le calcul de dT/dx
  317. *     \beta_y : coefficients pour le calcul de dT/dy
  318. *
  319.                         NLCENP=KRCENT.LECT(NPPRIM)
  320.                         IF (NLCENP.EQ.0) THEN
  321.                            WRITE(IOIMP,*) 'Erreur grave n°2'
  322.                            GOTO 9999
  323.                         ENDIF
  324. *       normale sortante pour IPD=1, rentrante pour IPD=2
  325.                         SIGNOR=(-1.D0)**(IPD+1)
  326.                         VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
  327.                         SURFFA=MPSURF.VPOCHA(NLFACE,1)
  328.                         CNX   =MPNORM.VPOCHA(NLFACE,1)
  329.                         CNY   =MPNORM.VPOCHA(NLFACE,2)
  330.                         BETAX =KCDTDX.VELCHE(IPP+1,IELEM)
  331.                         BETAY =KCDTDY.VELCHE(IPP+1,IELEM)
  332.                         RHOP  =MPROC.VPOCHA(NLCENP,1)
  333.                         UP    =MPVITC.VPOCHA(NLCENP,1)
  334.                         VP    =MPVITC.VPOCHA(NLCENP,2)
  335.                         TP    =MPTEMC.VPOCHA(NLCENP,1)
  336.                         FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*KAPPA
  337.      $                       *((CNX*BETAX)+(CNY*BETAY))
  338.                         DTDRHO=((((UP*UP)+(VP*VP))/(2.D0*CV))-TP)/RHOP
  339.                         DTDROU=-UP/(RHOP*CV)
  340.                         DTDROV=-VP/(RHOP*CV)
  341.                         DTDRET=1.D0/(RHOP*CV)
  342.                         IF (.NOT.LMUR) THEN
  343.                            CTGRT.POIPR2(IPP,IEL2)=NPPRIM
  344.                            RETRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRHO
  345.                            RETROU.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROU
  346.                            RETROV.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROV
  347.                            RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET
  348.                         ELSE
  349.                            CTGRT.POIPR1(IPP,IEL1)=NPPRIM
  350.                            RETRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRHO
  351.                            RETROU.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROU
  352.                            RETROV.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROV
  353.                            RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET
  354.                         ENDIF
  355.                      ENDIF
  356.  124              CONTINUE
  357.  122           CONTINUE
  358.                IF (.NOT.LMUR) THEN
  359.                   IEL2=IEL2+1
  360.                ELSE
  361.                   IEL1=IEL1+1
  362.                ENDIF
  363.             ENDIF
  364.  12      CONTINUE
  365.          NPP1=NPTEL-1
  366.          NPP2=NPTEL-1
  367.          NEL1=IEL1-1
  368.          NEL2=IEL2-1
  369.          SEGADJ RETRHO
  370.          SEGADJ RETROU
  371.          SEGADJ RETROV
  372.          SEGADJ RETRET
  373.          SEGADJ CTGRT
  374.          IF (ICLAU.EQ.0) THEN
  375.          CTGRT.LMATSI(**)=RETRHO
  376.          CTGRT.LMATSI(**)=RETROU
  377.          CTGRT.LMATSI(**)=RETROV
  378.          ENDIF
  379.          CTGRT.LMATSI(**)=RETRET
  380. * On accumule les matrices résultantes dans IJACO
  381.          CALL AJMTK(CTGRT,IJACO,IMPR,IRET)
  382.          IF (IRET.NE.0) GOTO 9999
  383.          SEGSUP RETRHO
  384.          SEGSUP RETROU
  385.          SEGSUP RETROV
  386.          SEGSUP RETRET
  387.          SEGSUP CTGRT
  388. *
  389.          SEGDES MCOGRT
  390.          SEGDES KCDTDY
  391.          SEGDES KCDTDX
  392.  1    CONTINUE
  393.       SEGDES ICOGRT
  394.       SEGDES MPTEMC
  395.       SEGDES MPVITC
  396.       SEGDES MPROC
  397.       SEGDES MPVOLU
  398.       SEGDES MPNORM
  399.       SEGDES MPSURF
  400.       SEGDES MELEFL
  401.       SEGDES KRFACE
  402.       SEGDES KRCENT
  403.       SEGDES NOMINC
  404.       IF (LCLIMQ) THEN
  405.          SEGDES KRQIMP
  406.       ENDIF
  407.       IF (LCLIMT) THEN
  408.          SEGDES KRTIMP
  409.       ENDIF
  410.       IF (LMIXT) THEN
  411.          SEGDES KRMIXT
  412.       ENDIF
  413. *
  414. * Normal termination
  415. *
  416.       IRET=0
  417.       RETURN
  418. *
  419. * Format handling
  420. *
  421. *
  422. * Error handling
  423. *
  424.  9999 CONTINUE
  425.       IRET=1
  426.       WRITE(IOIMP,*) 'An error was detected in subroutine ylap1b'
  427.       RETURN
  428. *
  429. * End of subroutine YLAP1B
  430. *
  431.       END
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales