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      subroutine wye(node1,node2,nodem,nelem,lakon,kon,ipkon,
19     &     nactdog,identity,ielprop,prop,iflag,v,xflow,f,
20     &     nodef,idirf,df,cp,r,physcon,numf,set,mi,ider,ttime,time,
21     &     iaxial,iplausi,dvi)
22!
23!     A wye split element(zeta calculation according to Idel'chik)
24!     Written by Yavor Dobrev
25!     For an explanation of the parameters see tee.f
26!
27!     author: Yannick Muller
28!
29      implicit none
30!
31      logical identity
32!
33      character*81 set(*)
34      character*8 lakon(*)
35!
36      integer nelem,nactdog(0:3,*),node1,node2,nodem,nodem1,numf,kon(*),
37     &ipkon(*),ielprop(*),nodef(*),idirf(*),index,iflag,inv,mi(2),
38     &nelem1,ichan_num,ider,icase,i,iaxial,iplausi
39!
40      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,R,Tt1,Tt2,pt1,pt2,
41     &cp,physcon(*),km1,kp1,kdkm1,pt2pt1,pt1pt2,pt2pt1_crit,tdkp1,
42     &A,A1,A2,A_s,calc_residual_wye,dh1,dh2,alpha,xflow1,xflow2,
43     &pi,zeta_fac,Ts0,pspt0,pspt2,M1,M2,Ts2,ttime,time,zeta,dvi
44!
45!
46!
47      index=ielprop(nelem)
48!
49      if (iflag.eq.0.d0) then
50         identity=.true.
51         if(nactdog(2,node1).ne.0)then
52            identity=.false.
53         elseif(nactdog(2,node2).ne.0)then
54            identity=.false.
55         elseif(nactdog(1,nodem).ne.0)then
56            identity=.false.
57         endif
58!
59      elseif (iflag.eq.1)then
60         if(v(1,nodem).ne.0.d0) then
61            xflow=v(1,nodem)
62            return
63         endif
64!
65         kappa=(cp/(cp-R))
66         kp1=kappa+1d0
67         km1=kappa-1d0
68         kdkm1=kappa/km1
69         tdkp1=2.d0/kp1
70!
71         if(nelem.eq.nint(prop(index+2))) then
72            A=prop(index+4)
73         elseif(nelem.eq.nint(prop(index+3)))then
74            A=prop(index+6)
75         endif
76!
77         pt1=v(2,node1)
78         pt2=v(2,node2)
79!
80         if(pt1.ge.pt2) then
81            inv=1
82            Tt1=v(0,node1)-physcon(1)
83            Tt2=v(0,node2)-physcon(1)
84         else
85            inv=-1
86            pt1=v(2,node2)
87            pt2=v(2,node1)
88            Tt1=v(0,node2)-physcon(1)
89            Tt2=v(0,node1)-physcon(1)
90         endif
91!
92         pt1pt2=pt1/pt2
93         pt2pt1=1/pt1pt2
94!
95         pt2pt1_crit=tdkp1**kdkm1
96!
97         if(pt2pt1.gt.pt2pt1_crit) then
98            xflow=inv*pt1*A*dsqrt(2.d0*kdkm1*pt2pt1**(2.d0/kappa)
99     &           *(1.d0-pt2pt1**(1.d0/kdkm1))/r)/dsqrt(Tt1)
100         else
101            xflow=inv*pt1*A*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
102     &           dsqrt(Tt1)
103         endif
104         xflow=xflow/50
105!
106      elseif (iflag.eq.2)then
107!
108         numf=6
109!
110         kappa=(cp/(cp-R))
111         pi=4.d0*datan(1.d0)
112!
113!        Inlet conditions are the same for both branches
114!
115!        Determining the previous element as
116!        the incoming mass flow is defined there
117!
118         nelem1=nint(prop(index+1))
119         nodem1=kon(ipkon(nelem1)+2)
120!
121!        Inlet conditions
122!
123         pt1=v(2,node1)
124         Tt1=v(0,node1)-physcon(1)
125         xflow1=v(1,nodem1)*iaxial
126         A1 = prop(index+4)
127!
128!        Outlet conditions
129!
130         Tt2=v(0,node2)
131         xflow2=v(1,nodem)*iaxial
132         pt2=v(2,node2)
133!
134         if(nelem.eq.nint(prop(index+2))) then
135!
136            A2 = A1
137            A_s = prop(index+6)
138!
139            dh1 = prop(index+9)
140            if(dh1.eq.0.d0) then
141               dh1 = dsqrt(4*A1/pi)
142            endif
143!
144            dh2 = dh1
145            alpha = 0
146            ichan_num = 1
147            zeta_fac = prop(index+13)
148!
149         elseif(nelem.eq.nint(prop(index+3))) then
150!
151            ichan_num = 2
152            A2 = prop(index+6)
153!
154            dh1 = prop(index+9)
155            if(dh1.eq.0.d0) then
156               dh1 = dsqrt(4*A1/pi)
157            endif
158!
159            dh2 = prop(index+10)
160            if(dh2.eq.0.d0) then
161               dh2 = dsqrt(4*A2/pi)
162            endif
163!
164            alpha = prop(index+8)
165            zeta_fac = prop(index+14)
166!
167         endif
168!
169!        Set the node numbers for the degrees of freedom
170         nodef(1)=node1
171         nodef(2)=node1
172         nodef(3)=nodem1
173         nodef(4)=nodem
174         nodef(5)=node2
175         nodef(6)=node2
176!
177!        Sets the types of the degrees of freedom
178         idirf(1)=2
179         idirf(2)=0
180         idirf(3)=1
181         idirf(4)=1
182         idirf(5)=2
183         idirf(6)=0
184!
185         if(ider.eq.0.d0) then
186!           Residual
187            f=calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2,
188     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,iflag
189     &,zeta)
190         else
191!           Derivatives
192            call calc_ider_wye(df,pt1,Tt1,xflow1,xflow2,pt2,
193     &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,iflag
194     &,zeta)
195         endif
196!
197      elseif(iflag.eq.3) then
198!
199         kappa=(cp/(cp-R))
200!     setting icase (always adiabatic)
201         icase=0;
202!
203         pi=4.d0*datan(1.d0)
204!
205!        Inlet conditions are the same for both branches
206!
207!        Determining the previous element as
208!        the incoming mass flow is defined there
209         nelem1=nint(prop(index+1))
210         nodem1=kon(ipkon(nelem1)+2)
211!
212!        Inlet conditions
213         pt1=v(2,node1)
214         Tt1=v(0,node1)-physcon(1)
215         xflow1=v(1,nodem1)*iaxial
216         A1 = prop(index+4)
217!
218!        Outlet conditions
219         Tt2=v(0,node2)
220         xflow2=v(1,nodem)*iaxial
221         pt2=v(2,node2)
222!
223         if(nelem.eq.nint(prop(index+2))) then
224!
225            A2 = A1
226            A_s = prop(index+6)
227!
228            dh1 = prop(index+9)
229            if(dh1.eq.0.d0) then
230               dh1 = dsqrt(4*A1/pi)
231            endif
232!
233            dh2 = dh1
234            alpha = 0
235            ichan_num = 1
236            zeta_fac = prop(index+13)
237!
238         elseif(nelem.eq.nint(prop(index+3))) then
239!
240            ichan_num = 2
241            A2 = prop(index+6)
242!
243            dh1 = prop(index+9)
244            if(dh1.eq.0.d0) then
245               dh1 = dsqrt(4*A1/pi)
246            endif
247!
248            dh2 = prop(index+10)
249            if(dh2.eq.0.d0) then
250               dh2 = dsqrt(4*A2/pi)
251            endif
252!
253            alpha = prop(index+8)
254            zeta_fac = prop(index+14)
255!
256         endif
257!
258!        Write the main information about the element
259!
260!        Flow velocity at inlet
261         call ts_calc(xflow1,Tt1,pt1,kappa,r,A1,Ts0,icase)
262         pspt0 = (Ts0/Tt1)**(kappa/(kappa-1))
263!        Calculate Mach numbers
264         call machpi(M1,pspt0,kappa,R)
265         call ts_calc(xflow2,Tt2,pt2,kappa,r,A2,Ts2,icase)
266!        Pressure ratio
267         pspt2 = (Ts2/Tt2)**(kappa/(kappa-1))
268         call machpi(M2,pspt2,kappa,R)
269!
270         write(1,*) ''
271         write(1,55) ' from node ',node1,
272     &        ' to node ', node2,':   air massflow rate= ',xflow
273!
274         write(1,56)'       Inlet node ',node1,':    Tt1= ',Tt1,
275     &        ' , Ts1= ',Ts0,' , Pt1= ',pt1,
276     &        ', M1= ',M1
277         write(1,*)'             Element ',nelem,lakon(nelem)
278     &        ,', Branch ',ichan_num
279!
280 55      format(1x,a,i6,a,i6,a,e11.4,a)
281 56      format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
282!
283!     Set ider to calculate the residual
284         ider = 0
285!
286!     Calculate the element one last time with enabled output
287         f=calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2,
288     &           Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,
289     &           kappa,R,ider,iflag,zeta)
290!
291         write(1,56)'      Outlet node ',node2,':   Tt2= ',Tt2,
292     &        ' , Ts2= ',Ts2,' , Pt2= ',pt2,
293     &        ', M2= ',M2
294!
295      endif
296!
297      xflow=xflow/iaxial
298      df(3)=df(3)*iaxial
299      df(4)=df(4)*iaxial
300!
301      return
302      end
303