1 : C C1DP SOURCE CHAT 05/01/12 21:44:47 5004 2 : SUBROUTINE C1DP(SIG,DSIG,YOUN,ANU,RT1,RT2,RDP,ADP,HDP, 3 : 1 XLAM1,XLAM2,XLAM3,IDAM,ANG,EPPLDP,ALPHA,KERRE) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8(A-H,O-Z) 6 : -INC CCOPTIO 7 : DIMENSION SIG(3),DSIG(3),SGG(3),DSGG(3),SFG(3),DSFG(3), 8 : 1 A(3),B(3),IDAM(3),EPC(3),SSI(3),SI(3),EPC1(3), 9 : 1 EPDP(3),EC(3),EPPLDP(3),SG(3),AT(4),BT(2),TENS(3) 10 : ITER=0 11 : Y=YOUN/(1.D0-ANU*ANU) 12 : SIREF=1.E-6*YOUN 13 : DO 10 ITYP=1,3 14 : EPDP(ITYP)=0.D0 15 : 10 EPPLDP(ITYP)=0.D0 16 : DL3=0.D0 17 : XT=0.D0 18 : IF(IIMPI.EQ.9) WRITE(IOIMP,9999) 19 : 9999 FORMAT(1X,'C1DP COUPLAGE COMP 1 DRUCKER',/) 20 : 9998 FORMAT(1X,'C1DP ',I4,'ITERATIONS INTERNES',/) 21 : C 22 : C------------------------------------------------ 23 : C COUPLAGE COMPRESSION 1 DRUCKER 24 : C------------------------------------------------ 25 : C---------------------------------------------- 26 : C DEFINITION D UNE CONTRAINTE DE REFERENCE 27 : C SI RMAX DU DRUCKER INFERIEURE A CETTE VALEUR 28 : C ALORS SIGMA=0 29 : C---------------------------------------------- 30 : C 31 : RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) 32 : IF(SIREF.GT.RMAX) THEN 33 : DO 20 ITYP=1,3 34 : 20 SIG(ITYP)=SIG(ITYP)+DSIG(ITYP) 35 : RDP=0.D0 36 : HDP=0.D0 37 : ADP=0.D0 37 : CALL CPHOIN(SIG,EPPLDP,YOUN,ANU) 37 : CALL NORME(SIG,DL3) 38 : DL3=DL3/YOUN*10.D0 39 : XLAM3=XLAM3+DL3 40 : DO 30 ITYP=1,3 41 : IDAM(ITYP)=0 42 : DSIG(ITYP)=0.D0 43 : 30 SIG(ITYP)=0.D0 44 : IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER 45 : RETURN 46 : ENDIF 47 : C----------------------------------------------- 48 : C ON SEPLACE DANS LE REPERE DE FISSURATION 49 : C----------------------------------------------- 50 : C 50 : CALL CHREP(SIG,SFG,ANG) 50 : CALL CHREP(DSIG,DSFG,ANG) 51 : C 52 : C----------------------------------------------------- 53 : C ESTIMATION DU PAS D'INCREMENT DE CONTRAINTES 54 : C----------------------------------------------------- 55 : C 55 : CALL SIJ(SFG,SI,SIETOI) 55 : CALL EPSDP(SI,SIETOI,ADP,EPC) 55 : CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA) 55 : CALL SCAL2(EPC,SSI,VAL) 56 : VAL=VAL-HDP 57 : IF(VAL.LT.0.D0) THEN 58 : KERRE=459 59 : RETURN 60 : ENDIF 60 : CALL SCAL(SSI,DSFG,VAL) 60 : CALL NORME(DSFG,VA1) 60 : CALL NORME(SSI,VA2) 61 : IF(VA1.EQ.0.D0) THEN 62 : IDAM(1)=0 63 : IDAM(2)=0 64 : IDAM(3)=0 65 : RETURN 66 : ENDIF 67 : RMIN=MIN((RDP/1.73),(RDP/(1.D0+2.D0*ADP))) 68 : X=VAL/VA1/VA2 69 : IF(X.GE.1.D0) THEN 70 : X=1.D0 71 : GOTO 41 72 : ENDIF 73 : X=1.D0/SQRT(1.0001D0-X*X)*RMIN/VA1/8.D0 74 : IF(X.GE.1.D0) X=1.D0 75 : C 76 : C---------------------------------------------------------- 77 : C ON ECOULE ET ON REGARDE LA VARIATION DE NORMALE 78 : C---------------------------------------------------------- 79 : C 80 : C 81 : C--------------------------------------- 82 : C RESOLUTION DU SYSTEME EN DL3 83 : C--------------------------------------- 84 : 41 X=X/2.D0 85 : EC(1)=1.D0 86 : EC(2)=0.D0 87 : EC(3)=0.D0 87 : CALL CPHOOB(EC,SG,YOUN,ANU) 88 : 100 X=X*2.D0 88 : CALL SIJ(SFG,SI,SIETOI) 88 : CALL EPSDP(SI,SIETOI,ADP,EPC) 88 : CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA) 88 : CALL SCAL2(EPC,SSI,VAL) 89 : VAL=VAL-HDP 90 : IF(VAL.LT.0.D0) THEN 91 : KERRE=459 92 : RETURN 93 : ENDIF 94 : IF((XT+X).GT.1.D0) X=1.D0-XT 95 : 110 RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) 96 : C 97 : C -------------------------------------------------- 98 : C CAS OU LE RAYON MAXIMAL DU DRUCKER-PRAGER 99 : C INFERIEUR A LA CONTRAINTE DE REFERENCE 100 : C -------------------------------------------------- 101 : C 102 : IF(SIREF.LT.RMAX) GOTO 131 103 : DO 120 ITYP=1,3 104 : 120 SFG(ITYP)=SFG(ITYP)+(1.D0-XT)*DSFG(ITYP) 105 : RDP=0.D0 106 : ADP=0.D0 107 : HDP=0.D0 107 : CALL NORME(SFG,DL3) 108 : DL3=DL3/YOUN*10.D0 109 : XLAM3=XLAM3+DL3 109 : CALL CPHOIN(SFG,EPPLDP,YOUN,ANU) 110 : DO 130 ITYP=1,3 111 : IDAM(ITYP)=0 112 : EPDP(ITYP)=EPDP(ITYP)+EPPLDP(ITYP) 113 : DSIG(ITYP)=0.D0 114 : 130 SIG(ITYP)=0.D0 114 : CALL CHREP(EPDP,EPPLDP,-ANG) 115 : IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER 116 : RETURN 117 : 131 ITER=ITER+1 118 : IF(ITER.GT.200) THEN 119 : KERRE=460 120 : RETURN 121 : ENDIF 122 : C 123 : C--------------------------------- 124 : C ON REALISE LE COUPLAGE 125 : C--------------------------------- 126 : C 127 : E=-(SFG(1)+X*DSFG(1))/SG(1) 128 : F=SSI(1)/SG(1) 129 : DO 140 ITYP=1,3 130 : DSGG(ITYP)=X*DSFG(ITYP) 131 : A(ITYP)=SFG(ITYP)+DSGG(ITYP)+E*SG(ITYP) 132 : 140 B(ITYP)=-SSI(ITYP)+F*SG(ITYP) 132 : CALL DLAMDP(A,B,DL3,RDP,ADP,HDP,ITEST) 133 : IF(ITEST.EQ.1) THEN 134 : X=X/2.D0 135 : GOTO 131 136 : ENDIF 137 : DL1=-E-F*DL3 138 : DO 145 ITYP=1,3 139 : 145 SGG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) 140 : C 141 : C---------------------------------------------- 142 : C SI UN DES DLAMDA < 0 FAUX COUPLAGE 143 : C---------------------------------------------- 144 : C 145 : IF(DL1.GT.0.D0) IDAM(1)=0 146 : IF(DL3.LT.0.D0) IDAM(3)=0 147 : IF(IDAM(1).EQ.0.OR.IDAM(3).EQ.0) THEN 148 : DL1=0.D0 149 : DL3=0.D0 150 : X=0.D0 151 : GOTO 3000 152 : ENDIF 153 : C 154 : C----------------------------------------------------- 155 : C ON VERIFIE QUE LA NORMALE SUIVANT LE DRUCKER 156 : C PRAGER TOURNE PEU 157 : C----------------------------------------------------- 158 : C 158 : CALL SIJ(SGG,SI,SIETOI) 158 : CALL EPSDP(SI,SIETOI,ADP,EPC1) 158 : CALL SCAL(EPC,EPC1,VAL) 158 : CALL NORME(EPC,VA1) 158 : CALL NORME(EPC1,VA2) 159 : CO=VAL/VA1/VA2 160 : IF(CO.LT.0.99) THEN 161 : X=X/2.D0 162 : GOTO 131 163 : ENDIF 164 : IF((DL1+XLAM1).GE.0.D0) GOTO 200 165 : C 166 : C------------------------------------------------ 167 : C CAS OU LA FISSURE EST COMPLETEMENT FERMEE 168 : C------------------------------------------------ 169 : C 170 : IDAM(1)=0 171 : DL1=-XLAM1 172 : IF(DSFG(1).EQ.0.D0) GOTO 160 173 : E=(DL1*SG(1)-SFG(1))/DSFG(1) 174 : F=SSI(1)/DSFG(1) 175 : DO 150 ITYP=1,3 176 : A(ITYP)=SFG(ITYP)-DL1*SG(ITYP)+E*DSFG(ITYP) 177 : 150 B(ITYP)=F*DSFG(ITYP)-SSI(ITYP) 177 : CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST) 178 : DO 151 ITYP=1,3 179 : SGG(ITYP)=A(ITYP)+DL31*B(ITYP) 180 : 151 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP) 180 : CALL SIJ(SGG,SI,SIETOI) 180 : CALL EPSDP(SI,SIETOI,ADP,EPC1) 180 : CALL SCAL(EPC,EPC1,VAL) 180 : CALL NORME(EPC,VA1) 180 : CALL NORME(EPC1,VA2) 181 : CO1=VAL/VA1/VA2 181 : CALL SIJ(DSGG,SI,SIETOI) 181 : CALL EPSDP(SI,SIETOI,ADP,EPC1) 181 : CALL SCAL(EPC,EPC1,VAL) 181 : CALL NORME(EPC,VA1) 181 : CALL NORME(EPC1,VA2) 182 : CO2=VAL/VA1/VA2 183 : XIN1=E+F*DL31 184 : XIN2=E+F*DL32 185 : IF(XIN1.GT.-1.D-10.AND.XIN1.LE.X.AND.CO1.GT.0.9) THEN 186 : DL3=DL31 187 : X=XIN1 188 : GOTO 200 189 : ENDIF 190 : IF(XIN2.GT.-1.D-10.AND.XIN2.LE.X.AND.CO2.GT.0.9) THEN 191 : DL3=DL32 192 : X=XIN2 193 : GOTO 200 194 : ENDIF 195 : IF(IIMPI.EQ.9) WRITE(IOIMP,10101) 196 : 10101 FORMAT(1X,'ERREUR DANS C1DP FERM',/) 197 : DL3=DL31 198 : X=XIN1 199 : GOTO 200 200 : 160 DL3=(SFG(1)-DL1*SG(1))/SSI(1) 201 : RD=RDP-HDP*DL3 202 : IF(RD.LT.0.D0) THEN 203 : RDP=0.D0 204 : GOTO 110 205 : ENDIF 206 : DO 170 ITYP=1,3 207 : A(ITYP)=SFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) 208 : 170 B(ITYP)=DSFG(ITYP) 208 : CALL XDP(A,B,RD,ADP,X,ITEST) 209 : 200 CONTINUE 210 : DO 205 ITYP=1,3 211 : DSGG(ITYP)=X*DSFG(ITYP) 212 : 205 SGG(ITYP)=SFG(ITYP)+DSGG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) 213 : C 214 : C----------------------------------------------- 215 : C ON REGARDE SI ON ENDOMMAGE PAS LES AUTRES 216 : C CRITERES 217 : C----------------------------------------------- 218 : C------------------------------------- 219 : C ON ENDOMMAGE LA TRACTION 2 220 : C------------------------------------- 221 : C 221 : CALL CTRAF(SGG(2),RT2,VCTR) 222 : IF(VCTR.LE.0.D0) GOTO 1000 223 : IDAM(1)=-1 224 : IDAM(2)=1 225 : IDAM(3)=1 226 : AT(1)=-DSFG(1) 227 : AT(2)=Y 228 : AT(4)=Y*ANU 229 : AT(3)=-DSFG(2) 230 : BT(1)=SFG(1) 231 : BT(2)=SFG(2)-RT2 231 : CALL SYLIN2(AT,BT,X1,DL11) 232 : BT(1)=-SSI(1) 233 : BT(2)=-SSI(2) 233 : CALL SYLIN2(AT,BT,X2,DL12) 234 : DO 230 ITYP=1,3 235 : A(ITYP)=SFG(ITYP)-DL11*SG(ITYP)+X1*DSFG(ITYP) 236 : 230 B(ITYP)=-DL12*SG(ITYP)+X2*DSFG(ITYP)-SSI(ITYP) 236 : CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST) 237 : DO 231 ITYP=1,3 238 : SGG(ITYP)=A(ITYP)+DL31*B(ITYP) 239 : 231 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP) 239 : CALL SIJ(SGG,SI,SIETOI) 239 : CALL EPSDP(SI,SIETOI,ADP,EPC1) 239 : CALL SCAL(EPC,EPC1,VAL) 239 : CALL NORME(EPC,VA1) 239 : CALL NORME(EPC1,VA2) 240 : CO1=VAL/VA1/VA2 240 : CALL SIJ(DSGG,SI,SIETOI) 240 : CALL EPSDP(SI,SIETOI,ADP,EPC1) 240 : CALL SCAL(EPC,EPC1,VAL) 240 : CALL NORME(EPC,VA1) 240 : CALL NORME(EPC1,VA2) 241 : CO2=VAL/VA1/VA2 242 : DLIN1=DL11+DL31*DL12 243 : DLIN2=DL11+DL32*DL12 244 : XIN1=X1+DL31*X2 245 : XIN2=X1+DL32*X2 246 : IF(XIN1.GT.-1.D-10.AND.DLIN1.LT.1.E-10.AND.XIN1.LE.X 247 : 1 .AND.CO1.GT.0.9) THEN 248 : DL3=DL31 249 : X=XIN1 250 : DL1=DLIN1 251 : GOTO 3000 252 : ENDIF 253 : IF(XIN2.GT.-1.D-10.AND.DLIN2.LT.1.E-10.AND.XIN2.LE.X 254 : 1 .AND.CO2.GT.0.9) THEN 255 : DL3=DL32 256 : X=XIN2 257 : DL1=DLIN2 258 : GOTO 3000 259 : ENDIF 260 : IF(IIMPI.EQ.9) WRITE(IOIMP,20202) 261 : 20202 FORMAT(1X,'ERREUR DANS C1DP TRAC ',/) 262 : DL3=DL31 263 : DL1=DLIN1 264 : X=XIN1 265 : GOTO 3000 266 : C 267 : C------------------------------------- 268 : C ON ENDOMMAGE LA COMPRESSION 2 269 : C------------------------------------- 270 : C 270 : 1000 CALL CCOAF(SGG(2),XLAM2,VCCO2) 271 : IF(VCCO2.LE.0.D0) GOTO 2000 272 : IDAM(1)=-1 273 : IDAM(3)=1 274 : IDAM(2)=-1 275 : AT(1)=-DSFG(1) 276 : AT(2)=Y 277 : AT(4)=Y*ANU 278 : AT(3)=-DSFG(2) 279 : BT(1)=SFG(1) 280 : BT(2)=SFG(2) 280 : CALL SYLIN2(AT,BT,X1,DL11) 281 : BT(1)=-SSI(1) 282 : BT(2)=-SSI(2) 282 : CALL SYLIN2(AT,BT,X2,DL12) 283 : DO 1010 ITYP=1,3 284 : A(ITYP)=SFG(ITYP)-DL11*SG(ITYP)+X1*DSFG(ITYP) 285 : 1010 B(ITYP)=-DL12*SG(ITYP)+X2*DSFG(ITYP)-SSI(ITYP) 285 : CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST) 286 : DO 1011 ITYP=1,3 287 : SGG(ITYP)=A(ITYP)+DL31*B(ITYP) 288 : 1011 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP) 288 : CALL SIJ(SGG,SI,SIETOI) 288 : CALL EPSDP(SI,SIETOI,ADP,EPC1) 288 : CALL SCAL(EPC,EPC1,VAL) 288 : CALL NORME(EPC,VA1) 288 : CALL NORME(EPC1,VA2) 289 : CO1=VAL/VA1/VA2 289 : CALL SIJ(DSGG,SI,SIETOI) 289 : CALL EPSDP(SI,SIETOI,ADP,EPC1) 289 : CALL SCAL(EPC,EPC1,VAL) 289 : CALL NORME(EPC,VA1) 289 : CALL NORME(EPC1,VA2) 290 : CO2=VAL/VA1/VA2 291 : DLIN1=DL11+DL31*DL12 292 : DLIN2=DL11+DL32*DL12 293 : XIN1=X1+DL31*X2 294 : XIN2=X1+DL32*X2 295 : IF(XIN1.GT.-1.D-10.AND.DLIN1.LT.1.E-10.AND.XIN1.LE.X 296 : 1 .AND.CO1.GT.0.9) THEN 297 : DL3=DL31 298 : X=XIN1 299 : DL1=DLIN1 300 : GOTO 3000 301 : ENDIF 302 : IF(XIN2.GT.-1.D-10.AND.DLIN2.LT.1.E-10.AND.XIN2.LE.X 303 : 1 .AND.CO2.GT.0.9) THEN 304 : DL3=DL32 305 : X=XIN2 306 : DL1=DLIN2 307 : GOTO 3000 308 : ENDIF 309 : IF(IIMPI.EQ.9) WRITE(IOIMP,30303) 310 : 30303 FORMAT(1X,'ERREUR DANS C1DP COMP ',/) 311 : DL3=DL31 312 : DL1=DLIN1 313 : X=XIN1 314 : GOTO 3000 315 : 2000 CONTINUE 316 : C 317 : C------------------------------------------------ 318 : C CAS OU AUCUN AUTRE CRITERE EST ENDOMMAGE 319 : C------------------------------------------------ 320 : C----------------------------------- 321 : C REMISE A JOUR DES VARIABLES 322 : C PUIS NOUVEAU PAS SI TOUT 323 : C L INCREMENT N EST PAS ECOULE 324 : C----------------------------------- 325 : C 326 : XT=XT+X 327 : RDP=RDP-DL3*HDP 328 : RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) 329 : IF(SIREF.GT.RMAX) THEN 330 : RDP=0.D0 331 : GOTO 110 332 : ENDIF 333 : DO 2010 ITYP=1,3 334 : SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) 335 : 2010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP) 336 : XLAM1=XLAM1+DL1 337 : XLAM3=XLAM3+DL3 338 : IF(XLAM1.LT.1.E-10) XLAM1=0.D0 339 : IF(XLAM1.LT.1.E-10) IDAM(1)=0 340 : C=1.D0-1.D-10 341 : IF(XT.LT.C) GOTO 2020 342 : DO 2015 ITYP=1,3 343 : IDAM(ITYP)=0 344 : 2015 DSIG(ITYP)=0.D0 344 : CALL CHREP(SFG,SIG,-ANG) 344 : CALL CHREP(EPDP,EPPLDP,-ANG) 345 : IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER 346 : RETURN 347 : 2020 IF(IDAM(1).EQ.-1) GOTO 100 348 : XLAM1=0.D0 349 : GAMDP=10.D0 350 : GAMTR1=10.D0 351 : GAMTR2=10.D0 352 : GAMCO2=10.D0 353 : DO 2030 ITYP=1,3 354 : DSFG(ITYP)=DSFG(ITYP)*(1.D0-XT) 355 : 2030 TENS(ITYP)=SFG(ITYP)+DSFG(ITYP) 355 : CALL CTRAF(TENS(1),RT1,VCTR1) 355 : CALL CTRAF(TENS(2),RT2,VCTR2) 355 : CALL CCOAF(TENS(2),XLAM2,VCCO2) 355 : CALL CDP(TENS,ADP,RDP,VCDP) 355 : IF(VCTR1.GT.0.D0) CALL GAMTAF(SFG(1),DSFG(1),RT1,GAMTR1) 355 : IF(VCTR2.GT.0.D0) CALL GAMTAF(SFG(2),DSFG(2),RT2,GAMTR2) 355 : IF(VCCO2.GT.0.D0) CALL GAMCAF(SFG(2),DSFG(2),GAMCO2) 355 : IF(VCDP.GT.0.D0) CALL GDP(SFG,DSFG,RDP,ADP,GAMDP) 356 : IDAM(1)=0 357 : IDAM(2)=0 358 : IDAM(3)=0 359 : GAM=MIN(GAMTR1,GAMTR2,GAMCO2,GAMDP) 360 : IF(GAM.GE.1.D0) THEN 361 : DO 2040 ITYP=1,3 362 : SFG(ITYP)=SFG(ITYP)+DSFG(ITYP) 363 : 2040 DSIG(ITYP)=0.D0 363 : CALL CHREP(SFG,SIG,-ANG) 364 : RETURN 365 : ENDIF 366 : IF(ABS(GAM-GAMTR1).LE.1.E-10) IDAM(1)=1 367 : IF(ABS(GAM-GAMTR2).LE.1.E-10) IDAM(2)=1 368 : IF(ABS(GAM-GAMCO2).LE.1.E-10) IDAM(2)=-1 369 : IF(ABS(GAM-GAMDP).LE.1.E-10) IDAM(3)=1 370 : DO 2050 ITYP=1,3 371 : SFG(ITYP)=SFG(ITYP)-GAM*DSFG(ITYP) 372 : 2050 DSFG(ITYP)=(1.D0-GAM)*DSFG(ITYP) 372 : CALL CHREP(SFG,SIG,-ANG) 372 : CALL CHREP(DSFG,DSIG,-ANG) 372 : CALL CHREP(EPDP,EPPLDP,-ANG) 373 : RETURN 374 : 3000 RDP=RDP-DL3*HDP 375 : RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) 376 : IF(SIREF.GT.RMAX) THEN 377 : RDP=0.D0 378 : GOTO 110 379 : ENDIF 380 : XT=XT+X 381 : DO 3010 ITYP=1,3 382 : SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) 383 : DSFG(ITYP)=(1.D0-XT)*DSFG(ITYP) 384 : 3010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP) 384 : CALL CHREP(SFG,SIG,-ANG) 384 : CALL CHREP(DSFG,DSIG,-ANG) 384 : CALL CHREP(EPDP,EPPLDP,-ANG) 385 : XLAM1=XLAM1+DL1 386 : XLAM3=XLAM3+DL3 387 : IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER 388 : RETURN 389 : END 390 :
© Cast3M 2003 - All rights reserved.
Disclaimer