CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
zeta_calc.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine zeta_calc (nelem, prop, ielprop, lakon, reynolds, zeta, isothermal, kon, ipkon, R, kappa, v, mi, iaxial)
 

Function/Subroutine Documentation

◆ zeta_calc()

subroutine zeta_calc ( integer  nelem,
real*8, dimension(*)  prop,
integer, dimension(*)  ielprop,
character*8, dimension(*)  lakon,
real*8  reynolds,
real*8  zeta,
logical  isothermal,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
real*8  R,
real*8  kappa,
real*8, dimension(0:mi(2),*)  v,
integer, dimension(*)  mi,
integer  iaxial 
)
35 !
36  implicit none
37 !
38  logical isothermal
39 !
40  character*8 lakon(*)
41 !
42  integer ielprop(*),nelem,iexp(2),i,j,ier,iwrit1,iexp3(2),
43  & iwrit2,nelem_ref,ipkon(*),kon(*),nelem0,nelem1,nelem2,node10,
44  & node20,nodem0,node11,node21,nodem1,node12,node22,nodem2,
45  & iexpbr1(2) /11,11/,icase,node0,node1,node2,mi(*),n0,n1,n2,
46  & n5,n6,n8,n10,n11,n12,n14,n15,n22,iaxial,index
47 !
48  real*8 zeta,prop(*),lzd,reynolds,ereo,fa2za1,zetap,zeta0,
49  & lambda,thau,a1,a2,dh,dl,a2za1,ldumm,dhdumm,ks,
50  & form_fact,zeta01,zeta02,alpha,rad,delta,a0,b0,azb,rzdh,
51  & a,c,rei,lam,ai,b1,c1,b2,c2,zeta1,re_val,k,ldre,
52  & zetah,cd,cdu,km,tt0,ts0,tt1,ts1,tt2,ts2,
53  & rho0,rho1,rho2,v0,v1,v2,a0a1,a0a2,zetlin,lam10,lam20,pi,
54  & alpha1,alpha2,r,kappa,ang1s,ang2s,cang1s,cang2s,
55  & v(0:mi(2),*),v1v0,v2v0,z1_60,z1_90,
56  & z2_60,z2_90,afakt,v2v0l,kb,ks2,a2a0,z90lim11,z90lim51,
57  & lam11,lam12,lam21,lam22,w2w0,w1w0,dh0,dh2,hq,z2d390,
58  & z1p090,z90,z60,pt0,pt2,pt1,m0,m1,m2,w0w1,w0w2,
59  & xflow0,xflow1,xflow2,qred_0, qred_1, qred_2,qred_crit
60 !
61 ! THICK EDGED ORIFICE IN STRAIGHT CONDUIT (L/DH > 0.015)
62 ! I.E. IDEL' CHIK (SECTION III PAGE 140)
63 !
64 ! I.E. IDEL'CHIK 'HANDBOOK OF HYDRAULIC RESISTANCE'
65 ! 2nd edition 1986,HEMISPHERE PUBLISHING CORP.
66 ! ISBN 0-899116-284-4
67 !
68 ! ***** long orifice *****
69 !
70 ! DIAGRAMS 4-19 p 175 - Reynolds R:epsilon^-_oRe
71 !
72  real*8 xre (14), yere (14)
73  data xre /25.,40.,60.0,100.,200.,400.,1000.,2000.,4000.,
74  & 10000.,20000.,100000.,200000.,1000000./
75  data yere /0.34,0.36,0.37,0.40,0.42,0.46,0.53,0.59,
76  & 0.64,0.74,0.81,0.94,0.95,0.98/
77 !
78 ! Diagram 4-19 p 175 - Reynolds | A1/A2 R: zeta_phi
79 !
80  real*8 zzeta (15,11)
81  data ((zzeta(i,j),i=1,15),j=1,11)
82  & /15.011 ,25.0,40.0,60.0,100.0,200.0,400.0,1000.0,2000.0,
83  & 4000.0,10000.0,20000.0,100000.0,200000.0,1000000.0,
84  & 0.00 ,1.94,1.38,1.14,0.89,0.69,0.64,0.39,0.30,0.22,0.15,
85  & 0.11,0.04,0.01,0.00,
86  & 0.20 ,1.78,1.36,1.05,0.85,0.67,0.57,0.36,0.26,0.20,0.13,
87  & 0.09,0.03,0.01,0.00,
88  & 0.30 ,1.57,1.16,0.88,0.75,0.57,0.43,0.30,0.22,0.17,0.10,
89  & 0.07,0.02,0.01,0.00,
90  & 0.40 ,1.35,0.99,0.79,0.57,0.40,0.28,0.19,0.14,0.10,0.06,
91  & 0.04,0.02,0.01,0.00,
92  & 0.50 ,1.10,0.75,0.55,0.34,0.19,0.12,0.07,0.05,0.03,0.02,
93  & 0.01,0.01,0.01,0.00,
94  & 0.60 ,0.85,0.56,0.30,0.19,0.10,0.06,0.03,0.02,0.01,0.01,
95  & 0.00,0.00,0.00,0.00,
96  & 0.70 ,0.58,0.37,0.23,0.11,0.06,0.03,0.02,0.01,0.00,0.00,
97  & 0.00,0.00,0.00,0.00,
98  & 0.80 ,0.40,0.24,0.13,0.06,0.03,0.02,0.01,0.00,0.00,0.00,
99  & 0.00,0.00,0.00,0.00,
100  & 0.90 ,0.20,0.13,0.08,0.03,0.01,0.00,0.00,0.00,0.00,0.00,
101  & 0.00,0.00,0.00,0.00,
102  & 0.95 ,0.03,0.03,0.02,0.00,0.00,0.00,0.00,0.00,0.00,0.00,
103  & 0.00,0.00,0.00,0.00/
104 !
105 ! Diagram 4-12 p 169 - l/Dh R: tau
106 !
107  real*8 xlzd (10), ytor (10)
108  data xlzd /0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.6,2.0,2.4/
109  data ytor /1.35,1.22,1.10,0.84,0.42,0.24,0.16,0.07,0.02,0.0/
110  data iexp /10, 1/
111 !
112 ! ***** wall orifice *****
113 !
114 ! THICK-WALLED ORIFICE IN LARGE WALL (L/DH > 0.015)
115 ! I.E. IDL'CHIK (page 174)
116 !
117 ! DIAGRAM 4-18 A - l/Dh R: zeta_o
118 !
119  real*8 xlqd(12)
120  DATA xlqd /
121  & 0.,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,10.0/
122  real*8 yzeta1(12)
123  DATA yzeta1 /
124  & 2.85,2.72,2.6,2.34,1.95,1.76,1.67,1.62,1.6,1.58,1.55,1.55/
125 !
126 ! DIAGRAM 4-19 p175 first line - Re (A1/A2=0) R: zeta_phi
127 !
128  real*8 xre2(14)
129  DATA xre2 /
130  & 25.,40.,60.,100.,200.,400.,1000.,2000.,4000.,10000.,
131  & 20000.,50000.,100000.,1000000./
132  real*8 yzeta2(14)
133  DATA yzeta2 /
134  & 1.94,1.38,1.14,.89,.69,.54,.39,.3,.22,.15,.11,.04,.01,0./
135 !
136 ! Diagram 4-18 p174 first case * (=multiplication) epsilon^-_oRe p 175
137 !
138  real*8 yere2(14)
139  DATA yere2 /
140  & 1.,1.05,1.09,1.15,1.23,1.37,1.56,1.71,1.88,2.17,2.38,2.56,
141  & 2.72,2.85/
142 !
143 ! ***** expansion *****
144 !
145 ! SUDDEN EXPANSION OF A STREAM WITH UNIFORM VELOCITY DISTRIBUTION
146 ! I.E. IDL'CHIK (page 160)
147 !
148 ! DIAGRAM 4-1 - Re | A1/A2 R:zeta
149 !
150  real*8 zzeta3(14,8)
151  DATA zzeta3 /
152  & 14.008, 10.000,15.0,20.0,30.0,40.0,50.0,100.0,200.0,500.0,
153  & 1000.0,2000.0,3000.0,3500.0,
154  & .01 ,3.10,3.20,3.00,2.40,2.15,1.95,1.70,1.65,1.70,2.00,
155  & 1.60,1.00,1.00,
156  & 0.1 ,3.10,3.20,3.00,2.40,2.15,1.95,1.70,1.65,1.70,2.00,
157  & 1.60,1.00,0.81,
158  & 0.2 ,3.10,3.20,2.80,2.20,1.85,1.65,1.40,1.30,1.30,1.60,
159  & 1.25,0.70,0.64,
160  & 0.3 ,3.10,3.10,2.60,2.00,1.60,1.40,1.20,1.10,1.10,1.30,
161  & 0.95,0.60,0.50,
162  & 0.4 ,3.10,3.00,2.40,1.80,1.50,1.30,1.10,1.00,0.85,1.05,
163  & 0.80,0.40,0.36,
164  & 0.5 ,3.10,2.80,2.30,1.65,1.35,1.15,0.90,0.75,0.65,0.90,
165  & 0.65,0.30,0.25,
166  & 0.6 ,3.10,2.70,2.15,1.55,1.25,1.05,0.80,0.60,0.40,0.60,
167  & 0.50,0.20,0.16/
168 !
169  DATA iexp3 /0,0/
170 !
171 ! ***** contraction *****
172 !
173 ! SUDDEN CONTRACTION WITH & WITHOUT CONICAL BELLMOUTH ENTRY
174 ! I.E. IDL'CHIK p 168
175 !
176 ! DIAGRAM 4-10 - Re | A1/A2 R: zeta
177 !
178  real*8 zzeta41(14,7)
179  DATA zzeta41 /
180  & 14.007 ,10.0,20.0,30.0,40.0,50.0,100.0,200.0,500.0,1000.0,
181  & 2000.0,4000.0,5000.0,10000.0,
182  &0.1 ,5.00,3.20,2.40,2.00,1.80,1.30,1.04,0.82,0.64,0.50,
183  & 0.80,0.75,0.50,
184  &0.2 ,5.00,3.10,2.30,1.84,1.62,1.20,0.95,0.70,0.50,0.40,
185  & 0.60,0.60,0.40,
186  &0.3 ,5.00,2.95,2.15,1.70,1.50,1.10,0.85,0.60,0.44,0.30,
187  & 0.55,0.55,0.35,
188  &0.4 ,5.00,2.80,2.00,1.60,1.40,1.00,0.78,0.50,0.35,0.25,
189  & 0.45,0.50,0.30,
190  &0.5 ,5.00,2.70,1.80,1.46,1.30,0.90,0.65,0.42,0.30,0.20,
191  & 0.40,0.42,0.25,
192  &0.6 ,5.00,2.60,1.70,1.35,1.20,0.80,0.56,0.35,0.24,0.15,
193  & 0.35,0.35,0.20/
194 !
195 ! Diagram 3-7 p128 - alpha | l/Dh R: zeta
196 !
197  real*8 zzeta42(10,7)
198  DATA zzeta42 /
199  & 10.007 ,0.,10.0,20.0,30.0,40.0,60.0,100.0,140.0,180.0,
200  & 0.025 ,0.50,0.47,0.45,0.43,0.41,0.40,0.42,0.45,0.50,
201  & 0.050 ,0.50,0.45,0.41,0.36,0.33,0.30,0.35,0.42,0.50,
202  & 0.075 ,0.50,0.42,0.35,0.30,0.26,0.23,0.30,0.40,0.50,
203  & 0.100 ,0.50,0.39,0.32,0.25,0.22,0.18,0.27,0.38,0.50,
204  & 0.150 ,0.50,0.37,0.27,0.20,0.16,0.15,0.25,0.37,0.50,
205  & 0.600 ,0.50,0.27,0.18,0.13,0.11,0.12,0.23,0.36,0.50/
206 !
207 ! ***** bends *****
208 !
209 ! SHARP ELBOW (R/DH=0) AT 0 < DELTA < 180
210 ! I.E. IDL'CHIK page 294
211 ! DIAGRAM 6-5 - a0/b0 R: C1
212 !
213  real*8 xaqb(12)
214  DATA xaqb /
215  & 0 .25,0.50,0.75,1.00,1.50,2.00,3.00,4.00,5.00,6.00,7.00,8.00/
216 !
217  real*8 yc(12)
218  DATA yc /
219  & 1.10,1.07,1.04,1.00,0.95,0.90,0.83,0.78,0.75,0.72,0.71,0.70/
220 !
221 ! DIAGRAM 6-5 - delta R: A
222 !
223  real*8 xdelta(10)
224  DATA xdelta /
225  & 20.0,30.0,45.0,60.0,75.0,90.0,110.,130.,150.,180./
226 !
227  real*8 ya(10)
228  DATA ya /
229  & 2.50,2.22,1.87,1.50,1.28,1.20,1.20,1.20,1.20,1.20/
230 !
231 ! SHARP BENDS 0.5 < R/DH < 1.5 AND 0 < DELTA < 180
232 ! I.E. IDL'CHIK page 289-290
233 ! DIAGRAM 6-1 (- delta from diagram 6-5) R: A1
234 !
235  real*8 ya1(10)
236  DATA ya1 /
237  & 0.31,0.45,0.60,0.78,0.90,1.00,1.13,1.20,1.28,1.40/
238 !
239 ! DIAGRAM 6-1 - R0/D0 R: B1
240 !
241  real*8 xrqdh(8)
242  DATA xrqdh /
243  & 0.50,0.60,0.70,0.80,0.90,1.00,1.25,1.50/
244 !
245  real*8 yb1(8)
246  DATA yb1 /
247  & 1.18,0.77,0.51,0.37,0.28,0.21,0.19,0.17/
248 !
249 ! DIAGRAM 6-1 (- a0/b0 from diagram 6-5) R: C1
250 !
251  real*8 yc1(12)
252  DATA yc1 /
253  & 1.30,1.17,1.09,1.00,0.90,0.85,0.85,0.90,095,0.98,1.00,1.00/
254 !
255 ! SMOOTH BENDS (R/DH > 1.5) AT 0 < DELTA < 180
256 ! I.E. IDL'CHIK
257 !
258 ! DIAGRAM 6-1 - R0/D0 R: B1 (continuation of XRQDH)
259 !
260  real*8 xrzdh(14)
261  DATA xrzdh /
262  & 1.00,2.00,4.00,6.00,8.00,10.0,15.0,20.0,25.0,30.0,35.0,40.0,
263  & 45.0,50.0/
264 !
265  real*8 yb2(14)
266  DATA yb2 /
267  & 0.21,0.15,0.11,0.09,0.07,0.07,0.06,0.05,0.05,0.04,0.04,0.03,
268  & 0.03,0.03/
269 !
270 ! (- a0/b0 from Diagram 6-5) R: C2
271 !
272  real*8 yc2(12)
273  DATA yc2 /
274  & 1.80,1.45,1.20,1.00,0.68,0.45,0.40,0.43,0.48,0.55,0.58,0.60/
275 !
276 ! D.S. MILLER 'INTERNAL FLOW SYSTEMS'
277 ! 1978,vol.5 B.H.R.A FLUID ENGINEERING SERIES
278 ! ISBN 0-900983-78-7
279 !
280 ! SMOOTH BENDS B.H.R.A HANDBOOK P.141
281 !
282  REAL*8 zzetao(14,15)
283  DATA((zzetao(i,j),i=1,14),j=1,8) /
284  & 14.015,0.5,0.6,0.8,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.,
285  & 10.00, 0.030,0.025,0.021,0.016,0.022,0.030,0.034,0.036,0.040,
286  & 0.042,0.043,0.044,0.044,
287  & 15.00, 0.036,0.035,0.025,0.025,0.033,0.042,0.045,0.050,0.055,
288  & 0.055,0.058,0.060,0.063,
289  & 20.00, 0.056,0.046,0.034,0.034,0.045,0.054,0.056,0.062,0.066,
290  & 0.067,0.072,0.075,0.080,
291  & 30.00, 0.122,0.094,0.063,0.056,0.063,0.071,0.075,0.082,0.087,
292  & 0.089,0.097,0.101,0.110,
293  & 40.00, 0.220,0.160,0.100,0.085,0.080,0.086,0.092,0.100,0.106,
294  & 0.122,0.121,0.126,0.136,
295  & 50.00, 0.340,0.245,0.148,0.117,0.097,0.100,0.108,0.116,0.123,
296  & 0.133,0.144,0.150,0.159,
297  & 60.00, 0.480,0.350,0.196,0.150,0.115,0.116,0.122,0.131,0.140,
298  & 0.153,0.164,0.171,0.181/
299  DATA((zzetao(i,j),i=1,14),j=9,15) /
300  & 70.00, 0.645,0.466,0.243,0.186,0.132,0.130,0.136,0.148,0.160,
301  & 0.172,0.185,0.191,0.200,
302  & 80.00, 0.827,0.600,0.288,0.220,0.147,0.142,0.150,0.166,0.180,
303  & 0.191,0.203,0.209,0.218,
304  & 90.00, 1.000,0.755,0.333,0.247,0.159,0.155,0.166,0.185,0.197,
305  & 0.209,0.220,0.227,0.236,
306  & 100.0, 1.125,0.863,0.375,0.264,0.167,0.166,0.183,0.202,0.214,
307  & 0.225,0.238,0.245,0.255,
308  & 120.0, 1.260,0.983,0.450,0.281,0.180,0.188,0.215,0.234,0.247,
309  & 0.260,0.273,0.282,0.291,
310  & 150.0, 1.335,1.060,0.536,0.289,0.189,0.214,0.251,0.272,0.297,
311  & 0.312,0.325,0.336,0.346,
312  & 180.0, 1.350,1.100,0.600,0.290,0.190,0.225,0.280,0.305,0.347,
313  & 0.364,0.378,0.390,0.400/
314 !
315  REAL*8 kre(22,4)
316  DATA kre /
317  & 22.004,1.e+3,2.e+3,3.e+3,4.e+3,5.e+3,6.e+3,7.e+3,8.e+3,9.e+3,
318  & 1.e+4,2.e+4,3.e+4,4.e+4,6.e+4,8.e+4,1.e+5,2.e+5,3.e+5,
319  & 5.e+5,7.e+5,1.e+6,
320  & 1.0, 3.88,3.06,2.77,2.60,2.49,2.40,2.33,2.27,2.22,2.18,
321  & 1.86,1.69,1.57,1.41,1.30,1.22,5*1.00,
322  & 1.5, 3.88,3.06,2.77,2.60,2.49,2.40,2.33,2.27,2.22,2.18,
323  & 1.90,1.76,1.67,1.54,1.46,1.40,1.22,1.12,3*1.00,
324  & 2.0, 3.88,3.06,2.77,2.60,2.49,2.40,2.33,2.27,2.22,2.18,
325  & 1.93,1.80,1.71,1.60,1.53,1.47,1.32,1.23,1.13,1.06,1.00/
326 !
327  integer iexp6(2)
328  DATA iexp6 /0,0/
329 !
330 ! Campbell, Slattery
331 ! "Flow in the entrance of a tube"
332 ! Journal of Basic Engineering, 1963
333 !
334 ! EXIT LOSS COEFFICIENT FOR LAMINAR FLOWS DEPENDING ON THE
335 ! ACTUAL VELOCITY DISTRIBUTION AT THE EXIT
336 !
337  real*8 xdre(12)
338  DATA xdre /
339  & 0.000,0.001,0.0035,0.0065,0.010,0.0150,0.020,
340  & 0.025,0.035,0.045,0.056,0.065/
341 !
342  real*8 zetaex(12)
343  DATA zetaex /
344  & 1.00,1.200,1.40,1.54,1.63,1.73,1.80,1.85,1.93,
345  & 1.97,2.00,2.00/
346 !
347 ! Branch Joint Genium
348 ! Branching Flow Part IV - TEES
349 ! Fluid Flow Division
350 ! Section 404.2 page 4 December 1986
351 ! Genium Publishing (see www.genium.com)
352 !
353 ! n.b: the values of this table have been scaled by a factor 64.
354 !
355  real*8 xang(11),yang(11)
356  data (xang(i),yang(i),i=1,11)
357  & /0.0d0,62.d0,
358  & 15.d0,62.d0,
359  & 30.d0,61.d0,
360  & 45.d0,61.d0,
361  & 60.d0,58.d0,
362  & 75.d0,52.d0,
363  & 90.d0,40.d0,
364  & 105.d0,36.d0,
365  & 120.d0,34.d0,
366  & 135.d0,33.d0,
367  & 150.d0,32.5d0/
368 !
369 ! Branch Joint Idelchik 1
370 ! Diagrams of resistance coefficients
371 ! I.E. IDEL'CHIK 'HANDBOOK OF HYDRAULIC RESISTANCE'
372 ! 2nd edition 1986,HEMISPHERE PUBLISHING CORP.
373 ! ISBN 0-899116-284-4
374 !
375  real*8 ta2a0(12),tafakt(12)
376  data (ta2a0(i),tafakt(i),i=1,12)
377  & /0.d0 ,1.d0 ,
378  & 0.16d0 ,1.d0 ,
379  & 0.20d0 ,0.99d0,
380  & 0.25d0 ,0.95d0,
381  & 0.29d0 ,0.90d0,
382  & 0.31d0 ,0.85d0,
383  & 0.33d0 ,0.80d0,
384  & 0.35d0 ,0.78d0,
385  & 0.4d0 ,0.75d0,
386  & 0.6d0 ,0.70d0,
387  & 0.8d0 ,0.65d0,
388  & 1.d0 ,0.60d0/
389 !
390 ! Branch Joint Idelchik 2
391 ! Diagrams of resistance coefficients p348-351 section VII
392 ! I.E. IDEL'CHIK 'HANDBOOK OF HYDRAULIC RESISTANCE'
393 ! 2nd edition 1986,HEMISPHERE PUBLISHING CORP.
394 ! ISBN 0-899116-284-4
395 !
396 ! page 352 diagram 7-9 - alpha | Fs/Fc
397 !
398  real*8 kbtab(6,7),kstab(6,6)
399  data ((kbtab(i,j),j=1,7),i=1,6)
400  & /6.007d0 ,0.d0,15.d0,30.d0,45.d0,60.d0 ,90.d0 ,
401  & 0.d0 ,0.d0, 0.d0, 0.d0, 0.d0, 0.d0 , 0.d0 ,
402  & 0.1d0 ,0.d0, 0.d0, 0.d0, 0.d0, 0.d0 , 0.d0 ,
403  & 0.2d0 ,0.d0, 0.d0, 0.d0, 0.d0, 0.d0 , 0.1d0 ,
404  & 0.33d0,0.d0, 0.d0, 0.d0, 0.d0, 0.d0 , 0.2d0 ,
405  & 0.5d0 ,0.d0, 0.d0, 0.d0, 0.d0, 0.1d0 , 0.25d0/
406 !
407 ! page 348-351 diagrams 7-5 to 7-8 - alpha | Fs/Fc
408 !
409  data ((kstab(i,j),j=1,6),i=1,6)
410  & /6.006d0 ,0.d0,15.d0 ,30.d0 ,45.d0 , 60.d0 ,
411  & 0.d0 ,0.d0, 0.d0 , 0.d0 , 0.d0 , 0.d0 ,
412  & 0.1d0 ,0.d0, 0.d0 , 0.d0 , 0.05d0, 0.d0 ,
413  & 0.2d0 ,0.d0, 0.d0 , 0.d0 , 0.14d0, 0.d0 ,
414  & 0.33d0,0.d0, 0.14d0, 0.17d0, 0.14d0, 0.1d0 ,
415  & 0.5d0 ,0.d0, 0.4d0 , 0.4d0 , 0.3d0 , 0.25d0/
416 !
417 ! page 352 diagram 7-9 R: zeta_c,st
418 !
419  real*8 z90tab(6,13)
420  data ((z90tab(i,j),j=1,13),i=1,6)/
421  &6.013,0. ,0.03,0.05,0.1 ,0.2 ,0.3 ,0.4 ,0.5 ,0.6 ,0.7 ,0.8 ,1.0 ,
422  & .06, .02, .05, .08, .08, .07, .01,-.15,1.e9,1.e9,1.e9,1.e9,1.e9,
423  & .10, .04, .08, .10, .20, .26, .20, .05,-.13,1.e9,1.e9,1.e9,1.e9,
424  & .20, .08, .12, .18, .25, .34, .32, .26, .16, .02,-.14,1.e9,1.e9,
425  & .33, .45, .50, .52, .59, .66, .64, .62, .58, .44, .27, .08,-.34,
426  & .50,1.00,1.04,1.06,1.16,1.25,1.25,1.22,1.10, .88, .70, .45,0. /
427 !
428 ! table to check the location of V2V0 in Z90TAB
429 !
430  real*8 z90limx (5),z90limy(5)
431  data z90limx
432  & /0.06d0,0.1d0,0.2d0,0.33,0.5d0 /
433 !
434  data z90limy
435  & /0.1d0,0.1d0,0.3d0,0.5d0,0.7d0/
436 !
437  data n0 /0/
438  data n1 /1/
439  data n2 /2/
440  data n5 /5/
441  data n6 /6/
442  data n8 /8/
443  data n10 /10/
444  data n11 /11/
445  data n12 /12/
446  data n14 /14/
447  data n15 /15/
448  data n22 /22/
449 !
450  pi=4.d0*datan(1.d0)
451  index=ielprop(nelem)
452 !
453  if((lakon(nelem)(2:5).eq.'REUS').or.
454  & (lakon(nelem)(2:5).eq.'LPUS')) then
455 !
456 ! user defined zeta
457 !
458  zeta=prop(index+4)
459 !
460  return
461 !
462  elseif((lakon(nelem)(2:5).eq.'REEN').or.
463  & (lakon(nelem)(2:5).eq.'LPEN')) then
464 !
465 ! entrance
466 !
467  zeta=prop(index+4)
468 !
469  return
470 !
471  elseif((lakon(nelem)(2:7).eq.'RELOID').or.
472  & (lakon(nelem)(2:7).eq.'LPLOID')) then
473 !
474 ! THICK EDGED ORIFICE IN STRAIGHT CONDUIT (L/DH > 0.015)
475 ! I.E. IDEL'CHIK p175
476 !
477 ! Input parameters
478 !
479 ! Inlet/outlet sections
480  a1=prop(index+1)
481  a2=prop(index+2)
482 ! Hydraulic diameter
483  dh=prop(index+3)
484  if((dh.eq.0.d0).and.(a1.le.a2)) then
485  dh=dsqrt(4d0*a1/pi)
486  elseif((dh.eq.0.d0).and.(a1.gt.a2)) then
487  dh=dsqrt(4d0*a2/pi)
488  endif
489 ! Length
490  dl=prop(index+4)
491 !
492  lzd=dl/dh
493  a2za1=min(a1/a2, 1.)
494 !
495  fa2za1=1.d0-a2za1
496 !
497  iwrit1= 0
498  if( lzd.gt.2.4 ) iwrit1= 1
499 !
500  ldumm=1.d0
501  dhdumm=-1.d0
502  ks=0.d0
503  form_fact=1.d0
504 !
505  call friction_coefficient(ldumm,dhdumm,ks,reynolds,
506  & form_fact,lambda)
507 !
508  call onedint(xlzd,ytor,n10,lzd,thau,n1,n1,n0,ier)
509  zeta0=((0.5+thau*dsqrt(fa2za1))+fa2za1) * fa2za1
510 !
511  if(reynolds.gt.1.e+05 ) then
512  zeta=zeta0 + lambda * dabs(lzd)
513  else
514  call onedint(xre,yere,n14,reynolds,ereo,n1,n1,n0,ier)
515 !
516  call twodint(zzeta,n15,n11,reynolds,
517  & a2za1,zetap,n1,iexp,ier)
518  zeta=zetap + ereo * zeta0 + lambda * dabs(lzd)
519  IF( a2za1.gt.0.95 ) iwrit1=1
520  endif
521 !
522  if(dabs(lzd) .le. 0.015 )then
523  write(*,*) '*WARNING in zeta_calc: L/DH outside valid'
524  write(*,*) ' range ie less than 0.015 !'
525  endif
526 !
527  if( iwrit1.eq.1 ) then
528  write(*,*)
529  & 'WARNING in zeta_calc: geometry data outside valid range'
530  write(*,*)
531  & ' l/dh greater than 2.4- extrapolated value(s) !'
532  endif
533 !
534  elseif((lakon(nelem)(2:7).eq.'REWAOR').or.
535  & (lakon(nelem)(2:7).eq.'LPWAOR'))then
536 !
537 ! THICK-WALLED ORIFICE IN LARGE WALL (L/DH > 0.015)
538 ! I.E. IDL'CHIK page 174
539 !
540 ! Input parameters
541 !
542 ! Inlet/outlet sections
543  a1=prop(index+1)
544  a2=prop(index+2)
545 ! Hydraulic diameter
546  dh=prop(index+3)
547  if((dh.eq.0.d0).and.(a1.le.a2)) then
548  dh=dsqrt(4d0*a1/pi)
549  elseif((dh.eq.0.d0).and.(a1.gt.a2)) then
550  dh=dsqrt(4d0*a2/pi)
551  endif
552 ! Length
553  dl=prop(index+4)
554 !
555  lzd=dl/dh
556  ldumm=1.d0
557  dhdumm=-1.d0
558  ks=0.d0
559  form_fact=1.d0
560 !
561  call friction_coefficient(ldumm,dhdumm,ks,reynolds,
562  & form_fact,lambda)
563  call onedint (xlqd,yzeta1,n12,lzd,zeta01,n1,n1,n0,ier)
564 !
565  iwrit1=0
566  if(lzd.gt.10.) iwrit1=1
567 !
568  if(reynolds.le.1.e+05) then
569 !
570  call onedint (xre2,yzeta2,n14,reynolds,zeta02,n1,n1,n10,ier)
571  call onedint (xre2,yere2,n14,reynolds,ereo,n1,n1,n0,ier)
572 !
573  zeta=zeta02+0.342*ereo*zeta01+lambda*lzd
574 !
575  elseif(reynolds.gt.1.e+05) then
576  zeta=zeta01+lambda*lzd
577  endif
578  if(lzd.le.0.015) then
579  write(*,*) '*WARNING in zeta_calc'
580  write(*,*)
581  & ' l/dh outside valid range i.e. less than 0.015 !'
582  endif
583  if(iwrit1.eq.1) then
584  write(*,*) '*WARNING in zeta_calc :extrapolated value(s)!'
585  endif
586 !
587  return
588 !
589  elseif((lakon(nelem)(2:5).eq.'REEL').or.
590  & (lakon(nelem)(2:5).eq.'LPEL')) then
591 !
592 ! SUDDEN EXPANSION OF A STREAM WITH UNIFORM VELOCITY DISTRIBUTION
593 ! I.E. IDL'CHIK page 160
594 !
595 ! Input parameters
596 !
597 ! Inlet/outlet sections
598  a1=prop(index+1)
599  a2=prop(index+2)
600 !
601  a2za1=a1/a2
602  iwrit1=0
603 !
604  if(reynolds.LE.10.) then
605  zeta=26.0/reynolds
606  elseif(reynolds.gt.10.and.reynolds.le.3.5e+03) then
607  call twodint(zzeta3,n14,n11,reynolds,a2za1,zeta,n1,iexp3,
608  & ier)
609  if(a2za1.lt.0.01.or.a2za1.gt.0.6) iwrit1=1
610  else
611  zeta=(1.-a2za1)**2
612  endif
613 !
614  if(iwrit1.eq.1) then
615  write(*,*) '*WARNING in zeta_calc: extrapolated value(s)!'
616  endif
617  return
618 !
619  elseif((lakon(nelem)(2:5).eq.'RECO').or.
620  & (lakon(nelem)(2:5).eq.'LPCO'))then
621 !
622 ! SUDDEN CONTRACTION WITH & WITHOUT CONICAL BELLMOUTH ENTRY
623 ! I.E. IDL'CHIK p 168
624 !
625 ! Input parameters
626 !
627 ! Inlet/outlet sections
628  a1=prop(index+1)
629  a2=prop(index+2)
630 ! Hydraulic diameter
631  dh=prop(index+3)
632  if((dh.eq.0.d0).and.(a1.le.a2)) then
633  dh=dsqrt(4d0*a1/pi)
634  elseif((dh.eq.0.d0).and.(a1.gt.a2)) then
635  dh=dsqrt(4d0*a2/pi)
636  endif
637 ! Length
638  dl=prop(index+4)
639 ! Angle
640  alpha=prop(index+5)
641 !
642  a2za1=a2/a1
643  iwrit1=0
644  dl=abs(dl)
645  lzd=dl/dh
646 !
647  if(dl.eq.0.d0) then
648  if(reynolds.le.10.) then
649  zeta=27.0/reynolds
650  elseif(reynolds.gt.10.and.reynolds.le.1.e+04) then
651  call twodint(zzeta41,n14,n11,reynolds,a2za1,zeta,n1,iexp,
652  & ier)
653  if(a2za1.le.0.1.or.a2za1.gt.0.6) iwrit1=1
654  elseif(reynolds.gt.1.e+04) then
655  zeta=0.5*(1.-a2za1)
656  endif
657  elseif(dl.gt.0.d0) then
658  call twodint(zzeta42,n10,n0,alpha,lzd,zeta0,n1,iexp,ier)
659  zeta=zeta0*(1.-a2za1)
660  if(lzd.lt.0.025 .or. lzd.gt.0.6) iwrit1=1
661  if(reynolds .le. 1.e+04) then
662  write(*,*) '*WARNING in zeta_calc: reynolds outside valid
663  & range i.e. < 10 000 !'
664  endif
665  endif
666 !
667  if( iwrit1.eq.1 ) then
668  WRITE(*,*) '*WARNING in zeta_calc: extrapolierte Werte!'
669  endif
670 !
671  return
672 !
673  elseif((lakon(nelem)(2:7).eq.'REBEID').or.
674  & (lakon(nelem)(2:7).eq.'LPBEID')) then
675 !
676 !
677 ! SHARP ELBOW (R/DH=0) AT 0 < DELTA < 180
678 ! I.E. IDL'CHIK page 294
679 !
680 ! SHARP BENDS 0.5 < R/DH < 1.5 AND 0 < DELTA < 180
681 ! I.E. IDL'CHIK page 289-290
682 !
683 ! SMOOTH BENDS (R/DH > 1.5) AT 0 < DELTA < 180
684 ! I.E. IDL'CHIK page 289-290
685 !
686 ! Input parameters
687 !
688 ! Inlet/outlet sections
689  a1=prop(index+1)
690  a2=prop(index+2)
691 ! Hydraulic diameter
692  dh=prop(index+3)
693  if((dh.eq.0.d0).and.(a1.le.a2)) then
694  dh=dsqrt(4d0*a1/pi)
695  elseif((dh.eq.0.d0).and.(a1.gt.a2)) then
696  dh=dsqrt(4d0*a2/pi)
697  endif
698 ! radius
699  rad=prop(index+4)
700 ! angle
701  delta=prop(index+5)
702 ! heigth/width (square section)
703  a0=prop(index+6)
704  b0=prop(index+7)
705 !
706  iwrit1=0
707  iwrit2=0
708  rzdh=rad/dh
709  if(a0.eq.0.) azb=1.0
710  if(a0.gt.0.) azb=a0/b0
711 !
712  if(rzdh.le.0.5) then
713  call onedint(xaqb,yc,n12,azb,c,n1,n1,n0,ier)
714  zeta1=0.95*(sin(delta*0.0087))**2+2.05*(sin(delta*0.0087))**4
715  call onedint(xdelta,ya,n10,delta,a,n1,n1,n10,ier)
716  zeta=c*a*zeta1
717  if(azb.le.0.25.or.azb.gt.8.0) iwrit2=1
718  if(reynolds.lt.4.e+04) then
719  if(reynolds.le.3.e+03) iwrit1=1
720  rei=max(2999.,reynolds)
721  ldumm=1.d0
722  dhdumm=-1.d0
723  ks=0.d0
724  form_fact=1.d0
725  call friction_coefficient(ldumm,dhdumm,ks,rei,form_fact
726  & ,lambda)
727  re_val=4.e+04
728  call friction_coefficient(ldumm,dhdumm,ks,re_val,form_fact
729  & , lam)
730  zeta=zeta*lambda/lam
731  endif
732 !
733  elseif(rzdh.gt.0.5.and.rzdh.lt.1.5) then
734  call onedint(xdelta,ya1,n10,delta,ai,n1,n1,n10,ier)
735  call onedint(xrqdh,yb1,n8,rzdh,b1,n1,n1,n10,ier)
736  call onedint(xaqb,yc1,n12,azb,c1,n1,n1,n10,ier)
737  rei=max(2.e5,reynolds)
738  ldumm=1.d0
739  dhdumm=-1.d0
740  ks=0.d0
741  form_fact=1.d0
742  call friction_coefficient(ldumm,dhdumm,ks,rei,form_fact
743  & , lambda)
744  zeta=ai*b1*c1+0.0175*delta*rzdh*lambda
745  if(azb.lt.0.25.or.azb.gt.8.0) iwrit2=1
746  if(reynolds.lt.2.e+05) then
747  IF(reynolds.lt.3.e+03) iwrit1=1
748  rei=max(2999.,reynolds)
749  call friction_coefficient(ldumm,dhdumm,ks,rei,form_fact
750  & ,lambda)
751  re_val=2.e+05
752  call friction_coefficient(ldumm,dhdumm,ks,re_val,form_fact
753  & , lam)
754  zeta=zeta*lambda/lam
755  endif
756 !
757  elseif(rzdh.ge.1.5.and.rzdh.lt.50.) then
758  call onedint(xdelta,ya1,n10,delta,ai,n1,n1,n10,ier)
759  call onedint(xaqb,yc2,n12,azb,c2,n1,n1,n10,ier)
760  call onedint(xrzdh,yb2,n8,rzdh,b2,n1,n1,n0,ier)
761  rei=max(2.e5,reynolds)
762  ldumm=1.d0
763  dhdumm=-1.d0
764  ks=0.d0
765  form_fact=1.d0
766  call friction_coefficient(ldumm,dhdumm,ks,rei,form_fact
767  & ,lambda)
768  zeta=ai*b2*c2+0.0175*delta*rzdh*lambda
769  if(azb.lt.0.25.or.azb.gt.8.0) iwrit2=1
770  if(reynolds.lt.2.e+05) then
771  if(reynolds.lt.3.e+03) iwrit1=1
772  rei=max(2999.,reynolds)
773  call friction_coefficient(ldumm,dhdumm,ks,rei,form_fact
774  & ,lambda)
775  re_val=2.e+05
776  call friction_coefficient(ldumm,dhdumm,ks,re_val,form_fact
777  & , lam)
778  zeta=zeta*lambda/lam
779  endif
780 !
781  elseif(rzdh.ge.50.) then
782  zeta=0.0175*rzdh*delta*lambda
783  if(reynolds.lt.2.e+04) then
784  write (*,*)'Reynolds outside valid range i.e. < 20 000!'
785  endif
786  endif
787 !
788  if(iwrit1.eq.1) then
789 !
790  write (*,*) 'Reynolds outside valid range i.e. < 3 000!'
791  endif
792 !
793  if(iwrit2.eq.1) then
794  write(*,*) '*WARNING in zeta_calc: extrapolated value(s)!'
795  endif
796  return
797 !
798  elseif((lakon(nelem)(2:7).eq.'REBEMI').or.
799  & (lakon(nelem)(2:7).eq.'LPBEMI')) then
800 !
801 ! SMOOTH BENDS B.H.R.A HANDBOOK
802 !
803 ! Input parameters
804 !
805 ! Inlet/outlet sections
806  a1=prop(index+1)
807  a2=prop(index+2)
808 ! Hydraulic diameter
809  dh=prop(index+3)
810 ! Radius:
811  rad=prop(index+4)
812 ! angle delta:
813  delta=prop(index+5)
814 !
815  rzdh=rad/dh
816 !
817  iwrit1=0
818  if( delta.lt.10. .or. delta.gt.180. .or.
819  & rzdh .lt.0.5 .or. rzdh. gt. 10. ) iwrit1=1
820 !
821  call twodint(zzetao,n14,n11,rzdh,delta,zeta0,n1,iexp6,ier)
822  call twodint(kre, n22,n11,reynolds,rzdh, k,n1,iexp6,ier)
823  zeta=zeta0 * k
824 !
825  if( reynolds.lt.1.e+3 .or. reynolds.gt.1.e+6 ) then
826  write (*,*)'Reynolds outside valid range <1.E+3 or >1.0E+6'
827  endif
828 !
829  if( iwrit1.eq.1 ) then
830  write (*,*)': geometry data outside valid range '
831  write (*,*)' - extrapolated value(s)!'
832  endif
833  RETURN
834 !
835  elseif((lakon(nelem)(2:7).eq.'REBEMA').or.
836  & (lakon(nelem)(2:7).eq.'LPBEMA')) then
837 !
838 ! Own tables and formula to be included
839 !
840  Write(*,*) '*WARNING in zeta_calc: ZETA implicitly equal 1'
841  zeta=1.d0
842 
843  RETURN
844 !
845  elseif((lakon(nelem)(2:5).eq.'REEX').or.
846  & (lakon(nelem)(2:5).eq.'LPEX')) then
847 !
848 ! EXIT LOSS COEFFICIENT FOR LAMINAR FLOWS DEPENDING ON THE
849 ! ACTUAL VELOCITY DISTRIBUTION AT THE EXIT
850 !
851 ! Input parameters
852 !
853 ! Inlet/outlet sections
854  a1=prop(index+1)
855  a2=prop(index+2)
856 ! Hydraulic diameter
857  dh=prop(index+3)
858  if((dh.eq.0.d0).and.(a1.le.a2)) then
859  dh=dsqrt(4d0*a1/pi)
860  elseif((dh.eq.0.d0).and.(a1.gt.a2)) then
861  dh=dsqrt(4d0*a2/pi)
862  endif
863 ! Reference element
864  nelem_ref=nint(prop(index+4))
865 !
866  if(lakon(nelem_ref)(2:5).ne.'GAPF') then
867  write(*,*) '*ERROR in zeta_calc :the reference element is no
868  &t of type GASPIPE'
869  call exit(201)
870  endif
871 !
872  if(lakon(nelem_ref)(2:6).eq.'GAPFI') then
873  isothermal=.true.
874  endif
875 ! Length of the previous pipe element
876  dl=abs(prop(ielprop(nelem_ref)+3))
877 !
878  if(reynolds .le. 2300.) then
879 ! (LAMINAR FLOW)
880  ldre=dl/dh/reynolds
881  call onedint (xdre,zetaex,n12,ldre,zeta,n1,n1,n0,ier)
882  elseif((reynolds.gt.2300).and.(reynolds.lt.3000)) then
883 ! (TRANSITION LAMINAR-TURBULENT)
884  ldre=dl/dh/2300.
885  call onedint (xdre,zetaex,n12,ldre,zetah,n1,n1,n0,ier)
886  zeta=zetah-(zetah-1.)*((reynolds-2300.)/700.)
887  else
888 ! (TURBULENT FLOW, RE.GT.3000)
889  zeta=1.
890  endif
891 !
892  RETURN
893 !
894  elseif((lakon(nelem)(2:7).eq.'RELOLI').or.
895  & (lakon(nelem)(2:7).eq.'LPLOLI')) then
896 !
897 ! 'METHOD OF LICHTAROWICZ'
898 ! "Discharge coeffcients for incompressible non-cavitating
899 ! flow through long orifices"
900 ! A. Lichtarowicz, R.K duggins and E. Markland
901 ! Journal Mechanical Engineering Science , vol 7, No. 2, 1965
902 !
903 ! TOTAL PRESSURE LOSS COEFFICIENT FOR LONG ORIFICES AND LOW REYNOLDS
904 ! NUMBERS ( RE < 2.E04 )
905 !
906 ! Input parameters
907 !
908 ! Inlet/outlet sections
909  a1=prop(index+1)
910  a2=prop(index+2)
911 ! Hydraulic diameter
912  dh=prop(index+3)
913  if((dh.eq.0.d0).and.(a1.le.a2)) then
914  dh=dsqrt(4d0*a1/pi)
915  elseif((dh.eq.0.d0).and.(a1.gt.a2)) then
916  dh=dsqrt(4d0*a2/pi)
917  endif
918 ! Length
919  dl=prop(index+4)
920 ! Isotermal
921 !
922  lzd=dabs(dl)/dh
923 !
924  cdu=0.827-0.0085*lzd
925  km=a1/a2
926  call cd_lichtarowicz(cd,cdu,reynolds,km,lzd)
927  if(reynolds.gt.2.e04) then
928  write(*,*)
929  & '*WARNING in zeta_calc: range of application exceeded !'
930  endif
931 !
932  zeta=1./cd**2
933 !
934  return
935 !
936 ! Branch
937 !
938  elseif((lakon(nelem)(2:5).eq.'REBR').or.
939  & (lakon(nelem)(2:5).eq.'LPBR')) then
940  nelem0=nint(prop(index+1))
941  nelem1=nint(prop(index+2))
942  nelem2=nint(prop(index+3))
943  a0=prop(index+4)
944  a1=prop(index+5)
945  a2=prop(index+6)
946  alpha1=prop(index+7)
947  alpha2=prop(index+8)
948 !
949 ! node definition
950 !
951  node10=kon(ipkon(nelem0)+1)
952  node20=kon(ipkon(nelem0)+3)
953  nodem0=kon(ipkon(nelem0)+2)
954 !
955  node11=kon(ipkon(nelem1)+1)
956  node21=kon(ipkon(nelem1)+3)
957  nodem1=kon(ipkon(nelem1)+2)
958 !
959  node12=kon(ipkon(nelem2)+1)
960  node22=kon(ipkon(nelem2)+3)
961  nodem2=kon(ipkon(nelem2)+2)
962 !
963 ! determining the nodes which are not in common
964 !
965  if(node10.eq.node11) then
966  node0=node10
967  node1=node21
968  if(node11.eq.node12) then
969  node2=node22
970  elseif(node11.eq.node22) then
971  node2=node12
972  endif
973  elseif(node10.eq.node21) then
974  node0=node10
975  node1=node11
976  if(node21.eq.node12) then
977  node0=node22
978  elseif(node21.eq.node22) then
979  node2=node12
980  endif
981  elseif(node20.eq.node11) then
982  node0=node20
983  node1=node21
984  if(node11.eq.node12) then
985  node2=node22
986  elseif(node11.eq.node22) then
987  node2=node12
988  endif
989  elseif(node20.eq.node21) then
990  node0=node20
991  node1=node11
992  if(node11.eq.node21) then
993  node2=node22
994  elseif(node21.eq.node22) then
995  node2=node12
996  endif
997  endif
998 !
999 ! density
1000 !
1001  if(lakon(nelem)(2:3).eq.'RE') then
1002 !
1003 ! for gases
1004 !
1005  qred_crit=dsqrt(kappa/r)*
1006  & (1+0.5d0*(kappa-1))**(-0.5d0*(kappa+1)/(kappa-1))
1007 !
1008  icase=0
1009 !
1010  tt0=v(0,node0)
1011  xflow0=iaxial*v(1,nodem0)
1012 c xflow0=v(1,nodem0)
1013  pt0=v(2,node0)
1014 !
1015  qred_0=dabs(xflow0)*dsqrt(tt0)/(a0*pt0)
1016  if(qred_0.gt.qred_crit)
1017  & then
1018  xflow0=qred_crit*(a0*pt0)/dsqrt(tt0)
1019  endif
1020 !
1021  call ts_calc(xflow0,tt0,pt0,kappa,r,a0,ts0,icase)
1022  m0=dsqrt(2/(kappa-1)*(tt0/ts0-1))
1023 !
1024  rho0=pt0/(r*tt0)*(tt0/ts0)**(-1/(kappa-1))
1025 !
1026  tt1=v(0,node1)
1027  xflow1=iaxial*v(1,nodem1)
1028 c xflow1=v(1,nodem1)
1029  pt1=v(2,node0)
1030 !
1031  qred_1=dabs(xflow1)*dsqrt(tt1)/(a1*pt1)
1032  if(qred_1.gt.qred_crit)
1033  & then
1034  xflow1=qred_crit*(a1*pt1)/dsqrt(tt1)
1035  endif
1036 !
1037  call ts_calc(xflow1,tt1,pt1,kappa,r,a1,ts1,icase)
1038  m1=dsqrt(2/(kappa-1)*(tt1/ts1-1))
1039 !
1040  rho1=pt1/(r*tt1)*(tt1/ts1)**(-1/(kappa-1))
1041 !
1042  tt2=v(0,node2)
1043  xflow2=iaxial*v(1,nodem2)
1044 c xflow2=v(1,nodem2)
1045  pt2=v(2,node0)
1046 !
1047  qred_2=dabs(xflow2)*dsqrt(tt2)/(a2*pt2)
1048  if(qred_2.gt.qred_crit) then
1049  xflow2=qred_crit*(a2*pt2)/dsqrt(tt2)
1050  endif
1051 !
1052  call ts_calc(xflow2,tt2,pt2,kappa,r,a2,ts2,icase)
1053  m2=dsqrt(2/(kappa-1)*(tt2/ts2-1))
1054  rho2=pt2/(r*tt2)*(tt2/ts2)**(-1/(kappa-1))
1055  else
1056 !
1057 ! for liquids the density is supposed to be constant
1058 ! across the element
1059 !
1060  rho0=1.d0
1061  rho1=1.d0
1062  rho2=1.d0
1063  endif
1064 !
1065 ! volumic flows (positive)
1066 !
1067  v0=dabs(xflow0/rho0)
1068  v1=dabs(xflow1/rho1)
1069  v2=dabs(xflow2/rho2)
1070 !
1071  v1v0=v1/v0
1072  v2v0=v2/v0
1073 !
1074  a0a1=a0/a1
1075  a0a2=a0/a2
1076  a2a0=1/a0a2
1077 !
1078  w0w1=1/(v1v0*a0a1)
1079  w0w2=1/(v2v0*a0a2)
1080 !
1081 ! Branch Joint Genium
1082 ! Branching Flow Part IV - TEES
1083 ! Fluid Flow Division
1084 ! Section 404.2 page 4 December 1986
1085 ! Genium Publishing (see www.genium.com)
1086 !
1087  if((lakon(nelem)(2:7).eq.'REBRJG').or.
1088  & (lakon(nelem)(2:7).eq.'LPBRJG')) then
1089 !
1090  ang1s=(1.41d0-0.00594*alpha1)*alpha1*pi/180
1091  ang2s=(1.41d0-0.00594*alpha2)*alpha2*pi/180
1092 !
1093  cang1s=dcos(ang1s)
1094  cang2s=dcos(ang2s)
1095 !
1096 ! linear part
1097 !
1098  zetlin=2.d0*(v1v0**2*a0a1*cang1s+v2v0**2*a0a2*cang2s)
1099 !
1100  if(nelem.eq.nelem1) then
1101  call onedint(xang,yang,n11,alpha1,lam10,n1,n2,n22,ier)
1102  zeta=lam10/64*(v1v0*a0a1)**2-zetlin+1d0
1103  zeta=zeta*(w0w1)**2
1104 !
1105  elseif(nelem.eq.nelem2) then
1106  call onedint(xang,yang,n11,alpha2,lam20,n1,n2,n22,ier)
1107  zeta=lam20/64*(v2v0*a0a2)**2-zetlin+1d0
1108  zeta=zeta*(w0w2)**2
1109  endif
1110  if(zeta.lt.0) then
1111  write(*,*) '*WARNING in zeta_calc: in Element',nelem,
1112  & 'TYPE= ',lakon(nelem)(2:5)
1113  write(*,*) ' Zeta value negative is set to 0.01'
1114  zeta=0.01d0
1115  endif
1116  return
1117 !
1118  elseif((lakon(nelem)(2:8).eq.'REBRJI1').or.
1119  & (lakon(nelem)(2:8).eq.'LPBRJI1')) then
1120 !
1121 ! Branch Joint Idelchik 1
1122 ! Diagrams of resistance coefficients p260-p266 section VII
1123 ! I.E. IDEL'CHIK 'HANDBOOK OF HYDRAULIC RESISTANCE'
1124 ! 2nd edition 1986,HEMISPHERE PUBLISHING CORP.
1125 ! ISBN 0-899116-284-4
1126 !
1127  a0a2=a0/a2
1128  if(alpha2.lt.60.d0) then
1129  if(nelem.eq.nelem1) then
1130  zeta=1.d0-v1v0**2
1131  & -2.d0*a0a2*v2v0**2*dcos(alpha2*pi/180)
1132  zeta=zeta*(w0w1)**2
1133  elseif(nelem.eq.nelem2) then
1134  zeta=1.d0-v1v0**2
1135  & -2.d0*a0a2*v2v0**2*dcos(alpha2*pi/180)
1136  & +(a0a2*v2v0)**2-v1v0**2
1137  zeta=zeta*(w0w2)**2
1138  endif
1139 !
1140  elseif(alpha2.eq.60.d0) then
1141 !
1142 ! proceeding as for alpha2<60 with cos(alpha2)=0.5
1143 !
1144  if(nelem.eq.nelem1) then
1145  zeta=1.d0-v1v0**2-a0a2*v2v0**2
1146  zeta=zeta*(w0w1)**2
1147  elseif(nelem.eq.nelem2) then
1148  zeta=1.d0-v1v0**2-a0a2*v2v0**2
1149  & +(a0a2*v2v0)**2-v1v0**2
1150  zeta=zeta*(w0w2)**2
1151  endif
1152 !
1153  elseif(alpha2.lt.90.d0) then
1154 !
1155 ! linear interpolation between alpha2=60 and alpha2=90
1156 !
1157  z1_60=1.d0-v1v0**2-a0a2*v2v0**2
1158  z1_90=(1.55d0-v2v0)*v2v0
1159  if(nelem.eq.nelem1) then
1160  zeta=z1_60+(z1_90-z1_60)*(alpha2-60.d0)/30
1161  zeta=zeta*(w0w1)**2
1162  elseif(nelem.eq.nelem2) then
1163  z2_60=z1_60+(a0a2*v2v0)**2-v1v0**2
1164  call onedint(ta2a0,tafakt,n12,a2a0,afakt,
1165  & n1,n1,n11,ier)
1166  z2_90=afakt*(1.d0+(a0a2*v2v0)**2-2.d0*v1v0**2)
1167  zeta=z2_60+(z2_90-z2_60)*(alpha2-60.d0)/30d0
1168  zeta=zeta*(w0w2)**2
1169  endif
1170 !
1171  elseif(alpha2.eq.90.d0) then
1172  if(nelem.eq.nelem1) then
1173  zeta=(1.55d0-v2v0)*v2v0
1174  zeta=zeta*(w0w1)**2
1175  elseif(nelem.eq.nelem2) then
1176  call onedint(ta2a0,tafakt,n12,a2a0,afakt,
1177  & n1,n1,n11,ier)
1178  zeta=afakt*(1.d0+(a0a2*v2v0)**2-2.d0*v1v0**2)
1179  zeta=zeta*(w0w2)**2
1180  endif
1181  endif
1182  if(zeta.lt.0) then
1183  write(*,*) '*WARNING in zeta_calc: in Element',nelem,
1184  & 'TYPE= ',lakon(nelem)(2:5)
1185  write(*,*) ' Zeta value negative is set to 0.01'
1186  zeta=0.01d0
1187  endif
1188  return
1189 !
1190  elseif((lakon(nelem)(2:8).eq.'REBRJI2').or.
1191  & (lakon(nelem)(2:8).eq.'LPBRJI2')) then
1192 !
1193 ! Branch Joint Idelchik 2
1194 ! Diagrams of resistance coefficients page 348-352
1195 ! I.E. IDEL'CHIK 'HANDBOOK OF HYDRAULIC RESISTANCE'
1196 ! 2nd edition 1986,HEMISPHERE PUBLISHING CORP.
1197 ! ISBN 0-899116-284-4 page 348-352
1198 !
1199  if(alpha2.lt.60.d0) then
1200  if(nelem.eq.nelem1) then
1201  zeta=1+a0a1*v1v0**2*(a0a1-2.)
1202  & -2d0*a0a2*v2v0**2*dcos(alpha2*pi/180)
1203 ! correction term
1204  call twodint(kstab,n6,n11,a2a0,alpha2,ks2,n1
1205  & ,iexpbr1,ier)
1206  zeta=zeta+ks2
1207  zeta=zeta*(w0w1)**2
1208  elseif(nelem.eq.nelem2) then
1209  zeta=1+a0a1*v1v0**2*(a0a1-2.)
1210  & -2d0*a0a2*v2v0**2*dcos(alpha2*pi/180)
1211  & -(a0a1*v1v0)**2+(a0a2*v2v0)**2
1212  call twodint(kbtab,n6,n11,a2a0,alpha2,kb,n1,
1213  & iexpbr1,ier)
1214  zeta=zeta+kb
1215  zeta=zeta*(w0w2)**2
1216  endif
1217 !
1218  elseif(alpha2.eq.60.d0) then
1219 ! as for alpha2 < 60 , with dcos(alpha2)=0.5
1220  if(nelem.eq.nelem1) then
1221  zeta=1+a0a1*v1v0**2*(a0a1-2.)-a0a2*v2v0**2
1222  call twodint(kstab,n6,n11,a2a0,alpha2,ks2,n1,
1223  & iexpbr1,ier)
1224  zeta=zeta+ks2
1225  zeta=zeta*(w0w1)**2
1226  elseif(nelem.eq.nelem2) then
1227  zeta=1+a0a1*v1v0**2*(a0a1-2.)-a0a2*v2v0**2
1228  & -(a0a1*v1v0)**2+(a0a2*v2v0)**2
1229  call twodint(kbtab,n6,n11,a2a0,alpha2,kb,n1,
1230  & iexpbr1,ier)
1231  zeta=zeta+kb
1232  zeta=zeta*(w0w2)**2
1233  endif
1234 !
1235  elseif(alpha2.lt.90.d0) then
1236 ! linear interpolation between alpha2=60 and alpha2=90
1237  z1_60=1+a0a1*v1v0**2*(a0a1-2.)-a0a2*v2v0**2
1238 ! correction term
1239  call twodint(kstab,n6,n11,a2a0,alpha2,ks2,n1,
1240  & iexpbr1,ier)
1241  z1_60=z1_60+ks2
1242  if(nelem.eq.nelem1) then
1243  call twodint(z90tab,n6,n11,a2a0,v2v0,z1_90,
1244  & n1,iexpbr1,ier)
1245  zeta=z1_60+(z1_90-z1_60)*(alpha2-60)/30
1246  zeta=zeta*(w0w1)**2
1247  elseif(nelem.eq.nelem2) then
1248  z2_60=z1_60-(a0a1*v1v0)**2+(a0a2*v2v0)**2
1249  call twodint(kbtab,n6,n11,a2a0,alpha2,kb,n1,
1250  & iexpbr1,ier)
1251  z2_60=z2_60+kb-ks2
1252  z2_90=1.+(a0a2*v2v0)**2-2*a0a1*v1v0**2+kb
1253  zeta=z2_60+(z2_90-z2_60)*(alpha2-60)/30
1254  zeta=zeta*(w0w2)**2
1255  endif
1256  elseif(alpha2.eq.90.d0) then
1257  if(nelem.eq.nelem2) then
1258  call twodint(kbtab,n6,n11,a2a0,alpha2,kb,n1,
1259  & iexpbr1,ier)
1260  zeta=1.+(a0a2*v2v0)**2-2*a0a1*v1v0**2+kb
1261  zeta=zeta*(w0w2)**2
1262  elseif(nelem.eq.nelem1) then
1263 ! table interpolation
1264  call twodint(z90tab,n6,n11,a2a0,v2v0,zeta,
1265  & n1,iexpbr1,ier)
1266  zeta=zeta*(w0w1)**2
1267 ! cheching whether the table eveluation in the eptrapolated domain
1268 ! (This procedure is guessed from the original table)
1269 !
1270  z90lim11=z90limx(1)
1271  z90lim51=z90limx(5)
1272  if((a2a0.ge.z90lim11)
1273  & .and.(a2a0.le.z90lim51))then
1274  call onedint(z90limx,z90limy,n5,a2a0,
1275  & v2v0l,n1,n1,n11,ier)
1276  if(v2v0.gt.v2v0l) then
1277  write(*,*) 'WARNING in zeta_calc: in element',
1278  & nelem
1279  write(*,*)
1280  & ' V2V0 in the extrapolated domain'
1281  write(*,*) ' for zeta table (branch 1)'
1282  endif
1283  endif
1284  endif
1285  endif
1286  if(zeta.lt.0) then
1287  write(*,*) '*WARNING in zeta_calc: in Element',nelem,
1288  & 'TYPE= ',lakon(nelem)(2:5)
1289  write(*,*) ' Zeta value negative is set to 0.01'
1290  zeta=0.01d0
1291  endif
1292  return
1293 !
1294  elseif((lakon(nelem)(2:7).eq.'REBRSG').or.
1295  & (lakon(nelem)(2:7).eq.'LPBRSG')) then
1296 !
1297 ! Branch Split Genium
1298 ! Branching Flow Part IV - TEES
1299 ! Fluid Flow Division
1300 ! Section 404.2 page 3 December 1986
1301 ! Genium Publishing (see www.genium.com)
1302 !
1303  if(nelem.eq.nelem1) then
1304 !
1305  ang1s=(1.41d0-0.00594*alpha1)*alpha1*pi/180
1306 !
1307  cang1s=dcos(ang1s)
1308 !
1309  if(alpha1.le.22.5) then
1310  lam11=0.0712*alpha1**0.7041+0.37
1311  lam12=0.0592*alpha1**0.7029+0.37
1312  else
1313  lam11=1.d0
1314  lam12=0.9d0
1315  endif
1316  zeta=lam11+(2.d0*lam12-lam11)*(v1v0*a0a1)**2
1317  & -2d0*lam12*v1v0*a0a1*cang1s
1318  zeta=zeta*(w0w1)**2
1319 !
1320  elseif(nelem.eq.nelem2) then
1321 !
1322  ang2s=(1.41d0-0.00594*alpha2)*alpha2*pi/180
1323 !
1324  cang2s=dcos(ang2s)
1325 !
1326  if(alpha2.le.22.5) then
1327  lam21=0.0712*alpha2**0.7041+0.37
1328  lam22=0.0592*alpha2**0.7029+0.37
1329  else
1330  lam21=1.d0
1331  lam22=0.9d0
1332  endif
1333 !
1334  zeta=lam21+(2.d0*lam22-lam21)*(v2v0*a0a2)**2
1335  & -2d0*lam22*v2v0*a0a2*cang2s
1336  zeta=zeta*(w0w2)**2
1337 !
1338  endif
1339  return
1340 !
1341  elseif((lakon(nelem)(2:8).eq.'REBRSI1').or.
1342  & (lakon(nelem)(2:8).eq.'LPBRSI1')) then
1343 !
1344 ! Branch Split Idelchik 1
1345 ! Diagrams of resistance coefficients p280,p282 section VII
1346 ! I.E. IDEL'CHIK 'HANDBOOK OF HYDRAULIC RESISTANCE'
1347 ! 2nd edition 1986,HEMISPHERE PUBLISHING CORP.
1348 ! ISBN 0-899116-284-4
1349 !
1350  w1w0=v1v0*a0a1
1351  w2w0=v2v0*a0a2
1352 !
1353  if(nelem.eq.nelem1) then
1354  zeta=0.4d0*(1-w1w0)**2
1355  zeta=zeta*(w0w1)**2
1356 !
1357  elseif(nelem.eq.nelem2) then
1358 !
1359  dh0=dsqrt(a0*4d0/pi)
1360  if(dh0.eq.0.d0) then
1361  dh0=dsqrt(4d0*a0/pi)
1362  endif
1363  dh2=dsqrt(a2*4d0/pi)
1364  if(dh2.eq.0.d0) then
1365  dh2=dsqrt(4d0*a2/pi)
1366  endif
1367 !
1368  hq=dh2/dh0
1369  if(alpha2.le.60.or.hq.le.2.d0/3.d0) then
1370  zeta=0.95d0*((w2w0-2d0*dcos(alpha2*pi/180))
1371  & *w2w0+1.d0)
1372  zeta=zeta*(w0w2)**2
1373  else
1374  z2d390=0.95d0*((w2w0-2d0*dcos(90.d0*pi/180))
1375  & *w2w0+1.d0)
1376  z1p090=0.95*(0.34d0+w2w0**2)
1377  z90=z2d390+(3*hq-2.d0)*(z1p090-z2d390)
1378  z60=0.95d0*((w2w0-2d0*dcos(60.d0*pi/180))
1379  & *w2w0+1.d0)
1380  zeta=z60+(alpha2/30.d0-2.d0)*(z90-z60)
1381  zeta=zeta*(w0w2)**2
1382  endif
1383  endif
1384  return
1385 !
1386  elseif((lakon(nelem)(2:8).eq.'REBRSI2').or.
1387  & (lakon(nelem)(2:8).eq.'LPBRSI2')) then
1388 !
1389 ! Branch Split Idelchik 2
1390 ! Diagrams of resistance coefficients p289,section VII
1391 ! I.E. IDEL'CHIK 'HANDBOOK OF HYDRAULIC RESISTANCE'
1392 ! 2nd edition 1986,HEMISPHERE PUBLISHING CORP.
1393 ! ISBN 0-899116-284-4
1394 !
1395  if(nelem.eq.nelem1) then
1396  w1w0=v1v0*a0a1
1397  w0w1=1/w1w0
1398  zeta=1.d0+0.3d0*w1w0**2
1399  zeta=zeta*(w0w1)**2
1400  elseif(nelem.eq.nelem2) then
1401  w2w0=v2v0*a0a2
1402  w0w2=1/w2w0
1403  zeta=1.d0+0.3d0*w2w0**2
1404  zeta=zeta*(w0w2)**2
1405  endif
1406  return
1407  endif
1408  if(zeta.lt.0) then
1409  write(*,*) '*WARNING in zeta_calc: in Element',nelem,
1410  & 'TYPE= ',lakon(nelem)(2:5)
1411  write(*,*) ' Zeta value negative is set to 0.01'
1412  zeta=0.01d0
1413  endif
1414  endif
1415 !
1416  return
max
#define max(a, b)
Definition: cascade.c:32
min
#define min(a, b)
Definition: cascade.c:31
c1
static double * c1
Definition: mafillvcompmain.c:30
ts_calc
subroutine ts_calc(xflow, Tt, Pt, kappa, r, A, Ts, icase)
Definition: ts_calc.f:20
friction_coefficient
subroutine friction_coefficient(l, d, ks, reynolds, form_fact, lambda)
Definition: friction_coefficient.f:26
twodint
subroutine twodint(T, LSP, IART, XA, YA, ZA, NA, IEXP, IER)
Definition: twodint.f:84
cd_lichtarowicz
subroutine cd_lichtarowicz(cd, cdu, reynolds, amod, bdh)
Definition: cd_lichtarowicz.f:27
v1
static double * v1
Definition: mafillsmmain_se.c:40
onedint
subroutine onedint(XE, YE, NE, XA, YA, NA, IART, IEXP, IER)
Definition: onedint.f:67
b1
static double * b1
Definition: mafillkmain.c:30
Hosted by OpenAircraft.com, (Michigan UAV, LLC)