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 real*8 function calc_residual_tee(pt1,Tt1,xflow1,xflow2, 19 &pt2,Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta) 20! 21! Function for calculating the residual of both branches of a tee 22! 23! author: Yannick Muller 24! 25 implicit none 26! 27 integer 28! Isothermal/adiabatic 29 &icase, 30! Residual/Derivatives 31 &ider, 32! Critical conditions at inlet or outlet 33 &icrit1, 34 &icrit2, 35! Where we are in the program 36 &iflag 37! 38 real*8 39! Residual 40 &f, 41! Kappa stuff 42 &kappa, 43 &R, 44! State variables 45 &pt1, 46 &pt2, 47 &Tt1, 48 &Tt2, 49 &Ts0, 50 &Ts1, 51 &Ts2, 52 &xflow1, 53 &xflow2, 54! 55 &pt2_lim, 56! 57 &zeta, 58 &zeta_fac, 59! 60! Areas 61 &A1, 62 &A2, 63! Reduced mass flows 64 &Q0, 65 &Q1, 66 &Q2, 67 &Q_crit, 68! Pressure ratios 69 &pspt_crit, 70 &pspt0, 71 &pspt1, 72 &pspt2, 73! Flow velocities 74 &w1, 75 &w2, 76 &w1w2, 77 &w2w1, 78! Mach numbers, 79 &M1, 80 &M2 81! 82 icrit1 = 0 83 icrit2 = 0 84! setting icase (always adiabatic) 85 icase=0; 86! 87! Critical values 88 Q_crit = dsqrt(kappa/R)* 89 & (1+0.5d0*(kappa-1))**(-0.5d0*(kappa+1)/(kappa-1)) 90 pspt_crit = (2/(kappa+1)) ** (kappa/(kappa-1)) 91! 92! These reduced mass flows are equivalent 93! (reduced mass flow at inlet) 94 Q0 = xflow1*dsqrt(Tt1)/pt1/A1 95 Q1 = xflow2*dsqrt(Tt1)/pt1/A2 96 if(Q1.ge.Q_crit) then 97 Q1 = Q_crit 98 icrit1 = 1 99 write(*,*)'*WARNING in Tee:' 100 write(*,*)'Critical conditions at 1' 101 endif 102! Reduced mass flow at outlet 103 Q2 = xflow2*dsqrt(Tt1)/pt2/A2 104 if(Q2.ge.Q_crit) then 105 Q2 = Q_crit 106 icrit2 = 1 107 write(*,*)'*WARNING in Tee:' 108 write(*,*)'Critical conditions at 2' 109 endif 110! 111! Flow velocity at inlet 112! Static temperature 113 Ts0=Tt1 114 call ts_calc(xflow1,Tt1,pt1,kappa,r,A1,Ts0,icase) 115! Pressure ratio 116 pspt0 = (Ts0/Tt1)**(kappa/(kappa-1)) 117! Velocity 118 call wpi(w1, pspt0, Q0, 119 & dsqrt(Tt1),kappa,R) 120! 121! Flow velocity at outlet 122! Static temperature 123 call ts_calc(xflow2,Tt1,pt1,kappa,r,A2,Ts1,icase) 124! Pressure ratio 125 pspt1 = (Ts1/Tt1)**(kappa/(kappa-1)) 126! Velocity 127 call wpi(w2, pspt1, Q1, 128 & dsqrt(Tt1),kappa,R) 129! 130! Velocity ratio 131c w2w1=w2/w1 132c w1w2=w1/w2 133 if(w2.eq.0.d0)then 134 w1w2=1d30 135 else 136 w1w2=w1/w2 137 endif 138 if(w1.eq.0.d0)then 139 w2w1=1d30 140 else 141 w2w1=w2/w1 142 endif 143! 144! Zeta calculation 145 zeta=1.d0+0.3d0*W2W1**2 146 zeta=zeta*(W1W2)**2 147! 148 zeta = zeta_fac*zeta 149! 150! Residual calculation 151 if(icrit2.ne.1) then 152 if(icrit1.ne.1) then 153 f = pt2 - pt1*pspt1**zeta 154 else 155 f = xflow2*dsqrt(Tt1)/pt1/A2-Q_crit 156 endif 157 else 158 f = xflow2*dsqrt(Tt1)/pt2/A2-Q_crit 159 endif 160! 161 if(iflag.eq.3) then 162! 163 write(1,57)' zeta= ',zeta 164 57 format(1x,a,f9.4) 165! 166 else if (iflag.eq.4) then 167 168! Calculate Mach numbers 169 call machpi(M1,pspt0,kappa,R) 170 call ts_calc(xflow2,Tt2,pt2,kappa,r,A2,Ts2,icase) 171! Pressure ratio 172 pspt2 = (Ts2/Tt2)**(kappa/(kappa-1)) 173 call machpi(M2,pspt2,kappa,R) 174 175! write(1,80)'Inlet: Tt1= ',Tt1, 176! & ', pt1= ',pt1,', M1= ',M1 177! 178! write(1,77)'mass flow = ',xflow2,', kappa = ',kappa, 179! & ', zeta= ',zeta 180! 181! write(1,80)'Outlet: Tt2= ',Tt2, 182! & ', pt2= ',pt2,', M2= ',M2 183! 184! 80 format(3x,a,f10.6,a,f10.2,a,f10.6) 185! 77 format(3x,a,f10.6,a,f10.2,a,f10.6) 186 endif 187! 188 calc_residual_tee=f 189! 190 return 191 end 192