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!     calculate the maximal admissible pressure ratio pt2/pt1
20!
21!     1) for icase=1 (sonic conditions are possible at 1):
22!                    assuming M1=1 M2 is calculated iteratively
23!                    using a dichotomy scheme
24!        for icase=2 (sonic conditions are possible at 2):
25!                    assuming M2=1 M1 is calculated iteratively
26!                    using a dichotomy scheme
27!
28!     2)the ratio of the critical pressure ratio
29!       Qred_1/Qred_2crit=Pt2/Pt1=D(M1)/D(M2_crit)
30!       is computed
31!       [D(M)=M*(1+0.5*(kappa-1)*M)**(-0.5*(kappa+1)/(kappa-1))]
32!       (general gas equation)
33!
34!     author: Yannick Muller/Guido Dhondt
35!
36!     Attention: this routine is only used for rotating gas pipes
37!                with rotational speed 0, i.e. static pipes with a possibly
38!                varying cross section. In such pipes the total pressure
39!                always decreases versus length in flow direction
40!                beta > 0: flow from node1 to node2
41!                beta < 0: flow from node2 to node1
42!
43      subroutine pt2zpt1_rot(pt2,pt1,kappa,r,l,pt2zpt1_c,crit,icase,
44     &     M1,M2,ca,cb,cc,alpha,beta,Qred1_crit,Qred2_crit,A1,A2)
45!
46      implicit none
47!
48      logical crit
49!
50      integer icase,i
51!
52      real*8 pt2,pt1,kappa,l,M1,pt2zpt1,pt2zpt1_c,km1,kp1,r,Qred1_crit,
53     &     fmin,Qred2_crit,f,fmax,M1_min,M1_max,Z1,Z1_min,Z1_max,Z2,
54     &     Z2_min,Z2_max,M2,M2_min,M2_max,tcdkm1,km1d2,ca,cb,cc,alpha,
55     &     beta,A1,A2,M1_root,M2_root
56!
57!
58!
59      crit=.false.
60!
61!     useful variables and constants
62!
63      km1=kappa-1.d0
64      kp1=kappa+1.d0
65      km1d2=km1/2.d0
66      tcdkm1=2.d0*cc/(kappa-1.d0)
67!
68      if(icase.eq.1) then
69!
70!        M1=1
71!        computing M2 using dichotomy method (dividing the interval
72!        with the funciton root iteratively by 2)
73!
74         i=1
75!
76!        because of
77!        (alpha+beta*Z2_min)/(alpha+beta))**(cb/beta)
78!        we have to avoid a root in the initial dichotomy interval
79!
80!        for icase.eq.1 we have:
81!     alpha+beta*Z is somewhere in the interval [0,1] negative;
82!
83!        if beta >= 0, alpha<0 must be true;
84!        since alpha+beta*Z is monotonically increasing, a root in the
85!        interval is only possible for alpha+beta>0
86!
87!        if beta <= 0, alpha+beta<0 must be true;
88!        since alpha+beta*Z is monotonically decreasing, a root in the
89!        interval is only possible for alpha>0
90!
91         if(beta.ge.0.d0) then
92            if(alpha+beta.gt.0.d0) then
93               M2_root=dsqrt(-alpha/beta)
94               M2_min=max(0.001d0,M2_root+0.001d0)
95               M2_max=0.999d0
96            else
97               M2_min=0.001d0
98               M2_max=0.999d0
99            endif
100         else
101            if(alpha.gt.0.d0) then
102               M2_root=dsqrt(-alpha/beta)
103               M2_min=max(0.001d0,M2_root+0.001d0)
104               M2_max=0.999d0
105            else
106               M2_min=0.001d0
107               M2_max=0.999d0
108            endif
109         endif
110
111c         M2_min=0.001d0
112c         M2_max=1
113         Z2_min=M2_min**2
114         Z2_max=M2_max**2
115!
116         fmin=dlog(Z2_min**ca*
117     &        ((alpha+beta*Z2_min)/(alpha+beta))**(cb/beta)*
118     &        ((1.d0+km1d2*Z2_min)/(1.d0+km1d2))**(tcdkm1))-l
119!
120         fmax=dlog(Z2_max**ca*
121     &        ((alpha+beta*Z2_max)/(alpha+beta))**(cb/beta)*
122     &        ((1.d0+km1d2*Z2_max)/(1.d0+km1d2))**(tcdkm1))-l
123!
124         if(fmin*fmax.gt.0.d0) then
125            pt2zpt1_c=1.d30
126            Qred2_crit=1.d30
127            crit=.false.
128            return
129         endif
130!
131         do
132            i=i+1
133            M2=(M2_min+M2_max)*0.5d0
134            Z2=M2**2
135!
136            f=dlog(Z2**ca*
137     &           ((alpha+beta*Z2)/(alpha+beta))**(cb/beta)*
138     &           ((1.d0+km1d2*Z2)/(1.d0+km1d2))**(tcdkm1))-l
139!
140            if(abs(f).le.1d-6) then
141               exit
142            endif
143            if(i.gt.50) then
144               exit
145            endif
146!
147            if(fmin*f.le.0.d0) then
148               M2_max=M2
149               fmax=f
150            else
151               M2_min=M2
152               fmin=f
153            endif
154         enddo
155         write(*,*) 'pt2zpt1_rot M2 ',M2
156!
157         pt2zpt1_c=(0.5d0*kp1)**(-0.5d0*kp1/km1)
158     &        *(1.d0+0.5d0*km1*Z2)**(0.5d0*kp1/km1)*A1/(A2*M2)
159!
160         Qred2_crit=M2*dsqrt(kappa/r)
161     &        *(1.d0+0.5d0*km1*Z2)**(-0.5d0*kp1/km1)
162!
163         pt2zpt1=pt2/pt1
164         write(*,*) 'pt2zpt1_rot pt2/pt1 ',pt2zpt1_c,pt2zpt1
165         if(beta.ge.0.d0) then
166            if(pt2zpt1.le.pt2zpt1_c) then
167               crit=.true.
168            endif
169         else
170            if(pt2zpt1.ge.pt2zpt1_c) then
171               crit=.true.
172            endif
173         endif
174!
175      elseif (icase.eq.2) then
176!
177!        M2=1
178!        computing M1 using dichotomy method (dividing the interval
179!        with the function root iteratively by 2)
180!
181         i=1
182!
183!        because of
184!        (alpha+beta*Z2_min)/(alpha+beta))**(cb/beta)
185!        we have to avoid a root in the initial dichotomy interval
186!
187!        for icase.eq.2 we have:
188!     alpha+beta*Z is somewhere in the interval [0,1] positive;
189
190!        if beta >= 0 alpha+beta*Z is monotonically increasing, and a root
191!        in the interval is only possible for alpha<0
192
193!        if beta <= 0 alpha+beta*Z is monotonically decreasing, and a root
194!        in the interval is only possible for alpha+beta<0
195!
196!
197         if(beta.ge.0.d0) then
198            if(alpha.lt.0.d0) then
199               M1_root=dsqrt(-alpha/beta)
200               M1_min=max(0.001d0,M1_root+0.001d0)
201               M1_max=0.999d0
202            else
203               M1_min=0.001d0
204               M1_max=0.999d0
205            endif
206         else
207            if(alpha+beta.lt.0.d0) then
208               M1_root=dsqrt(-alpha/beta)
209               M1_min=max(0.001d0,M1_root+0.001d0)
210               M1_max=0.999d0
211            else
212               M1_min=0.001d0
213               M1_max=0.999d0
214            endif
215         endif
216c         M1_min=0.001d0
217c         M1_max=1
218         Z1_min=M1_min**2
219         Z1_max=M1_max**2
220!
221         fmin=dlog((1.d0/Z1_min)**ca*
222     &        ((alpha+beta)/(alpha+beta*Z1_min))**(cb/beta)*
223     &        ((1.d0+km1d2)/(1.d0+km1d2*Z1_min))**(tcdkm1))-l
224!
225         fmax=dlog((1.d0/Z1_max)**ca*
226     &        ((alpha+beta)/(alpha+beta*Z1_max))**(cb/beta)*
227     &        ((1.d0+km1d2)/(1.d0+km1d2*Z1_max))**(tcdkm1))-l
228!
229         if(fmin*fmax.gt.0.d0) then
230            pt2zpt1_c=1.d30
231            Qred1_crit=1.d30
232            crit=.false.
233            return
234         endif
235!
236         do
237            i=i+1
238            M1=(M1_min+M1_max)*0.5d0
239            Z1=M1**2
240!
241            f=dlog((1.d0/Z1)**ca*
242     &           ((alpha+beta)/(alpha+beta*Z1))**(cb/beta)*
243     &           ((1.d0+km1d2)/(1.d0+km1d2*Z1))**(tcdkm1))-l
244!
245            if(abs(f).le.1d-6) then
246               exit
247            endif
248            if(i.gt.50) then
249               exit
250            endif
251!
252            if(fmin*f.le.0.d0) then
253               M1_max=M1
254               fmax=f
255            else
256               M1_min=M1
257               fmin=f
258            endif
259         enddo
260!
261         pt2zpt1_c=M1*(0.5d0*kp1)**(0.5d0*kp1/km1)
262     &        *(1.d0+0.5d0*km1*Z1)**(-0.5d0*kp1/km1)*A1/A2
263!
264         Qred1_crit=M1*dsqrt(kappa/r)
265     &        *(1.d0+0.5d0*km1*Z1)**(-0.5d0*kp1/km1)
266!
267         pt2zpt1=pt2/pt1
268         if(beta.ge.0.d0) then
269            if(pt2zpt1.le.pt2zpt1_c) then
270               crit=.true.
271            endif
272         else
273            if(pt2zpt1.ge.pt2zpt1_c) then
274               crit=.true.
275            endif
276         endif
277!
278      endif
279!
280      return
281      end
282
283
284