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