Sources Esope | Cast3M

tufikt
C TUFIKT    SOURCE    CHAT      05/01/13    03:55:50     5004
      SUBROUTINE TUFIKT(XMAT,CONT,CARA,VARI,REG,XPREC,KERRE)
C=======================================================================
C   CALCUL DE LA MATRICE DE RIGIDITE TANGENTE DE L'ELEMENT
C    DE TUYAU FISSURE DANS LE REPERE GLOBAL
C   ENTREES:
C     MAT   = LE MATERIAU
C     CONT  = LES CONTRAINTES
C     CARA  = LES CARACTERISTIQUES
C     VARI  = LES VARIABLES INTERNES
C     XPREC = PRECISION POUR TESTER SI ON EST ELASTIQUE
C   SORTIE:
C     REG(12,12)=  MATRICE DE RIGIDITE  GLOBALE
C     KERRE = CODE D'ERREUR CORRESPONDANT A LA POSITION RELATIVE DES
C             DEUX VECTEURS V1 ET V2 DEFINISSANT L'ELEMENT TUFI:
C             0 SI PAS DE PROBLEME, 1 SINON.
C             1 SI V2 EST NUL OU SI V1 ET V2 SONT COLINEAIRES.
C             2 SI V1 EST NUL.
C             3 SI LE TUYAU FISSURE EST TRP EPAI R/T  INF A 2
C             4 CONTRAINTES A L EXTERIEUR DE LA SURFACE
C=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
      DIMENSION XMAT(*),CONT(*),CARA(*),VARI(*)
      DIMENSION REL(12,12),V1(3),V2(3),BPSS(3,3)
      DIMENSION REG(12,*),BV1(3),BV2(3),RIL(6,6),RAL(12,12)
C
      DATA EPS/1.D-10/
      DATA XZER,UN,UNS2/0.D0,1.D0,.5D0/
C
C   D6 = COEFFICIENT INTERVENANT POUR LES TERMES DE PENALISATION.
C
      DATA DEUX,D6/2.D0,1.D06/
C
C     ON RECUPERE LES INFORMATIONS
C
      XP=CONT(1)
      XM=-CONT(6)
      YOU=XMAT(1)
      XNU=XMAT(2)
      SIGMA0=XMAT(5)
      RAYON=CARA(1)
      EPAI=CARA(2)
      V1(1)=CARA(3)
      V1(2)=CARA(4)
      V1(3)=CARA(5)
      V2(1)=CARA(6)
      V2(2)=CARA(7)
      V2(3)=CARA(8)
      TETA1=VARI(3)/2.D0
      CISA  = YOU /(DEUX *(UN + XNU))
      KERRE = 0
C=======================================================================
C   LES CARACTERISTIQUES REPRESENTENT:
C      EPAI   = EPAISSEUR DU TUYAU
C      RAYON  = RAYON EXTERIEUR DU TUYAU
C      TETA1  = DEMI-ANGLE D'OUVERTURE DE LA FISSURE EN DEGRE
C      V1     = VECTEUR ORIENTANT LES NOEUDS 1 ET 2 (AXE DU TUYAU)
C      V2     = VECTEUR ORIENTANT L'AXE DE LA FISSURE
C      YOU    = MODULE D'YOUNG
C      CISA   = MODULE DE CISAILLEMENT
C   TABLEAU DE TRAVAIL:
C      BPSS   = MATRICE DE PASSAGE
C   ON OBTIENT:
C      REG(12,12) = MATRICE DE RIGIDITE  AXES GLOBAUX
C
C=======================================================================
C
C     CONVERSION DE TETA1 EN RADIAN
C
      TETA = (TETA1 * XPI)/180.D0
C
C     ON MET DANS 'RAYON' LE RAYON MOYEN DU TUYAU.
C
      RAYMO = RAYON - EPAI*UNS2
      XM0=4.D0*SIGMA0*RAYMO*RAYMO*EPAI
      XP0=2.D0*SIGMA0*XPI*RAYMO*EPAI
