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 subroutine cd_ms_ms(p1,p2,T1,rad,d,xl,kappa,r,reynolds,u,vid,cd) 20! 21! This subroutine enables to calculate the discharge coefficient for an 22! orifice (shap edged , rotating..) following the results obtained 23! by Mcgreehan and Schotsch 24! The decription of the method can be found in : 25! "Flow characteristics of long orifices with rotation and 26! corner radiusing" 27! ASME 87-GT-162 28! 29! author: Yannick Muller 30! 31 implicit none 32! 33 real*8 p1,p2,T1,rad,d,xl,kappa,r,reynolds,u,cd,qlim,q, 34 & c1,c2,c3,fakt,aux,rzd,lkorr,qkorr,rv,vid 35! 36 qlim=10.d0 37! 38! taking in account the influence of the Reynolds number 39! 40 cd=0.5885d0+372.d0/reynolds 41 cd=min(cd,1.d0) 42! 43! taking in account the edge radius 44! 45 rzd=rad/d 46 aux=exp(-(3.5d0*rzd+5.5d0)*rzd) 47 fakt=aux+0.008d0*(1.d0-aux) 48 cd=1.d0-fakt*(1.d0-cd) 49 cd=min(max(cd,0.d0),1.d0) 50! 51! taking in account the lenght of the orifice 52! 53 lkorr=xl-rad 54 q=lkorr/d 55 qkorr=min(q,qlim) 56 fakt=(1.d0+1.3d0*exp(-1.606d0*qkorr**2.d0))* 57 & (0.435d0+0.021d0*qkorr)/(2.3d0*0.435d0) 58 cd=1.d0-fakt*(1.d0-cd) 59 cd=min(max(cd,0.d0),1.d0) 60! 61! taking in account the tangential velocity 62! 63 if(u.ne.0d0) then 64 vid=dsqrt(2.d0*kappa/(kappa-1.d0)*r*T1* 65 & (1.d0-(p2/p1)**((kappa-1.d0)/kappa))) 66 rv=u/vid*(0.6d0/cd)**3 67 c1=exp(-rv**1.2d0) 68 c2=0.5d0*rv**0.6d0*dsqrt(0.6d0/cd) 69 c3=exp(-0.5d0*rv**0.9d0) 70 cd=cd*(c1+c2*c3) 71 cd=min(max(cd,0.d0),1.d0) 72! 73 endif 74! 75! 76 return 77 end 78