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! this subroutine enables to calculate a compressibility correction factor 20! following the results that can be found in: 21! 22! S.L.Bragg 23! "Effect of conpressibility on the discharge coefficient of orifices 24! and convergent nozzles" 25! Journal of Mechanical engineering vol 2 No 1 1960 26! 27! author: Yannick Muller 28! 29 subroutine cd_bragg(cd,p2p1,cdbragg,itype) 30! 31 implicit none 32! 33! itype is used in the proprietary version of cd_bragg 34! 35 integer nx,ny,idx,idy,i,j,itype 36! 37 real*8 cd,p2p1,cdbragg,z1,z2,z3,z4,et,xi 38! 39 real*8 cd_tab (12) 40 data cd_tab 41 & /0.457d0,0.500d0,0.550d0,0.600d0,0.650d0,0.700d0, 42 & 0.750d0,0.800d0,0.850d0,0.900d0,0.950d0,1.000d0/ 43! 44 real*8 p2p1_tab (19) 45 data p2p1_tab 46 & /0.00d0,0.10d0,0.15d0,0.20d0,0.25d0,0.30d0,0.35d0,0.40d0, 47 & 0.45d0,0.50d0,0.55d0,0.60d0,0.65d0,0.70d0,0.75d0,0.80d0, 48 & 0.85d0,0.90d0,1.00d0/ 49! 50 real*8 cd_bragg_tab(19,12) 51 data ((cd_bragg_tab(i,j),i=1,19),j=1,12) 52 & /0.754d0,0.735d0,0.724d0,0.712d0,0.701d0,0.688d0,0.672d0, 53 & 0.655d0,0.633d0,0.610d0,0.590d0,0.570d0,0.549d0,0.530d0, 54 & 0.514d0,0.500d0,0.488d0,0.477d0,0.454d0, 55! 56 & 0.789d0,0.770d0,0.760d0,0.749d0,0.747d0,0.733d0,0.709d0, 57 & 0.691d0,0.672d0,0.650d0,0.628d0,0.606d0,0.588d0,0.572d0, 58 & 0.558d0,0.543d0,0.531d0,0.520d0,0.500d0, 59! 60 & 0.833d0,0.815d0,0.805d0,0.796d0,0.783d0,0.771d0,0.758d0, 61 & 0.740d0,0.720d0,0.700d0,0.675d0,0.655d0,0.638d0,0.621d0, 62 & 0.607d0,0.592d0,0.580d0,0.569d0,0.550d0, 63! 64 & 0.870d0,0.855d0,0.846d0,0.828d0,0.827d0,0.815d0,0.801d0, 65 & 0.786d0,0.769d0,0.749d0,0.725d0,0.704d0,0.685d0,0.670d0, 66 & 0.654d0,0.641d0,0.630d0,0.619d0,0.600d0, 67! 68 & 0.902d0,0.890d0,0.882d0,0.875d0,0.867d0,0.855d0,0.842d0, 69 & 0.830d0,0.811d0,0.792d0,0.773d0,0.751d0,0.732d0,0.718d0, 70 & 0.700d0,0.689d0,0.678d0,0.668d0,0.650d0, 71! 72 & 0.929d0,0.920d0,0.912d0,0.908d0,0.900d0,0.890d0,0.880d0, 73 & 0.869d0,0.852d0,0.835d0,0.815d0,0.794d0,0.778d0,0.761d0, 74 & 0.749d0,0.736d0,0.725d0,0.716d0,0.700d0, 75! 76 & 0.952d0,0.946d0,0.940d0,0.936d0,0.930d0,0.921d0,0.913d0, 77 & 0.903d0,0.889d0,0.873d0,0.854d0,0.836d0,0.820d0,0.808d0, 78 & 0.796d0,0.785d0,0.775d0,0.766d0,0.750d0, 79! 80 & 0.970d0,0.966d0,0.962d0,0.958d0,0.953d0,0.948d0,0.941d0, 81 & 0.935d0,0.923d0,0.909d0,0.890d0,0.874d0,0.860d0,0.849d0, 82 & 0.838d0,0.829d0,0.820d0,0.812d0,0.800d0, 83! 84 & 0.983d0,0.9805d0,0.98d0,0.978d0,0.975d0,0.970d0,0.965d0, 85 & 0.958d0,0.950d0,0.949d0,0.926d0,0.911d0,0.900d0,0.890d0, 86 & 0.881d0,0.874d0,0.867d0,0.860d0,0.850d0, 87! 88 & 0.992d0,0.991d0,0.990d0,0.989d0,0.988d0,0.985d0,0.981d0, 89 & 0.980d0,0.973d0,0.967d0,0.956d0,0.943d0,0.935d0,0.928d0, 90 & 0.920d0,0.915d0,0.910d0,0.907d0,0.900d0, 91! 92 & 0.999d0,0.999d0,0.998d0,0.998d0,0.998d0,0.997d0,0.995d0, 93 & 0.992d0,0.990d0,0.988d0,0.981d0,0.975d0,0.970d0,0.964d0, 94 & 0.960d0,0.958d0,0.954d0,0.952d0,0.950d0, 95! 96 & 1.000d0,1.000d0,1.000d0,1.000d0,1.000d0,1.000d0,1.000d0, 97 & 1.000d0,1.000d0,1.000d0,1.000d0,1.000d0,1.000d0,1.000d0, 98 & 1.000d0,1.000d0,1.000d0,1.000d0,1.000d0/ 99! 100 nx=19 101 ny=12 102! 103 call ident(p2p1_tab,p2p1,nx,idx) 104 call ident(cd_tab,cd,ny,idy) 105! 106 if (idx.eq.0) then 107 if(idy.eq.0) then 108 cdbragg=cd_bragg_tab(1,1) 109 else 110 if(idy.eq.ny) then 111 cdbragg=cd_bragg_tab(1,ny) 112 else 113 cdbragg=cd_bragg_tab(1,idy)+(cd_bragg_tab(1,idy+1) 114 & -cd_bragg_tab(1,idy)) 115 & *(cd-cd_tab(idy))/(cd_tab(idy+1)-cd_tab(idy)) 116 endif 117 endif 118! 119 elseif(idx.ge.nx) then 120 if(idy.le.0) then 121 cdbragg=cd_bragg_tab(nx,1) 122 else 123 if(idy.ge.ny) then 124 cdbragg=cd_bragg_tab(nx,ny) 125 else 126 cdbragg=cd_bragg_tab(nx,idy)+ 127 & (cd_bragg_tab(nx,idy+1)-cd_bragg_tab(nx,idy)) 128 & *(cd-cd_tab(idy))/(cd_tab(idy+1)-cd_tab(idy)) 129 endif 130 endif 131 else 132 if(idy.le.0) then 133! 134 cdbragg=cd_bragg_tab(idx,1)+(cd_bragg_tab(idx+1,1) 135 & -cd_bragg_tab(idx,1)) 136 & *(p2p1-p2p1_tab(idx))/(p2p1_tab(idx+1) 137 & -p2p1_tab(idx)) 138 elseif(idy.ge.ny) then 139 cdbragg=cd_bragg_tab(idx,ny)+(cd_bragg_tab(idx+1,ny) 140 & -cd_bragg_tab(idx,ny)) 141 & *(p2p1-p2p1_tab(idx))/(p2p1_tab(idx+1) 142 & -p2p1_tab(idx)) 143 else 144 xi=(p2p1-p2p1_tab(idx))/(p2p1_tab(idx+1) 145 & -p2p1_tab(idx)) 146 et=(cd-cd_tab(idy))/(cd_tab(idy+1)-cd_tab(idy)) 147 z1=cd_bragg_tab(idx,idy) 148 z2=cd_bragg_tab(idx+1,idy) 149 z3=cd_bragg_tab(idx,idy+1) 150 z4=cd_bragg_tab(idx+1,idy+1) 151 cdbragg=(1-xi)*(1-et)*z1+(1-xi)*et*z3 152 & +xi*(1-et)*z2+xi*et*z4 153 endif 154 endif 155! 156 return 157 end 158