C
C     ON CALCULE A
C
      RSURT= RAYMO / EPAI
      IF(RSURT.LE.10.0D0.AND.RSURT.GE.4.9D0) THEN
         AXX= (.125D0*RSURT - .25D0 ) **.25D0
      ELSE IF(RSURT.GT.10.D0.AND.RSURT.LE.35.D0) THEN
         AXX= (.4D0*RSURT - 3.D0 ) **.25D0
      ELSE
         KERRE=3
      ENDIF
      IF(KERRE.EQ.3) GOTO 666
C
C     NORMALISATION DES VECTEURS V1 ET V2.
C
      XNV1 = XZER
      XNV2 = XZER
      DO 5 I=1,3
         XNV1 = XNV1 + (V1(I)*V1(I))
         XNV2 = XNV2 + (V2(I)*V2(I))
 5    CONTINUE
C
C     MISE A ZERO DE REG
C
      CALL ZERO(REG,12,12)
      IF (XNV1.LT.EPS) THEN
         KERRE = 2
         GOTO 666
      ENDIF
      IF (XNV2.LT.EPS) THEN
         KERRE = 1
         GOTO 666
      ENDIF
      XNV1 = UN/SQRT(XNV1)
      XNV2 = UN/SQRT(XNV2)
      DO 10 I=1,3
         BV1(I) = V1(I)*XNV1
         BV2(I) = V2(I)*XNV2
 10   CONTINUE
C
C     ORTHOGONALISATION ET RENORMALISATION DE V2.
C
      PS = (BV1(1)*BV2(1))+(BV1(2)*BV2(2))+(BV1(3)*BV2(3))
C
C     TEST SUR LA COLINEARITE DE V1 ET V2.
C
      IF(ABS(PS).GE.(.99D0)) THEN
         KERRE = 1
         GOTO 666
      ELSE
         XNV2 = UN/SQRT(UN-(PS**2))
         DO 15 I=1,3
            BV2(I) = XNV2*(BV2(I)-(PS*BV1(I)))
 15      CONTINUE
      ENDIF
C
C     CALCUL DE LA MATRICE DE PASSAGE
C      LOCAL = BPSS * GLOBAL
C    BV1 REPRESENTE L'AXE LOCAL DES Z
C    BV2 REPRESENTE L'AXE LOCAL DES Y
C
      DO 20 I=1,3
         BPSS(3,I) =BV1(I)
         BPSS(2,I) =BV2(I)
 20   CONTINUE
      BPSS(1,1) = (BPSS(2,2)*BPSS(3,3))-(BPSS(2,3)*BPSS(3,2))
      BPSS(1,2) = (BPSS(2,3)*BPSS(3,1))-(BPSS(2,1)*BPSS(3,3))
      BPSS(1,3) = (BPSS(2,1)*BPSS(3,2))-(BPSS(2,2)*BPSS(3,1))
C
C     CALCUL DES TERMES DE LA MATRICE DE RIGIDITE     REPERE LOCAL
C
      TESPI = TETA/XPI
      CALL ZERO(RAL,6,6)
C
C     ON TESTE SUR L'IMPORTANCE DE L'ANGLE D'OUVERTURE DE LA FISSURE. SI
C     CE DERNIER EST INFERIEUR A UN DEGRE, ON ANNULE TOUS LES COEFFI-
C     CIENTS DE COUPLAGE ET ON PENALISE TOUS LES DEGRES DE LIBERTE. DANS
C     LE CAS CONTRAIRE, SEULS SONT PERMIS LA ROTATION AUTOUR DE L'AXE
C     LOCAL DES X ET L'ALLONGEMENT SUIVANT L'AXE LOCAL DES Z ,CES DEUX
C     DEGRES DE LIBERTE ETANT COUPLES.
C
      IF (TETA1.LE.(0.5D0)) THEN
         RAL(3,3) = DEUX * YOU * XPI * EPAI * RAYMO      * D6
         RAL(4,4) =        YOU * XPI * EPAI * (RAYMO**2) * D6
      ELSE
C
         CALL TUFIFP(TESPI,AXX,FP,FM,FMP,FOP,FOM)
         DENOM = FP*FM - FMP*FMP
         COEF = (YOU * XPI * EPAI)/(DENOM * (TETA**2))
         RAL(3,3) =  COEF * DEUX * FM
         RAL(3,4) = -COEF * RAYMO * FMP
         RAL(4,3) =  RAL(3,4)
         RAL(4,4) =  COEF * (RAYMO**2) * FP * UNS2
      ENDIF
