1 : C INCRE7 SOURCE CHAT 05/01/13 00:35:30 5004 2 : SUBROUTINE INCRE7(SIG,VAR,DSPT,EPSVPT,VARPT,XMAT,NSTRS0, 3 : & MFR,NVARI,NCOMAT,IFOUR) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8(A-H,O-Z) 6 : DIMENSION SIG(*),VAR(*),DSPT(*) 7 : DIMENSION EPSVPT(*),VARPT(*),XMAT(*) 8 : DIMENSION XX(6),YY(6),EVPT(6),EPPT(6),SIGPT(6) 9 : DIMENSION SIG0(6),EPS0(6) 10 : PARAMETER (AMAX = 1.0D20 , AMIN = 1.D-10) 11 : DETIER = 2.0D0/3.0D0 12 : ROOT = SQRT(3.0D0/2.0D0) 12 : CALL ZDANUL(VARPT,NVARI) 13 : IF (MFR.EQ.5) THEN 14 : NSTRS=6 15 : SIG0(1)=SIG(1) 16 : SIG0(2)=SIG(2) 17 : SIG0(3)=0.D0 18 : SIG0(4)=SIG(3) 19 : SIG0(5)=SIG(4) 20 : SIG0(6)=SIG(5) 21 : ELSE 22 : NSTRS=NSTRS0 23 : DO 10 I=1,NSTRS 24 : SIG0(I)=SIG(I) 25 : 10 CONTINUE 26 : ENDIF 27 : TRACE =(SIG0(1)+SIG0(2)+SIG0(3))/3.0D0 28 : C-------------------------------------------------------------------| 29 : C COEFFICIENTS MATERIAU | 30 : C-------------------------------------------------------------------| 31 : YOU =XMAT( 1) 32 : XNU =XMAT( 2) 33 : RP0 =XMAT( 5) 34 : QP =XMAT( 6) 35 : BP =XMAT( 7) 36 : CP1 =XMAT( 8) 37 : CP2 =XMAT( 9) 38 : DP1 =XMAT(10) 39 : DP2 =XMAT(11) 40 : XK =XMAT(12) 41 : XN =XMAT(13) 42 : RV0 =XMAT(14) 43 : QV =XMAT(15) 44 : BV =XMAT(16) 45 : CV1 =XMAT(17) 46 : CV2 =XMAT(18) 47 : DV1 =XMAT(19) 48 : DV2 =XMAT(20) 49 : CVP1=XMAT(21) 50 : CVP2=XMAT(22) 51 : 52 : IF (CP1.EQ.0.D0) THEN 53 : DCP1 = 0.D0 54 : ELSE 55 : DCP1 = DP1 / CP1 56 : ENDIF 57 : 58 : IF (CP2.EQ.0.D0) THEN 59 : DCP2 = 0.D0 60 : ELSE 61 : DCP2 = DP2 / CP2 62 : ENDIF 63 : 64 : IF (CV1.EQ.0.D0) THEN 65 : DCV1 = 0.D0 66 : ELSE 67 : DCV1 = DV1 / CV1 68 : ENDIF 69 : 70 : IF (CV2.EQ.0.D0) THEN 71 : DCV2 = 0.D0 72 : ELSE 73 : DCV2 = DV2 / CV2 74 : ENDIF 75 : 76 : C-------------------------------------------------------------------| 77 : C PARTIE VISCOPLASTIQUE | 78 : C-------------------------------------------------------------------I 79 : C 80 : C *** CALCUL DE XX=DEVIATEUR(SIGMA-XP) 81 : C 82 : TRACE =(SIG0(1)+SIG0(2)+SIG0(3))/3.0D0 83 : DO 71 I=1,NSTRS 84 : A = 0.0D0 85 : IF (I.LE.3) A=1.0D0 86 : XX(I) = SIG0(I)-A*TRACE - VAR(5*NSTRS+I)-VAR(6*NSTRS+I) 87 : 71 CONTINUE 88 : C 89 : C *** CRITERE VISCOPLASTIQUE 90 : C 91 : AJ2 = PROCON(XX,XX,NSTRS) 92 : AJ2 = SQRT(1.5D0*AJ2) 93 : CRIT= AJ2-VAR(10*NSTRS+4)-RV0 94 : C 95 : C *** ECOULEMENT 96 : C 97 : IF(CRIT.GT.0.D0) THEN 98 : VPT=(CRIT/XK)**XN 99 : ELSE 100 : VPT=0.D0 101 : ENDIF 102 : VARPT(10*NSTRS+3)=VPT 103 : DO 80 I=1,NSTRS 104 : IF(VPT.EQ.0.D0) THEN 105 : EVPT(I)=0.D0 106 : VARPT(7*NSTRS+I)=0.D0 107 : VARPT(8*NSTRS+I)=0.D0 108 : VARPT(9*NSTRS+I)=0.D0 109 : ELSE 110 : EVPT(I)=1.5*XX(I)*VPT/AJ2 111 : VARPT(9*NSTRS+I)=EVPT(I) 112 : VARPT(7*NSTRS+I)=EVPT(I)-(1.5*DCV1*VAR(5*NSTRS+I)*VPT) 113 : VARPT(8*NSTRS+I)=EVPT(I)-(1.5*DCV2*VAR(6*NSTRS+I)*VPT) 114 : ENDIF 115 : 80 CONTINUE 116 : C-------------------------------------------------------------------| 117 : C PARTIE PLASTIQUE | 118 : C-------------------------------------------------------------------I 119 : C 120 : C *** CALCUL DE XX=DEVIATEUR(SIGMA-XP)/J2(SIGMA-X) 121 : C 122 : DO 70 I=1,NSTRS 123 : A = 0.0D0 124 : IF (I.LE.3) A=1.0D0 125 : XX(I) = SIG0(I)-A*TRACE - VAR(I)-VAR(NSTRS+I) 126 : 70 CONTINUE 127 : AJ2 = PROCON (XX,XX,NSTRS) 128 : AJ2 = SQRT(1.5D0*AJ2) 129 : IF(AJ2.NE.0.D0) THEN 130 : DO 81 I=1,NSTRS 131 : XX(I)=1.5*XX(I)/AJ2 132 : 81 CONTINUE 133 : ENDIF 134 : C 135 : C *** CRITERE PLASTIQUE 136 : C 137 : CRIT=AJ2-VAR(10*NSTRS+2)-RP0 138 : C 139 : C *** CALCUL DE PPT1=XX:(E:(EPSPT-EVPT)-2/3*SIGVPPT) 140 : C 141 : IF(CRIT.LE.0.D0) THEN 142 : PPT=0.D0 143 : VARPT(10*NSTRS+1)=PPT 144 : DO 85 I=1,NSTRS 145 : VARPT(2*NSTRS+I)=0.D0 146 : VARPT(3*NSTRS+I)=0.D0 147 : VARPT(4*NSTRS+I)=0.D0 148 : EPPT(I)=0.D0 149 : 85 CONTINUE 150 : ELSE 151 : IF(IFOUR.EQ.-2) THEN 152 : X2MU=YOU/(1.0+XNU) 153 : XCO = X2MU/(1.0-XNU) 154 : SIGPT(1)=DSPT(1)-XCO*(EVPT(1)+XNU*EVPT(2)) 155 : SIGPT(2)=DSPT(2)-XCO*(EVPT(2)+XNU*EVPT(1)) 156 : SIGPT(4)=DSPT(4)-X2MU*EVPT(4) 157 : SIGPT(3)=0.0D0 158 : ELSEIF (IFOUR.EQ.2.OR.IFOUR.EQ.0.OR.IFOUR.EQ.1.OR. 159 : & IFOUR.EQ.-1.OR.IFOUR.EQ.-3) THEN 160 : X2MU = YOU/(1+XNU) 161 : ALMB = X2MU*XNU/(1-2*XNU) 162 : IF (MFR.EQ.5) THEN 163 : X2MU=YOU/(1.0+XNU) 164 : XCO = X2MU/(1.0-XNU) 165 : SIGPT(1)=DSPT(1)-XCO*(EVPT(1)+XNU*EVPT(2)) 166 : SIGPT(2)=DSPT(2)-XCO*(EVPT(2)+XNU*EVPT(1)) 167 : SIGPT(3)=DSPT(3)-X2MU*EVPT(3) 168 : SIGPT(4)=DSPT(4)-X2MU*EVPT(4) 169 : SIGPT(5)=DSPT(5)-X2MU*EVPT(5) 170 : ELSE 171 : TRACE=EVPT(1)+EVPT(2)+EVPT(3) 172 : DO 20 I=1,NSTRS,1 173 : A = 0.0 174 : IF (I.LE.3) A=1.0 175 : SIGPT(I) = DSPT(I)-ALMB*A*TRACE-X2MU*EVPT(I) 176 : 20 CONTINUE 177 : ENDIF 178 : C= Modes de calcul 1D contraintes planes suivant z (DYCZ et GYCZ) 179 : ELSE IF (IFOUR.EQ.4.OR.IFOUR.EQ.8) THEN 180 : XCO = YOU/(1.-XNU*XNU) 181 : SIGPT(1) = DSPT(1)-XCO*(EVPT(1)+XNU*EVPT(2)) 182 : SIGPT(2) = DSPT(2)-XCO*(EVPT(2)+XNU*EVPT(1)) 183 : SIGPT(3) = 0.D0 184 : C= Modes de calcul 1D contraintes planes suivant y (CYDZ et CYGZ) et 185 : C= 1D axisymetrique contraintes axiales nulles (AXCZ) 186 : ELSE IF (IFOUR.EQ.5.OR.IFOUR.EQ.10.OR.IFOUR.EQ.13) THEN 187 : XCO = YOU/(1.-XNU*XNU) 188 : SIGPT(1) = DSPT(1)-XCO*(EVPT(1)+XNU*EVPT(3)) 189 : SIGPT(2) = 0.D0 190 : SIGPT(3) = DSPT(3)-XCO*(EVPT(3)+XNU*EVPT(1)) 191 : C= Modes de calcul 1D contraintes planes suivant y et z (CYCZ) 192 : ELSE IF (IFOUR.EQ.6) THEN 193 : SIGPT(1) = DSPT(1)-YOU*EVPT(1) 194 : SIGPT(2) = 0.D0 195 : SIGPT(3) = 0.D0 196 : C= Autres modes de calcul 1D deformations planes (DYDZ GYDZ DYGZ GYGZ) 197 : C= 1D axisymetrique (AXDZ AXGZ) et 1D spherique 198 : ELSE IF (IFOUR.EQ. 3.OR.IFOUR.EQ. 7.OR.IFOUR.EQ. 9.OR.IFOUR.EQ.11 199 : . .OR.IFOUR.EQ.12.OR.IFOUR.EQ.14.OR.IFOUR.EQ.15) THEN 200 : X2MU = YOU/(1.+XNU) 201 : TRACE = EVPT(1)+EVPT(2)+EVPT(3) 202 : XCO = X2MU*XNU*TRACE/(1.-XNU) 203 : SIGPT(1) = DSPT(1)-X2MU*EVPT(1)-XCO 204 : SIGPT(2) = DSPT(2)-X2MU*EVPT(2)-XCO 205 : SIGPT(3) = DSPT(3)-X2MU*EVPT(3)-XCO 206 : ENDIF 207 : DO 30 I=1,NSTRS 208 : SIGPT(I)=SIGPT(I)-DETIER*( 209 : & CVP1*VARPT(7*NSTRS+I)-CVP2*VARPT(8*NSTRS+I)) 210 : 30 CONTINUE 211 : PPT1=PROCON(SIGPT,XX,NSTRS) 212 : C 213 : C *** CALCUL DE 214 : C H1=2/3*XX:(C1P(NP-3/2*(D1P/C1P)*X1P)+C2P(NP-3/2*(D2P/C2P)*X2P)) 215 : C 216 : DO 40 I=1,NSTRS 217 : YY(I)=CP1*(XX(I)-1.5*DCP1*VAR(I))+ 218 : & CP2*(XX(I)-1.5*DCP2*VAR(NSTRS+I)) 219 : 40 CONTINUE 220 : H1=DETIER*PROCON(XX,YY,NSTRS) 221 : C 222 : C *** CALCUL DE H2=XX:(E:XX) 223 : C 224 : IF(IFOUR.EQ.-2) THEN 225 : X2MU=YOU/(1.0+XNU) 226 : XCO = X2MU/(1.0-XNU) 227 : YY(1)=XCO*(XX(1)+XNU*XX(2)) 228 : YY(2)=XCO*(XX(2)+XNU*XX(1)) 229 : YY(4)=X2MU*XX(4) 230 : YY(3)=0.0D0 231 : ELSEIF (IFOUR.EQ.2.OR.IFOUR.EQ.0.OR.IFOUR.EQ.1.OR. 232 : & IFOUR.EQ.-1.OR.IFOUR.EQ.-3) THEN 233 : X2MU = YOU/(1+XNU) 234 : ALMB = X2MU*XNU/(1-2*XNU) 235 : IF (MFR.EQ.5) THEN 236 : X2MU=YOU/(1.0+XNU) 237 : XCO = X2MU/(1.0-XNU) 238 : YY(1)=XCO*(XX(1)+XNU*XX(2)) 239 : YY(2)=XCO*(XX(2)+XNU*XX(1)) 240 : YY(3)=X2MU*XX(3) 241 : YY(4)=X2MU*XX(4) 242 : YY(5)=X2MU*XX(5) 243 : ELSE 244 : TRACE=XX(1)+XX(2)+XX(3) 245 : DO 50 I=1,NSTRS,1 246 : A = 0.0 247 : IF (I.LE.3) A=1.0 248 : YY(I) = ALMB*A*TRACE+X2MU*XX(I) 249 : 50 CONTINUE 250 : ENDIF 251 : C= Modes de calcul 1D contraintes planes suivant z (DYCZ et GYCZ) 252 : ELSE IF (IFOUR.EQ.4.OR.IFOUR.EQ.8) THEN 253 : XCO = YOU/(1.-XNU*XNU) 254 : YY(1) = XCO*(XX(1)+XNU*XX(2)) 255 : YY(2) = XCO*(XX(2)+XNU*XX(1)) 256 : YY(3) = 0.D0 257 : C= Modes de calcul 1D contraintes planes suivant y (CYDZ et CYGZ) et 258 : C= 1D axisymetrique contraintes axiales nulles (AXCZ) 259 : ELSE IF (IFOUR.EQ.5.OR.IFOUR.EQ.10.OR.IFOUR.EQ.13) THEN 260 : XCO = YOU/(1.-XNU*XNU) 261 : YY(1) = XCO*(XX(1)+XNU*XX(3)) 262 : YY(2) = 0.D0 263 : YY(3) = XCO*(XX(3)+XNU*XX(1)) 264 : C= Modes de calcul 1D contraintes planes suivant y et z (CYCZ) 265 : ELSE IF (IFOUR.EQ.6) THEN 266 : YY(1) = YOU*XX(1) 267 : YY(2) = 0.D0 268 : YY(3)=0.D0 269 : C= Autres modes de calcul 1D deformations planes (DYDZ GYDZ DYGZ GYGZ) 270 : C= 1D axisymetrique (AXDZ AXGZ) et 1D spherique 271 : ELSE IF (IFOUR.EQ. 3.OR.IFOUR.EQ. 7.OR.IFOUR.EQ. 9.OR.IFOUR.EQ.11 272 : . .OR.IFOUR.EQ.12.OR.IFOUR.EQ.14.OR.IFOUR.EQ.15) THEN 273 : X2MU = YOU/(1.+XNU) 274 : TRACE = XX(1)+XX(2)+XX(3) 275 : XCO = X2MU*XNU*TRACE/(1.-XNU) 276 : YY(1) = X2MU*XX(1)+XCO 277 : YY(2) = X2MU*XX(2)+XCO 278 : YY(3) = X2MU*XX(3)+XCO 279 : ENDIF 280 : H2=PROCON(XX,YY,NSTRS) 281 : C 282 : C *** CALCUL DE H3=BP*(RP-RP0-QP) 283 : C 284 : H3=BP*(VAR(10*NSTRS+2)-QP) 285 : C 286 : C *** MODULE PLASTIQUE 287 : C 288 : H=H1+H2-H3 289 : C 290 : C *** CALCUL DE PPT,EPPT,ALPHAPT 291 : C 292 : PPT=PPT1/H 293 : IF(PPT.LT.0.D0) PPT=0.D0 294 : VARPT(10*NSTRS+1)=PPT 295 : DO 60 I=1,NSTRS 296 : EPPT(I)=XX(I)*PPT 297 : VARPT(2*NSTRS+I)=EPPT(I)-(1.5*DCP1*VAR(I)*PPT) 298 : VARPT(3*NSTRS+I)=EPPT(I)-(1.5*DCP2*VAR(NSTRS+I)*PPT) 299 : VARPT(4*NSTRS+I)=EPPT(I) 300 : 60 CONTINUE 301 : ENDIF 302 : C 303 : IF (MFR.EQ.5) THEN 304 : EPSVPT(1)=EVPT(1)+EPPT(1) 305 : EPSVPT(2)=EVPT(2)+EPPT(2) 306 : EPSVPT(3)=EVPT(4)+EPPT(4) 307 : EPSVPT(4)=EVPT(5)+EPPT(5) 308 : EPSVPT(5)=EVPT(6)+EPPT(6) 309 : ELSE 310 : DO 11 I=1,NSTRS 311 : EPSVPT(I)=EVPT(I)+EPPT(I) 312 : 11 CONTINUE 313 : ENDIF 314 : C 315 : VARPT(10*NSTRS+5)=VPT+PPT 316 : RETURN 317 : END 318 : 319 : 320 : 321 : 322 : 323 : 324 : 325 : 326 : 327 : 328 : 329 :
© Cast3M 2003 - All rights reserved.
Disclaimer