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