Télécharger thetor.eso

Retour à la liste

thetor
C THETOR    SOURCE    BP208322  17/03/01    21:18:30     9325           
 
C=======================================================================
C=                            T H E T O R                              =
C=                            -----------                              =
C=                                                                     =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=   Calcul des contraintes d'origine THERMIQUE dans le cas des mate-  =
C=   riaux ORTHOTROPES. Sous-programme appele par THETAP (thetap.eso). =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   MFR     (E)   Numero de la FORMULATION utilisee                   =
C=   MELE    (E)   Numero de l'element fini dans NOMTP (cf. CCHAMP)    =
C=   VALMAT  (E)   Tableau des caracteristiques du materiau            =
C=   LHOOK   (E)   Dimension de la matrice de Hooke                    =
C=   NSTRS   (E)   Nombre de composantes de contraintes                =
C=   TEMP    (E)   Temperature au point de Gauss considere             =
C=   THIF    (E)   | Temperatures des differentes couches              =
C=   THM     (E)   | dans les cas des elements coques                  =
C=   THSU    (E)   | (couches INFErieure, MOYEnne, SUPErieure)         =
C=   E3      (E)   Excentrement dans la cas d'elements coques          =
C=   SD      (E)   Section dans le cas d'elements barres               =
C=   EPAIST  (E)   Epaisseur dans le cas d'elements coques             =
C=   RES     (S)   Tableau des valeurs des deformation thermiques      =
C=   KERRE   (S)   Indicateur d'erreur si non nul                      =
C=   DHOOK   (S)   Matrice de Hooke au point de Gauss considere        =
C=   TXR,XLOC,XGLOB,ROTS   (E/S)   Tableaux de travail                 =
C=======================================================================
 
      SUBROUTINE THETOR (MFR,MELE,VALMAT,LHOOK,NSTRS,TEMP,THIF,THM,
     .                   THSU,E3,SD,EPAIST,TXR,XLOC,XGLOB,ROTS,DHOOK,
     .                   RES,KERRE)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC CCHAMP
-INC CCREEL
 
      DIMENSION XLOC(3,3),XGLOB(3,3),ROTS(NSTRS,*),DHOOK(LHOOK,*),
     .          TXR(IDIM,*),RES(NSTRS),VALMAT(*)
 
      PARAMETER (UN=1.,DEUX=2.)
 
      KERRE=0
 
C  1 - Elements MASSIFS et POREUX
C ================================
      IF (MFR.EQ.1.OR.MFR.EQ.33.OR.MFR.EQ.31) THEN
C =====
C  1.1 - Bidimensionnel PLAN Contraintes Planes
C =====
        IF (IFOUR.EQ.-2)THEN
          YG1=VALMAT(1)
          YG2=VALMAT(2)
          XNU12=VALMAT(3)
          ALP1=VALMAT(4)
          ALP2=VALMAT(5)
          XLOC(1,1)=VALMAT(6)
          XLOC(2,1)=VALMAT(7)
          XLOC(1,2)=-XLOC(2,1)
          XLOC(2,2)=XLOC(1,1)
          DO K=1,IDIM
            DO J=1,IDIM
              AA=XZero
              DO I=1,IDIM
                AA=AA+TXR(J,I)*XLOC(I,K)
              ENDDO
              XGLOB(J,K)=AA
            ENDDO
          ENDDO
          CC=XGLOB(1,1)*XGLOB(1,1)
          SS=XGLOB(2,1)*XGLOB(2,1)
          CS=XGLOB(1,1)*XGLOB(2,1)
          AA=TEMP/(UN-(XNU12*XNU12*YG2/YG1))
          st1=(YG1*ALP1+XNU12*YG2*ALP2)*AA
          st2=(YG2*ALP2+XNU12*YG2*ALP1)*AA
          RES(1)=CC*st1+SS*st2
          RES(2)=CC*st2+SS*st1
          RES(3)=XZero
          RES(4)=CS*(st1-st2)
