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