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