C =====
C  1.2 - Bidimensionnel PLAN (DP/DPGE) et Axisymetrie
C =====
        ELSE IF (IFOUR.EQ.-1.OR.IFOUR.EQ.-3.OR.IFOUR.EQ.0) THEN
          YG1=VALMAT(1)
          YG2=VALMAT(2)
          YG3=VALMAT(3)
          XNU12=VALMAT(4)
          XNU23=VALMAT(5)
          XNU13=VALMAT(6)
          ALP1=VALMAT(7)
          ALP2=VALMAT(8)
          ALP3=VALMAT(9)
          XLOC(1,1)=VALMAT(10)
          XLOC(2,1)=VALMAT(11)
          XLOC(1,2)=-XLOC(2,1)
          XLOC(2,2)=XLOC(1,1)
          DO K=1,IDIM
            DO J=1,IDIM
              AA=XZero
              DO I=1,IDIM
                AA=AA+TXR(J,I)*XLOC(I,K)
              ENDDO
              XGLOB(J,K)=AA
            ENDDO
          ENDDO
          CC=XGLOB(1,1)*XGLOB(1,1)
          SS=XGLOB(2,1)*XGLOB(2,1)
          CS=XGLOB(1,1)*XGLOB(2,1)
          XNU21=(YG2/YG1)*XNU12
          XNU32=(YG3/YG2)*XNU23
          XNU31=(YG3/YG1)*XNU13
          AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
     .         -DEUX*XNU21*XNU32*XNU13)
          AUX1=YG1/AUX
          AUX2=YG2/AUX
          AUX3=YG3/AUX
          DHOOK(1,1)=(UN-XNU23*XNU32)*AUX1
          DHOOK(1,2)=(XNU21+XNU31*XNU23)*AUX1
          DHOOK(2,1)=DHOOK(1,2)
          DHOOK(1,3)=(XNU31+XNU21*XNU32)*AUX1
          DHOOK(3,1)=DHOOK(1,3)
          DHOOK(2,2)=(UN-XNU13*XNU31)*AUX2
          DHOOK(2,3)=(XNU32+XNU12*XNU31)*AUX2
          DHOOK(3,2)=DHOOK(2,3)
          DHOOK(3,3)=(UN-XNU12*XNU21)*AUX3
          st1=DHOOK(1,1)*ALP1+DHOOK(1,2)*ALP2+DHOOK(1,3)*ALP3
          st2=DHOOK(2,1)*ALP1+DHOOK(2,2)*ALP2+DHOOK(2,3)*ALP3
          st3=DHOOK(3,1)*ALP1+DHOOK(3,2)*ALP2+DHOOK(3,3)*ALP3
          RES(1)=(CC*st1+SS*st2)*TEMP
          RES(2)=(CC*st2+SS*st1)*TEMP
          RES(3)=st3*TEMP
          RES(4)=CS*(st1-st2)*TEMP
          IF (IFOUR.EQ.1) THEN
            RES(5)=XZero
            RES(6)=XZero
          ENDIF
C =====
C  1.3 - Tridimensionnel ou 2D Fourier
C =====
        ELSE IF (IFOUR.EQ.2.OR.IFOUR.EQ.1) THEN
          YG1=VALMAT(1)
          YG2=VALMAT(2)
          YG3=VALMAT(3)
          XNU12=VALMAT(4)
          XNU23=VALMAT(5)
          XNU13=VALMAT(6)
          ALP1=VALMAT(7)
          ALP2=VALMAT(8)
          ALP3=VALMAT(9)
          XLOC(1,1)=VALMAT(10)
          XLOC(2,1)=VALMAT(11)
          XLOC(3,1)=VALMAT(12)
          XLOC(1,2)=VALMAT(13)
          XLOC(2,2)=VALMAT(14)
          XLOC(3,2)=VALMAT(15)
          CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR)
          DO K=1,3
            DO J=1,IDIM
              AA=XZero
              DO I=1,IDIM
                AA=AA+TXR(J,I)*XLOC(I,K)
              ENDDO
              XGLOB(J,K)=AA
              ROTS(J,K)=AA*AA
            ENDDO
          ENDDO
