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 calcheatnet(nelem,lakon,ipkon,kon,v,ielprop,prop, 20 & ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,ipobody,ibody,xbody, 21 & mi,nacteq,bc,qat,nalt) 22! 23! user subroutine film 24! 25! 26! INPUT: 27! 28! nelem actual element number 29! lakon(i) contains the label of element i 30! ipkon(i) points to the location in field kon preceding 31! the topology of element i 32! kon(*) contains the topology of all elements. The 33! topology of element i starts at kon(ipkon(i)+1) 34! and continues until all nodes are covered. The 35! number of nodes depends on the element label 36! v(0..4,1..nk) actual solution field in all nodes; 37! for structural nodes: 38! 0: temperature 39! 1: displacement in global x-direction 40! 2: displacement in global y-direction 41! 3: displacement in global z-direction 42! 4: static pressure 43! for network nodes 44! 0: total temperature (at end nodes) 45! = static temperature for liquids 46! 1: mass flow (at middle nodes) 47! 2: total pressure (at end nodes) 48! = static pressure for liquids 49! 3: static temperature (at end nodes; only for gas) 50! ielprop(i) points to the location in field prop preceding 51! the properties of element i 52! prop(*) contains the properties of all network elements. The 53! properties of element i start at prop(ielprop(i)+1) 54! and continues until all properties are covered. The 55! appropriate amount of properties depends on the 56! element label. The kind of properties, their 57! number and their order corresponds 58! to the description in the user's manual, 59! cf. the sections "Fluid Section Types" 60! ielmat(j,i) contains the material number for element i 61! and layer j 62! ntmat_ maximum number of temperature data points for 63! any material property for any material 64! shcon(0,j,i) temperature at temperature point j of material i 65! shcon(1,j,i) specific heat at constant pressure at the 66! temperature point j of material i 67! shcon(2,j,i) dynamic viscosity at the temperature point j of 68! material i 69! shcon(3,1,i) specific gas constant of material i 70! nshcon(i) number of temperature data points for the specific 71! heat of material i 72! rhcon(0,j,i) temperature at density temperature point j of 73! material i 74! rhcon(1,j,i) density at the density temperature point j of 75! material i 76! nrhcon(i) number of temperature data points for the density 77! of material i 78! ipobody(1,i) points to an entry in fields ibody and xbody 79! containing the body load applied to element i, 80! if any, else 0 81! ipobody(2,i) index referring to the line in field ipobody 82! containing a pointer to the next body load 83! applied to element i, else 0 84! ibody(1,i) code identifying the kind of body load i: 85! 1=centrifugal, 2=gravity, 3=generalized gravity 86! ibody(2,i) amplitude number for load i 87! ibody(3,i) load case number for load i 88! xbody(1,i) size of body load i 89! xbody(2..4,i) for centrifugal loading: point on the axis, 90! for gravity loading with known gravity vector: 91! normalized gravity vector 92! xbody(5..7,i) for centrifugal loading: normalized vector on the 93! rotation axis 94! mi(1) max # of integration points per element (max 95! over all elements) 96! mi(2) max degree of freedom per node (max over all 97! nodes) in fields like v(0:mi(2))... 98! mi(3) max # of layers in any element 99! nacteq(j,i) contains the number of an equation expressed at 100! network node i (i=1..,ntg, ntg is the number of 101! network nodes) 102! j=0: energy conservation 103! j=1: mass conservation 104! j=2: element equation 105! j=3: geometric equation 106! 107! OUTPUT: 108! 109! bc(*) right hand side of the network equation system 110! qat sum of energy contributions in the network 111! nalt number of energy contributions: used to calculate 112! a typical energy flow (used in the convergence 113! criteria) 114! 115 implicit none 116! 117 character*8 lakon(*) 118! 119 integer nelem,ipkon(*),kon(*),nshcon(*),nrhcon(*),ipobody(2,*), 120 & ibody(3,*),mi(*),ntmat_,ielprop(*),ielmat(mi(3),*),node1, 121 & node2,nodem,imat,nalt,ieq,nacteq(0:3,*),index 122! 123 real*8 shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),xbody(7,*),heat, 124 & v(0:mi(2),*),prop(*),xflow,Uout,Tg1,Tg2,Rin,Rout,R1,R2,r,cp,rho, 125 & dvi,Uin,gastemp,cp_cor,bc(*),qat,om 126! 127! 128! 129 heat=0.d0 130! 131 nodem=kon(ipkon(nelem)+2) 132 xflow=v(1,nodem) 133 node1=kon(ipkon(nelem)+1) 134 node2=kon(ipkon(nelem)+3) 135 if((node1.eq.0).or.(node2.eq.0)) return 136! 137 if(lakon(nelem)(2:3).eq.'VO') then 138! 139 index=ielprop(nelem) 140! 141 if(xflow.gt.0d0) then 142 R1=prop(index+2) 143 R2=prop(index+1) 144 Rout=R2 145 Rin=R1 146 else 147 R1=prop(index+2) 148 R2=prop(index+1) 149 Rout=R1 150 Rin=R2 151 endif 152! 153! computing temperature corrected Cp=Cp(T) coefficient 154! 155 Tg1=v(0,node1) 156 Tg2=v(0,node2) 157 if((lakon(nelem)(2:3).ne.'LP').and. 158 & (lakon(nelem)(2:3).ne.'LI')) then 159 gastemp=(Tg1+Tg2)/2.d0 160 else 161 if(xflow.gt.0) then 162 gastemp=Tg1 163 else 164 gastemp=Tg2 165 endif 166 endif 167! 168 imat=ielmat(1,nelem) 169 call materialdata_tg(imat,ntmat_,gastemp, 170 & shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho) 171! 172 call cp_corrected(cp,Tg1,Tg2,cp_cor) 173! 174c Uout=prop(index+5)*Rout 175c Uin=prop(index+5)*Rin 176! 177! free and forced vortices with temperature 178! change in the relative system of coordinates 179! 180 if((lakon(nelem)(2:5).eq.'VOFR') .and. 181 & (nint(prop(index+8)).eq.(-1))) then 182! 183 Uout=prop(index+7)*Rout 184 Uin=prop(index+7)*Rin 185! 186 heat=0.5d0*Cp/Cp_cor*(Uout**2-Uin**2)*xflow 187! 188 elseif (((lakon(nelem)(2:5).eq.'VOFO') 189 & .and.(nint(prop(index+6)).eq.(-1)))) then 190! 191 Uout=prop(index+5)*Rout 192 Uin=prop(index+5)*Rin 193! 194 heat=0.5d0*Cp/Cp_cor*(Uout**2-Uin**2)*xflow 195! 196! forced vortices with temperature change in the absolute system 197! 198 elseif((lakon(nelem)(2:5).eq.'VOFO') 199 & .and.((nint(prop(index+6)).eq.1))) then 200! 201 Uout=prop(index+5)*Rout 202 Uin=prop(index+5)*Rin 203 heat=Cp/Cp_cor*(Uout**2-Uin**2)*xflow 204! 205 endif 206 elseif(lakon(nelem)(2:5).eq.'GAPR') then 207! 208! heat production in a rotating pipe (in the relative 209! system) 210! 211 index=ielprop(nelem) 212! 213 R1=prop(index+8) 214 R2=prop(index+9) 215 om=prop(index+10) 216 if(xflow.gt.0.d0) then 217 Rin=R1 218 Rout=R2 219 else 220 Rin=R2 221 Rout=R1 222 endif 223 Uin=Rin*om 224 Uout=Rout*om 225! 226 heat=(Uout**2-Uin**2)*xflow/2.d0 227! 228 elseif(lakon(nelem)(2:2).eq.'U') then 229! 230! insert here the heat generated in user defined network elements 231! 232! START insert 233! 234 heat=0.d0 235! 236! END insert 237! 238 else 239 return 240 endif 241! 242! including the resulting additional heat flux in the energy equation 243! 244 if(xflow.gt.0d0)then 245 ieq=nacteq(0,node2) 246 if(ieq.ne.0) then 247 bc(ieq)=bc(ieq)+heat 248 qat=qat+dabs(heat) 249 nalt=nalt+1 250 endif 251 else 252 ieq=nacteq(0,node1) 253 if(ieq.ne.0) then 254 bc(ieq)=bc(ieq)+heat 255 qat=qat+dabs(heat) 256 nalt=nalt+1 257 endif 258 endif 259! 260 return 261 end 262 263