1 : C PRVIK2 SOURCE PV 11/03/07 21:17:56 6885 2 : SUBROUTINE PRVIK2(SIG0,NSTRS,DEPST,VAR0,XMAT,NCOMAT,XCAR,ICARA, 3 : 1 NVARI,SIGF,VARF,DEFP,MFR1,DDAUX,CMATE,VALMAT,VALCAR,N2EL, 4 : 2 N2PTEL,NBPGAU,IFOU,IB,IGAU,EPAIST,MELE,NPINT,NBGMAT, 5 : 3 NELMAT,SECT,LHOOK,TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU, 6 : 4 CRIGI,DSIGT,KERRE,DT) 7 : comm ****************************************************************** 8 : comm * execute le calcul annonce au niveau du p_gauss * 9 : comm * loi viscoplastique : le comportement anelastique est decrit * 10 : comm * par un mecanisme "plastique" (ecrouissage cinematique) * 11 : comm * et un mecanisme "visqueux" (Maxwell) * 12 : comm * en parallele, ce au-dela du seuil de plasticite * 13 : comm * * 14 : comm * il faut donc gerer 2 jeux de variables cachees tensorielles * 15 : comm * * 16 : comm * cette loi est derivee des travaux deja realises a propos du * 17 : comm * comportement mecanique du polyethylene * 18 : comm * (these Kichenin, Ecole Polytechnique, 1992) * 19 : comm * algorithme Nguyen Quoc Son, mise au point J Frelat * 20 : comm * $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$* 21 : comm * RECOMMANDATION : * 22 : comm * la chasse au bugs est bienvenue MAIS PAS * 23 : comm * les ameliorations ou nouveaux developpements dans * 24 : comm * l inspiration de VISK2 (rapport SEMT/LAMS 98 Kichenin, Charras)* 25 : comm * Il est preferable d une part de creer une nouvelle loi, d autre* 26 : comm * part de RECOPIER dans un AUTRE fichier les lignes voulues * 27 : comm * MERCI. JK * 28 : comm ****************************************************************** 29 : IMPLICIT INTEGER(I-N) 30 : IMPLICIT REAL*8(A-H,O-Z) 31 : -INC CCOPTIO 32 : DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*), 33 : & VARF(*),DEFP(*) 34 : DIMENSION RSIG0(6),RDSIGT(6),RDEFP(6),RSIGF(6),DSIGT(*) 35 : DIMENSION VALMAT(*),VALCAR(*) 36 : DIMENSION TXR(IDIM,*),CRIGI(12) 37 : DIMENSION DDAUX(LHOOK,*),DDHOMU(LHOOK,*) 38 : DIMENSION XLOC(IDIM,*),XGLOB(IDIM,*) 39 : DIMENSION D1HOOK(LHOOK,*),ROTHOO(LHOOK,*) 40 : CHARACTER*8 CMATE 41 : real*8 lambdav 42 : logical plast 43 : dimension sigpca(6),sigppa(6),sig2(6) 44 : c.. 45 : c.. VOIR PB IFOUR IFOURB IFOU 46 : IF (IFOUR.NE.-2.AND.IFOUR.NE.-1.AND.IFOUR.NE.-3.AND.IFOUR.NE.0) 47 : & THEN 48 : iuv1 = 7 49 : ELSE 50 : iuv1 = 5 51 : ENDIF 52 : c.. 53 : if((1.+xmat(2)).eq.0.or.(1-2.*xmat(2)).eq.0) then 54 : interr(5)=1 55 : kerre = 82 56 : return 57 : endif 58 : * suppose fondalement isotropie du comportement 59 : x2mu= xmat(1)/(1.+xmat(2)) 60 : x2muv=xmat(8)/(1.+xmat(2)) 61 : seuil=xmat(5)/1.73205080757 62 : ecrou=xmat(6) 63 : if (x2mu.le.0.or.seuil.lt.0.or.ecrou.lt.0) then 64 : interr(5) = 4 65 : kerre = 82 66 : return 67 : endif 68 : c.. initialise ecrouissage mecanisme plastique 69 : c.. l indice 1 est la deformation plastique equivalente 70 : varf(1) = var0(1) 71 : varf(2) = var0(2) 72 : varf(3) = var0(3) 73 : varf(4) = var0(4) 74 : varf(5) = var0(5) 75 : if (ifour.gt.0) then 76 : varf(6) = var0(6) 77 : varf(7) = var0(7) 78 : endif 79 : c mecanisme visqueux 80 : svxx=var0(iuv1+1) 81 : svyy=var0(iuv1+2) 82 : svzz=var0(iuv1+3) 83 : svxy=var0(iuv1+4) 84 : if (ifour.gt.0) then 85 : svxz=var0(iuv1+5) 86 : svyz=var0(iuv1+6) 87 : endif 88 : c.. 89 : c 90 : C CALCUL DE L'INCREMENT DE CONTRAINTES "elastique" 91 : C 91 : CALL CALSIG(DEPST,DDAUX,NSTRS,CMATE,VALMAT,VALCAR, 92 : 1 N2EL,N2PTEL,MFR1,IFOU,IB,IGAU,EPAIST,NBPGAU, 93 : 2 MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,XLOC, 94 : 3 XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD) 95 : * 96 : IF(IRTD.NE.1) THEN 97 : KERRE=69 98 : GOTO 1000 99 : ENDIF 100 : 101 : 102 : plast = .false. 103 : c.. composante plastique + increment elastique : cinematique admissible 104 : sigpca(1)= sig0(1) - svxx + DSIGT(1) 105 : sigpca(2)= sig0(2) - svyy + DSIGT(2) 106 : sigpca(3)= sig0(3) - svzz + DSIGT(3) 107 : sigpca(4)= sig0(4) - svxy + DSIGT(4) 108 : if (ifour.gt.0)then 109 : sigpca(5)= sig0(5) - svxz + DSIGT(5) 110 : sigpca(6)= sig0(6) - svyz + DSIGT(6) 111 : endif 112 : trsigp=0.33333*(sigpca(1)+sigpca(2)+sigpca(3)) 113 : c.. 114 : c.. evalue cette composante (deviateur) dans le convexe definissant l'ecrouissage 115 : sig2(1)= sigpca(1)-trsigp-varf(2) 116 : sig2(2)= sigpca(2)-trsigp-varf(3) 117 : sig2(3)= sigpca(3)-trsigp-varf(4) 118 : sig2(4)= sigpca(4)-varf(5) 119 : if (ifour.gt.0)then 120 : sig2(5)= sigpca(5)-varf(6) 121 : sig2(6)= sigpca(6)-varf(7) 122 : endif 123 : c.. calcule la norme de Von Mises 124 : fsig2 = sig2(1)*sig2(1) + sig2(2)*sig2(2) + sig2(3)*sig2(3) 125 : fsig2 = fsig2 + 2.d0*sig2(4)*sig2(4) 126 : if(ifour.gt.0)then 127 : fsig2 = fsig2 + 2.d0*sig2(5)*sig2(5) + 2.d0*sig2(6)*sig2(6) 128 : endif 129 : fsig2 = dsqrt(0.5d0*fsig2) 130 : 131 : 132 : c.. tests 133 : c a l interieur du convexe ou cas pathologique 134 : if ((fsig2-seuil).le.0.or.fsig2.eq.0) then 135 : plast=.false. 136 : coef = 0. 137 : defp(1) = 0. 138 : defp(2) = 0. 139 : defp(3) = 0. 140 : defp(4) = 0. 141 : if (ifour.gt.0) then 142 : defp(5) = 0. 143 : defp(6) = 0. 144 : endif 145 : goto 100 146 : endif 147 : c.. sinon le seuil est franchi 148 : plast = .true. 149 : c.. et il est necessaire de faire un calcul d ecoulement 150 : coef=(fsig2-seuil)/(x2mu+ecrou)/fsig2 151 : c.. increment de la deformation plastique 152 : defp(1) = coef * sig2(1) 153 : defp(2) = coef * sig2(2) 154 : defp(3) = coef * sig2(3) 155 : defp(4) = coef * sig2(4) * 2. 156 : if (ifour.gt.0) then 157 : defp(5) = coef * sig2(5) * 2. 158 : defp(6) = coef * sig2(6) * 2. 159 : endif 160 : 100 continue 161 : c.. increment de la variable d ecrouissage 162 : varf(1) = varf(1) + coef*fsig2*1.414213562 163 : varf(2) = varf(2) + ecrou*defp(1) 164 : varf(3) = varf(3) + ecrou*defp(2) 165 : varf(4) = varf(4) + ecrou*defp(3) 166 : varf(5) = varf(5) + ecrou*defp(4) / 2. 167 : if (ifour.gt.0) then 168 : varf(6) = varf(6) + ecrou*defp(5) / 2. 169 : varf(7) = varf(7) + ecrou*defp(6) / 2. 170 : endif 171 : c.. contrainte plastiquement admissible 172 : sigppa(1) = sigpca(1) - x2mu*defp(1) 173 : sigppa(2) = sigpca(2) - x2mu*defp(2) 174 : sigppa(3) = sigpca(3) - x2mu*defp(3) 175 : sigppa(4) = sigpca(4) - x2mu*defp(4) / 2. 176 : if(ifour.gt.0) then 177 : sigppa(5) = sigpca(5) - x2mu*defp(5) / 2. 178 : sigppa(6) = sigpca(6) - x2mu*defp(6) / 2. 179 : endif 180 : 181 : fsigpa = sigppa(1)*sigppa(1) + sigppa(2)*sigppa(2) 182 : & + sigppa(3)*sigppa(3) + 2.d0*sigppa(4)*sigppa(4) 183 : if (ifour.gt.0) then 184 : fsigpa = fsigpa + 2.d0*sigppa(5)*sigppa(5) 185 : & + 2.d0*sigppa(6)*sigppa(6) 186 : endif 187 : fsigpa = dsqrt(0.5d0*fsigpa) 188 : c.. 189 : 190 : 300 CONTINUE 191 : 192 : c.. 193 : c.. traite le mecanisme visqueux 194 : fsv1 = 0. 195 : comm *................................................................* 196 : comm * viscosite non-lineaire ... si devis est le parametre cache * 197 : comm * increment de deformation-visqueuse, et sigvis le parametre * 198 : comm * contrainte-visqueuse initial, la loi de comportement s ecrit * 199 : comm * eta*dev = (sigvis + loiv*(du-dev))*|sigvis+loiv*(du-dev)|**m * 200 : comm * (m = n-1) le mecanisme visqueux etant somme toute secondaire * 201 : comm * l approche bricolee consiste a evaluer * 202 : comm * fsv1 = |sigvis+loiv*du|**m, dev=0 * 203 : comm * puis ecrire sigv-final = (sigvis + loiv*(du-dev)) * coefv * 204 : comm * le coefv etant evalue a partir de fsv1 et valant * 205 : comm * coefv = 1./(1+x2muv/eta) quand m=0 * 206 : comm *................................................................* 207 : c.. on suppose le comportement visqueux homogene et isotrope 208 : lambdav=xmat(8)*xmat(2)/(1+xmat(2))/(1-2.*xmat(2)) 209 : x2muv=xmat(8)/(1.+xmat(2)) 210 : if(x2muv.lt.0) then 211 : interr(5)=5 212 : kerre = 82 213 : return 214 : endif 215 : m = nint(xmat(9)) - 1 216 : c.. coefficient applique a la branche visqueuse : mecanisme actif 217 : c.. au-dela du seuil, inactif sinon 218 : if (x2mu.gt.0) then 219 : xambda2 = x2muv/x2mu 220 : endif 221 : if (xambda2.le.0) xambda2 = 0. 222 : c.. viscosite 223 : if (xmat(7).le.0.) then 224 : c.. mecanisme visqueux inactif, 225 : varf(iuv1+1)=0. 226 : varf(iuv1+2)=0. 227 : varf(iuv1+3)=0. 228 : varf(iuv1+4)=0. 229 : if (ifour.gt.0)then 230 : varf(iuv1+5)=0. 231 : varf(iuv1+6)=0. 232 : endif 233 : goto 400 234 : endif 235 : c.. le mecanisme n'est efficace qu a partir d un pas de temps assez grand 236 : c.. et si l on est dans le domaine plastique 237 : if (dt.gt.1.e-15) then 238 : eta=xmat(7)/dt 239 : x2muv2=x2muv+eta 240 : coedev=1./x2muv2 241 : coesph=lambdav/x2muv2/(3.*lambdav+x2muv2) 242 : c.. 243 : c.. incremente "elastiquement" la contrainte-visqueuse 244 : svxx=svxx+x2muv * defp(1) 245 : svyy=svyy+x2muv * defp(2) 246 : svzz=svzz+x2muv * defp(3) 247 : svxy=svxy+x2muv * defp(4) 248 : if(ifour.gt.0)then 249 : svxz=svxz+x2muv * defp(5) 250 : svyz=svyz+x2muv * defp(6) 251 : endif 252 : c.. 253 : c.. deduit l'increment de deformation-visqueuse 254 : trsv=svxx+svyy+svzz 255 : svxx=eta*(coedev*svxx-coesph*trsv) 256 : svyy=eta*(coedev*svyy-coesph*trsv) 257 : svzz=eta*(coedev*svzz-coesph*trsv) 258 : evxy=eta*coedev*svxy 259 : if(ifour.gt.0)then 260 : svxz=eta*coedev*svxz 261 : svyz=eta*coedev*svyz 262 : endif 263 : c.. 264 : fsv1 = svxx*svxx + svyy*svyy + svzz*svzz + 2.*svxy*svxy 265 : if(ifour.gt.0) then 266 : fsv1= fsv1 + 2.*svxz*svxz+2.*svyz*svyz 267 : fsv1 = fsv1/eta 268 : endif 269 : c 270 : c.. coef multiplicateur 271 : if (fsv1.gt.0.and.m.gt.0) fsv1 = (dsqrt(0.5*fsv1))**m 272 : if (fsv1.gt.0.and.m.lt.0) fsv1 = (1./dsqrt(0.5*fsv1))**(-1.*m) 273 : if (fsv1.gt.0.and.m.eq.0) fsv1 = 1. 274 : * if(igau.eq.4) write(6,*) ' fsv1 ', fsv1 275 : if (fsv1.eq.0) then 276 : c.. composante visqueuse nulle 277 : varf(iuv1+1) = 0. 278 : varf(iuv1+2) = 0. 279 : varf(iuv1+3) = 0. 280 : varf(iuv1+4) = 0. 281 : if (ifour.gt.0) then 282 : varf(iuv1+5) = 0. 283 : varf(iuv1+6) = 0. 284 : endif 285 : goto 400 286 : endif 287 : c.. coeff reducteur 288 : coefv = fsv1/(1. + (m+1)*x2muv*fsv1/eta) 289 : varf(iuv1+1)=svxx*coefv 290 : varf(iuv1+2)=svyy*coefv 291 : varf(iuv1+3)=svzz*coefv 292 : varf(iuv1+4)=svxy*coefv 293 : if (ifour.gt.0)then 294 : varf(iuv1+5)=svxz*coefv 295 : varf(iuv1+6)=svyz*coefv 296 : endif 297 : 298 : c.. quand le pas de temps est petit la variable cachee est inchangee 299 : else 300 : varf(iuv1+1)=var0(iuv1+1) 301 : varf(iuv1+2)=var0(iuv1+2) 302 : varf(iuv1+3)=var0(iuv1+3) 303 : varf(iuv1+4)=var0(iuv1+4) 304 : if (ifour.gt.0)then 305 : varf(iuv1+5)=var0(iuv1+5) 306 : varf(iuv1+6)=var0(iuv1+6) 307 : endif 308 : endif 309 : comm *................................................................* 310 : comm * calcul de la contrainte finale * 311 : comm *................................................................* 312 : 400 continue 313 : sigf(1) = sigppa(1) + varf(iuv1+1) 314 : sigf(2) = sigppa(2) + varf(iuv1+2) 315 : sigf(3) = sigppa(3) + varf(iuv1+3) 316 : sigf(4) = sigppa(4) + varf(iuv1+4) 317 : if(ifour.gt.0)then 318 : sigf(5) = sigppa(5) + varf(iuv1+5) 319 : sigf(6) = sigppa(6) + varf(iuv1+6) 320 : endif 321 : return 322 : 323 : 1000 continue 324 : end 325 : 326 : 327 :
© Cast3M 2003 - All rights reserved.
Disclaimer