Télécharger la plaquette Cast3M


fr en

Télécharger xlap1b.eso

Retour à la liste

Numérotation des lignes : 

xlap1b

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

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