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