C
C     INITIALISATION DES TERMES DE PENALISATION
C
      RAL(1,1) = DEUX * CISA * XPI * EPAI * RAYMO      * D6
      RAL(2,2) = RAL(1,1)
      RAL(5,5) =         YOU * XPI * EPAI * (RAYMO**2) * D6
      RAL(6,6) = DEUX * CISA * XPI * EPAI * (RAYMO**2) * D6
C
C   TEST POUR SAVOIR SI L ETAT DE CONTRAINTES EST ELASTIQUE
C    SI OUI ON RENVOIT LA RAIDEUR INITIALE
C    SI ON EST SUR LA SURFACE ON RENVOIT LA RAIDEUR TANGENTE
C    SI ON EST A L EXTERIEUR ON MET KERRE A 4
C
      TES=CRIT1(XM,XP,TETA,XM0,XP0)
      IF(TES.LT.(-XPREC)) THEN
         DO 40 I=1,6
            DO 401 J=1,6
               RAL(I+6,J+6) =  RAL(I,J)
               RAL(I+6,J  ) = -RAL(I,J)
               RAL(I  ,J+6) = -RAL(I,J)
 401        CONTINUE
 40      CONTINUE
C
         CALL POUROT(REG,12,BPSS,RAL)
         GOTO 666
      ENDIF
      IF(TES.GT.XPREC) THEN
         KERRE=4
         GOTO 666
      ENDIF
C     ON CALCULE LA MATRICE TANGENTE EN MULTIPLIANT LA
C     MATRICE DE RAIDEUR PAR
C                  T                               T
C       I-((DF/DSIG)(DF/DSIG)K)/((DF/DSIG)K(DF/DSIG))
C
C
C     ON TESTE POUR SAVOIR SUR QUELLE SURFACE ON TRAVAILLE
C
      CRI=2.D0*XP0*(1.D0-TETA/XPI)*XM
      CRI=CRI+XM0*SIN(TETA)*XP
      IF(CRI.GE.XZER) THEN
         DFSDM=1.D0/XM0
         DFSDP=XPI*SIN(0.5D0*(TETA+XPI*XP/XP0))/2.D0/XP0
      ELSE
         DFSDM=-1.D0/XM0
         DFSDP=-XPI*SIN(0.5D0*(TETA-XPI*XP/XP0))/2.D0/XP0
      ENDIF
      DM2=DFSDM*DFSDM
      DP2=DFSDP*DFSDP
      DMP=DFSDM*DFSDP
C
C     ON FORME LE TERME CORRECTEUR
C
      DEN=RAL(3,3)*DM2
      DEN=DEN+2.D0*RAL(3,4)*DMP
      DEN=DEN+RAL(4,4)*DP2
      CALL ZERO(RIL,6,6)
      RIL(3,3)=1.D0-(RAL(3,3)*DM2+RAL(3,4)*DMP)/DEN
      RIL(3,4)=-(RAL(3,4)*DM2+RAL(4,4)*DMP)/DEN
      RIL(4,3)=-(RAL(3,3)*DMP+RAL(3,4)*DP2)/DEN
      RIL(4,4)=1.D0-(RAL(3,4)*DMP+RAL(4,4)*DP2)/DEN
      RIL(1,1)=1.D0
      RIL(2,2)=1.D0
      RIL(5,5)=1.D0
      RIL(6,6)=1.D0
      CALL MULMAT(REL,RAL,RIL,6,6,6)
C
C     CALCUL DE LA MATRICE (12,12) A PARTIR DE LA MATRICE (6,6)
C
      DO 30 I=1,6
         DO 301 J=1,6
            REL(I+6,J+6) =  REL(I,J)
            REL(I+6,J  ) = -REL(I,J)
            REL(I  ,J+6) = -REL(I,J)
 301     CONTINUE
 30   CONTINUE
C
      CALL POUROT(REG,12,BPSS,REL)
 666  CONTINUE
      RETURN
      END