C BONE SOURCE CHAT 05/01/12 21:42:05 5004 SUBROUTINE BONE(SIGR,SIGF,DSTRN,IPLA,IFISU,SIG1,SIG2 A ,NSTRS,D,D1,IFOUB,SIGP,EPST,SIR,SIRL,ENDO B ,ITHHER,T0,TF,BETJEF,VISCO,BETFLU B ,NECH0,NECH1,NECH2,NECH3) C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION SIGR(4),DSIGT(4),SIGF(4),S1(4),DSTRN(4),SIGT(4) DIMENSION EPSR(9),S2(4),V2(4),V1(4) DIMENSION S(6),SI1(6),EPST(4),SIGRV(4) DIMENSION ST0(6),V02(4),SI0(6),V01(4),SIGP(4) DIMENSION SIGIST(4),EXU(9),CODU(9,9),COD(8),TZU(19) DIMENSION STRNI(9),SOM(9),DEFMAX(9),DDSTRN(9) DIMENSION H(4,4),DSTFL(4),S3(4) DIMENSION DEP(4,4),D1(4,4),D(4,4),SIRL(8,4),SIR(9,4) DIMENSION DEPE(4,4),DH(4,4),HI(4,4),CODL(8,8) * SEGMENT BETJEF REAL*8 AA,BETA,RB,ALPHA,EX,XNU,GFC,GFT,CAR,ETA,TDEF, & TCON,DPSTF1,DPSTF2,TETA,PDT,TP00 INTEGER ICT,ICC,IMOD,IVIS,ITER, & ISIM,IBB,IGAU,IZON ENDSEGMENT SEGMENT VISCO REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV ENDSEGMENT SEGMENT BETFLU REAL*8 DATCOU,DATCUR,DATSEC,E28,PGTZO,PGDUR,TAU1,TAU2, & TP0,TZER INTEGER ITYPE,IMD,NBRC,NCOE,NTZERO,NTPS,IFOR ENDSEGMENT SEGMENT NECH0 REAL*8 DT,DC,ALFG,S0,ENDO ENDSEGMENT SEGMENT NECH1 REAL*8 DLMT ENDSEGMENT SEGMENT NECH2 REAL*8 ATR,GTR,ALPH0 ENDSEGMENT SEGMENT NECH3 REAL*8 RBT,ALFAT,YOUNT,GFCT,GFTT,ALPH ENDSEGMENT C TAU1=0.05 NTPS=18 NTZERO=18 NCOE=8 MC=NBRC+1 NC=NCOE+1 C C-------------------------------------------------------------------------- C ************************************************************ C * APPLICATION DES CRITERES DE PLASTICITE ET FISSURATION * C * CRITERE BETON * C * * C * IFISU=INDICE DE FISSURATION * C * =0 PAS DE FISSURE * C * =1 UNE FISSURE * C * =2 RUINE DANS DIRECTION DE TRACTION * C * * C * IPLA =INDICE DE PLASTICITE * C * =0 ELASTIQUE * C * =1 ECROUISSAGE POSITIF * C * =2 RUPTURE PAR COMPRESSION DANS 1 DIRECTION * C * * C * * C * * C ************************************************************ C--------------------------------------------------------------------- AA = 1.D0/3.D0 CALL ZERO(DSIGT,4,1) CALL ZERO(SIGT,4,1) CALL ZERO(SIGIST,4,1) CALL ZERO(SIGRV,4,1) CALL ZERO(CODU,9,9) CALL ZERO(CODL,8,8) CALL ZERO(COD,8,1) CALL ZERO(S1,4,1) CALL ZERO(S2,4,1) CALL ZERO(V1,4,1) CALL ZERO(V2,4,1) CALL ZERO(D,4,4) CALL ZERO(DEP,4,4) CALL ZERO(H,4,4) CALL ZERO(HI,4,4) DAM = 0.D0 C C ****************** INITIALISATION POUR VISCOPLASTIQUE *************** C IF (IVIS.EQ.1.OR.IVIS.EQ.4) THEN DO 5 I=1,NSTRS SIGRV(I)= SIGR(I) 5 CONTINUE ENDIF IF ((IVIS.EQ.1.AND.(IPLA.GT.0.OR.IFISU.GT.0)).OR. & (IVIS.EQ.4.AND.(IPLA.GT.0.OR.IFISU.GT.0))) THEN C DO 6 I=1,NSTRS SIGR(I)= SIGP(I) 6 CONTINUE ENDIF DO 7 I=1,NSTRS SIGR(I)= SIGR(I)/(1.D0-ENDO) 7 CONTINUE C C ****************** INITIALISATION POUR VISCOELASTIQUE *************** C IF (IVIS.EQ.2) THEN SEL1=0.D0 SEL2=0.D0 EX=0.D0 TPS1 = TP0 TPS2 = TP0+PDT ENDIF C CALL INFICH(CODU,CODL,COD,BETJEF,BETFLU) CALL MODBET(TPS1,TPS2,SEL1,SEL2,EXU,EXUL,EX,CODU,CODL & ,COD,BETJEF,BETFLU) C WRITE(*,*)'EX=',EX,'a T=',TPS2 IF (IFOR.EQ.1) THEN CALL CALIS1(SIGIST,NSTRS,DSTRN,IFOUB,SIR,CODU,CODL, & COD,BETJEF,BETFLU) ENDIF IF (IFOR.EQ.2) THEN CALL CALIS2(SIGIST,NSTRS,DSTRN,IFOUB,SIRL,CODU,CODL, & COD,BETJEF,BETFLU) ENDIF C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C INITIALISATION DE L'ENDOMMAGEMENT C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF(IVIS.EQ.3.OR.IVIS.EQ.4) THEN IF((IFISU.EQ.0.D0).AND.(IPLA.EQ.0.D0))THEN ENDO=0.D0 DO 300 I =1, NSTRS STRNI(I)=0.D0 300 CONTINUE ENDIF C DO 400 I=1, NSTRS SOM(I)=SOM(I)+DSTRN(I) DEFMAX(I)=MAX(DABS(SOM(I)),DABS(STRNI(I))) DDSTRN(I)=DABS(SOM(I))-DABS(DEFMAX(I)) 400 CONTINUE C DO 500 I=1, NSTRS STRNI(I)=DEFMAX(I) 500 CONTINUE C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C TEST SUR LE TYPE DE CHARGEMENT C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C DO 200 I=1,NSTRS IF (DDSTRN(I).LT.0.D0) THEN DSTRN(I)=(1.D0-ENDO)*DSTRN(I) ENDIF 200 CONTINUE C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C CALCUL DE LA MATRICE D'INTERACTION THERMO-MECA [H] C ET INVERSE DE [H] C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C CALL INTERF(ITHHER,T0,TF,NSTRS,H,BETJEF,NECH2) C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C CALCUL DE LA DEFORMATION D'INTERACTION C THERMO-MECA [DSTFL]AU PAS N C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C CALL TRANSI(ITHHER,T0,TF,NSTRS,HI,BETJEF,NECH2) C DO 250 I=1,NSTRS DSTFL(I)=0.D0 DO 250 J=1,NSTRS DSTFL(I)=DSTFL(I)+HI(I,J)*SIGF(J) 250 CONTINUE C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C CALCUL DE LA DEFORMATION DU TIRE ELASTIQUE C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C DO 260 I=1,NSTRS DSTRN(I)=DSTRN(I)-DSTFL(I) 260 CONTINUE C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C CALL INVMA2(H,NSTRS,ISING) ENDIF C C ************************** CAS CONTRAINTE PLANE ********************* C IF (IFOUB.EQ.-2) THEN IF (IVIS.EQ.3.OR.IVIS.EQ.4) THEN CALL CREMAT(DH,EX,XNU,3,IFOUB) CALL PRODMM(H,DH,DEP,NSTRS,NSTRS,NSTRS) ELSE CALL CREMAT(DEP,EX,XNU,3,IFOUB) ENDIF IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN CALL PROMA2(DEP,DSTRN,3,DSIGT) DO 2 I=1,NSTRS SIGT(I)=SIGR(I)+DSIGT(I)+SIGIST(I) 2 CONTINUE ENDIF IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN WRITE(*,*)'ATTENTION IMOD=2 ou 4 -> DEFORMATIONS PLANES' STOP ENDIF ENDIF C C ***** CAS DEFORMATION PLANE ou AXISYMETRIQUE *************** C IF (IFOUB.EQ.-1.OR.IFOUB.EQ.0) THEN IF (IVIS.EQ.3.OR.IVIS.EQ.4) THEN CALL CREMAT(DH,EX,XNU,4,IFOUB) CALL PRODMM(H,DH,DEP,NSTRS,NSTRS,NSTRS) ELSE CALL CREMAT(DEP,EX,XNU,4,IFOUB) ENDIF CALL PROMA2(DEP,DSTRN,4,DSIGT) DO 3 I=1,NSTRS SIGT(I)=SIGR(I)+DSIGT(I)+SIGIST(I) 3 CONTINUE IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN WRITE(*,*)'ATTENTION IMOD=1 ou 3 -> CONTRAINTES PLANES' STOP ENDIF ENDIF C C ************************ Normalisation ********************* C DO 1 I=1,NSTRS IF(DABS(SIGR(I)).LT.1.D-6) SIGR(I)=0.D0 IF(DABS(SIGT(I)).LT.1.D-6) SIGT(I)=0.D0 1 CONTINUE C--------------------------------------------------------------------- DO 11 I=1,NSTRS S1(I)=SIGT(I) 11 CONTINUE CALL PRINC(S1,V1,NSTRS) C C **************** TEST DU TYPE DE SOLLICITATION ********************* C IZON=1 IF(V1(1).LT.0.D0.AND.V1(2).LT.0.D0) IZON=0 IF(V1(1).GE.0.D0.AND.V1(2).GE.0.D0) IZON=2 C--------------------------------------------------------------------- C IF (IFISU.NE.0) GOTO 15 IF (IZON.EQ.0) GOTO 30 IF (IZON.EQ.1.OR.IZON.EQ.2) GOTO 15 C C ******************** COMPORTEMENT DE TRACTION ********************** C 15 CONTINUE C ##################################### C * POINT DEJA FISSURE * C * COMPORTEMENT PLASTIQUE * C ##################################### C IF(IFISU.NE.0) THEN IF(IZON.EQ.2) THEN CALL BEHAV1(SIGR,DSTRN,DPSTF1,DPSTF2,IPLA, A SIG1,SIG2,IFISU,S1,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, A BETJEF,VISCO,NECH0,NECH1) GOTO 40 ENDIF IF(IZON.EQ.1) THEN CALL BEHAV3(SIGR,DSTRN,DPSTF1,DPSTF2,IPLA, A SIG1,SIG2,IFISU,S1,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, A BETJEF,VISCO,NECH0,NECH1) GOTO 40 ENDIF IF(IZON.EQ.0) THEN WRITE(*,*)'Dans elasbet un point deja fissure est teste en *bicompression' WRITE(*,*)'Element num:',IBB ,'au point de Gauss num:',IGAU STOP ENDIF ENDIF C--------------------------------------------------------------------- C ##################################### C * POINT FISSURE 1ER FOIS * C * COMPORTEMENT PLASTIQUE * C ##################################### C Ft=ALPHA*RB IF(V1(1).GT.Ft) THEN IF(IZON.EQ.2) THEN CALL BEHAV1(SIGR,DSTRN,DPSTF1,DPSTF2,IPLA, A SIG1,SIG2,IFISU,S1,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, A BETJEF,VISCO,NECH0,NECH1) GOTO 40 ENDIF IF(IZON.EQ.1) THEN CALL BEHAV3(SIGR,DSTRN,DPSTF1,DPSTF2,IPLA, A SIG1,SIG2,IFISU,S1,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, A BETJEF,VISCO,NECH0,NECH1) GOTO 40 ENDIF ENDIF IF(IZON.EQ.1.AND.IPLA.EQ.0) THEN IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN CALL CRIVON(S1,SEQ,NSTRS,BETJEF) ENDIF IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN CALL CRIDRU(S1,SEQ,NSTRS,BETJEF) ENDIF FCRI=SEQ-RB*AA IF(FCRI.GT.0.D0) THEN CALL BEHAV3(SIGR,DSTRN,DPSTF1,DPSTF2,IPLA, A SIG1,SIG2,IFISU,S1,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, A BETJEF,VISCO,NECH0,NECH1) GOTO 40 ENDIF ENDIF IF(IZON.EQ.1.AND.IPLA.GE.1) THEN IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN CALL CRIVON(S1,SEQ,NSTRS,BETJEF) ENDIF IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN CALL CRIDRU(S1,SEQ,NSTRS,BETJEF) ENDIF FCRI=SEQ-SIG2 IF(FCRI.GT.0.D0) THEN CALL BEHAV3(SIGR,DSTRN,DPSTF1,DPSTF2,IPLA, A SIG1,SIG2,IFISU,S1,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, A BETJEF,VISCO,NECH0,NECH1) GOTO 40 ENDIF ENDIF GOTO 40 C C *************** COMPORTEMENT DE BICOMPRESSION ********************** C 30 CONTINUE C--------------------------------------------------------------------- IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN CALL CRIVON(S1,SEQ,NSTRS,BETJEF) ENDIF IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN CALL CRIDRU(S1,SEQ,NSTRS,BETJEF) ENDIF FCRI=SEQ-RB*AA IF (FCRI.GT.0.D0.OR.IPLA.NE.0) THEN CALL BEHAV2(SIGR,DSTRN,DPSTF1,DPSTF2,IPLA, A SIG1,SIG2,IFISU,S1,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, A BETJEF,VISCO,NECH0,NECH1) GOTO 40 ENDIF C--------------------------------------------------------------------- 40 CONTINUE C--------------------------------------------------------------------- C IF (IVIS.EQ.3.OR.IVIS.EQ.4) THEN IF(IVIS.EQ.3) THEN CALL REELLE(S1,DEP,DPSTF1,DPSTF2,DEPE,S2,DAM, & NSTRS,IFISU,IPLA,BETJEF,NECH0,NECH1) ENDO = DAM DO 98 I=1,NSTRS S1(I)= S2(I) 98 CONTINUE ENDIF IF(IVIS.EQ.4.AND.(IPLA.GT.0.OR.IFISU.GT.0)) THEN CALL REELLE(S1,DEP,DPSTF1,DPSTF2,DEPE,SIGP,DAM, & NSTRS,IFISU,IPLA,BETJEF,NECH0,NECH1) ENDO = DAM CALL VISPL1(SIGP,DSIGT,NSTRS,S2,SIGRV,DSTRN, & DPSTF1,DPSTF2,BETJEF,VISCO,NECH0,NECH1) DO 99 I=1,NSTRS S1(I)= S2(I) 99 CONTINUE ENDIF ENDIF C C MODIF DEP ET CONTRAINTES C DO 121 I=1,NSTRS DO 121 J=1,NSTRS D(I,J)=DEP(I,J) 121 CONTINUE C--------------------------------------------------------------------- 50 CONTINUE C--------------------------------------------------------------------- DO 42 I=1,NSTRS C IF(DABS(S1(I)).LT.1.D-6) S1(I)=0.D0 SIGF(I)= S1(I) C WRITE(*,*)'SIGF(',I,')=',SIGF(I) 42 CONTINUE C C 70 FORMAT (4(1X,E12.5)) RETURN END