cbp   en 2D Fourier, vrai TXR = [TXR(2x2) [0] ; [0] 1]   
          IF (IFOUR.EQ.1) THEN
            XGLOB(3,1)=XLOC(3,1)
            XGLOB(3,2)=XLOC(3,2)
            XGLOB(3,3)=XLOC(3,3)        
          ENDIF          
          DO i=1,3
            ROTS(i,4)=2.*XGLOB(i,1)*XGLOB(i,2)
            ROTS(i,5)=2.*XGLOB(i,2)*XGLOB(i,3)
            ROTS(i,6)=2.*XGLOB(i,1)*XGLOB(i,3)
            ROTS(4,i)=XGLOB(1,i)*XGLOB(2,i)
            ROTS(5,i)=XGLOB(2,i)*XGLOB(3,i)
            ROTS(6,i)=XGLOB(1,i)*XGLOB(3,i)
          ENDDO
          DO i=4,6
            IL1=i-3
            IL2=IL1+1
            IF (IL2.GT.3) IL2=IL2-3
            DO j=4,6
              IC1=j-3
              IC2=IC1+1
              IF (IC2.GT.3) IC2=IC2-3
              ROTS(i,j)= XGLOB(IL1,IC1)*XGLOB(IL2,IC2)
     .                  +XGLOB(IL1,IC2)*XGLOB(IL2,IC1)
            ENDDO
          ENDDO
          DO i=1,6
            AA=ROTS(6,i)
            ROTS(6,i)=ROTS(5,i)
            ROTS(5,i)=AA
          ENDDO
          DO i=1,6
            AA=ROTS(i,6)
            ROTS(i,6)=ROTS(i,5)
            ROTS(i,5)=AA
          ENDDO
          XNU21=(YG2/YG1)*XNU12
          XNU32=(YG3/YG2)*XNU23
          XNU31=(YG3/YG1)*XNU13
          AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
     .         -DEUX*XNU21*XNU32*XNU13)
          AUX1=YG1/AUX
          AUX2=YG2/AUX
          AUX3=YG3/AUX
          DHOOK(1,1)=(UN-XNU23*XNU32)*AUX1
          DHOOK(1,2)=(XNU21+XNU31*XNU23)*AUX1
          DHOOK(2,1)=DHOOK(1,2)
          DHOOK(1,3)=(XNU31+XNU21*XNU32)*AUX1
          DHOOK(3,1)=DHOOK(1,3)
          DHOOK(2,2)=(UN-XNU13*XNU31)*AUX2
          DHOOK(2,3)=(XNU32+XNU12*XNU31)*AUX2
          DHOOK(3,2)=DHOOK(2,3)
          DHOOK(3,3)=(UN-XNU12*XNU21)*AUX3
          st1=DHOOK(1,1)*ALP1+DHOOK(1,2)*ALP2+DHOOK(1,3)*ALP3
          st2=DHOOK(2,1)*ALP1+DHOOK(2,2)*ALP2+DHOOK(2,3)*ALP3
          st3=DHOOK(3,1)*ALP1+DHOOK(3,2)*ALP2+DHOOK(3,3)*ALP3
          RES(1)=(st1*ROTS(1,1)+st2*ROTS(1,2)+st3*ROTS(1,3))*TEMP
          RES(2)=(st1*ROTS(2,1)+st2*ROTS(2,2)+st3*ROTS(2,3))*TEMP
          RES(3)=(st1*ROTS(3,1)+st2*ROTS(3,2)+st3*ROTS(3,3))*TEMP
          RES(4)=(st1*ROTS(4,1)+st2*ROTS(4,2)+st3*ROTS(4,3))*TEMP
          RES(5)=(st1*ROTS(5,1)+st2*ROTS(5,2)+st3*ROTS(5,3))*TEMP
          RES(6)=(st1*ROTS(6,1)+st2*ROTS(6,2)+st3*ROTS(6,3))*TEMP
