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 defplas(elconloc,elas,emec,ithermal,icmd, 20 & beta,stre,ckl,vj) 21! 22! calculates stiffness and stresses for the deformation plasticity 23! material law 24! 25! icmd=3: calculates stress at mechanical strain 26! else: calculates stress and stiffness matrix at mechanical strain 27! 28 implicit none 29! 30 integer cauchy 31! 32 integer ithermal(*),icmd,i,j,k,l,m,n,ii,istart,iend,nt,kk(84) 33! 34 real*8 elconloc(*),elas(*),emec(*),beta(*),s(6),al, 35 & ee,un,s0,xn,stre(*),eq,c0,c1,c2,c3,dkl(3,3),ekl(3,3), 36 & q,dq,pp,el(6),ckl(3,3),vj 37! 38 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, 39 & 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, 40 & 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, 41 & 2,3,2,3/) 42! 43 cauchy=1 44! 45 istart=1 46 iend=1 47! 48! determining linear elastic material constants 49! 50 ee=elconloc(1) 51 un=elconloc(2) 52 s0=elconloc(3) 53 xn=elconloc(4) 54 al=elconloc(5) 55! 56 do i=1,6 57 el(i)=emec(i) 58 enddo 59! 60! major loop 61! 62 do ii=istart,iend 63! 64 c0=(el(1)+el(2)+el(3))/3.d0 65! 66 el(1)=el(1)-c0 67 el(2)=el(2)-c0 68 el(3)=el(3)-c0 69! 70! equivalent deviatoric strain 71! 72 eq=dsqrt(2.d0/3.d0*(el(1)*el(1)+el(2)*el(2)+ 73 & el(3)*el(3)+2.d0*(el(4)*el(4)+ 74 & el(5)*el(5)+el(6)*el(6)))) 75! 76! initial value of the Mises equivalent stress (q) 77! 78 c1=3.d0*ee*eq/(2.d0*(1.d0+un)) 79! 80 if(c1.le.s0) then 81 q=c1 82 else 83 q=(s0**(xn-1)*ee*eq/al)**(1.d0/xn) 84 endif 85! 86! determining the Mises equivalent stress q 87! 88 c1=2.d0*(1.d0+un)/3.d0 89 do 90 c2=al*(q/s0)**(xn-1.d0) 91 dq=(ee*eq-(c1+c2)*q)/(c1+xn*c2) 92 if((dabs(dq).lt.q*1.d-4).or.(dabs(dq).lt.1.d-10)) exit 93 q=q+dq 94 enddo 95! 96 if(icmd.ne.3) then 97! 98! calculating the tangent stiffness matrix 99! 100! initialization of the Delta Dirac function 101! 102 do i=1,3 103 do j=1,3 104 dkl(i,j)=0.d0 105 enddo 106 enddo 107 do i=1,3 108 dkl(i,i)=1.d0 109 enddo 110! 111 ekl(1,1)=el(1) 112 ekl(2,2)=el(2) 113 ekl(3,3)=el(3) 114 ekl(1,2)=el(4) 115 ekl(1,3)=el(5) 116 ekl(2,3)=el(6) 117 ekl(2,1)=ekl(1,2) 118 ekl(3,1)=ekl(1,3) 119 ekl(3,2)=ekl(2,3) 120! 121 if(eq.lt.1.d-10) then 122 c1=ee/(1.d0+un) 123 c2=0.d0 124 else 125 c1=2.d0/(3.d0*eq) 126 c2=c1*(1.d0/eq-1.d0/(eq+(xn-1.d0)*c2*q/ee)) 127 c1=c1*q 128 endif 129 c3=(ee/(1.d0-2.d0*un)-c1)/3.d0 130! 131 nt=0 132 do i=1,21 133 k=kk(nt+1) 134 l=kk(nt+2) 135 m=kk(nt+3) 136 n=kk(nt+4) 137 nt=nt+4 138 elas(i)=c1*((dkl(k,m)*dkl(l,n)+dkl(k,n)*dkl(l,m))/2.d0 139 & -c2*ekl(k,l)*ekl(m,n)) 140 & +c3*dkl(k,l)*dkl(m,n) 141 enddo 142! 143! conversion of the stiffness matrix from spatial coordinates 144! coordinates into material coordinates 145! 146 call stiff2mat(elas,ckl,vj,cauchy) 147! 148 endif 149! 150! calculating the stress 151! 152 pp=-ee*c0/(1.d0-2.d0*un) 153! 154 if(eq.lt.1.d-10) then 155 c1=0.d0 156 else 157 c1=2.d0*q/(3.d0*eq) 158 endif 159! 160 do i=1,6 161 s(i)=el(i)*c1 162 enddo 163 do i=1,3 164 s(i)=s(i)-pp 165 enddo 166! 167 do i=1,6 168 stre(i)=s(i) 169 enddo 170! 171! converting the stress into the material frame of 172! reference 173! 174 call str2mat(stre,ckl,vj,cauchy) 175! 176 enddo 177! 178 return 179 end 180