1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine incplas_lin(elconloc,plconloc,xstate,xstateini, 20 & elas,emec,ithermal,icmd,beta,stre,vj,kode, 21 & ielas,amat,t1l,dtime,time,ttime,iel,iint,nstate_,mi, 22 & eloc,pgauss,nmethod,pnewdt,depvisc) 23! 24! calculates stiffness and stresses for the incremental plasticity 25! material law 26! 27! icmd=3: calculates stress at mechanical strain 28! else: calculates stress at mechanical strain and the stiffness 29! matrix 30! 31! This routine is meant for small strains. Reference: 32! Dhondt, Guido; The Finite Element Method for Three-dimensional 33! Thermomechanical Applications (Wiley, 2004), Section 5.3 34! 35 implicit none 36! 37 character*80 amat 38! 39 integer ithermal(*),icmd,i,k,l,m,n,kode,ivisco,ielastic,kel(4,21), 40 & niso,nkin,ielas,iel,iint,nstate_,mi(*),id,leximp,lend,layer, 41 & kspt,kstep,kinc,iloop,nmethod,user_hardening,user_creep 42! 43 real*8 elconloc(*),elas(21),emec(6),beta(6),stre(6), 44 & vj,plconloc(802),stbl(6),stril(6),xitril(6), 45 & ee,un,um,al,cop,dxitril,xn(3,3),epl(6),c1,c2,c3,c4,c7, 46 & c8,ftrial,xiso(200),yiso(200),xkin(200),ykin(200), 47 & fiso,dfiso,fkin,dfkin,fiso0,fkin0,ep,t1l,dtime, 48 & epini,a1,dsvm,xxa,xxn,dkl(3,3),el(6),tracee,traces, 49 & dcop,time,ttime,eloc(6),xstate(nstate_,mi(1),*), 50 & xstateini(nstate_,mi(1),*),decra(5),deswa(5),serd, 51 & esw(2),ec(2),p,qtild,predef(1),dpred(1),timeabq(2),pgauss(3), 52 & dtemp,pnewdt,um2,depvisc 53! 54! 55 kel=reshape((/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3, 56 & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3, 57 & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3, 58 & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/)) 59 dkl=reshape((/1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/), 60 & (/3,3/)) 61! 62 leximp=1 63 lend=2 64 user_creep=0 65! 66! localizing the plastic fields 67! 68 do i=1,6 69 epl(i)=xstateini(1+i,iint,iel) 70 stbl(i)=xstateini(7+i,iint,iel) 71 enddo 72 epini=xstateini(1,iint,iel) 73! 74 ee=elconloc(1) 75 un=elconloc(2) 76 um2=ee/(1.d0+un) 77 al=um2*un/(1.d0-2.d0*un) 78 um=um2/2.d0 79! 80 ep=epini 81! 82! check whether the user activated a viscous calculation 83! (*VISCO instead of *STATIC) 84! 85 if((nmethod.ne.1).or.(ithermal(1).eq.3)) then 86 ivisco=1 87 else 88 ivisco=0 89 endif 90! 91! check for user subroutines 92! 93 if((plconloc(801).lt.0.8d0).and.(plconloc(802).lt.0.8d0)) then 94 user_hardening=1 95 else 96 user_hardening=0 97 endif 98 if((kode.eq.-52).and.(ivisco.eq.1)) then 99 if(elconloc(3).lt.0.d0) then 100 user_creep=1 101 else 102 user_creep=0 103 xxa=elconloc(3)*(ttime+time)**elconloc(5) 104 if(xxa.lt.1.d-20) xxa=1.d-20 105 xxn=elconloc(4) 106 a1=xxa*dtime 107 endif 108 endif 109! 110! trial elastic strain 111! 112 do i=1,6 113 el(i)=emec(i)-epl(i) 114 enddo 115! 116! trial stress 117! 118 tracee=el(1)+el(2)+el(3) 119 do i=1,6 120 stre(i)=um2*el(i)-beta(i) 121 enddo 122 do i=1,3 123 stre(i)=stre(i)+al*tracee 124 enddo 125! 126! trial deviatoric stress 127! 128 traces=(stre(1)+stre(2)+stre(3))/3.d0 129 do i=1,3 130 stril(i)=stre(i)-traces 131 enddo 132 do i=4,6 133 stril(i)=stre(i) 134 enddo 135! 136! subtracting the back stress -> trial radius vector 137! 138 do i=1,6 139 xitril(i)=stril(i)-stbl(i) 140 enddo 141! 142! size of trial radius vector 143! 144 dxitril=0.d0 145 do i=1,3 146 dxitril=dxitril+xitril(i)*xitril(i) 147 enddo 148 do i=4,6 149 dxitril=dxitril+2.d0*xitril(i)*xitril(i) 150 enddo 151 dxitril=dsqrt(dxitril) 152! 153! restoring the hardening curves for the actual temperature 154! plconloc contains the true stresses. By multiplying by 155! the Jacobian, yiso and ykin are Kirchhoff stresses, as 156! required by the hyperelastic theory (cf. Simo, 1988). 157! 158 niso=int(plconloc(801)) 159 nkin=int(plconloc(802)) 160 if(niso.ne.0) then 161 do i=1,niso 162 xiso(i)=plconloc(2*i-1) 163 yiso(i)=plconloc(2*i) 164 enddo 165 endif 166 if(nkin.ne.0) then 167 do i=1,nkin 168 xkin(i)=plconloc(399+2*i) 169 ykin(i)=plconloc(400+2*i) 170 enddo 171 endif 172! 173! if no viscous calculation is performed a pure creep calculatino 174! (without plasticity) is reduced to an elastic calculation 175! 176 ielastic=0 177 if(ivisco.eq.0) then 178 if(niso.eq.2) then 179 if((dabs(yiso(1)).lt.1.d-10).and. 180 & (dabs(yiso(2)).lt.1.d-10)) then 181 ielastic=1 182 endif 183 endif 184 endif 185! 186! check for yielding 187! 188 if(user_hardening.eq.1) then 189 call uhardening(amat,iel,iint,t1l,epini,ep,dtime, 190 & fiso,dfiso,fkin,dfkin) 191 else 192 if(niso.ne.0) then 193 call ident(xiso,ep,niso,id) 194 if(id.eq.0) then 195 fiso=yiso(1) 196 elseif(id.eq.niso) then 197 fiso=yiso(niso) 198 else 199 dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 200 fiso=yiso(id)+dfiso*(ep-xiso(id)) 201 endif 202 elseif(nkin.ne.0) then 203 fiso=ykin(1) 204 else 205 fiso=0.d0 206 endif 207 endif 208! 209 c1=2.d0/3.d0 210 c2=dsqrt(c1) 211! 212 ftrial=dxitril-c2*fiso 213 if((ftrial.le.1.d-10).or.(ielas.eq.1).or.(ielastic.eq.1) 214 & .or.(dtime.lt.1.d-30)) then 215! 216! updating the plastic fields 217! 218 do i=1,6 219 xstate(1+i,iint,iel)=epl(i) 220 xstate(7+i,iint,iel)=stbl(i) 221 enddo 222 xstate(1,iint,iel)=ep 223! 224 if(icmd.ne.3) then 225 elas(1)=al+um2 226 elas(2)=al 227 elas(3)=al+um2 228 elas(4)=al 229 elas(5)=al 230 elas(6)=al+um2 231 elas(7)=0.d0 232 elas(8)=0.d0 233 elas(9)=0.d0 234 elas(10)=um 235 elas(11)=0.d0 236 elas(12)=0.d0 237 elas(13)=0.d0 238 elas(14)=0.d0 239 elas(15)=um 240 elas(16)=0.d0 241 elas(17)=0.d0 242 elas(18)=0.d0 243 elas(19)=0.d0 244 elas(20)=0.d0 245 elas(21)=um 246 endif 247! 248 return 249 endif 250! 251! plastic deformation 252! 253! calculating the consistency parameter 254! 255 iloop=0 256 cop=0.d0 257! 258 loop: do 259 iloop=iloop+1 260 ep=epini+c2*cop 261! 262 if(user_hardening.eq.1) then 263 call uhardening(amat,iel,iint,t1l,epini,ep,dtime, 264 & fiso,dfiso,fkin,dfkin) 265 else 266 if(niso.ne.0) then 267 call ident(xiso,ep,niso,id) 268 if(id.eq.0) then 269 fiso=yiso(1) 270 dfiso=0.d0 271 elseif(id.eq.niso) then 272 fiso=yiso(niso) 273 dfiso=0.d0 274 else 275 dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 276 fiso=yiso(id)+dfiso*(ep-xiso(id)) 277 endif 278 elseif(nkin.ne.0) then 279 fiso=ykin(1) 280 dfiso=0.d0 281 else 282 fiso=0.d0 283 dfiso=0.d0 284 endif 285! 286 if(nkin.ne.0) then 287 call ident(xkin,ep,nkin,id) 288 if(id.eq.0) then 289 fkin=ykin(1) 290 dfkin=0.d0 291 elseif(id.eq.nkin) then 292 fkin=ykin(nkin) 293 dfkin=0.d0 294 else 295 dfkin=(ykin(id+1)-ykin(id))/(xkin(id+1)-xkin(id)) 296 fkin=ykin(id)+dfkin*(ep-xkin(id)) 297 endif 298 elseif(niso.ne.0) then 299 fkin=yiso(1) 300 dfkin=0.d0 301 else 302 fkin=0.d0 303 dfkin=0.d0 304 endif 305 endif 306! 307 if(dabs(cop).lt.1.d-10) then 308 fiso0=fiso 309 fkin0=fkin 310 endif 311! 312 if((kode.eq.-51).or.(ivisco.eq.0)) then 313 dcop=(dxitril-c2*(fiso+fkin-fkin0)-um2*cop)/ 314 & (um2+c1*(dfiso+dfkin)) 315 else 316 if(user_creep.eq.1) then 317 if(ithermal(1).eq.0) then 318 write(*,*) '*ERROR in incplas: no temperature defined' 319 call exit(201) 320 endif 321 timeabq(1)=time 322 timeabq(2)=ttime+time 323! 324! Eqn. (5.139) 325! 326 qtild=(dxitril-um2*cop)/c2-fiso-(fkin-fkin0) 327! 328! the Von Mises stress must be positive 329! 330 if(qtild.lt.1.d-10) qtild=1.d-10 331 ec(1)=epini 332 call creep(decra,deswa,xstateini(1,iint,iel),serd,ec, 333 & esw,p,qtild,t1l,dtemp,predef,dpred,timeabq,dtime, 334 & amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt, 335 & kstep,kinc) 336 dsvm=1.d0/decra(5) 337 dcop=(decra(1)-c2*cop)/ 338 & (c2*(decra(5)*(dfiso+um*(3.d0+dfkin/um))+1.d0)) 339 else 340 qtild=(dxitril-um2*cop)/c2-fiso-(fkin-fkin0) 341! 342! the Von Mises stress must be positive 343! 344 if(qtild.lt.1.d-10) qtild=1.d-10 345 decra(1)=a1*qtild**xxn 346 decra(5)=xxn*decra(1)/qtild 347 dsvm=1.d0/decra(5) 348 dcop=(decra(1)-c2*cop)/ 349 & (c2*(decra(5)*(dfiso+um*(3.d0+dfkin/um))+1.d0)) 350 endif 351 endif 352 cop=cop+dcop 353! 354 if((dabs(dcop).lt.cop*1.d-4).or. 355 & (dabs(dcop).lt.1.d-10)) exit 356! 357! check for endless loops or a negative consistency 358! parameter 359! 360 if((iloop.gt.15).or.(cop.le.0.d0)) then 361 pnewdt=0.25d0 362 return 363 endif 364! 365 enddo loop 366! 367! updating the equivalent plastic strain 368! 369 ep=epini+c2*cop 370! 371! updating the back stress 372! 373 c7=c2*(fkin-fkin0)/dxitril 374 do i=1,6 375 stbl(i)=stbl(i)-c7*xitril(i) 376 enddo 377! 378! updating the plastic strain tensor 379! 380 c7=cop/dxitril 381 do i=1,6 382 epl(i)=epl(i)+c7*xitril(i) 383 enddo 384! 385! trial elastic strain 386! 387 do i=1,6 388 el(i)=emec(i)-epl(i) 389 enddo 390! 391! trial stress 392! 393 tracee=el(1)+el(2)+el(3) 394 do i=1,6 395 stre(i)=um2*el(i)-beta(i) 396 enddo 397 do i=1,3 398 stre(i)=stre(i)+al*tracee 399 enddo 400! 401! calculating the local stiffness matrix 402! 403 if(icmd.ne.3) then 404 c7=um2**2*cop/dxitril 405! 406 if((kode.eq.-51).or.(ivisco.eq.0)) then 407 c8=um2**2/(um2+c1*(dfiso+dfkin)) 408 else 409 c8=um2**2/(um2+c1*(dfiso+dfkin+1.d0/decra(5))) 410 endif 411! 412 c1=um-c7/2.d0 413 c2=al+c7/3.d0 414 c3=(c7-c8)/(dxitril**2) 415! 416 xn(1,1)=xitril(1) 417 xn(2,2)=xitril(2) 418 xn(3,3)=xitril(3) 419 xn(1,2)=xitril(4) 420 xn(1,3)=xitril(5) 421 xn(2,3)=xitril(6) 422 xn(2,1)=xn(1,2) 423 xn(3,1)=xn(1,3) 424 xn(3,2)=xn(2,3) 425! 426 do i=1,21 427 k=kel(1,i) 428 l=kel(2,i) 429 m=kel(3,i) 430 n=kel(4,i) 431 elas(i)=c1*(dkl(k,m)*dkl(l,n)+dkl(k,n)*dkl(l,m)) 432 & +c2*dkl(k,l)*dkl(m,n) 433 & +c3*xn(k,l)*xn(m,n) 434 enddo 435 endif 436! 437! updating the plastic fields 438! 439 do i=1,6 440 xstate(1+i,iint,iel)=epl(i) 441 xstate(7+i,iint,iel)=stbl(i) 442 enddo 443 xstate(1,iint,iel)=ep 444! 445! maximum difference in the equivalent viscoplastic strain 446! in this increment based on the viscoplastic strain rate at 447! the start and the end of the increment 448! 449 if(ivisco.eq.1) then 450 depvisc=max(depvisc,abs(decra(1)-dtime*xstateini(14,iint,iel))) 451 xstate(14,iint,iel)=decra(1)/dtime 452 endif 453! 454! maximum difference in the equivalent viscoplastic strain 455! in this increment based on the viscoplastic strain rate at 456! the start and the end of the increment 457! 458 if(ivisco.eq.1) then 459 depvisc=max(depvisc,abs(decra(1)-dtime*xstateini(14,iint,iel))) 460 xstate(14,iint,iel)=decra(1)/dtime 461 endif 462c depvisc=max(depvisc,ep-epini) 463! 464 return 465 end 466