C =====
C  1.4 - Unidimensionnel (1D)
C        Les axes d'orthotropie sont confondus avec le repere global.
C =====
        ELSE IF (IFOUR.GE.3.AND.IFOUR.LE.15.AND.MFR.EQ.1) THEN
          IF (IFOUR.EQ.6) THEN
            RES(1)=VALMAT(1)*VALMAT(7)*TEMP
            RES(2)=0.
            RES(3)=0.
          ELSE IF (IFOUR.EQ.4.OR.IFOUR.EQ.8) THEN
            YG1=VALMAT(1)
            YG2=VALMAT(2)
            XNU12=VALMAT(4)
            XNU21=(YG2/YG1)*XNU12
            AUX=Un/(Un-XNU12*XNU21)
            DHOOK(1,1)=YG1*AUX
            DHOOK(2,1)=XNU21*YG1*AUX
            DHOOK(1,2)=DHOOK(2,1)
            DHOOK(2,2)=YG2*AUX
            ALP1=VALMAT(7)
            ALP2=VALMAT(8)
            RES(1)=(DHOOK(1,1)*ALP1+DHOOK(1,2)*ALP2)*TEMP
            RES(2)=(DHOOK(2,1)*ALP1+DHOOK(2,2)*ALP2)*TEMP
            RES(3)=0.
          ELSE IF (IFOUR.EQ.5.OR.IFOUR.EQ.10.OR.IFOUR.EQ.13) THEN
            YG1=VALMAT(1)
            YG3=VALMAT(3)
            XNU13=VALMAT(6)
            XNU31=(YG3/YG1)*XNU13
            AUX=Un/(Un-XNU13*XNU31)
            DHOOK(1,1)=YG1*AUX
            DHOOK(3,1)=XNU31*YG1*AUX
            DHOOK(1,3)=DHOOK(3,1)
            DHOOK(3,3)=YG3*AUX
            ALP1=VALMAT(7)
            ALP3=VALMAT(8)
            RES(1)=(DHOOK(1,1)*ALP1+DHOOK(1,3)*ALP3)*TEMP
            RES(2)=0.
            RES(3)=(DHOOK(3,1)*ALP1+DHOOK(3,3)*ALP3)*TEMP
          ELSE
            YG1=VALMAT(1)
            YG2=VALMAT(2)
            YG3=VALMAT(3)
            XNU12=VALMAT(4)
            XNU23=VALMAT(5)
            XNU13=VALMAT(6)
            XNU21=(YG2/YG1)*XNU12
            XNU32=(YG3/YG2)*XNU23
            XNU31=(YG3/YG1)*XNU13
            AUX=Un/(Un-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
     .              -Deux*XNU12*XNU23*XNU31)
            AUX1=AUX*YG1
            DHOOK(1,1)=(Un-XNU23*XNU32)*AUX1
            DHOOK(1,2)=(XNU21+XNU31*XNU23)*AUX1
            DHOOK(2,1)=DHOOK(1,2)
            DHOOK(1,3)=(XNU31+XNU21*XNU32)*AUX1
            DHOOK(3,1)=DHOOK(1,3)
            AUX2=AUX*YG2
            DHOOK(2,2)=(Un-XNU13*XNU31)*AUX2
            DHOOK(2,3)=(XNU32+XNU12*XNU31)*AUX2
            DHOOK(3,2)=DHOOK(2,3)
            DHOOK(3,3)=(Un-XNU12*XNU21)*AUX*YG3
            ALP1=VALMAT(7)*TEMP
            ALP2=VALMAT(8)*TEMP
            ALP3=VALMAT(9)*TEMP
            RES(1)=DHOOK(1,1)*ALP1+DHOOK(1,2)*ALP2+DHOOK(1,3)*ALP3
            RES(2)=DHOOK(2,1)*ALP1+DHOOK(2,2)*ALP2+DHOOK(2,3)*ALP3
            RES(3)=DHOOK(3,1)*ALP1+DHOOK(3,2)*ALP2+DHOOK(3,3)*ALP3
          ENDIF
C =====
C  1.5 - Cas non prevus
C =====
        ELSE
          KERRE=19
        ENDIF
 
C  2 - Elements JOINT 3D : JOT4 et JOI4
C ======================================
      ELSE IF (MELE.EQ.87.OR.MELE.EQ.88) THEN
        IF (IFOUR.EQ.2) THEN
          ZKN=VALMAT(3)
          ALPN=VALMAT(4)
          RES(1)=XZero
          RES(2)=XZero
          RES(3)=ALPN * ZKN * TEMP
        ELSE
          KERRE=19
        ENDIF
 
C  3 - Elements COQUES MINCES sauf COQ2
C ======================================
      ELSE IF (MELE.EQ.27.OR.MELE.EQ.28.OR.MELE.EQ.49.OR.
     .         MELE.EQ.93) THEN
        IF (IFOUR.GT.0) THEN
          YG1=VALMAT(1)
          YG2=VALMAT(2)
          XNU12=VALMAT(3)
          ALP1=VALMAT(4)
          ALP2=VALMAT(5)
          V1X=VALMAT(6)
          V1Y=VALMAT(7)
          CC=V1X*V1X
          SS=V1Y*V1Y
          CS=V1X*V1Y
          X1NUNU=UN-(XNU12*XNU12 * YG2 /YG1)
          THMO=(THSU+THIF+4.D0*THM)/6.D0
          THDI=THSU-THIF
          EXCEN=XZero
          EP2=EPAIST*EPAIST/12.D0
          FACM=THMO*EPAIST
          FACF=THMO*EPAIST*EXCEN+THDI*EP2
          st1=(YG1*ALP1+XNU12*YG2*ALP2)/X1NUNU
          st2=(YG2*ALP2+XNU12*YG2*ALP1)/X1NUNU
          RES(4)=CC*st1+SS*st2
          RES(5)=CC*st2+SS*st1
          RES(6)=CS*(st1-st2)
          RES(1)=RES(4)*FACM
          RES(2)=RES(5)*FACM
          RES(3)=RES(6)*FACM
          RES(4)=RES(4)*FACF
          RES(5)=RES(5)*FACF
          RES(6)=RES(6)*FACF
        ELSE
          KERRE=19
        ENDIF
 
