Télécharger la plaquette Cast3M


fr en

Télécharger hdebit.eso

Retour à la liste

Numérotation des lignes : 

hdebit

  1. C HDEBIT    SOURCE    CB215821  24/04/12    21:16:14     11897          
  2.       SUBROUTINE HDEBIT
  3. C-----------------------------------------------------------------------
  4. C Calcul du débit à travers chaque face dans le cas d'une
  5. C formulation DARCY en éléments finis mixte hybride.
  6. C-----------------------------------------------------------------------
  7. C
  8. C Une face appartient en général à deux éléments. On ne définit qu'une
  9. C normale à la face. Le ddl aux faces est VN * S * signe(normale)
  10. C avec signe=1 si normale sortante et -1 sinon, S=Surface de la face et
  11. C VN=produit scalaire entre la vitesse à la face et la normale sortante.
  12. C
  13. C On obtient le débit en effectuant le produit du champoint de trace de
  14. C charge TH par une matrice - exprimée en fonction de la matrice de
  15. C rigidité de sous type DARCY et de la matrice DIV (matrice uniligne
  16. C composé de 1) - d'une part, par le MCHAML des orientations de normale
  17. C d'autre part.
  18. C
  19. C Comme nous utilisons des élémennts finis non conforme, il ne faut pas
  20. C assembler les quantitées calculées. Comme le flux est continu au signe
  21. C près, on ne stocke que le flux dans le sens de la normale calculée.
  22. C
  23. C          -1 t       (n)        -1    (n)
  24. C  Q  =  M1  * DIV * H       - M1  * TH
  25. C
  26. C
  27. C avec M1    : Matrice massse hybride
  28. C      DIV   : Matrice uniligne representant la divergence
  29. C
  30. C
  31. C---------------------------
  32. C Phrase d'appel (GIBIANE) :
  33. C---------------------------
  34. C
  35. C CHP1 = HDEB MMODEL MRIGI1 CHP2 CHP3  (MRIGI2 CHPO4) (CHP5);
  36. C
  37. C------------------------
  38. C Opérandes et résultat :
  39. C------------------------
  40. C
  41. C CHP1    : CHAMPOINT contenant le débit porté par la normale à la face
  42. C MMODEL  : Objet MMODEL spécifiant la formulation
  43. C MRIGI1  : RIGIDITE masses hybrides élémentaires de sous-type DARCY.
  44. C MRIGI2  : Matrices masses hybrides elementaires de sous-type MASSE.
  45. C CHP2    : CHPOINT centre contenant la charge H (ou concentration)
  46. C CHP3    : CHPOINT contenant la trace de charge TH
  47. C              (ou les traces de concentration)
  48. C CHPO4   : CHPO centre de composante FX,FY(,FZ),
  49. C                                  densite de force moyenne par élément
  50. C CHPO5   : CHPOIN face contenant le flux de la vitesse convective
  51. C
  52. C
  53. C-----------------------------------------------------------------------
  54. C
  55. C Langage : ESOPE + FORTRAN77
  56. C
  57. C Auteurs : 09/93 F.DABBENE - Cas permanent
  58. C           09/94 X.NOUVELLON - Extension au cas transitoire
  59. C           02/96 L.V.BENET - Prise en compte de forces de volume
  60. C           06/98 F.AURIOL - formulation prenant en compte les valeurs
  61. C                           des charges ou des concentrations au centre
  62. C                           des elements. Prise en compte de la vitesse
  63. C                           convective
  64. C
  65. C-----------------------------------------------------------------------
  66.       IMPLICIT INTEGER(I-N)
  67.       IMPLICIT REAL*8 (A-H,O-Z)
  68. *
  69.  
  70. -INC PPARAM
  71. -INC CCOPTIO
  72. -INC SMCHAML
  73. -INC SMCHPOI
  74. -INC SMELEME
  75. -INC SMMODEL
  76. -INC SMRIGID
  77. -INC SMTABLE
  78. *
  79.       SEGMENT IPMAHY
  80.          INTEGER MAHYBR(NSOUS)
  81.       ENDSEGMENT
  82. *
  83.       LOGICAL LOGRE,LOGIN
  84.       CHARACTER*4 NOMTOT(1)
  85.       CHARACTER*4 NOMTO3(3)
  86.       CHARACTER*4 NOMTH
  87.       CHARACTER*6 CHAR6
  88.       CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  89. *
  90. * Initialisations
  91. *
  92.       IVALIN = 0
  93.       XVALIN = 0.D0
  94.       LOGIN  = .TRUE.
  95.       IOBIN  = 0
  96.       TAPIND = 'MOT     '
  97. *
  98. * Lecture du MMODEL et test de la formulation
  99. *
  100.       CALL LIROBJ('MMODEL ',IPMODE,1,IRET)
  101.       CALL ACTOBJ('MMODEL ',IPMODE,1)
  102.       IF (IERR.NE.0) RETURN
  103.       MMODEL = IPMODE
  104.       SEGACT MMODEL
  105.       N1=KMODEL(/1)
  106.       DO 7 I=1,N1
  107.          IMODEL=KMODEL(I)
  108.          SEGACT IMODEL
  109.          IF(FORMOD(1).NE.'DARCY')THEN
  110.              MOTERR(1:16)  = 'DARCY           '
  111.              CALL ERREUR(719)
  112.              RETURN
  113.          ENDIF
  114.     7 CONTINUE
  115. *
  116. * Lecture de la TABLE domaine
  117. *
  118.       IPTABL = 0
  119.       CALL LEKMOD(MMODEL,IPTABL,IRET)
  120.       IF (IERR.NE.0) RETURN
  121.       CHARIN = 'MAILLAGE'
  122.       TYPOBJ = 'MAILLAGE'
  123.       CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  124.       IF (IERR.NE.0) RETURN
  125.       IPGEOM = IOBRE
  126.       CALL LEKTAB(IPTABL,'ELTFA',IOBRE)
  127.       IF (IERR.NE.0) RETURN
  128.       IELTFA = IOBRE
  129.       CALL LEKTAB(IPTABL,'CENTRE',IOBRE)
  130.       IF (IERR.NE.0) RETURN
  131.       ICENTR = IOBRE
  132.       CALL LEKTAB(IPTABL,'FACE',IOBRE)
  133.       IF (IERR.NE.0) RETURN
  134.       IPFACE = IOBRE
  135. *
  136. * récup. MCHAML orientation normale
  137. *
  138.       CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  139.       IF (IERR.NE.0) RETURN
  140. *
  141. * récup. CHPO orientation normale
  142. *
  143.       CALL LEKTAB(IPTABL,'XXNORMAF',INORM)
  144.       IF (IERR.NE.0) RETURN
  145. *
  146. * récup. CHPO surface des faces
  147. *
  148.       CALL LEKTAB(IPTABL,'XXSURFAC',ISURF)
  149.       IF (IERR.NE.0) RETURN
  150. *
  151. * Lecture de RIGIDITE  obligatoire
  152. *
  153.       CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET)
  154.       IF (IERR.NE.0) RETURN
  155.       RI1 = IPRIGI
  156. *
  157. * Lecture du CHAMPOINT des charges (ou concentrations au centre)
  158. *
  159.       CALL LIROBJ('CHPOINT ',ITPP,1,IRET)
  160.       CALL ACTOBJ('CHPOINT ',ITPP,1)
  161. *
  162.       ICHP1 = ITPP
  163. *
  164. * Lecture du CHAMPOINT de concentration aux facees
  165. *
  166.       CALL LIROBJ('CHPOINT ',ITPN,1,IRET)
  167.       CALL ACTOBJ('CHPOINT ',ITPN,1)
  168.       IF (IERR.NE.0) RETURN
  169. *
  170. * Lecture eventuelle du 1er CHAMPOINT optionnel
  171. *
  172.       ICHP2 = 0
  173.       CALL LIROBJ('CHPOINT',ICHP2,0,IRET2)
  174.       IF(IRET2 .EQ. 1) CALL ACTOBJ('CHPOINT',ICHP2,1)
  175.       IF (IERR.NE.0) RETURN
  176. *
  177. * Lecture eventuelle du 2em CHAMPOINT optionnel
  178. *
  179.       ICHP3 = 0
  180.       CALL LIROBJ('CHPOINT ',ICHP3,0,IRET3)
  181.       IF(IRET3 .EQ. 1) CALL ACTOBJ('CHPOINT ',ICHP3,1)
  182.       IF (IERR.NE.0) RETURN
  183. *
  184. * Lecture de la RIGIDITE complementaire
  185. *
  186.       IPRIGI = 0
  187.       CALL LIROBJ('RIGIDITE',IPRIGI,0,IRETR)
  188.       IF (IERR.NE.0) RETURN
  189.       RI2 = IPRIGI
  190. *
  191. * Test de la formulation contenu dans MMODEL
  192. *
  193.       SEGACT MMODEL
  194.       NSOUS  = KMODEL(/1)
  195.       SEGINI IPMAHY
  196.       IDARCY = 0
  197.       IPT1= IPGEOM
  198.       IPT2= IPGEOM
  199.       DO 10 ISOUS=1,NSOUS
  200.          IF(NSOUS.GT.1) THEN
  201.               SEGACT IPT2
  202.               IPT1= IPT2.LISOUS(ISOUS)
  203.          ENDIF
  204.          IMODEL = KMODEL(ISOUS)
  205.          SEGACT IMODEL
  206.          LETYPE = FORMOD(1)
  207.          IF (LETYPE.EQ.'DARCY') THEN
  208.             IDARCY = IDARCY + 1
  209.             MAHYBR(ISOUS) = IPT1
  210.          ENDIF
  211.  10   CONTINUE
  212.       IF (IDARCY.EQ.0) THEN
  213.          MOTERR = LETYPE
  214.          CALL ERREUR(193)
  215.          GOTO 100
  216.       ENDIF
  217. *
  218. * Récuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  219. *
  220.       MELEME = IELTFA
  221.       SEGACT MELEME
  222.       LZONES = LISOUS(/1)
  223.       IF (LZONES.EQ.0) LZONES = 1
  224.       IPT4   = IPGEOM
  225.       SEGACT IPT4
  226.       DO 30 ISOUS=1,NSOUS
  227.          IMAMEL = MAHYBR(ISOUS)
  228.          IF (IMAMEL.NE.0) THEN
  229.             DO 20 ISZ=1,LZONES
  230.                IF (LZONES.EQ.1) THEN
  231.                   IPT2 = IPT4
  232.                   IPT3 = MELEME
  233.                ELSE
  234.                   IPT2 = IPT4.LISOUS(ISZ)
  235.                   IPT3 = LISOUS(ISZ)
  236.                ENDIF
  237.                IF (IPT2.EQ.IMAMEL) THEN
  238.                   MAHYBR(ISOUS) = IPT3
  239.                   GOTO 30
  240.                ENDIF
  241.  20         CONTINUE
  242.             IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  243.                INTERR(1) = ISOUS
  244.                CALL ERREUR(698)
  245.                GOTO 100
  246.             ENDIF
  247.          ENDIF
  248.  30   CONTINUE
  249. *
  250. * Test et tri du sous-type des matrices de rigiditée récupérées
  251. *
  252.       IPFORC=0
  253.       IPDARC=0
  254.       SEGACT RI1
  255.       LETYPE =  RI1.MTYMAT
  256.       IF(LETYPE.EQ.'DARCY')THEN
  257.         IPDARC=RI1
  258.       ELSE IF(LETYPE.EQ.'MASSE')THEN
  259.         IPFORC=RI1
  260.       ENDIF
  261.       IF(RI2.NE.0)THEN
  262.         SEGACT RI2
  263.         LETYPE =  RI2.MTYMAT
  264.         IF(LETYPE.EQ.'MASSE')THEN
  265.            IPFORC=RI2
  266.         ELSE IF(LETYPE.EQ.'DARCY')THEN
  267.            IPDARC=RI2
  268.         ENDIF
  269.         IF(IPFORC.EQ.0)THEN
  270.           MOTERR(1:8)  = 'RIGIDITE'
  271.           MOTERR(9:16) = 'MASSE'
  272.           CALL ERREUR(79)
  273.           SEGDES RI1
  274.           SEGDES RI2
  275.           GOTO 100
  276.         ENDIF
  277.       ENDIF
  278.       IF(IPDARC.EQ.0)THEN
  279.         MOTERR(1:8)  = 'RIGIDITE'
  280.         MOTERR(9:16) = 'DARCY'
  281.         CALL ERREUR(79)
  282.         SEGDES RI1
  283.         IF(IPFORC.NE.0)SEGDES RI2
  284.         GOTO 100
  285.       ENDIF
  286. *
  287. * Controle des pointeurs de MELEME de la rigidité
  288. *
  289.       DO 40 ISOUS=1,NSOUS
  290.          IMAMEL = MAHYBR(ISOUS)
  291.          IF (IMAMEL.NE.0) THEN
  292.             RI1=IPDARC
  293.             IPTEST = RI1.IRIGEL(1,ISOUS)
  294.             CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  295.             IF(.NOT. LOGRE)THEN
  296.                MOTERR(1:8)  = 'DARCY'
  297.                MOTERR(9:16) = ' ELTFA  '
  298.                INTERR(1)    = ISOUS
  299.                CALL ERREUR(698)
  300.                SEGDES RI1
  301.                GOTO 100
  302.             ENDIF
  303.             IF(IPFORC.NE.0)THEN
  304.               RI2=IPFORC
  305.               IPTEST = RI2.IRIGEL(1,ISOUS)
  306.                CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  307.                IF(.NOT. LOGRE)THEN
  308.                  MOTERR(1:8)  = 'MASSE'
  309.                  MOTERR(9:16) = ' ELTFA  '
  310.                  INTERR(1)    = ISOUS
  311.                  CALL ERREUR(698)
  312.                  SEGDES RI1
  313.                  SEGDES RI2
  314.                  GOTO 100
  315.                ENDIF
  316.             ENDIF
  317.          ENDIF
  318.  40   CONTINUE
  319.       SEGDES RI1
  320.       IF(IPFORC.NE.0)SEGDES RI2
  321. *
  322. *  Test du  CHAMPOINT de valeurs aux faces
  323. *
  324.       IFORC=0
  325.       MCHPOI=ITPN
  326.       SEGACT MCHPOI
  327.       NBSOUS=IPCHP(/1)
  328.       ICHP0= ITPN
  329.       IITH=0
  330.       IVCONV=0
  331.       IF(NBSOUS.NE.1)THEN
  332. C     c'est le cas en presence de multiplicateurs de Lagrange
  333. C     on va extraire la composante TH
  334.          NOMTH='TH  '
  335.          CALL EXCOPP(ITPN,NOMTH,NIFOUR,ITPTH,NOMTH,NIFOUR,0)
  336.          IF(IERR.NE.0)GO TO 100
  337.          ICHP0= ITPTH
  338.          IITH=1
  339.       ENDIF
  340.       NOMTOT(1)='    '
  342.       NBCOMP=-1
  343.       CALL QUEPOI(ICHP0,IPFACE,INDIC,NBCOMP,NOMTOT)
  344.       IF(IERR.NE.0)GO TO 100
  345.       IF(INDIC.LT.0)THEN
  346.          IRETOU=0
  347.          IPOINT=0
  348.          CHARIN='        '
  349.          CALL SKNAME(ITPN,CHARIN,IRETOU,IPOINT)
  350.          MOTERR(1:8)= CHARIN(1:8)
  351.          IF(IRETOU.EQ.1)CALL ERREUR(-301)
  352.          CALL ERREUR(609)
  353.          GO TO 100
  354.        ENDIF
  355.        IF(IITH.EQ.0)THEN
  356.           IF(NBSOUS.EQ.1)THEN
  357.              MCHPOI=ICHP0
  358.              SEGACT MCHPOI
  359.              MSOUPO=IPCHP(1)
  360.              SEGACT MSOUPO
  361.              NC=NOCOMP(/2)
  362.              IF(NC.EQ.1)THEN
  363.                 IF(NOCOMP(1).EQ.'TH  ')IITH=1
  364.              ENDIF
  365.           ENDIF
  366.        ENDIF
  367. *
  368. * Tri et test des CHAMPOINTs ARGUMENTS optionnels
  369. *
  370.        IF(IPFORC.NE.0) THEN
  371.           IF(IITH.NE.1)THEN
  372.              IRETOU=0
  373.              IPOINT=0
  374.              CHARIN='        '
  375.              CALL SKNAME(ITPN,CHARIN,IRETOU,IPOINT)
  376.              MOTERR(1:4)= 'TH  '
  377.              MOTERR(5:12)= CHARIN(1:8)
  378.              CALL ERREUR(77)
  379.           ENDIF
  380.           IF(IRET2.EQ.1)THEN
  381.              IFORC=ICHP2
  382.              NOMTOT(1)='    '
  384.              NBCOMP=-1
  385.              CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTOT)
  386.              IF(IERR.NE.0)GO TO 100
  387.              IF(INDIC.LT.0)THEN
  388.                  CALL ERREUR(348)
  389.                  GO TO 100
  390.              ENDIF
  391.           ELSE
  392.              CALL ERREUR(641)
  393.              GO TO 100
  394.           ENDIF
  395.           NBCOMP    = 0
  396.           INDIC     = 1
  397.           NOMTO3(1) = ' '
  398.           NOMTO3(2) = ' '
  399.           NOMTO3(3) = ' '
  400.           CALL QUEPOI(IFORC,ICENTR,INDIC,NBCOMP,NOMTO3)
  401.           IF(NOMTO3(1).NE.'FX')THEN
  402.              MOTERR(1:4)= NOMTO3(1)
  403.              CALL ERREUR(197)
  404.           ENDIF
  405.           IF(IRET3.NE.0)THEN
  406.              NBCOMP    = -1
  407.              INDIC     = 0
  408.              NOMTOT(1) = ' '
  409.              CALL QUEPOI(ICHP3,IPFACE,INDIC,NBCOMP,NOMTOT)
  410.              IF(IERR.NE.0)GO TO 100
  411.              IF(INDIC.LT.0)THEN
  412.                  CALL ERREUR(348)
  413.                  GO TO 100
  414.              ENDIF
  415.              IVCONV=ICHP3
  416.           ENDIF
  417. C     Cas ou il n y a pas de rigidite RI2
  418.       ELSE
  419.          NOMTOT(1)='    '
  421.          NBCOMP=-1
  422.          CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTOT)
  423.          IF(IERR.NE.0)GO TO 100
  424.          IF(INDIC.LT.0)THEN
  425.              CALL ERREUR(348)
  426.              GO TO 100
  427.          ENDIF
  428.          MCHPOI=ICHP0
  429.          SEGACT MCHPOI
  430.          MSOUPO=IPCHP(1)
  431.          SEGACT MSOUPO
  432.          NC=NOCOMP(/2)
  433.          MCHPO1= ICHP1
  434.          SEGACT MCHPO1
  435.          MSOUP1=MCHPO1.IPCHP(1)
  436.          SEGACT MSOUP1
  437.          NC1=MSOUP1.NOCOMP(/2)
  438.          IF(NC1.NE.NC)THEN
  439.             MOTERR(1:40)=
  440.      &         'pas le meme nombre de composantes       '
  441.             CALL ERREUR(-301)
  442.             CALL ERREUR(21)
  443.             GO TO 100
  444.          ENDIF
  445.          IF(IITH.EQ.1) THEN
  446.          IF((MSOUP1.NOCOMP(1)).NE.'H   ')THEN
  447.             MOTERR(1:40)=
  448.      &         'le nom de la composante doit etre H     '
  449.             CALL ERREUR(-301)
  450.             CALL ERREUR(21)
  451.             GO TO 100
  452.          ENDIF
  453.          IF(IRET2.NE.0)THEN
  454.          NBCOMP    = -1
  455.          INDIC     = 0
  456.          NOMTOT(1) = ' '
  457.          CALL QUEPOI(ICHP2,IPFACE,INDIC,NBCOMP,NOMTOT)
  458.          IF(IERR.NE.0)GO TO 100
  459.          IF(INDIC.LT.0)THEN
  460.              CALL ERREUR(348)
  461.              GO TO 100
  462.          ENDIF
  463.          IVCONV=ICHP2
  464.          ENDIF
  465.          ELSE
  466.             NNC1=0
  467.             DO 50 I=1,NC
  468.                DO 35 J=1,NC1
  469.                   IF(NOCOMP(I).EQ.MSOUP1.NOCOMP(J))THEN
  470.                      NNC1=NNC1+1
  471.                      GO TO 32
  472.                   ENDIF
  473.    35          CONTINUE
  474.                MOTERR(1:40)=
  475.      &            'les noms de composantes sont differents '
  476.                CALL ERREUR(-301)
  477.                CALL ERREUR(21)
  478.                GO TO 100
  479.    32          CONTINUE
  480.    50       CONTINUE
  481.             IF(NNC1.NE.NC1)THEN
  482.                MOTERR(1:40)=
  483.      &            'un des CHPO a plus de comp. que l autre '
  484.                CALL ERREUR(21)
  485.                GO TO 100
  486.             ENDIF
  487.             IVCONV=ICHP2
  488.          ENDIF
  489.       ENDIF
  490. *
  491. * Construction du CHAMPOINT des debits aux faces
  492. *
  493.       SEGDES IPMAHY
  494.       CALL HDEBI2(IPMAHY,IPFACE,IPDARC,ICHP0,IORIE,
  495.      S     ICHP1,IVCONV,IITH,INORM,IORIE,ISURF,IPFORC,IFORC,IRET)
  496.  
  497.  
  498.       IF (IRET.NE.0) THEN
  499.         CALL ACTOBJ('CHPOINT ',IRET,1)
  500.         CALL ECROBJ('CHPOINT ',IRET)
  501.       ENDIF
  502. *
  503. * Ménage
  504. *
  505.  100  CONTINUE
  506.       SEGSUP IPMAHY
  507.       END
  508.  
  509.  
  510.  
  511.  
  512.  

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