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 umat_aniso_creep(amat,iel,iint,kode,elconloc,emec, 20 & emec0,beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime, 21 & icmd,ielas,mi,nstate_,xstateini,xstate,stre,stiff,iorien, 22 & pgauss,orab,nmethod,pnewdt,depvisc) 23! 24! calculates stiffness and stresses for a user defined material 25! law 26! 27! icmd=3: calculates stress at mechanical strain 28! else: calculates stress at mechanical strain and the stiffness 29! matrix 30! 31! INPUT: 32! 33! amat material name 34! iel element number 35! iint integration point number 36! 37! kode material type (-100-#of constants entered 38! under *USER MATERIAL): can be used for materials 39! with varying number of constants 40! 41! elconloc(21) user defined constants defined by the keyword 42! card *USER MATERIAL (max. 21, actual # = 43! -kode-100), interpolated for the 44! actual temperature t1l 45! 46! emec(6) Lagrange mechanical strain tensor (component order: 47! 11,22,33,12,13,23) at the end of the increment 48! (thermal strains are subtracted) 49! emec0(6) Lagrange mechanical strain tensor at the start of the 50! increment (thermal strains are subtracted) 51! beta(6) residual stress tensor (the stress entered under 52! the keyword *INITIAL CONDITIONS,TYPE=STRESS) 53! 54! xokl(3,3) deformation gradient at the start of the increment 55! voj Jacobian at the start of the increment 56! xkl(3,3) deformation gradient at the end of the increment 57! vj Jacobian at the end of the increment 58! 59! ithermal 0: no thermal effects are taken into account: for 60! creep this does not make sense. 61! >0: thermal effects are taken into account (triggered 62! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) 63! t1l temperature at the end of the increment 64! dtime time length of the increment 65! time step time at the end of the current increment 66! ttime total time at the start of the current step 67! 68! icmd not equal to 3: calculate stress and stiffness 69! 3: calculate only stress 70! ielas 0: no elastic iteration: irreversible effects 71! are allowed 72! 1: elastic iteration, i.e. no irreversible 73! deformation allowed 74! 75! mi(1) max. # of integration points per element in the 76! model 77! nstate_ max. # of state variables in the model 78! 79! xstateini(nstate_,mi(1),# of elements) 80! state variables at the start of the increment 81! xstate(nstate_,mi(1),# of elements) 82! state variables at the end of the increment 83! 84! stre(6) Piola-Kirchhoff stress of the second kind 85! at the start of the increment 86! 87! iorien number of the local coordinate axis system 88! in the integration point at stake (takes the value 89! 0 if no local system applies) 90! pgauss(3) global coordinates of the integration point 91! orab(7,*) description of all local coordinate systems. 92! If a local coordinate system applies the global 93! tensors can be obtained by premultiplying the local 94! tensors with skl(3,3). skl is determined by calling 95! the subroutine transformatrix: 96! call transformatrix(orab(1,iorien),pgauss,skl) 97! 98! 99! OUTPUT: 100! 101! xstate(nstate_,mi(1),# of elements) 102! updated state variables at the end of the increment 103! stre(6) Piola-Kirchhoff stress of the second kind at the 104! end of the increment 105! stiff(21): consistent tangent stiffness matrix in the material 106! frame of reference at the end of the increment. In 107! other words: the derivative of the PK2 stress with 108! respect to the Lagrangian strain tensor. The matrix 109! is supposed to be symmetric, only the upper half is 110! to be given in the same order as for a fully 111! anisotropic elastic material (*ELASTIC,TYPE=ANISO). 112! Notice that the matrix is an integral part of the 113! fourth order material tensor, i.e. the Voigt notation 114! is not used. 115! 116 implicit none 117! 118 integer exitcriterion 119! 120 character*80 amat 121! 122 integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien, 123 & i,j,ipiv(6),info,neq,lda,ldb,j1,j2,j3,j4,j5,j6,j7,j8, 124 & nrhs,iplas,kel(4,21),iloop,leximp,lend,layer,kspt,kstep, 125 & kinc,ii,nmethod 126! 127 real*8 ep0(6),epqini,ep(6),b,Pn(6),dg,ddg,c(21),x(21),cm1(21), 128 & stri(6),htri,sg(6),r(13),ee(6),dd,gl(6,6),gr(6,6),c0,c1,c2, 129 & skl(3,3),gcreep,gm1,ya(3,3,3,3),dsg,detc,strinv,depvisc, 130 & depq,svm,dsvm,dg1,dg2,fu,fu1,fu2,expon,ec(2),pnewdt, 131 & timeabq(2),r1(13),ep1(6),gl1(6,6),sg1(6),ckl(3,3), 132 & elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6), 133 & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), 134 & time,ttime,decra(5),deswa(5),serd,esw(2),p,predef(1),dpred(1), 135 & dtemp,xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*) 136! 137 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, 138 & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3, 139 & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3, 140 & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/)) 141! 142 leximp=1 143 lend=3 144! 145 if(ithermal(1).eq.0) then 146 write(*,*)'*ERROR in umat_aniso_creep: no temperature defined;' 147 write(*,*) ' a creep calculation without temperature' 148 write(*,*) ' does not make sense' 149 write(*,*) 150 call exit(201) 151 endif 152! 153 iloop=0 154 exitcriterion=0 155! 156 c0=dsqrt(2.d0/3.d0) 157 c1=2.d0/3.d0 158 c2=-1.d0/3.d0 159! 160! elastic constants 161! 162 if(iorien.gt.0) then 163! 164 call transformatrix(orab(1,iorien),pgauss,skl) 165! 166 call orthotropic(elconloc,ya) 167! 168 do j=1,21 169 j1=kel(1,j) 170 j2=kel(2,j) 171 j3=kel(3,j) 172 j4=kel(4,j) 173 c(j)=0.d0 174 do j5=1,3 175 do j6=1,3 176 do j7=1,3 177 do j8=1,3 178 c(j)=c(j)+ya(j5,j6,j7,j8)* 179 & skl(j1,j5)*skl(j2,j6)*skl(j3,j7)*skl(j4,j8) 180 enddo 181 enddo 182 enddo 183 enddo 184 enddo 185! 186 else 187 do i=1,9 188 c(i)=elconloc(i) 189 enddo 190 endif 191! 192! state variables 193! 194! equivalent plastic strain 195! 196 epqini=xstateini(1,iint,iel) 197! 198! plastic strain 199! 200 do i=1,6 201 ep0(i)=xstateini(1+i,iint,iel) 202 enddo 203! elastic strains 204! 205 do i=1,6 206 ee(i)=emec(i)-ep0(i) 207 enddo 208! 209! global trial stress tensor 210! 211 if(iorien.gt.0) then 212 stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+ 213 & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6)) 214 & -beta(1) 215 stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+ 216 & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6)) 217 & -beta(2) 218 stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+ 219 & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6)) 220 & -beta(3) 221 stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+ 222 & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6)) 223 & -beta(4) 224 stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+ 225 & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6)) 226 & -beta(5) 227 stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+ 228 & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6)) 229 & -beta(6) 230 else 231 stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1) 232 stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1) 233 stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1) 234 stri(4)=2.d0*c(7)*ee(4)-beta(4) 235 stri(5)=2.d0*c(8)*ee(5)-beta(5) 236 stri(6)=2.d0*c(9)*ee(6)-beta(6) 237 endif 238! 239! stress radius (only deviatoric part of stress enters) 240! 241 strinv=(stri(1)+stri(2)+stri(3))/3.d0 242 do i=1,3 243 sg(i)=stri(i)-strinv 244 enddo 245 do i=4,6 246 sg(i)=stri(i) 247 enddo 248 dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+ 249 & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6))) 250! 251! evaluation of the yield surface 252! 253 ec(1)=epqini 254! 255 htri=dsg 256! 257! check whether plasticity occurs 258! 259 if(htri.gt.1.d-10) then 260 iplas=1 261 else 262 iplas=0 263 endif 264! 265 if((iplas.eq.0).or.(ielas.eq.1).or.(dtime.lt.1.d-30).or. 266 & ((nmethod.eq.1).and.(ithermal(1).ne.3))) then 267! 268! elastic stress 269! 270 do i=1,6 271 stre(i)=stri(i) 272 enddo 273! 274! updating the state variables 275! 276 xstate(1,iint,iel)=epqini 277 do i=1,6 278 xstate(1+i,iint,iel)=ep0(i) 279 enddo 280! 281! elastic stiffness 282! 283 if(icmd.ne.3) then 284 if(iorien.gt.0) then 285 do i=1,21 286 stiff(i)=c(i) 287 enddo 288 else 289 stiff(1)=c(1) 290 stiff(2)=c(2) 291 stiff(3)=c(3) 292 stiff(4)=c(4) 293 stiff(5)=c(5) 294 stiff(6)=c(6) 295 stiff(7)=0.d0 296 stiff(8)=0.d0 297 stiff(9)=0.d0 298 stiff(10)=c(7) 299 stiff(11)=0.d0 300 stiff(12)=0.d0 301 stiff(13)=0.d0 302 stiff(14)=0.d0 303 stiff(15)=c(8) 304 stiff(16)=0.d0 305 stiff(17)=0.d0 306 stiff(18)=0.d0 307 stiff(19)=0.d0 308 stiff(20)=0.d0 309 stiff(21)=c(9) 310 endif 311 endif 312! 313 return 314 endif 315! 316! plastic deformation 317! 318 neq=6 319 nrhs=1 320 lda=6 321 ldb=6 322! 323! initializing the state variables 324! 325 do i=1,6 326 ep(i)=ep0(i) 327 enddo 328 dg=0.d0 329! 330! determining the inverse of c 331! 332 if(iorien.gt.0) then 333! 334! solve gl:C=gr 335! 336 gl(1,1)=c(1) 337 gl(1,2)=c(2) 338 gl(2,2)=c(3) 339 gl(1,3)=c(4) 340 gl(2,3)=c(5) 341 gl(3,3)=c(6) 342 gl(1,4)=c(7) 343 gl(2,4)=c(8) 344 gl(3,4)=c(9) 345 gl(4,4)=c(10) 346 gl(1,5)=c(11) 347 gl(2,5)=c(12) 348 gl(3,5)=c(13) 349 gl(4,5)=c(14) 350 gl(5,5)=c(15) 351 gl(1,6)=c(16) 352 gl(2,6)=c(17) 353 gl(3,6)=c(18) 354 gl(4,6)=c(19) 355 gl(5,6)=c(20) 356 gl(6,6)=c(21) 357 do i=1,6 358 do j=1,i-1 359 gl(i,j)=gl(j,i) 360 enddo 361 enddo 362 do i=1,6 363 do j=1,6 364 gr(i,j)=0.d0 365 enddo 366 gr(i,i)=1.d0 367 enddo 368 nrhs=6 369 call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info) 370 if(info.ne.0) then 371 write(*,*) '*ERROR in sc.f: linear equation solver' 372 write(*,*) ' exited with error: info = ',info 373 call exit(201) 374 endif 375 nrhs=1 376 cm1(1)=gr(1,1) 377 cm1(2)=gr(1,2) 378 cm1(3)=gr(2,2) 379 cm1(4)=gr(1,3) 380 cm1(5)=gr(2,3) 381 cm1(6)=gr(3,3) 382 cm1(7)=gr(1,4)/2.d0 383 cm1(8)=gr(2,4)/2.d0 384 cm1(9)=gr(3,4)/2.d0 385 cm1(10)=gr(4,4)/4.d0 386 cm1(11)=gr(1,5)/2.d0 387 cm1(12)=gr(2,5)/2.d0 388 cm1(13)=gr(3,5)/2.d0 389 cm1(14)=gr(4,5)/4.d0 390 cm1(15)=gr(5,5)/4.d0 391 cm1(16)=gr(1,6)/2.d0 392 cm1(17)=gr(2,6)/2.d0 393 cm1(18)=gr(3,6)/2.d0 394 cm1(19)=gr(4,6)/4.d0 395 cm1(20)=gr(5,6)/4.d0 396 cm1(21)=gr(6,6)/4.d0 397 else 398 detc=c(1)*(c(3)*c(6)-c(5)*c(5))- 399 & c(2)*(c(2)*c(6)-c(4)*c(5))+ 400 & c(4)*(c(2)*c(5)-c(4)*c(3)) 401 cm1(1)=(c(3)*c(6)-c(5)*c(5))/detc 402 cm1(2)=(c(5)*c(4)-c(2)*c(6))/detc 403 cm1(3)=(c(1)*c(6)-c(4)*c(4))/detc 404 cm1(4)=(c(2)*c(5)-c(3)*c(4))/detc 405 cm1(5)=(c(2)*c(4)-c(1)*c(5))/detc 406 cm1(6)=(c(1)*c(3)-c(2)*c(2))/detc 407 cm1(7)=1.d0/(4.d0*c(7)) 408 cm1(8)=1.d0/(4.d0*c(8)) 409 cm1(9)=1.d0/(4.d0*c(9)) 410 endif 411! 412! first attempt: root search with Newton-Raphson 413! 414 loop: do 415! 416 iloop=iloop+1 417! 418! elastic strains 419! 420 do i=1,6 421 ee(i)=emec(i)-ep(i) 422 enddo 423! 424! global trial stress tensor 425! 426 if(iorien.gt.0) then 427 stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+ 428 & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6)) 429 & -beta(1) 430 stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+ 431 & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6)) 432 & -beta(2) 433 stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+ 434 & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6)) 435 & -beta(3) 436 stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+ 437 & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6)) 438 & -beta(4) 439 stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+ 440 & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6)) 441 & -beta(5) 442 stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+ 443 & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6)) 444 & -beta(6) 445 else 446 stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1) 447 stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1) 448 stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1) 449 stri(4)=2.d0*c(7)*ee(4)-beta(4) 450 stri(5)=2.d0*c(8)*ee(5)-beta(5) 451 stri(6)=2.d0*c(9)*ee(6)-beta(6) 452 endif 453! 454! stress radius (only deviatoric part of stress enters) 455! 456 strinv=(stri(1)+stri(2)+stri(3))/3.d0 457 do i=1,3 458 sg(i)=stri(i)-strinv 459 enddo 460 do i=4,6 461 sg(i)=stri(i) 462 enddo 463 dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+ 464 & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6))) 465! 466! evaluation of the yield surface 467! 468 ec(1)=epqini 469 decra(1)=c0*dg 470 timeabq(1)=time 471 timeabq(2)=ttime+time 472 call creep(decra,deswa,xstateini(1,iint,iel),serd,ec, 473 & esw,p,svm,t1l,dtemp,predef,dpred,timeabq,dtime, 474 & amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt, 475 & kstep,kinc) 476! 477! if the creep routine returns an increased value of decra(1) 478! it means that there is a lower cut-off for decra(1); 479! if the routine stays in a range lower than this cut-off, 480! it will never leave it and the exit conditions are 481! assumed to be satisfied. 482! 483 if(decra(1).gt.c0*dg) then 484 dg=decra(1)/c0 485 if(iloop.gt.1) exitcriterion=1 486 endif 487! 488 htri=dsg-c0*svm 489! 490 do i=1,6 491 sg(i)=sg(i)/dsg 492 enddo 493! 494! determining the residual matrix 495! 496 do i=1,6 497 r(i)=ep0(i)-ep(i)+dg*sg(i) 498 enddo 499! 500! check convergence 501! 502 if(exitcriterion.eq.1) exit 503 if((dabs(htri).le.1.d-3).and. 504 & ((iloop.gt.1).and.((dabs(ddg).lt.1.d-10).or. 505 & (dabs(ddg).lt.1.d-3*dabs(dg))))) then 506 dd=0.d0 507 do i=1,6 508 dd=dd+r(i)*r(i) 509 enddo 510 dd=sqrt(dd) 511 if(dd.le.1.d-10) then 512 exit 513 endif 514 endif 515! 516! determining b.x 517! 518 b=dg/dsg 519! 520 x(1)=b*(c1-sg(1)*sg(1)) 521 x(2)=b*(c2-sg(1)*sg(2)) 522 x(3)=b*(c1-sg(2)*sg(2)) 523 x(4)=b*(c2-sg(1)*sg(3)) 524 x(5)=b*(c2-sg(2)*sg(3)) 525 x(6)=b*(c1-sg(3)*sg(3)) 526 x(7)=-b*sg(1)*sg(4) 527 x(8)=-b*sg(2)*sg(4) 528 x(9)=-b*sg(3)*sg(4) 529 x(10)=b*(.5d0-sg(4)*sg(4)) 530 x(11)=-b*sg(1)*sg(5) 531 x(12)=-b*sg(2)*sg(5) 532 x(13)=-b*sg(3)*sg(5) 533 x(14)=-b*sg(4)*sg(5) 534 x(15)=b*(.5d0-sg(5)*sg(5)) 535 x(16)=-b*sg(1)*sg(6) 536 x(17)=-b*sg(2)*sg(6) 537 x(18)=-b*sg(3)*sg(6) 538 x(19)=-b*sg(4)*sg(6) 539 x(20)=-b*sg(5)*sg(6) 540 x(21)=b*(.5d0-sg(6)*sg(6)) 541! 542! filling the LHS 543! 544 if(iorien.gt.0) then 545 gl(1,1)=cm1(1)+x(1) 546 gl(1,2)=cm1(2)+x(2) 547 gl(2,2)=cm1(3)+x(3) 548 gl(1,3)=cm1(4)+x(4) 549 gl(2,3)=cm1(5)+x(5) 550 gl(3,3)=cm1(6)+x(6) 551 gl(1,4)=cm1(7)+x(7) 552 gl(2,4)=cm1(8)+x(8) 553 gl(3,4)=cm1(9)+x(9) 554 gl(4,4)=cm1(10)+x(10) 555 gl(1,5)=cm1(11)+x(11) 556 gl(2,5)=cm1(12)+x(12) 557 gl(3,5)=cm1(13)+x(13) 558 gl(4,5)=cm1(14)+x(14) 559 gl(5,5)=cm1(15)+x(15) 560 gl(1,6)=cm1(16)+x(16) 561 gl(2,6)=cm1(17)+x(17) 562 gl(3,6)=cm1(18)+x(18) 563 gl(4,6)=cm1(19)+x(19) 564 gl(5,6)=cm1(20)+x(20) 565 gl(6,6)=cm1(21)+x(21) 566 do i=1,6 567 do j=1,i-1 568 gl(i,j)=gl(j,i) 569 enddo 570 enddo 571 else 572 gl(1,1)=cm1(1)+x(1) 573 gl(1,2)=cm1(2)+x(2) 574 gl(2,2)=cm1(3)+x(3) 575 gl(1,3)=cm1(4)+x(4) 576 gl(2,3)=cm1(5)+x(5) 577 gl(3,3)=cm1(6)+x(6) 578 gl(1,4)=x(7) 579 gl(2,4)=x(8) 580 gl(3,4)=x(9) 581 gl(4,4)=cm1(7)+x(10) 582 gl(1,5)=x(11) 583 gl(2,5)=x(12) 584 gl(3,5)=x(13) 585 gl(4,5)=x(14) 586 gl(5,5)=cm1(8)+x(15) 587 gl(1,6)=x(16) 588 gl(2,6)=x(17) 589 gl(3,6)=x(18) 590 gl(4,6)=x(19) 591 gl(5,6)=x(20) 592 gl(6,6)=cm1(9)+x(21) 593 do i=1,6 594 do j=1,i-1 595 gl(i,j)=gl(j,i) 596 enddo 597 enddo 598 endif 599! 600! filling the RHS 601! 602 do i=1,6 603 gr(i,1)=sg(i) 604 enddo 605! 606! solve gl:(P:n)=gr 607! 608 call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info) 609 if(info.ne.0) then 610 write(*,*) '*ERROR in sc.f: linear equation solver' 611 write(*,*) ' exited with error: info = ',info 612 call exit(201) 613 endif 614! 615 do i=1,6 616 Pn(i)=gr(i,1) 617 enddo 618! 619! calculating the creep contribution 620! 621 gcreep=c1/decra(5) 622! 623! calculating the correction to the consistency parameter 624! 625 gm1=Pn(1)*sg(1)+Pn(2)*sg(2)+Pn(3)*sg(3)+ 626 & (Pn(4)*sg(4)+Pn(5)*sg(5)+Pn(6)*sg(6)) 627 gm1=1.d0/(gm1+gcreep) 628 ddg=gm1*(htri-(Pn(1)*r(1)+Pn(2)*r(2)+Pn(3)*r(3)+ 629 & (Pn(4)*r(4)+Pn(5)*r(5)+Pn(6)*r(6)))) 630! 631! updating the residual matrix 632! 633 do i=1,6 634 r(i)=r(i)+ddg*sg(i) 635 enddo 636! 637! update the plastic strain 638! 639 gr(1,1)=r(1) 640 gr(2,1)=r(2) 641 gr(3,1)=r(3) 642 gr(4,1)=r(4) 643 gr(5,1)=r(5) 644 gr(6,1)=r(6) 645! 646 call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info) 647 if(info.ne.0) then 648 write(*,*) '*ERROR in sc.f: linear equation solver' 649 write(*,*) ' exited with error: info = ',info 650 call exit(201) 651 endif 652! 653 if(iorien.gt.0) then 654 ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+cm1(4)*gr(3,1)+ 655 & (cm1(7)*gr(4,1)+cm1(11)*gr(5,1)+cm1(16)*gr(6,1)) 656 ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+cm1(5)*gr(3,1)+ 657 & (cm1(8)*gr(4,1)+cm1(12)*gr(5,1)+cm1(17)*gr(6,1)) 658 ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+cm1(6)*gr(3,1)+ 659 & (cm1(9)*gr(4,1)+cm1(13)*gr(5,1)+cm1(18)*gr(6,1)) 660 ep(4)=ep(4)+cm1(7)*gr(1,1)+cm1(8)*gr(2,1)+cm1(9)*gr(3,1)+ 661 & (cm1(10)*gr(4,1)+cm1(14)*gr(5,1)+cm1(19)*gr(6,1)) 662 ep(5)=ep(5)+cm1(11)*gr(1,1)+cm1(12)*gr(2,1)+cm1(13)*gr(3,1)+ 663 & (cm1(14)*gr(4,1)+cm1(15)*gr(5,1)+cm1(20)*gr(6,1)) 664 ep(6)=ep(6)+cm1(16)*gr(1,1)+cm1(17)*gr(2,1)+cm1(18)*gr(3,1)+ 665 & (cm1(19)*gr(4,1)+cm1(20)*gr(5,1)+cm1(21)*gr(6,1)) 666 else 667 ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+cm1(4)*gr(3,1) 668 ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+cm1(5)*gr(3,1) 669 ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+cm1(6)*gr(3,1) 670 ep(4)=ep(4)+cm1(7)*gr(4,1) 671 ep(5)=ep(5)+cm1(8)*gr(5,1) 672 ep(6)=ep(6)+cm1(9)*gr(6,1) 673 endif 674! 675! update the consistency parameter 676! 677 dg=dg+ddg 678! 679! end of major loop 680! 681 if((iloop.gt.15).or.(dg.le.0.d0)) then 682 iloop=1 683 dg=0.d0 684 do i=1,6 685 ep(i)=ep0(i) 686 enddo 687! 688! second attempt: root search through interval division 689! 690 do 691 if(iloop.gt.100) then 692c NOTE: write statements cause problems for 693c parallellized execution 694c write(*,*) 695c & '*WARNING in umat_aniso_creep: material loop' 696c write(*,*) ' did not converge in integration' 697c write(*,*) ' point',iint,'in element',iel,';' 698c write(*,*) ' the increment size is reduced' 699c write(*,*) 700 pnewdt=0.25d0 701 return 702 endif 703! 704! elastic strains 705! 706 do i=1,6 707 ee(i)=emec(i)-ep(i) 708 enddo 709! 710! global trial stress tensor 711! 712 if(iorien.gt.0) then 713 stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+ 714 & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6)) 715 & -beta(1) 716 stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+ 717 & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6)) 718 & -beta(2) 719 stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+ 720 & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6)) 721 & -beta(3) 722 stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+ 723 & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6)) 724 & -beta(4) 725 stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+ 726 & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6)) 727 & -beta(5) 728 stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+ 729 & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6)) 730 & -beta(6) 731 else 732 stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1) 733 stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1) 734 stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1) 735 stri(4)=2.d0*c(7)*ee(4)-beta(4) 736 stri(5)=2.d0*c(8)*ee(5)-beta(5) 737 stri(6)=2.d0*c(9)*ee(6)-beta(6) 738 endif 739! 740! stress radius (only deviatoric part of stress enters) 741! 742 strinv=(stri(1)+stri(2)+stri(3))/3.d0 743 do i=1,3 744 sg(i)=stri(i)-strinv 745 enddo 746 do i=4,6 747 sg(i)=stri(i) 748 enddo 749 dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+ 750 & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6))) 751! 752! evaluation of the yield surface 753! 754 ec(1)=epqini 755 decra(1)=c0*dg 756 timeabq(1)=time 757 timeabq(2)=ttime+time 758 call creep(decra,deswa,xstateini(1,iint,iel),serd,ec, 759 & esw,p,svm,t1l,dtemp,predef,dpred,timeabq,dtime, 760 & amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt, 761 & kstep,kinc) 762 if(decra(1).gt.c0*dg) then 763 dg=decra(1)/c0 764 if(abs(iloop).gt.2) exitcriterion=1 765 endif 766! 767! needed in case decra(1) was changed in subroutine creep, 768! for instance because it is too small 769! 770 dg=decra(1)/c0 771! 772 htri=dsg-c0*svm 773! 774 do i=1,6 775 sg(i)=sg(i)/dsg 776 enddo 777! 778! determining the residual matrix 779! 780 do i=1,6 781 r(i)=ep0(i)-ep(i)+dg*sg(i) 782 enddo 783! 784! check convergence 785! 786 if(exitcriterion.eq.1) exit loop 787 if((dabs(htri).le.1.d-3).and. 788 & ((iloop.gt.2).and.((dabs(ddg).lt.1.d-10).or. 789 & (dabs(ddg).lt.1.d-3*dabs(dg))))) then 790 dd=0.d0 791 do i=1,6 792 dd=dd+r(i)*r(i) 793 enddo 794 dd=sqrt(dd) 795 if(dd.le.1.d-10) then 796 exit loop 797 endif 798 endif 799 if(iloop.gt.100) then 800c write(*,*) 801c & '*ERROR: no convergence in umat_aniso_creep' 802c write(*,*) ' iloop>100' 803c write(*,*) 'htri,dd ',htri,dd 804 exit loop 805 endif 806! 807! determining b.x 808! 809 b=dg/dsg 810! 811 x(1)=b*(c1-sg(1)*sg(1)) 812 x(2)=b*(c2-sg(1)*sg(2)) 813 x(3)=b*(c1-sg(2)*sg(2)) 814 x(4)=b*(c2-sg(1)*sg(3)) 815 x(5)=b*(c2-sg(2)*sg(3)) 816 x(6)=b*(c1-sg(3)*sg(3)) 817 x(7)=-b*sg(1)*sg(4) 818 x(8)=-b*sg(2)*sg(4) 819 x(9)=-b*sg(3)*sg(4) 820 x(10)=b*(.5d0-sg(4)*sg(4)) 821 x(11)=-b*sg(1)*sg(5) 822 x(12)=-b*sg(2)*sg(5) 823 x(13)=-b*sg(3)*sg(5) 824 x(14)=-b*sg(4)*sg(5) 825 x(15)=b*(.5d0-sg(5)*sg(5)) 826 x(16)=-b*sg(1)*sg(6) 827 x(17)=-b*sg(2)*sg(6) 828 x(18)=-b*sg(3)*sg(6) 829 x(19)=-b*sg(4)*sg(6) 830 x(20)=-b*sg(5)*sg(6) 831 x(21)=b*(.5d0-sg(6)*sg(6)) 832! 833! filling the LHS 834! 835 if(iorien.gt.0) then 836 gl(1,1)=cm1(1)+x(1) 837 gl(1,2)=cm1(2)+x(2) 838 gl(2,2)=cm1(3)+x(3) 839 gl(1,3)=cm1(4)+x(4) 840 gl(2,3)=cm1(5)+x(5) 841 gl(3,3)=cm1(6)+x(6) 842 gl(1,4)=cm1(7)+x(7) 843 gl(2,4)=cm1(8)+x(8) 844 gl(3,4)=cm1(9)+x(9) 845 gl(4,4)=cm1(10)+x(10) 846 gl(1,5)=cm1(11)+x(11) 847 gl(2,5)=cm1(12)+x(12) 848 gl(3,5)=cm1(13)+x(13) 849 gl(4,5)=cm1(14)+x(14) 850 gl(5,5)=cm1(15)+x(15) 851 gl(1,6)=cm1(16)+x(16) 852 gl(2,6)=cm1(17)+x(17) 853 gl(3,6)=cm1(18)+x(18) 854 gl(4,6)=cm1(19)+x(19) 855 gl(5,6)=cm1(20)+x(20) 856 gl(6,6)=cm1(21)+x(21) 857 do i=1,6 858 do j=1,i-1 859 gl(i,j)=gl(j,i) 860 enddo 861 enddo 862 else 863 gl(1,1)=cm1(1)+x(1) 864 gl(1,2)=cm1(2)+x(2) 865 gl(2,2)=cm1(3)+x(3) 866 gl(1,3)=cm1(4)+x(4) 867 gl(2,3)=cm1(5)+x(5) 868 gl(3,3)=cm1(6)+x(6) 869 gl(1,4)=x(7) 870 gl(2,4)=x(8) 871 gl(3,4)=x(9) 872 gl(4,4)=cm1(7)+x(10) 873 gl(1,5)=x(11) 874 gl(2,5)=x(12) 875 gl(3,5)=x(13) 876 gl(4,5)=x(14) 877 gl(5,5)=cm1(8)+x(15) 878 gl(1,6)=x(16) 879 gl(2,6)=x(17) 880 gl(3,6)=x(18) 881 gl(4,6)=x(19) 882 gl(5,6)=x(20) 883 gl(6,6)=cm1(9)+x(21) 884 do i=1,6 885 do j=1,i-1 886 gl(i,j)=gl(j,i) 887 enddo 888 enddo 889 endif 890! 891! filling the RHS 892! 893 do i=1,6 894 gr(i,1)=sg(i) 895 enddo 896! 897! solve gl:(P:n)=gr 898! 899 call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info) 900 if(info.ne.0) then 901 write(*,*) '*ERROR in sc.f: linear equation solver' 902 write(*,*) ' exited with error: info = ',info 903 call exit(201) 904 endif 905! 906 do i=1,6 907 Pn(i)=gr(i,1) 908 enddo 909! 910! calculating the creep contribution 911! 912 gcreep=c1/decra(5) 913! 914! calculating the correction to the consistency parameter 915! 916 gm1=Pn(1)*sg(1)+Pn(2)*sg(2)+Pn(3)*sg(3)+ 917 & (Pn(4)*sg(4)+Pn(5)*sg(5)+Pn(6)*sg(6)) 918 gm1=1.d0/(gm1+gcreep) 919 fu=(htri-(Pn(1)*r(1)+Pn(2)*r(2)+Pn(3)*r(3)+ 920 & (Pn(4)*r(4)+Pn(5)*r(5)+Pn(6)*r(6)))) 921! 922 if(iloop.eq.1) then 923c write(*,*) 'iloop,dg,fu ',iloop,dg,fu 924 dg1=0.d0 925 fu1=fu 926 iloop=2 927 dg=1.d-10 928 ddg=dg 929 do i=1,6 930 ep1(i)=ep(i) 931 r1(i)=r(i) 932 sg1(i)=sg(i) 933 do j=1,6 934 gl1(i,j)=gl(i,j) 935 enddo 936 enddo 937 elseif((iloop.eq.2).or.(iloop.lt.0)) then 938 if(fu*fu1.lt.0.d0) then 939c write(*,*) 'iloop,dg,fu ',iloop,dg,fu 940 if(iloop.eq.2) then 941 iloop=3 942 else 943 iloop=-iloop+1 944 endif 945 fu2=fu 946 dg2=dg 947 dg=(dg1+dg2)/2.d0 948 ddg=(dg2-dg1)/2.d0 949 do i=1,6 950 ep(i)=ep1(i) 951 r(i)=r1(i) 952 sg(i)=sg1(i) 953 do j=1,6 954 gl(i,j)=gl1(i,j) 955 enddo 956 enddo 957 else 958c write(*,*) 'iloop,dg,fu ',iloop,dg,fu 959c dg1=dg 960c fu1=fu 961 if(iloop.eq.2) then 962 if(dabs(fu).gt.dabs(fu1)) exitcriterion=1 963 dg1=dg 964 fu1=fu 965 ddg=dg*9.d0 966 dg=dg*10.d0 967 else 968 dg1=dg 969 fu1=fu 970 dg=dg+ddg 971 iloop=iloop-1 972 endif 973 if(dg.gt.10.1d0) then 974 write(*,*) 975 & '*ERROR: no convergence in umat_aniso_creep' 976 write(*,*) ' dg>10.' 977 call exit(201) 978 endif 979 do i=1,6 980 ep1(i)=ep(i) 981 r1(i)=r(i) 982 sg1(i)=sg(i) 983 do j=1,6 984 gl1(i,j)=gl(i,j) 985 enddo 986 enddo 987 endif 988 else 989c write(*,*) 'iloop,dg,fu ',iloop,dg,fu 990 if(fu*fu1.ge.0.d0) then 991 dg1=dg 992 fu1=fu 993 dg=(dg1+dg2)/2.d0 994 ddg=(dg2-dg1)/2.d0 995 do i=1,6 996 ep1(i)=ep(i) 997 r1(i)=r(i) 998 sg1(i)=sg(i) 999 do j=1,6 1000 gl1(i,j)=gl(i,j) 1001 enddo 1002 enddo 1003 iloop=-iloop-1 1004 else 1005 dg2=dg 1006 fu2=fu 1007 dg=(dg1+dg2)/2.d0 1008 ddg=(dg2-dg1)/2.d0 1009 do i=1,6 1010 ep(i)=ep1(i) 1011 r(i)=r1(i) 1012 sg(i)=sg1(i) 1013 do j=1,6 1014 gl(i,j)=gl1(i,j) 1015 enddo 1016 enddo 1017 iloop=iloop+1 1018 endif 1019 endif 1020! 1021! updating the residual matrix 1022! 1023 do i=1,6 1024 r(i)=r(i)+ddg*sg(i) 1025 enddo 1026! 1027! update the plastic strain 1028! 1029 gr(1,1)=r(1) 1030 gr(2,1)=r(2) 1031 gr(3,1)=r(3) 1032 gr(4,1)=r(4) 1033 gr(5,1)=r(5) 1034 gr(6,1)=r(6) 1035! 1036 call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb, 1037 & info) 1038 if(info.ne.0) then 1039 write(*,*) '*ERROR in sc.f: linear equation solver' 1040 write(*,*) ' exited with error: info = ',info 1041 call exit(201) 1042 endif 1043! 1044 if(iorien.gt.0) then 1045 ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+ 1046 & cm1(4)*gr(3,1)+ 1047 & (cm1(7)*gr(4,1)+cm1(11)*gr(5,1)+ 1048 & cm1(16)*gr(6,1)) 1049 ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+ 1050 & cm1(5)*gr(3,1)+ 1051 & (cm1(8)*gr(4,1)+cm1(12)*gr(5,1)+ 1052 & cm1(17)*gr(6,1)) 1053 ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1) 1054 & +cm1(6)*gr(3,1)+ 1055 & (cm1(9)*gr(4,1)+cm1(13)*gr(5,1)+ 1056 & cm1(18)*gr(6,1)) 1057 ep(4)=ep(4)+cm1(7)*gr(1,1)+cm1(8)*gr(2,1)+ 1058 & cm1(9)*gr(3,1)+ 1059 & (cm1(10)*gr(4,1)+cm1(14)*gr(5,1)+ 1060 & cm1(19)*gr(6,1)) 1061 ep(5)=ep(5)+cm1(11)*gr(1,1)+cm1(12)*gr(2,1)+ 1062 & cm1(13)*gr(3,1)+ 1063 & (cm1(14)*gr(4,1)+cm1(15)*gr(5,1)+ 1064 & cm1(20)*gr(6,1)) 1065 ep(6)=ep(6)+cm1(16)*gr(1,1)+cm1(17)*gr(2,1)+ 1066 & cm1(18)*gr(3,1)+ 1067 & (cm1(19)*gr(4,1)+cm1(20)*gr(5,1)+ 1068 & cm1(21)*gr(6,1)) 1069 else 1070 ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+ 1071 & cm1(4)*gr(3,1) 1072 ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+ 1073 & cm1(5)*gr(3,1) 1074 ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+ 1075 & cm1(6)*gr(3,1) 1076 ep(4)=ep(4)+cm1(7)*gr(4,1) 1077 ep(5)=ep(5)+cm1(8)*gr(5,1) 1078 ep(6)=ep(6)+cm1(9)*gr(6,1) 1079 endif 1080! 1081! end of major loop 1082! 1083 enddo 1084! 1085 endif 1086! 1087 enddo loop 1088! 1089! storing the stress 1090! 1091 do i=1,6 1092 stre(i)=stri(i) 1093 enddo 1094! 1095! calculating the tangent stiffness matrix 1096! 1097 if(icmd.ne.3) then 1098! 1099! determining p 1100! 1101 gr(1,1)=1.d0 1102 gr(1,2)=0.d0 1103 gr(2,2)=1.d0 1104 gr(1,3)=0.d0 1105 gr(2,3)=0.d0 1106 gr(3,3)=1.d0 1107 gr(1,4)=0.d0 1108 gr(2,4)=0.d0 1109 gr(3,4)=0.d0 1110 gr(4,4)=1.d0 1111 gr(1,5)=0.d0 1112 gr(2,5)=0.d0 1113 gr(3,5)=0.d0 1114 gr(4,5)=0.d0 1115 gr(5,5)=1.d0 1116 gr(1,6)=0.d0 1117 gr(2,6)=0.d0 1118 gr(3,6)=0.d0 1119 gr(4,6)=0.d0 1120 gr(5,6)=0.d0 1121 gr(6,6)=1.d0 1122 do i=1,6 1123 do j=1,i-1 1124 gr(i,j)=gr(j,i) 1125 enddo 1126 enddo 1127 nrhs=6 1128! 1129 call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info) 1130 if(info.ne.0) then 1131 write(*,*) '*ERROR in sc.f: linear equation solver' 1132 write(*,*) ' exited with error: info = ',info 1133 call exit(201) 1134 endif 1135! 1136 stiff(1)=gr(1,1)-gm1*Pn(1)*Pn(1) 1137 stiff(2)=gr(1,2)-gm1*Pn(1)*Pn(2) 1138 stiff(3)=gr(2,2)-gm1*Pn(2)*Pn(2) 1139 stiff(4)=gr(1,3)-gm1*Pn(1)*Pn(3) 1140 stiff(5)=gr(2,3)-gm1*Pn(2)*Pn(3) 1141 stiff(6)=gr(3,3)-gm1*Pn(3)*Pn(3) 1142 stiff(7)=(gr(1,4)-gm1*Pn(1)*Pn(4))/2.d0 1143 stiff(8)=(gr(2,4)-gm1*Pn(2)*Pn(4))/2.d0 1144 stiff(9)=(gr(3,4)-gm1*Pn(3)*Pn(4))/2.d0 1145 stiff(10)=(gr(4,4)-gm1*Pn(4)*Pn(4))/4.d0 1146 stiff(11)=(gr(1,5)-gm1*Pn(1)*Pn(5))/2.d0 1147 stiff(12)=(gr(2,5)-gm1*Pn(2)*Pn(5))/2.d0 1148 stiff(13)=(gr(3,5)-gm1*Pn(3)*Pn(5))/2.d0 1149 stiff(14)=(gr(4,5)-gm1*Pn(4)*Pn(5))/4.d0 1150 stiff(15)=(gr(5,5)-gm1*Pn(5)*Pn(5))/4.d0 1151 stiff(16)=(gr(1,6)-gm1*Pn(1)*Pn(6))/2.d0 1152 stiff(17)=(gr(2,6)-gm1*Pn(2)*Pn(6))/2.d0 1153 stiff(18)=(gr(3,6)-gm1*Pn(3)*Pn(6))/2.d0 1154 stiff(19)=(gr(4,6)-gm1*Pn(4)*Pn(6))/4.d0 1155 stiff(20)=(gr(5,6)-gm1*Pn(5)*Pn(6))/4.d0 1156 stiff(21)=(gr(6,6)-gm1*Pn(6)*Pn(6))/4.d0 1157 endif 1158! 1159! updating the state variables 1160! 1161 xstate(1,iint,iel)=epqini+c0*dg 1162 do i=1,6 1163 xstate(1+i,iint,iel)=ep(i) 1164 enddo 1165! 1166! maximum equivalent viscoplastic strain in this increment 1167! 1168c depvisc=max(depvisc,c0*dg) 1169 depvisc=0.d0 1170! 1171 return 1172 end 1173