1 : C TNUMAC SOURCE CHAT 12/06/07 21:16:09 7389 2 : 3 : C======================================================================= 4 : C= T N U M A C = 5 : C= ----------- = 6 : C= = 7 : C= Fonction : = 8 : C= ---------- = 9 : C= Calcul de la matrice de CONDUCTIVITE THERMIQUE (conduction) pour = 10 : C= les elements finis MASSIFs a integration NUMERIQUE = 11 : C= = 12 : C= Parametres : (E)=Entree (S)=Sortie = 13 : C= ------------ = 14 : C= NEF (E) Numero de l'ELEMENT FINI dans NOMTP (cf. CCHAMP) = 15 : C= IMAIL (E) Numero du segment IMODEL dans le segment MMODEL = 16 : C= IPMODE (E) Pointeur sur un segment IMODEL suppose ACTIF = 17 : C= IPCHEM (E) Pointeur sur un segment MCHELM de CARACTERISTIQUES = 18 : C= IPRIGI (E/S) Pointeur sur l'objet RIGIDITE (CONDUCTIVITE) = 19 : C= = 20 : C= Variables locales : = 21 : C= ------------------- = 22 : C= NBNN Nombre de NOEUDS dans l'element considere = 23 : C= NbElem Nombre d'ELEMENTS dans le maillage elementaire = 24 : C= NbPGau Nombre de points de Gauss de l'element fini = 25 : C= XE(3,NBNN) Coordonnees GLOBALES des noeuds de l'ELEMENT = 26 : C= SHP(6,NBNN) Tableau de TRAVAIL = 27 : C= FORME(NBNN) Fonctions de FORME = 28 : C= NDIM Nombre de lignes de la matrice GRADIENT = 29 : C= GRAD(NDIM,NBNN) Matrice GRADIENT des fonctions de forme = 30 : C= CEL(NBNN,NBNN) Matrice de CONDUCTIVITE ELEMENTAIRE = 31 : C= CMAT(NDIM,NDIM) Matrice de CONDUCTIBILITE du materiau = 32 : C= VALMAT Valeurs des composantes de CONDUCTIBILITE et des = 33 : C= cos.directeurs des axes orthotropie/repere LOCAL = 34 : C= XGLOB Cos.directeurs des 3 axes d'ortho./repere GLOBAL = 35 : C= XLOC Cos.directeurs des 3 axes d'ortho./repere LOCAL = 36 : C= TXR Cos.directeurs des axes LOCAUX/repere GLOBAL = 37 : C= = 38 : C= Denis ROBERT, le 7 janvier 1988. P. DOWLATYARI, octobre 1990. = 39 : C======================================================================= 40 : 41 : SUBROUTINE TNUMAC (NEF,IPMAIL,IPINTE,IMATE,IVAMAT,NVAMAT, 42 : & IPMATR,NLIGR) 43 : 44 : IMPLICIT INTEGER(I-N) 45 : IMPLICIT REAL*8 (A-H,O-Z) 46 : 47 : -INC CCOPTIO 48 : -INC CCREEL 49 : -INC CCHAMP 50 : 51 : -INC SMCHAML 52 : -INC SMCOORD 53 : -INC SMELEME 54 : -INC SMINTE 55 : -INC SMRIGID 56 : 57 : SEGMENT MMAT1 58 : REAL*8 VALMAT(NV1) 59 : REAL*8 CEL(NLIGR,NLIGR),XE(3,NBNN) 60 : REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN) 61 : REAL*8 CMAT(NDIM,NDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM) 62 : ENDSEGMENT 63 : 64 : SEGMENT MAXE 65 : REAL*8 TXR(IDIM,IDIM),XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM) 66 : ENDSEGMENT 67 : 68 : SEGMENT MPTVAL 69 : INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU) 70 : CHARACTER*16 TYVAL(NCOSOU) 71 : ENDSEGMENT 72 : 73 : C 1 - INITIALISATIONS ET VERIFICATIONS 74 : C ====================================== 75 : MELEME = IPMAIL 76 : c* SEGACT,MELEME 77 : NBNN = NUM(/1) 78 : NBELEM = NUM(/2) 79 : C ===== 80 : MINTE = IPINTE 81 : c* SEGACT,MINTE 82 : NBPGAU = POIGAU(/1) 83 : C ===== 84 : MPTVAL = IVAMAT 85 : c* SEGACT,MPTVAL 86 : C ===== 87 : XMATRI = IPMATR 88 : c* SEGACT,XMATRI*MOD 89 : C ===== 90 : C Recuperation des fonctions de forme et de leurs derivees au 91 : C centre de gravite de l'element pour le calcul des axes locaux 92 : C d'orthotropie ou d'anisotropie 93 : C ===== 94 : IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN 95 : NLG = NUMGEO(NEF) 95 : CALL RESHPT(1,NBNN,NLG,NEF,0,IPINT1,IOK) 96 : C*OF IF (IOK.EQ.0) RETURN 97 : MINTE1 = IPINT1 98 : SEGACT,MINTE1 99 : NBSH = MINTE1.SHPTOT(/2) 100 : ENDIF 101 : C ===== 102 : C Initialisation des segments de travail 103 : C ===== 104 : IF (IFOMOD.EQ.1) THEN 105 : NDIM = 3 106 : ELSE 107 : NDIM = IDIM 108 : ENDIF 109 : NV1 = NVAMAT 110 : SEGINI,MMAT1 111 : IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN 112 : SEGINI,MAXE 113 : ENDIF 114 : 115 : C 2 - BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL 116 : C ============================================================ 117 : DO iElt = 1, NBELEM 118 : C ===== 119 : C 2.1 - Mise a zero de la matrice de CONDUCTIVITE de l'element iElt 120 : C ===== 120 : CALL ZERO(CEL,NBNN,NBNN) 121 : C ===== 122 : C 2.2 - Recuperation des coordonnees GLOBALES des noeuds de l'element 123 : C ===== 123 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE) 124 : C ===== 125 : C 2.3 - Calculs des axes locaux d'orthotropie ou d'anisotropie 126 : C ===== 127 : IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN 127 : CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR) 128 : IF (nbsh.EQ.-1) THEN 128 : CALL ERREUR(525) 129 : GOTO 9990 130 : ENDIF 131 : ENDIF 132 : C ===== 133 : C 2.4 - Boucle sur les points de Gauss de l'element iElt 134 : C ===== 135 : iFois=0 136 : DO iGau = 1, NBPGAU 137 : C ======= 138 : C 2.4.1 - Calcul du jacobien, des fonctions de forme et de leurs 139 : C derivees (matrice gradient) au point de Gauss iGau 140 : C ======= 140 : CALL TCOND5(iGau,NBNN,NDIM,XE,SHPTOT,SHP,GRAD,DJAC) 141 : IF (IERR.NE.0) GOTO 9990 142 : IF (DJAC.LT.XZERO) iFois=iFois+1 143 : DJAC=ABS(DJAC) 144 : C ======= 145 : C 2.4.2 - Erreur si le jacobien est nul en ce point de Gauss 146 : C ======= 147 : IF (DJAC.LT.XPETIT) THEN 148 : INTERR(1) = iElt 148 : CALL ERREUR(259) 149 : GOTO 9990 150 : ENDIF 151 : C ======= 152 : C 2.4.3 - Recuperation de la ou des valeurs de conductibilite au point 153 : C de Gauss iGau (tableau VALMAT) 154 : C ======= 155 : DJAC=DJAC*POIGAU(iGau) 156 : DO i = 1, NVAMAT 157 : IF (IVAL(i).NE.0) THEN 158 : MELVAL=IVAL(i) 159 : IBMN = MIN(iElt,VELCHE(/2)) 160 : IGMN = MIN(iGau,VELCHE(/1)) 161 : VALMAT(i) = VELCHE(IGMN,IBMN) 162 : ELSE 163 : VALMAT(i) = XZERO 164 : ENDIF 165 : ENDDO 166 : C ======= 167 : C 2.4.4 - Cas d'un materiau ISOTROPE de conductibilite K 168 : C Calcul de la contribution du point de Gauss a la matrice 169 : C CONDUCTIVITE elementaire de cet element fini 170 : C Ajout du terme XK*transposee(GRAD)*GRAD 171 : C Seul cas valide en dimension 1 172 : C ======= 173 : IF (IMATE.EQ.1) THEN 174 : XK = VALMAT(1)*DJAC 174 : CALL NTNST(GRAD,XK,NBNN,NDIM,CEL) 175 : C ======= 176 : C 2.4.5 - Cas d'un materiau ORTHOTROPE ou ANISOTROPE 177 : C Calcul de la contribution du point de Gauss a la matrice 178 : C CONDUCTIVITE elementaire de cet element fini 179 : C ======= 180 : ELSE 181 : C 2.4.5.1 - Calcul de la matrice de conductibilite (CMAT1) dans le 182 : C repere de "tropie" et de la matrice des cosinus directeurs 183 : C (XLOC) des directions de "tropie" dans le repere LOCAL 184 : C Cas des elements MASSIFS 2D 185 : IF (IDIM.EQ.2) THEN 186 : C= Cas d'un materiau ORTHOTROPE 187 : IF (IMATE.EQ.2) THEN 188 : CMAT1(1,1) = VALMAT(1) 189 : CMAT1(2,1) = XZERO 190 : CMAT1(1,2) = XZERO 191 : CMAT1(2,2) = VALMAT(2) 192 : XLOC(1,1) = VALMAT(3) 193 : XLOC(2,1) = VALMAT(4) 194 : XLOC(1,2) = -VALMAT(4) 195 : XLOC(2,2) = VALMAT(3) 196 : ELSE IF (IMATE.EQ.3) THEN 197 : C= Cas d'un materiau ANISOTROPE 198 : CMAT1(1,1) = VALMAT(1) 199 : CMAT1(2,1) = VALMAT(3) 200 : CMAT1(1,2) = VALMAT(3) 201 : CMAT1(2,2) = VALMAT(2) 202 : XLOC(1,1) = VALMAT(4) 203 : XLOC(2,1) = VALMAT(5) 204 : XLOC(1,2) = -VALMAT(5) 205 : XLOC(2,2) = VALMAT(4) 206 : ENDIF 207 : C 2.4.5.2 - Calcul de la matrice de conductibilite (CMAT1) dans le 208 : C repere de "tropie" et de la matrice des cosinus directeurs 209 : C (XLOC) des directions de "tropie" dans le repere LOCAL 210 : C Cas des elements MASSIFS 3D 211 : ELSE 212 : C= Cas d'un materiau ORTHOTROPE 213 : IF (IMATE.EQ.2) THEN 214 : CMAT1(1,1) = VALMAT(1) 215 : CMAT1(2,1) = XZERO 216 : CMAT1(3,1) = XZERO 217 : CMAT1(1,2) = XZERO 218 : CMAT1(2,2) = VALMAT(2) 219 : CMAT1(3,2) = XZERO 220 : CMAT1(1,3) = XZERO 221 : CMAT1(2,3) = XZERO 222 : CMAT1(3,3) = VALMAT(3) 223 : XLOC(1,1) = VALMAT(4) 224 : XLOC(2,1) = VALMAT(5) 225 : XLOC(3,1) = VALMAT(6) 226 : XLOC(1,2) = VALMAT(7) 227 : XLOC(2,2) = VALMAT(8) 228 : XLOC(3,2) = VALMAT(9) 229 : C= Cas d'un materiau ANISOTROPE 230 : ELSE IF (IMATE.EQ.3) THEN 231 : CMAT1(1,1) = VALMAT(1) 232 : CMAT1(2,1) = VALMAT(4) 233 : CMAT1(3,1) = VALMAT(5) 234 : CMAT1(1,2) = VALMAT(4) 235 : CMAT1(2,2) = VALMAT(2) 236 : CMAT1(3,2) = VALMAT(6) 237 : CMAT1(1,3) = VALMAT(5) 238 : CMAT1(2,3) = VALMAT(6) 239 : CMAT1(3,3) = VALMAT(3) 240 : XLOC(1,1) = VALMAT(7) 241 : XLOC(2,1) = VALMAT(8) 242 : XLOC(3,1) = VALMAT(9) 243 : XLOC(1,2) = VALMAT(10) 244 : XLOC(2,2) = VALMAT(11) 245 : XLOC(3,2) = VALMAT(12) 246 : ENDIF 247 : C= Calcul de la troisieme direction de "tropie" 247 : CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IOK) 248 : ENDIF 249 : C 2.4.5.3 - Calcul de la matrice (XGLOB) des cosinus directeurs des 250 : C axes de "tropie" par rapport au repere GLOBAL 250 : C*OF CALL ZERO(XGLOB,IDIM,IDIM) 251 : DO j=1,IDIM 252 : DO k=1,IDIM 253 : Cte = XZERO 254 : DO i=1,IDIM 255 : Cte=Cte+TXR(j,i)*XLOC(i,k) 256 : ENDDO 257 : XGLOB(k,j) = Cte 258 : ENDDO 259 : ENDDO 260 : C 2.4.5.4 - Calcul de la matrice de conductibilite (CMAT) dans le 261 : C repere GLOBAL (traitement particulier du mode 'FOURIER') 261 : CALL ZERO(CMAT,NDIM,NDIM) 262 : IF (IFOMOD.EQ.1) THEN 262 : CALL ZERO(CMAT2,IDIM,IDIM) 262 : CALL PRODT(CMAT2,CMAT1,XGLOB,IDIM,IDIM) 263 : CMAT(1,1)=CMAT2(1,1) 264 : IF (IMATE.EQ.2) THEN 265 : CMAT(2,2)=VALMAT(5) 266 : ELSE 267 : CMAT(2,2)=VALMAT(6) 268 : ENDIF 269 : CMAT(1,3)=CMAT2(1,2) 270 : CMAT(3,1)=CMAT2(1,2) 271 : CMAT(3,3)=CMAT2(2,2) 272 : ELSE 272 : CALL PRODT(CMAT,CMAT1,XGLOB,IDIM,IDIM) 273 : ENDIF 274 : C 2.4.5.5 - Contribution a la matrice de CONDUCTIVITE de cet element 275 : C Ajout du terme transposee(GRAD)*CMAT*GRAD 275 : CALL BDBST(GRAD,DJAC,CMAT,NBNN,NDIM,CEL) 276 : ENDIF 277 : ENDDO 278 : C ===== 279 : C 2.5 - Erreur si, en un point de Gauss, le jacobien change de signe. 280 : C ===== 281 : IF (iFois.NE.0.AND.iFois.NE.NBPGAU) THEN 282 : INTERR(1) = iElt 282 : CALL ERREUR(195) 283 : GOTO 9990 284 : ENDIF 285 : C ===== 286 : C 2.6 - Stockage de la matrice de CONDUCTIVITE pour cet element fini 287 : C ===== (remplissage de XMATRI) 287 : CALL REMPMT(CEL,NLIGR,RE(1,1,ielt)) 288 : C 289 : ENDDO 290 : 291 : C 3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS 292 : C ==================================================== 293 : 9990 CONTINUE 294 : SEGSUP,MMAT1 295 : IF (IMATE.EQ.2.OR.IMATE.EQ.3) THEN 296 : SEGDES,MINTE1 297 : SEGSUP,MAXE 298 : ENDIF 299 : 300 : RETURN 301 : END 302 : 303 : 304 :
© Cast3M 2003 - All rights reserved.
Disclaimer