Esope source files

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