1 : C TUFIKT SOURCE CHAT 05/01/13 03:55:50 5004 2 : SUBROUTINE TUFIKT(XMAT,CONT,CARA,VARI,REG,XPREC,KERRE) 3 : C======================================================================= 4 : C CALCUL DE LA MATRICE DE RIGIDITE TANGENTE DE L'ELEMENT 5 : C DE TUYAU FISSURE DANS LE REPERE GLOBAL 6 : C ENTREES: 7 : C MAT = LE MATERIAU 8 : C CONT = LES CONTRAINTES 9 : C CARA = LES CARACTERISTIQUES 10 : C VARI = LES VARIABLES INTERNES 11 : C XPREC = PRECISION POUR TESTER SI ON EST ELASTIQUE 12 : C SORTIE: 13 : C REG(12,12)= MATRICE DE RIGIDITE GLOBALE 14 : C KERRE = CODE D'ERREUR CORRESPONDANT A LA POSITION RELATIVE DES 15 : C DEUX VECTEURS V1 ET V2 DEFINISSANT L'ELEMENT TUFI: 16 : C 0 SI PAS DE PROBLEME, 1 SINON. 17 : C 1 SI V2 EST NUL OU SI V1 ET V2 SONT COLINEAIRES. 18 : C 2 SI V1 EST NUL. 19 : C 3 SI LE TUYAU FISSURE EST TRP EPAI R/T INF A 2 20 : C 4 CONTRAINTES A L EXTERIEUR DE LA SURFACE 21 : C======================================================================= 22 : IMPLICIT INTEGER(I-N) 23 : IMPLICIT REAL*8(A-H,O-Z) 24 : C Include contenant quelques constantes dont XPI : 25 : -INC CCREEL 26 : DIMENSION XMAT(*),CONT(*),CARA(*),VARI(*) 27 : DIMENSION REL(12,12),V1(3),V2(3),BPSS(3,3) 28 : DIMENSION REG(12,*),BV1(3),BV2(3),RIL(6,6),RAL(12,12) 29 : C 30 : DATA EPS/1.D-10/ 31 : DATA XZER,UN,UNS2/0.D0,1.D0,.5D0/ 32 : C 33 : C D6 = COEFFICIENT INTERVENANT POUR LES TERMES DE PENALISATION. 34 : C 35 : DATA DEUX,D6/2.D0,1.D06/ 36 : C 37 : C ON RECUPERE LES INFORMATIONS 38 : C 39 : XP=CONT(1) 40 : XM=-CONT(6) 41 : YOU=XMAT(1) 42 : XNU=XMAT(2) 43 : SIGMA0=XMAT(5) 44 : RAYON=CARA(1) 45 : EPAI=CARA(2) 46 : V1(1)=CARA(3) 47 : V1(2)=CARA(4) 48 : V1(3)=CARA(5) 49 : V2(1)=CARA(6) 50 : V2(2)=CARA(7) 51 : V2(3)=CARA(8) 52 : TETA1=VARI(3)/2.D0 53 : CISA = YOU /(DEUX *(UN + XNU)) 54 : KERRE = 0 55 : C======================================================================= 56 : C LES CARACTERISTIQUES REPRESENTENT: 57 : C EPAI = EPAISSEUR DU TUYAU 58 : C RAYON = RAYON EXTERIEUR DU TUYAU 59 : C TETA1 = DEMI-ANGLE D'OUVERTURE DE LA FISSURE EN DEGRE 60 : C V1 = VECTEUR ORIENTANT LES NOEUDS 1 ET 2 (AXE DU TUYAU) 61 : C V2 = VECTEUR ORIENTANT L'AXE DE LA FISSURE 62 : C YOU = MODULE D'YOUNG 63 : C CISA = MODULE DE CISAILLEMENT 64 : C TABLEAU DE TRAVAIL: 65 : C BPSS = MATRICE DE PASSAGE 66 : C ON OBTIENT: 67 : C REG(12,12) = MATRICE DE RIGIDITE AXES GLOBAUX 68 : C 69 : C======================================================================= 70 : C 71 : C CONVERSION DE TETA1 EN RADIAN 72 : C 73 : TETA = (TETA1 * XPI)/180.D0 74 : C 75 : C ON MET DANS 'RAYON' LE RAYON MOYEN DU TUYAU. 76 : C 77 : RAYMO = RAYON - EPAI*UNS2 78 : XM0=4.D0*SIGMA0*RAYMO*RAYMO*EPAI 79 : XP0=2.D0*SIGMA0*XPI*RAYMO*EPAI 80 : C 81 : C ON CALCULE A 82 : C 83 : RSURT= RAYMO / EPAI 84 : IF(RSURT.LE.10.0D0.AND.RSURT.GE.4.9D0) THEN 85 : AXX= (.125D0*RSURT - .25D0 ) **.25D0 86 : ELSE IF(RSURT.GT.10.D0.AND.RSURT.LE.35.D0) THEN 87 : AXX= (.4D0*RSURT - 3.D0 ) **.25D0 88 : ELSE 89 : KERRE=3 90 : ENDIF 91 : IF(KERRE.EQ.3) GOTO 666 92 : C 93 : C NORMALISATION DES VECTEURS V1 ET V2. 94 : C 95 : XNV1 = XZER 96 : XNV2 = XZER 97 : DO 5 I=1,3 98 : XNV1 = XNV1 + (V1(I)*V1(I)) 99 : XNV2 = XNV2 + (V2(I)*V2(I)) 100 : 5 CONTINUE 101 : C 102 : C MISE A ZERO DE REG 103 : C 103 : CALL ZERO(REG,12,12) 104 : IF (XNV1.LT.EPS) THEN 105 : KERRE = 2 106 : GOTO 666 107 : ENDIF 108 : IF (XNV2.LT.EPS) THEN 109 : KERRE = 1 110 : GOTO 666 111 : ENDIF 112 : XNV1 = UN/SQRT(XNV1) 113 : XNV2 = UN/SQRT(XNV2) 114 : DO 10 I=1,3 115 : BV1(I) = V1(I)*XNV1 116 : BV2(I) = V2(I)*XNV2 117 : 10 CONTINUE 118 : C 119 : C ORTHOGONALISATION ET RENORMALISATION DE V2. 120 : C 121 : PS = (BV1(1)*BV2(1))+(BV1(2)*BV2(2))+(BV1(3)*BV2(3)) 122 : C 123 : C TEST SUR LA COLINEARITE DE V1 ET V2. 124 : C 125 : IF(ABS(PS).GE.(.99D0)) THEN 126 : KERRE = 1 127 : GOTO 666 128 : ELSE 129 : XNV2 = UN/SQRT(UN-(PS**2)) 130 : DO 15 I=1,3 131 : BV2(I) = XNV2*(BV2(I)-(PS*BV1(I))) 132 : 15 CONTINUE 133 : ENDIF 134 : C 135 : C CALCUL DE LA MATRICE DE PASSAGE 136 : C LOCAL = BPSS * GLOBAL 137 : C BV1 REPRESENTE L'AXE LOCAL DES Z 138 : C BV2 REPRESENTE L'AXE LOCAL DES Y 139 : C 140 : DO 20 I=1,3 141 : BPSS(3,I) =BV1(I) 142 : BPSS(2,I) =BV2(I) 143 : 20 CONTINUE 144 : BPSS(1,1) = (BPSS(2,2)*BPSS(3,3))-(BPSS(2,3)*BPSS(3,2)) 145 : BPSS(1,2) = (BPSS(2,3)*BPSS(3,1))-(BPSS(2,1)*BPSS(3,3)) 146 : BPSS(1,3) = (BPSS(2,1)*BPSS(3,2))-(BPSS(2,2)*BPSS(3,1)) 147 : C 148 : C CALCUL DES TERMES DE LA MATRICE DE RIGIDITE REPERE LOCAL 149 : C 150 : TESPI = TETA/XPI 150 : CALL ZERO(RAL,6,6) 151 : C 152 : C ON TESTE SUR L'IMPORTANCE DE L'ANGLE D'OUVERTURE DE LA FISSURE. SI 153 : C CE DERNIER EST INFERIEUR A UN DEGRE, ON ANNULE TOUS LES COEFFI- 154 : C CIENTS DE COUPLAGE ET ON PENALISE TOUS LES DEGRES DE LIBERTE. DANS 155 : C LE CAS CONTRAIRE, SEULS SONT PERMIS LA ROTATION AUTOUR DE L'AXE 156 : C LOCAL DES X ET L'ALLONGEMENT SUIVANT L'AXE LOCAL DES Z ,CES DEUX 157 : C DEGRES DE LIBERTE ETANT COUPLES. 158 : C 159 : IF (TETA1.LE.(0.5D0)) THEN 160 : RAL(3,3) = DEUX * YOU * XPI * EPAI * RAYMO * D6 161 : RAL(4,4) = YOU * XPI * EPAI * (RAYMO**2) * D6 162 : ELSE 163 : C 163 : CALL TUFIFP(TESPI,AXX,FP,FM,FMP,FOP,FOM) 164 : DENOM = FP*FM - FMP*FMP 165 : COEF = (YOU * XPI * EPAI)/(DENOM * (TETA**2)) 166 : RAL(3,3) = COEF * DEUX * FM 167 : RAL(3,4) = -COEF * RAYMO * FMP 168 : RAL(4,3) = RAL(3,4) 169 : RAL(4,4) = COEF * (RAYMO**2) * FP * UNS2 170 : ENDIF 171 : C 172 : C INITIALISATION DES TERMES DE PENALISATION 173 : C 174 : RAL(1,1) = DEUX * CISA * XPI * EPAI * RAYMO * D6 175 : RAL(2,2) = RAL(1,1) 176 : RAL(5,5) = YOU * XPI * EPAI * (RAYMO**2) * D6 177 : RAL(6,6) = DEUX * CISA * XPI * EPAI * (RAYMO**2) * D6 178 : C 179 : C TEST POUR SAVOIR SI L ETAT DE CONTRAINTES EST ELASTIQUE 180 : C SI OUI ON RENVOIT LA RAIDEUR INITIALE 181 : C SI ON EST SUR LA SURFACE ON RENVOIT LA RAIDEUR TANGENTE 182 : C SI ON EST A L EXTERIEUR ON MET KERRE A 4 183 : C 184 : TES=CRIT1(XM,XP,TETA,XM0,XP0) 185 : IF(TES.LT.(-XPREC)) THEN 186 : DO 40 I=1,6 187 : DO 401 J=1,6 188 : RAL(I+6,J+6) = RAL(I,J) 189 : RAL(I+6,J ) = -RAL(I,J) 190 : RAL(I ,J+6) = -RAL(I,J) 191 : 401 CONTINUE 192 : 40 CONTINUE 193 : C 193 : CALL POUROT(REG,12,BPSS,RAL) 194 : GOTO 666 195 : ENDIF 196 : IF(TES.GT.XPREC) THEN 197 : KERRE=4 198 : GOTO 666 199 : ENDIF 200 : C ON CALCULE LA MATRICE TANGENTE EN MULTIPLIANT LA 201 : C MATRICE DE RAIDEUR PAR 202 : C T T 203 : C I-((DF/DSIG)(DF/DSIG)K)/((DF/DSIG)K(DF/DSIG)) 204 : C 205 : C 206 : C ON TESTE POUR SAVOIR SUR QUELLE SURFACE ON TRAVAILLE 207 : C 208 : CRI=2.D0*XP0*(1.D0-TETA/XPI)*XM 209 : CRI=CRI+XM0*SIN(TETA)*XP 210 : IF(CRI.GE.XZER) THEN 211 : DFSDM=1.D0/XM0 212 : DFSDP=XPI*SIN(0.5D0*(TETA+XPI*XP/XP0))/2.D0/XP0 213 : ELSE 214 : DFSDM=-1.D0/XM0 215 : DFSDP=-XPI*SIN(0.5D0*(TETA-XPI*XP/XP0))/2.D0/XP0 216 : ENDIF 217 : DM2=DFSDM*DFSDM 218 : DP2=DFSDP*DFSDP 219 : DMP=DFSDM*DFSDP 220 : C 221 : C ON FORME LE TERME CORRECTEUR 222 : C 223 : DEN=RAL(3,3)*DM2 224 : DEN=DEN+2.D0*RAL(3,4)*DMP 225 : DEN=DEN+RAL(4,4)*DP2 225 : CALL ZERO(RIL,6,6) 226 : RIL(3,3)=1.D0-(RAL(3,3)*DM2+RAL(3,4)*DMP)/DEN 227 : RIL(3,4)=-(RAL(3,4)*DM2+RAL(4,4)*DMP)/DEN 228 : RIL(4,3)=-(RAL(3,3)*DMP+RAL(3,4)*DP2)/DEN 229 : RIL(4,4)=1.D0-(RAL(3,4)*DMP+RAL(4,4)*DP2)/DEN 230 : RIL(1,1)=1.D0 231 : RIL(2,2)=1.D0 232 : RIL(5,5)=1.D0 233 : RIL(6,6)=1.D0 233 : CALL MULMAT(REL,RAL,RIL,6,6,6) 234 : C 235 : C CALCUL DE LA MATRICE (12,12) A PARTIR DE LA MATRICE (6,6) 236 : C 237 : DO 30 I=1,6 238 : DO 301 J=1,6 239 : REL(I+6,J+6) = REL(I,J) 240 : REL(I+6,J ) = -REL(I,J) 241 : REL(I ,J+6) = -REL(I,J) 242 : 301 CONTINUE 243 : 30 CONTINUE 244 : C 244 : CALL POUROT(REG,12,BPSS,REL) 245 : 666 CONTINUE 246 : RETURN 247 : END 248 : 249 : 250 : 251 : 252 :
© Cast3M 2003 - All rights reserved.
Disclaimer