Télécharger tnumac.eso

Retour à la liste

tnumac
C TNUMAC    SOURCE    CB215821  17/06/06    21:15:11     9449           
 
C=======================================================================
C=                            T N U M A C                              =
C=                            -----------                              =
C=                                                                     =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=  Calcul de la matrice de CONDUCTIVITE THERMIQUE (conduction) pour   =
C=  les elements finis MASSIFs a integration NUMERIQUE                 =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   NEF      (E)   Numero de l'ELEMENT FINI dans NOMTP (cf. CCHAMP)   =
C=   IMAIL    (E)   Numero du segment IMODEL dans le segment MMODEL    =
C=   IPMODE   (E)   Pointeur sur un segment IMODEL suppose ACTIF       =
C=   IPCHEM   (E)   Pointeur sur un segment MCHELM de CARACTERISTIQUES =
C=   IPRIGI  (E/S)  Pointeur sur l'objet RIGIDITE (CONDUCTIVITE)       =
C=                                                                     =
C=  Variables locales :                                                =
C=  -------------------                                                =
C=   NBNN            Nombre de NOEUDS dans l'element considere         =
C=   NbElem          Nombre d'ELEMENTS dans le maillage elementaire    =
C=   NbPGau          Nombre de points de Gauss de l'element fini       =
C=   XE(3,NBNN)      Coordonnees GLOBALES des noeuds de l'ELEMENT      =
C=   SHP(6,NBNN)     Tableau de TRAVAIL                                =
C=   FORME(NBNN)     Fonctions de FORME                                =
C=   NDIM            Nombre de lignes de la matrice GRADIENT           =
C=   GRAD(NDIM,NBNN) Matrice GRADIENT des fonctions de forme           =
C=   CEL(NBNN,NBNN)  Matrice de CONDUCTIVITE ELEMENTAIRE               =
C=   CMAT(NDIM,NDIM) Matrice de CONDUCTIBILITE du materiau             =
C=   VALMAT          Valeurs des composantes de CONDUCTIBILITE et des  =
C=                   cos.directeurs des axes orthotropie/repere LOCAL  =
C=   XGLOB           Cos.directeurs des 3 axes d'ortho./repere GLOBAL  =
C=   XLOC            Cos.directeurs des 3 axes d'ortho./repere LOCAL   =
C=   TXR             Cos.directeurs des axes LOCAUX/repere GLOBAL      =
C=                                                                     =
C=  Denis ROBERT, le 7 janvier 1988. P. DOWLATYARI, octobre 1990.      =
C=======================================================================
 
      SUBROUTINE TNUMAC (NEF,IPMAIL,IPINTE,IMATE,IVAMAT,NVAMAT,
     &                   IPMATR,NLIGR)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC CCHAMP
 
-INC SMCHAML
-INC SMCOORD
-INC SMELEME
-INC SMINTE
-INC SMRIGID
 
      SEGMENT MMAT1
        REAL*8 VALMAT(NV1)
        REAL*8 CEL(NLIGR,NLIGR),XE(3,NBNN)
        REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN)
        REAL*8 CMAT(NDIM,NDIM),CMAT1(IDIM,IDIM),CMAT2(IDIM,IDIM)
      ENDSEGMENT
 
      SEGMENT MAXE
        REAL*8 TXR(IDIM,IDIM),XLOC(3,3),XGLOB(3,3)
      ENDSEGMENT
 
      SEGMENT MPTVAL
        INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU)
        CHARACTER*16 TYVAL(NCOSOU)
      ENDSEGMENT
 
C  1 - INITIALISATIONS ET VERIFICATIONS
C ======================================
      MELEME = IPMAIL
c*      SEGACT,MELEME
      NBNN   = NUM(/1)
      NBELEM = NUM(/2)
C =====
      MINTE = IPINTE
c*      SEGACT,MINTE
      NBPGAU = POIGAU(/1)
C =====
      MPTVAL = IVAMAT
c*      SEGACT,MPTVAL
C =====
      XMATRI = IPMATR
c*      SEGACT,XMATRI*MOD
C =====
C Recuperation des fonctions de forme et de leurs derivees au
C        centre de gravite de l'element pour le calcul des axes locaux
C        d'orthotropie ou d'anisotropie
C =====
      IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN
        NLG = NUMGEO(NEF)
        CALL RESHPT(1,NBNN,NLG,NEF,0,IPINT1,IOK)
C*OF    IF (IOK.EQ.0) RETURN
        MINTE1 = IPINT1
        SEGACT,MINTE1
        NBSH = MINTE1.SHPTOT(/2)
      ENDIF
