1 : C BETOKT SOURCE PV 11/03/07 21:15:10 6885 2 : SUBROUTINE BETOKT(XMAT,TRAC,VARF,SIGF,NCURVT,NCURVC,XKTAT,KERRE) 3 : C 4 : C Matrice tangente pour le modele endommageable en cyclique 5 : C D. Combescure 18/09/95 ELSA-Ispra 6 : C 7 : C ENTREES: 8 : C ------- 9 : C XMAT = Caracteristiques materiaux 10 : C 11 : C VARF = Variables internes 12 : C 13 : C SIGF = Contraintes 14 : C 15 : C SORTIES: 16 : C ------- 17 : C XKTAT = Matrice tangente 18 : C 19 : C================================================================ 20 : IMPLICIT INTEGER(I-N) 21 : IMPLICIT REAL*8(A-H,O-Z) 22 : -INC CCOPTIO 23 : DIMENSION XMAT(*) 24 : DIMENSION TRAC(*) 25 : C 26 : DIMENSION VARF(14),SIGF(4),EPSI(2),FG(2) 27 : DIMENSION XNOR1(4),XNOR2(4),XCOU1(3,3),XCOU2(3,3) 28 : DIMENSION GAMBDA(2),EEE(2,2),SIGP(2),EEE2(2,2) 29 : DIMENSION UUU(3,2),VVV(2,3),UUU2(3,2),VVV2(2,3),UUU3(3,2) 30 : DIMENSION HHH1(3,3),HHH2(3,3),HHH3(3,3) 31 : DIMENSION XKTAT(4,4),RESIST(6),XKISO(2) 32 : PARAMETER (UN=1.D0,UNDEMI=0.5D0,EPSILO=1.D-20,XZER=0.D0) 33 : PARAMETER (PRECIS=1.D-4) 34 : KERRE=0 35 : C================================================================ 36 : C LECTURE DES CARACTERISTIQUES MAT. ET DES VAR. INT. 37 : C=============================================================== 38 : C Pour la lecture des courbes 39 : IPT = 1 40 : IPC = IPT+2*NCURVT 41 : C 42 : YOUN = XMAT(1) 43 : XNU = XMAT(2) 44 : XGGG = (UNDEMI*YOUN/(UN+XNU)) 45 : XHH1 = XMAT(3) 46 : DEGRAD = XMAT(8)/XMAT(9) 47 : RESIST(1) = XMAT(4) 48 : RESIST(2) = XMAT(4) 49 : RESIST(3) = XMAT(5) 50 : RESIST(4) = XMAT(5) 51 : RESIST(5) = XMAT(6) 52 : RESIST(6) = XMAT(7) 53 : XLCAT = XMAT(10) 54 : XLCAC = XMAT(11) 55 : CPHI = VARF(5) + UN 56 : SPHI = VARF(6) 57 : EPSI(1) = VARF(3) 58 : EPSI(2) = VARF(4) 59 : ICAS1 = nint(VARF(11)) 60 : ICAS2 = nint(VARF(12)) 61 : C============================================================== 62 : C RAIDEUR ELASTIQUE 63 : C=============================================================== 63 : CALL ZERO(XKTAT,4,4) 64 : XKTAT(1,1) = YOUN/(UN - (XNU**2)) 65 : XKTAT(1,2) = XNU*XKTAT(1,1) 66 : XKTAT(2,1) = XKTAT(1,2) 67 : XKTAT(2,2) = XKTAT(1,1) 68 : XKTAT(4,4) = YOUN/(2.D0*(UN+XNU)) 69 : C============================================================== 70 : IF (ICAS1.EQ.0) THEN 71 : GOTO 9999 72 : ENDIF 73 : C 73 : CALL VALEUR(EPSI(1),TRAC(IPT),NCURVT,XLCAT, 74 : &XKISO(1),XPEN1,KERRE) 74 : CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC, 75 : &XKISO(2),XPEN2,KERRE) 76 : GAMBDA(1) = VARF(13) 77 : GAMBDA(2) = VARF(14) 78 : IF (ICAS1.LE.4) THEN 79 : FG(1) = RESIST(ICAS1) 80 : ELSE 81 : FG(1) = XKISO(ICAS1-4)*RESIST(ICAS1) 82 : ENDIF 83 : IF (ICAS2.LE.4) THEN 84 : FG(2) = RESIST(ICAS2) 85 : ELSE 86 : FG(2) = XKISO(ICAS2-4)*RESIST(ICAS2) 87 : ENDIF 88 : C============================================================== 89 : C CALCUL DE L'INVERSE 90 : C============================================================== 90 : CALL ZERO(HHH1,3,3) 91 : HHH1(1,1) = UN/YOUN 92 : HHH1(1,2) = -XNU*HHH1(1,1) 93 : HHH1(2,1) = HHH1(1,2) 94 : HHH1(2,2) = HHH1(1,1) 95 : HHH1(3,3) = (2.D0*(UN+XNU)/YOUN) 96 : C 96 : CALL CONPRI (SIGF,CPHI,SPHI,SIGP,CPHI,SPHI) 97 : C 98 : IF (ICAS1.EQ.1) THEN 99 : ISGNTC = 1 100 : ELSE IF (ICAS1.EQ.2) THEN 101 : ISGNTC = 1 102 : ELSE IF (ICAS1.EQ.3) THEN 103 : ISGNTC = -1 104 : ELSE IF (ICAS1.EQ.4) THEN 105 : ISGNTC = -1 106 : ELSE IF (ICAS1.EQ.5) THEN 107 : ISGNTC = 1 108 : ELSE IF (ICAS1.EQ.6) THEN 109 : ISGNTC = -1 110 : ENDIF 111 : IF (ICAS1.LE.4) THEN 112 : IF (SIGP(2).GE.SIGP(1)) THEN 113 : IF ((ICAS1.EQ.1).OR.(ICAS1.EQ.3)) THEN 114 : ISGN2 = -1 115 : ELSE 116 : ISGN2 = 1 117 : ENDIF 118 : ELSE 119 : IF ((ICAS1.EQ.1).OR.(ICAS1.EQ.3)) THEN 120 : ISGN2 = 1 121 : ELSE 122 : ISGN2 = -1 123 : ENDIF 124 : ENDIF 125 : ELSE 126 : ISGN2 = 1 127 : ENDIF 128 : IF (ICAS1.LE.4) THEN 128 : CALL COURBU(SIGF,FG(1),ISGNTC,ISGN2,1,XCOU1) 128 : CALL NORMBL(SIGF,FG(1),ISGNTC,ISGN2,1,XNOR1) 129 : ELSE 130 : PHI = SQRT ((UNDEMI*(SIGF(1)-SIGF(2)))**2+SIGF(4)**2) 131 : IF ((PHI.LT.ABS(PRECIS*FG(1))).AND.(ICAS1.EQ.ICAS2)) THEN 132 : C Cas de l'apex ou on calcule directement H-1 133 : HHH1(1,1) = XGGG*((4.D0*PHI)+(YOUN*GAMBDA(1)))/ 134 : & (2.D0*(UN-XNU)*((2.D0*PHI)+(XGGG*GAMBDA(1)))) 135 : HHH1(2,2)=HHH1(1,1) 136 : HHH1(1,2) = XGGG*((4.D0*XNU*PHI)+(YOUN*GAMBDA(1)))/ 137 : & (2.D0*(UN-XNU)*((2.D0*PHI)+(XGGG*GAMBDA(1)))) 138 : HHH1(2,1)=HHH1(1,2) 139 : HHH1(3,3) = XGGG*PHI/(PHI+XGGG*GAMBDA(1)) 139 : CALL NORMBL(SIGF,FG(1),ISGNTC,1,2,XNOR1) 140 : GOTO 1010 141 : ELSE 141 : CALL COURBU(SIGF,FG(1),ISGNTC,1,2,XCOU1) 141 : CALL NORMBL(SIGF,FG(1),ISGNTC,1,2,XNOR1) 142 : ENDIF 143 : ENDIF 144 : C 145 : C 146 : IF (ICAS1.EQ.ICAS2) THEN 147 : GOTO 1000 148 : ELSE 149 : IF (ICAS2.LE.4) THEN 150 : FG(2) = RESIST(ICAS2) 151 : ELSE 152 : FG(2) = XKISO(ICAS2-4)*RESIST(ICAS2) 153 : ENDIF 154 : IF (ICAS2.EQ.1) THEN 155 : ISGNTCB = 1 156 : ELSE IF (ICAS2.EQ.2) THEN 157 : ISGNTCB = 1 158 : ELSE IF (ICAS2.EQ.3) THEN 159 : ISGNTCB = -1 160 : ELSE IF (ICAS2.EQ.4) THEN 161 : ISGNTCB = -1 162 : ELSE IF (ICAS2.EQ.5) THEN 163 : ISGNTCB = 1 164 : ELSE IF (ICAS2.EQ.6) THEN 165 : ISGNTCB = -1 166 : ENDIF 167 : IF (ICAS2.LE.4) THEN 168 : IF (SIGP(2).GE.SIGP(1)) THEN 169 : IF ((ICAS2.EQ.1).OR.(ICAS2.EQ.3)) THEN 170 : ISGN2B = -1 171 : ELSE 172 : ISGN2B = 1 173 : ENDIF 174 : ELSE 175 : IF ((ICAS2.EQ.1).OR.(ICAS2.EQ.3)) THEN 176 : ISGN2B = 1 177 : ELSE 178 : ISGN2B = -1 179 : ENDIF 180 : ENDIF 181 : ELSE 182 : ISGN2B = 1 183 : ENDIF 184 : C 185 : IF (ICAS2.LE.4) THEN 185 : CALL COURBU(SIGF,FG(2),ISGNTCB,ISGN2B,1,XCOU2) 185 : CALL NORMBL(SIGF,FG(2),ISGNTCB,ISGN2B,1,XNOR2) 186 : ELSE 186 : CALL COURBU(SIGF,FG(2),ISGNTCB,1,2,XCOU2) 186 : CALL NORMBL(SIGF,FG(2),ISGNTCB,1,2,XNOR2) 187 : ENDIF 188 : C 189 : GOTO 2000 190 : ENDIF 191 : C 192 : C 193 : C 194 : 1000 CONTINUE 195 : C 196 : C Cas 1 mecanisme 197 : C 198 : C Calcul de l'inverse de H et de H 198 : CALL PLMATR(HHH1,UN,XCOU1,GAMBDA(1),3,3,HHH1) 198 : CALL INVALM(HHH1,3,3,IRD,1.D-8) 199 : IF (IRD.NE.0) THEN 200 : WRITE (IOIMP,*)'BETOKT: Determinant nul' 201 : GOTO 9999 202 : ENDIF 203 : C Point d'entree dans le cas de l'apex 204 : 1010 CONTINUE 205 : C Calcul de U,V et E-1 205 : CALL ZERO(UUU,3,2) 205 : CALL ZERO(VVV,2,3) 205 : CALL ZERO(UUU2,3,2) 205 : CALL ZERO(VVV2,2,3) 206 : UUU(1,1) = XNOR1(1) 207 : UUU(2,1) = XNOR1(2) 208 : UUU(3,1) = XNOR1(4) 209 : VVV(1,1) = XNOR1(1) 210 : VVV(1,2) = XNOR1(2) 211 : VVV(1,3) = XNOR1(4) 212 : C 213 : IF (ICAS1.LE.4) THEN 214 : EEE(1,1) = XHH1 215 : ELSE IF (ICAS1.EQ.5) THEN 216 : EEE(1,1) = ISGNTC*XPEN1*RESIST(ICAS1) 217 : ELSE IF (ICAS1.EQ.6) THEN 218 : EEE(1,1) = ISGNTC*XPEN2*RESIST(ICAS1) 219 : ENDIF 220 : C 220 : CALL MATMAT(HHH1,UUU,3,3,2,UUU2) 220 : CALL MATMAT(VVV,HHH1,2,3,3,VVV2) 221 : C 221 : CALL MATMAT(VVV,UUU2,2,3,2,HHH2) 222 : C 223 : EEE(1,1) = EEE(1,1) + HHH2(1,1) 224 : IF (ABS(EEE(1,1)).GT.EPSILO) THEN 225 : EEE(1,1) = -UN/EEE(1,1) 226 : ELSE 227 : WRITE (IOIMP,*)'BETOKT: Determinant nul' 228 : GOTO 9999 229 : ENDIF 229 : CALL MATMAT(UUU2,VVV2,3,2,3,HHH3) 229 : CALL PLMATR(HHH1,UN,HHH3,EEE(1,1),3,3,HHH1) 230 : C 231 : GOTO 3000 232 : C 233 : 2000 CONTINUE 234 : C 235 : C Cas 2 mecanismes 236 : C 237 : C Calcul de H-1 et H 237 : CALL PLMATR(HHH1,UN,XCOU1,GAMBDA(1),3,3,HHH1) 237 : CALL PLMATR(HHH1,UN,XCOU2,GAMBDA(2),3,3,HHH1) 237 : CALL INVALM(HHH1,3,3,IRD,1.D-8) 238 : IF (IRD.NE.0) THEN 239 : WRITE (IOIMP,*)'BETOKT: Determinant nul' 240 : GOTO 9999 241 : ENDIF 242 : C Calcul de U ,V et E 243 : UUU(1,1) = XNOR1(1) 244 : UUU(2,1) = XNOR1(2) 245 : UUU(3,1) = XNOR1(4) 246 : UUU(1,2) = XNOR2(1) 247 : UUU(2,2) = XNOR2(2) 248 : UUU(3,2) = XNOR2(4) 249 : VVV(1,1) = XNOR1(1) 250 : VVV(1,2) = XNOR1(2) 251 : VVV(1,3) = XNOR1(4) 252 : VVV(2,1) = XNOR2(1) 253 : VVV(2,2) = XNOR2(2) 254 : VVV(2,3) = XNOR2(4) 255 : C 256 : IF (ICAS1.LE.4) THEN 257 : EEE(1,1) = XHH1 258 : IF (ICAS2.EQ.6) THEN 259 : EEE(2,1) = DEGRAD*ISGNTCB*RESIST(6) 260 : &*ISGNTC*RESIST(ICAS1)*XPEN2 261 : ELSE 262 : EEE(2,1) = XZER 263 : ENDIF 264 : ELSE IF (ICAS1.EQ.5) THEN 265 : EEE(1,1) = ISGNTC*XPEN1*RESIST(ICAS1) 266 : EEE(2,1) = XZER 267 : ELSE IF (ICAS1.EQ.6) THEN 268 : EEE(1,1) = ISGNTC*XPEN2*RESIST(ICAS1) 269 : EEE(2,1) = XZER 270 : ENDIF 271 : IF (ICAS2.LE.4) THEN 272 : EEE(2,2) = XHH1 273 : IF (ICAS1.EQ.6) THEN 274 : EEE(1,2) = DEGRAD*ISGNTC*RESIST(6) 275 : &*ISGNTCB*RESIST(ICAS2)*XPEN2 276 : ELSE 277 : EEE(1,2) = XZER 278 : ENDIF 279 : ELSE IF (ICAS2.EQ.5) THEN 280 : EEE(2,2) = ISGNTCB*XPEN1*RESIST(ICAS2) 281 : EEE(1,2) = XZER 282 : ELSE IF (ICAS2.EQ.6) THEN 283 : EEE(2,2) = ISGNTCB*XPEN2*RESIST(ICAS2) 284 : EEE(1,2) = XZER 285 : ENDIF 286 : C 286 : CALL MATMAT(HHH1,UUU,3,3,2,UUU2) 286 : CALL MATMAT(VVV,HHH1,2,3,3,VVV2) 287 : C 287 : CALL MATMAT(VVV,UUU2,2,3,2,HHH2) 288 : C 288 : CALL PLMATR(EEE,UN,HHH2,UN,2,2,EEE) 289 : XDET = EEE(1,1)*EEE(2,2) - EEE(1,2)*EEE(2,1) 290 : IF (ABS(XDET).GT.EPSILO) THEN 291 : XDET = UN/XDET 292 : EEE2(1,1) = XDET*EEE(2,2) 293 : EEE2(2,2) = XDET*EEE(1,1) 294 : EEE2(2,1) = -XDET*EEE(1,2) 295 : EEE2(1,2) = -XDET*EEE(2,1) 296 : ELSE 297 : WRITE (IOIMP,*)'BETOKT: Determinant nul' 298 : GOTO 9999 299 : ENDIF 299 : CALL MATMAT(UUU2,EEE2,3,2,2,UUU3) 299 : CALL MATMAT(UUU3,VVV2,3,2,3,HHH3) 299 : CALL PLMATR(HHH1,UN,HHH3,-UN,3,3,HHH1) 300 : C 301 : GOTO 3000 302 : C 303 : C============================================================== 304 : 3000 CONTINUE 304 : CALL ZERO(XKTAT,4,4) 305 : XKTAT(1,1) = HHH1(1,1) 306 : XKTAT(1,2) = HHH1(1,2) 307 : XKTAT(2,1) = HHH1(2,1) 308 : XKTAT(2,2) = HHH1(2,2) 309 : XKTAT(4,4) = HHH1(3,3) 310 : XKTAT(4,1) = HHH1(3,1) 311 : XKTAT(4,2) = HHH1(3,2) 312 : XKTAT(1,4) = HHH1(1,3) 313 : XKTAT(2,4) = HHH1(2,3) 314 : 9999 CONTINUE 315 : C=============================================================== 316 : C Fin de la routine 317 : C=============================================================== 318 : RETURN 319 : END 320 : 321 : 322 :
© Cast3M 2003 - All rights reserved.
Disclaimer