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! author: Yannick Muller 19! 20 real*8 function calc_residual_wye(pt1,Tt1,xflow1,xflow2, 21 & pt2,Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa, 22 & R,ider,iflag,zeta) 23! 24 implicit none 25! 26 integer icase,ichan_num,icrit1,icrit2,ider,iflag 27! 28 real*8 29! In- and Output 30 & f, 31 & R, 32! 33! Kappa stuff 34 & kappa, 35! 36 & pt1, 37 & pt2, 38 & Tt1, 39 & Tt2, 40 & xflow1, 41 & xflow2, 42! 43 & zeta, 44! 45 & A1, 46 & A2, 47 & A_s, 48! 49 & Ts1, 50 & Ts2, 51 & Ts_s, 52 & dh1, 53 & dh2, 54 & alpha, 55 & Q_crit, 56 & pspt_crit, 57 & Q0, 58 & Q1, 59 & Q2, 60 & pspt0, 61 & pspt1, 62 & w1, 63 & w2, 64 & w1w2, 65 & w2w1, 66 & VsV1, 67 & pi, 68 & z2d390, 69 & z1p090, 70 & z60, 71 & z90, 72 & hq, 73 & Ts0, 74 & zeta_fac, 75 & M1, 76 & M2, 77 & pspt2 78! 79 real*8 Table_A(2,11) 80! 81 icrit1 = 0 82 icrit2 = 0 83! 84! setting icase (always adiabatic) 85 icase=0; 86! 87 pi=4.d0*datan(1.d0) 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 Q0 = xflow1*dsqrt(Tt1)/pt1/A1 93 Q1 = xflow2*dsqrt(Tt1)/pt1/A2 94 if(Q1.ge.Q_crit) then 95 Q1 = Q_crit 96 icrit1 = 1 97 write(*,*)'*WARNING in Wye:' 98 write(*,*)'Critical conditions at 1' 99 endif 100 Q2 = xflow2*dsqrt(Tt1)/pt2/A2 101 if(Q2.ge.Q_crit) then 102 Q2 = Q_crit 103 icrit2 = 1 104 write(*,*)'*WARNING in Wye:' 105 write(*,*)'Critical conditions at 2' 106 endif 107! 108! Flow velocity at inlet 109 Ts0=Tt1 110 111 call ts_calc(xflow1,Tt1,pt1,kappa,r,A1,Ts0,icase) 112 pspt0 = (Ts0/Tt1)**(kappa/(kappa-1)) 113 call wpi(w1, pspt0, Q0, 114 & dsqrt(Tt1),kappa,R) 115! 116! Flow velocity at outlet 117 call ts_calc(xflow2,Tt1,pt1,kappa,r,A2,Ts1,icase) 118 pspt1 = (Ts1/Tt1)**(kappa/(kappa-1)) 119 call wpi(w2, pspt1, Q1, 120 & dsqrt(Tt1),kappa,R) 121! 122 if(w2.eq.0.d0)then 123 w1w2=1d30 124 else 125 w1w2=w1/w2 126 endif 127 if(w1.eq.0.d0)then 128 w2w1=1d30 129 else 130 w2w1=w2/w1 131 endif 132! 133! Zeta calculation 134! Main branch 135 if(ichan_num.eq.1) then 136! Zeta as in Calculix and old ACC tool 137 zeta=0.4d0*(1-W2W1)**2 138 zeta=zeta*(W1W2)**2 139! 140! Branch 141 elseif(ichan_num.eq.2) then 142 hq=dh2/dh1 143! 144! Interpolation as in CalculiX 145 if((alpha.le.60.and.hq.le.1.d0).or. 146 & (alpha.le.90.d0.and.hq.le.2.d0/3.d0)) then 147 zeta=0.95d0*((W2W1-2d0*dcos(alpha*pi/180)) 148 & *W2W1+1.d0) 149 zeta=zeta*(W1W2)**2 150 elseif(alpha.le.90.d0.and.hq.le.1.d0) then 151 z2d390=0.95d0*((W2W1-2d0*dcos(90.d0*pi/180)) 152 & *W2W1+1.d0) 153! 154 z1p090=0.95d0*(1.d0+0.3d0*W2W1**2) 155! 156 z90=z2d390+(3*hq-2.d0)*(z1p090-z2d390) 157! 158 Z60=0.95d0*((W2W1-2d0*dcos(60.d0*pi/180)) 159 & *W2W1+1.d0) 160! 161 zeta=z60+(alpha/30.d0-2.d0)*(z90-z60) 162 zeta=zeta*(W1W2)**2 163 elseif(alpha.le.60.d0.and.hq.gt.1.d0) then 164! 165! extrapolation for hq>1 (not part of Idelchik 166! but the previous definition results sometimes 167! into zeta<0 values 168! 169 write(*,*)'WARNING in calc_residual_wye:' 170 write(*,*)' Branch element is 171 &outside valid range defined by Idelchik' 172 zeta=0.95d0*((W2W1-2d0*dcos(alpha*pi/180)) 173 & *W2W1+1.d0) 174 zeta=zeta*(W1W2)**2 175 elseif(alpha.le.90.d0.and.hq.gt.1.d0) then 176! 177! extrapolation for hq>1 (not part of Idelchik 178! but the previous definition results sometimes 179! into zeta<0 values 180! 181 write(*,*)'WARNING in calc_residual_wye:' 182 write(*,*)' Branch element is 183 &outside valid range defined by Idelchik' 184 Z60=0.95d0*((W2W1-2d0*dcos(60.d0*pi/180)) 185 & *W2W1+1.d0) 186! 187 z1p090=0.95d0*(1.d0+0.3d0*W2W1**2) 188! 189 zeta=Z60+(z1p090-Z60)/(90.d0-60.d0)*(alpha-60) 190! 191 zeta=zeta*(W1W2)**2 192! interpolation between exit loss zeta=1 193! and last point of definition range (hq=1) 194 zeta=zeta+(1-zeta)/(50.d0-1.d0)*(hq-1) 195 else 196! 197! values are absolutely out of def. range 198! 199 write(*,*)'ERROR in wye.f:' 200 write(*,*)' Branch element is 201 &outside valid range' 202 call exit(201) 203 endif 204 endif 205! 206 zeta = zeta*zeta_fac 207! 208 if(icrit2.ne.1) then 209 if(icrit1.ne.1) then 210 f = pt2 - pt1*pspt1**zeta 211 else 212 f = xflow2*dsqrt(Tt1)/pt1/A2-Q_crit 213 endif 214 else 215 f = xflow2*dsqrt(Tt1)/pt2/A2-Q_crit 216 endif 217! 218 calc_residual_wye=f 219 220 if(iflag.eq.4) then 221! Calculate Mach numbers 222 call machpi(M1,pspt0,kappa,R) 223 call ts_calc(xflow2,Tt2,pt2,kappa,r,A2,Ts2,icase) 224! Pressure ratio 225 pspt2 = (Ts2/Tt2)**(kappa/(kappa-1)) 226 call machpi(M2,pspt2,kappa,R) 227! 228 80 format(3x,a,f10.6,a,f10.2,a,f10.6) 229 77 format(3x,a,f10.6,a,f10.2,a,f10.6) 230 endif 231! 232 return 233 end 234