1 : C FIBMAZ SOURCE CHAT 12/04/06 21:15:11 7348 2 : SUBROUTINE FIBMAZ (XMAT,DEPS,SIG0,VAR0,SIGF,VARF) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : DIMENSION XMAT(*),DEPS(3),SIG0(3),VAR0(5),SIGF(3),VARF(5) 6 : C--------------------------------------------------------------- 7 : C MODELE DE MAZARS (MAZARS_FIB) SPECIALISE POUR LE MODELE A FIBRE 8 : C Didier Combescure et Pierre Pegon Nov. 93 9 : C--------------------------------------------------------------- 10 : C 11 : C 12 : C variables en entree 13 : C 14 : C XMAT : Caracteristique du materiaux 15 : C 16 : C EPS : Deformations initiales 17 : C 18 : C DEPS : Increment de deformations 19 : C 20 : C SIG0 : Contraintes initiales 21 : C 22 : C VAR0 : Variables internes initiales 23 : C 24 : C 25 : C variables en sortie 26 : C 27 : C SIGF : Contraintes finales 28 : C 29 : C VARF : Variables internes finales 30 : C 31 : C declaration des variables 32 : C 33 : C 34 : REAL*8 EPS33(3,3),EPSIPP(3),EPSILC(3),VALP33(3,3) 35 : REAL*8 SIGP(3),SIGPC(3) 36 : PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12) 37 : 38 : C 39 : C 40 : C recuperation des variables initiales dans les tableaux 41 : C 42 : C 43 : YOUN = XMAT(1) 44 : XNU = XMAT(2) 45 : EPSD0= XMAT(5) 46 : ACOM = XMAT(6) 47 : BCOM = XMAT(7) 48 : ATRA = XMAT(8) 49 : BTRA = XMAT(9) 50 : BETA = XMAT(10) 51 : DINI = VAR0(2) 52 : EPSX = VAR0(3) 53 : EPSY = VAR0(4) 54 : EPSZ = VAR0(5) 55 : 56 : C 57 : C calcul des seuils 58 : C 59 : C 60 : C calcul de la deformation totale 61 : C 62 : DO 100 ISTRS=1,3 63 : VARF(ISTRS+2)=VAR0(ISTRS+2)+DEPS(ISTRS) 64 : 100 CONTINUE 65 : EPS33(1,1) = VARF(3) 66 : EPS33(1,2) = (VARF(4)/2.D0) 67 : EPS33(2,1) = EPS33(1,2) 68 : EPS33(1,3) = (VARF(5)/2.D0) 69 : EPS33(3,1) = EPS33(1,3) 70 : EPS33(2,2) = ((-1.D0*XNU)*VARF(3)) 71 : EPS33(3,3) = EPS33(2,2) 72 : EPS33(3,2) = XZERO 73 : EPS33(2,3) = XZERO 74 : 75 : C 76 : C calcul des deformations principales 77 : C 78 : C on diagonalise 79 : C 79 : CALL JACOB3 (EPS33,3,EPSIPP,VALP33) 80 : 81 : C 82 : C on calcule les contraintes ppales 83 : C 84 : C 85 : SIGP(1)=((UN-XNU)*EPSIPP(1)+ 86 : 1 (XNU*EPSIPP(2))+ 87 : 2 (XNU*EPSIPP(3))) 88 : SIGP(2)=((UN-XNU)*EPSIPP(2)+ 89 : 1 (XNU*EPSIPP(1))+ 90 : 2 (XNU*EPSIPP(3))) 91 : SIGP(3)=((UN-XNU)*EPSIPP(3)+ 92 : 1 (XNU*EPSIPP(2))+ 93 : 2 (XNU*EPSIPP(1))) 94 : DO 200 ISTRS=1,3 95 : SIGP(ISTRS)=(SIGP(ISTRS)*YOUN/((UN+XNU)*(UN-2.D0*XNU))) 96 : 200 CONTINUE 97 : 98 : C 99 : C on calcule le epsilontild 100 : C 101 : EPSTIL=MAX( XZERO , EPSIPP(1) )**2 + 102 : 1 MAX( XZERO , EPSIPP(2) )**2 + 103 : 2 MAX( XZERO , EPSIPP(3) )**2 104 : EPSTIL=SQRT (EPSTIL) 105 : C 106 : C on est toujours dans le cas istep=0 107 : C 108 : EPSTIM=EPSTIL 109 : VARF(1)=EPSTIL 110 : C 111 : C on calcule l'endommagement et les contraintes 112 : C 113 : C 114 : C on calcule ALFAT ALFAC DT et DC puis D 115 : C dans le cas ou le seuil initial est depasse 116 : C 117 : C on calcule le signe des contraintes elastiques 118 : C 119 : DO 300 ISTRS=1,3 120 : SIGPC(ISTRS)=MIN(XZERO,SIGP(ISTRS)) 121 : 300 CONTINUE 122 : TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3) 123 : C 124 : C ainsi que GAMMA et EPSTIM si BETA est different de zero 125 : C 126 : IF (BETA.GT.1.D-5) THEN 127 : IF (ABS(TRSIGC).GT.XPETIT) THEN 128 : GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+ 129 : 1 SIGPC(3)*SIGPC(3) 130 : GAMMA= SQRT(GAMMA)/ABS(TRSIGC) 131 : EPSTIM=EPSTIM*GAMMA 132 : END IF 133 : END IF 134 : C 135 : IF ( EPSTIM .GT. EPSD0) THEN 136 : C 137 : C calcul de l'endommagement 138 : C 139 : C on calcule les deformations dues aux contraintes negatives 140 : C 141 : DO 400 ISTRS=1,3 142 : EPSILC(ISTRS)=(SIGPC(ISTRS)*(UN+XNU)-TRSIGC*XNU)/YOUN 143 : 400 CONTINUE 144 : C 145 : C on en deduit ALFAT et ALFAC 146 : C 147 : ALFAC = MAX(XZERO,EPSIPP(1))*EPSILC(1) + 148 : 1 MAX(XZERO,EPSIPP(2))*EPSILC(2) + 149 : 2 MAX(XZERO,EPSIPP(3))*EPSILC(3) 150 : ALFAC = ALFAC/(EPSTIL*EPSTIL) 151 : ALFAT = UN - ALFAC 152 : C 153 : C amelioration de la reponse en cisaillement pour beta > 1. 154 : C 155 : IF (BETA .GT. UN) THEN 156 : IF ( ALFAT .GT. XPETIT ) THEN 157 : ALFAT=ALFAT**BETA 158 : END IF 159 : IF ( ALFAC .GT. XPETIT ) THEN 160 : ALFAC=ALFAC**BETA 161 : END IF 162 : END IF 163 : 164 : C 165 : C on calcule DT et DC puis D 166 : C dans le cas ou le seuil initial est depasse 167 : C 168 : DT=UN - EPSD0*(UN-ATRA)/EPSTIM - 169 : 1 ATRA/EXP(BTRA*(EPSTIM-EPSD0)) 170 : DC=UN - EPSD0*(UN-ACOM)/EPSTIM - 171 : 1 ACOM/EXP(BCOM*(EPSTIM-EPSD0)) 172 : D = ALFAT*DT + ALFAC*DC 173 : C 174 : C on borne la valeur de D a 0.99999 175 : C 176 : D=MIN ( D , UN-1.D-05 ) 177 : ELSE 178 : D=XZERO 179 : END IF 180 : C 181 : C on teste la croissance de D 182 : C 183 : D=MAX ( D , DINI ) 184 : C 185 : C on le stocke dans les variables internes finales 186 : C 187 : VARF(2)= D 188 : C 189 : C on calcule les contraintes finales 190 : C 191 : SIGF(1) = YOUN*VARF(3) 192 : SIGF(2) = YOUN*VARF(4)/2.D0/(UN+XNU) 193 : SIGF(3) = YOUN*VARF(5)/2.D0/(UN+XNU) 194 : DO 500 ISTRS=1,3 195 : SIGF(ISTRS)=SIGF(ISTRS)*(UN-D) 196 : 500 CONTINUE 197 : 198 : RETURN 199 : END 200 : 201 : 202 :
© Cast3M 2003 - All rights reserved.
Disclaimer