Télécharger la plaquette Cast3M


fr en

Télécharger zlap13.eso

Retour à la liste

Numérotation des lignes : 

zlap13

  1. C ZLAP13    SOURCE    CB215821  20/11/25    13:45:02     10792          
  2.       SUBROUTINE ZLAP13(PROPHY,IROC,ITEMC,
  3.      &     ITIMP,IRIMP,
  4.      &     MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
  5.      &     ICHFLU,DT)
  6. C***********************************************************************
  7. C NOM         : ZLAP13
  8. C DESCRIPTION : Calcul des flux diffusifs Navier-Stokes3D multi-espèces.
  9. C
  10. C
  11. C LANGAGE     : ESOPE
  12. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  13. C               mél : gounand@semt2.smts.cea.fr
  14. C***********************************************************************
  15. C APPELES (UTIL)   : KRIPAD : MELEME -> (num. globale->locale)
  16. C                    LICHT  : Lecture des pointeurs (maillages, valeurs)
  17. C                             d'un objet de type MCHPOI.
  18. C                    ERREUR : Gestion des erreurs par GIBI.
  19. C APPELE PAR       : ZLAP11 : Chapeau de l'opérateur Gibiane 'LAPN'
  20. C                             option 'VF'.
  21. C***********************************************************************
  22. C ENTREES            : PROPHY (type PROPHY) : propriétés des espèces
  23. C                      PROPH2 (type PROPH2) : précond. de PROPHY
  24. C                      IROC   (type MCHPOI) : masse volumique par
  25. C                                             élément.
  26. C                      ITEMC  (type MCHPOI) : température par
  27. C                                             élément.
  28. C                      ITIMP  (type MCHPOI) : CL de
  29. C                        Dirichlet sur la température.
  30. C                      IRIMP  (type MCHPOI) : CL de
  31. C                        Dirichlet sur la densité.
  32. C                      MELEMC (type MELEME) : maillage des centres des
  33. C                                             éléments.
  34. C                      MELEMF (type MELEME) : maillage des faces des
  35. C                                             éléments.
  36. C                      MELEFL (type MELEME) : connectivités face-(centre
  37. C                                             gauche, centre droit).
  38. C                      ISURF (type MCHPOI) : surface des faces.
  39. C                      INORM (type MCHPOI) : normale aux faces.
  40. C                      IDIAM (type MCHPOI) : diamètre des éléments.
  41. C ENTREES/SORTIES    : -
  42. C SORTIES            : ICHFLU (type MCHPOI) : flux diffusif aux
  43. C                                             interfaces.
  44. C                      DT     (type REAL*8) : pas de temps de stabilité
  45. C                                             (Fourier)
  46. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  47. C***********************************************************************
  48. C VERSION    : v1, 05/03/2002, version initiale
  49. C HISTORIQUE : v1, 05/03/2002, création
  50. C HISTORIQUE :
  51. C HISTORIQUE :
  52. C***********************************************************************
  53. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  54. C en cas de modification de ce sous-programme afin de faciliter
  55. C la maintenance !
  56. C***********************************************************************
  57.       IMPLICIT INTEGER(I-N)
  58.  
  59. -INC PPARAM
  60. -INC CCOPTIO
  61. -INC CCREEL
  62. -INC SMCHPOI
  63. -INC SMELEME
  64. -INC SMCOORD
  65. -INC SMLENTI
  66. -INC SMLMOTS
  67. C
  68.       POINTEUR IROC.MCHPOI,ITEMC.MCHPOI,ICDIFF.MCHPOI,IGRYKF.MCHPOI
  69.       POINTEUR MPROC.MPOVAL,MPTEMC.MPOVAL,
  70.      &     MPTIMP.MPOVAL,MPRIMP.MPOVAL,MPCDIF.MPOVAL,MPGRYK.MPOVAL,
  71.      &     MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL,
  72.      &     MPFLUX.MPOVAL,MPCDIN.MPOVAL
  73. C
  74.       POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME
  75.       POINTEUR MLCENT.MLENTI,MLETIM.MLENTI,MLERIM.MLENTI
  76. C
  77.       INTEGER NESP
  78.       SEGMENT PROPHY
  79.          CHARACTER*4 NOMESP(NESP+1)
  80.          REAL*8 CV(NESP+1)
  81.          REAL*8 R(NESP+1)
  82.          REAL*8 H0K(NESP+1)
  83.          POINTEUR CDIFF(NESP+1).MCHPOI
  84.          POINTEUR YK(NESP+1).MCHPOI
  85.          POINTEUR GRADYK(NESP+1).MCHPOI
  86.          POINTEUR CGRYK(NESP+1).MCHELM
  87.          POINTEUR CLYK(NESP+1).MCHPOI
  88.       ENDSEGMENT
  89.       SEGMENT LIPOVA
  90.       POINTEUR MPGRAD(0).MPOVAL
  91.       ENDSEGMENT
  92.       SEGMENT LIPOV2
  93.       POINTEUR MPDIFF(0).MPOVAL
  94.       ENDSEGMENT
  95. C
  96.       INTEGER ITIMP,IRIMP
  97.      &       ,ISURF,INORM,IDIAM,ICHFLU
  98.      &       ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED
  99.      &       ,NLCEG,NLCED,NLFTI,NLFRI
  100.      &       , ICOORX, IGEOM
  101.  
  102.       REAL*8 DT, UNSDT
  103.      &       ,XG,YG,ZG,XFMXG,YFMYG,ZFMZG,DRG
  104.      &       ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH
  105.      &       ,XF,YF,ZF
  106.      &       ,CNX,CNY,CNZ,ORIENT,DIAM,DIAM2,CELL,SURF
  107.       INTEGER IESP
  108.       REAL*8 CPK,DK,DYKDX,DYKDY,DYKDZ,HKF,CPN,DN,HNF
  109.       REAL*8 RHOD,RHOF,RHOG
  110.       REAL*8 TEMD,TEMF,TEMG
  111. C
  112.       CHARACTER*8 TYPE
  113. C
  114. C**** Initialisation de 1/DT
  115. C
  116.       UNSDT = 0.0D0
  117. C
  118. C**** KRIPAD pour la correspondance global/local de centre
  119. C
  120.       CALL KRIPAD(MELEMC,MLCENT)
  121. C
  122. C     EN KRIPAD
  123. C     SEGACT MELEMC
  124. C     SEGACT MLCENT
  125. C
  126.       CALL LICHT(IROC,MPROC,TYPE,IGEOM)
  127.       CALL LICHT(ITEMC,MPTEMC,TYPE,IGEOM)
  128.       CALL LICHT(ISURF,MPSURF,TYPE,IGEOM)
  129.       CALL LICHT(INORM,MPNORM,TYPE,IGEOM)
  130.       CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM)
  131.       CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM)
  132. C
  133. C     EN LICHT
  134. C     SEGACT*MOD MPROC
  135. C     SEGACT*MOD MPTEMC
  136. C     SEGACT*MOD MPSURF
  137. C     SEGACT*MOD MPNORM
  138. C     SEGACT*MOD MPDIAM
  139. C     SEGACT*MOD MPFLUX
  140. C
  141.       SEGINI LIPOVA
  142.       SEGINI LIPOV2
  143.       SEGACT PROPHY
  144.       NESP=PROPHY.CV(/1)-1
  145.       DO IESP=1,NESP+1
  146.          ICDIFF=PROPHY.CDIFF(IESP)
  147. * In licht SEGACT LIPOV2.MPDIFF(*)
  148.          CALL LICHT(ICDIFF,MPCDIF,TYPE,IGEOM)
  149.          LIPOV2.MPDIFF(**)=MPCDIF
  150.          IGRYKF=PROPHY.GRADYK(IESP)
  151. * In licht SEGACT LIPOVA.MPGRAD(*)
  152.          CALL LICHT(IGRYKF,MPGRYK,TYPE,IGEOM)
  153.          LIPOVA.MPGRAD(**)=MPGRYK
  154.       ENDDO
  155.  
  156.       IF(ITIMP .NE. 0)THEN
  157.          CALL LICHT(ITIMP,MPTIMP,TYPE,IGEOM)
  158. C        SEGACT*MOD MPVIMP
  159.          CALL KRIPAD(IGEOM,MLETIM)
  160. C        SEGACT IGEOM
  161. C        SEGACT MLETIM
  162.          MELEME = IGEOM
  163.          SEGDES MELEME
  164.       ENDIF
  165.       IF(IRIMP .NE. 0)THEN
  166.          CALL LICHT(IRIMP,MPRIMP,TYPE,IGEOM)
  167. C        SEGACT*MOD MPVIMP
  168.          CALL KRIPAD(IGEOM,MLERIM)
  169. C        SEGACT IGEOM
  170. C        SEGACT MLERIM
  171.          MELEME = IGEOM
  172.          SEGDES MELEME
  173.       ENDIF
  174.  
  175. C
  176.       SEGACT MELEFL
  177.       SEGACT MELEMF
  178.       NFAC = MELEMF.NUM(/2)
  179. C
  180. C**** Boucle sur les faces
  181. C
  182.       DO NLCF = 1, NFAC, 1
  183. C
  184. C******* NLCF  = numero local du centre de facel
  185. C        NGCF  = numero global du centre de facel
  186. C        NLCF1 = numero local du centre de face
  187. C        NGCEG = numero global du centre ELT "gauche"
  188. C        NLCEG = numero local du centre ELT "gauche"
  189. C        NGCED = numero global du centre ELT "droite"
  190. C        NLCED = numero local du centre ELT "droite"
  191. C
  192.          NGCF = MELEMF.NUM(1,NLCF)
  193.          NGCF1 = MELEFL.NUM(2,NLCF)
  194.          IF(NGCF .NE. NGCF1)THEN
  195.             MOTERR(1:40)= 'FACEL et FACE = ?                       '
  196.             CALL ERREUR(5)
  197.             GOTO 9999
  198.          ENDIF
  199. C
  200.          NGCEG = MELEFL.NUM(1,NLCF)
  201.          NGCED = MELEFL.NUM(3,NLCF)
  202.          NLCEG = MLCENT.LECT(NGCEG)
  203.          NLCED = MLCENT.LECT(NGCED)
  204. C
  205. C******* On controlle si sur NGCF on impose de CL
  206. C
  207. C        NLFTI  = numero local du centre de face sul le maillage des
  208. C                 "températures" "imposées"
  209. C
  210.          IF(ITIMP .NE. 0)THEN
  211.             NLFTI = MLETIM.LECT(NGCF)
  212.          ELSE
  213.             NLFTI = 0
  214.          ENDIF
  215. C
  216.          IF(IRIMP .NE. 0)THEN
  217.             NLFRI = MLERIM.LECT(NGCF)
  218.          ELSE
  219.             NLFRI = 0
  220.          ENDIF
  221. C
  222.          IF(NGCEG .NE. NGCED)THEN
  223. C
  224. C********** Parametres geometriques
  225. C
  226.             ICOORX = ((IDIM + 1) * (NGCF - 1))+1
  227.             XF = MCOORD.XCOOR(ICOORX)
  228.             YF = MCOORD.XCOOR(ICOORX+1)
  229.             ZF = MCOORD.XCOOR(ICOORX+2)
  230. C
  231.             ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  232.             XG = MCOORD.XCOOR(ICOORX)
  233.             YG = MCOORD.XCOOR(ICOORX+1)
  234.             ZG = MCOORD.XCOOR(ICOORX+2)
  235.             XFMXG = XF - XG
  236.             YFMYG = YF - YG
  237.             ZFMZG = ZF - ZG
  238.             DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
  239. C
  240.             ICOORX = ((IDIM + 1) * (NGCED - 1))+1
  241.             XD = MCOORD.XCOOR(ICOORX)
  242.             YD = MCOORD.XCOOR(ICOORX+1)
  243.             ZD = MCOORD.XCOOR(ICOORX+2)
  244.             XFMXD = XF - XD
  245.             YFMYD = YF - YD
  246.             ZFMZD = ZF - ZD
  247.             DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
  248. C
  249. C********** F=G -> DRG = 0 -> ALPHA = 0
  250. C
  251.             ALPHA=DRG/(DRG+DRD)
  252.             UMALPH= 1.0D0 - ALPHA
  253. C
  254. C********** Les valeurs à l'interface
  255. C
  256. C           DRG=0 -> F=G
  257. C
  258. C
  259. C********** Température
  260. C
  261.             IF( NLFTI .GT. 0) THEN
  262.                TEMF = MPTIMP.VPOCHA(NLFTI,1)
  263.             ELSE
  264. C************* Gauche
  265.                TEMG = MPTEMC.VPOCHA(NLCEG,1)
  266. C************* Droite
  267.                TEMD = MPTEMC.VPOCHA(NLCED,1)
  268. C************* Face
  269.                TEMF = UMALPH * TEMG + ALPHA * TEMD
  270.             ENDIF
  271. C
  272. C********** Rho
  273. C
  274.             IF( NLFRI .GT. 0) THEN
  275.                RHOF = MPRIMP.VPOCHA(NLFRI,1)
  276.             ELSE
  277. C************* Gauche
  278.                RHOG = MPROC.VPOCHA(NLCEG,1)
  279. C************* Droite
  280.                RHOD = MPROC.VPOCHA(NLCED,1)
  281. C************* Face
  282.                RHOF = UMALPH * RHOG + ALPHA * RHOD
  283.             ENDIF
  284.          ELSE
  285. C
  286. C********** MURS
  287. C
  288. C           Etat a gauche = Etat droite
  289. C
  290.             ALPHA=0.0D0
  291.             UMALPH=1.0D0
  292. C
  293. C********** Parametres geometriques
  294. C
  295.             ICOORX = ((IDIM + 1) * (NGCF - 1))+1
  296.             XF = MCOORD.XCOOR(ICOORX)
  297.             YF = MCOORD.XCOOR(ICOORX+1)
  298.             ZF = MCOORD.XCOOR(ICOORX+2)
  299. C
  300.             ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  301.             XG = MCOORD.XCOOR(ICOORX)
  302.             YG = MCOORD.XCOOR(ICOORX+1)
  303.             ZG = MCOORD.XCOOR(ICOORX+2)
  304.             XFMXG = XF - XG
  305.             YFMYG = YF - YG
  306.             ZFMZG = ZF - ZG
  307. C
  308. C********** Température
  309. C
  310.             IF( NLFTI .GT. 0) THEN
  311.                TEMF = MPTIMP.VPOCHA(NLFTI,1)
  312.             ELSE
  313.                TEMF = MPTEMC.VPOCHA(NLCEG,1)
  314.             ENDIF
  315. C
  316. C********** Rho
  317. C
  318.             IF( NLFRI .GT. 0) THEN
  319.                RHOF = MPRIMP.VPOCHA(NLFRI,1)
  320.             ELSE
  321.                RHOF = MPROC.VPOCHA(NLCEG,1)
  322.             ENDIF
  323.          ENDIF
  324. C
  325. C******* On calcule le sign du pruduit scalare
  326. C        (Normales de Castem) * (vecteur "gauche" -> "centre")
  327. C
  328.          CNX = MPNORM.VPOCHA(NLCF,1)
  329.          CNY = MPNORM.VPOCHA(NLCF,2)
  330.          CNZ = MPNORM.VPOCHA(NLCF,3)
  331.          ORIENT = CNX * XFMXG + CNY * YFMYG + CNZ * ZFMZG
  332.          ORIENT = SIGN(1.0D0,ORIENT)
  333.          IF(ORIENT .NE. 1.0D0)THEN
  334.             MOTERR(1:40)=
  335.      &           'LAPN , subroutine zlap13.eso.           '
  336.             WRITE(IOIMP,*) MOTERR(1:40)
  337.             MOTERR(1:40)=
  338.      &           'Orientation normales.                   '
  339.             WRITE(IOIMP,*) MOTERR(1:40)
  340.             CALL ERREUR(5)
  341.             GOTO 9999
  342.          ENDIF
  343. C
  344. C******* Les flux aux interfaces
  345. C
  346.          SURF = MPSURF.VPOCHA(NLCF,1)
  347.          DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
  348.      &        ALPHA*MPDIAM.VPOCHA(NLCED,1)
  349.          DIAM2=DIAM*DIAM
  350.          MPCDIN=LIPOV2.MPDIFF(NESP+1)
  351.          CPN=PROPHY.CV(NESP+1)+PROPHY.R(NESP+1)
  352.          HNF=(CPN*TEMF)+PROPHY.H0K(NESP+1)
  353.          DO IESP=1,NESP
  354.             MPCDIF=LIPOV2.MPDIFF(IESP)
  355.             DK=UMALPH*MPCDIF.VPOCHA(NLCEG,1) +
  356.      $           ALPHA*MPCDIF.VPOCHA(NLCED,1)
  357.             DN=UMALPH*MPCDIN.VPOCHA(NLCEG,1) +
  358.      $           ALPHA*MPCDIN.VPOCHA(NLCED,1)
  359.             MPGRYK=LIPOVA.MPGRAD(IESP)
  360.             DYKDX=MPGRYK.VPOCHA(NLCF,1)
  361.             DYKDY=MPGRYK.VPOCHA(NLCF,2)
  362.             DYKDZ=MPGRYK.VPOCHA(NLCF,3)
  363.             CPK=PROPHY.CV(IESP)+PROPHY.R(IESP)
  364.             HKF=(CPK*TEMF)+PROPHY.H0K(IESP)
  365. * Contribution pour l'espèce IESP
  366.             MPFLUX.VPOCHA(NLCF,IDIM+1+IESP)=
  367.      $           RHOF*DK*((DYKDX*CNX)+(DYKDY*CNY)+(DYKDZ*CNZ))
  368.      $           *SURF*(-1.D0)
  369.  
  370. * Contribution pour l'énergie totale
  371.             MPFLUX.VPOCHA(NLCF,IDIM+1)=
  372.      $           MPFLUX.VPOCHA(NLCF,IDIM+1)+
  373.      $           (RHOF*((HKF*DK)-(HNF*DN))
  374.      $           *((DYKDX*CNX)+(DYKDY*CNY)+(DYKDZ*CNZ))*SURF*(-1.D0))
  375. * Le pas de temps
  376.             CELL=(4.0D0*DK)/DIAM2
  377.             IF(CELL .GT. UNSDT)THEN
  378.                UNSDT=CELL
  379.             ENDIF
  380.          ENDDO
  381.       ENDDO
  382. C
  383. C
  384.       DT = 1.0D0 / (UNSDT + XPETIT)
  385. C
  386.       SEGDES MELEFL
  387.       SEGDES MELEMF
  388.       SEGDES MELEMC
  389.       SEGDES MPSURF
  390.       SEGDES MPNORM
  391.       SEGDES MPDIAM
  392.       SEGSUP MLCENT
  393. C
  394.       SEGDES MPROC
  395.       SEGDES MPTEMC
  396.       SEGDES PROPHY
  397.       SEGDES LIPOVA.MPGRAD(*)
  398.       SEGSUP LIPOVA
  399.       SEGDES LIPOV2.MPDIFF(*)
  400.       SEGSUP LIPOV2
  401.       SEGDES MPFLUX
  402. C
  403.       IF(ITIMP .NE. 0) THEN
  404.          SEGDES MPTIMP
  405.          SEGSUP MLETIM
  406.       ENDIF
  407. C
  408.       IF(IRIMP .NE. 0) THEN
  409.          SEGDES MPRIMP
  410.          SEGSUP MLERIM
  411.       ENDIF
  412. C
  413.  9999 CONTINUE
  414.       RETURN
  415.       END
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  

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