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! 20! center of gravity of the projection of the vertices for 21! visibility purposes 22! exact integration for one triangle: routine cubtri 23! if the surfaces are far enough away, one-point integration 24! is used 25! 26 subroutine radmatrix(ntr,adrad,aurad,bcr,sideload,nelemload, 27 & xloadact,lakon,vold,ipkon,kon,co,nloadtr,tarea,tenv,physcon, 28 & erad,adview,auview,ithermal,iinc,iit,fenv,istep,dtime,ttime, 29 & time,iviewfile,xloadold,reltime,nmethod,mi,iemchange,nam, 30 & iamload,jqrad,irowrad,nzsrad) 31! 32 implicit none 33! 34 character*8 lakonl,lakon(*) 35 character*20 sideload(*) 36! 37 integer ntr,nelemload(2,*),nope,nopes,mint2d,i,j,k,l, 38 & node,ifaceq(8,6),ifacet(6,4),iviewfile,mi(*), 39 & ifacew(8,5),nelem,ig,index,konl(20),iflag, 40 & ipkon(*),kon(*),nloadtr(*),i1,ithermal(*),iinc,iit, 41 & istep,jltyp,nfield,nmethod,iemchange,nam, 42 & iamload(2,*),irowrad(*),jqrad(*),nzsrad 43! 44 real*8 adrad(*),aurad(*),bcr(ntr,1),xloadact(2,*),h(2), 45 & xl2(3,8),coords(3),dxsj2,temp,xi,et,weight,xsj2(3), 46 & vold(0:mi(2),*),co(3,*),shp2(7,8),xs2(3,7), 47 & areamean,tarea(*),tenv(*),erad(*),fenv(*),physcon(*), 48 & adview(*),auview(*),tl2(8),tvar(2),field, 49 & dtime,ttime,time,areaj,xloadold(2,*),reltime, 50 & fform 51! 52 data ifaceq /4,3,2,1,11,10,9,12, 53 & 5,6,7,8,13,14,15,16, 54 & 1,2,6,5,9,18,13,17, 55 & 2,3,7,6,10,19,14,18, 56 & 3,4,8,7,11,20,15,19, 57 & 4,1,5,8,12,17,16,20/ 58 data ifacet /1,3,2,7,6,5, 59 & 1,2,4,5,9,8, 60 & 2,3,4,6,10,9, 61 & 1,4,3,8,10,7/ 62 data ifacew /1,3,2,9,8,7,0,0, 63 & 4,5,6,10,11,12,0,0, 64 & 1,2,5,4,7,14,10,13, 65 & 2,3,6,5,8,15,11,14, 66 & 4,6,3,1,12,15,9,13/ 67 data iflag /2/ 68! 69 external fform 70! 71 include "gauss.f" 72! 73 tvar(1)=time 74 tvar(2)=ttime+time 75! 76! filling acr and bcr 77! 78 do i1=1,ntr 79 i=nloadtr(i1) 80 nelem=nelemload(1,i) 81 lakonl=lakon(nelem) 82! 83! calculate the mean temperature of the face 84! 85 read(sideload(i)(2:2),'(i1)') ig 86! 87! number of nodes and integration points in the face 88! 89 if(lakonl(4:4).eq.'2') then 90 nope=20 91 nopes=8 92 elseif(lakonl(4:4).eq.'8') then 93 nope=8 94 nopes=4 95 elseif(lakonl(4:5).eq.'10') then 96 nope=10 97 nopes=6 98 elseif(lakonl(4:4).eq.'4') then 99 nope=4 100 nopes=3 101 elseif(lakonl(4:5).eq.'15') then 102 nope=15 103 else 104 nope=6 105 endif 106! 107 if(lakonl(4:5).eq.'8R') then 108 mint2d=1 109 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) 110 & then 111 if(lakonl(7:7).eq.'A') then 112 mint2d=2 113 else 114 mint2d=4 115 endif 116 elseif(lakonl(4:4).eq.'2') then 117 mint2d=9 118 elseif(lakonl(4:5).eq.'10') then 119 mint2d=3 120 elseif(lakonl(4:4).eq.'4') then 121 mint2d=1 122 endif 123! 124 if(lakonl(4:4).eq.'6') then 125 mint2d=1 126 if(ig.le.2) then 127 nopes=3 128 else 129 nopes=4 130 endif 131 endif 132 if(lakonl(4:5).eq.'15') then 133 if(ig.le.2) then 134 mint2d=3 135 nopes=6 136 else 137 mint2d=4 138 nopes=8 139 endif 140 endif 141! 142! connectivity of the element 143! 144 index=ipkon(nelem) 145 if(index.lt.0) then 146 write(*,*) '*ERROR in radmatrix: element ',nelem 147 write(*,*) ' is not defined' 148 call exit(201) 149 endif 150 do k=1,nope 151 konl(k)=kon(index+k) 152 enddo 153! 154! coordinates of the nodes belonging to the face 155! 156 if((nope.eq.20).or.(nope.eq.8)) then 157 do k=1,nopes 158 tl2(k)=vold(0,konl(ifaceq(k,ig))) 159! 160 do j=1,3 161 xl2(j,k)=co(j,konl(ifaceq(k,ig)))+ 162 & vold(j,konl(ifaceq(k,ig))) 163 enddo 164 enddo 165 elseif((nope.eq.10).or.(nope.eq.4)) then 166 do k=1,nopes 167 tl2(k)=vold(0,konl(ifacet(k,ig))) 168 do j=1,3 169 xl2(j,k)=co(j,konl(ifacet(k,ig)))+ 170 & vold(j,konl(ifacet(k,ig))) 171 enddo 172 enddo 173 else 174 do k=1,nopes 175 tl2(k)=vold(0,konl(ifacew(k,ig))) 176 do j=1,3 177 xl2(j,k)=co(j,konl(ifacew(k,ig)))+ 178 & vold(j,konl(ifacew(k,ig))) 179 enddo 180 enddo 181 endif 182! 183! integration to obtain the center of gravity and the mean 184! temperature; radiation coefficient 185! 186 areamean=0.d0 187 tarea(i1)=0.d0 188! 189 read(sideload(i)(2:2),'(i1)') jltyp 190 jltyp=jltyp+10 191 if(sideload(i)(5:6).ne.'NU') then 192 erad(i1)=xloadact(1,i) 193! 194! if an amplitude was defined for the emissivity it is 195! assumed that the emissivity changes with the step, so 196! acr has to be calculated anew in every iteration 197! 198 if(nam.gt.0) then 199 if(iamload(1,i).ne.0) then 200 iemchange=1 201 endif 202 endif 203 else 204 erad(i1)=0.d0 205 endif 206! 207 do l=1,mint2d 208 if((lakonl(4:5).eq.'8R').or. 209 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 210 xi=gauss2d1(1,l) 211 et=gauss2d1(2,l) 212 weight=weight2d1(l) 213 elseif((lakonl(4:4).eq.'8').or. 214 & (lakonl(4:6).eq.'20R').or. 215 & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then 216 xi=gauss2d2(1,l) 217 et=gauss2d2(2,l) 218 weight=weight2d2(l) 219 elseif(lakonl(4:4).eq.'2') then 220 xi=gauss2d3(1,l) 221 et=gauss2d3(2,l) 222 weight=weight2d3(l) 223 elseif((lakonl(4:5).eq.'10').or. 224 & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then 225 xi=gauss2d5(1,l) 226 et=gauss2d5(2,l) 227 weight=weight2d5(l) 228 elseif((lakonl(4:4).eq.'4').or. 229 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 230 xi=gauss2d4(1,l) 231 et=gauss2d4(2,l) 232 weight=weight2d4(l) 233 endif 234! 235 if(nopes.eq.8) then 236 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 237 elseif(nopes.eq.4) then 238 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 239 elseif(nopes.eq.6) then 240 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 241 else 242 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 243 endif 244! 245 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 246 & xsj2(3)*xsj2(3)) 247! 248 temp=0.d0 249 do j=1,nopes 250 temp=temp+tl2(j)*shp2(4,j) 251 enddo 252! 253 tarea(i1)=tarea(i1)+temp*dxsj2*weight 254 areamean=areamean+dxsj2*weight 255! 256 if(sideload(i)(5:6).eq.'NU') then 257 areaj=dxsj2*weight 258 do k=1,3 259 coords(k)=0.d0 260 enddo 261 do j=1,nopes 262 do k=1,3 263 coords(k)=coords(k)+xl2(k,j)*shp2(4,j) 264 enddo 265 enddo 266 call radiate(h(1),tenv(i1),temp,istep, 267 & iinc,tvar,nelem,l,coords,jltyp,field,nfield, 268 & sideload(i),node,areaj,vold,mi,iemchange) 269 if(nmethod.eq.1) h(1)=xloadold(1,i)+ 270 & (h(1)-xloadold(1,i))*reltime 271 erad(i1)=erad(i1)+h(1) 272 endif 273! 274 enddo 275 tarea(i1)=tarea(i1)/areamean-physcon(1) 276 if(sideload(i)(5:6).eq.'NU') then 277 erad(i1)=erad(i1)/mint2d 278 endif 279! 280! updating the right hand side 281! 282 bcr(i1,1)=physcon(2)*(erad(i1)*tarea(i1)**4+ 283 & (1.d0-erad(i1))*fenv(i1)*tenv(i1)**4) 284c write(*,*) 'radmatrix ',bcr(i1,1) 285! 286 enddo 287! 288! adrad and aurad is recalculated only if the viewfactors 289! or the emissivity changed 290! 291 if(((ithermal(1).eq.3).and.(iviewfile.ge.0)).or. 292 & ((iit.eq.-1).and.(iviewfile.ne.-2)).or.(iemchange.eq.1).or. 293 & ((iit.eq.0).and.(abs(nmethod).eq.1))) then 294! 295 do i=1,ntr 296 erad(i)=1.d0-erad(i) 297 enddo 298! 299! diagonal entries 300! 301 do i=1,ntr 302 adrad(i)=1.d0-erad(i)*adview(i) 303 enddo 304! 305! lower triangular entries 306! 307 do i=1,nzsrad 308 aurad(i)=-erad(irowrad(i))*auview(i) 309 enddo 310! 311! upper triangular entries 312! 313 do i=1,ntr 314 do j=nzsrad+jqrad(i),nzsrad+jqrad(i+1)-1 315 aurad(j)=-erad(i)*auview(j) 316 enddo 317 enddo 318! 319 do i=1,ntr 320 erad(i)=1.d0-erad(i) 321 enddo 322 endif 323! 324 return 325 end 326