Esope source files

   1 : C CAPAC3    SOURCE    FANDEUR   11/07/19    21:15:27     7042           
   2 : 
   3 : C=======================================================================
   4 : C=                            C A P A C 3                              =
   5 : C=                            -----------                              =
   6 : C=                                                                     =
   7 : C=  Fonction :                                                         =
   8 : C=  ----------                                                         =
   9 : C=  Calcul de la matrice de CAPACITE CALORIFIQUE pour des elements de  =
  10 : C=  COQUE TRIANGLE (COQ3) a integration semi-analytique                =
  11 : C=                                                                     =
  12 : C=  Parametres :  (E)=Entree  (S)=Sortie                               =
  13 : C=  ------------                                                       =
  14 : C=   NEF      (E)    Numero de l'ELEMENT FINI dans NOMTP               =
  15 : C=   IPMAIL   (E)    Numero du segment IMODEL dans le segment MMODEL   =
  16 : C=   CLAT     (E)    Chaleur latente du changement de phase            =
  17 : C=   IPRIGI   (E/S)  Matrice de CAPACITE (RIGIDITE) resultat (ACTIF)   =
  18 : C=======================================================================
  19 : 
  20 :       SUBROUTINE CAPAC3 (NEF,IPMAIL,IPINTE,IVAMAT,NVAMAT,IVAPHA,NVAPHA,
  21 :      &                   IPMATR,NLIGR)
  22 : 
  23 :       IMPLICIT INTEGER(I-N)
  24 :       IMPLICIT REAL*8 (A-H,O-Z)
  25 : 
  26 : -INC CCOPTIO
  27 : -INC CCREEL
  28 : -INC CCHAMP
  29 : 
  30 : -INC SMCHAML
  31 : -INC SMCOORD
  32 : -INC SMELEME
  33 : -INC SMINTE
  34 : -INC SMRIGID
  35 : 
  36 :       SEGMENT MMAT1
  37 :         REAL*8 XE(3,NBNN),FORME(NBNN)
  38 :         REAL*8 CAPSS(NBNN,NBNN),CAPV(NLIGR,NLIGR)
  39 :       ENDSEGMENT
  40 : 
  41 :       SEGMENT MPTVAL
  42 :         INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU)
  43 :         CHARACTER*16 TYVAL(NCOSOU)
  44 :       ENDSEGMENT
  45 : 
  46 : C= Quelques constantes numeriques
  47 :       PARAMETER (X1s15=0.066666666666666666666666666667D0)
  48 :       PARAMETER (X2s15=0.133333333333333333333333333333D0)
  49 :       PARAMETER (X8s15=0.533333333333333333333333333333D0)
  50 :       PARAMETER (X1s30=0.033333333333333333333333333333D0)
  51 : 
  52 :       REAL*8 VACOMP(3)
  53 : 
  54 : C  1 - INITIALISATIONS ET VERIFICATIONS
  55 : C ======================================
  56 :       MELEME = IPMAIL
  57 : c*      SEGACT,MELEME
  58 :       NBNN   = NUM(/1)
  59 :       NBELEM = NUM(/2)
  60 :       NBNN2 = 2*NBNN
  61 : c*      NBNN3 = 3*NBNN
  62 : C =====
  63 :       MINTE = IPINTE
  64 : c*      SEGACT,MINTE
  65 :       NBPGAU = POIGAU(/1)
  66 : C =====
  67 :       MPTVAL = IVAMAT
  68 : c*      SEGACT,MPTVAL
  69 : C- Test sur la constance du champ d'epaisseur : supprime
  70 : c*      IPMELV = IVAL(3)
  70 : c*      CALL QUELCH(IPMELV,IOK)
  71 : c*      IF (IOK.NE.0) THEN
  71 : c*        CALL ERREUR(566)
  72 : c*        GOTO 9990
  73 : c*      ENDIF
  74 : C =====
  75 : c*      IF (IVAPHA.NE.0) THEN
  76 : c*        MPTVAL = IVAPHA
  77 : c*        SEGACT,MPTVAL
  78 : c*      ENDIF
  79 : C =====
  80 :       XMATRI = IPMATR
  81 : c*      SEGACT,XMATRI*MOD
  82 : c*      NLIGRP = NBNN3 = NLIGR
  83 : c*      NLIGRD = NBNN3 = NLIGR
  84 : C =====
  85 :       SEGINI,MMAT1
  86 : 
  87 : C  2 - BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL
  88 : C ============================================================
  89 :       DO iElt = 1, NBELEM
  90 : C =====
  91 : C  2.1 - Recuperation des coordonnees GLOABLES des noeuds de l'element
  92 : C =====
  92 :         CALL ZERO(CAPV,NLIGR,NLIGR)
  93 : C =====
  94 : C  2.2 - Recuperation des coordonnees GLOABLES des noeuds de l'element
  95 : C =====
  95 :         CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE)
  96 : C =====
  97 : C  2.3 - Boucle sur les points de Gauss de l'element iElt
  98 : C =====
  99 :         DO iGau = 1, NBPGAU
 100 : 
 101 : C- Calcul du volume associe a ce point de Gauss (jacobien)
 102 :           S1=XZero
 103 :           S2=XZero
 104 :           S3=XZero
 105 :           S4=XZero
 106 :           S5=XZero
 107 :           S6=XZero
 108 :           DO iNoe = 1, NBNN
 109 :             S1=S1+SHPTOT(2,iNoe,iGau)*XE(2,iNoe)
 110 :             S2=S2+SHPTOT(3,iNoe,iGau)*XE(3,iNoe)
 111 :             S3=S3+SHPTOT(3,iNoe,iGau)*XE(2,iNoe)
 112 :             S4=S4+SHPTOT(2,iNoe,iGau)*XE(3,iNoe)
 113 :             S5=S5+SHPTOT(3,iNoe,iGau)*XE(1,iNoe)
 114 :             S6=S6+SHPTOT(2,iNoe,iGau)*XE(1,iNoe)
 115 :           ENDDO
 116 :           SurfX=S1*S2-S3*S4
 117 :           SurfY=S4*S5-S2*S6
 118 :           SurfZ=S6*S3-S5*S1
 119 :           DJAC = ABS(SurfX*SurfX+SurfY*SurfY+SurfZ*SurfZ)
 120 : C- Verification que le volume n'est pas nul en ce point de Gauss
 121 :           IF (DJAC.LT.XPETIT) THEN
 122 :             INTERR(1) = iElt
 122 :             CALL ERREUR(259)
 123 :             GOTO 9990
 124 :           ENDIF
 125 :           DJAC = SQRT(DJAC)
 126 : C- Calcul du terme Rho.Cp.Vol en ce point de Gauss
 127 : C- (avec prise en compte d'un eventuel changement de phase)
 127 :           CALL CAPA5(iGau,iElt,IVAMAT,NVAMAT,IVAPHA,NVAPHA, VACOMP)
 128 :           CAPA = DJAC * POIGAU(iGau) * VACOMP(1)
 129 : C- Calcul de la contribution du point de Gauss a la matrice
 130 : C- CAPACITE elementaire pour cet element fini
 130 :           CALL ZERO(CAPSS,NBNN,NBNN)
 131 :           do iou=1,nbnn
 132 :             forme(iou)=shptot(1,iou,igau)
 133 :           enddo
 133 :           CALL NTNST(FORME,CAPA,NBNN,1,CAPSS)
 134 : C- Ajout de termes specifiques dus a l'integration (analytique)
 135 : C- suivant l'epaisseur de l'element de type COQUE
 136 : C =======
 137 :           EP = VACOMP(3)
 138 : C- Erreur si l'epaisseur est est nulle
 139 : c*          IF (EP.LE.XPetit) THEN
 139 : c*            CALL ERREUR(517)
 140 : c*            GOTO 9990
 141 : c*          ENDIF
 142 :           C1 =  X2s15*EP
 143 :           C2 =  X1s15*EP
 144 :           C3 = -X1s30*EP
 145 :           C4 =  X8s15*EP
 146 :           C5 =  C2
 147 :           C6 =  C1
 148 :           DO j=1,NBNN
 149 :             j1 = j + NBNN
 150 :             j2 = j + NBNN2
 151 :             DO i=1,NBNN
 152 :               i1 = i + NBNN
 153 :               i2 = i + NBNN2
 154 :               Cte = CAPSS(i,j)
 155 :               CAPV( i, j) = CAPV( i, j) + C1*Cte
 156 :               CAPV(i1, j) = CAPV(i1, j) + C2*Cte
 157 :               CAPV(i2, j) = CAPV(i2, j) + C3*Cte
 158 :               CAPV(i1,j1) = CAPV(i1,j1) + C4*Cte
 159 :               CAPV(i2,j1) = CAPV(i2,j1) + C5*Cte
 160 :               CAPV(i2,j2) = CAPV(i2,j2) + C6*Cte
 161 :             ENDDO
 162 :           ENDDO
 163 :         ENDDO
 164 : C =====
 165 : C  2.4 - Stockage de la matrice de CAPACITE pour cet element fini
 166 : C        (remplissage de XMATRI)
 167 : C =====
 167 :         CALL REMPMT(CAPV,NLIGR,RE(1,1,iElt))
 168 :       ENDDO
 169 : 
 170 : C  3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS
 171 : C ====================================================
 172 :  9990 CONTINUE
 173 :       SEGSUP,MMAT1
 174 : c*      SEGDES,MELEME,MINTE,XMATRI
 175 : c*      MPTVAL = IVAMAT
 176 : c*      SEGDES,MPTVAL
 177 : c*      IF (IVAPHA.NE.0) THEN
 178 : c*        MPTVAL = IVAPHA
 179 : c*        SEGDES,MPTVAL
 180 : c*      ENDIF
 181 : 
 182 :       RETURN
 183 :       END
 184 : 
 185 :