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_undo_nlgeom_lin_el(amat,iel,iint,kode, 20 & elconloc,emec,emec0,beta,xokl,voj,xkl,vj,ithermal,t1l, 21 & dtime,time,ttime,icmd,ielas,mi,nstate_,xstateini,xstate, 22 & stre,stiff,iorien,pgauss,orab,eloc,nlgeom_undo) 23! 24! calculates stiffness and stresses for a linear elastic isotropic 25! material with special modification of the strain tensor 26! 27! The strain tensor that enters the routine is the Lagrange strain 28! tensor. This tensor is modified into an infinitesimal strain tensor 29! for large rotations. Application: singular stress/strain field at a 30! crack tip 31! 32! Procedure: The Biot strain U-I is 33! used to calculate the stresses by a 34! linear elastic isotropic relationship 35! 36! 37! icmd=3: calculates stress at mechanical strain 38! else: calculates stress at mechanical strain and the stiffness 39! matrix 40! 41! INPUT: 42! 43! amat material name 44! iel element number 45! iint integration point number 46! 47! kode material type (-100-#of constants entered 48! under *USER MATERIAL): can be used for materials 49! with varying number of constants 50! 51! elconloc(21) user defined constants defined by the keyword 52! card *USER MATERIAL (max. 21, actual # = 53! -kode-100), interpolated for the 54! actual temperature t1l 55! 56! emec(6) Lagrange mechanical strain tensor (component order: 57! 11,22,33,12,13,23) at the end of the increment 58! (thermal strains are subtracted) 59! emec0(6) Lagrange mechanical strain tensor at the start of the 60! increment (thermal strains are subtracted) 61! beta(6) residual stress tensor (the stress entered under 62! the keyword *INITIAL CONDITIONS,TYPE=STRESS) 63! 64! xokl(3,3) deformation gradient at the start of the increment 65! voj Jacobian at the start of the increment 66! xkl(3,3) deformation gradient at the end of the increment 67! vj Jacobian at the end of the increment 68! 69! ithermal 0: no thermal effects are taken into account 70! >0: thermal effects are taken into account (triggered 71! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) 72! t1l temperature at the end of the increment 73! dtime time length of the increment 74! time step time at the end of the current increment 75! ttime total time at the start of the current step 76! 77! icmd not equal to 3: calculate stress and stiffness 78! 3: calculate only stress 79! ielas 0: no elastic iteration: irreversible effects 80! are allowed 81! 1: elastic iteration, i.e. no irreversible 82! deformation allowed 83! 84! mi(1) max. # of integration points per element in the 85! model 86! nstate_ max. # of state variables in the model 87! 88! xstateini(nstate_,mi(1),# of elements) 89! state variables at the start of the increment 90! xstate(nstate_,mi(1),# of elements) 91! state variables at the end of the increment 92! 93! stre(6) Piola-Kirchhoff stress of the second kind 94! at the start of the increment 95! 96! iorien number of the local coordinate axis system 97! in the integration point at stake (takes the value 98! 0 if no local system applies) 99! pgauss(3) global coordinates of the integration point 100! orab(7,*) description of all local coordinate systems. 101! If a local coordinate system applies the global 102! tensors can be obtained by premultiplying the local 103! tensors with skl(3,3). skl is determined by calling 104! the subroutine transformatrix: 105! call transformatrix(orab(1,iorien),pgauss,skl) 106! eloc(6) Lagrange total strain tensor (component order: 107! 11,22,33,12,13,23) at the end of the increment 108! (thermal strains are subtracted) 109! 110! 111! OUTPUT: 112! 113! xstate(nstate_,mi(1),# of elements) 114! updated state variables at the end of the increment 115! stre(6) Piola-Kirchhoff stress of the second kind at the 116! end of the increment 117! stiff(21): consistent tangent stiffness matrix in the material 118! frame of reference at the end of the increment. In 119! other words: the derivative of the PK2 stress with 120! respect to the Lagrangian strain tensor. The matrix 121! is supposed to be symmetric, only the upper half is 122! to be given in the same order as for a fully 123! anisotropic elastic material (*ELASTIC,TYPE=ANISO). 124! Notice that the matrix is an integral part of the 125! fourth order material tensor, i.e. the Voigt notation 126! is not used. 127! emec(6) linear mechanical strain tensor for large 128! rotations (component order: 129! 11,22,33,12,13,23) at the end of the increment 130! (thermal strains are subtracted) 131! eloc(6) linear total strain tensor for large 132! rotations (component order: 133! 11,22,33,12,13,23) at the end of the increment 134! nlgeom_undo 0: Lagrange strain goes out 135! 1: linear strain for large rotations goes out 136! 137 implicit none 138! 139 character*80 amat 140! 141 integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien, 142 & i,j,k,nlgeom_undo,n,matz,ier,nconstants,mattyp 143! 144 real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6), 145 & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), 146 & time,ttime,ukl(3,3),ca(3,3),elag(3,3),elin(3,3), 147 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*), 148 & e,un,al,um,am1,am2,eme1(6),eloc(6),w(3),z(3,3),fv1(3), 149 & fv2(3),v1,v2,v3,c(6),c2(6),dd,d(6),u(6),ur(3,3) 150! 151 d=(/1.d0,1.d0,1.d0,0.d0,0.d0,0.d0/) 152! 153 nlgeom_undo=1 154! 155! a) calculate the Cauchy-Green tensor C from the mechanical Lagrange strain 156! 157 do i=1,6 158 c(i)=2*emec(i) 159 enddo 160 do i=1,3 161 c(i)=c(i)+1.d0 162 enddo 163! 164! b) storing C as matrix 165! 166 ca(1,1)=c(1) 167 ca(2,2)=c(2) 168 ca(3,3)=c(3) 169 ca(1,2)=c(4) 170 ca(1,3)=c(5) 171 ca(2,3)=c(6) 172 ca(2,1)=c(4) 173 ca(3,1)=c(5) 174 ca(3,2)=c(6) 175! 176! c) calculating the principal stretches 177! 178 n=3 179 matz=1 180! 181 call rs(n,n,ca,w,matz,z,fv1,fv2,ier) 182! 183 if(ier.ne.0) then 184 write(*,*) ' 185 & *ERROR calculating the eigenvalues/vectors in ' 186 write(*,*) ' umat_undo_nlgeom_lin_el' 187 call exit(201) 188 endif 189 w(1)=dsqrt(w(1)) 190 w(2)=dsqrt(w(2)) 191 w(3)=dsqrt(w(3)) 192! 193! d) calculating the invariants of U 194! 195 v1=w(1)+w(2)+w(3) 196 v2=w(1)*w(2)+w(2)*w(3)+w(3)*w(1) 197 v3=w(1)*w(2)*w(3) 198! 199! e) calculating C.C 200! 201 c2(1)=c(1)*c(1)+c(4)*c(4)+c(5)*c(5) 202 c2(2)=c(4)*c(4)+c(2)*c(2)+c(6)*c(6) 203 c2(3)=c(5)*c(5)+c(6)*c(6)+c(3)*c(3) 204 c2(4)=c(1)*c(4)+c(4)*c(2)+c(5)*c(6) 205 c2(5)=c(1)*c(5)+c(4)*c(6)+c(5)*c(3) 206 c2(6)=c(4)*c(5)+c(2)*c(6)+c(6)*c(3) 207! 208! f) calculating the right stretch tensor U 209! (cf. Simo and Hughes, Computational Inelasticity) 210! 211 dd=v1*v2-v3 212 do i=1,6 213 u(i)=(-c2(i)+(v1*v1-v2)*c(i)+v1*v3*d(i))/dd 214 enddo 215! 216! g) calculating the Biot strain = U-I 217! 218 do i=1,3 219 emec(i)=u(i)-1.d0 220 enddo 221 do i=4,6 222 emec(i)=u(i) 223 enddo 224! 225 nconstants=-kode-100 226! 227! calculating the stress and the linear elastic material data 228! 229 call linel(nconstants,mattyp,beta,emec,stre,stiff,elconloc, 230 & iorien,orab,pgauss) 231! 232c do i=1,6 233c write(*,*) 'umat...lin_el',time,iel,iint,elin(i),stre(i) 234c enddo 235! 236! for a user material (umat) the material is considered to be 237! fully anisotropic; rewriting the isotropic or orthotropic 238! data in an anisotropic format 239! 240 if(icmd.ne.3) then 241 if(mattyp.eq.1) then 242! 243! isotropic 244! 245 e=elconloc(1) 246 un=elconloc(2) 247 um=e/(1.d0+un) 248 al=un*um/(1.d0-2.d0*un) 249 um=um/2.d0 250! 251 stiff(1)=al+2.d0*um 252 stiff(2)=al 253 stiff(3)=al+2.d0*um 254 stiff(4)=al 255 stiff(5)=al 256 stiff(6)=al+2.d0*um 257 stiff(7)=0.d0 258 stiff(8)=0.d0 259 stiff(9)=0.d0 260 stiff(10)=um 261 stiff(11)=0.d0 262 stiff(12)=0.d0 263 stiff(13)=0.d0 264 stiff(14)=0.d0 265 stiff(15)=um 266 stiff(16)=0.d0 267 stiff(17)=0.d0 268 stiff(18)=0.d0 269 stiff(19)=0.d0 270 stiff(20)=0.d0 271 stiff(21)=um 272 elseif(mattyp.eq.2) then 273! 274! orthotropic 275! 276 stiff(10)=stiff(7) 277 stiff(15)=stiff(8) 278 stiff(21)=stiff(9) 279 do i=7,9 280 stiff(i)=0.d0 281 enddo 282 do i=11,14 283 stiff(i)=0.d0 284 enddo 285 do i=16,20 286 stiff(i)=0.d0 287 enddo 288 endif 289 endif 290! 291! calculating the stresses 292! 293c e=elconloc(1) 294c un=elconloc(2) 295c al=un*e/(1.d0+un)/(1.d0-2.d0*un) 296c um=e/2.d0/(1.d0+un) 297c am1=al+2.d0*um 298c am2=2.d0*um 299c! 300c stre(1)=am1*emec(1)+al*(emec(2)+emec(3))-beta(1) 301c stre(2)=am1*emec(2)+al*(emec(1)+emec(3))-beta(2) 302c stre(3)=am1*emec(3)+al*(emec(1)+emec(2))-beta(3) 303c stre(4)=am2*emec(4)-beta(4) 304c stre(5)=am2*emec(5)-beta(5) 305c stre(6)=am2*emec(6)-beta(6) 306cc write(*,*) 'umat_undo..stre',stre(1) 307c! 308cc do i=1,6 309cc write(*,*) 'umat...lin_el',time,iel,iint,eloc(i),stre(i) 310cc enddo 311c! 312c if(icmd.ne.3) then 313c! 314c! calculating the stiffness matrix 315c! 316c stiff(1)=al+2.d0*um 317c stiff(2)=al 318c stiff(3)=al+2.d0*um 319c stiff(4)=al 320c stiff(5)=al 321c stiff(6)=al+2.d0*um 322c stiff(7)=0.d0 323c stiff(8)=0.d0 324c stiff(9)=0.d0 325c stiff(10)=um 326c stiff(11)=0.d0 327c stiff(12)=0.d0 328c stiff(13)=0.d0 329c stiff(14)=0.d0 330c stiff(15)=um 331c stiff(16)=0.d0 332c stiff(17)=0.d0 333c stiff(18)=0.d0 334c stiff(19)=0.d0 335c stiff(20)=0.d0 336c stiff(21)=um 337c endif 338! 339 return 340 end 341