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! author: Yannick Muller 20! 21 real*8 function calc_residual_cross_split(pt1,Tt1,xflow1,xflow2, 22 &pt2,Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac, 23 &kappa,R,ider,iflag) 24! 25 implicit none 26! 27 integer icase,ichan_num,ider,icrit1,icrit2,iflag,ier 28! 29 real*8 30! In- and Output 31 &f,R, 32! 33! Kappa stuff 34 &kappa,km1,kp1, 35! 36 &pt1,pt2,Tt1,Tt2,xflow1,xflow2, 37! 38 &pt2_lim, 39! 40 &zeta, 41! 42 &A1,A2, 43! 44 &Ts0,Ts1,Ts2,dh1,dh2,alpha,Q_crit,pspt_crit,Q0,Q1,Q2,pspt0, 45 &pspt1,pspt2,w1,w2,w1w2,w2w1,pi,z2d390,z1p090,z60,z90,hq,M1,M2, 46 &zeta_fac,xflow_s,Q_s,Ts_s,pspt_s,w_s,wsw1,A_s,AsA1,VsV1 47! 48 real*8 Table_zeta(2,10) 49! 50 pi=4.d0*datan(1.d0) 51! 52 icrit1=0 53 icrit2=0 54! 55! setting icase (always adiabatic) 56! 57 icase=0; 58! 59 km1=kappa-1.d0 60 kp1=kappa+1.d0 61 Q_crit=dsqrt(kappa/R)* 62 & (1+0.5d0*(kappa-1))**(-0.5d0*(kappa+1)/(kappa-1)) 63 pspt_crit=(2.d0/(KAPPA+1.d0))**(KAPPA/(KAPPA-1.d0)) 64! 65 Q0=xflow1*dsqrt(Tt1)/pt1/A1 66 Q1=xflow2*dsqrt(Tt1)/pt1/A2 67 if(Q1.ge.Q_crit) then 68 Q1=Q_crit 69 icrit1=1 70 write(*,*)'*WARNING in Cross Split:' 71 write(*,*)'Critical conditions at 1' 72 endif 73 Q2=xflow2*dsqrt(Tt1)/pt2/A2 74 if(Q2.ge.Q_crit) then 75 Q2=Q_crit 76 icrit2=1 77 write(*,*)'*WARNING in Cross Split:' 78 write(*,*)'Critical conditions at 2' 79 endif 80! 81! Flow velocity at inlet 82 call ts_calc(xflow1,Tt1,pt1,kappa,r,A1,Ts0,icase) 83 pspt0=(Ts0/Tt1)**(kappa/(kappa-1)) 84 call wpi(w1, pspt0, Q0, 85 & dsqrt(Tt1),kappa,R) 86! 87! Flow velocity at outlet 88 call ts_calc(xflow2,Tt1,pt1,kappa,r,A2,Ts1,icase) 89 pspt1=(Ts1/Tt1)**(kappa/(kappa-1)) 90 call wpi(w2, pspt1, Q1, 91 & dsqrt(Tt2),kappa,R) 92! 93 w2w1=w2/w1 94 w1w2=w1/w2 95! 96! Main branch 97 if(ichan_num.eq.1) then 98! 99! Zeta as in Calculix 100 zeta=0.4d0*(1-W2W1)**2 101! 102 zeta=zeta*(W1W2)**2 103! 104! First branch 105 elseif((ichan_num.eq.2).or.(ichan_num.eq.3)) then 106 hq=dh2/dh1 107 if(alpha.le.60.or.hq.le.2.d0/3.d0) then 108 zeta=0.95d0*((W2W1-2d0*dcos(alpha*pi/180)) 109 & *W2W1+1.d0) 110 zeta=zeta*(W1W2)**2 111 else 112 z2d390=0.95d0*((W2W1-2d0*dcos(90.d0*pi/180)) 113 & *W2W1+1.d0) 114 z1p090=0.95d0*(0.34d0+W2W1**2) 115 z90=z2d390+(3*hq-2.d0)*(z1p090-z2d390) 116 Z60=0.95d0*((W2W1-2d0*dcos(60.d0*pi/180)) 117 & *W2W1+1.d0) 118 zeta=z60+(alpha/30.d0-2.d0)*(z90-z60) 119 zeta=zeta*(W1W2)**2 120 endif 121! 122 endif 123! 124! zeta_fac for side branches are all =1 125! main branch can be set by the user in ACC Designer 126 zeta=zeta*zeta_fac 127! 128 if(icrit2.ne.1) then 129 if(icrit1.ne.1) then 130 f=pt2-pt1*pspt1**zeta 131 else 132 f=xflow2*dsqrt(Tt1)/pt1/A2-Q_crit 133 endif 134 else 135 f=xflow2*dsqrt(Tt1)/pt2/A2-Q_crit 136 endif 137! 138 if(iflag.eq.3) then 139! 140 write(1,57)' zeta= ',zeta 141 57 format(1x,a,f9.4) 142! 143 else if (iflag.eq.4) then 144! 145! Calculate Mach numbers 146 call machpi(M1,pspt0,kappa,R) 147 call ts_calc(xflow2,Tt2,pt2,kappa,r,A2,Ts2,icase) 148! Pressure ratio 149 pspt2=(Ts2/Tt2)**(kappa/(kappa-1)) 150 call machpi(M2,pspt2,kappa,R) 151 152 write(1,80)'Inlet: Tt1= ',Tt1, 153 & ', pt1= ',pt1,', M1= ',M1 154 155 write(1,77)'mass flow = ',xflow2,', kappa = ',kappa, 156 & ', zeta= ',zeta 157 158 write(1,80)'Outlet: Tt2= ',Tt2, 159 & ', pt2= ',pt2,', M2= ',M2 160 161 80 format(3x,a,f10.6,a,f10.2,a,f10.6) 162 77 format(3x,a,f10.6,a,f10.2,a,f10.6) 163 164 endif 165 166! 167 calc_residual_cross_split=f 168! 169 return 170 end 171