C C1DP SOURCE CHAT 05/01/12 21:44:47 5004 SUBROUTINE C1DP(SIG,DSIG,YOUN,ANU,RT1,RT2,RDP,ADP,HDP, 1 XLAM1,XLAM2,XLAM3,IDAM,ANG,EPPLDP,ALPHA,KERRE) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO DIMENSION SIG(3),DSIG(3),SGG(3),DSGG(3),SFG(3),DSFG(3), 1 A(3),B(3),IDAM(3),EPC(3),SSI(3),SI(3),EPC1(3), 1 EPDP(3),EC(3),EPPLDP(3),SG(3),AT(4),BT(2),TENS(3) ITER=0 Y=YOUN/(1.D0-ANU*ANU) SIREF=1.E-6*YOUN DO 10 ITYP=1,3 EPDP(ITYP)=0.D0 10 EPPLDP(ITYP)=0.D0 DL3=0.D0 XT=0.D0 IF(IIMPI.EQ.9) WRITE(IOIMP,9999) 9999 FORMAT(1X,'C1DP COUPLAGE COMP 1 DRUCKER',/) 9998 FORMAT(1X,'C1DP ',I4,'ITERATIONS INTERNES',/) C C------------------------------------------------ C COUPLAGE COMPRESSION 1 DRUCKER C------------------------------------------------ C---------------------------------------------- C DEFINITION D UNE CONTRAINTE DE REFERENCE C SI RMAX DU DRUCKER INFERIEURE A CETTE VALEUR C ALORS SIGMA=0 C---------------------------------------------- C RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) IF(SIREF.GT.RMAX) THEN DO 20 ITYP=1,3 20 SIG(ITYP)=SIG(ITYP)+DSIG(ITYP) RDP=0.D0 HDP=0.D0 ADP=0.D0 CALL CPHOIN(SIG,EPPLDP,YOUN,ANU) CALL NORME(SIG,DL3) DL3=DL3/YOUN*10.D0 XLAM3=XLAM3+DL3 DO 30 ITYP=1,3 IDAM(ITYP)=0 DSIG(ITYP)=0.D0 30 SIG(ITYP)=0.D0 IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER RETURN ENDIF C----------------------------------------------- C ON SEPLACE DANS LE REPERE DE FISSURATION C----------------------------------------------- C CALL CHREP(SIG,SFG,ANG) CALL CHREP(DSIG,DSFG,ANG) C C----------------------------------------------------- C ESTIMATION DU PAS D'INCREMENT DE CONTRAINTES C----------------------------------------------------- C CALL SIJ(SFG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC) CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA) CALL SCAL2(EPC,SSI,VAL) VAL=VAL-HDP IF(VAL.LT.0.D0) THEN KERRE=459 RETURN ENDIF CALL SCAL(SSI,DSFG,VAL) CALL NORME(DSFG,VA1) CALL NORME(SSI,VA2) IF(VA1.EQ.0.D0) THEN IDAM(1)=0 IDAM(2)=0 IDAM(3)=0 RETURN ENDIF RMIN=MIN((RDP/1.73),(RDP/(1.D0+2.D0*ADP))) X=VAL/VA1/VA2 IF(X.GE.1.D0) THEN X=1.D0 GOTO 41 ENDIF X=1.D0/SQRT(1.0001D0-X*X)*RMIN/VA1/8.D0 IF(X.GE.1.D0) X=1.D0 C C---------------------------------------------------------- C ON ECOULE ET ON REGARDE LA VARIATION DE NORMALE C---------------------------------------------------------- C C C--------------------------------------- C RESOLUTION DU SYSTEME EN DL3 C--------------------------------------- 41 X=X/2.D0 EC(1)=1.D0 EC(2)=0.D0 EC(3)=0.D0 CALL CPHOOB(EC,SG,YOUN,ANU) 100 X=X*2.D0 CALL SIJ(SFG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC) CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA) CALL SCAL2(EPC,SSI,VAL) VAL=VAL-HDP IF(VAL.LT.0.D0) THEN KERRE=459 RETURN ENDIF IF((XT+X).GT.1.D0) X=1.D0-XT 110 RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) C C -------------------------------------------------- C CAS OU LE RAYON MAXIMAL DU DRUCKER-PRAGER C INFERIEUR A LA CONTRAINTE DE REFERENCE C -------------------------------------------------- C IF(SIREF.LT.RMAX) GOTO 131 DO 120 ITYP=1,3 120 SFG(ITYP)=SFG(ITYP)+(1.D0-XT)*DSFG(ITYP) RDP=0.D0 ADP=0.D0 HDP=0.D0 CALL NORME(SFG,DL3) DL3=DL3/YOUN*10.D0 XLAM3=XLAM3+DL3 CALL CPHOIN(SFG,EPPLDP,YOUN,ANU) DO 130 ITYP=1,3 IDAM(ITYP)=0 EPDP(ITYP)=EPDP(ITYP)+EPPLDP(ITYP) DSIG(ITYP)=0.D0 130 SIG(ITYP)=0.D0 CALL CHREP(EPDP,EPPLDP,-ANG) IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER RETURN 131 ITER=ITER+1 IF(ITER.GT.200) THEN KERRE=460 RETURN ENDIF C C--------------------------------- C ON REALISE LE COUPLAGE C--------------------------------- C E=-(SFG(1)+X*DSFG(1))/SG(1) F=SSI(1)/SG(1) DO 140 ITYP=1,3 DSGG(ITYP)=X*DSFG(ITYP) A(ITYP)=SFG(ITYP)+DSGG(ITYP)+E*SG(ITYP) 140 B(ITYP)=-SSI(ITYP)+F*SG(ITYP) CALL DLAMDP(A,B,DL3,RDP,ADP,HDP,ITEST) IF(ITEST.EQ.1) THEN X=X/2.D0 GOTO 131 ENDIF DL1=-E-F*DL3 DO 145 ITYP=1,3 145 SGG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) C C---------------------------------------------- C SI UN DES DLAMDA < 0 FAUX COUPLAGE C---------------------------------------------- C IF(DL1.GT.0.D0) IDAM(1)=0 IF(DL3.LT.0.D0) IDAM(3)=0 IF(IDAM(1).EQ.0.OR.IDAM(3).EQ.0) THEN DL1=0.D0 DL3=0.D0 X=0.D0 GOTO 3000 ENDIF C C----------------------------------------------------- C ON VERIFIE QUE LA NORMALE SUIVANT LE DRUCKER C PRAGER TOURNE PEU C----------------------------------------------------- C CALL SIJ(SGG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC1) CALL SCAL(EPC,EPC1,VAL) CALL NORME(EPC,VA1) CALL NORME(EPC1,VA2) CO=VAL/VA1/VA2 IF(CO.LT.0.99) THEN X=X/2.D0 GOTO 131 ENDIF IF((DL1+XLAM1).GE.0.D0) GOTO 200 C C------------------------------------------------ C CAS OU LA FISSURE EST COMPLETEMENT FERMEE C------------------------------------------------ C IDAM(1)=0 DL1=-XLAM1 IF(DSFG(1).EQ.0.D0) GOTO 160 E=(DL1*SG(1)-SFG(1))/DSFG(1) F=SSI(1)/DSFG(1) DO 150 ITYP=1,3 A(ITYP)=SFG(ITYP)-DL1*SG(ITYP)+E*DSFG(ITYP) 150 B(ITYP)=F*DSFG(ITYP)-SSI(ITYP) CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST) DO 151 ITYP=1,3 SGG(ITYP)=A(ITYP)+DL31*B(ITYP) 151 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP) CALL SIJ(SGG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC1) CALL SCAL(EPC,EPC1,VAL) CALL NORME(EPC,VA1) CALL NORME(EPC1,VA2) CO1=VAL/VA1/VA2 CALL SIJ(DSGG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC1) CALL SCAL(EPC,EPC1,VAL) CALL NORME(EPC,VA1) CALL NORME(EPC1,VA2) CO2=VAL/VA1/VA2 XIN1=E+F*DL31 XIN2=E+F*DL32 IF(XIN1.GT.-1.D-10.AND.XIN1.LE.X.AND.CO1.GT.0.9) THEN DL3=DL31 X=XIN1 GOTO 200 ENDIF IF(XIN2.GT.-1.D-10.AND.XIN2.LE.X.AND.CO2.GT.0.9) THEN DL3=DL32 X=XIN2 GOTO 200 ENDIF IF(IIMPI.EQ.9) WRITE(IOIMP,10101) 10101 FORMAT(1X,'ERREUR DANS C1DP FERM',/) DL3=DL31 X=XIN1 GOTO 200 160 DL3=(SFG(1)-DL1*SG(1))/SSI(1) RD=RDP-HDP*DL3 IF(RD.LT.0.D0) THEN RDP=0.D0 GOTO 110 ENDIF DO 170 ITYP=1,3 A(ITYP)=SFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) 170 B(ITYP)=DSFG(ITYP) CALL XDP(A,B,RD,ADP,X,ITEST) 200 CONTINUE DO 205 ITYP=1,3 DSGG(ITYP)=X*DSFG(ITYP) 205 SGG(ITYP)=SFG(ITYP)+DSGG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) C C----------------------------------------------- C ON REGARDE SI ON ENDOMMAGE PAS LES AUTRES C CRITERES C----------------------------------------------- C------------------------------------- C ON ENDOMMAGE LA TRACTION 2 C------------------------------------- C CALL CTRAF(SGG(2),RT2,VCTR) IF(VCTR.LE.0.D0) GOTO 1000 IDAM(1)=-1 IDAM(2)=1 IDAM(3)=1 AT(1)=-DSFG(1) AT(2)=Y AT(4)=Y*ANU AT(3)=-DSFG(2) BT(1)=SFG(1) BT(2)=SFG(2)-RT2 CALL SYLIN2(AT,BT,X1,DL11) BT(1)=-SSI(1) BT(2)=-SSI(2) CALL SYLIN2(AT,BT,X2,DL12) DO 230 ITYP=1,3 A(ITYP)=SFG(ITYP)-DL11*SG(ITYP)+X1*DSFG(ITYP) 230 B(ITYP)=-DL12*SG(ITYP)+X2*DSFG(ITYP)-SSI(ITYP) CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST) DO 231 ITYP=1,3 SGG(ITYP)=A(ITYP)+DL31*B(ITYP) 231 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP) CALL SIJ(SGG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC1) CALL SCAL(EPC,EPC1,VAL) CALL NORME(EPC,VA1) CALL NORME(EPC1,VA2) CO1=VAL/VA1/VA2 CALL SIJ(DSGG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC1) CALL SCAL(EPC,EPC1,VAL) CALL NORME(EPC,VA1) CALL NORME(EPC1,VA2) CO2=VAL/VA1/VA2 DLIN1=DL11+DL31*DL12 DLIN2=DL11+DL32*DL12 XIN1=X1+DL31*X2 XIN2=X1+DL32*X2 IF(XIN1.GT.-1.D-10.AND.DLIN1.LT.1.E-10.AND.XIN1.LE.X 1 .AND.CO1.GT.0.9) THEN DL3=DL31 X=XIN1 DL1=DLIN1 GOTO 3000 ENDIF IF(XIN2.GT.-1.D-10.AND.DLIN2.LT.1.E-10.AND.XIN2.LE.X 1 .AND.CO2.GT.0.9) THEN DL3=DL32 X=XIN2 DL1=DLIN2 GOTO 3000 ENDIF IF(IIMPI.EQ.9) WRITE(IOIMP,20202) 20202 FORMAT(1X,'ERREUR DANS C1DP TRAC ',/) DL3=DL31 DL1=DLIN1 X=XIN1 GOTO 3000 C C------------------------------------- C ON ENDOMMAGE LA COMPRESSION 2 C------------------------------------- C 1000 CALL CCOAF(SGG(2),XLAM2,VCCO2) IF(VCCO2.LE.0.D0) GOTO 2000 IDAM(1)=-1 IDAM(3)=1 IDAM(2)=-1 AT(1)=-DSFG(1) AT(2)=Y AT(4)=Y*ANU AT(3)=-DSFG(2) BT(1)=SFG(1) BT(2)=SFG(2) CALL SYLIN2(AT,BT,X1,DL11) BT(1)=-SSI(1) BT(2)=-SSI(2) CALL SYLIN2(AT,BT,X2,DL12) DO 1010 ITYP=1,3 A(ITYP)=SFG(ITYP)-DL11*SG(ITYP)+X1*DSFG(ITYP) 1010 B(ITYP)=-DL12*SG(ITYP)+X2*DSFG(ITYP)-SSI(ITYP) CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST) DO 1011 ITYP=1,3 SGG(ITYP)=A(ITYP)+DL31*B(ITYP) 1011 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP) CALL SIJ(SGG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC1) CALL SCAL(EPC,EPC1,VAL) CALL NORME(EPC,VA1) CALL NORME(EPC1,VA2) CO1=VAL/VA1/VA2 CALL SIJ(DSGG,SI,SIETOI) CALL EPSDP(SI,SIETOI,ADP,EPC1) CALL SCAL(EPC,EPC1,VAL) CALL NORME(EPC,VA1) CALL NORME(EPC1,VA2) CO2=VAL/VA1/VA2 DLIN1=DL11+DL31*DL12 DLIN2=DL11+DL32*DL12 XIN1=X1+DL31*X2 XIN2=X1+DL32*X2 IF(XIN1.GT.-1.D-10.AND.DLIN1.LT.1.E-10.AND.XIN1.LE.X 1 .AND.CO1.GT.0.9) THEN DL3=DL31 X=XIN1 DL1=DLIN1 GOTO 3000 ENDIF IF(XIN2.GT.-1.D-10.AND.DLIN2.LT.1.E-10.AND.XIN2.LE.X 1 .AND.CO2.GT.0.9) THEN DL3=DL32 X=XIN2 DL1=DLIN2 GOTO 3000 ENDIF IF(IIMPI.EQ.9) WRITE(IOIMP,30303) 30303 FORMAT(1X,'ERREUR DANS C1DP COMP ',/) DL3=DL31 DL1=DLIN1 X=XIN1 GOTO 3000 2000 CONTINUE C C------------------------------------------------ C CAS OU AUCUN AUTRE CRITERE EST ENDOMMAGE C------------------------------------------------ C----------------------------------- C REMISE A JOUR DES VARIABLES C PUIS NOUVEAU PAS SI TOUT C L INCREMENT N EST PAS ECOULE C----------------------------------- C XT=XT+X RDP=RDP-DL3*HDP RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) IF(SIREF.GT.RMAX) THEN RDP=0.D0 GOTO 110 ENDIF DO 2010 ITYP=1,3 SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) 2010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP) XLAM1=XLAM1+DL1 XLAM3=XLAM3+DL3 IF(XLAM1.LT.1.E-10) XLAM1=0.D0 IF(XLAM1.LT.1.E-10) IDAM(1)=0 C=1.D0-1.D-10 IF(XT.LT.C) GOTO 2020 DO 2015 ITYP=1,3 IDAM(ITYP)=0 2015 DSIG(ITYP)=0.D0 CALL CHREP(SFG,SIG,-ANG) CALL CHREP(EPDP,EPPLDP,-ANG) IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER RETURN 2020 IF(IDAM(1).EQ.-1) GOTO 100 XLAM1=0.D0 GAMDP=10.D0 GAMTR1=10.D0 GAMTR2=10.D0 GAMCO2=10.D0 DO 2030 ITYP=1,3 DSFG(ITYP)=DSFG(ITYP)*(1.D0-XT) 2030 TENS(ITYP)=SFG(ITYP)+DSFG(ITYP) CALL CTRAF(TENS(1),RT1,VCTR1) CALL CTRAF(TENS(2),RT2,VCTR2) CALL CCOAF(TENS(2),XLAM2,VCCO2) CALL CDP(TENS,ADP,RDP,VCDP) IF(VCTR1.GT.0.D0) CALL GAMTAF(SFG(1),DSFG(1),RT1,GAMTR1) IF(VCTR2.GT.0.D0) CALL GAMTAF(SFG(2),DSFG(2),RT2,GAMTR2) IF(VCCO2.GT.0.D0) CALL GAMCAF(SFG(2),DSFG(2),GAMCO2) IF(VCDP.GT.0.D0) CALL GDP(SFG,DSFG,RDP,ADP,GAMDP) IDAM(1)=0 IDAM(2)=0 IDAM(3)=0 GAM=MIN(GAMTR1,GAMTR2,GAMCO2,GAMDP) IF(GAM.GE.1.D0) THEN DO 2040 ITYP=1,3 SFG(ITYP)=SFG(ITYP)+DSFG(ITYP) 2040 DSIG(ITYP)=0.D0 CALL CHREP(SFG,SIG,-ANG) RETURN ENDIF IF(ABS(GAM-GAMTR1).LE.1.E-10) IDAM(1)=1 IF(ABS(GAM-GAMTR2).LE.1.E-10) IDAM(2)=1 IF(ABS(GAM-GAMCO2).LE.1.E-10) IDAM(2)=-1 IF(ABS(GAM-GAMDP).LE.1.E-10) IDAM(3)=1 DO 2050 ITYP=1,3 SFG(ITYP)=SFG(ITYP)-GAM*DSFG(ITYP) 2050 DSFG(ITYP)=(1.D0-GAM)*DSFG(ITYP) CALL CHREP(SFG,SIG,-ANG) CALL CHREP(DSFG,DSIG,-ANG) CALL CHREP(EPDP,EPPLDP,-ANG) RETURN 3000 RDP=RDP-DL3*HDP RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP))) IF(SIREF.GT.RMAX) THEN RDP=0.D0 GOTO 110 ENDIF XT=XT+X DO 3010 ITYP=1,3 SFG(ITYP)=SFG(ITYP)+X*DSFG(ITYP)-DL1*SG(ITYP)-DL3*SSI(ITYP) DSFG(ITYP)=(1.D0-XT)*DSFG(ITYP) 3010 EPDP(ITYP)=DL3*EPC(ITYP)+EPDP(ITYP) CALL CHREP(SFG,SIG,-ANG) CALL CHREP(DSFG,DSIG,-ANG) CALL CHREP(EPDP,EPPLDP,-ANG) XLAM1=XLAM1+DL1 XLAM3=XLAM3+DL3 IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER RETURN END