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 film(h,sink,temp,kstep,kinc,time,noel,npt, 20 & coords,jltyp,field,nfield,loadtype,node,area,vold,mi, 21 & ipkon,kon,lakon,iponoel,inoel,ielprop,prop,ielmat, 22 & shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon, 23 & ipobody,xbody,ibody,heatnod,heatfac) 24! 25! user subroutine film 26! 27! 28! INPUT: 29! 30! sink most recent sink temperature 31! temp current temperature value 32! kstep step number 33! kinc increment number 34! time(1) current step time 35! time(2) current total time 36! noel element number 37! npt integration point number 38! coords(1..3) global coordinates of the integration point 39! jltyp loading face kode: 40! 11 = face 1 41! 12 = face 2 42! 13 = face 3 43! 14 = face 4 44! 15 = face 5 45! 16 = face 6 46! field currently not used 47! nfield currently not used (value = 1) 48! loadtype load type label 49! node network node (only for forced convection) 50! area area covered by the integration point 51! vold(0..4,1..nk) actual solution field in all nodes; 52! for structural nodes: 53! 0: temperature 54! 1: displacement in global x-direction 55! 2: displacement in global y-direction 56! 3: displacement in global z-direction 57! 4: static pressure 58! for network nodes 59! 0: total temperature (at end nodes) 60! = static temperature for liquids 61! 1: mass flow (at middle nodes) 62! 2: total pressure (at end nodes) 63! = static pressure for liquids 64! 3: static temperature (at end nodes; only for gas) 65! mi(1) max # of integration points per element (max 66! over all elements) 67! mi(2) max degree of freedom per node (max over all 68! nodes) in fields like v(0:mi(2))... 69! mi(3) max # of layers in any element 70! ipkon(i) points to the location in field kon preceding 71! the topology of element i 72! kon(*) contains the topology of all elements. The 73! topology of element i starts at kon(ipkon(i)+1) 74! and continues until all nodes are covered. The 75! number of nodes depends on the element label 76! lakon(i) contains the label of element i 77! iponoel(i) the network elements to which node i belongs 78! are stored in inoel(1,iponoel(i)), 79! inoel(1,inoel(2,iponoel(i)))...... until 80! inoel(2,inoel(2,inoel(2......)=0 81! inoel(1..2,*) field containing the network elements 82! ielprop(i) points to the location in field prop preceding 83! the properties of element i 84! prop(*) contains the properties of all network elements. The 85! properties of element i start at prop(ielprop(i)+1) 86! and continues until all properties are covered. The 87! appropriate amount of properties depends on the 88! element label. The kind of properties, their 89! number and their order corresponds 90! to the description in the user's manual, 91! cf. the sections "Fluid Section Types" 92! ielmat(j,i) contains the material number for element i 93! and layer j 94! shcon(0,j,i) temperature at temperature point j of material i 95! shcon(1,j,i) specific heat at constant pressure at the 96! temperature point j of material i 97! shcon(2,j,i) dynamic viscosity at the temperature point j of 98! material i 99! shcon(3,1,i) specific gas constant of material i 100! nshcon(i) number of temperature data points for the specific 101! heat of material i 102! rhcon(0,j,i) temperature at density temperature point j of 103! material i 104! rhcon(1,j,i) density at the density temperature point j of 105! material i 106! nrhcon(i) number of temperature data points for the density 107! of material i 108! ntmat_ maximum number of temperature data points for 109! any material property for any material 110! ncocon(1,i) number of conductivity constants for material i 111! ncocon(2,i) number of temperature data points for the 112! conductivity coefficients of material i 113! cocon(0,j,i) temperature at conductivity temperature point 114! j of material i 115! cocon(k,j,i) conductivity coefficient k at conductivity 116! temperature point j of material i 117! ipobody(1,i) points to an entry in fields ibody and xbody 118! containing the body load applied to element i, 119! if any, else 0 120! ipobody(2,i) index referring to the line in field ipobody 121! containing a pointer to the next body load 122! applied to element i, else 0 123! ibody(1,i) code identifying the kind of body load i: 124! 1=centrifugal, 2=gravity, 3=generalized gravity 125! ibody(2,i) amplitude number for load i 126! ibody(3,i) load case number for load i 127! xbody(1,i) size of body load i 128! xbody(2..4,i) for centrifugal loading: point on the axis, 129! for gravity loading with known gravity vector: 130! normalized gravity vector 131! xbody(5..7,i) for centrifugal loading: normalized vector on the 132! rotation axis 133! 134! OUTPUT: 135! 136! h(1) magnitude of the film coefficient 137! h(2) not used; please do NOT assign any value 138! sink (updated) sink temperature (need not be 139! defined for forced convection) 140! heatnod extra heat flow going to the network node 141! (zero if not specified) 142! heatfac extra heat flow going to the structural face 143! (zero if not specified) 144! 145 implicit none 146! 147 character*8 lakon(*) 148 character*20 loadtype 149! 150 integer kstep,kinc,noel,npt,jltyp,nfield,node,mi(*),ipkon(*), 151 & kon(*),iponoel(*),inoel(2,*),ielprop(*),ielmat(mi(3),*),ntmat_, 152 & nshcon(*),nrhcon(*),ncocon(2,*),nodem,indexprop,indexe, 153 & iel1,iel2,ielup,iit,imat,icase,itherm,ipobody(2,*), 154 & ibody(3,*) 155! 156 real*8 h(2),sink,time(2),coords(3),temp,field(nfield),area, 157 & vold(0:mi(2),*),prop(*),shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*), 158 & cocon(0:6,ntmat_,*),rho,r,pt,Pr,xl,Tt,Ts,xflow,Tsold,Re,um, 159 & xks,xkappa,xlambda,f,cp,A,D,form_fact,xbody(7,*),heatnod, 160 & heatfac 161! 162! 163! 164! 165 if((node.eq.0).or.(iponoel(node).eq.0)) then 166! 167! simple example: constant film coefficient 168! 169 h(1)=200.d0 170 else 171! 172! complicated example: forced convection with a film coefficient 173! determined by the Gnielinski equation (pipe application for 174! turbulent flow), 175! which runs like: 176! 177! NuD=f/8*(ReD-1000)*pr/(1+12.7*sqrt(f/8)*(Pr**(2/3)-1)) 178! 179! NuD=h*D/xlambda: Nusselt number 180! ReD=massflow*D/(um*A): Reynolds number 181! Pr=um*cp/xlambda: Prandl number 182! h: film coefficient 183! D: hydraulic diameter 184! A: area of the pipe cross section 185! f: friction coefficient (cf. User's manual for formula) 186! xlambda: heat conduction coefficient 187! um: dynamic viscosity 188! cp: heat capacity 189! 190! The convection node comes from the calling program: "node" 191! 192! The elements connected to this node are determined by: 193! iel1=inoel(1,iponoel(node)), 194! iel2=inoel(1,inoel(2,iponoel(node))), 195! .... until 196! inoel(2,inoel(2,inoel(2,.......inoel(2,iponoel(node)))))=0 197! 198! Let us assumed that "node" belongs to exactly two network 199! elements, i.e. inoel(2,inoel(2,iponoel(node)))=0 200! 201! Let us look at the upstream element. To determine the 202! upstream element one looks at the sign of the mass flow 203! in the middle node of the elements. 204! 205! The middle node nodem of element iel1 is given by 206! kon(ipkon(iel1)+2), the end nodes by 207! kon(ipkon(iel1)+1) and kon(ipkon(iel1)+3). The mass 208! flow in the element is given by vold(1,nodem). If this value 209! is positive and kon(ipkon(iel1)+3)=node or if is negative 210! and kon(ipkon(iel1)+1)=node iel1 is the upstream element 211! 212! Let us assume this element is a gas pipe element, i.e. 213! lakon(iel1)='DGAPFA '. The properties of such an element 214! are the cross sectional area, the hydraulic diameter.... (in 215! that order), i.e. prop(ielprop(iel1)+1)=A, 216! prop(ielprop(iel1)+2)=D,... 217! 218! The material number of element iel1 is given by ielmat(iel1) 219! with this information and the gas temperature (sink=Tfluid) the 220! gas material properties can be determined: xlambda by calling 221! materialdata_cond, um by calling materialdata_dvi and cp by 222! calling materialdata_cp 223! 224! This yields all data needed to determine the film coefficient 225! with the Gnielinski equation 226! 227! For laminar flow (Re < 3000) the following laminar flow 228! equation is used: 229! 230! h=4.36*xlambda/D 231! 232! elements belonging to the network node 233! 234 iel1=inoel(1,iponoel(node)) 235 if(inoel(2,iponoel(node)).ne.0) then 236 iel2=inoel(1,inoel(2,iponoel(node))) 237! 238! check whether the node belongs to maximum two elements 239! 240 if(inoel(2,inoel(2,iponoel(node))).ne.0) then 241 write(*,*) 'ERROR in film: the network node' 242 write(*,*) ' belongs to more than 2 elements' 243 call exit(201) 244 endif 245 else 246! 247! node (must be an end node) belongs only to one 248! element (pure thermal network without entry or 249! exit element) 250! 251 iel2=iel1 252 endif 253! 254! only adiabatic gas pipes considered (heat exchange only 255! takes place at the end nodes, not within the pipe) 256! 257 icase=0 258! 259! check whether iel1 is the upstream element 260! 261 indexe=ipkon(iel1) 262 nodem=kon(indexe+2) 263 if(((vold(1,nodem).ge.0.d0).and.(kon(indexe+3).eq.node)).or. 264 & ((vold(1,nodem).le.0.d0).and.(kon(indexe+1).eq.node))) then 265 ielup=iel1 266 else 267 ielup=iel2 268 endif 269! 270! check whether the upstream element is an adiabatic gas 271! pipe element or a White-Colebrook liquid pipe 272! 273 if((lakon(ielup)(1:6).ne.'DGAPFA').and. 274 & (lakon(ielup)(1:7).ne.'DLIPIWC')) then 275 write(*,*) 'ERROR in film: upstream element',ielup 276 write(*,*) ' is no adiabatic' 277 write(*,*) ' gas pipe nor White-Colebrook' 278 write(*,*) ' liquid pipe' 279 call exit(201) 280 endif 281! 282 indexprop=ielprop(ielup) 283 A=prop(indexprop+1) 284 D=prop(indexprop+2) 285 xl=prop(indexprop+3) 286 xks=prop(indexprop+4) 287 form_fact=prop(indexprop+5) 288! 289! material of upstream element 290! 291 imat=ielmat(1,ielup) 292! 293 xflow=dabs(vold(1,nodem)) 294 Tt=vold(0,node) 295 pt=vold(2,node) 296! 297 if(lakon(ielup)(2:2).eq.'G') then 298! 299! compressible flow 300! 301 iit=0 302 Ts=Tt 303 do 304 iit=iit+1 305 Tsold=Ts 306 call materialdata_tg(imat,ntmat_,Ts,shcon,nshcon,cp,r, 307 & um,rhcon,nrhcon,rho) 308 call ts_calc(xflow,Tt,pt,xkappa,r,A,Ts,icase) 309 if((dabs(Ts-Tsold).le.1.d-5*Ts).or.(iit.eq.10)) exit 310 enddo 311! 312 call materialdata_tg(imat,ntmat_,Ts,shcon,nshcon,cp,r, 313 & um,rhcon,nrhcon,rho) 314 else 315! 316! incompressible flow 317! 318 Ts=Tt 319 call materialdata_cp(imat,ntmat_,Ts,shcon,nshcon,cp) 320! 321 itherm=2 322 call materialdata_dvi(shcon,nshcon,imat,um,Ts,ntmat_, 323 & itherm) 324! 325 endif 326! 327 call materialdata_cond(imat,ntmat_,Ts,cocon,ncocon,xlambda) 328! 329 Re=xflow*D/(um*A) 330! 331 if(Re.lt.3000.d0) then 332! 333! laminar flow 334! 335 h(1)=4.36d0*xlambda/D 336 else 337! 338! turbulent flow 339! 340 if(Re.gt.5.e6) then 341 write(*,*) '*ERROR in film: Reynolds number ',Re 342 write(*,*) ' is outside valid range' 343 call exit(201) 344 endif 345! 346 Pr=um*cp/xlambda 347! 348 if((Pr.lt.0.5d0).or.(Pr.gt.2000.d0)) then 349 write(*,*) '*ERROR in film: Prandl number ',Pr 350 write(*,*) ' is outside valid range' 351 call exit(201) 352 endif 353! 354 call friction_coefficient(xl,D,xks,Re,form_fact,f) 355! 356 h(1)=f/8.d0*(Re-1000.d0)*Pr/ 357 & (1.d0+12.7d0*dsqrt(f/8.d0)*(Pr**(2.d0/3.d0)-1.d0)) 358 & *xlambda/D 359 endif 360 endif 361! 362 return 363 end 364 365