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