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_compression_only( 20 & amat,iel,iint,kode,elconloc,emec,emec0, 21 & beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime, 22 & icmd,ielas,mi,nstate_,xstateini,xstate,stre,stiff, 23 & iorien,pgauss,orab,pnewdt,ipkon) 24! 25! calculates stiffness and stresses for a user defined material 26! law 27! 28! icmd=3: calculates stress at mechanical strain 29! else: calculates stress at mechanical strain and the stiffness 30! matrix 31! 32! INPUT: 33! 34! amat material name 35! iel element number 36! iint integration point number 37! 38! kode material type (-100-#of constants entered 39! under *USER MATERIAL): can be used for materials 40! with varying number of constants 41! 42! elconloc(21) user defined constants defined by the keyword 43! card *USER MATERIAL (max. 21, actual # = 44! -kode-100), interpolated for the 45! actual temperature t1l 46! 47! emec(6) Lagrange mechanical strain tensor (component order: 48! 11,22,33,12,13,23) at the end of the increment 49! (thermal strains are subtracted) 50! emec0(6) Lagrange mechanical strain tensor at the start of the 51! increment (thermal strains are subtracted) 52! beta(6) residual stress tensor (the stress entered under 53! the keyword *INITIAL CONDITIONS,TYPE=STRESS) 54! 55! xokl(3,3) deformation gradient at the start of the increment 56! voj Jacobian at the start of the increment 57! xkl(3,3) deformation gradient at the end of the increment 58! vj Jacobian at the end of the increment 59! 60! ithermal 0: no thermal effects are taken into account 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! pnewdt to be specified by the user if the material 116! routine is unable to return the stiffness matrix 117! and/or the stress due to divergence within the 118! routine. pnewdt is the factor by which the time 119! increment is to be multiplied in the next 120! trial and should exceed zero but be less than 1. 121! Default is -1 indicating that the user routine 122! has converged. 123! ipkon(*) ipkon(iel) points towards the position in field 124! kon prior to the first node of the element's 125! topology. If ipkon(iel) is set to -1, the 126! element is removed from the mesh 127! 128 implicit none 129! 130 character*80 amat 131! 132 integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien, 133 & ipkon(*),kel(4,21),n,matz,ier,i,j,kal(2,6),j1,j2,j3,j4 134! 135 real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6), 136 & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), 137 & time,ttime,pnewdt,xstate(nstate_,mi(1),*),e(3,3),we(3), 138 & xstateini(nstate_,mi(1),*),wc(3),z(3,3),fv1(3),fv2(3),eps, 139 & pi,young,fla(3),xm1(3,3),xm2(3,3),xm3(3,3),dfla(3),a(21),b(21), 140 & d(3,3),c(3,3),xmm1(21),xmm2(21),xmm3(21),ca(3),cb(3) 141! 142 kal=reshape((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/)) 143! 144 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, 145 & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3, 146 & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3, 147 & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/)) 148! 149 d=reshape((/1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/), 150 & (/3,3/)) 151! 152! compression-only material. The constants are: 153! 1. Young's modulus 154! 2. maximum tension allowed (in absolute value) 155! 156 pi=4.d0*datan(1.d0) 157! 158 if(elconloc(1).lt.1.d-30) then 159 write(*,*) '*ERROR in umat_compression_only: the Young modulus' 160 write(*,*) ' is too small' 161 call exit(201) 162 endif 163! 164 if(elconloc(2).lt.1.d-30) then 165 write(*,*) '*ERROR in umat_compression_only: maximum tension' 166 write(*,*) ' value is too small' 167 call exit(201) 168 endif 169 young=elconloc(1) 170 eps=elconloc(2)*pi/young 171! 172! calculation of the eigenvalues and eigenvectors of the 173! Lagrange strain 174! 175 e(1,1)=emec(1) 176 e(2,2)=emec(2) 177 e(3,3)=emec(3) 178 e(1,2)=emec(4) 179 e(1,3)=emec(5) 180 e(2,3)=emec(6) 181 e(2,1)=emec(4) 182 e(3,1)=emec(5) 183 e(3,2)=emec(6) 184c if(iint.eq.1) write(*,*) 'emec' 185c if(iint.eq.1) write(*,100) (emec(i),i=1,6) 186! 187 n=3 188 matz=1 189 ier=0 190! 191 call rs(n,n,e,we,matz,z,fv1,fv2,ier) 192! 193 if(ier.ne.0) then 194 write(*,*) ' 195 & *ERROR calculating the eigenvalues/vectors in umat_compression' 196 call exit(201) 197 endif 198! 199! calculating the eigenvalues of the Cauchy tensor 200! at the end of the increment (eigenvalues of the Cauchy tensor) 201! 202 do i=1,3 203 wc(i)=2.d0*we(i)+1.d0 204c if(iint.eq.1) 205c & write(*,*) 'eigenvalues',i,we(i),z(1,i),z(2,i),z(3,i) 206 enddo 207! 208! check for 3 equal eigenvalues 209! 210 if((dabs(we(3)-we(2)).lt.1.d-10).and. 211 & (dabs(we(2)-we(1)).lt.1.d-10)) then 212 fla(1)=young*we(1)*(0.5d0+datan(-we(1)/eps)/pi) 213 do i=1,3 214 stre(i)=fla(1)-beta(i) 215 enddo 216 do i=4,6 217 stre(i)=0.d0-beta(i) 218 enddo 219c if(iint.eq.1) write(*,*) '1 diff stre' 220c if(iint.eq.1) write(*,100) (stre(i),i=1,6) 221! 222 if(icmd.ne.3) then 223 dfla(1)=young*((0.5d0+datan(-we(1)/eps)/pi)/2.d0 224 & -we(1)/(2.d0*pi*eps*(1.d0+(we(1)/eps)**2))) 225! 226 do j=1,21 227 j1=kel(1,j) 228 j2=kel(2,j) 229 j3=kel(3,j) 230 j4=kel(4,j) 231 stiff(j)=dfla(1)*(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1)) 232 enddo 233c if(iint.eq.1) write(*,*) '1 diff stiff' 234c if(iint.eq.1) write(*,100) (stiff(i),i=1,21) 235 endif 236! 237! 2 equal eigenvalues: w2=w3 238! 239 elseif(dabs(we(3)-we(2)).lt.1.d-10) then 240! 241! calculating the structural tensor for the unequal 242! eigenvalue of the Cauchy and the Lagrange tensor 243! 244 do i=1,3 245 do j=1,3 246 xm1(i,j)=z(i,1)*z(j,1) 247 enddo 248 enddo 249! 250! calculating the modified Young's modulus (due to the 251! compression-only modification) 252! 253 do i=1,2 254 fla(i)=young*we(i)*(0.5d0+datan(-we(i)/eps)/pi) 255 enddo 256! 257! calculating the stresses 258! 259 do j=1,6 260 j1=kal(1,j) 261 j2=kal(2,j) 262 stre(j)=fla(1)*xm1(j1,j2)+fla(2)*(d(j1,j2)-xm1(j1,j2)) 263 & -beta(j) 264 enddo 265c if(iint.eq.1) write(*,*) '2 diff stre (w2=w3)' 266c if(iint.eq.1) write(*,100) (stre(i),i=1,6) 267! 268 if(icmd.ne.3) then 269! 270! calculating the stiffness matrix 271! 272! derivative of the modified young's modulus w.r.t. 273! the Cauchy eigenvalues 274! 275 do i=1,2 276 dfla(i)=young*((0.5d0+datan(-we(i)/eps)/pi)/2.d0 277 & -we(i)/(2.d0*pi*eps*(1.d0+(we(i)/eps)**2))) 278 enddo 279! 280! stiffness matrix 281! 282 do j=1,21 283 j1=kel(1,j) 284 j2=kel(2,j) 285 j3=kel(3,j) 286 j4=kel(4,j) 287! 288! calculating the auxiliary field xmm1 289! 290 xmm1(j)=xm1(j1,j2)*xm1(j3,j4) 291! 292! calculating the auxiliary fields a and b 293! 294 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 295 b(j)=(d(j3,j1)*xm1(j4,j2)+d(j4,j1)*xm1(j3,j2)+ 296 & d(j4,j2)*xm1(j1,j3)+d(j3,j2)*xm1(j1,j4))/2.d0 297! 298! calculating the stiffness matrix 299! 300 stiff(j)=2.d0*(dfla(1)*xmm1(j) 301 & +dfla(2)*(a(j)+xmm1(j)-b(j)) 302 & +(fla(1)-fla(2))*(b(j)-2.d0*xmm1(j)) 303 & /(2.d0*(we(1)-we(2)))) 304 enddo 305c if(iint.eq.1) write(*,*) '2 diff stiff (w2=w3)' 306c if(iint.eq.1) write(*,100) (stiff(i),i=1,21) 307 endif 308! 309! 2 equal eigenvalues: w1=w2 310! 311 elseif(dabs(we(2)-we(1)).lt.1.d-10) then 312! 313! calculating the structural tensor for the unequal 314! eigenvalue of the Cauchy and the Lagrange tensor 315! 316 do i=1,3 317 do j=1,3 318 xm3(i,j)=z(i,3)*z(j,3) 319 enddo 320 enddo 321! 322! calculating the modified Young's modulus (due to the 323! compression-only modification) 324! 325 do i=2,3 326 fla(i)=young*we(i)*(0.5d0+datan(-we(i)/eps)/pi) 327 enddo 328! 329! calculating the stresses 330! 331 do j=1,6 332 j1=kal(1,j) 333 j2=kal(2,j) 334 stre(j)=fla(3)*xm3(j1,j2)+fla(2)*(d(j1,j2)-xm3(j1,j2)) 335 & -beta(j) 336 enddo 337c if(iint.eq.1) write(*,*) '2 diff stre (w1=w2)' 338c if(iint.eq.1) write(*,100) (stre(i),i=1,6) 339! 340 if(icmd.ne.3) then 341! 342! calculating the stiffness matrix 343! 344! derivative of the modified young's modulus w.r.t. 345! the Cauchy eigenvalues 346! 347 do i=2,3 348 dfla(i)=young*((0.5d0+datan(-we(i)/eps)/pi)/2.d0 349 & -we(i)/(2.d0*pi*eps*(1.d0+(we(i)/eps)**2))) 350 enddo 351! 352! stiffness matrix 353! 354 do j=1,21 355 j1=kel(1,j) 356 j2=kel(2,j) 357 j3=kel(3,j) 358 j4=kel(4,j) 359! 360! calculating the auxiliary field xmm3 361! 362 xmm3(j)=xm3(j1,j2)*xm3(j3,j4) 363! 364! calculating the auxiliary fields a and b 365! 366 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 367 b(j)=(d(j3,j1)*xm3(j4,j2)+d(j4,j1)*xm3(j3,j2)+ 368 & d(j4,j2)*xm3(j1,j3)+d(j3,j2)*xm3(j1,j4))/2.d0 369! 370! calculating the stiffness matrix 371! 372 stiff(j)=2.d0*(dfla(3)*xmm3(j) 373 & +dfla(2)*(a(j)+xmm3(j)-b(j)) 374 & +(fla(2)-fla(3))*(b(j)-2.d0*xmm3(j)) 375 & /(2.d0*(we(2)-we(3)))) 376 enddo 377c if(iint.eq.1) write(*,*) '2 diff stiff (w1=w2)' 378c if(iint.eq.1) write(*,100) (stiff(i),i=1,21) 379 endif 380! 381! 2 equal eigenvalues: w1=w3 382! 383 elseif(dabs(we(3)-we(1)).lt.1.d-10) then 384! 385! calculating the structural tensor for the unequal 386! eigenvalue of the Cauchy and the Lagrange tensor 387! 388 do i=1,3 389 do j=1,3 390 xm3(i,j)=z(i,3)*z(j,3) 391 enddo 392 enddo 393! 394! calculating the modified Young's modulus (due to the 395! compression-only modification) 396! 397 do i=2,3 398 fla(i)=young*we(i)*(0.5d0+datan(-we(i)/eps)/pi) 399 enddo 400! 401! calculating the stresses 402! 403 do j=1,6 404 j1=kal(1,j) 405 j2=kal(2,j) 406 stre(j)=fla(2)*xm3(j1,j2)+fla(3)*(d(j1,j2)-xm3(j1,j2)) 407 & -beta(j) 408 enddo 409c if(iint.eq.1) write(*,*) '2 diff stre (w1=w2)' 410c if(iint.eq.1) write(*,100) (stre(i),i=1,6) 411! 412 if(icmd.ne.3) then 413! 414! calculating the stiffness matrix 415! 416! derivative of the modified young's modulus w.r.t. 417! the Cauchy eigenvalues 418! 419 do i=2,3 420 dfla(i)=young*((0.5d0+datan(-we(i)/eps)/pi)/2.d0 421 & -we(i)/(2.d0*pi*eps*(1.d0+(we(i)/eps)**2))) 422 enddo 423! 424! stiffness matrix 425! 426 do j=1,21 427 j1=kel(1,j) 428 j2=kel(2,j) 429 j3=kel(3,j) 430 j4=kel(4,j) 431! 432! calculating the auxiliary field xmm3 433! 434 xmm3(j)=xm3(j1,j2)*xm3(j3,j4) 435! 436! calculating the auxiliary fields a and b 437! 438 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 439 b(j)=(d(j3,j1)*xm3(j4,j2)+d(j4,j1)*xm3(j3,j2)+ 440 & d(j4,j2)*xm3(j1,j3)+d(j3,j2)*xm3(j1,j4))/2.d0 441! 442! calculating the stiffness matrix 443! 444 stiff(j)=2.d0*(dfla(2)*xmm3(j) 445 & +dfla(3)*(a(j)+xmm3(j)-b(j)) 446 & +(fla(3)-fla(2))*(b(j)-2.d0*xmm3(j)) 447 & /(2.d0*(we(3)-we(2)))) 448 enddo 449c if(iint.eq.1) write(*,*) '2 diff stiff (w1=w2)' 450c if(iint.eq.1) write(*,100) (stiff(i),i=1,21) 451 endif 452! 453! 3 different eigenvalues 454! 455 else 456! 457! calculating the structural tensors of the Cauchy and the 458! Lagrange tensor 459! 460 do i=1,3 461 do j=1,3 462 xm1(i,j)=z(i,1)*z(j,1) 463 xm2(i,j)=z(i,2)*z(j,2) 464 xm3(i,j)=z(i,3)*z(j,3) 465 enddo 466 enddo 467! 468! calculating the modified Young's modulus (due to the 469! compression-only modification) 470! 471 do i=1,3 472 fla(i)=young*we(i)*(0.5d0+datan(-we(i)/eps)/pi) 473 enddo 474! 475! calculating the stresses 476! 477 do j=1,6 478 j1=kal(1,j) 479 j2=kal(2,j) 480 stre(j)=fla(1)*xm1(j1,j2)+fla(2)*xm2(j1,j2) 481 & +fla(3)*xm3(j1,j2)-beta(j) 482 enddo 483c if(iint.eq.1) write(*,*) '3 diff stre' 484c if(iint.eq.1) write(*,100) (stre(i),i=1,6) 485! 486 if(icmd.ne.3) then 487! 488! calculating the stiffness matrix 489! 490! derivative of the modified young's modulus w.r.t. 491! the Cauchy eigenvalues 492! 493 do i=1,3 494 dfla(i)=young*((0.5d0+datan(-we(i)/eps)/pi)/2.d0 495 & -we(i)/(2.d0*pi*eps*(1.d0+(we(i)/eps)**2))) 496 enddo 497! 498 cb(1)=1.d0/(4.d0*(we(1)-we(2))*(we(1)-we(3))) 499 ca(1)=-(wc(2)+wc(3))*cb(1) 500 cb(2)=1.d0/(4.d0*(we(2)-we(3))*(we(2)-we(1))) 501 ca(2)=-(wc(3)+wc(1))*cb(2) 502 cb(3)=1.d0/(4.d0*(we(3)-we(1))*(we(3)-we(2))) 503 ca(3)=-(wc(1)+wc(2))*cb(3) 504! 505! Cauchy tensor 506! 507 do i=1,3 508 do j=1,3 509 c(i,j)=2.d0*e(i,j) 510 enddo 511 c(i,i)=c(i,i)+1.d0 512 enddo 513! 514! stiffness matrix 515! 516 do j=1,21 517 j1=kel(1,j) 518 j2=kel(2,j) 519 j3=kel(3,j) 520 j4=kel(4,j) 521! 522! calculating the auxiliary fields xmm1,xmm2,xmm3 523! 524 xmm1(j)=xm1(j1,j2)*xm1(j3,j4) 525 xmm2(j)=xm2(j1,j2)*xm2(j3,j4) 526 xmm3(j)=xm3(j1,j2)*xm3(j3,j4) 527! 528! calculating the auxiliary fields a and b 529! (Dhondt, G., The Finite Element Method for Three- 530! Dimensional Thermomechanical Applications, Wiley 2004, 531! formulae (4.203) and (4.204)) 532! 533 a(j)=(d(j3,j1)*d(j4,j2)+d(j3,j2)*d(j4,j1))/2.d0 534 & -xmm1(j)-xmm2(j)-xmm3(j) 535 b(j)=(d(j3,j1)*c(j4,j2)+d(j4,j1)*c(j3,j2)+ 536 & d(j4,j2)*c(j1,j3)+d(j3,j2)*c(j1,j4))/2.d0 537 & -2.d0*(wc(1)*xmm1(j)+wc(2)*xmm2(j)+wc(3)*xmm3(j)) 538! 539! calculating the stiffness matrix 540! 541 stiff(j)=2.d0*(dfla(1)*xmm1(j)+dfla(2)*xmm2(j)+ 542 & dfla(3)*xmm3(j) 543 & +fla(1)*(cb(1)*b(j)+ca(1)*a(j)) 544 & +fla(2)*(cb(2)*b(j)+ca(2)*a(j)) 545 & +fla(3)*(cb(3)*b(j)+ca(3)*a(j))) 546 enddo 547c if(iint.eq.1) write(*,*) '3 diff stiff' 548c if(iint.eq.1) write(*,100) (stiff(i),i=1,21) 549 endif 550 endif 551! 552c 100 format(6(1x,e11.4)) 553! 554 return 555 end 556 557 558