1 : C RICJ3 SOURCE CHAT 11/03/16 21:30:40 6902 2 : C sub ricjoi3d 3 : C RICJOI3D SOURCE CHAT 05/01/12 22:24:09 5004 4 : SUBROUTINE RICJ3(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI, 5 : & DEPST,IFOUR,TETA1,TETA2,XMAT, 6 : & NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE) 7 : C----------------------------------------------------------------------- 8 : C Modele RICJOI en version element joint 3D 9 : C----------------------------------------------------------------------- 10 : C Développé par : Benjamin Richard 11 : C 12 : C Supervisé par : Frederic Ragueneau 13 : C Lucas Hector Adelaide 14 : C Christian Cremona 15 : C Jean Louis Tailhan 16 : C----------------------------------------------------------------------- 17 : 18 : C 19 : C-----DECLARATION GENERALE---------------------------------------------- 20 : C 21 : IMPLICIT INTEGER(I-N) 22 : IMPLICIT REAL*8(A-H,O-Z) 23 : C 24 : C-----APPEL AUX INCLUDES------------------------------------------------ 25 : C 26 : -INC CCREEL 27 : C 28 : C-----PARAMETRES 29 : C 30 : PARAMETER(XDILAT = 80.0D0) 31 : C 32 : C-----DIMENSION--------------------------------------------------------- 33 : C 34 : DIMENSION SIG0(*),EPIN0(*),VAR0(*),DEPST(*),XMAT(*) 35 : DIMENSION DD(NSTRS,NSTRS) 36 : DIMENSION SIGF(*),DEFP(*),VARF(*) 37 : 38 : FVARLIM = 2.0D-5 39 : C 40 : C-----COMPTEUR POUR LE NOMBRE DE PASSAGE DAS LA BOUCLE 41 : C 42 : COUN = VAR0(17) 43 : C 44 : C-----ETAT PRECEDENT---------------------------------------------------- 45 : C 46 : D = VAR0(1 ) 47 : EPSF1 = VAR0(2 ) 48 : EPSF2 = VAR0(3 ) 49 : XEC1 = VAR0(4 ) 50 : XEC2 = VAR0(5 ) 51 : EPST1 = VAR0(6 ) 52 : EPST2 = VAR0(7 ) 53 : EPSN = VAR0(8 ) 54 : CIN1 = VAR0(9 ) 55 : CIN2 = VAR0(10) 56 : EPSR = VAR0(11) 57 : IFLAG = nint(VAR0(12)) 58 : FVAR = VAR0(13) 59 : SIGM = VAR0(14) 60 : XKS = VAR0(15) 61 : XKN = VAR0(16) 62 : crit2tol=1.001d0 63 : C 64 : C-----INITIALISATION SPECIFIQUE POUR LES VARIABLES INTERNES DIFF. DE ZERO 65 : C 66 : IF (COUN.EQ.0.0D0) THEN 67 : FVAR = 0.2D0 68 : SIGM = XMAT(13) 69 : XKS = XMAT(1) 70 : XKN = XMAT(2) 71 : ENDIF 72 : C 73 : C-----CONSTANTES MATERIAUX---------------------------------------------- 74 : C 75 : XKS0 = XMAT(1) 76 : XKN0 = XMAT(2) 77 : 78 : XAD = XMAT(5) 79 : XY0 = XMAT(6) 80 : XAL = XMAT(7) 81 : XGA = XMAT(8) 82 : XAA = XMAT(9) 83 : 84 : Q1COE = XMAT(10) 85 : Q2COE = XMAT(11) 86 : Q3COE = XMAT(12) 87 : SYCOE = XMAT(13) 88 : NCOEF = nint(XMAT(14) ) 89 : KCOEF = nint(XMAT(15)) 90 : 91 : TCCOR = XMAT(16) 92 : NGONF = nint(XMAT(17)) 93 : 93 : CALL APP3(TCCOR,YR,15) 94 : 95 : 96 : C 97 : C-----ACTUALISATION DE LA DEFORMATION----------------------------------- 98 : C 99 : EPST1 = EPST1 + DEPST(1) 100 : EPST2 = EPST2 + DEPST(2) 101 : EPSN = EPSN + DEPST(3) 102 : 103 : IF (NGONF.EQ.1) THEN 104 : GOTO 123 105 : ENDIF 106 : 107 : XSTE = XKN*EPSN 108 : XSAD = XDILAT*0.5D0*(-DABS(XSTE)+XSTE) 109 : 110 : C 111 : C-----CALCUL DE L ENERGIE----------------------------------------------- 112 : C 113 : EPSNP = 0.5D0*(DABS(EPSN)+EPSN) 114 : YN = 0.5D0*XKN*EPSNP*EPSNP 115 : YD1 = 0.5D0*XKS*EPST1*EPST1 116 : YD2 = 0.5D0*XKS*EPST2*EPST2 117 : 118 : C 119 : C-----CALCUL DU SEUIL D ENDOMMAGEMENT 1--------------------------------- 120 : C 121 : XFD1 = XAL*YN+YD1-(XY0+XEC1+YR) 122 : 123 : C 124 : C-----TEST SUR LE SEUIL D ENDOMMAGEMENT--------------------------------- 125 : C 126 : IF (XFD1.GT.0.0D0) THEN 127 : Dn = 1.0D0-1.0D0/(1+XAD*(XAL*YN+YD1-XY0-YR)) 128 : XEC1 = XAL*YN+YD1-XY0-YR 129 : IF (Dn.GT.D) THEN 130 : D = Dn 131 : ENDIF 132 : ENDIF 133 : IF (D.GE.1.0D0) THEN 134 : D = 0.999D0 135 : ENDIF 136 : 137 : C 138 : C-----CALCUL DU SEUIL D ENDOMMAGEMENT 2--------------------------------- 139 : C 140 : XFD2 = XAL*YN+YD2-(XY0+XEC2+YR) 141 : 142 : C 143 : C-----TEST SUR LE SEUIL D ENDOMMAGEMENT--------------------------------- 144 : C 145 : IF (XFD2.GT.0.0D0) THEN 146 : Dn = 1.0D0-1.0D0/(1+XAD*(XAL*YN+YD2-XY0-YR)) 147 : XEC2 = XAL*YN+YD2-XY0-YR 148 : IF (Dn.GT.D) THEN 149 : D = Dn 150 : ENDIF 151 : ENDIF 152 : IF (D.GE.1.0D0) THEN 153 : D = 0.999D0 154 : ENDIF 155 : 156 : C 157 : C-----FROTTEMENT 1------------------------------------------------------- 158 : C 159 : SPI1 = D*XKS*(EPST1-EPSF1) 160 : 161 : CRITF1= DABS(SPI1-CIN1)+XSAD 162 : 163 : IF (CRITF1.GT.(0.0D0)) THEN 164 : XF0 = CRITF1 165 : XDLAM = 1.0D0 166 : 167 : DO I=1,500 168 : 169 : SEUIL = (DABS(SPI1-CIN1)+XSAD)/XF0 170 : 171 : IF ((SEUIL.LE.1.0D-5).OR.(XDLAM.LE.1.0D-10)) THEN 172 : GOTO 900 173 : ELSE 174 : XSG = DABS(SPI1-CIN1)/(SPI1-CIN1) 175 : XDLAM = (DABS(SPI1-CIN1)+XSAD)/ 176 : & (XKS*D-XSG*XGA*(-XSG+XAA*CIN1)) 177 : CIN1 = CIN1-XGA*XDLAM*(-XSG+XAA*CIN1) 178 : SPI1 = SPI1-D*XKS*XDLAM*XSG 179 : ENDIF 180 : 181 : EPSF1 = -SPI1/(D*XKS)+EPST1 182 : 183 : ENDDO 184 : ENDIF 185 : 186 : 900 CONTINUE 187 : 188 : C 189 : C-----FROTTEMENT 2------------------------------------------------------- 190 : C 191 : SPI2 = D*XKS*(EPST2-EPSF2) 192 : 193 : CRITF2= DABS(SPI2-CIN2)+XSAD 194 : 195 : IF (CRITF2.GT.(0.0D0)) THEN 196 : XF0 = CRITF2 197 : XDLAM = 1.0D0 198 : 199 : DO I=1,500 200 : 201 : SEUIL = (DABS(SPI2-CIN2)+XSAD)/XF0 202 : 203 : IF ((SEUIL.LE.1.0D-5).OR.(XDLAM.LE.1.0D-10)) THEN 204 : GOTO 901 205 : ELSE 206 : XSG = DABS(SPI2-CIN2)/(SPI2-CIN2) 207 : XDLAM = (DABS(SPI2-CIN2)+XSAD)/ 208 : & (XKS*D-XSG*XGA*(-XSG+XAA*CIN2)) 209 : CIN2 = CIN2-XGA*XDLAM*(-XSG+XAA*CIN2) 210 : SPI2 = SPI2-D*XKS*XDLAM*XSG 211 : ENDIF 212 : 213 : EPSF2 = -SPI2/(D*XKS)+EPST2 214 : 215 : ENDDO 216 : ENDIF 217 : 218 : 901 CONTINUE 219 : 220 : C 221 : C-----ANELASTICITE-------------------------------------------------------- 222 : C 223 : IF (EPSN.GT.0.0D0) THEN 224 : SIGN = (1-D)*XKN*EPSN 225 : ELSE 226 : SIGN = XKN*(EPSN-EPSR) 227 : ENDIF 228 : 229 : C 230 : C-----CONTRAINTE DUE AUX EFFETS INELASTIQUE DE LA ROUILLE----------------- 231 : C 232 : sigr=sign 233 : IF (EPSN.LE.0.0D0) THEN 234 : SIGR = -1.0D0*SIGN 235 : ENDIF 236 : 237 : C 238 : C-----CRITERE D ANELASTICITE DE LA ROUILLE-------------------------------- 239 : C 240 : IF (IFLAG.EQ.1) THEN 241 : GOTO 20 242 : ELSE 243 : CRIT2 = 2.0D0*Q1COE*FVAR*DCOSH(Q2COE*SIGR/2.0D0/SYCOE)- 244 : & (1.0D0+(Q3COE*FVAR)**2.0D0) 245 : ENDIF 246 : 247 : * print*,'crit2=',crit2 248 : 249 : C 250 : C SI LE CRITERE EST DEPASSE ON DEMARRE LE RETURN MAPPING 251 : C 252 : IF ((CRIT2.GT.0.0D0).AND. 253 : 1 (EPSN .LE.0.0D0).AND. 254 : 2 (TCCOR.GT.0.0D0).AND. 255 : 3 (FVAR.GT.FVARLIM).AND. 256 : 4 (IFLAG.EQ.0).AND. 257 : 5 ((EPSN*DEPST(3)).GT.0.0D0)) THEN 258 : CRIT20 = CRIT2 259 : DLAM2 = 1.0D0 260 : GOTO 10 261 : ELSE IF (FVAR.EQ.FVARLIM) THEN 262 : IFLAG = 1 263 : GOTO 20 264 : ELSE 265 : GOTO 30 266 : END IF 267 : C 268 : C-----DEBUT DES ITERATIONS INTERNES 269 : C 270 : 10 DO I=1,10000 271 : C 272 : C EVALUATION DU CRITERE 273 : C 274 : CRIT2 = 2.0D0*Q1COE*FVAR*DCOSH(Q2COE*SIGR/2.0D0/SYCOE)- 275 : & (1.0D0+(Q3COE*FVAR)**2.0D0) 276 : C 277 : C TEST DE CONVERGENCE SUR LA VALEUR DU CRITERE RELATIVE 278 : C 279 : IF ((DABS(CRIT2/CRIT20).LE.CRIT2TOL).OR. 280 : & (DLAM2.LT.0.0D0).OR. 281 : & (CRIT2.LT.0.0D0)) THEN 282 : GOTO 20 283 : END IF 284 : C 285 : C CALCUL DES DERIVEES DU CRITERES 286 : C 287 : C dFr/dSIGR : DERIV21 288 : C dFr/dSIGM : DERIV22 289 : C dFr/dFVAR : DERIV23 290 : C 291 : DERIV21 = DSINH(Q2COE*SIGR/2.0D0/SIGM)* 292 : & Q1COE*Q2COE*FVAR/SIGM 293 : 294 : DERIV22 = DSINH(Q2COE*SIGR/2.0D0/SIGM)* 295 : & (-1.0D0)*Q1COE*Q2COE*FVAR*SIGR/SIGM**2.0D0 296 : 297 : DERIV23 = 2.0D0*Q1COE*DCOSH(Q2COE*SIGR/2.0D0/SIGM)- 298 : & 2.0D0*(Q3COE**2.0D0)*FVAR 299 : C 300 : C CALUL DU MODULE TANGENT COHERENTS ENTRE DSIGM ET DEPSM 301 : C 302 : ETAN = XKN/NCOEF*(SYCOE/SIGM)**(NCOEF-1.0D0) 303 : ESTA = 1.0D0/ETAN - 1.0D0/XKN 304 : ESTA = 1.0D0/ESTA 305 : C 306 : C CALCUL DU MULTIPLICATEUR DE LAGRANGE 307 : C 308 : DLAM2 = -CRIT2/(XKN*(DERIV21)**2.0D0- 309 : & DERIV23*KCOEF*FVAR*(1.0D0-FVAR)*DERIV21- 310 : & ESTA/((1.0D0-FVAR)*SIGM)*DERIV21*SIGR*DERIV22) 311 : C 312 : C ACTUALISATION DES VARIABLES FORCES 313 : C 314 : SIGR = SIGR + DLAM2 * XKN * DERIV21 315 : FVAR = FVAR - DLAM2 * KCOEF * FVAR * (1.0D0-FVAR) * DERIV21 316 : SIGM = SIGM - DLAM2 * SIGR * DERIV21 * 317 : & ESTA/(1.0D0-FVAR)/SIGM 318 : C 319 : C SECURITE SUR LA VALEUR DE FVAR (CAS OU ELLE DEVIENT TROP FAIBLE) 320 : C 321 : IF (FVAR.LE.FVARLIM) THEN 322 : FVAR = FVARLIM 323 : MODR = nint(DLAM2 * XKN * DERIV21/DEPST(2) ) 324 : END IF 325 : C 326 : C -------------------------------------------- 327 : C RETURN MAPPING : FIN DES ITERATIONS INTERNES 328 : C -------------------------------------------- 329 : C 330 : END DO 331 : C 332 : C BALISE DE SORTIE DU RETURN MAPPING AVANT ACUALISATION 333 : C 334 : 20 CONTINUE 335 : C 336 : C CAS OU LA FRACTION VOLUMIQUE DES PORES EST EGALE A 0.02 337 : C 338 : IF ((FVAR.EQ.FVARLIM).AND. 339 : & (IFLAG.EQ.1)) THEN 340 : SIGR = SIGR - MODR * DEPST(3) 341 : END IF 342 : C 343 : C CALCUL DES QUANTITES RESULTATS (EPSR, VKACT ET XKN) 344 : C 345 : SIGN = -SIGR 346 : EPSR = -(EPSN - SIGR/XKN) 347 : C 348 : C BALISE DE SORTIE DU RETURN MAPPING APRES ACUALISATION 349 : C 350 : 30 CONTINUE 351 : C 352 : C ACTUALISATION DES PARAMETRES ELASTIQUES 353 : C 354 : XKNACT = 4.0D0*XKN0*XKS0*(1.0D0-FVAR)/ 355 : & (4.0D0*XKS0+3.0D0*XKN0*FVAR) 356 : XKSACT = XKS0*(1.0D0-FVAR)/ 357 : & (1.0D0+((6.0D0*XKN0+12.0D0*XKS0)/ 358 : & (9.0D0*XKN0+8.0D0*XKS0))*FVAR) 359 : 360 : XKN = XKNACT 361 : XKS = XKSACT 362 : 363 : C 364 : C-----CALCUL DES CONTRAINTES FINALES------------------------------------ 365 : C 366 : SIGF(1) = XKS*(1-D)*EPST1 + SPI1 367 : SIGF(2) = XKS*(1-D)*EPST2 + SPI2 368 : 369 : IF (EPSN.GT.0.0D0) THEN 370 : SIGF(3) = XKN*(1-D)*EPSN 371 : ELSE 372 : SIGF(3) = XKN*(EPSN-EPSR) 373 : ENDIF 374 : C 375 : C-----STOCKAGE EN SORTIE DES VARIABLES INTERNES------------------------- 376 : C 377 : VARF(1 ) = D 378 : VARF(2 ) = EPSF1 379 : VARF(3 ) = EPSF2 380 : VARF(4 ) = XEC1 381 : VARF(5 ) = XEC2 382 : VARF(9 ) = CIN1 383 : VARF(10) = CIN2 384 : VARF(11) = EPSR 385 : VARF(12) = IFLAG 386 : VARF(13) = FVAR 387 : VARF(14) = SIGM 388 : VARF(15) = XKS 389 : VARF(16) = XKN 390 : VARF(17) = COUN + 1 391 : 392 : 123 CONTINUE 393 : 394 : IF (NGONF.EQ.1) THEN 395 : SIGF(1) = XMAT(1)*EPST1 396 : SIGF(2) = XMAT(1)*EPST2 397 : SIGF(3) = XMAT(2)*EPSN 398 : ENDIF 399 : 400 : VARF(6 ) = EPST1 401 : VARF(7 ) = EPST2 402 : VARF(8 ) = EPSN 403 : 404 : RETURN 405 : END 406 : 407 : 408 : 409 : 410 : 411 : 412 : 413 : 414 : 415 :
© Cast3M 2003 - All rights reserved.
Disclaimer