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! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12! GNU General Public License for more details. 13! 14! You should have received a copy of the GNU General Public License 15! along with this program; if not, write to the Free Software 16! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 17! 18 subroutine wye(node1,node2,nodem,nelem,lakon,kon,ipkon, 19 & nactdog,identity,ielprop,prop,iflag,v,xflow,f, 20 & nodef,idirf,df,cp,r,physcon,numf,set,mi,ider,ttime,time, 21 & iaxial,iplausi,dvi) 22! 23! A wye split element(zeta calculation according to Idel'chik) 24! Written by Yavor Dobrev 25! For an explanation of the parameters see tee.f 26! 27! author: Yannick Muller 28! 29 implicit none 30! 31 logical identity 32! 33 character*81 set(*) 34 character*8 lakon(*) 35! 36 integer nelem,nactdog(0:3,*),node1,node2,nodem,nodem1,numf,kon(*), 37 &ipkon(*),ielprop(*),nodef(*),idirf(*),index,iflag,inv,mi(2), 38 &nelem1,ichan_num,ider,icase,i,iaxial,iplausi 39! 40 real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,R,Tt1,Tt2,pt1,pt2, 41 &cp,physcon(*),km1,kp1,kdkm1,pt2pt1,pt1pt2,pt2pt1_crit,tdkp1, 42 &A,A1,A2,A_s,calc_residual_wye,dh1,dh2,alpha,xflow1,xflow2, 43 &pi,zeta_fac,Ts0,pspt0,pspt2,M1,M2,Ts2,ttime,time,zeta,dvi 44! 45! 46! 47 index=ielprop(nelem) 48! 49 if (iflag.eq.0.d0) then 50 identity=.true. 51 if(nactdog(2,node1).ne.0)then 52 identity=.false. 53 elseif(nactdog(2,node2).ne.0)then 54 identity=.false. 55 elseif(nactdog(1,nodem).ne.0)then 56 identity=.false. 57 endif 58! 59 elseif (iflag.eq.1)then 60 if(v(1,nodem).ne.0.d0) then 61 xflow=v(1,nodem) 62 return 63 endif 64! 65 kappa=(cp/(cp-R)) 66 kp1=kappa+1d0 67 km1=kappa-1d0 68 kdkm1=kappa/km1 69 tdkp1=2.d0/kp1 70! 71 if(nelem.eq.nint(prop(index+2))) then 72 A=prop(index+4) 73 elseif(nelem.eq.nint(prop(index+3)))then 74 A=prop(index+6) 75 endif 76! 77 pt1=v(2,node1) 78 pt2=v(2,node2) 79! 80 if(pt1.ge.pt2) then 81 inv=1 82 Tt1=v(0,node1)-physcon(1) 83 Tt2=v(0,node2)-physcon(1) 84 else 85 inv=-1 86 pt1=v(2,node2) 87 pt2=v(2,node1) 88 Tt1=v(0,node2)-physcon(1) 89 Tt2=v(0,node1)-physcon(1) 90 endif 91! 92 pt1pt2=pt1/pt2 93 pt2pt1=1/pt1pt2 94! 95 pt2pt1_crit=tdkp1**kdkm1 96! 97 if(pt2pt1.gt.pt2pt1_crit) then 98 xflow=inv*pt1*A*dsqrt(2.d0*kdkm1*pt2pt1**(2.d0/kappa) 99 & *(1.d0-pt2pt1**(1.d0/kdkm1))/r)/dsqrt(Tt1) 100 else 101 xflow=inv*pt1*A*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/ 102 & dsqrt(Tt1) 103 endif 104 xflow=xflow/50 105! 106 elseif (iflag.eq.2)then 107! 108 numf=6 109! 110 kappa=(cp/(cp-R)) 111 pi=4.d0*datan(1.d0) 112! 113! Inlet conditions are the same for both branches 114! 115! Determining the previous element as 116! the incoming mass flow is defined there 117! 118 nelem1=nint(prop(index+1)) 119 nodem1=kon(ipkon(nelem1)+2) 120! 121! Inlet conditions 122! 123 pt1=v(2,node1) 124 Tt1=v(0,node1)-physcon(1) 125 xflow1=v(1,nodem1)*iaxial 126 A1 = prop(index+4) 127! 128! Outlet conditions 129! 130 Tt2=v(0,node2) 131 xflow2=v(1,nodem)*iaxial 132 pt2=v(2,node2) 133! 134 if(nelem.eq.nint(prop(index+2))) then 135! 136 A2 = A1 137 A_s = prop(index+6) 138! 139 dh1 = prop(index+9) 140 if(dh1.eq.0.d0) then 141 dh1 = dsqrt(4*A1/pi) 142 endif 143! 144 dh2 = dh1 145 alpha = 0 146 ichan_num = 1 147 zeta_fac = prop(index+13) 148! 149 elseif(nelem.eq.nint(prop(index+3))) then 150! 151 ichan_num = 2 152 A2 = prop(index+6) 153! 154 dh1 = prop(index+9) 155 if(dh1.eq.0.d0) then 156 dh1 = dsqrt(4*A1/pi) 157 endif 158! 159 dh2 = prop(index+10) 160 if(dh2.eq.0.d0) then 161 dh2 = dsqrt(4*A2/pi) 162 endif 163! 164 alpha = prop(index+8) 165 zeta_fac = prop(index+14) 166! 167 endif 168! 169! Set the node numbers for the degrees of freedom 170 nodef(1)=node1 171 nodef(2)=node1 172 nodef(3)=nodem1 173 nodef(4)=nodem 174 nodef(5)=node2 175 nodef(6)=node2 176! 177! Sets the types of the degrees of freedom 178 idirf(1)=2 179 idirf(2)=0 180 idirf(3)=1 181 idirf(4)=1 182 idirf(5)=2 183 idirf(6)=0 184! 185 if(ider.eq.0.d0) then 186! Residual 187 f=calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2, 188 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,iflag 189 &,zeta) 190 else 191! Derivatives 192 call calc_ider_wye(df,pt1,Tt1,xflow1,xflow2,pt2, 193 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,iflag 194 &,zeta) 195 endif 196! 197 elseif(iflag.eq.3) then 198! 199 kappa=(cp/(cp-R)) 200! setting icase (always adiabatic) 201 icase=0; 202! 203 pi=4.d0*datan(1.d0) 204! 205! Inlet conditions are the same for both branches 206! 207! Determining the previous element as 208! the incoming mass flow is defined there 209 nelem1=nint(prop(index+1)) 210 nodem1=kon(ipkon(nelem1)+2) 211! 212! Inlet conditions 213 pt1=v(2,node1) 214 Tt1=v(0,node1)-physcon(1) 215 xflow1=v(1,nodem1)*iaxial 216 A1 = prop(index+4) 217! 218! Outlet conditions 219 Tt2=v(0,node2) 220 xflow2=v(1,nodem)*iaxial 221 pt2=v(2,node2) 222! 223 if(nelem.eq.nint(prop(index+2))) then 224! 225 A2 = A1 226 A_s = prop(index+6) 227! 228 dh1 = prop(index+9) 229 if(dh1.eq.0.d0) then 230 dh1 = dsqrt(4*A1/pi) 231 endif 232! 233 dh2 = dh1 234 alpha = 0 235 ichan_num = 1 236 zeta_fac = prop(index+13) 237! 238 elseif(nelem.eq.nint(prop(index+3))) then 239! 240 ichan_num = 2 241 A2 = prop(index+6) 242! 243 dh1 = prop(index+9) 244 if(dh1.eq.0.d0) then 245 dh1 = dsqrt(4*A1/pi) 246 endif 247! 248 dh2 = prop(index+10) 249 if(dh2.eq.0.d0) then 250 dh2 = dsqrt(4*A2/pi) 251 endif 252! 253 alpha = prop(index+8) 254 zeta_fac = prop(index+14) 255! 256 endif 257! 258! Write the main information about the element 259! 260! Flow velocity at inlet 261 call ts_calc(xflow1,Tt1,pt1,kappa,r,A1,Ts0,icase) 262 pspt0 = (Ts0/Tt1)**(kappa/(kappa-1)) 263! Calculate Mach numbers 264 call machpi(M1,pspt0,kappa,R) 265 call ts_calc(xflow2,Tt2,pt2,kappa,r,A2,Ts2,icase) 266! Pressure ratio 267 pspt2 = (Ts2/Tt2)**(kappa/(kappa-1)) 268 call machpi(M2,pspt2,kappa,R) 269! 270 write(1,*) '' 271 write(1,55) ' from node ',node1, 272 & ' to node ', node2,': air massflow rate= ',xflow 273! 274 write(1,56)' Inlet node ',node1,': Tt1= ',Tt1, 275 & ' , Ts1= ',Ts0,' , Pt1= ',pt1, 276 & ', M1= ',M1 277 write(1,*)' Element ',nelem,lakon(nelem) 278 & ,', Branch ',ichan_num 279! 280 55 format(1x,a,i6,a,i6,a,e11.4,a) 281 56 format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4) 282! 283! Set ider to calculate the residual 284 ider = 0 285! 286! Calculate the element one last time with enabled output 287 f=calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2, 288 & Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac, 289 & kappa,R,ider,iflag,zeta) 290! 291 write(1,56)' Outlet node ',node2,': Tt2= ',Tt2, 292 & ' , Ts2= ',Ts2,' , Pt2= ',pt2, 293 & ', M2= ',M2 294! 295 endif 296! 297 xflow=xflow/iaxial 298 df(3)=df(3)*iaxial 299 df(4)=df(4)*iaxial 300! 301 return 302 end 303