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