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 jouleheating(ipkon,lakon,kon,co,elcon,nelcon, 20 & mi,ne,sti,ielmat,nelemload,sideload,xload,nload,nload_, 21 & iamload,nam,idefload,ncmat_,ntmat_, 22 & alcon,nalcon,ithermal,vold,t1,nmethod) 23! 24! determines the effect of Joule heating 25! 26 implicit none 27! 28 character*8 lakon(*) 29 character*20 label,sideload(*) 30! 31 integer ipkon(*),nelem,kon(*),mi(*),nope,indexe,j,k,null, 32 & mint3d,jj,iflag,ne,nelemload(2,*),iamload(2,*),nload,nload_, 33 & ielmat(mi(3),*),konl(20),idefload(*),iamplitude,isector,nam, 34 & one,nelcon(2,*),nalcon(2,*),ithermal(*),i1,ncmat_,ntmat_, 35 & imat,nmethod 36! 37 real*8 co(3,*),xl(3,20),xi,et,ze,xsj,shp(4,20),weight,xload(2,*), 38 & sti(6,mi(1),*),alpha(6),heat,elcon(0:ncmat_,ntmat_,*),volume, 39 & elconloc(21),t1l,alcon(0:6,ntmat_,*),vold(0:mi(2),*),t1(*) 40! 41 include "gauss.f" 42! 43 data iflag /2/ 44! 45 null=0 46 one=1 47! 48 do nelem=1,ne 49 if(ipkon(nelem).lt.0) cycle 50! 51! only elements belonging to the A,V-domain experience 52! Joule heating 53! 54 if(int(elcon(2,1,ielmat(1,nelem))).ne.2) cycle 55! 56 imat=ielmat(1,nelem) 57 indexe=ipkon(nelem) 58! 59 if(lakon(nelem)(1:5).eq.'C3D8I') then 60 nope=11 61 elseif(lakon(nelem)(4:4).eq.'2') then 62 nope=20 63 elseif(lakon(nelem)(4:4).eq.'8') then 64 nope=8 65 elseif(lakon(nelem)(4:5).eq.'10') then 66 nope=10 67 elseif(lakon(nelem)(4:4).eq.'4') then 68 nope=4 69 elseif(lakon(nelem)(4:5).eq.'15') then 70 nope=15 71 elseif(lakon(nelem)(4:5).eq.'6') then 72 nope=6 73 endif 74! 75 do j=1,nope 76 konl(j)=kon(indexe+j) 77 do k=1,3 78 xl(k,j)=co(k,konl(j)) 79 enddo 80 enddo 81! 82 heat=0.d0 83 volume=0.d0 84! 85 if(lakon(nelem)(4:5).eq.'8R') then 86 mint3d=1 87 elseif((lakon(nelem)(4:4).eq.'8').or. 88 & (lakon(nelem)(4:6).eq.'20R')) then 89 mint3d=8 90 elseif(lakon(nelem)(4:4).eq.'2') then 91 mint3d=27 92 elseif(lakon(nelem)(4:5).eq.'10') then 93 mint3d=4 94 elseif(lakon(nelem)(4:4).eq.'4') then 95 mint3d=1 96 elseif(lakon(nelem)(4:5).eq.'15') then 97 mint3d=9 98 elseif(lakon(nelem)(4:5).eq.'6') then 99 mint3d=2 100 endif 101! 102! loop over the integration points 103! 104 do jj=1,mint3d 105 if(lakon(nelem)(4:5).eq.'8R') then 106 xi=gauss3d1(1,jj) 107 et=gauss3d1(2,jj) 108 ze=gauss3d1(3,jj) 109 weight=weight3d1(jj) 110 elseif((lakon(nelem)(4:4).eq.'8').or. 111 & (lakon(nelem)(4:6).eq.'20R')) 112 & then 113 xi=gauss3d2(1,jj) 114 et=gauss3d2(2,jj) 115 ze=gauss3d2(3,jj) 116 weight=weight3d2(jj) 117 elseif(lakon(nelem)(4:4).eq.'2') then 118 xi=gauss3d3(1,jj) 119 et=gauss3d3(2,jj) 120 ze=gauss3d3(3,jj) 121 weight=weight3d3(jj) 122 elseif(lakon(nelem)(4:5).eq.'10') then 123 xi=gauss3d5(1,jj) 124 et=gauss3d5(2,jj) 125 ze=gauss3d5(3,jj) 126 weight=weight3d5(jj) 127 elseif(lakon(nelem)(4:4).eq.'4') then 128 xi=gauss3d4(1,jj) 129 et=gauss3d4(2,jj) 130 ze=gauss3d4(3,jj) 131 weight=weight3d4(jj) 132 elseif(lakon(nelem)(4:5).eq.'15') then 133 xi=gauss3d8(1,jj) 134 et=gauss3d8(2,jj) 135 ze=gauss3d8(3,jj) 136 weight=weight3d8(jj) 137 else 138 xi=gauss3d7(1,jj) 139 et=gauss3d7(2,jj) 140 ze=gauss3d7(3,jj) 141 weight=weight3d7(jj) 142 endif 143! 144 if(nope.eq.20) then 145 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 146 elseif(nope.eq.8) then 147 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 148 elseif(nope.eq.10) then 149 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 150 elseif(nope.eq.4) then 151 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 152 elseif(nope.eq.15) then 153 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 154 else 155 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 156 endif 157! 158! calculating the temperature 159! 160 t1l=0.d0 161 if(ithermal(1).eq.1) then 162 if(lakon(nelem)(4:5).eq.'8 ') then 163 do i1=1,nope 164 t1l=t1l+t1(konl(i1))/8.d0 165 enddo 166 elseif(lakon(nelem)(4:6).eq.'20 ')then 167 call linscal(t1,konl,nope,jj,t1l,one) 168 elseif(lakon(nelem)(4:6).eq.'10T') then 169 call linscal10(t1,konl,t1l,null,shp) 170 else 171 do i1=1,nope 172 t1l=t1l+shp(4,i1)*t1(konl(i1)) 173 enddo 174 endif 175 elseif(ithermal(1).ge.2) then 176 if(lakon(nelem)(4:5).eq.'8 ') then 177 do i1=1,nope 178 t1l=t1l+vold(0,konl(i1))/8.d0 179 enddo 180 elseif(lakon(nelem)(4:6).eq.'20 ')then 181 call linscal(vold,konl,nope,jj,t1l,mi(2)) 182 elseif(lakon(nelem)(4:6).eq.'10T') then 183 call linscal10(vold,konl,t1l,mi(2),shp) 184 else 185 do i1=1,nope 186 t1l=t1l+shp(4,i1)*vold(0,konl(i1)) 187 enddo 188 endif 189 endif 190! 191! material data (electric conductivity and 192! magnetic permeability) 193! 194 call materialdata_em(elcon,nelcon,alcon,nalcon, 195 & imat,ntmat_,t1l,elconloc,ncmat_,alpha) 196! 197 if(nmethod.ne.2) then 198! 199! time dependent current: 200! heat=sigma*E**2 201! 202 heat=heat+weight*xsj*alpha(1)* 203 & (sti(1,jj,nelem)**2+ 204 & sti(2,jj,nelem)**2+ 205 & sti(3,jj,nelem)**2) 206 else 207! 208! alernating current: 209! heat=sigma*(E_real**2+E_imaginary**2) 210! 211 heat=heat+weight*xsj*alpha(1)* 212 & (sti(1,jj,nelem)**2+ 213 & sti(2,jj,nelem)**2+ 214 & sti(3,jj,nelem)**2+ 215 & sti(1,jj,nelem+ne)**2+ 216 & sti(2,jj,nelem+ne)**2+ 217 & sti(3,jj,nelem+ne)**2) 218 endif 219 volume=volume+weight*xsj 220 enddo 221! 222 heat=heat/volume 223! 224! adding the Joule heating to the distributed loading 225! 226 label='BF ' 227 iamplitude=0 228 isector=0 229 call loadadd(nelem,label,heat,nelemload,sideload, 230 & xload,nload,nload_,iamload,iamplitude,nam,isector, 231 & idefload) 232 enddo 233! 234 return 235 end 236 237