C =====
C Initialisation des segments de travail
C =====
      IF (IFOMOD.EQ.1) THEN
        NDIM = 3
      ELSE
        NDIM = IDIM
      ENDIF
      NV1 = NVAMAT
      SEGINI,MMAT1
      IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN
        SEGINI,MAXE
      ENDIF
 
C  2 - BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL
C ============================================================
      DO iElt = 1, NBELEM
C =====
C  2.1 - Mise a zero de la matrice de CONDUCTIVITE de l'element iElt
C =====
        CALL ZERO(CEL,NBNN,NBNN)
C =====
C  2.2 - Recuperation des coordonnees GLOBALES des noeuds de l'element
C =====
        CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE)
C =====
C  2.3 - Calculs des axes locaux d'orthotropie ou d'anisotropie
C =====
        IF (IMATE.EQ.2 .OR. IMATE.EQ.3) THEN
          CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR)
          IF (nbsh.EQ.-1) THEN
            CALL ERREUR(525)
            GOTO 9990
          ENDIF
        ENDIF
C =====
C  2.4 - Boucle sur les points de Gauss de l'element iElt
C =====
        iFois=0
        DO iGau = 1, NBPGAU
C =======
C  2.4.1 - Calcul du jacobien, des fonctions de forme et de leurs
C          derivees (matrice gradient) au point de Gauss iGau
C =======
          CALL TCOND5(iGau,NBNN,NDIM,XE,SHPTOT,SHP,GRAD,DJAC)
          IF (IERR.NE.0) GOTO 9990
          IF (DJAC.LT.XZERO) iFois=iFois+1
          DJAC=ABS(DJAC)
C =======
C  2.4.2 - Erreur si le jacobien est nul en ce point de Gauss
C =======
          IF (DJAC.LT.XPETIT) THEN
            INTERR(1) = iElt
            CALL ERREUR(259)
            GOTO 9990
          ENDIF
C =======
C  2.4.3 - Recuperation de la ou des valeurs de conductibilite au point
C          de Gauss iGau (tableau VALMAT)
C =======
          DJAC=DJAC*POIGAU(iGau)
          DO i = 1, NVAMAT
            IF (IVAL(i).NE.0) THEN
              MELVAL=IVAL(i)
              IBMN = MIN(iElt,VELCHE(/2))
              IGMN = MIN(iGau,VELCHE(/1))
              VALMAT(i) = VELCHE(IGMN,IBMN)
            ELSE
              VALMAT(i) = XZERO
            ENDIF
          ENDDO
C =======
C  2.4.4 - Cas d'un materiau ISOTROPE de conductibilite K
C          Calcul de la contribution du point de Gauss a la matrice
C          CONDUCTIVITE elementaire de cet element fini
C          Ajout du terme XK*transposee(GRAD)*GRAD
C          Seul cas valide en dimension 1
C =======
          IF (IMATE.EQ.1) THEN
            XK = VALMAT(1)*DJAC
            CALL NTNST(GRAD,XK,NBNN,NDIM,CEL)
C =======
C  2.4.5 - Cas d'un materiau ORTHOTROPE ou ANISOTROPE
C          Calcul de la contribution du point de Gauss a la matrice
C          CONDUCTIVITE elementaire de cet element fini
C =======
          ELSE
C  2.4.5.1 - Calcul de la matrice de conductibilite (CMAT1) dans le
C            repere de "tropie" et de la matrice des cosinus directeurs
C            (XLOC) des directions de "tropie" dans le repere LOCAL
C            Cas des elements MASSIFS 2D
            IF (IDIM.EQ.2) THEN
C= Cas d'un materiau ORTHOTROPE
              IF (IMATE.EQ.2) THEN
                CMAT1(1,1) = VALMAT(1)
                CMAT1(2,1) = XZERO
                CMAT1(1,2) = XZERO
                CMAT1(2,2) = VALMAT(2)
                XLOC(1,1) =  VALMAT(3)
                XLOC(2,1) =  VALMAT(4)
                XLOC(1,2) = -VALMAT(4)
                XLOC(2,2) =  VALMAT(3)
              ELSE IF (IMATE.EQ.3) THEN
C= Cas d'un materiau ANISOTROPE
                CMAT1(1,1) = VALMAT(1)
                CMAT1(2,1) = VALMAT(3)
                CMAT1(1,2) = VALMAT(3)
                CMAT1(2,2) = VALMAT(2)
                XLOC(1,1) =  VALMAT(4)
                XLOC(2,1) =  VALMAT(5)
                XLOC(1,2) = -VALMAT(5)
                XLOC(2,2) =  VALMAT(4)
              ENDIF
