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_ideal_gas(amat,iel,iint,kode,elconloc,emec,emec0, 20 & beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime, 21 & icmd,ielas,mi, 22 & nstate_,xstateini,xstate,stre,stiff,iorien,pgauss,orab) 23! 24! calculates stiffness and stresses for an ideal gas 25! For this material there is just one material constant equal to 26! initial density x specific gas constant 27! The user should list this constant (which does not depend on 28! the temperature) 29! underneath the *USER MATERIAL,CONSTANTS=1 card. The name of the 30! material has to start with IDEAL_GAS, e.g. IDEAL_GAS_AIR or 31! IDEAL_GAS_NITROGEN etc. 32! 33! This routine should only be used with nlgeom=yes 34! 35! icmd=3: calculates stress at mechanical strain 36! else: calculates stress at mechanical strain and the stiffness 37! matrix 38! 39! INPUT: 40! 41! amat material name 42! iel element number 43! iint integration point number 44! 45! kode material type (-100-#of constants entered 46! under *USER MATERIAL): can be used for materials 47! with varying number of constants 48! 49! elconloc(21) user defined constants defined by the keyword 50! card *USER MATERIAL (max. 21, actual # = 51! -kode-100), interpolated for the 52! actual temperature t1l 53! 54! emec(6) Lagrange mechanical strain tensor (component order: 55! 11,22,33,12,13,23) at the end of the increment 56! (thermal strains are subtracted) 57! emec0(6) Lagrange mechanical strain tensor at the start of the 58! increment (thermal strains are subtracted) 59! beta(6) residual stress tensor (the stress entered under 60! the keyword *INITIAL CONDITIONS,TYPE=STRESS) 61! 62! xokl(3,3) deformation gradient at the start of the increment 63! voj Jacobian at the start of the increment 64! xkl(3,3) deformation gradient at the end of the increment 65! vj Jacobian at the end of the increment 66! 67! ithermal 0: no thermal effects are taken into account 68! >0: thermal effects are taken into account (triggered 69! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) 70! t1l temperature at the end of the increment 71! dtime time length of the increment 72! time step time at the end of the current increment 73! ttime total time at the start of the current step 74! 75! icmd not equal to 3: calculate stress and stiffness 76! 3: calculate only stress 77! ielas 0: no elastic iteration: irreversible effects 78! are allowed 79! 1: elastic iteration, i.e. no irreversible 80! deformation allowed 81! 82! mi(1) max. # of integration points per element in the 83! model 84! nstate_ max. # of state variables in the model 85! 86! xstateini(nstate_,mi(1),# of elements) 87! state variables at the start of the increment 88! xstate(nstate_,mi(1),# of elements) 89! state variables at the end of the increment 90! 91! stre(6) Piola-Kirchhoff stress of the second kind 92! at the start of the increment 93! 94! iorien number of the local coordinate axis system 95! in the integration point at stake (takes the value 96! 0 if no local system applies) 97! pgauss(3) global coordinates of the integration point 98! orab(7,*) description of all local coordinate systems. 99! If a local coordinate system applies the global 100! tensors can be obtained by premultiplying the local 101! tensors with skl(3,3). skl is determined by calling 102! the subroutine transformatrix: 103! call transformatrix(orab(1,iorien),pgauss,skl) 104! 105! 106! OUTPUT: 107! 108! xstate(nstate_,mi(1),# of elements) 109! updated state variables at the end of the increment 110! stre(6) Piola-Kirchhoff stress of the second kind at the 111! end of the increment 112! stiff(21): consistent tangent stiffness matrix in the material 113! frame of reference at the end of the increment. In 114! other words: the derivative of the PK2 stress with 115! respect to the Lagrangian strain tensor. The matrix 116! is supposed to be symmetric, only the upper half is 117! to be given in the same order as for a fully 118! anisotropic elastic material (*ELASTIC,TYPE=ANISO). 119! Notice that the matrix is an integral part of the 120! fourth order material tensor, i.e. the Voigt notation 121! is not used. 122! 123 implicit none 124! 125 character*80 amat 126! 127 integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien, 128 & i, 129 & kk(84),k,l,m,n,nt 130! 131 real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6), 132 & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), 133 & time,ttime,xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*), 134 & rho0r,c(3,3),cinv(3,3),v3,didc(3,3) 135! 136 kk=(/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3, 137 & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,3,3,1,3, 138 & 1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,1,2,2,3,1,3,2,3, 139 & 2,3,2,3/) 140! 141 rho0r=elconloc(1) 142! 143! calculation of the Green deformation tensor for the total 144! strain and the thermal strain 145! 146 do i=1,3 147 c(i,i)=emec(i)*2.d0+1.d0 148 enddo 149 c(1,2)=2.d0*emec(4) 150 c(1,3)=2.d0*emec(5) 151 c(2,3)=2.d0*emec(6) 152! 153! calculation of the third invariant of C 154! 155 v3=c(1,1)*(c(2,2)*c(3,3)-c(2,3)*c(2,3)) 156 & -c(1,2)*(c(1,2)*c(3,3)-c(1,3)*c(2,3)) 157 & +c(1,3)*(c(1,2)*c(2,3)-c(1,3)*c(2,2)) 158! 159! inversion of c 160! 161 cinv(1,1)=(c(2,2)*c(3,3)-c(2,3)*c(2,3))/v3 162 cinv(2,2)=(c(1,1)*c(3,3)-c(1,3)*c(1,3))/v3 163 cinv(3,3)=(c(1,1)*c(2,2)-c(1,2)*c(1,2))/v3 164 cinv(1,2)=(c(1,3)*c(2,3)-c(1,2)*c(3,3))/v3 165 cinv(1,3)=(c(1,2)*c(2,3)-c(2,2)*c(1,3))/v3 166 cinv(2,3)=(c(1,2)*c(1,3)-c(1,1)*c(2,3))/v3 167 cinv(2,1)=cinv(1,2) 168 cinv(3,1)=cinv(1,3) 169 cinv(3,2)=cinv(2,3) 170! 171! changing the meaning of v3 172! 173 v3=v3*rho0r 174! 175! stress at mechanical strain 176! 177 stre(1)=v3*cinv(1,1) 178 stre(2)=v3*cinv(2,2) 179 stre(3)=v3*cinv(3,3) 180 stre(4)=v3*cinv(1,2) 181 stre(5)=v3*cinv(1,3) 182 stre(6)=v3*cinv(2,3) 183! 184! tangent 185! 186 if(icmd.ne.3) then 187! 188! 189! second derivative of the c-invariants w.r.t. c(k,l) 190! and c(m,n) 191! 192 if(icmd.ne.3) then 193 nt=0 194 do i=1,21 195 k=kk(nt+1) 196 l=kk(nt+2) 197 m=kk(nt+3) 198 n=kk(nt+4) 199 nt=nt+4 200 stiff(i)=v3*(cinv(m,n)*cinv(k,l)- 201 & (cinv(k,m)*cinv(n,l)+cinv(k,n)*cinv(m,l))/2.d0) 202 enddo 203 endif 204 endif 205! 206 return 207 end 208