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 calcstressheatfluxfem(kon,lakon,ipkon, 20 & ielmat,ntmat_,vold,matname,mi,shcon,nshcon, 21 & iturbulent,compressible,ipvar,var,sti,qfx,cocon, 22 & ncocon,ne,isti,iqfx,ithermal,rhcon,nrhcon,vcon,nk) 23! 24! calculating the viscous stresses and the heat flow 25! in the integration points (CBS-method) 26! 27 implicit none 28! 29 character*8 lakon(*) 30 character*80 matname(*),amat 31! 32 integer kon(*),nelem,mi(*),ne,konl(20),ipkon(*),j,i1,i2,j1,ii, 33 & jj,indexe,isti,iqfx,kk,ielmat(mi(3),*),nshcon(*),ntmat_,nope, 34 & imat,ncocon(2,*),mint3d,k1,ipvar(*),index,iturbulent,nk, 35 & compressible,ithermal(*),nrhcon(*) 36! 37 real*8 shp(4,20),dvi,cond,cocon(0:6,ntmat_,*),div,arg2,c1,c2, 38 & shcon(0:3,ntmat_,*),voldl(0:mi(2),20),vold(0:mi(2),*),temp, 39 & tt(3,3),rho,t(3,3),vkl(3,3),xkin,unt,umt,f2,a1,vort,xtuf, 40 & var(*),sti(6,mi(1),*),dtem(3),qfx(3,mi(1),*),y, 41 & rhcon(0:1,ntmat_,*),vcon(nk,0:mi(2)),vconl(0:mi(2),8) 42! 43 if(iturbulent.gt.0) then 44 a1=0.31d0 45 endif 46! 47 do nelem=1,ne 48! 49 if(ipkon(nelem).lt.0) cycle 50 if(lakon(nelem)(1:1).ne.'F') cycle 51 indexe=ipkon(nelem) 52! 53 imat=ielmat(1,nelem) 54 amat=matname(imat) 55! 56 if(lakon(nelem)(4:4).eq.'4') then 57 nope=4 58 mint3d=1 59 elseif(lakon(nelem)(4:4).eq.'6') then 60 nope=6 61 mint3d=2 62 elseif(lakon(nelem)(4:5).eq.'8R') then 63 nope=8 64 mint3d=1 65 elseif(lakon(nelem)(4:4).eq.'8') then 66 nope=8 67 mint3d=8 68 endif 69! 70 do j=1,nope 71 konl(j)=kon(indexe+j) 72 enddo 73! 74! storing the local temperature and velocity 75! 76 do i1=1,nope 77 do i2=0,3 78 voldl(i2,i1)=vold(i2,konl(i1)) 79 enddo 80 enddo 81! 82! storing the local turbulent kinetic energy and turbulence 83! frequency 84! 85 if(iturbulent.gt.0) then 86 do i1=1,nope 87 do i2=5,mi(2) 88 voldl(i2,i1)=vold(i2,konl(i1)) 89 enddo 90 vconl(4,i1)=vcon(konl(i1),4) 91 enddo 92 endif 93! 94! computation of the matrix: loop over the Gauss points 95! 96 index=ipvar(nelem) 97 do kk=1,mint3d 98! 99! copying the shape functions and their derivatives from field var 100! 101 do jj=1,nope 102 do ii=1,4 103 index=index+1 104 shp(ii,jj)=var(index) 105 enddo 106 enddo 107 index=index+2 108 y=var(index) 109! 110! calculating of 111! the velocity gradient vkl 112! the temperature gradient dtem 113! 114 if(isti.gt.0) then 115 do i1=1,3 116 do j1=1,3 117 vkl(i1,j1)=0.d0 118 enddo 119 enddo 120 do i1=1,nope 121 do j1=1,3 122 do k1=1,3 123 vkl(j1,k1)=vkl(j1,k1)+shp(k1,i1)*voldl(j1,i1) 124 enddo 125 enddo 126 enddo 127 if(compressible.eq.1) div=vkl(1,1)+vkl(2,2)+vkl(3,3) 128 endif 129! 130 if(iqfx.gt.0) then 131 do i1=1,3 132 dtem(i1)=0.d0 133 enddo 134 do i1=1,nope 135 do j1=1,3 136 dtem(j1)=dtem(j1)+shp(j1,i1)*voldl(0,i1) 137 enddo 138 enddo 139 endif 140! 141! calculating the temperature 142! 143 temp=0.d0 144! 145 do i1=1,nope 146 temp=temp+shp(4,i1)*voldl(0,i1) 147 enddo 148! 149! determining the dissipative stress 150! 151 if(isti.gt.0) then 152 call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon, 153 & dvi) 154 do i1=1,3 155 do j1=i1,3 156 t(i1,j1)=vkl(i1,j1)+vkl(j1,i1) 157 enddo 158 if(compressible.eq.1) t(i1,i1)=t(i1,i1)-2.d0*div/3.d0 159 enddo 160! 161! calculating the stress 162! 163 if(iturbulent.gt.0) then 164! 165! calculation of the density (liquid or gas) 166! 167 if(compressible.eq.1) then 168 rho=0.d0 169 do i1=1,nope 170 rho=rho+shp(4,i1)*vconl(4,i1) 171 enddo 172 else 173 call materialdata_rho(rhcon,nrhcon,imat,rho, 174 & temp,ntmat_,ithermal) 175 endif 176! 177! calculation of the turbulent kinetic energy, turbulence 178! frequency and turbulent kinematic viscosity 179! 180 xkin=0.d0 181 xtuf=0.d0 182 do i1=1,nope 183 xkin=xkin+shp(4,i1)*voldl(5,i1) 184 xtuf=xtuf+shp(4,i1)*voldl(6,i1) 185 enddo 186! 187! adding the turbulent stress 188! 189! factor F2 190! 191 c1=dsqrt(xkin)/(0.09d0*xtuf*y) 192 c2=500.d0*dvi/(y*y*xtuf*rho) 193! 194! kinematic turbulent viscosity 195! 196 if(iturbulent.eq.4) then 197! 198! vorticity 199! 200 vort=dsqrt((vkl(3,2)-vkl(2,3))**2+ 201 & (vkl(1,3)-vkl(3,1))**2+ 202 & (vkl(2,1)-vkl(1,2))**2) 203 arg2=max(2.d0*c1,c2) 204 f2=dtanh(arg2*arg2) 205 unt=a1*xkin/max(a1*xtuf,vort*f2) 206 else 207 unt=xkin/xtuf 208 endif 209! 210 umt=unt*rho 211! 212! calculating the turbulent stress 213! 214 do i1=1,3 215 do j1=i1,3 216 tt(i1,j1)=umt*t(i1,j1) 217 enddo 218 tt(i1,i1)=tt(i1,i1)-2.d0*rho*xkin/3.d0 219 enddo 220! 221! adding the viscous stress 222! 223 do i1=1,3 224 do j1=i1,3 225 t(i1,j1)=dvi*t(i1,j1)+tt(i1,j1) 226 enddo 227 enddo 228 else 229! 230 do i1=1,3 231 do j1=i1,3 232 t(i1,j1)=dvi*t(i1,j1) 233 enddo 234 enddo 235 endif 236! 237! storing the total stress 238! 239 sti(1,kk,nelem)=t(1,1) 240 sti(2,kk,nelem)=t(2,2) 241 sti(3,kk,nelem)=t(3,3) 242 sti(4,kk,nelem)=t(1,2) 243 sti(5,kk,nelem)=t(1,3) 244 sti(6,kk,nelem)=t(2,3) 245 endif 246! 247! storing the heat flow 248! 249 if(iqfx.gt.0) then 250 call materialdata_cond(imat,ntmat_,temp,cocon, 251 & ncocon,cond) 252 qfx(1,kk,nelem)=-cond*dtem(1) 253 qfx(2,kk,nelem)=-cond*dtem(2) 254 qfx(3,kk,nelem)=-cond*dtem(3) 255 endif 256! 257 enddo 258 enddo 259! 260 return 261 end 262 263