C  2.4.5.2 - Calcul de la matrice de conductibilite (CMAT1) dans le
C            repere de "tropie" et de la matrice des cosinus directeurs
C            (XLOC) des directions de "tropie" dans le repere LOCAL
C            Cas des elements MASSIFS 3D
            ELSE
C= Cas d'un materiau ORTHOTROPE
              IF (IMATE.EQ.2) THEN
                CMAT1(1,1) = VALMAT(1)
                CMAT1(2,1) = XZERO
                CMAT1(3,1) = XZERO
                CMAT1(1,2) = XZERO
                CMAT1(2,2) = VALMAT(2)
                CMAT1(3,2) = XZERO
                CMAT1(1,3) = XZERO
                CMAT1(2,3) = XZERO
                CMAT1(3,3) = VALMAT(3)
                XLOC(1,1) = VALMAT(4)
                XLOC(2,1) = VALMAT(5)
                XLOC(3,1) = VALMAT(6)
                XLOC(1,2) = VALMAT(7)
                XLOC(2,2) = VALMAT(8)
                XLOC(3,2) = VALMAT(9)
C= Cas d'un materiau ANISOTROPE
              ELSE IF (IMATE.EQ.3) THEN
                CMAT1(1,1) = VALMAT(1)
                CMAT1(2,1) = VALMAT(4)
                CMAT1(3,1) = VALMAT(5)
                CMAT1(1,2) = VALMAT(4)
                CMAT1(2,2) = VALMAT(2)
                CMAT1(3,2) = VALMAT(6)
                CMAT1(1,3) = VALMAT(5)
                CMAT1(2,3) = VALMAT(6)
                CMAT1(3,3) = VALMAT(3)
                XLOC(1,1) = VALMAT(7)
                XLOC(2,1) = VALMAT(8)
                XLOC(3,1) = VALMAT(9)
                XLOC(1,2) = VALMAT(10)
                XLOC(2,2) = VALMAT(11)
                XLOC(3,2) = VALMAT(12)
              ENDIF
C= Calcul de la troisieme direction de "tropie"
              CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IOK)
            ENDIF
C  2.4.5.3 - Calcul de la matrice (XGLOB) des cosinus directeurs des
C            axes de "tropie" par rapport au repere GLOBAL
C*OF        CALL ZERO(XGLOB,IDIM,IDIM)
            DO j=1,IDIM
              DO k=1,IDIM
                Cte = XZERO
                DO i=1,IDIM
                  Cte=Cte+TXR(j,i)*XLOC(i,k)
                ENDDO
                XGLOB(k,j) = Cte
              ENDDO
            ENDDO
C  2.4.5.4 - Calcul de la matrice de conductibilite (CMAT) dans le
C            repere GLOBAL (traitement particulier du mode 'FOURIER')
            CALL ZERO(CMAT,NDIM,NDIM)
            IF (IFOMOD.EQ.1) THEN
              CALL ZERO(CMAT2,IDIM,IDIM)
              CALL PRODT3(CMAT2,CMAT1,XGLOB,IDIM,IDIM)
              CMAT(1,1)=CMAT2(1,1)
              IF (IMATE.EQ.2) THEN
                CMAT(2,2)=VALMAT(5)
              ELSE
                CMAT(2,2)=VALMAT(6)
              ENDIF
              CMAT(1,3)=CMAT2(1,2)
              CMAT(3,1)=CMAT2(1,2)
              CMAT(3,3)=CMAT2(2,2)
            ELSE
              CALL PRODT3(CMAT,CMAT1,XGLOB,IDIM,IDIM)
            ENDIF
C  2.4.5.5 - Contribution a la matrice de CONDUCTIVITE de cet element
C            Ajout du terme transposee(GRAD)*CMAT*GRAD
            CALL BDBST(GRAD,DJAC,CMAT,NBNN,NDIM,CEL)
          ENDIF
        ENDDO
C =====
C  2.5 - Erreur si, en un point de Gauss, le jacobien change de signe.
C =====
        IF (iFois.NE.0.AND.iFois.NE.NBPGAU) THEN
          INTERR(1) = iElt
          CALL ERREUR(195)
          GOTO 9990
        ENDIF
C =====
C  2.6 - Stockage de la matrice de CONDUCTIVITE pour cet element fini
C =====  (remplissage de XMATRI)
        CALL REMPMT(CEL,NLIGR,RE(1,1,ielt))
C
      ENDDO
 
C  3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS
C ====================================================
 9990   CONTINUE
      SEGSUP,MMAT1
      IF (IMATE.EQ.2.OR.IMATE.EQ.3) THEN
        SEGDES,MINTE1
        SEGSUP,MAXE
      ENDIF
 
      RETURN
      END
 
 
 
 
 
 
 

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