C  4 - Element COQ2
C ==================
      ELSE IF (MELE.EQ.44) THEN
        YG1=VALMAT(1)
        YG2=VALMAT(2)
        XNU12=VALMAT(3)
        ALP1=VALMAT(4)
        ALP2=VALMAT(5)
        V1X=VALMAT(6)
        V1Y=VALMAT(7)
        CC=V1X*V1X
        SS=V1Y*V1Y
        CS=V1X*V1Y
        X1NUNU=UN-(XNU12*XNU12 * YG2 / YG1)
        THMO=(THSU+THIF+4.D0*THM)/6.D0
        THDI=THSU-THIF
        EXCEN=XZero
        EP2=EPAIST*EPAIST/12.D0
        FACM=THMO*EPAIST
        FACF=THMO*EPAIST*EXCEN+THDI*EP2
        st1=(YG1*ALP1+XNU12*YG2*ALP2)/X1NUNU
        st2=(YG2*ALP2+XNU12*YG2*ALP1)/X1NUNU
        IF (IFOUR.GT.0) THEN
          RES(4)=CC*st1+SS*st2
          RES(5)=CC*st2+SS*st1
          RES(6)=CS*(st1-st2)
          RES(1)=RES(4)*FACM
          RES(2)=RES(5)*FACM
          RES(3)=RES(6)*FACM
          RES(4)=RES(4)*FACF
          RES(5)=RES(5)*FACF
          RES(6)=RES(6)*FACF
        ELSE IF (IFOUR.EQ.0) THEN
          RES(3)=CC*st1+SS*st2
          RES(4)=CC*st2+SS*st1
          RES(1)=RES(3)*FACM
          RES(2)=RES(4)*FACM
          RES(3)=RES(3)*FACF
          RES(4)=RES(4)*FACF
        ELSE
          KERRE=19
        ENDIF
 
C  5 - Elements COQUES EPAISSES
C ==============================
      ELSE IF (MELE.EQ.41.OR.MELE.EQ.56) THEN
        YG1=VALMAT(1)
        YG2=VALMAT(2)
        XNU12=VALMAT(3)
        ALP1=VALMAT(4)
        ALP2=VALMAT(5)
        V1X=VALMAT(6)
        V1Y=VALMAT(7)
        CC=V1X*V1X
        SS=V1Y*V1Y
        CS=V1X*V1Y
        X1NUNU=UN-(XNU12*XNU12 * YG2 / YG1)
        THG=0.5*E3*((E3+UN)*THSU+(E3-UN)*THIF)+(UN-E3*E3)*THM
        st1=(YG1*ALP1+XNU12*YG2*ALP2)/X1NUNU
        st2=(YG2*ALP2+XNU12*YG2*ALP1)/X1NUNU
        RES(1)=(CC*st1+SS*st2)*THG
        RES(2)=(CC*st2+SS*st1)*THG
        RES(3)=CS*(st1-st2)*THG
        RES(4)=XZero
        RES(5)=XZero
 
C  6 - Elements POUTRES et TUYAUX
C ================================
      ELSE IF (MELE.EQ.29.OR.MELE.EQ.42.OR.MELE.EQ.84) THEN
        YOU=VALMAT(1)
        ALPH=VALMAT(2)
        RES(1)=YOU*ALPH*SD*TEMP
 
C  7 - Elements BARRES, CERCES et POINT
C ======================================
      ELSEIF(MELE.EQ.46.OR.MELE.EQ.95.OR.MELE.EQ.45)THEN
        IF (MELE.EQ.45.AND.IFOUR.NE.-3) THEN
          KERRE=86
        ELSE IF (MELE.EQ.95.AND.IFOUR.NE.0) THEN
          KERRE=86
        ELSE
          YOU=VALMAT(1)
          ALPH=VALMAT(2)
          RES(1)=YOU*ALPH*SD*TEMP
        ENDIF
 
C  8 - Autres cas non prevus
C ===========================
      ELSE
        KERRE=19
      ENDIF
 
      RETURN
      END
 
 
 
 
 
 
 

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