1 : C MICRO1 SOURCE CHAT 05/01/13 01:44:59 5004 2 : SUBROUTINE MICRO1(WRK0,NMATT,WRK1,NSTRS,DSIGT, 3 : 1 NVARI,IDECAL,KERRE,LIMPR) 4 : IMPLICIT INTEGER(I-N) 5 : -INC CCOPTIO 6 : INTEGER NSTRS,NVARI,KERRE,ISTRS,NMATT,ITER,NSIG,IDECAL 7 : LOGICAL CONV,LCP,LIMPR 8 : REAL*8 SEQ,TRDFDS,YOUNG,XNU,DEUXMU,LAMB,R0,DR,DR0,RINF,RI 9 : REAL*8 FPC,MP,EPSE,R,TEST,TRSIG,F0,DELTL,DELTLT 10 : REAL*8 DEPS3,EPS3,DFDELTL,DFDEPS3,DS3DELTL,DS3DEPS3 11 : REAL*8 A,B,C,A1,B1,C1,A2,B2,C2 12 : CHARACTER*8 CMATE 13 : REAL*8 DSIGT(6) 14 : * 15 : SEGMENT WRK0 16 : REAL*8 XMAT(NMATT) 17 : ENDSEGMENT 18 : * 19 : SEGMENT WRK1 20 : REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) 21 : REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) 22 : REAL*8 DEFP(NSTRS),XCAR(ICARA) 23 : ENDSEGMENT 24 : * 25 : SEGMENT WRKK1 26 : REAL*8 SIGE0(NSTRS),DFDS(NSTRS) 27 : END SEGMENT 28 : * 29 : SEGINI WRKK1 30 : YOUNG=XMAT(1) 31 : XNU=XMAT(2) 32 : FPC=XMAT(7) 33 : MP=XMAT(8) 34 : EPSE=VAR0(1) 35 : RI=VAR0(2) 36 : R0=0.4D0*FPC 37 : RINF=0.6d0*FPC 38 : LCP=IFOUR.EQ.-2 39 : DEUXMU=YOUNG/(1.D0+XNU) 40 : LAMB=XNU*DEUXMU/(1.D0-2.D0*XNU) 41 : TEST=1.D-6*R0 42 : DR0=0.D0 43 : KERRE=0 44 : C 45 : C calcul de la matrice elastique 46 : C 47 : CMATE = 'ISOTROPE' 48 : KCAS=1 48 : CALL DOHMAS(XMAT,CMATE,IFOUR,NSTRS,KCAS,DDHOOK,IRTD) 49 : IF ( IRTD .EQ. 0) THEN 50 : KERRE=123 51 : return 52 : end if 53 : C 54 : C calcul de l'increment de contrainte 55 : C 55 : CALL MATVE1 (DDHOOK,DEPST,NSTRS,NSTRS,DSIGT,1) 56 : * print*,'IFOUR=',IFOUR 57 : * print*,'DDHOOK' 58 : * do i=1,3 59 : * print*,(ddhook(i,j),j=1,3) 60 : * end do 61 : * do i=1,3 62 : * print*,'DEPST(',I,')=',DEPST(i) 63 : * print*,'DSIGT(',I,')=',DSIGT(I) 64 : * end do 65 : * 66 : * calcul des parametres du critere de NADAI 67 : * 68 : C- 69 : C- A,B,C : PARAMETRES DU CRITERE DE PLASTICITE (NADAI) 70 : C- ILS SONT DIFFERENTS SUIVANT LE CHARGEMENT 71 : C- 72 : 73 : CONV=.FALSE. 73 : CALL ZDANUL (DEFP,NSTRS) 74 : ITER=0 75 : EPS3=0.d0 76 : DELTLT=0.D0 77 : R=RI 78 : A1=0.4d0 79 : B1=(1.d0 - A1) / (1.D0 + A1) 80 : C1=(2.D0 * A1) / (1.D0 + A1) 81 : A2=1.16d0 82 : B2=(A2 - 1.D0) / (2.D0 * A2 - 1.D0) 83 : C2=A2 / (2.D0 * A2 - 1.D0) 84 : * 85 : * NSIG est la dimension du tenseur des contraintes: 86 : * NSIG=3 sauf 2 en contraintes planes 87 : * 88 : IF (LCP)THEN 89 : NSIG=2 90 : ELSE 91 : NSIG=3 92 : END IF 93 : * 94 : * on recupere les contraintes anelastiques due a l'endommagement 95 : * 96 : DO ISTRS=1,NSTRS 97 : SIGE0(ISTRS)=VAR0(ISTRS+IDECAL) 98 : END DO 99 : * 100 : * prediction elastique de la contrainte effective 101 : * contrainte + contrainte anelastique + increment de contrainte elastique 102 : * 103 : * l'increment de contrainte est decoupe: 104 : * 1) point de contact avec la surface seuil ??? 105 : * 2) decoupage du reste de facon a ce que le pas en contrainte 106 : * equivalente soit inferieur a disons pour le moment E*1.e-04 107 : * 108 : * en deux coups les gros 109 : * 110 : * 111 : * DS=0.d0 112 : * do istrs=1,nstrs 113 : * ds=ds + dsigt(istrs)**2 114 : * end do 115 : * nstep=int(1.e+4*sqrt(ds)/youn) 116 : DO ISTRS=1,NSTRS 117 : SIGE0(ISTRS)=SIG0(ISTRS)+SIGE0(ISTRS)+DSIGT(ISTRS) 118 : END DO 119 : if(limpr)then 120 : DO ISTRS=1,NSTRS 121 : PRINT*,'SIG0(',ISTRS,')=',SIGE0(ISTRS) 122 : PRINT*,'DSIGT(',ISTRS,')=',DSIGT(ISTRS) 123 : END DO 124 : endif 125 : 126 : DO WHILE (ITER .LE. 100 .AND. .NOT. CONV) 127 : ITER=ITER+1 128 : * 129 : * on calcule la contrainte equivalente de Mises (SEQ) 130 : * et la trace des contraintes (TRSIG) 131 : * 132 : * on utilise la fonction modulo: 133 : * MOD( 1,3)= 1 134 : * MOD( 2,3)= 2 135 : * MOD( 3,3)= 0 136 : * MOD( 4,3)= 1 137 : * MOD( 5,3)= 2 138 : * MOD( 6,3)= 0 139 : * 140 : SEQ=0.D0 141 : TRSIG=0.d0 142 : * ici NSIG change en 3 le 6/11/95 (clb) 143 : DO ISTRS=1,3 144 : TRSIG=TRSIG + SIGE0(ISTRS) 145 : SEQ=SEQ+SIGE0(ISTRS)*(SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1)) 146 : END DO 147 : DO ISTRS=4,NSTRS 148 : SEQ=SEQ + 3.D0*SIGE0(ISTRS)**2 149 : END DO 150 : if (limpr) then 151 : print*,'VMIS**2=',SEQ,'TRSIG=',TRSIG 152 : end if 153 : SEQ=SQRT(SEQ) 154 : * 155 : * en fonction de la valeur de la trace de sigma, on determine les 156 : * coef b et c a applique a la surface seuil 157 : * 158 : IF (MOD(ITER-1,3).EQ.0 )then 159 : IF ( TRSIG .GT. (R+R0)*(C2-C1)/(B2-B1))THEN 160 : B=B1 161 : C=C1 162 : ELSE 163 : B=B2 164 : C=C2 165 : END IF 166 : ENDIF 167 : * 168 : * fonction seuil 169 : * 170 : F0=(SEQ + B*TRSIG)/C - R0 -R 171 : if (limpr) then 172 : PRINT*,'B=',B,'C=',C,'F0=',F0 173 : PRINT*,'SIGMAEQ=',(SEQ + B*TRSIG)/C 174 : end if 175 : * 176 : * le critere de convergence est verifie si la valeur absolue 177 : * de la fonction seuil est plus petite que le critere ou si 178 : * la fonction seuil est negative pour la premiere iteration 179 : * 180 : CONV=(ABS(F0).LE.TEST .AND. 181 : 1 (.NOT. LCP .OR. ABS(SIGE0(3)) .LE. TEST ) ) 182 : 2 .OR. (F0 .LE. 0.D0 .AND. ITER .EQ. 1) 183 : * 184 : * dans le cas des contraintes planes on rajoute aussi 185 : * que la valeur absolue de sigma33 est inferieure au critere 186 : * 187 : IF (.NOT. CONV) THEN 188 : * 189 : * calcul de dF/dSigma 190 : * 191 : DO ISTRS=1,3 192 : DFDS(ISTRS)=(2.D0*SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1) 193 : & -SIGE0(MOD(ISTRS+1,3)+1) )/(2.D0*C*SEQ) +B/C 194 : END DO 195 : DO ISTRS=4,NSTRS 196 : DFDS(ISTRS)=3.D0*SIGE0(ISTRS)/(2.d0*C*SEQ) 197 : END DO 198 : TRDFDS=0.D0 199 : DO ISTRS=1,NSIG 200 : TRDFDS=TRDFDS+DFDS(ISTRS) 201 : END DO 202 : if (limpr) then 203 : * on verifie que les valeur ppales sont bonnes 204 : do istrs=1,nstrs 205 : print*,'dfds(',istrs,')=',dfds(istrs) 206 : end do 207 : print*,'trdfds=',trdfds 208 : end if 209 : * 210 : * calcul de dsigma/deltl=C * depsp 211 : * attention !! dans le cas des contraintes planes:epsp33 n'est 212 : * pas une inconnue directe mais c'est depse33 213 : * 214 : * 215 : * R=MP*p -----> DR=MP * Dp =MP*DELTL 216 : * R=Rinf(1-exp(-bR)) ---> dR=b(Rinf-R)deltl 217 : * 218 : * DR = (RINF - R)*MP 219 : * DR=MP 220 : * R=MP*(p)**0.5 ---> DR=0.5 * MP *P**(-05)dp 221 : IF ((EPSE+DELTLT).GE. 1.D-20 )THEN 222 : DR=0.5*MP*(EPSE+DELTLT)**(-0.5) 223 : ELSE 224 : DR=0.D0 225 : END IF 226 : DFDELTL=-DR 227 : DO ISTRS=1,NSIG 228 : DFDELTL=DFDELTL-DFDS(ISTRS)*(DFDS(ISTRS)*DEUXMU+TRDFDS*LAMB) 229 : END DO 230 : DO ISTRS=4,NSTRS 231 : DFDELTL=DFDELTL - 2.d0*DEUXMU * DFDS(ISTRS)**2 232 : END DO 233 : IF (LCP) THEN 234 : DFDELTL=DFDELTL-DFDS(3)*TRDFDS*LAMB 235 : DFDEPS3=LAMB*TRDFDS + (LAMB+DEUXMU)*DFDS(3) 236 : DS3DELTL=-LAMB*TRDFDS 237 : DS3DEPS3=LAMB + DEUXMU 238 : DELTL=(SIGE0(3)*DFDEPS3 - F0*DS3DEPS3)/ 239 : & (DS3DEPS3*DFDELTL - DS3DELTL*DFDEPS3) 240 : DEPS3=-(SIGE0(3)-DS3DELTL*DELTL)/DS3DEPS3 241 : EPS3=EPS3+DEPS3 242 : DELTLT=DELTLT+DELTL 243 : DO ISTRS=1,NSIG 244 : SIGE0(ISTRS)=SIGE0(ISTRS)+LAMB*DEPS3 245 : END DO 246 : SIGE0(3)=SIGE0(3)+(LAMB+DEUXMU)*DEPS3-LAMB*DELTL*TRDFDS 247 : ELSE 248 : DELTL=-F0/DFDELTL 249 : DELTLT=DELTLT+DELTL 250 : END IF 251 : DO ISTRS=1,NSIG 252 : SIGE0(ISTRS)=SIGE0(ISTRS) 253 : & -DELTL*(DEUXMU*DFDS(ISTRS)+LAMB*TRDFDS) 254 : END DO 255 : DO ISTRS=4,NSTRS 256 : SIGE0(ISTRS)=SIGE0(ISTRS)-DELTL*DEUXMU*DFDS(ISTRS) 257 : END DO 258 : ELSE 259 : DELTL=0.D0 260 : END IF 261 : R=R + DR*DELTL 262 : DO ISTRS=1,3 263 : DEFP(ISTRS)=DEFP(ISTRS)+DELTL*DFDS(ISTRS) 264 : END DO 265 : DO ISTRS=4,NSTRS 266 : DEFP(ISTRS)=DEFP(ISTRS)+2.D0*DELTL*DFDS(ISTRS) 267 : END DO 268 : if (limpr) then 269 : print*,'ITER=',ITER,'F0=',F0,'DELTL=', DELTL 270 : print*,'DELTLT=',DELTLT,'R=',R 271 : do istrs=1,nstrs 272 : print*,'SIGMA(',ISTRS,')=', SIGE0(ISTRS) 273 : end do 274 : end if 275 : END DO 276 : * EN CAS DE NON CONVERGENCE: ON MET KERRE a 1 277 : IF (.NOT. CONV) THEN 278 : KERRE=1 279 : END IF 280 : * PRINT*,'CONVERGENCE OBTENUE APRES ',ITER, 'ITERATIONS' 281 : * 282 : * on remplit les variables finales 283 : * 284 : VARF(1)=EPSE+DELTLT 285 : VARF(2)=R 286 : * 287 : * ON REND LES CONTRAINTES 288 : * 289 : DO ISTRS=1,NSTRS 290 : SIGF(ISTRS)=SIGE0(ISTRS) 291 : * PRINT*,'SIGF(',ISTRS,')=',SIGF(ISTRS) 292 : END DO 293 : * 294 : * ON REND L'INCREMENT DE DEFORMATIONS PLASTIQUES 295 : * 296 : 297 : SEGSUP WRKK1 298 : RETURN 299 : END 300 : * 301 : 302 : 303 :
© Cast3M 2003 - All rights reserved.
Disclaimer