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

Go to the source code of this file.

Functions/Subroutines

subroutine orifice (node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, cp, R, physcon, dvi, numf, set, co, vold, mi, ttime, time, iaxial)
 

Function/Subroutine Documentation

◆ orifice()

subroutine orifice ( integer  node1,
integer  node2,
integer  nodem,
integer  nelem,
character*8, dimension(*)  lakon,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
integer, dimension(0:3,*)  nactdog,
logical  identity,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
integer  iflag,
real*8, dimension(0:mi(2),*)  v,
real*8  xflow,
real*8  f,
integer, dimension(*)  nodef,
integer, dimension(*)  idirf,
real*8, dimension(*)  df,
real*8  cp,
real*8  R,
real*8, dimension(*)  physcon,
real*8  dvi,
integer  numf,
character*81, dimension(*)  set,
real*8, dimension(3,*)  co,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(*)  mi,
real*8  ttime,
real*8  time,
integer  iaxial 
)
23 !
24 ! orifice element
25 !
26 ! author: Yannick Muller
27 !
28  implicit none
29 !
30  logical identity
31  character*8 lakon(*)
32  character*81 set(*)
33 !
34  integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
35  & ielprop(*),nodef(*),idirf(*),index,iflag,
36  & inv,ipkon(*),kon(*),number,kgas,nelemswirl,
37  & nodea,nodeb,iaxial,mi(*),i,itype
38 !
39 c real*4 ofvidg
40 !
41  real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,a,d,dl,
42  & p1,p2,t1,aeff,c1,c2,c3,cd,cp,physcon(*),p2p1,km1,dvi,
43  & kp1,kdkm1,tdkp1,km1dk,x,y,ca1,cb1,ca2,cb2,dt1,alambda,
44  & rad,beta,reynolds,theta,k_phi,c2u_new,u,pi,xflow_oil,
45  & ps1pt1,uref,cd_chamf,angle,vid,cdcrit,t2,radius,
46  & initial_radius,co(3,*),vold(0:mi(2),*),offset,ttime,time,
47  & x_tab(15), y_tab(15),x_tab2(15),y_tab2(15),curve,xmach
48 !
49 c external ofvidg
50 !
51  pi=4.d0*datan(1.d0)
52  if(iflag.eq.0) then
53  identity=.true.
54 !
55  if(nactdog(2,node1).ne.0)then
56  identity=.false.
57  elseif(nactdog(2,node2).ne.0)then
58  identity=.false.
59  elseif(nactdog(1,nodem).ne.0)then
60  identity=.false.
61  endif
62 !
63  elseif(iflag.eq.1)then
64 !
65  index=ielprop(nelem)
66  kappa=(cp/(cp-r))
67  a=prop(index+1)
68  d=prop(index+2)
69  dl=prop(index+3)
70 !
71  if(lakon(nelem)(2:5).eq.'ORFL') then
72  nodea=nint(prop(index+1))
73  nodeb=nint(prop(index+2))
74  offset=prop(index+4)
75  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
76  & co(1,nodea)-vold(1,nodea))**2)-offset
77  initial_radius=dsqrt((co(1,nodeb)-co(1,nodea))**2)-offset
78  a=pi*radius**2
79  d=2*radius
80  endif
81 !
82  p1=v(2,node1)
83  p2=v(2,node2)
84  if(p1.ge.p2) then
85  inv=1
86  t1=v(0,node1)-physcon(1)
87  else
88  inv=-1
89  p1=v(2,node2)
90  p2=v(2,node1)
91  t1=v(0,node2)-physcon(1)
92  endif
93 !
94  cd=1.d0
95 !
96  p2p1=p2/p1
97  km1=kappa-1.d0
98  kp1=kappa+1.d0
99  kdkm1=kappa/km1
100  tdkp1=2.d0/kp1
101  c2=tdkp1**kdkm1
102  aeff=a*cd
103  if(p2p1.gt.c2) then
104  xflow=inv*p1*aeff*dsqrt(2.d0*kdkm1*p2p1**(2.d0/kappa)
105  & *(1.d0-p2p1**(1.d0/kdkm1))/r)/dsqrt(t1)
106  else
107  xflow=inv*p1*aeff*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
108  & dsqrt(t1)
109  endif
110 !
111  elseif(iflag.eq.2)then
112 !
113  numf=4
114  alambda=10000.d0
115  index=ielprop(nelem)
116  kappa=(cp/(cp-r))
117  a=prop(index+1)
118 !
119  p1=v(2,node1)
120  p2=v(2,node2)
121  if(p1.ge.p2) then
122  inv=1
123  xflow=v(1,nodem)*iaxial
124  t1=v(0,node1)-physcon(1)
125  nodef(1)=node1
126  nodef(2)=node1
127  nodef(3)=nodem
128  nodef(4)=node2
129  else
130  inv=-1
131  p1=v(2,node2)
132  p2=v(2,node1)
133  xflow=-v(1,nodem)*iaxial
134  t1=v(0,node2)-physcon(1)
135  nodef(1)=node2
136  nodef(2)=node2
137  nodef(3)=nodem
138  nodef(4)=node1
139  endif
140 !
141  idirf(1)=2
142  idirf(2)=0
143  idirf(3)=1
144  idirf(4)=2
145 !
146 ! calculation of the dynamic viscosity
147 !
148 !
149  if(dabs(dvi).lt.1e-30) then
150  write(*,*) '*ERROR in orifice: '
151  write(*,*) ' no dynamic viscosity defined'
152  write(*,*) ' dvi= ',dvi
153  call exit(201)
154 c kgas=0
155 c call dynamic_viscosity(kgas,T1,dvi)
156  endif
157 !
158  if((lakon(nelem)(4:5).ne.'BT').and.
159  & (lakon(nelem)(4:5).ne.'PN').and.
160  & (lakon(nelem)(4:5).ne.'C1').and.
161  & (lakon(nelem)(4:5).ne.'FL') ) then
162  d=prop(index+2)
163  dl=prop(index+3)
164 ! circumferential velocity of the rotating hole (same as disc @ given radius)
165  u=prop(index+7)
166  nelemswirl=nint(prop(index+8))
167  if(nelemswirl.eq.0) then
168  uref=0.d0
169  else
170 ! swirl generating element
171 !
172 ! preswirl nozzle
173  if(lakon(nelemswirl)(2:5).eq.'ORPN') then
174  uref=prop(ielprop(nelemswirl)+5)
175 ! rotating orifices
176  else if((lakon(nelemswirl)(2:5).eq.'ORMM').or.
177  & (lakon(nelemswirl)(2:5).eq.'ORMA').or.
178  & (lakon(nelemswirl)(2:5).eq.'ORPM').or.
179  & (lakon(nelemswirl)(2:5).eq.'ORPA')) then
180  uref=prop(ielprop(nelemswirl)+7)
181 ! forced vortex
182  elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
183  uref=prop(ielprop(nelemswirl)+7)
184 ! free vortex
185  elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
186  uref=prop(ielprop(nelemswirl)+9)
187 ! Moehring
188  elseif(lakon(nelemswirl)(2:4).eq.'MRG') then
189  uref=prop(ielprop(nelemswirl)+10)
190 ! RCAVO
191  elseif((lakon(nelemswirl)(2:4).eq.'ROR').or.
192  & (lakon(nelemswirl)(2:4).eq.'ROA'))then
193  uref=prop(ielprop(nelemswirl)+6)
194 ! RCAVI
195  elseif(lakon(nelemswirl)(2:4).eq.'RCV') then
196  uref=prop(ielprop(nelemswirl)+5)
197 !
198  else
199  write(*,*) '*ERROR in orifice:'
200  write(*,*) ' element',nelemswirl
201  write(*,*) ' refered by element',nelem
202  write(*,*) ' is not a swirl generating element'
203  endif
204  endif
205 ! write(*,*) 'nelem',nelem, u, uref
206  u=u-uref
207  angle=prop(index+5)
208 !
209  endif
210 !
211 ! calculate the discharge coefficient using Bragg's Method
212 ! "Effect of Compressibility on the discharge coefficient
213 ! of orifices and convergent nozzles"
214 ! journal of mechanical Engineering
215 ! vol2 No 1 1960
216 !
217  if((lakon(nelem)(2:5).eq.'ORBG')) then
218 !
219  p2p1=p2/p1
220  cdcrit=prop(index+2)
221 !
222  itype=2
223  call cd_bragg(cdcrit,p2p1,cd,itype)
224 !
225  elseif(lakon(nelem)(2:5).eq.'ORMA') then
226 !
227 ! calculate the discharge coefficient using own table data and
228 ! using Dr.Albers method for rotating cavities
229 !
230  call cd_own_albers(p1,p2,dl,d,cd,u,t1,r,kappa)
231 !
232 ! outlet circumferential velocity of the fluid is equal to the circumferential velocity of the hole
233 ! as the holes are perpendicular to the rotating surface and rotating with it
234 ! prop(index+7)
235 !
236 ! chamfer correction
237 !
238  if(angle.gt.0.d0)then
239  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
240  cd=cd*cd_chamf
241  endif
242 !
243  elseif(lakon(nelem)(2:5).eq.'ORMM') then
244 !
245 ! calculate the discharge coefficient using McGreehan and Schotsch method
246 !
247  rad=prop(index+4)
248 !
249  reynolds=dabs(xflow)*d/(dvi*a)
250 !
251 ! outlet circumferential velocity of the fluid is equal to the circumferential velocity of the hole
252 ! as the holes are perpendicular to the rotating surface and rotating with it
253 ! prop(index+7)
254 !
255  call cd_ms_ms(p1,p2,t1,rad,d,dl,kappa,r,reynolds,u,vid,cd)
256 !
257  if(cd.ge.1) then
258  write(*,*) ''
259  write(*,*) '**WARNING**'
260  write(*,*) 'in RESTRICTOR ',nelem
261  write(*,*) 'Calculation using'
262  write(*,*) ' McGreehan and Schotsch method:'
263  write(*,*) ' Cd=',cd,'>1 !'
264  write(*,*) 'Calcultion will proceed will Cd=1'
265  write(*,*) 'l/d=',dl/d,'r/d=',rad/d,'u/vid=',u/vid
266  cd=1.d0
267  endif
268 !
269 ! chamfer correction
270 !
271  if(angle.gt.0.d0) then
272  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
273  cd=cd*cd_chamf
274  endif
275 !
276  elseif (lakon(nelem)(2:5).eq.'ORPA') then
277 !
278 ! calculate the discharge coefficient using Parker and Kercher method
279 ! and using Dr. Albers method for rotating cavities
280 !
281  rad=prop(index+4)
282 !
283  beta=prop(index+6)
284 !
285  reynolds=dabs(xflow)*d/(dvi*a)
286 !
287  call cd_pk_albers(rad,d,dl,reynolds,p2,p1,beta,kappa,
288  & cd,u,t1,r)
289 !
290 ! outlet circumferential velocity of the fluid is equal to the circumferential velocity of the hole
291 ! as the holes are perpendicular to the rotating surface and rotating with it
292 ! prop(index+7)
293 !
294 ! chamfer correction
295 !
296  if(angle.gt.0.d0) then
297  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
298  cd=cd*cd_chamf
299  endif
300 !
301  elseif(lakon(nelem)(2:5).eq.'ORPM') then
302 !
303 ! calculate the discharge coefficient using Parker and Kercher method
304 ! and using Mac Grehan and Schotsch method for rotating cavities
305 !
306  rad=prop(index+4)
307 !
308  beta=prop(index+6)
309  reynolds=dabs(xflow)*d/(dvi*a)
310 !
311  call cd_pk_ms(rad,d,dl,reynolds,p2,p1,beta,kappa,cd,
312  & u,t1,r)
313 !
314 ! outlet circumferential velocity of the fluid is equal to the circumferential velocity of the hole
315 ! as the holes are perpendicular to the rotating surface and rotating with it
316 ! prop(index+7)
317 !
318 ! chamfer correction
319 !
320  if(angle.gt.0.d0) then
321  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
322  cd=cd*cd_chamf
323  endif
324 !
325  elseif(lakon(nelem)(2:5).eq.'ORC1') then
326 !
327  d=dsqrt(a*4/pi)
328  reynolds=dabs(xflow)*d/(dvi*a)
329  cd=1.d0
330 !
331  elseif(lakon(nelem)(2:5).eq.'ORBT') then
332 !
333 ! calculate the discharge coefficient of bleed tappings (OWN tables)
334 !
335  ps1pt1=prop(index+2)
336  curve=nint(prop(index+3))
337  number=nint(prop(index+4))
338 !
339  if(number.ne.0.d0)then
340  do i=1,number
341  x_tab(i)=prop(index+2*i+3)
342  y_tab(i)=prop(index+2*i+4)
343  enddo
344  endif
345 !
346  call cd_bleedtapping(p2,p1,ps1pt1,number,curve,x_tab,y_tab,
347  & cd)
348 !
349  elseif(lakon(nelem)(2:5).eq.'ORPN') then
350 !
351 ! calculate the discharge coefficient of preswirl nozzle (OWN tables)
352 !
353  d=dsqrt(4*a/pi)
354  reynolds=dabs(xflow)*d/(dvi*a)
355  curve=nint(prop(index+4))
356  number=nint(prop(index+6))
357  if(number.ne.0.d0)then
358  do i=1,number
359  x_tab2(i)=prop(index+2*i+5)
360  y_tab2(i)=prop(index+2*i+6)
361  enddo
362  endif
363  call cd_preswirlnozzle(p2,p1,number,curve,x_tab2,y_tab2
364  & ,cd)
365 !
366  theta=prop(index+2)
367  k_phi=prop(index+3)
368 !
369  if(p2/p1.gt.(2/(kappa+1.d0))**(kappa/(kappa-1.d0))) then
370  c2u_new=k_phi*cd*sin(theta*pi/180.d0)*r*
371  & dsqrt(2.d0*kappa/(r*(kappa-1)))*
372  & dsqrt(t1*(1.d0-(p2/p1)**((kappa-1)/kappa)))
373 !
374  else
375  c2u_new=k_phi*cd*sin(theta*pi/180.d0)*r*
376  & dsqrt(2.d0*kappa/(r*(kappa-1)))*
377  & dsqrt(t1*(1.d0-2/(kappa+1)))
378  endif
379  prop(index+5)=c2u_new
380 !
381  elseif(lakon(nelem)(2:5).eq.'ORFL') then
382  nodea=nint(prop(index+1))
383  nodeb=nint(prop(index+2))
384 c iaxial=nint(prop(index+3))
385  offset=prop(index+4)
386  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
387  & co(1,nodea)-vold(1,nodea))**2)-offset
388 !
389  initial_radius=dsqrt((co(1,nodeb)-co(1,nodea))**2)-offset
390 !
391 c if(iaxial.ne.0) then
392 c A=pi*radius**2/iaxial
393 c else
394  a=pi*radius**2
395 c endif
396  d=2*radius
397  reynolds=dabs(xflow)*d/(dvi*a)
398  cd=1.d0
399 !
400  endif
401 !
402  if(cd.gt.1.d0) then
403  Write(*,*) '*WARNING:'
404  Write(*,*) 'In RESTRICTOR',nelem
405  write(*,*) 'Cd greater than 1'
406  write (*,*) 'Calculation will proceed using Cd=1'
407  cd=1.d0
408  endif
409 !
410  p2p1=p2/p1
411  km1=kappa-1.d0
412  kp1=kappa+1.d0
413  kdkm1=kappa/km1
414  tdkp1=2.d0/kp1
415  c2=tdkp1**kdkm1
416  aeff=a*cd
417  dt1=dsqrt(t1)
418 !
419  if(p2p1.gt.c2) then
420  c1=dsqrt(2.d0*kdkm1/r)*aeff
421  km1dk=1.d0/kdkm1
422  y=p2p1**km1dk
423  x=dsqrt(1.d0-y)
424  ca1=-c1*x/(kappa*p1*y)
425  cb1=c1*km1dk/(2.d0*p1)
426  ca2=-ca1*p2p1-xflow*dt1/(p1*p1)
427  cb2=-cb1*p2p1
428  f=xflow*dt1/p1-c1*p2p1**(1.d0/kappa)*x
429  if(cb2.le.-(alambda+ca2)*x) then
430  df(1)=-alambda
431  elseif(cb2.ge.(alambda-ca2)*x) then
432  df(1)=alambda
433  else
434  df(1)=ca2+cb2/x
435  endif
436  df(2)=xflow/(2.d0*p1*dt1)
437  df(3)=inv*dt1/p1
438  if(cb1.le.-(alambda+ca1)*x) then
439  df(4)=-alambda
440  elseif(cb1.ge.(alambda-ca1)*x) then
441  df(4)=alambda
442  else
443  df(4)=ca1+cb1/x
444  endif
445  else
446  c3=dsqrt(kappa/r)*(tdkp1)**(kp1/(2.d0*km1))*aeff
447  f=xflow*dt1/p1-c3
448  df(1)=-xflow*dt1/(p1)**2
449  df(2)=xflow/(2*p1*dt1)
450  df(3)=inv*dt1/p1
451  df(4)=0.d0
452  endif
453 !
454 ! output
455 !
456  elseif(iflag.eq.3) then
457 !
458  pi=4.d0*datan(1.d0)
459  p1=v(2,node1)
460  p2=v(2,node2)
461  if(p1.ge.p2) then
462  inv=1
463  xflow=v(1,nodem)*iaxial
464  t1=v(0,node1)-physcon(1)
465  t2=v(0,node2)-physcon(1)
466  else
467  inv=-1
468  p1=v(2,node2)
469  p2=v(2,node1)
470  xflow=-v(1,nodem)*iaxial
471  t1=v(0,node2)-physcon(1)
472  t2=v(0,node1)-physcon(1)
473  endif
474 !
475 ! calculation of the dynamic viscosity
476 !
477  if(dabs(dvi).lt.1e-30) then
478  write(*,*) '*ERROR in orifice: '
479  write(*,*) ' no dynamic viscosity defined'
480  write(*,*) ' dvi= ',dvi
481  call exit(201)
482 c kgas=0
483 c call dynamic_viscosity(kgas,T1,dvi)
484  endif
485 !
486  index=ielprop(nelem)
487  kappa=(cp/(cp-r))
488  a=prop(index+1)
489 !
490  if((lakon(nelem)(4:5).ne.'BT').and.
491  & (lakon(nelem)(4:5).ne.'PN').and.
492  & (lakon(nelem)(4:5).ne.'C1')) then
493  d=prop(index+2)
494  dl=prop(index+3)
495  u=prop(index+7)
496  nelemswirl=nint(prop(index+8))
497  if(nelemswirl.eq.0) then
498  uref=0.d0
499  else
500 ! swirl generating element
501 !
502 ! preswirl nozzle
503  if(lakon(nelemswirl)(2:5).eq.'ORPN') then
504  uref=prop(ielprop(nelemswirl)+5)
505 ! rotating orifices
506  else if((lakon(nelemswirl)(2:5).eq.'ORMM').or.
507  & (lakon(nelemswirl)(2:5).eq.'ORMA').or.
508  & (lakon(nelemswirl)(2:5).eq.'ORPM').or.
509  & (lakon(nelemswirl)(2:5).eq.'ORPA')) then
510  uref=prop(ielprop(nelemswirl)+7)
511 ! forced vortex
512  elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
513  uref=prop(ielprop(nelemswirl)+7)
514 !
515 ! free vortex
516  elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
517  uref=prop(ielprop(nelemswirl)+9)
518 ! Moehring
519  elseif(lakon(nelemswirl)(2:4).eq.'MRG') then
520  uref=prop(ielprop(nelemswirl)+10)
521 ! RCAVO
522  elseif((lakon(nelemswirl)(2:4).eq.'ROR').or.
523  & (lakon(nelemswirl)(2:4).eq.'ROA'))then
524  uref=prop(ielprop(nelemswirl)+6)
525 ! RCAVI
526  elseif(lakon(nelemswirl)(2:4).eq.'RCV') then
527  uref=prop(ielprop(nelemswirl)+5)
528  else
529  write(*,*) '*ERROR in orifice:'
530  write(*,*) ' element',nelemswirl
531  write(*,*) 'refered by element',nelem
532  write(*,*) 'is not a swirl generating element'
533  endif
534  endif
535 ! write(*,*) 'nelem',nelem, u, uref
536  u=u-uref
537  angle=prop(index+5)
538 !
539  endif
540 !
541 ! calculate the discharge coefficient using Bragg's Method
542 ! "Effect of Compressibility on the discharge coefficient
543 ! of orifices and convergent nozzles"
544 ! journal of mechanical Engineering
545 ! vol2 No 1 1960
546 !
547  if((lakon(nelem)(2:5).eq.'ORBG')) then
548 !
549  p2p1=p2/p1
550  d=dsqrt(a*4/pi)
551  reynolds=dabs(xflow)*d/(dvi*a)
552  cdcrit=prop(index+2)
553 !
554  itype=2
555  call cd_bragg(cdcrit,p2p1,cd,itype)
556 !
557  elseif(lakon(nelem)(2:5).eq.'ORMA') then
558 !
559 ! calculate the discharge coefficient using own table data and
560 ! using Dr.Albers method for rotating cavities
561 !
562  reynolds=dabs(xflow)*d/(dvi*a)
563 !
564  call cd_own_albers(p1,p2,dl,d,cd,u,t1,r,kappa)
565 !
566 ! chamfer correction
567 !
568  if(angle.gt.0.d0)then
569  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
570  cd=cd*cd_chamf
571  endif
572 !
573  elseif(lakon(nelem)(2:5).eq.'ORMM') then
574 !
575 ! calculate the discharge coefficient using McGreehan and Schotsch method
576 !
577  rad=prop(index+4)
578 !
579  reynolds=dabs(xflow)*d/(dvi*a)
580 !
581  call cd_ms_ms(p1,p2,t1,rad,d,dl,kappa,r,reynolds,u,vid,cd)
582 !
583  if(cd.ge.1) then
584  write(*,*) ''
585  write(*,*) '**WARNING**'
586  write(*,*) 'in RESTRICTOR ',nelem
587  write(*,*) 'Calculation using'
588  write(*,*) ' McGreehan and Schotsch method:'
589  write(*,*) ' Cd=',cd,'>1 !'
590  write(*,*) 'Calcultion will proceed will Cd=1'
591  write(*,*) 'l/d=',dl/d,'r/d=',rad/d,'u/vid=',u/vid
592  cd=1.d0
593  endif
594 !
595 ! chamfer correction
596 !
597  if(angle.gt.0.d0) then
598  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
599  cd=cd*cd_chamf
600  endif
601 !
602  elseif (lakon(nelem)(2:5).eq.'ORPA') then
603 !
604 ! calculate the discharge coefficient using Parker and Kercher method
605 ! and using Dr. Albers method for rotating cavities
606 !
607  rad=prop(index+4)
608 !
609  beta=prop(index+6)
610 !
611  reynolds=dabs(xflow)*d/(dvi*a)
612 !
613  call cd_pk_albers(rad,d,dl,reynolds,p2,p1,beta,kappa,
614  & cd,u,t1,r)
615 !
616 ! chamfer correction
617 !
618  if(angle.gt.0.d0) then
619  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
620  cd=cd*cd_chamf
621  endif
622 !
623  elseif(lakon(nelem)(2:5).eq.'ORPM') then
624 !
625 ! calculate the discharge coefficient using Parker and Kercher method
626 ! and using Mac Grehan and Schotsch method for rotating cavities
627 !
628  rad=prop(index+4)
629 !
630  beta=prop(index+6)
631  reynolds=dabs(xflow)*d/(dvi*a)
632 !
633  call cd_pk_ms(rad,d,dl,reynolds,p2,p1,beta,kappa,cd,
634  & u,t1,r)
635 !
636 ! chamfer correction
637 !
638  if(angle.gt.0.d0) then
639  call cd_chamfer(dl,d,p1,p2,angle,cd_chamf)
640  cd=cd*cd_chamf
641  endif
642 !
643  elseif(lakon(nelem)(2:5).eq.'ORC1') then
644 !
645  d=dsqrt(a*4/pi)
646  reynolds=dabs(xflow)*d/(dvi*a)
647  cd=1.d0
648 !
649  elseif(lakon(nelem)(2:5).eq.'ORBT') then
650 !
651 ! calculate the discharge coefficient of bleed tappings (OWN tables)
652 !
653  d=dsqrt(a*pi/4)
654  reynolds=dabs(xflow)*d/(dvi*a)
655  ps1pt1=prop(index+2)
656  curve=nint(prop(index+3))
657  number=nint(prop(index+4))
658  reynolds=dabs(xflow)*d/(dvi*a)
659  if(number.ne.0.d0)then
660  do i=1,number
661  x_tab(i)=prop(index+2*i+3)
662  y_tab(i)=prop(index+2*i+4)
663  enddo
664  endif
665 !
666  call cd_bleedtapping(p2,p1,ps1pt1,number,curve,x_tab,y_tab,
667  & cd)
668 !
669  elseif(lakon(nelem)(2:5).eq.'ORPN') then
670 !
671 ! calculate the discharge coefficient of preswirl nozzle (OWN tables)
672 !
673  d=dsqrt(4*a/pi)
674  reynolds=dabs(xflow)*d/(dvi*a)
675  curve=nint(prop(index+4))
676  number=nint(prop(index+6))
677 !
678  if(number.ne.0.d0)then
679  do i=1,number
680  x_tab2(i)=prop(index+2*i+5)
681  y_tab2(i)=prop(index+2*i+6)
682  enddo
683  endif
684 !
685  call cd_preswirlnozzle(p2,p1,number,curve,x_tab2,y_tab2,cd)
686 !
687  theta=prop(index+2)
688  k_phi=prop(index+3)
689 !
690  if(p2/p1.gt.(2/(kappa+1.d0))**(kappa/(kappa-1.d0))) then
691  c2u_new=k_phi*cd*sin(theta*pi/180.d0)*r*
692  & dsqrt(2.d0*kappa/(r*(kappa-1)))*
693  & dsqrt(t1*(1.d0-(p2/p1)**((kappa-1)/kappa)))
694 !
695  else
696  c2u_new=k_phi*cd*sin(theta*pi/180.d0)*r*
697  & dsqrt(2.d0*kappa/(r*(kappa-1)))*
698  & dsqrt(t1*(1.d0-2/(kappa+1)))
699  endif
700  prop(index+5)=c2u_new
701  endif
702 !
703  if(cd.gt.1.d0) then
704  Write(*,*) '*WARNING:'
705  Write(*,*) 'In RESTRICTOR',nelem
706  write(*,*) 'Cd greater than 1'
707  write(*,*) 'Calculation will proceed using Cd=1'
708  cd=1.d0
709  endif
710  xflow_oil=0
711 !
712  if(iflag.eq.3)then
713 !
714  xmach=dsqrt(((p1/p2)**((kappa-1.d0)/kappa)-1.d0)*2.d0/
715  & (kappa-1.d0))
716  write(1,*) ''
717  write(1,55) ' from node ',node1,
718  & ' to node ', node2,' : air massflow rate = '
719  &,inv*xflow,' ',
720  & ', oil massflow rate = ',xflow_oil,' '
721 !
722  if(inv.eq.1) then
723  write(1,56)' Inlet node ',node1,' : Tt1 = ',t1,
724  & ' , Ts1 = ',t1,' , Pt1 = ',p1, ' '
725 !
726  write(1,*)' Element ',nelem,lakon(nelem)
727  write(1,57)' dyn.visc = ',dvi,' , Re = '
728  & ,reynolds,', M = ',xmach
729  if(lakon(nelem)(2:5).eq.'ORC1') then
730  write(1,58)' CD = ',cd
731  else if((lakon(nelem)(2:5).eq.'ORMA').or.
732  & (lakon(nelem)(2:5).eq.'ORMM').or.
733  & (lakon(nelem)(2:5).eq.'ORPM').or.
734  & (lakon(nelem)(2:5).eq.'ORPA'))then
735  write(1,59)' CD = ',cd,' , C1u = ',u,
736  & ' , C2u = ', prop(index+7), ' '
737  endif
738 ! special for bleed tappings
739  if(lakon(nelem)(2:5).eq.'ORBT') then
740  write(1,63) ' P2/P1 = ',p2/p1,
741  &' , ps1pt1 = ', ps1pt1, ' , DAB = ',(1-p2/p1)/(1-ps1pt1),
742  &' , curve N° = ', curve,' , cd = ',cd
743 ! special for preswirlnozzles
744  elseif(lakon(nelem)(2:5).eq.'ORPN') then
745  write(1,62) ' cd = ', cd,
746  &' , C2u = ',c2u_new,' '
747 ! special for recievers
748  endif
749 !
750  write(1,56)' Outlet node ',node2,' : Tt2 = ',t2,
751  & ' , Ts2 = ',t2,' , Pt2 = ',p2,' '
752 !
753  else if(inv.eq.-1) then
754  write(1,56)' Inlet node ',node2,': Tt1 = ',t1,
755  & ' , Ts1 = ',t1,' , Pt1 = ',p1, ' '
756  &
757  write(1,*)' element R ',set(numf)(1:30)
758  write(1,57)' dyn.visc. = ',dvi,' , Re ='
759  & ,reynolds,', M = ',xmach
760  if(lakon(nelem)(2:5).eq.'ORC1') then
761  write(1,58)' CD = ',cd
762  else if((lakon(nelem)(2:5).eq.'ORMA').or.
763  & (lakon(nelem)(2:5).eq.'ORMM').or.
764  & (lakon(nelem)(2:5).eq.'ORPM').or.
765  & (lakon(nelem)(2:5).eq.'ORPA'))then
766  write(1,59)' CD = ',cd,' , C1u = ',u,
767  & ' , C2u = ', prop(index+7), ' '
768  endif
769 ! special for bleed tappings
770  if(lakon(nelem)(2:5).eq.'ORBT') then
771  write(1,63) ' P2/P1 = ',p2/p1,
772  &' , ps1pt1 = ', ps1pt1, ' , DAB = ',(1-p2/p1)/(1-ps1pt1),
773  &' , curve N° = ', curve,' , cd = ',cd
774 ! special for preswirlnozzles
775  elseif(lakon(nelem)(2:5).eq.'ORPN') then
776  write(1,*) ' cd = ', cd,' , C2u = '
777  &,c2u_new,' '
778  endif
779 !
780  write(1,56)' Outlet node ',node1,': Tt2 = ',t2,
781  & ' , Ts2 = ',t2,' , Pt2 = ',p2, ' '
782  endif
783  endif
784 !
785 !
786  endif
787 !
788  55 format(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a)
789  56 format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a)
790  57 format(1x,a,e11.4,a,e11.4,a,e11.4)
791  58 format(1x,a,e11.4)
792  59 format(1x,a,e11.4,a,e11.4,a,e11.4,a)
793  62 format(1x,a,e11.4,a,e11.4,a,e11.4,a)
794  63 format(1x,a,e11.4,a,e11.4,a,e11.4,a,i2,a,e11.4)
795 !
796  xflow=xflow/iaxial
797  df(3)=df(3)*iaxial
798 !
799  return
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
cd_pk_ms
subroutine cd_pk_ms(rad, d, xl, reynolds, p2, p1, beta, kappa, cd, u, T1, R)
Definition: cd_pk_ms.f:21
cd_ms_ms
subroutine cd_ms_ms(p1, p2, T1, rad, d, xl, kappa, r, reynolds, u, vid, cd)
Definition: cd_ms_ms.f:20
cd_preswirlnozzle
subroutine cd_preswirlnozzle(ps2, pt1, number, curve, x_tab, y_tab, cd)
Definition: cd_preswirlnozzle.f:25
cd_pk_albers
subroutine cd_pk_albers(rad, d, xl, reynolds, p2, p1, beta, kappa, cd, u, T1, R)
Definition: cd_pk_albers.f:23
c1
static double * c1
Definition: mafillvcompmain.c:30
cd_bragg
subroutine cd_bragg(cd, p2p1, cdbragg, itype)
Definition: cd_bragg.f:30
cd_own_albers
subroutine cd_own_albers(p1, p2, xl, d, cd, u, T1, R, kappa)
Definition: cd_own_albers.f:22
cd_bleedtapping
subroutine cd_bleedtapping(ps2, ps1, ps1pt1, nummer, curve, x_tab, y_tab, cd)
Definition: cd_bleedtapping.f:24
e11
static double * e11
Definition: radflowload.c:42
cd_chamfer
subroutine cd_chamfer(l, d, p_up, p_down, angle, cd)
Definition: cd_chamfer.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)