1 : C MSHEAR SOURCE PV 11/03/07 21:17:33 6885 2 : SUBROUTINE MSHEAR(DDEP,FOR0,FORF,VAR0,VARF,NVAR,DDINL, 3 : & DTRANP,DTRANM,BETA,NPELA,TRFA,DOCP,DOCM,EXPN, 4 : & CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP, 5 : & CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,KERRE) 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8(A-H,O-Z) 8 : C 9 : C====================================================================== 10 : C MUR_SHEAR, Elisa, Armelle, Alessandra et Pierre, ELSA/ISPRA 11 : C 01/2002 --> ../2004 ... 12 : C====================================================================== 13 : C 14 : C MODELE GLOBAL DE MUR EN CISAILLEMENT 15 : C (Sur des elements de poutre TIMO - Effort tranchant/Cisail.) 16 : C 17 : C====================================================================== 18 : C======================================================================= 19 : C 20 : C LISTE D'ECHANGE 21 : C --------------- 22 : C 23 : C DDEP = Incrément de déplacement 24 : C FOR0 = Effort initial 25 : C FORF = Effort final 26 : C VAR0 = Variables internes initiales 27 : C VARF = Variables internes finales 28 : C NVAR = Nb de variable interne 29 : C DDINL = Increment de deformation inelastique 30 : C 31 : C DTRANP = zone de transition (sens positif) 32 : C DTRANM = (sens négatif) 33 : C 34 : C BETA = Coefficient BETA (equivalence nrj-deplacement) 35 : C NPELA = Nb de pts par unite de zone elastique pour les 36 : C segments non lineaires 37 : C+2003 38 : C TRFA = Coefficient d'entrainement des origine (>0,<1) 39 : C DOCP = Deplacement d'ouverture des fissure + (>0) 40 : C DOCM = Deplacement d'ouverture des fissure - (<0) 41 : C+2003 42 : C+2004 43 : C EXPN = Exposant n de la fonction de cumul du domage 44 : C+2004 45 : C 46 : C CURFP,NCURFP = courbe de charge + (x>0,y>0) 47 : C CURKP,NCURKP = courbe de raideur + (x>0,y>0) 48 : C CURLP,NCURLP = courbe de domaine elastique + (x>0,y>0) 49 : C CURFM,NCURFM = courbe de charge - (x<0,y<0) 50 : C CURKM,NCURKM = courbe de raideur - (x<0,y>0) 51 : C CURLM,NCURLM = courbe de domaine elastique - (x<0,y>0) 52 : C 53 : C VARIABLES INTERNES (0=initiales, F=finales) 54 : C ------------------ 55 : C 56 : C : Numero du cas 57 : C DEP : Deplacement total 58 : C EPLUS : Energie "positive" 59 : C EMOIN : Energie "negative" 60 : C DPLUS : Deplacement de reference + 61 : C DMOIN : Deplacement de reference - 62 : C FCINI : Force de reference pour la surface cinematique 63 : C FCAMP : Taille de la surface cinematique 64 : C KCINE : Raideur de la surface cinematique 65 : C APLUS : Coefficient multiplicateur de l'amplitude + 66 : C AMOIN : Coefficient multiplicateur de l'amplitude - 67 : C+2003 68 : C OPLUS : Position de l'origine + <0 69 : C OMOIN : Position de l'origine - <0 70 : C+2003 71 : C 72 : C 73 : C======================================================================= 74 : C 75 : PARAMETER (XZER=0.D0,UN=1.D0,XPETIT=1.D-7) 76 : DIMENSION CURFP(2,NCURFP),CURKP(2,NCURKP),CURLP(2,NCURLP) 77 : DIMENSION CURFM(2,NCURFM),CURKM(2,NCURKM),CURLM(2,NCURLM) 78 : C 79 : REAL*8 KCINE0,KCINEF 80 : C 81 : DIMENSION VAR0(NVAR),VARF(NVAR) 82 : LOGICAL LSWICH 83 : C 84 : KERRE=0 85 : C 86 : C Lecture variables internes 87 : C 88 : ICAS = nint(VAR0(1)) 89 : DEP0 = VAR0(2) 90 : EPLUS0= VAR0(3) 91 : EMOIN0= VAR0(4) 92 : DPLUS0= VAR0(5) 93 : DMOIN0= VAR0(6) 94 : FCINI0= VAR0(7) 95 : FCAMP0= VAR0(8) 96 : KCINE0= VAR0(9) 97 : APLUS0= VAR0(10) 98 : AMOIN0= VAR0(11) 99 : C+2003 100 : OPLUS0= VAR0(12) 101 : OMOIN0= VAR0(13) 102 : C+2003 103 : C+2004 104 : E= VAR0(12) 105 : OMOIN0= VAR0(13) 106 : C+2004 107 : C 108 : C Initialisation des variables finales 109 : C 110 : DDINL=XZER 111 : DEPF =DEP0+DDEP 112 : DO IE1=1,NVAR 113 : VARF(IE1)=VAR0(IE1) 114 : ENDDO 115 : VARF(2)=DEPF 116 : FORF=FOR0 117 : * 118 : * WRITE(6,*)ICAS,DEP0,DEPF 119 : * 120 : C 121 : C Cas de l'increment nul 122 : C 123 : IF(DDEP.EQ.XZER)RETURN 124 : C================================================================== 125 : C 126 : C DRIVER 127 : C 128 : C================================================================== 129 : GOTO(1,2,3,4,5,6,7),ICAS+1 130 : C========================================================== 131 : C A: CAS DU MODELE DANS SON ETAT ELASTIQUE ICAS= 0 132 : C========================================================== 133 : 1 CONTINUE 134 : C-->Condition de branchement sur B0 ou C0 135 : IF(DEPF.GT.CURFP(1,2))THEN 136 : ICAS=1 137 : DEP0 = CURFP(1,2) 138 : FOR0 = CURFP(2,2) 139 : KCINE0 = FOR0/DEP0 140 : DMOIN0 = CURFM(1,2) 141 : VARF(6)= DMOIN0 142 : GOTO 2 143 : ELSEIF(DEPF.LT.CURFM(1,2))THEN 144 : ICAS=4 145 : DEP0 = CURFM(1,2) 146 : FOR0 = CURFM(2,2) 147 : KCINE0 = FOR0/DEP0 148 : DPLUS0 = CURFP(1,2) 149 : VARF(5)= DPLUS0 150 : GOTO 5 151 : ENDIF 152 : C-->Regime normal: On calcul l'etat elastique 153 : IF(DEPF.GE.0)THEN 154 : KCINEF=CURFP(2,2)/CURFP(1,2) 155 : FORF=KCINEF*DEPF 156 : ELSE 157 : KCINEF=CURFM(2,2)/CURFM(1,2) 158 : FORF=KCINEF*DEPF 159 : ENDIF 160 : RETURN 161 : C========================================================== 162 : C B0: CAS DU CHARGEMENT MONOTONE DIRECTION + ICAS= 1 163 : C========================================================== 164 : 2 CONTINUE 165 : C-->Condition de branchement sur C1 166 : IF(DDEP.LT.XZER)THEN 167 : ICAS=2 168 : GOTO 3 169 : ENDIF 170 : C-->Regime normal: On procede sur la courbe de charge + 171 : C 2003DPLUSF=DEPF 172 : DPLUSF=DEPF-OPLUS0 172 : C 2003CALL YOFXCU(DEPF ,CURFP,NCURFP, FORF,DYSDX,KERRE) 172 : CALL YOFXCU(DPLUSF,CURFP,NCURFP, FORF,DYSDX,KERRE) 173 : FORF=FORF*(UN-APLUS0) 174 : VARF(1)=ICAS 175 : VARF(5)=DPLUSF 176 : C-->On calcule et on sauve les caracteristiques de la surface 177 : C Cinematique en tenant compte de la zone de transition 177 : C 2003CALL MSHEA1(DEPF,ICAS,FORF,DPLUSF,DMOIN0, DTRANP,DTRANM, 177 : CALL MSHEA1(DPLUSF,ICAS,FORF,DPLUSF,DMOIN0, DTRANP,DTRANM, 178 : > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP, 179 : > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM, 180 : > FCINIF,FCAMPF,KCINEF) 181 : C2003> FCINIF,FCAMPF,KCINEF,ALPHP) 182 : VARF(7) =FCINIF 183 : VARF(8) =FCAMPF 184 : VARF(9) =KCINEF 185 : C-->Increment de d.p. 186 : DDINL=DDINL+DEPF-FORF/KCINEF-(DEP0-FOR0/KCINE0) 187 : C+2003 188 : IF(DPLUSF.GE.DOCP)THEN 189 : IF(OMOIN0.EQ.XZER)THEN 190 : OMOINF=TRFA*(DPLUSF-DOCP) 191 : ELSE 192 : OMOINF=OMOIN0+TRFA*(DPLUSF-DPLUS0) 193 : ENDIF 194 : VARF(13)=OMOINF 195 : ENDIF 196 : C+2003 197 : RETURN 198 : C========================================================== 199 : C C1: BRANCHE DE UNLOADING ELASTIQUE - ICAS= 2 200 : C========================================================== 201 : 3 CONTINUE 202 : C-->Condition de branchement sur B0, B2 ou C2 203 : FORF=FOR0+KCINE0*DDEP 204 : IF(FORF.GT.FCINI0)THEN 205 : C -->On sort par le haut (B0 ou B2) 206 : DEP0=DEP0+(FCINI0-FOR0)/KCINE0 207 : DDEP=DEPF-DEP0 208 : FOR0=FCINI0 209 : C -->On retourne sur B0... si on vient de B0 210 : IF(ABS(DEP0-DPLUS0).LT.XPETIT*DPLUS0)THEN 211 : ICAS=1 212 : GOTO 2 213 : C -->On retourne sur B2... si on ne vient PAS de B0 214 : ELSE 215 : ICAS=6 216 : GOTO 7 217 : ENDIF 218 : ELSEIF(FORF.LT.FCINI0-FCAMP0)THEN 219 : C --> On continue sur C2 220 : DEP0=DEP0+(FCINI0-FCAMP0-FOR0)/KCINE0 221 : DDEP=DEPF-DEP0 222 : FOR0=FCINI0-FCAMP0 223 : C 224 : C+2002 On repointe "tangent" a la courbe (si possible) 224 : CALL YOFXCU(DMOIN0,CURFM,NCURFM, FMOIN0,DYSDX,KERRE) 225 : C 2003 STRIAL=(FMOIN0-FOR0)/(DMOIN0-DEP0) 226 : STRIAL=(FMOIN0-FOR0)/(DMOIN0+OMOIN0-DEP0) 227 : DO IE1=1,NCURFM 228 : IF(DMOIN0.GT.CURFM(1,IE1))THEN 229 : C 2003 SNEXT=(CURFM(2,IE1)-FOR0)/(CURFM(1,IE1)-DEP0) 230 : SNEXT=(CURFM(2,IE1)-FOR0)/(CURFM(1,IE1)+OMOIN0-DEP0) 231 : IF(SNEXT.GT.STRIAL)THEN 232 : DMOIN0=CURFM(1,IE1) 233 : STRIAL=SNEXT 234 : ELSE 235 : GOTO 31 236 : ENDIF 237 : ENDIF 238 : ENDDO 239 : 31 CONTINUE 240 : VARF(6)=DMOIN0 241 : C+2002 242 : ICAS=3 243 : GOTO 4 244 : ENDIF 245 : C-->Regime normal: On decharge elastiquement 246 : VARF(1)=ICAS 247 : RETURN 248 : C========================================================== 249 : C C2: BRANCHE DE UNLOADING ANELASTIQUE - ICAS= 3 250 : C========================================================== 251 : 4 CONTINUE 252 : C-->Condition de branchement sur B1 253 : IF(DDEP.GT.XZER)THEN 254 : ICAS=5 255 : GOTO 6 256 : ENDIF 257 : C-->Valeur de ALPHP 258 : C 2003ALPHP=(DEP0-FOR0/KCINE0-DTRANM)/(DTRANP-DTRANM) 259 : C 2003ALPHP=MAX(XZER,ALPHP) 260 : C 2003ALPHP=MIN(UN ,ALPHP) 261 : C+2003 261 : CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0, 261 : CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0, 262 : > FOR0,FCAMP0,KCINE0,ICAS, 263 : > ALPHP,ALPHM) 264 : C+2003 265 : C-->Determination de la taille du sous increment 266 : C WARNING DDEP et XINCR sont negatifs 267 : XINCR=(CURFM(2,2)-CURFP(2,2))/NPELA 268 : NINCR=INT(DDEP/XINCR)+1 269 : XINCR=DDEP/NINCR 270 : C-->Determination initiale du point "parfait" que l'on vise (DMOIN0,FMOIN0) 270 : CALL YOFXCU(DMOIN0,CURFM,NCURFM, FMOIN0,DYSDX,KERRE) 271 : C-->Loop sur les sous-increments 272 : LSWICH=.FALSE. 273 : DO IE1=1,NINCR 274 : C-->Traitement du branchement possible sur C0 275 : C 2003 IF(DEP0+XINCR.LT.DMOIN0)THEN 276 : C IF(DEP0+XINCR.LT.DMOIN0+OMOIN0)THEN 277 : C 2003 IF(DEP0+XINCR.LT.(UN-XPETIT)*DMOIN0)THEN 278 : IF(DEP0+XINCR.LT.(UN-XPETIT)*(DMOIN0+OMOIN0))THEN 279 : LSWICH=.TRUE. 280 : C 2003 XINCR =DMOIN0-DEP0 281 : XINCR =DMOIN0+OMOIN0-DEP0 282 : ENDIF 283 : C-->On incremente explicitement tout ce qui doit etre incremente 284 : C WARNING:force effective pointee=FMOIN0*(UN-AMOIN0) 285 : C 2003 FORF=FOR0+XINCR/(DMOIN0-DEP0)*(FMOIN0*(UN-AMOIN0)-FOR0) 286 : FORF=FOR0+XINCR/(DMOIN0+OMOIN0-DEP0)*(FMOIN0*(UN-AMOIN0)-FOR0) 287 : DEP0=DEP0+XINCR 288 : DENE=ABS(FCAMP0/2*(XINCR-(FORF-FOR0)/KCINE0)) 289 : DDINL=DDINL+XINCR+FOR0/KCINE0 290 : FOR0=FORF 291 : C 292 : EPLUS0=EPLUS0+ ALPHP *DENE 293 : C 2003 EMOIN0=EMOIN0+(UN-ALPHP)*DENE 294 : EMOIN0=EMOIN0+ ALPHM *DENE 295 : * APLUS0=MIN(UN,APLUS0+BETA* ALPHP* DENE) 296 : * AMOIN0=MIN(UN,AMOIN0+BETA*(UN-ALPHP)*DENE) 297 : C 2003 APLUS0=TANH(BETA*EPLUS0) 298 : C 2003 AMOIN0=TANH(BETA*EMOIN0) 298 : CALL MSHEA3(BETA,EXPN, EPLUS0,EMOIN0, APLUS0,AMOIN0) 299 : C-->On calcule les caracteristiques de la surface 300 : C Cinematique en tenant compte de la zone de transition 300 : CALL MSHEA1(DEP0,ICAS,FOR0,DPLUS0,DMOIN0, 301 : > DTRANP,DTRANM, 302 : > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP, 303 : > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM, 304 : > FCINI0,FCAMP0,KCINE0) 305 : C2003> FCINI0,FCAMP0,KCINE0,ALPHP) 306 : C+2003 306 : CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0, 306 : CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0, 307 : > FOR0,FCAMP0,KCINE0,ICAS, 308 : > ALPHP,ALPHM) 309 : C+2003 310 : DDINL=DDINL-FOR0/KCINE0 311 : C-->Switch sur C0 312 : IF(LSWICH)THEN 313 : ICAS=4 314 : VARF(3)=EPLUS0 315 : VARF(4)=EMOIN0 316 : VARF(10)=APLUS0 317 : VARF(11)=AMOIN0 318 : GOTO 5 319 : ENDIF 320 : C-->Fin normale loop sur les sous-increments 321 : ENDDO 322 : C 323 : FORF=FOR0 324 : C 325 : VARF(1)=ICAS 326 : VARF(3)=EPLUS0 327 : VARF(4)=EMOIN0 328 : VARF(7)=FCINI0 329 : VARF(8)=FCAMP0 330 : VARF(9)=KCINE0 331 : VARF(10)=APLUS0 332 : VARF(11)=AMOIN0 333 : C 334 : RETURN 335 : C========================================================== 336 : C B0: CAS DU CHARGEMENT MONOTONE DIRECTION - ICAS= 4 337 : C========================================================== 338 : 5 CONTINUE 339 : C-->Condition de branchement sur B1 340 : IF(DDEP.GT.XZER)THEN 341 : ICAS=5 342 : GOTO 6 343 : ENDIF 344 : C-->Regime normal: On procede sur la courbe de charge - 345 : C 2003DMOINF=DEPF 346 : DMOINF=DEPF-OMOIN0 346 : C 2003CALL YOFXCU(DEPF ,CURFM,NCURFM, FORF,DYSDX,KERRE) 346 : CALL YOFXCU(DMOINF,CURFM,NCURFM, FORF,DYSDX,KERRE) 347 : FORF=FORF*(UN-AMOIN0) 348 : VARF(1)=ICAS 349 : VARF(6)=DMOINF 350 : C-->On calcule et on sauve les caracteristiques de la surface 351 : C Cinematique en tenant compte de la zone de transition 351 : C 2003CALL MSHEA1(DEPF,ICAS,FORF,DPLUS0,DMOINF, DTRANP,DTRANM, 351 : CALL MSHEA1(DEPF,ICAS,FORF,DPLUS0,DMOINF, DTRANP,DTRANM, 352 : > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP, 353 : > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM, 354 : > FCINIF,FCAMPF,KCINEF) 355 : C2003> FCINIF,FCAMPF,KCINEF,ALPHP) 356 : VARF(7) =FCINIF 357 : VARF(8) =FCAMPF 358 : VARF(9) =KCINEF 359 : C-->Increment de d.p. 360 : DDINL=DDINL+DEPF-FORF/KCINEF-(DEP0-FOR0/KCINE0) 361 : C+2003 362 : IF(DMOINF.LE.DOCM)THEN 363 : IF(OPLUS0.EQ.XZER)THEN 364 : OPLUSF=TRFA*(DMOINF-DOCM) 365 : ELSE 366 : OPLUSF=OPLUS0+TRFA*(DMOINF-DMOIN0) 367 : ENDIF 368 : VARF(12)=OPLUSF 369 : ENDIF 370 : C+2003 371 : RETURN 372 : C========================================================== 373 : C B1: BRANCHE DE RELOADING ELASTIQUE + ICAS= 5 374 : C========================================================== 375 : 6 CONTINUE 376 : C-->Condition de branchement sur C0, C2 ou B2 377 : FORF=FOR0+KCINE0*DDEP 378 : IF(FORF.LT.FCINI0)THEN 379 : C -->On sort par le bas (C0 ou C2) 380 : DEP0=DEP0+(FCINI0-FOR0)/KCINE0 381 : DDEP=DEPF-DEP0 382 : FOR0=FCINI0 383 : C -->On retourne sur C0... si on vient de C0 384 : IF(ABS(DEP0-DMOIN0).LT.ABS(XPETIT*DMOIN0))THEN 385 : ICAS=4 386 : GOTO 5 387 : C -->On retourne sur C2... si on ne vient PAS de C0 388 : ELSE 389 : ICAS=3 390 : GOTO 4 391 : ENDIF 392 : ELSEIF(FORF.GT.FCINI0+FCAMP0)THEN 393 : C --> On continue sur B2 394 : DEP0=DEP0+(FCINI0+FCAMP0-FOR0)/KCINE0 395 : DDEP=DEPF-DEP0 396 : FOR0=FCINI0+FCAMP0 397 : C 398 : C2002 On repointe "tangent" a la courbe (si possible) 398 : CALL YOFXCU(DPLUS0,CURFP,NCURFP, FPLUS0,DYSDX,KERRE) 399 : C 2003 STRIAL=(FPLUS0-FOR0)/(DPLUS0-DEP0) 400 : STRIAL=(FPLUS0-FOR0)/(DPLUS0+OPLUS0-DEP0) 401 : DO IE1=1,NCURFP 402 : IF(DPLUS0.LT.CURFP(1,IE1))THEN 403 : C 2003 SNEXT=(CURFP(2,IE1)-FOR0)/(CURFP(1,IE1)-DEP0) 404 : SNEXT=(CURFP(2,IE1)-FOR0)/(CURFP(1,IE1)+OPLUS0-DEP0) 405 : IF(SNEXT.GT.STRIAL)THEN 406 : DPLUS0=CURFP(1,IE1) 407 : STRIAL=SNEXT 408 : ELSE 409 : GOTO 61 410 : ENDIF 411 : ENDIF 412 : ENDDO 413 : 61 CONTINUE 414 : VARF(5)=DPLUS0 415 : C2002 416 : ICAS=6 417 : GOTO 7 418 : ENDIF 419 : C-->Regime normal: On decharge elastiquement 420 : VARF(1)=ICAS 421 : RETURN 422 : C========================================================== 423 : C B2: BRANCHE DE UNLOADING ANELASTIQUE - ICAS= 6 424 : C========================================================== 425 : 7 CONTINUE 426 : C-->Condition de branchement sur C1 427 : IF(DDEP.LT.XZER)THEN 428 : ICAS=2 429 : GOTO 3 430 : ENDIF 431 : C-->Valeur de ALPHP 432 : C 2003ALPHP=(DEP0-FOR0/KCINE0-DTRANM)/(DTRANP-DTRANM) 433 : C 2003ALPHP=MAX(XZER,ALPHP) 434 : C 2003ALPHP=MIN(UN ,ALPHP) 435 : C+2003 435 : CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0, 435 : CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0, 436 : > FOR0,FCAMP0,KCINE0,ICAS, 437 : > ALPHP,ALPHM) 438 : C+2003 439 : 440 : C-->Determination de la taille du sous increment 441 : C WARNING DDEP et XINCR sont negatifs 442 : XINCR=(CURFP(2,2)-CURFM(2,2))/NPELA 443 : NINCR=INT(DDEP/XINCR)+1 444 : XINCR=DDEP/NINCR 445 : C-->Determination initiale du point "parfait" que l'on vise (DPLUS0,FPLUS0) 445 : CALL YOFXCU(DPLUS0,CURFP,NCURFP, FPLUS0,DYSDX,KERRE) 446 : C-->Loop sur les sous-increments 447 : LSWICH=.FALSE. 448 : DO IE1=1,NINCR 449 : C-->Traitement du branchement possible sur B0 450 : C 2003 IF(DEP0+XINCR.GT.DPLUS0)THEN 451 : C IF(DEP0+XINCR.GT.DPLUS0+OPLUS0)THEN 452 : C 2003 IF(DEP0+XINCR.GT.(UN-XPETIT)*DPLUS0)THEN 453 : IF(DEP0+XINCR.GT.(UN-XPETIT)*(DPLUS0+OPLUS0))THEN 454 : LSWICH=.TRUE. 455 : C 2003 XINCR =DPLUS0-DEP0 456 : XINCR =DPLUS0+OPLUS0-DEP0 457 : ENDIF 458 : C-->On incremente explicitement tout ce qui doit etre incremente 459 : C WARNING:l'increment de d.p. est calcule en 2 steps 460 : C WARNING:force effective pointee=FPLUS0*(UN-APLUS0) 461 : C 2003 FORF=FOR0+XINCR/(DPLUS0-DEP0)*(FPLUS0*(UN-APLUS0)-FOR0) 462 : FORF=FOR0+XINCR/(DPLUS0+OPLUS0-DEP0)*(FPLUS0*(UN-APLUS0)-FOR0) 463 : DEP0=DEP0+XINCR 464 : DENE=ABS(FCAMP0/2*(XINCR-(FORF-FOR0)/KCINE0)) 465 : DDINL=DDINL+XINCR+FOR0/KCINE0 466 : FOR0=FORF 467 : C 468 : EPLUS0=EPLUS0+ ALPHP *DENE 469 : C 2003 EMOIN0=EMOIN0+(UN-ALPHP)*DENE 470 : EMOIN0=EMOIN0+ ALPHM *DENE 471 : * APLUS0=MIN(UN,APLUS0+BETA* ALPHP* DENE) 472 : * AMOIN0=MIN(UN,AMOIN0+BETA*(UN-ALPHP)*DENE) 473 : C 2003 APLUS0=TANH(BETA*EPLUS0) 474 : C 2003 AMOIN0=TANH(BETA*EMOIN0) 474 : CALL MSHEA3(BETA,EXPN, EPLUS0,EMOIN0, APLUS0,AMOIN0) 475 : C-->On calcule les caracteristiques de la surface 476 : C Cinematique en tenant compte de la zone de transition 476 : CALL MSHEA1(DEP0,ICAS,FOR0,DPLUS0,DMOIN0, 477 : > DTRANP,DTRANM, 478 : > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP, 479 : > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM, 480 : > FCINI0,FCAMP0,KCINE0) 481 : C2003> FCINI0,FCAMP0,KCINE0,ALPHP) 482 : C+2003 482 : CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0, 482 : CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0, 483 : > FOR0,FCAMP0,KCINE0,ICAS, 484 : > ALPHP,ALPHM) 485 : C+2003 486 : DDINL=DDINL-FOR0/KCINE0 487 : C-->Switch sur B0 488 : IF(LSWICH)THEN 489 : ICAS=1 490 : VARF(3)=EPLUS0 491 : VARF(4)=EMOIN0 492 : VARF(6)=DMOIN0 493 : VARF(10)=APLUS0 494 : VARF(11)=AMOIN0 495 : GOTO 2 496 : ENDIF 497 : C-->Fin normale loop sur les sous-increments 498 : ENDDO 499 : C 500 : FORF=FOR0 501 : C 502 : VARF(1)=ICAS 503 : VARF(3)=EPLUS0 504 : VARF(4)=EMOIN0 505 : VARF(7)=FCINI0 506 : VARF(8)=FCAMP0 507 : VARF(9)=KCINE0 508 : VARF(10)=APLUS0 509 : VARF(11)=AMOIN0 510 : C 511 : RETURN 512 : C 513 : END 514 : 515 : 516 :
© Cast3M 2003 - All rights reserved.
Disclaimer