Télécharger la plaquette Cast3M


fr en

Télécharger cli281.eso

Retour à la liste

Numérotation des lignes : 

cli281

  1. C CLI281    SOURCE    CB215821  20/11/25    13:20:56     10792          
  2.       SUBROUTINE CLI281(NSP,MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,
  3.      &     ICHPSU,LRECP,LRECV,IROC,IVITC,IPC,IYN,ICHLIM,ICHRES,ICHRLI)
  4. C************************************************************************
  5. C
  6. C PROJET            :  CASTEM 2000
  7. C
  8. C NOM               :  CLI281
  9. C
  10. C DESCRIPTION       :  Subroutine appellée par CLIM22
  11. C                       calcul de RESIDU et CLIM at the board
  12. C                       OPTION: 'INOU' 2D
  13. C
  14. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  15. C
  16. C AUTEUR            :  S.Kudriakov, DEN/DM2S/SFME/LTMF
  17. C
  18. C************************************************************************
  19. C
  20. C APPELES (Calcul)  :
  21. C
  22. C************************************************************************
  23. C
  24. C HISTORIQUE (Anomalies et modifications éventuelles)
  25. C
  26. C HISTORIQUE :
  27. C
  28. C************************************************************************
  29. C
  30.       IMPLICIT INTEGER(I-N)
  31.       IMPLICIT REAL*8(A-H,O-Z)
  32.  
  33. -INC PPARAM
  34. -INC CCOPTIO
  35. -INC SMLMOTS
  36. -INC SMELEME
  37.       POINTEUR MELEFC.MELEME
  38. -INC SMLENTI
  39.       POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI
  40. -INC SMCHPOI
  41.       POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL,
  42.      &     MPVC.MPOVAL, MPPC.MPOVAL, MPYN.MPOVAL, MPLIM.MPOVAL,
  43.      &     MPRES.MPOVAL,  MPRLI.MPOVAL
  44.  
  45.       INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
  46.      &     ,IYN,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC
  47.      &     ,NGF,NGC,NLF,NLC,NLCB,LRECP,LRECV,I,NSP,NESP
  48.       REAL*8 VOLU,SURF,RC,PC,UXC,UYC,UZC,GAMC,CNX,CNY,CNZ,CTX,CTY,CTZ
  49.      &     ,CT2X,CT2Y,CT2Z,RF,PF,UXF,UYF,UZF,TOP,BOT
  50.      &     ,UNC,UNF,UTF,UT2F,SF,GAMF,ECIN,PSRF,HTF,GM1
  51.      &     ,CELLT,UT2C,UTC,TEMP
  52.       CHARACTER*(8) TYPE
  53. C------------------------------------------------------------
  54. -INC SMLREEL
  55.        POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  56. C-------------------------------------------------------
  57. C**********  Les CP's and CV's   ***********************
  58. C-------------------------------------------------------
  59.       SEGMENT GCONST
  60.          REAL*8 GC(NSP)
  61.       ENDSEGMENT
  62.       POINTEUR CP.GCONST, CV.GCONST
  63. C-------------------------------------------------------------
  64. C******* Les fractionines massiques **************************
  65. C-------------------------------------------------------------
  66.       SEGMENT FRAMAS
  67.          REAL*8 YET(NSP)
  68.       ENDSEGMENT
  69.       POINTEUR YC.FRAMAS, YF.FRAMAS
  70. C-------------------------------------------------------------
  71. C**********  Segments for the flux-vector  *******************
  72. C-------------------------------------------------------------
  73.        SEGMENT FUNEL
  74.          REAL*8 FU(4+NSP)
  75.        ENDSEGMENT
  76.        POINTEUR flux2D.funel, flux3D.funel
  77.        SEGINI FLUX2D
  78.        SEGINI FLUX3D
  79. C------------------------------------------------------
  80. C**** KRIPAD pour la correspondance global/local
  81. C------------------------------------------------------
  82.       CALL KRIPAD(MELEMC,MLEMC)
  83.       CALL KRIPAD(MELECB,MLEMCB)
  84.       CALL KRIPAD(MELEMF,MLEMF)
  85. C------------------------------------------------------
  86. C**** CHPOINTs de la table DOMAINE
  87. C------------------------------------------------------
  88.       CALL LICHT(INORM,MPNORM,TYPE,ICEL)
  89.       CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL)
  90.       CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL)
  91. C------------------------------------------------------
  92. C**** CHPOINTs des variables
  93. C------------------------------------------------------
  94.       CALL LICHT(IROC,MPRC,TYPE,ICEL)
  95.       CALL LICHT(IVITC,MPVC,TYPE,ICEL)
  96.       CALL LICHT(IPC,MPPC,TYPE,ICEL)
  97.       CALL LICHT(IYN,MPYN,TYPE,ICEL)
  98.       CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL)
  99.       CALL LICHT(ICHRES,MPRES,TYPE,ICEL)
  100.       CALL LICHT(ICHRLI,MPRLI,TYPE,ICEL)
  101. C---------------------------------------------------------
  102. C**** Boucle sur le face pour le calcul des invariants de
  103. C     Riemann et du flux
  104. C---------------------------------------------------------
  105.       SEGACT MELEFC
  106.       NFAC=MELEFC.NUM(/2)
  107. c      write(*,*) 'hell', NFAC
  108.       UZC=0.0D0
  109.       CNZ=0.0D0
  110.       CTZ=0.0D0
  111.       CT2X=0.0D0
  112.       CT2Y=0.0D0
  113.       CT2Z=0.0D0
  114.       DO 1 IFAC=1,NFAC,1
  115.          NGF=MELEFC.NUM(1,IFAC)
  116.          NGC=MELEFC.NUM(2,IFAC)
  117.          NLF=MLEMF.LECT(NGF)
  118.          NLC=MLEMC.LECT(NGC)
  119.          NLCB=MLEMCB.LECT(NGF)
  120.          VOLU=MPVOL.VPOCHA(NLC,1)
  121.          SURF=MPSURF.VPOCHA(NLF,1)
  122. C----------------------------------------------
  123. C        In CASTEM les normales sont sortantes
  124. C----------------------------------------------
  125.          CNX=MPNORM.VPOCHA(NLF,1)
  126.          CNY=MPNORM.VPOCHA(NLF,2)
  127.          IF(IDIM.EQ.2)THEN
  128.             CTX=-1.0D0*CNY
  129.             CTY=CNX
  130.          ELSE
  131.             CNZ=MPNORM.VPOCHA(NLF,3)
  132.             CTX=MPNORM.VPOCHA(NLF,4)
  133.             CTY=MPNORM.VPOCHA(NLF,5)
  134.             CTZ=MPNORM.VPOCHA(NLF,6)
  135.             CT2X=MPNORM.VPOCHA(NLF,7)
  136.             CT2Y=MPNORM.VPOCHA(NLF,8)
  137.             CT2Z=MPNORM.VPOCHA(NLF,9)
  138.          ENDIF
  139. C----------------------------------------
  140.        SEGINI CP, CV
  141.        MLRECP = LRECP
  142.        MLRECV = LRECV
  143.        SEGACT MLRECP, MLRECV
  144.        DO 10 I=1,(NSP-1)
  145.            CP.GC(I)=MLRECP.PROG(I)
  146.            CV.GC(I)=MLRECV.PROG(I)
  147.  10    CONTINUE
  148.        CP.GC(NSP)=MLRECP.PROG(NSP)
  149.        CV.GC(NSP)=MLRECV.PROG(NSP)
  150. C----------------------------
  151. C        Variables au centre
  152. C----------------------------
  153.          RC=MPRC.VPOCHA(NLC,1)
  154.          PC=MPPC.VPOCHA(NLC,1)
  155.          UXC=MPVC.VPOCHA(NLC,1)
  156.          UYC=MPVC.VPOCHA(NLC,2)
  157.          IF(IDIM.EQ.3)UZC=MPVC.VPOCHA(NLC,3)
  158.          SEGINI YC
  159.          SEGACT MPYN
  160.          DO 100 I=1,(NSP-1)
  161.            YC.YET(I)=MPYN.VPOCHA(NLC,I)
  162.  100     CONTINUE
  163. c-------------------------------------------------------------
  164. c  Computing GAMMA at the cell-center
  165. c-------------------------------------------------------------
  166.          top=0.0D0
  167.          bot=0.0D0
  168.          do 102  i=1,(nsp-1)
  169.           top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp))
  170.           bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp))
  171.  102     continue
  172.          top=cp.gc(nsp)+top
  173.          bot=cv.gc(nsp)+bot
  174.          GAMC=top/bot
  175. C-------------------------------------------------------------
  176. C   We decide whether the boundary conditions are 'inlet' or
  177. C   'outlet' by looking at the UNC (UNC > 0.0 -- 'outlet')
  178. C-------------------------------------------------------------
  179.         UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
  180. c        write(*,*) 'before', UNC, NLC,NGC
  181.         IF(UNC .GE. 0.0D0) THEN
  182. C--------------------------------------
  183. C 'OUTP' -- outlet with given pressure
  184. C--------------------------------------
  185. C-----------------------------------------
  186. C        Variables à la face
  187. C-----------------------------------------
  188.          PF=MPLIM.VPOCHA(NLCB,NSP+2)
  189. C---------------------------------------
  190. C******* On calcule UN, UT, UT2
  191. C---------------------------------------
  192.          UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
  193.          UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z)
  194. C-----------------------------------------------
  195. C******* Densite, vitesse, pression sur le bord
  196. C-----------------------------------------------
  197.          MPRLI.VPOCHA(NLCB,1)=RC
  198.          MPRLI.VPOCHA(NLCB,2)=UXC
  199.          MPRLI.VPOCHA(NLCB,3)=UYC
  200.          IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZC
  201.          MPRLI.VPOCHA(NLCB,IDIM+2)=PF
  202.          do 104  i=1,(nsp-1)
  203.            MPRLI.VPOCHA(NLCB,IDIM+2+I)=YC.YET(I)
  204.  104     continue
  205. C---------------------------------------------------
  206. C******* Probleme de Riemann entre l'etat gauche
  207. C        RC,UNC,UTC,UT2C,PC et l'etat droite
  208. C        RC,UNC,UTC,UT2C,PF
  209. C        On utilise AUSM+
  210. C        Flux dans le repaire normale
  211. C---------------------------------------------------
  212.          NESP=NSP-1
  213.          IF(IDIM.EQ.2)THEN
  214.             CALL FAUSMP(NESP,
  215.      &           GAMC,RC,PC,UNC,UTC,
  216.      &           GAMC,RC,PF,UNC,UTC,
  217.      &           YC.YET,YC.YET,
  218.      &           FLUX2D.FU,
  219.      &           CELLT)
  220. C-------------------------------------------------------
  221. C******* Residuum (son SPG a le meme ordre que MELEFC)
  222. C-------------------------------------------------------
  223.             MPRES.VPOCHA(IFAC,1)=-1*FLUX2D.FU(1)*SURF/VOLU
  224.              MPRES.VPOCHA(IFAC,2)=-1*((FLUX2D.FU(2)*CNX)+
  225.      &          (FLUX2D.FU(3)*CTX))*SURF/VOLU
  226.              MPRES.VPOCHA(IFAC,3)=-1*((FLUX2D.FU(2)*CNY)+
  227.      &          (FLUX2D.FU(3)*CTY))*SURF/VOLU
  228.              MPRES.VPOCHA(IFAC,4)=-1*FLUX2D.FU(4)*SURF/VOLU
  229.              do 105  i=1,(nsp-1)
  230.               MPRES.VPOCHA(IFAC,4+I)=-1*FLUX2D.FU(4+I)*SURF/VOLU
  231.  105       continue
  232.          ELSE
  233.             CALL FAUSM3(NESP,
  234.      &           GAMC,RC,PC,UNC,UTC,UT2C,
  235.      &           GAMC,RC,PF,UNC,UTC,UT2C,
  236.      &           YC.YET,YC.YET,
  237.      &           FLUX3D.FU,
  238.      &           CELLT)
  239. C------------------------------------------------------
  240. C******* Residuum (son SPG a le meme ordre que MELEFC)
  241. C------------------------------------------------------
  242.             MPRES.VPOCHA(IFAC,1)=-1*FLUX3D.FU(1)*SURF/VOLU
  243.             MPRES.VPOCHA(IFAC,2)=-1*((FLUX3D.FU(2)*CNX)+
  244.      &       (FLUX3D.FU(3)*CTX)+(FLUX3D.FU(4)*CT2X))*SURF/VOLU
  245.             MPRES.VPOCHA(IFAC,3)=-1*((FLUX3D.FU(2)*CNY)+
  246.      &       (FLUX3D.FU(3)*CTY)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  247.             MPRES.VPOCHA(IFAC,4)=-1*((FLUX3D.FU(2)*CNZ)+
  248.      &       (FLUX3D.FU(3)*CTZ)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  249.             MPRES.VPOCHA(IFAC,5)=-1*FLUX3D.FU(5)*SURF/VOLU
  250.             do 106  i=1,(nsp-1)
  251.               MPRES.VPOCHA(IFAC,5+I)=-1*FLUX3D.FU(5+I)*SURF/VOLU
  252.  106        continue
  253.          ENDIF
  254.          TEMP=-1000.0
  255.         ELSE
  256. c        write(*,*) 'insuuuuuuuu', NLCB
  257. C-------------------------------------------------------
  258. C  'INSU' -- 'inlet' b.c. with given Ht and S (entropy)
  259. C-------------------------------------------------------
  260.          UZF=0.0D0
  261.          UT2F=0.0D0
  262. C----------------------------
  263. C        Variables à la face
  264. C----------------------------
  265.          HTF=MPLIM.VPOCHA(NLCB,1)
  266.          SF=MPLIM.VPOCHA(NLCB,2)
  267.          SEGINI YF
  268.          TEMP=1000.0
  269.          DO 101 I=1,(NSP-1)
  270.            YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I)
  271.  101     CONTINUE
  272.          UTF=0.0D0
  273. c-------------------------------------------------------------
  274. c  Computing GAMMA at the face-center
  275. c-------------------------------------------------------------
  276.          top=0.0D0
  277.          bot=0.0D0
  278.          do 103  i=1,(nsp-1)
  279.           top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp))
  280.           bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp))
  281.  103     continue
  282.          top=cp.gc(nsp)+top
  283.          bot=cv.gc(nsp)+bot
  284.          GAMF=top/bot
  285.          GM1=GAMF-1.0D0
  286. C-----------------------------------------
  287. C changing direction of the normal vector
  288. C-----------------------------------------
  289.          CNX=-1.0D0*CNX
  290.          CNY=-1.0D0*CNY
  291.          IF(IDIM.EQ.2)THEN
  292.             CTX=-1.0D0*CNY
  293.             CTY=CNX
  294.          ELSE
  295.             CNZ=-1.0D0*CNZ
  296.             CTX=-1.0D0*CTX
  297.             CTY=-1.0D0*CTY
  298.             CTZ=-1.0D0*CTZ
  299.             CT2X=-1.0D0*CT2X
  300.             CT2Y=-1.0D0*CT2Y
  301.             CT2Z=-1.0D0*CT2Z
  302.          ENDIF
  303. C---------------------------------------
  304. C******* On calcule UN, UT, UT2, ASON, S
  305. C---------------------------------------
  306.          UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
  307.          UNF=UNC
  308.          UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
  309. C----------------------------------
  310.          UXF=UNF*CNX+UTF*CTX+UT2F*CT2X
  311.          UYF=UNF*CNY+UTF*CTY+UT2F*CT2Y
  312.          UZF=UNF*CNZ+UTF*CTZ+UT2F*CT2Z
  313. C----------------------------------
  314.          ECIN=0.5D0*((UXF*UXF)+(UYF*UYF)+(UZF*UZF))
  315.          PSRF=(GM1/GAMF)*(HTF-ECIN)
  316.          RF=PSRF/SF
  317.          RF=RF**(1.0D0/GM1)
  318.          PF=SF*(RF**GAMF)
  319. C-----------------------------------------------
  320. C******* Densite, vitesse, pression sur le bord
  321. C-----------------------------------------------
  322.          MPRLI.VPOCHA(NLCB,1)=RF
  323.          MPRLI.VPOCHA(NLCB,2)=UXF
  324.          MPRLI.VPOCHA(NLCB,3)=UYF
  325.          IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZF
  326.          MPRLI.VPOCHA(NLCB,IDIM+2)=PF
  327.          do 1040  i=1,(nsp-1)
  328.            MPRLI.VPOCHA(NLCB,IDIM+2+I)=YF.YET(I)
  329.  1040    continue
  330. C---------------------------------------------------
  331. C******* Probleme de Riemann entre l'etat gauche
  332. C        RF,UNC,UTF,UT2F,PF et l'etat droite
  333. C        RC,UNC,UTC,UT2C,PC
  334. C        On utilise AUSM+
  335. C        Flux dans le repaire normale
  336. C---------------------------------------------------
  337.          NESP=NSP-1
  338.          IF(IDIM.EQ.2)THEN
  339.             CALL FAUSMP(NESP,
  340.      &           GAMF,RF,PF,UNC,UTF,
  341.      &           GAMC,RC,PC,UNC,UTC,
  342.      &           YF.YET,YC.YET,
  343.      &           FLUX2D.FU,
  344.      &           CELLT)
  345. C-------------------------------------------------------
  346. C******* Residuum (son SPG a le meme ordre que MELEFC)
  347. C-------------------------------------------------------
  348.             MPRES.VPOCHA(IFAC,1)=FLUX2D.FU(1)*SURF/VOLU
  349.          MPRES.VPOCHA(IFAC,2)=((FLUX2D.FU(2)*CNX)+(FLUX2D.FU(3)*CTX))
  350.      &           *SURF/VOLU
  351.          MPRES.VPOCHA(IFAC,3)=((FLUX2D.FU(2)*CNY)+(FLUX2D.FU(3)*CTY))
  352.      &           *SURF/VOLU
  353.             MPRES.VPOCHA(IFAC,4)=FLUX2D.FU(4)*SURF/VOLU
  354.             do 1050  i=1,(nsp-1)
  355.               MPRES.VPOCHA(IFAC,4+I)=FLUX2D.FU(4+I)*SURF/VOLU
  356.  1050       continue
  357.          ELSE
  358.             CALL FAUSM3(NESP,
  359.      &           GAMF,RF,PF,UNC,UTF,UT2F,
  360.      &           GAMC,RC,PC,UNC,UTC,UT2C,
  361.      &           YF.YET,YC.YET,
  362.      &           FLUX3D.FU,
  363.      &           CELLT)
  364. C------------------------------------------------------
  365. C******* Residuum (son SPG a le meme ordre que MELEFC)
  366. C------------------------------------------------------
  367.             MPRES.VPOCHA(IFAC,1)=FLUX3D.FU(1)*SURF/VOLU
  368.          MPRES.VPOCHA(IFAC,2)=((FLUX3D.FU(2)*CNX)+(FLUX3D.FU(3)*CTX)+
  369.      &           (FLUX3D.FU(4)*CT2X))*SURF/VOLU
  370.          MPRES.VPOCHA(IFAC,3)=((FLUX3D.FU(2)*CNY)+(FLUX3D.FU(3)*CTY)+
  371.      &           (FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  372.          MPRES.VPOCHA(IFAC,4)=((FLUX3D.FU(2)*CNZ)+(FLUX3D.FU(3)*CTZ)+
  373.      &           (FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  374.             MPRES.VPOCHA(IFAC,5)=FLUX3D.FU(5)*SURF/VOLU
  375.             do 1060  i=1,(nsp-1)
  376.               MPRES.VPOCHA(IFAC,5+I)=FLUX3D.FU(5+I)*SURF/VOLU
  377.  1060       continue
  378.          ENDIF
  379.         ENDIF
  380.  1      CONTINUE
  381. C
  382.       SEGDES MELEFC
  383. C
  384.       SEGDES MLEMC
  385.       SEGDES MLEMCB
  386.       SEGDES MLEMF
  387. C
  388.       SEGDES MPNORM
  389.       SEGDES MPVOL
  390.       SEGDES MPSURF
  391.       SEGDES MPRC
  392.       SEGDES MPPC
  393.       SEGDES MPVC
  394.       SEGDES MPYN
  395.       SEGDES MPLIM
  396.       SEGDES MPRES
  397.       SEGDES MPRLI
  398.       SEGDES MLRECP
  399.       SEGDES MLRECV
  400.       SEGDES YC
  401.       IF(TEMP .GT. 0.0D0) SEGDES YF
  402.       SEGDES FLUX2D
  403.       SEGDES FLUX3D
  404. C
  405.       RETURN
  406.       END
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  

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