Télécharger la plaquette Cast3M


fr en

Télécharger ylap2b.eso

Retour à la liste

Numérotation des lignes : 

ylap2b

  1. C YLAP2B    SOURCE    CB215821  20/11/25    13:44:14     10792          
  2.       SUBROUTINE YLAP2B(ICOGRT,MPROC,MPVITC,MPTEMC,
  3.      $     MPVOLU,MPNORM,MPSURF,MELEFL,
  4.      $     KRFACE,KRCENT,LCLIMT,KRTIMP,LCLIMQ,KRQIMP,
  5.      $     NOMINC,
  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         : YLAP2B
  13. C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
  14. C               VF 3D.
  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 w,   \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       : YLAP2A : Calcul de la matrice jacobienne du
  32. C                    résidu du laplacien VF 3D.
  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                      NOMINC (type MLMOTS) : noms des inconnues.
  64. C                      KAPPA (type réel) : conductivité thermique (SI)
  65. C                      CV (type réel) : chaleur spécifique à volume
  66. C                                       constant (SI).
  67. C ENTREES/SORTIES    : IJACO (type MATRIK) : matrice jacobienne du
  68. C                        résidu du laplacien VF 3D.
  69. C SORTIES            : -
  70. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  71. C***********************************************************************
  72. C VERSION    : v1, 28/08/2001, version initiale
  73. C HISTORIQUE : v1, 28/08/2001, création
  74. C HISTORIQUE :
  75. C HISTORIQUE :
  76. C***********************************************************************
  77. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  78. C en cas de modification de ce sous-programme afin de faciliter
  79. C la maintenance !
  80. C***********************************************************************
  81.  
  82. -INC PPARAM
  83. -INC CCOPTIO
  84. -INC SMCHPOI
  85.       POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL,MPTEMC.MPOVAL
  86.       POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
  87. -INC SMCHAML
  88.       POINTEUR ICOGRT.MCHELM,JCOGRT.MCHAML
  89.       POINTEUR KCDTDX.MELVAL,KCDTDY.MELVAL,KCDTDZ.MELVAL
  90. -INC SMELEME
  91.       POINTEUR MELEFL.MELEME
  92.       POINTEUR MCOGRT.MELEME
  93. -INC SMLENTI
  94.       POINTEUR KRTIMP.MLENTI,KRQIMP.MLENTI
  95.       POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
  96. -INC SMLMOTS
  97.       POINTEUR NOMINC.MLMOTS
  99.       POINTEUR IJACO.MATRIK
  100. *
  101. * Objet matrice élémentaire simplifié
  102. *
  103.       SEGMENT GMATSI
  104.       INTEGER POIPR1(NPP1,NEL1)
  105.       INTEGER POIDU1(1,NEL1)
  106.       INTEGER POIPR2(NPP2,NEL2)
  107.       INTEGER POIDU2(2,NEL2)
  108.       POINTEUR LMATSI(0).MATSIM
  109.       ENDSEGMENT
  110. *  Contributions de la part du gradient de température (CTGRT)
  111.       POINTEUR CTGRT.GMATSI
  112.       SEGMENT MATSIM
  113.       CHARACTER*8 NOMPRI,NOMDUA
  114.       REAL*8 VALMA1(1,NPP1,NEL1)
  115.       REAL*8 VALMA2(2,NPP2,NEL2)
  116.       ENDSEGMENT
  117.       POINTEUR RETRHO.MATSIM
  118.       POINTEUR RETROU.MATSIM
  119.       POINTEUR RETROV.MATSIM
  120.       POINTEUR RETROW.MATSIM
  121.       POINTEUR RETRET.MATSIM
  122. *
  123.       REAL*8 KAPPA,CV
  124. *
  125.       INTEGER IMPR,IRET
  126. *
  127.       LOGICAL LCLIMT,LCLIMQ
  128.       LOGICAL LMUR
  129.       LOGICAL LCTRB1
  130. *
  131.       INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
  132.       INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
  133.       INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
  134.       INTEGER NLCENP,NLCEND,NLFACE,NLCLQ,NLCLT
  135.       INTEGER NPTEL
  136. *
  137.       REAL*8 BETAX,BETAY,BETAZ,CNX,CNY,CNZ
  138.       REAL*8 SIGNOR,SURFFA,VOLUEL
  139.       REAL*8 RHOP,UP,VP,WP,TP
  140.       REAL*8 FACTOR
  141.       REAL*8 DTDRHO,DTDROU,DTDROV,DTDROW,DTDRET
  142. C
  143. *
  144. * Executable statements
  145. *
  146.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap2b.eso'
  147. *   On calcule la partie de d Res_{\rho e_t} / d var
  148. *   (var prend successivement les valeurs :
  149. *      \rho,   \rho u,    \rho v,  \rho w,   \rho e_t)
  150. *   faisant intervenir les coefficients pour le calcul des gradients de
  151. *   température (ICOGRT).
  152. *   C'est la partie contenant le terme : - \vect{q} \pscal \vect{n}
  153. *   Les noms de matrices élémentaires (type MATSIM) associées sont :
  154. *   RETRHO, RETROU, RETROV, RETROW, RETRET
  155.       IF (LCLIMT) THEN
  156.          SEGACT KRTIMP
  157.       ENDIF
  158.       IF (LCLIMQ) THEN
  159.          SEGACT KRQIMP
  160.       ENDIF
  161.       SEGACT NOMINC
  162.       SEGACT KRCENT
  163.       SEGACT KRFACE
  164.       SEGACT MELEFL
  165.       SEGACT MPSURF
  166.       SEGACT MPNORM
  167.       SEGACT MPVOLU
  168.       SEGACT MPROC
  169.       SEGACT MPVITC
  170.       SEGACT MPTEMC
  171.       SEGACT ICOGRT
  172.       NSOUCH=ICOGRT.IMACHE(/1)
  173.       DO 1 ISOUCH=1,NSOUCH
  174.          MCOGRT=ICOGRT.IMACHE(ISOUCH)
  175.          JCOGRT=ICOGRT.ICHAML(ISOUCH)
  176.          SEGACT JCOGRT
  177.          KCDTDX=JCOGRT.IELVAL(1)
  178.          KCDTDY=JCOGRT.IELVAL(2)
  179.          KCDTDZ=JCOGRT.IELVAL(3)
  180.          SEGDES JCOGRT
  181.          SEGACT KCDTDX
  182.          SEGACT KCDTDY
  183.          SEGACT KCDTDZ
  184.          SEGACT MCOGRT
  185.          NELEM=MCOGRT.NUM(/2)
  186.          NPTEL=MCOGRT.NUM(/1)
  187.          NPP1=NPTEL-1
  188.          NPP2=NPTEL-1
  189.          NEL1=NELEM
  190.          NEL2=NELEM
  191.          IEL1=1
  192.          IEL2=1
  193.          SEGINI RETRHO
  194.          SEGINI RETROU
  195.          SEGINI RETROV
  196.          SEGINI RETROW
  197.          SEGINI RETRET
  198.          SEGINI CTGRT
  199.          RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  200.          RETRHO.NOMPRI(5:8)='    '
  201.          RETRHO.NOMDUA(1:4)=NOMINC.MOTS(5)
  202.          RETRHO.NOMDUA(5:8)='    '
  203.          RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
  204.          RETROU.NOMPRI(5:8)='    '
  205.          RETROU.NOMDUA(1:4)=NOMINC.MOTS(5)
  206.          RETROU.NOMDUA(5:8)='    '
  207.          RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
  208.          RETROV.NOMPRI(5:8)='    '
  209.          RETROV.NOMDUA(1:4)=NOMINC.MOTS(5)
  210.          RETROV.NOMDUA(5:8)='    '
  211.          RETROW.NOMPRI(1:4)=NOMINC.MOTS(4)
  212.          RETROW.NOMPRI(5:8)='    '
  213.          RETROW.NOMDUA(1:4)=NOMINC.MOTS(5)
  214.          RETROW.NOMDUA(5:8)='    '
  215.          RETRET.NOMPRI(1:4)=NOMINC.MOTS(5)
  216.          RETRET.NOMPRI(5:8)='    '
  217.          RETRET.NOMDUA(1:4)=NOMINC.MOTS(5)
  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. C
  270. C AB : we do not check the BC anymore
  271. C
  272.                      NLCENP=KRCENT.LECT(NPPRIM)
  273.                      IF (NLCENP .EQ. 0) THEN
  274. *     Lorsque une contribution est nulle, on fixe artificiellement le
  275. *     point primal égal au point dual.
  276.                         IF (.NOT.LMUR) THEN
  277.                            CTGRT.POIPR2(IPP,IEL2)=NPDUAL
  278.                            RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
  279.                            RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
  280.                            RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
  281.                            RETROW.VALMA2(IPD,IPP,IEL2)=0.D0
  282.                            RETRET.VALMA2(IPD,IPP,IEL2)=0.D0
  283.                         ELSE
  284.                            CTGRT.POIPR1(IPP,IEL1)=NPDUAL
  285.                            RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
  286.                            RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
  287.                            RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
  288.                            RETROW.VALMA1(IPD,IPP,IEL1)=0.D0
  289.                            RETRET.VALMA1(IPD,IPP,IEL1)=0.D0
  290.                         ENDIF
  291.                      ELSE
  292. *     Les contributions valent :
  293. *     (d Res_{\rho e_t})_d / (d var)_p =
  294. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f) * \kappa
  295. *     * [ ((n_x * \beta_x) + (n_y * \beta_y) + (n_z * \beta_z)) *
  296. *         ((dT)_p / (d var)_p)]
  297. *     avec :
  298. *          (dT)_p / (d \rho)_p = (1 / \rho_p) *
  299. *                ( (((u_p)^2 + (v_p)^2 + (w_p)^2) / (2 * c_v)) - T )
  300. *          (dT)_p / (d \rho u)_p = - u_p / (\rho_p * c_v)
  301. *          (dT)_p / (d \rho v)_p = - v_p / (\rho_p * c_v)
  302. *          (dT)_p / (d \rho w)_p = - w_p / (\rho_p * c_v)
  303. *          (dT)_p / (d \rho e_t)_p =   1 / (\rho_p * c_v)
  304. *     \beta_x : coefficients pour le calcul de dT/dx
  305. *     \beta_y : coefficients pour le calcul de dT/dy
  306. *     \beta_z : coefficients pour le calcul de dT/dz
  307. *
  308. *       normale sortante pour IPD=1, rentrante pour IPD=2
  309.                         SIGNOR=(-1.D0)**(IPD+1)
  310.                         VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
  311.                         SURFFA=MPSURF.VPOCHA(NLFACE,1)
  312.                         CNX   =MPNORM.VPOCHA(NLFACE,1)
  313.                         CNY   =MPNORM.VPOCHA(NLFACE,2)
  314.                         CNZ   =MPNORM.VPOCHA(NLFACE,3)
  315.                         BETAX =KCDTDX.VELCHE(IPP+1,IELEM)
  316.                         BETAY =KCDTDY.VELCHE(IPP+1,IELEM)
  317.                         BETAZ =KCDTDZ.VELCHE(IPP+1,IELEM)
  318.                         RHOP  =MPROC.VPOCHA(NLCENP,1)
  319.                         UP    =MPVITC.VPOCHA(NLCENP,1)
  320.                         VP    =MPVITC.VPOCHA(NLCENP,2)
  321.                         WP    =MPVITC.VPOCHA(NLCENP,3)
  322.                         TP    =MPTEMC.VPOCHA(NLCENP,1)
  323.                         FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*KAPPA
  324.      $                       *((CNX*BETAX)+(CNY*BETAY)+(CNZ*BETAZ))
  325.                         DTDRHO=((((UP*UP)+(VP*VP)+(WP*WP))
  326.      $                       /(2.D0*CV))-TP)/RHOP
  327.                         DTDROU=-UP/(RHOP*CV)
  328.                         DTDROV=-VP/(RHOP*CV)
  329.                         DTDROW=-WP/(RHOP*CV)
  330.                         DTDRET=1.D0/(RHOP*CV)
  331.                         IF (.NOT.LMUR) THEN
  332.                            CTGRT.POIPR2(IPP,IEL2)=NPPRIM
  333.                            RETRHO.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRHO
  334.                            RETROU.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROU
  335.                            RETROV.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROV
  336.                            RETROW.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDROW
  337.                            RETRET.VALMA2(IPD,IPP,IEL2)=FACTOR*DTDRET
  338.                         ELSE
  339.                            CTGRT.POIPR1(IPP,IEL1)=NPPRIM
  340.                            RETRHO.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRHO
  341.                            RETROU.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROU
  342.                            RETROV.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROV
  343.                            RETROW.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDROW
  344.                            RETRET.VALMA1(IPD,IPP,IEL1)=FACTOR*DTDRET
  345.                         ENDIF
  346.                      ENDIF
  347.  124              CONTINUE
  348.  122           CONTINUE
  349.                IF (.NOT.LMUR) THEN
  350.                   IEL2=IEL2+1
  351.                ELSE
  352.                   IEL1=IEL1+1
  353.                ENDIF
  354.             ENDIF
  355.  12      CONTINUE
  356.          NPP1=NPTEL-1
  357.          NPP2=NPTEL-1
  358.          NEL1=IEL1-1
  359.          NEL2=IEL2-1
  360.          SEGADJ RETRHO
  361.          SEGADJ RETROU
  362.          SEGADJ RETROV
  363.          SEGADJ RETROW
  364.          SEGADJ RETRET
  365.          SEGADJ CTGRT
  366.          CTGRT.LMATSI(**)=RETRHO
  367.          CTGRT.LMATSI(**)=RETROU
  368.          CTGRT.LMATSI(**)=RETROV
  369.          CTGRT.LMATSI(**)=RETROW
  370.          CTGRT.LMATSI(**)=RETRET
  371. * On accumule les matrices résultantes dans IJACO
  372.          CALL AJMTK(CTGRT,IJACO,IMPR,IRET)
  373.          IF (IRET.NE.0) GOTO 9999
  374.          SEGSUP RETRHO
  375.          SEGSUP RETROU
  376.          SEGSUP RETROV
  377.          SEGSUP RETROW
  378.          SEGSUP RETRET
  379.          SEGSUP CTGRT
  380. *
  381.          SEGDES MCOGRT
  382.          SEGDES KCDTDZ
  383.          SEGDES KCDTDY
  384.          SEGDES KCDTDX
  385.  1    CONTINUE
  386.       SEGDES ICOGRT
  387.       SEGDES MPTEMC
  388.       SEGDES MPVITC
  389.       SEGDES MPROC
  390.       SEGDES MPVOLU
  391.       SEGDES MPNORM
  392.       SEGDES MPSURF
  393.       SEGDES MELEFL
  394.       SEGDES KRFACE
  395.       SEGDES KRCENT
  396.       SEGDES NOMINC
  397.       IF (LCLIMQ) THEN
  398.          SEGDES KRQIMP
  399.       ENDIF
  400.       IF (LCLIMT) THEN
  401.          SEGDES KRTIMP
  402.       ENDIF
  403. *
  404. * Normal termination
  405. *
  406.       IRET=0
  407.       RETURN
  408. *
  409. * Format handling
  410. *
  411. *
  412. * Error handling
  413. *
  414.  9999 CONTINUE
  415.       IRET=1
  416.       WRITE(IOIMP,*) 'An error was detected in subroutine ylap2b'
  417.       RETURN
  418. *
  419. * End of subroutine YLAP2B
  420. *
  421.       END
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  

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