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!
20!     this subroutine enables to calculate thecorrection factor of the discharge
21!     coefficient of a labyrinth with one fin as a function of the ratio b/s and the
22!     pressure ratio Pdownstream/Pupstream
23!     the results are interpolated
24!
25!     the relevant data can be found in:
26!     "Air System Correlations Part 1: Labyrinth Seals"
27!     H.Zimmermann and K.H. Wollf
28!     ASME 98-GT-206
29!     fig 12 p 7
30!
31!     author: Yannick Muller
32!
33      subroutine cd_lab_correction(p1p2,s,b,cd_correction)
34!
35      implicit none
36!
37      integer nx,ny,idx,idy
38!
39      real*8 s,b,cd_correction,z1,z2,z3,z4,xi,et,szb,p1p2
40!
41      real*8 puszpds_tab(7)
42      data puszpds_tab
43     &    /1.d0,1.2d0,1.4d0,1.6d0,1.8d0,2.d0,2.5d0/
44!
45      real*8 szb_tab(9)
46      data szb_tab
47     &    /0.25d0,0.5d0,1.d0,1.5d0,2d0,2.5d0,3d0,3.5d0,4d0/
48!
49      real*8 cd_correction_tab(9,7)
50      data cd_correction_tab
51     &  /1.05d0,1.07d0,1.03d0,0.98d0,0.95d0,0.94d0,0.95d0,0.95d0,0.95d0,
52     &   1.15d0,1.07d0,1.02d0,0.95d0,0.92d0,0.91d0,0.91d0,0.92d0,0.92d0,
53     &   1.15d0,1.05d0,0.98d0,0.91d0,0.88d0,0.86d0,0.86d0,0.87d0,0.87d0,
54     &   1.15d0,1.04d0,0.95d0,0.87d0,0.85d0,0.84d0,0.83d0,0.83d0,0.83d0,
55     &   1.15d0,1.03d0,0.91d0,0.85d0,0.81d0,0.80d0,0.80d0,0.80d0,0.80d0,
56     &   1.15d0,1.01d0,0.90d0,0.82d0,0.79d0,0.79d0,0.77d0,0.77d0,0.77d0,
57     &   1.10d0,1.00d0,0.88d0,0.79d0,0.75d0,0.74d0,0.73d0,0.72d0,0.70d0/
58!
59      szb=s/b
60!
61      nx=9
62      ny=7
63!
64!      p1p2=1/p2p1
65!      if ((p1p2.ge.2.5d0).or.(szb.ge.4d0))then
66!         write(*,*) '*WARNING in cd_lab_correction'
67!         write(*,*) 'p1p2>2.5 or szb>4'
68!         write(*,*) 'check input file'
69!         write(*,*) 'calculation will proceed using cd_lab_correction=1'
70!         cd_correction=1.d0
71!        return
72!      endif
73!
74      call ident(puszpds_tab,p1p2,ny,idy)
75      call ident(szb_tab,szb,nx,idx)
76!
77      if (idx.eq.0) then
78         if(idy.eq.0) then
79            cd_correction=cd_correction_tab(1,1)
80         else
81            if(idy.eq.ny) then
82               cd_correction=cd_correction_tab(1,ny)
83            else
84               cd_correction=cd_correction_tab(1,idy)
85     &           +(cd_correction_tab(1,idy+1)-cd_correction_tab(1,idy))
86     &              *(szb-szb_tab(idx))/(szb_tab(idx+1)-szb_tab(idx))
87            endif
88         endif
89!
90      elseif(idx.ge.nx) then
91         if(idy.le.0) then
92            cd_correction=cd_correction_tab(nx,1)
93         else
94            if(idy.ge.ny) then
95               cd_correction=cd_correction_tab(nx,ny)
96            else
97            cd_correction=cd_correction_tab(nx,idy)
98     &     +(cd_correction_tab(nx,idy+1)-cd_correction_tab(nx,idy))
99     &     *(szb-szb_tab(idx))/(szb_tab(idx+1)-szb_tab(idx))
100            endif
101         endif
102      else
103         if(idy.le.0) then
104!
105            cd_correction=cd_correction_tab(idx,1)
106     &          +(cd_correction_tab(idx+1,1)-cd_correction_tab(idx,1))
107     &           *(p1p2-puszpds_tab(idy))/(puszpds_tab(idy+1)
108     &           -puszpds_tab(idy))
109         elseif(idy.ge.ny) then
110            cd_correction=cd_correction_tab(idx,ny)
111     &         +(cd_correction_tab(idx+1,ny)-cd_correction_tab(idx,ny))
112     &           *(p1p2-puszpds_tab(idy))/(puszpds_tab(idy+1)
113     &           -puszpds_tab(idy))
114         else
115            et=(p1p2-puszpds_tab(idy))/(puszpds_tab(idy+1)
116     &           -puszpds_tab(idy))
117            xi=(szb-szb_tab(idx))/(szb_tab(idx+1)-szb_tab(idx))
118            z1=cd_correction_tab(idx,idy)
119            z2=cd_correction_tab(idx+1,idy)
120            z3=cd_correction_tab(idx,idy+1)
121            z4=cd_correction_tab(idx+1,idy+1)
122            cd_correction=(1-xi)*(1-et)*z1+(1-xi)*et*z3
123     &           +xi*(1-et)*z2+xi*et*z4
124         endif
125      endif
126!
127!      if (cd_correction.ge.1.d0)then
128!         cd_correction=1.d0
129!      endif
130!
131      return
132      end
133