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 user_network_element_p1(node1,node2,nodem,nelem,lakon, 20 & kon,ipkon,nactdog,identity,ielprop,prop,iflag,v,xflow,f, 21 & nodef,idirf,df,cp,R,physcon,dvi,numf,set,co,vold,mi,ttime, 22 & time,iaxial,iplausi) 23! 24! UP1 element: mass flow = c * dsqrt(pt2-pt1) 25! 26! f:= c * dsqrt(pt2-pt1) - mass flow 27! 28! for this element it is known that the flow direction 29! has to correspond to the pressure drop 30! 31! INPUT: 32! 33! node1 first node in element topology 34! node2 third node in element topology 35! nodem second node in element topology (middle node) 36! nelem element number 37! lakon(i) label of element i 38! kon connectivity list of all elements; the topology 39! of element i starts at kon(ipkon(i)) 40! ipkon(i) pointer of element i into list kon 41! nactdog(j,i) global degree of freedom in the network 42! equation system of local dof j (0,1 or 2 for 43! networks) of node i. If nactdog(j,i)=0 the 44! variable is known 45! ielprop(i) pointer for element i into field prop. The 46! properties for element i start at 47! prop(ielprop(i)+1,...) 48! prop field of all properties 49! iflag indicates what information should be returned 50! by the routine: 51! 0: identity 52! 1: xflow 53! 2: numf, nodef, idirf, f, df 54! 3: none 55! v(0..mi(2),i) values at node i in the current network 56! iteration (0=total temperature, 57! 1=mass flow, 2=total pressure for network nodes, 58! 0=temperature, 1..3=displacements for structural 59! nodes) 60! cp specific heat at constant pressure corresponding 61! to a mean static temperature across the element 62! R specific gas constant 63! physcon(1..) physical constants (e.g. physcon(1) is absolute 64! zero in the unit systemof the user; cf. the 65! user's manual for the other entries) 66! dvi dynamical viscosity corresponding 67! to a mean static temperature across the element 68! set(i) set name corresponding to set i 69! co(1..3,i) coordinates of node i in the global system 70! vold(0..mi(2),i) values at node i at the start of the current network 71! iterations (0=total temperature, 72! 1=mass flow, 2=total pressure for network nodes, 73! 0=temperature, 1..3=displacements for structural 74! nodes) 75! mi(2) max degree of freedom per node (max over all 76! nodes) in fields like v(0:mi(2))...; the other 77! values of mi are not relevant here 78! ttime total time at the end of the current 79! thermo-mechanical increment. To reach the end 80! of this increment several thermo-mechanical 81! iterations are performed. For each of these 82! iterations a loop of network iterations is 83! performed 84! time step time a the end of the current thermo- 85! mechanical increment 86! iaxial number of times the current structure fits into 87! 360 degrees 88! iplausi flag telling whether any plausibility checks 89! have been violated up to entry in this routine 90! 0: plausibility checks not satisfied 91! 1: plausibility checks (if any) are satisfied 92! 93! 94! OUTPUT: 95! 96! identity if .true. the user_network_element routine is 97! not needed (all variables known) 98! xflow mass flow 99! f value of the element equation 100! nodef nodes corresponding to the variables in the 101! element equation 102! idirf degrees of freedom corresponding to the variables 103! in the element equation 104! df derivatives of the element equation w.r.t. its 105! variables 106! numf number of variables in the element equation 107! iplausi flag telling whether any plausibility checks 108! were violated at return time from this routine 109! 0: plausibility checks not satisfied 110! 1: plausibility checks (if any) are satisfied 111! only feasible change within this routine is from 112! 1 to 0. 113! 114 implicit none 115! 116 logical identity 117 character*8 lakon(*) 118 character*81 set(*) 119! 120 integer nelem,nactdog(0:3,*),node1,node2,nodem,numf, 121 & ielprop(*),nodef(*),idirf(*),iflag,ipkon(*),kon(*), 122 & iaxial,mi(*),inv,index,icase,iplausi 123! 124 real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),R,cp,physcon(*),dvi, 125 & co(3,*),vold(0:mi(2),*),ttime,time,xmach,xflow_oil,Tt1,Tt2, 126 & reynolds,pt1,pt2,kappa,c,a,d,pi,T1,T2 127! 128! 129! 130 pi=4.d0*datan(1.d0) 131 if(iflag.eq.0) then 132! 133! called by envtemp.f: 134! 135! check whether element equation is needed (this is the 136! case if identity=.false. 137! 138 identity=.true. 139! 140 if(nactdog(2,node1).ne.0)then 141 identity=.false. 142 elseif(nactdog(2,node2).ne.0)then 143 identity=.false. 144 elseif(nactdog(1,nodem).ne.0)then 145 identity=.false. 146 endif 147! 148 elseif(iflag.eq.1)then 149 if(v(1,nodem).ne.0.d0) then 150 xflow=v(1,nodem) 151 return 152 endif 153! 154! called by initialnet.f: 155! 156! calculation of the mass flow if everything else is known 157! 158 index=ielprop(nelem) 159 c=prop(index+2) 160! 161 pt1=v(2,node1) 162 pt2=v(2,node2) 163 if(pt1.ge.pt2) then 164 inv=1 165 else 166 inv=-1 167 pt1=v(2,node2) 168 pt2=v(2,node1) 169 endif 170! 171 xflow=c*dsqrt(pt1-pt2) 172! 173 elseif(iflag.eq.2)then 174! 175! called by resultnet.f and mafillnet.f 176! 177 numf=3 178 index=ielprop(nelem) 179! 180 pt1=v(2,node1) 181 pt2=v(2,node2) 182 if(pt1.ge.pt2) then 183! 184! flow direction corresponds to element orientation 185! 186 inv=1 187 xflow=v(1,nodem)*iaxial 188 nodef(1)=node1 189 nodef(2)=nodem 190 nodef(3)=node2 191 else 192! 193! flow direction does not correspond to element orientation 194! 195 inv=-1 196 pt1=v(2,node2) 197 pt2=v(2,node1) 198 xflow=-v(1,nodem)*iaxial 199 nodef(1)=node2 200 nodef(2)=nodem 201 nodef(3)=node1 202 endif 203! 204 idirf(1)=2 205 idirf(2)=1 206 idirf(3)=2 207! 208! nodef and idirf show that the variables are (if inv=1): 209! 1. total pressure in node 1 210! 2. mass flow 211! 3. total pressure in node 2 212! 213 c=prop(index+2) 214! 215 f=c*dsqrt(pt1-pt2)-xflow 216 df(1)=c/(2.d0*dsqrt(pt1-pt2)) 217 df(2)=-1.d0*inv 218 df(3)=-df(1) 219! 220! output 221! 222 elseif(iflag.eq.3) then 223! 224! called by flowoutput.f 225! 226! storage in the .net-file 227! 228 index=ielprop(nelem) 229 a=prop(index+1) 230 kappa=(cp/(cp-R)) 231 icase=1 232! 233 pt1=v(2,node1) 234 pt2=v(2,node2) 235 if(pt1.ge.pt2) then 236 inv=1 237 xflow=v(1,nodem)*iaxial 238 Tt1=v(0,node1)-physcon(1) 239 Tt2=v(0,node2)-physcon(1) 240 else 241 inv=-1 242 pt1=v(2,node2) 243 pt2=v(2,node1) 244 xflow=-v(1,nodem)*iaxial 245 Tt1=v(0,node2)-physcon(1) 246 Tt2=v(0,node1)-physcon(1) 247 endif 248! 249! calculation of the static temperatures 250! 251 if(lakon(nelem)(3:3).eq.'P') then 252! 253! "pipe"-like element: total and static temperatures 254! differ 255! 256 call ts_calc(xflow,Tt1,pt1,kappa,r,A,T1,icase) 257 call ts_calc(xflow,Tt2,pt2,kappa,r,A,T2,icase) 258 else 259! 260! "chamber-connecting"-like element: total and static 261! temperatures are equal 262! 263 T1=Tt1 264 T2=Tt2 265 endif 266! 267! calculation of the dynamic viscosity 268! 269 if(dabs(dvi).lt.1d-30) then 270 write(*,*) '*ERROR in orifice: ' 271 write(*,*) ' no dynamic viscosity defined' 272 write(*,*) ' dvi= ',dvi 273 call exit(201) 274 endif 275! 276 index=ielprop(nelem) 277 kappa=(cp/(cp-R)) 278 a=prop(index+1) 279 d=dsqrt(a*4/Pi) 280 reynolds=dabs(xflow)*d/(dvi*a) 281 xmach=dsqrt(((pt1/pt2)**((kappa-1.d0)/kappa)-1.d0)*2.d0/ 282 & (kappa-1.d0)) 283! 284 xflow_oil=0.d0 285! 286 write(1,*) '' 287 write(1,55) ' from node ',node1, 288 & ' to node ', node2,' : air massflow rate = ' 289 &,inv*xflow,' ', 290 & ', oil massflow rate = ',xflow_oil,' ' 291! 292 if(inv.eq.1) then 293 write(1,56)' Inlet node ',node1,' : Tt1 = ',Tt1, 294 & ' , Ts1 = ',T1,' , Pt1 = ',pt1, ' ' 295! 296 write(1,*)' Element ',nelem,lakon(nelem) 297 write(1,57)' dyn.visc = ',dvi,' , Re = ' 298 & ,reynolds,', M = ',xmach 299! 300 write(1,56)' Outlet node ',node2,' : Tt2 = ',Tt2, 301 & ' , Ts2 = ',T2,' , Pt2 = ',pt2,' ' 302! 303 else if(inv.eq.-1) then 304 write(1,56)' Inlet node ',node2,': Tt1 = ',Tt1, 305 & ' , Ts1 = ',T1,' , Pt1 = ',pt1, ' ' 306 & 307 write(1,*)' element R ',set(numf)(1:30) 308 write(1,57)' dyn.visc. = ',dvi,' , Re =' 309 & ,reynolds,', M = ',xmach 310! 311 write(1,56)' Outlet node ',node1,': Tt2 = ',Tt2, 312 & ' , Ts2 = ',T2,' , Pt2 = ',pt2, ' ' 313 endif 314! 315 endif 316! 317 55 format(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a) 318 56 format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a) 319 57 format(1x,a,e11.4,a,e11.4,a,e11.4) 320! 321! degree of freedom 2 is the mass flow dof (cf. field idirf) 322! 323 xflow=xflow/iaxial 324 df(2)=df(2)*iaxial 325! 326 return 327 end 328