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 linel(kode,mattyp,beta,emec,stre,elas,elconloc, 20 & iorien,orab,pgauss) 21! 22! calculates stresses for linear elastic materials 23! 24 implicit none 25! 26 integer mattyp,j1,j2,j3,j4,j5,j6,j7,j8,j,jj,kel(4,21), 27 & iorien,i,kode 28! 29 real*8 beta(6),elas(21),stre(6),fxx,fyy,fzz,fxy,fxz,fyz, 30 & elconloc(*),emax,ya(3,3,3,3),orab(7,*),skl(3,3),e,un, 31 & um,um2,al,am1,pgauss(3),emec(6) 32! 33 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, 34 & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3, 35 & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3, 36 & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/)) 37! 38! engineering strain 39! 40 fxx=emec(1) 41 fyy=emec(2) 42 fzz=emec(3) 43 fxy=2.d0*emec(4) 44 fxz=2.d0*emec(5) 45 fyz=2.d0*emec(6) 46! 47 if(kode.eq.2) then 48! 49! isotropic 50! 51 do i=1,2 52 elas(i)=elconloc(i) 53 enddo 54! 55 e=elas(1) 56 un=elas(2) 57 um2=e/(1.d0+un) 58 al=un*um2/(1.d0-2.d0*un) 59 um=um2/2.d0 60 am1=al+um2 61! 62 stre(1)=am1*fxx+al*(fyy+fzz)-beta(1) 63 stre(2)=am1*fyy+al*(fxx+fzz)-beta(2) 64 stre(3)=am1*fzz+al*(fxx+fyy)-beta(3) 65 stre(4)=um*fxy-beta(4) 66 stre(5)=um*fxz-beta(5) 67 stre(6)=um*fyz-beta(6) 68! 69 mattyp=1 70! 71 elseif((kode.eq.9).or.(kode.eq.21)) then 72! 73 if((kode.eq.9).and.(iorien.eq.0)) then 74! 75! orthotropic 76! 77 do i=1,9 78 elas(i)=elconloc(i) 79 enddo 80 do i=10,21 81 elas(i)=0.d0 82 enddo 83! 84 stre(1)=elas(1)*fxx+elas(2)*fyy+ 85 & elas(4)*fzz-beta(1) 86 stre(2)=elas(2)*fxx+elas(3)*fyy+ 87 & elas(5)*fzz-beta(2) 88 stre(3)=elas(4)*fxx+elas(5)*fyy+ 89 & elas(6)*fzz-beta(3) 90 stre(4)=elas(7)*fxy-beta(4) 91 stre(5)=elas(8)*fxz-beta(5) 92 stre(6)=elas(9)*fyz-beta(6) 93! 94 mattyp=2 95! 96 else 97! 98 do i=1,21 99 elas(i)=elconloc(i) 100 enddo 101! 102 mattyp=3 103! 104 if(iorien.ne.0) then 105! 106! calculating the transformation matrix 107! 108 call transformatrix(orab(1,iorien),pgauss,skl) 109! 110! transforming the elastic coefficients 111! 112 if(kode.eq.9) then 113 call orthotropic(elas,ya) 114 else 115 call anisotropic(elas,ya) 116 endif 117! 118 do jj=1,21 119 j1=kel(1,jj) 120 j2=kel(2,jj) 121 j3=kel(3,jj) 122 j4=kel(4,jj) 123 elas(jj)=0.d0 124 do j5=1,3 125 do j6=1,3 126 do j7=1,3 127 do j8=1,3 128 elas(jj)=elas(jj)+ya(j5,j6,j7,j8)* 129 & skl(j1,j5)*skl(j2,j6)*skl(j3,j7)* 130 & skl(j4,j8) 131 enddo 132 enddo 133 enddo 134 enddo 135 enddo 136! 137! determining the type: orthotropic or anisotropic 138! 139 emax=0.d0 140 do j=1,21 141 emax=max(emax,dabs(elas(j))) 142 enddo 143 do j=7,9 144 if(dabs(elas(j)).gt.emax*1.d-10) then 145 emax=-1.d0 146 exit 147 endif 148 enddo 149 if(emax.ge.0.d0) then 150 do j=11,14 151 if(dabs(elas(j)).gt.emax*1.d-10) then 152 emax=-1.d0 153 exit 154 endif 155 enddo 156 endif 157 if(emax.ge.0.d0) then 158 do j=16,20 159 if(dabs(elas(j)).gt.emax*1.d-10) then 160 emax=-1.d0 161 exit 162 endif 163 enddo 164 endif 165 if(emax.ge.0.d0) then 166 elas(7)=elas(10) 167 elas(8)=elas(15) 168 elas(9)=elas(21) 169! 170 do j=10,21 171 elas(j)=0.d0 172 enddo 173c elas(10)=0.d0 174c elas(15)=0.d0 175c elas(21)=0.d0 176! 177 mattyp=2 178 endif 179 endif 180! 181 if(mattyp.eq.2) then 182! 183! orthotropic 184! 185 stre(1)=elas(1)*fxx+elas(2)*fyy+ 186 & elas(4)*fzz-beta(1) 187 stre(2)=elas(2)*fxx+elas(3)*fyy+ 188 & elas(5)*fzz-beta(2) 189 stre(3)=elas(4)*fxx+elas(5)*fyy+ 190 & elas(6)*fzz-beta(3) 191 stre(4)=elas(7)*fxy-beta(4) 192 stre(5)=elas(8)*fxz-beta(5) 193 stre(6)=elas(9)*fyz-beta(6) 194 else 195! 196! fully anisotropic 197! 198 stre(1)=elas(1)*fxx+ 199 & elas(2)*fyy+ 200 & elas(4)*fzz+ 201 & elas(7)*fxy+ 202 & elas(11)*fxz+ 203 & elas(16)*fyz-beta(1) 204 stre(2)=elas(2)*fxx+ 205 & elas(3)*fyy+ 206 & elas(5)*fzz+ 207 & elas(8)*fxy+ 208 & elas(12)*fxz+ 209 & elas(17)*fyz-beta(2) 210 stre(3)=elas(4)*fxx+ 211 & elas(5)*fyy+ 212 & elas(6)*fzz+ 213 & elas(9)*fxy+ 214 & elas(13)*fxz+ 215 & elas(18)*fyz-beta(3) 216 stre(4)=elas(7)*fxx+ 217 & elas(8)*fyy+ 218 & elas(9)*fzz+ 219 & elas(10)*fxy+ 220 & elas(14)*fxz+ 221 & elas(19)*fyz-beta(4) 222 stre(5)=elas(11)*fxx+ 223 & elas(12)*fyy+ 224 & elas(13)*fzz+ 225 & elas(14)*fxy+ 226 & elas(15)*fxz+ 227 & elas(20)*fyz-beta(5) 228 stre(6)=elas(16)*fxx+ 229 & elas(17)*fyy+ 230 & elas(18)*fzz+ 231 & elas(19)*fxy+ 232 & elas(20)*fxz+ 233 & elas(21)*fyz-beta(6) 234! 235 endif 236 endif 237 endif 238! 239 return 240 end 241 242