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! Author: Sven Kassbohm 20! Permission to use this routine was given by the author on 21! January 17, 2015. 22! 23 subroutine umat_ciarlet_el(amat,iel,iint,kode,elconloc,emec,emec0, 24 & beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime, 25 & icmd,ielas,mi, 26 & nstate_,xstateini,xstate,stre,stiff,iorien,pgauss,orab) 27! 28! calculates stiffness and stresses for a user defined material 29! law 30! 31! icmd=3: calculates stress at mechanical strain 32! else: calculates stress at mechanical strain and the stiffness 33! matrix 34! 35! INPUT: 36! 37! amat material name 38! iel element number 39! iint integration point number 40! 41! kode material type (-100-#of constants entered 42! under *USER MATERIAL): can be used for materials 43! with varying number of constants 44! 45! elconloc(21) user defined constants defined by the keyword 46! card *USER MATERIAL (max. 21, actual # = 47! -kode-100), interpolated for the 48! actual temperature t1l 49! 50! emec(6) Lagrange mechanical strain tensor (component order: 51! 11,22,33,12,13,23) at the end of the increment 52! (thermal strains are subtracted) 53! emec0(6) Lagrange mechanical strain tensor at the start of the 54! increment (thermal strains are subtracted) 55! beta(6) residual stress tensor (the stress entered under 56! the keyword *INITIAL CONDITIONS,TYPE=STRESS) 57! 58! xokl(3,3) deformation gradient at the start of the increment 59! voj Jacobian at the start of the increment 60! xkl(3,3) deformation gradient at the end of the increment 61! vj Jacobian at the end of the increment 62! 63! ithermal 0: no thermal effects are taken into account 64! >0: thermal effects are taken into account (triggered 65! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) 66! t1l temperature at the end of the increment 67! dtime time length of the increment 68! time step time at the end of the current increment 69! ttime total time at the start of the current increment 70! 71! icmd not equal to 3: calculate stress and stiffness 72! 3: calculate only stress 73! ielas 0: no elastic iteration: irreversible effects 74! are allowed 75! 1: elastic iteration, i.e. no irreversible 76! deformation allowed 77! 78! mi(1) max. # of integration points per element in the 79! model 80! nstate_ max. # of state variables in the model 81! 82! xstateini(nstate_,mi(1),# of elements) 83! state variables at the start of the increment 84! xstate(nstate_,mi(1),# of elements) 85! state variables at the end of the increment 86! 87! stre(6) Piola-Kirchhoff stress of the second kind 88! at the start of the increment 89! 90! iorien number of the local coordinate axis system 91! in the integration point at stake (takes the value 92! 0 if no local system applies) 93! pgauss(3) global coordinates of the integration point 94! orab(7,*) description of all local coordinate systems. 95! If a local coordinate system applies the global 96! tensors can be obtained by premultiplying the local 97! tensors with skl(3,3). skl is determined by calling 98! the subroutine transformatrix: 99! call transformatrix(orab(1,iorien),pgauss,skl) 100! 101! 102! OUTPUT: 103! 104! xstate(nstate_,mi(1),# of elements) 105! updated state variables at the end of the increment 106! stre(6) Piola-Kirchhoff stress of the second kind at the 107! end of the increment 108! stiff(21): consistent tangent stiffness matrix in the material 109! frame of reference at the end of the increment. In 110! other words: the derivative of the PK2 stress with 111! respect to the Lagrangian strain tensor. The matrix 112! is supposed to be symmetric, only the upper half is 113! to be given in the same order as for a fully 114! anisotropic elastic material (*ELASTIC,TYPE=ANISO). 115! Notice that the matrix is an integral part of the 116! fourth order material tensor, i.e. the Voigt notation 117! is not used. 118! 119 implicit none 120! 121 character*80 amat 122 integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien 123 real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6), 124 & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), 125 & time,ttime 126 real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*) 127! 128! Ela-Pla-modifications: 129! 130 real*8 C1(3,3), C1i(3,3), detC1, dS2dE(3,3,3,3), S2PK1(3,3) 131 real*8 lamda, mu,un,e 132! 133! change on 190315: input are Young's modulus and Poisson's 134! coefficient instead of the Lame constants 135! 136 e=elconloc(1) 137 un=elconloc(2) 138 mu=e/(1.d0+un) 139 lamda=mu*un/(1.d0-2.d0*un) 140 mu=mu/2.d0 141 142c lamda=elconloc(1) 143c mu=elconloc(2) 144 145C ciarlet_doc 146C 147C Compute pullback C1(3,3) of metric tensor g. 148c C1 = matmul(transpose(xkl),xkl) 149 150 151! 152! change on 190315: the Cauchy tensor should be 153! computed based on the mechanical Lagrange strain in 154! order to discard the thermal strains; the deformation 155! gradient still contains the thermal strain effect 156! 157! right Cauchy-Green tensor (eloc contains the Lagrange strain, 158! including thermal strain) 159! 160 C1(1,1)=2.d0*emec(1)+1.d0 161 C1(2,2)=2.d0*emec(2)+1.d0 162 C1(3,3)=2.d0*emec(3)+1.d0 163 C1(1,2)=2.d0*emec(4) 164 C1(2,1)=C1(1,2) 165 C1(1,3)=2.d0*emec(5) 166 C1(3,1)=C1(1,3) 167 C1(2,3)=2.d0*emec(6) 168 C1(3,2)=C1(2,3) 169C 170C 2PK-Stress at C: 171 call S2_Ciarlet(C1,lamda,mu, S2PK1,C1i,detC1) 172C Tangent/linearization: 173 call dS2_Ciarlet_dE(lamda,mu,detC1,C1i,dS2dE) 174 175C Dhondt-Voigt-Order, see umat_lin_iso_el.f 176 call voigt_33mat_to_6vec(S2PK1,stre) 177 call voigt_tetrad_to_matrix(dS2dE,stiff) 178 179C ciarlet_cod 180 181 return 182 end 183 184 185 186 subroutine voigt_33mat_to_6vec(mat,vec) 187 implicit none 188C in: 189 double precision mat(3,3) 190C out: 191 double precision vec(6) 192C 11,22,33,12,13,23 193 vec(1)=mat(1,1) 194 vec(2)=mat(2,2) 195 vec(3)=mat(3,3) 196 vec(4)=mat(1,2) 197 vec(5)=mat(1,3) 198 vec(6)=mat(2,3) 199 200 return 201 end 202 203 204 subroutine voigt_tetrad_to_matrix(C,stiff) 205 implicit none 206C in: 207 double precision C(3,3,3,3) 208C out: 209 double precision stiff(21) 210 211C indices according to 212C 1) file anisotropic.f and 213C 2) html-documentation 214C -> Making C symmetric in the **last two** indices: 215C 216C stiff stores: 217C 1111, 1122, 2222, 1133 218C 2233, 3333, 1112, 2212 219C 220C 1111 221 stiff(1) = C(1,1,1,1) 222C 1122 = 2211 223 stiff(2) = C(1,1,2,2) 224C 2222 225 stiff(3) = C(2,2,2,2) 226C 1133 = 3311 227 stiff(4) = C(1,1,3,3) 228C 2233 = 3322 229 stiff(5) = C(2,2,3,3) 230C 3333 231 stiff(6) = C(3,3,3,3) 232C 1112 = 1121 = 1211 = 2111 233 stiff(7) = 0.5d0*(C(1,1,1,2)+C(1,1,2,1)) 234C 2212 = 1222 = 2122 = 2221 235 stiff(8) = 0.5d0*(C(2,2,1,2)+C(2,2,2,1)) 236C 237C stiff stores: 238C 3312, 1212, 1113, 2213 239C 3313, 1213, 1313, 1123 240C 241C 3312 = 1233 = 2133 = 3321 242 stiff(9) = 0.5d0*(C(3,3,1,2)+C(3,3,2,1)) 243C 1212 = 1221 = 2112 = 2121 244 stiff(10) = 0.5d0*(C(1,2,1,2)+C(1,2,2,1)) 245C 1113 = 1131 = 1311 = 3111 246 stiff(11) = 0.5d0*(C(1,1,1,3)+C(1,1,3,1)) 247C 2213 = 1322 = 2231 = 3122 248 stiff(12) = 0.5d0*(C(2,2,1,3)+C(2,2,3,1)) 249C 250C 3313 = 1333 = 3133 = 3331 251 stiff(13) = 0.5d0*(C(3,3,1,3)+C(3,3,3,1)) 252C 1213 = 1231 = 1312 = 1321 = 2113 = 2131 = 3112 = 3121: 253 stiff(14) = 0.5d0*(C(1,2,1,3)+C(1,2,3,1)) 254C 1313 = 1331 = 3113 = 3131 255 stiff(15) = 0.5d0*(C(1,3,1,3)+C(1,3,3,1)) 256C 1123 = 1132 = 2311 = 3211 257 stiff(16) = 0.5d0*(C(1,1,2,3)+C(1,1,3,2)) 258C 259C stiff stores: 260C 2223, 3323, 1223, 1323 261C 2323 262C 263C 2223 = 2232 = 2322 = 3222 264 stiff(17) = 0.5d0*(C(2,2,2,3)+C(2,2,3,2)) 265C 3323 = 2333 = 3233 = 3332 266 stiff(18) = 0.5d0*(C(3,3,2,3)+C(3,3,3,2)) 267C 1223 = 1232 = 2123 = 2132 = 2312 = 2321 = 3212 = 3221 268 stiff(19) = 0.5d0*(C(1,2,2,3)+C(1,2,3,2)) 269C 1323 = 1332 = 2313 = 2331 = 3123 = 3132 = 3213 = 3231 270 stiff(20) = 0.5d0*(C(1,3,2,3)+C(1,3,3,2)) 271C 2323 = 2332 272 stiff(21) = 0.5d0*(C(2,3,2,3)+C(2,3,3,2)) 273 274 return 275 end 276 277 subroutine S2_Ciarlet(Cb,lamda,mu, S2,Cbi,detC) 278C computes 2PK stresses for Ciarlet-model 279C 280 implicit none 281C input: 282 double precision Cb(3,3),lamda,mu 283C output: 284 double precision S2(3,3), Cbi(3,3), detC 285C local: 286 double precision id(3,3) 287 double precision f1 288 integer I,J 289 logical OK_FLAG 290 291C Set id: 292 do I=1,3 293 do J=1,3 294 id(I,J)=0.d0 295 enddo 296 id(I,I)=1.d0 297 enddo 298 299C Inverse of Cb and det(C): 300 call M33INV_DET(Cb, Cbi, detC, OK_FLAG) 301C 302C 303 f1 = 0.5d0*lamda*(detC-1.d0) 304 do I=1,3 305 do J=1,3 306 S2(I,J)=f1*Cbi(I,J) + mu * ( id(I,J)-Cbi(I,J) ) 307 enddo 308 enddo 309 return 310 end 311 312 subroutine dS2_Ciarlet_dE(lamda,mu,detC,Cbi,A) 313C computes linearization of 2PK stresses with E 314C for Ciarlet-model 315 implicit none 316C input: 317 double precision lamda,mu,detC,Cbi(3,3) 318C output: 319 double precision A(3,3,3,3) 320C local: 321 double precision f1,f2 322 integer I,J,K,L 323 324 f1 = lamda*detC 325 f2 = 2.d0*mu - lamda*(detC-1.d0) 326 do I=1,3 327 do J=1,3 328 do K=1,3 329 do L=1,3 330 A(I,J,K,L) = 331 & f1 * Cbi(K,L)*Cbi(I,J) + 332 & f2 * Cbi(I,K)*Cbi(L,J) 333 enddo 334 enddo 335 enddo 336 enddo 337 return 338 end 339 340 subroutine M33INV_DET(A, AINV, DET, OK_FLAG) 341C 342C 343C !****************************************************************** 344C ! M33INV_DET - Compute the inverse of a 3x3 matrix. 345C ! 346C ! A = input 3x3 matrix to be inverted 347C ! AINV = output 3x3 inverse of matrix A 348C ! OK_FLAG = (output) .TRUE. if the input matrix could be inverted, 349C ! and .FALSE. if the input matrix is singular. 350C !****************************************************************** 351 352 IMPLICIT NONE 353 354 DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: A 355 DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: AINV 356 LOGICAL, INTENT(OUT) :: OK_FLAG 357 358 DOUBLE PRECISION, PARAMETER :: EPS = 1.0D-10 359 DOUBLE PRECISION :: DET 360 DOUBLE PRECISION, DIMENSION(3,3) :: COFACTOR 361 362 363 DET = A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2) 364 & - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1) 365 & + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1) 366 367 IF (ABS(DET) .LE. EPS) THEN 368 AINV = 0.0D0 369 OK_FLAG = .FALSE. 370 RETURN 371 END IF 372 373 COFACTOR(1,1) = +(A(2,2)*A(3,3)-A(2,3)*A(3,2)) 374 COFACTOR(1,2) = -(A(2,1)*A(3,3)-A(2,3)*A(3,1)) 375 COFACTOR(1,3) = +(A(2,1)*A(3,2)-A(2,2)*A(3,1)) 376 COFACTOR(2,1) = -(A(1,2)*A(3,3)-A(1,3)*A(3,2)) 377 COFACTOR(2,2) = +(A(1,1)*A(3,3)-A(1,3)*A(3,1)) 378 COFACTOR(2,3) = -(A(1,1)*A(3,2)-A(1,2)*A(3,1)) 379 COFACTOR(3,1) = +(A(1,2)*A(2,3)-A(1,3)*A(2,2)) 380 COFACTOR(3,2) = -(A(1,1)*A(2,3)-A(1,3)*A(2,1)) 381 COFACTOR(3,3) = +(A(1,1)*A(2,2)-A(1,2)*A(2,1)) 382 383 AINV = TRANSPOSE(COFACTOR) / DET 384 385 OK_FLAG = .TRUE. 386 387 RETURN 388 END 389 390 391