CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
labyrinth.f File Reference

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ labyrinth()

subroutine labyrinth ( integer  node1,
integer  node2,
integer  nodem,
integer  nelem,
character*8, dimension(*)  lakon,
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, dimension(3,*)  co,
real*8  dvi,
integer  numf,
real*8, dimension(0:mi(2),*)  vold,
character*81, dimension(*)  set,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
integer, dimension(*)  mi,
real*8  ttime,
real*8  time,
integer  iaxial 
)
22 !
23 ! labyrinth element
24 !
25 ! author: Yannick Muller
26 !
27  implicit none
28 !
29  logical identity
30  character*8 lakon(*)
31  character*81 set(*)
32 !
33  integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
34  & ielprop(*),nodef(*),idirf(*),index,iflag,mi(*),
35  & inv,kgas,n,iaxial,nodea,nodeb,ipkon(*),kon(*),i,itype
36 !
37  real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,a,d,
38  & p1,p2,t1,aeff,c1,c2,c3,cd,cp,physcon(*),p2p1,km1,dvi,
39  & kp1,kdkm1,tdkp1,km1dk,x,y,ca1,cb1,ca2,cb2,dt1,alambda,
40  & rad,reynolds,pi,ppkrit,co(3,*),ttime,time,
41  & carry_over,dlc,hst,e,szt,num,denom,t,s,b,h,cdu,
42  & cd_radius,cst,dh,cd_honeycomb,cd_lab,bdh,
43  & pt0zps1,cd_1spike,cdbragg,rzdh,
44  & cd_correction,p1p2,xflow_oil,t2,vold(0:mi(2),*)
45 !
46  itype=1
47  pi=4.d0*datan(1.d0)
48  e=2.718281828459045d0
49 !
50  index=ielprop(nelem)
51 !
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  kappa=(cp/(cp-r))
66 !
67 ! Usual Labyrinth
68 !
69  if(lakon(nelem)(2:5).ne.'LABF') then
70  t=prop(index+1)
71  s=prop(index+2)
72  d=prop(index+4)
73  n=nint(prop(index+5))
74  b=prop(index+6)
75  h=prop(index+7)
76  dlc=prop(index+8)
77  rad=prop(index+9)
78  x=prop(index+10)
79  hst=prop(index+11)
80 !
81  a=pi*d*s
82 !
83 ! "flexible" labyrinth for thermomechanical coupling
84 !
85  elseif(lakon(nelem)(2:5).eq.'LABF') then
86  nodea=nint(prop(index+1))
87  nodeb=nint(prop(index+2))
88  t=prop(index+4)
89  d=prop(index+5)
90  n=nint(prop(index+6))
91  b=prop(index+7)
92  h=prop(index+8)
93  dlc=prop(index+9)
94  rad=prop(index+10)
95  x=prop(index+11)
96  hst=prop(index+12)
97 
98 !
99 ! gap definition
100  s=dsqrt((co(1,nodeb)+vold(1,nodeb)-
101  & co(1,nodea)-vold(1,nodea))**2)
102  a=pi*d*s
103  endif
104 !
105  p1=v(2,node1)
106  p2=v(2,node2)
107  if(p1.ge.p2) then
108  inv=1
109  t1=v(0,node1)-physcon(1)
110  else
111  inv=-1
112  p1=v(2,node2)
113  p2=v(2,node1)
114  t1=v(0,node2)-physcon(1)
115  endif
116 !
117  cd=1.d0
118  aeff=a*cd
119  p2p1=p2/p1
120 !
121 !************************
122 ! one fin
123 !*************************
124  if(n.eq.1) then
125 !
126  km1=kappa-1.d0
127  kp1=kappa+1.d0
128  kdkm1=kappa/km1
129  tdkp1=2.d0/kp1
130  c2=tdkp1**kdkm1
131 !
132 ! subcritical
133 !
134  if(p2p1.gt.c2) then
135  xflow=inv*p1*aeff*dsqrt(2.d0*kdkm1*p2p1**(2.d0/kappa)
136  & *(1.d0-p2p1**(1.d0/kdkm1))/r)/dsqrt(t1)
137 !
138 ! critical
139 !
140  else
141  xflow=inv*p1*aeff*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
142  & dsqrt(t1)
143  endif
144  endif
145 !
146 !***********************
147 ! straight labyrinth and stepped labyrinth
148 ! method found in "Air system Correlations Part1 Labyrinth Seals"
149 ! H.Zimmermann and K.H. Wolff
150 ! ASME 98-GT-206
151 !**********************
152 !
153  if(n.ge.2) then
154 !
155  call lab_straight_ppkrit(n,ppkrit)
156 !
157 ! subcritical case
158 !
159  if(p2p1.gt.ppkrit) then
160  xflow=inv*p1*aeff/dsqrt(t1)*dsqrt((1.d0-p2p1**2.d0)
161  & /(r*(n-log(p2p1)/log(e))))
162 !
163 ! critical case
164 !
165  else
166  xflow=inv*p1*aeff/dsqrt(t1)*dsqrt(2.d0/r)*ppkrit
167  endif
168  endif
169 !
170  elseif(iflag.eq.2)then
171  numf=4
172  alambda=10000.d0
173 !
174  p1=v(2,node1)
175  p2=v(2,node2)
176  if(p1.ge.p2) then
177  inv=1
178  xflow=v(1,nodem)*iaxial
179  t1=v(0,node1)-physcon(1)
180  t2=v(0,node2)-physcon(1)
181  nodef(1)=node1
182  nodef(2)=node1
183  nodef(3)=nodem
184  nodef(4)=node2
185  else
186  inv=-1
187  p1=v(2,node2)
188  p2=v(2,node1)
189  xflow=-v(1,nodem)*iaxial
190  t1=v(0,node2)-physcon(1)
191  t2=v(0,node1)-physcon(1)
192  nodef(1)=node2
193  nodef(2)=node2
194  nodef(3)=nodem
195  nodef(4)=node1
196  endif
197 !
198  idirf(1)=2
199  idirf(2)=0
200  idirf(3)=1
201  idirf(4)=2
202 !
203 ! Usual labyrinth
204 !
205  if(lakon(nelem)(2:5).ne. 'LABF') then
206  kappa=(cp/(cp-r))
207  t=prop(index+1)
208  s=prop(index+2)
209  d=prop(index+4)
210  n=nint(prop(index+5))
211  b=prop(index+6)
212  h=prop(index+7)
213  dlc=prop(index+8)
214  rad=prop(index+9)
215  x=prop(index+10)
216  hst=prop(index+11)
217  a=pi*d*s
218 !
219 ! Flexible labyrinth for coupled calculations
220 !
221  elseif(lakon(nelem)(2:5).eq.'LABF') then
222  nodea=nint(prop(index+1))
223  nodeb=nint(prop(index+2))
224 c iaxial=nint(prop(index+3))
225  t=prop(index+4)
226  d=prop(index+5)
227  n=nint(prop(index+6))
228  b=prop(index+7)
229  h=prop(index+8)
230  dlc=prop(index+9)
231  rad=prop(index+10)
232  x=prop(index+11)
233  hst=prop(index+12)
234 !
235 ! gap definition
236  s=dsqrt((co(1,nodeb)+vold(1,nodeb)-
237  & co(1,nodea)-vold(1,nodea))**2)
238  a=pi*d*s
239  endif
240 !
241  p2p1=p2/p1
242  dt1=dsqrt(t1)
243 !
244  aeff=a
245 !
246 ! honeycomb stator correction
247 !
248  cd_honeycomb=1.d0
249  if(dlc.ne.0.d0)then
250  call cd_lab_honeycomb(s,dlc,cd_honeycomb)
251  cd_honeycomb=1+cd_honeycomb/100
252  endif
253 !
254 ! inlet radius correction
255 !
256  cd_radius=1.d0
257  if((rad.ne.0.d0).and.(n.ne.1d0)) then
258  call cd_lab_radius(rad,s,hst,cd_radius)
259  endif
260 !
261 ! carry over factor (only for straight throught labyrinth)
262 !
263  if((n.ge.2).and.(hst.eq.0.d0)) then
264  cst=n/(n-1.d0)
265  szt=s/t
266  carry_over=cst/dsqrt(cst-szt/(szt+0.02))
267  aeff=aeff*carry_over
268  endif
269 !
270 ! calculation of the dynamic viscosity
271 !
272  if(dabs(dvi).lt.1e-30) then
273  write(*,*) '*ERROR in labyrinth: '
274  write(*,*) ' no dynamic viscosity defined'
275  write(*,*) ' dvi= ',dvi
276  call exit(201)
277  endif
278 !
279 ! calculation of the number of reynolds for a gap
280 !
281  reynolds=dabs(xflow)*2.d0*s/(dvi*a*cd_honeycomb/cd_radius)
282 !
283 !**************************************
284 ! single fin labyrinth
285 ! the resolution procedure is the same as for the restrictor
286 !**************************************
287 !
288  if(n.eq.1)then
289 !
290 ! single fin labyrinth
291 !
292 ! incompressible basis cd , reynolds correction,and radius correction
293 !
294 ! "Flow Characteristics of long orifices with rotation and corner radiusing"
295 ! W.F. Mcgreehan and M.J. Schotsch
296 ! ASME 87-GT-162
297 !
298  dh=2*s
299  bdh=b/dh
300  rzdh=rad/dh
301 !
302  call cd_mcgreehan_schotsch(rzdh,bdh,reynolds,cdu)
303 !
304 ! compressibility correction factor
305 !
306 ! S.L.Bragg
307 ! "Effect of conpressibility on the discharge coefficient of orifices and convergent nozzles"
308 ! Journal of Mechanical engineering vol 2 No 1 1960
309 !
310  call cd_bragg(cdu,p2p1,cdbragg,itype)
311  cd=cdbragg
312  aeff=aeff*cd
313 !
314  km1=kappa-1.d0
315  kp1=kappa+1.d0
316  kdkm1=kappa/km1
317  tdkp1=2.d0/kp1
318  c2=tdkp1**kdkm1
319 !
320  if(p2p1.gt.c2) then
321  c1=dsqrt(2.d0*kdkm1/r)*aeff
322  km1dk=1.d0/kdkm1
323  y=p2p1**km1dk
324  x=dsqrt(1.d0-y)
325  ca1=-c1*x/(kappa*p1*y)
326  cb1=c1*km1dk/(2.d0*p1)
327  ca2=-ca1*p2p1-xflow*dt1/(p1*p1)
328  cb2=-cb1*p2p1
329  f=xflow*dt1/p1-c1*p2p1**(1.d0/kappa)*x
330  if(cb2.le.-(alambda+ca2)*x) then
331  df(1)=-alambda
332  elseif(cb2.ge.(alambda-ca2)*x) then
333  df(1)=alambda
334  else
335  df(1)=ca2+cb2/x
336  endif
337  df(2)=xflow/(2.d0*p1*dt1)
338  df(3)=inv*dt1/p1
339  if(cb1.le.-(alambda+ca1)*x) then
340  df(4)=-alambda
341  elseif(cb1.ge.(alambda-ca1)*x) then
342  df(4)=alambda
343  else
344  df(4)=ca1+cb1/x
345  endif
346  else
347  c3=dsqrt(kappa/r)*(tdkp1)**(kp1/(2.d0*km1))*aeff
348  f=xflow*dt1/p1-c3
349  df(1)=-xflow*dt1/(p1)**2
350  df(2)=xflow/(2*p1*dt1)
351  df(3)=inv*dt1/p1
352  df(4)=0.d0
353  endif
354  endif
355 !
356 !****************************************
357 ! straight labyrinth & stepped labyrinth
358 ! method found in "Air system Correlations Part1 Labyrinth Seals"
359 ! H.Zimmermann and K.H. Wolff
360 ! ASME 98-GT-206
361 !****************************************
362 !
363  if(n.ge.2) then
364  num=(1.d0-p2p1**2)
365  denom=r*(n-log(p2p1)/log(e))
366 !
367 ! straight labyrinth
368 !
369  if((hst.eq.0.d0).and.(n.ne.1)) then
370  call cd_lab_straight(n,p2p1,s,b,reynolds,cd_lab)
371  aeff=aeff*cd_lab*cd_honeycomb*cd_radius
372 !
373 ! Stepped Labyrinth
374 !
375  else
376 ! corrective term for the first spike
377  p1p2=p1/p2
378  pt0zps1=(p1p2)**(1/prop(index+4))
379  call cd_lab_1spike (pt0zps1,s,b,cd_1spike)
380 !
381 ! corrective term for cd_lab_1spike
382 !
383  call cd_lab_correction (p1p2,s,b,cd_correction)
384 !
385 ! calculation of the discharge coefficient of the stepped labyrinth
386 !
387  cd=cd_1spike*cd_correction
388  cd_lab=cd
389 !
390  aeff=aeff*cd_lab*cd_radius*cd_honeycomb
391  endif
392 !
393  call lab_straight_ppkrit(n,ppkrit)
394 !
395 ! subcritical case
396 !
397  if(p2p1.gt.ppkrit) then
398 !
399  f=xflow*dt1/p1-dsqrt(num/denom)*aeff
400 !
401  df(1)=xflow*dt1/p1**2.d0-aeff/2.d0
402  & *dsqrt(denom/num)*(2.d0*(p2**2.d0/p1**3.d0)/denom)
403  & +num/denom**2.d0*r/p1
404  df(2)=xflow/(2.d0*p1*dt1)
405  df(3)=inv*dt1/p1
406  df(4)=-aeff/2.d0*dsqrt(denom/num)*(-2.d0*(p2/p1**2.d0)
407  & /denom)+num/denom**2.d0*r/p2
408 !
409 ! critical case
410 !
411  else
412  c2=dsqrt(2/r)*aeff*ppkrit
413 !
414  f=xflow*dt1/p1-c2
415  df(1)=-xflow*dt1/(p1**2)
416  df(2)=xflow/(2.d0*p1*dt1)
417  df(3)=inv*dt1/p1
418  df(4)=0.d0
419  endif
420  endif
421 !
422 ! output
423 !
424  elseif(iflag.eq.3)then
425 !
426 
427  p1=v(2,node1)
428  p2=v(2,node2)
429  if(p1.ge.p2) then
430  inv=1
431  xflow=v(1,nodem)*iaxial
432  t1=v(0,node1)-physcon(1)
433  t2=v(0,node2)-physcon(1)
434  nodef(1)=node1
435  nodef(2)=node1
436  nodef(3)=nodem
437  nodef(4)=node2
438  else
439  inv=-1
440  p1=v(2,node2)
441  p2=v(2,node1)
442  xflow=-v(1,nodem)*iaxial
443  t1=v(0,node2)-physcon(1)
444  t2=v(0,node2)-physcon(1)
445  nodef(1)=node2
446  nodef(2)=node2
447  nodef(3)=nodem
448  nodef(4)=node1
449  endif
450 !
451  kappa=(cp/(cp-r))
452  t=prop(index+1)
453  s=prop(index+2)
454  d=prop(index+3)
455  n=nint(prop(index+4))
456  b=prop(index+5)
457  h=prop(index+6)
458  dlc=prop(index+7)
459  rad=prop(index+8)
460  x=prop(index+9)
461  hst=prop(index+10)
462 !
463  p2p1=p2/p1
464  dt1=dsqrt(t1)
465 !
466  pi=4.d0*datan(1.d0)
467  a=pi*d*s
468  aeff=a
469  e=2.718281828459045d0
470 !
471 ! honeycomb stator correction
472 !
473  if(dlc.ne.0.d0)then
474  call cd_lab_honeycomb(s,dlc,cd_honeycomb)
475  aeff=aeff*(1.d0+cd_honeycomb/100.d0)
476  else
477  cd_honeycomb=0
478  endif
479 !
480 ! inlet radius correction
481 !
482  if((rad.ne.0.d0).and.(n.ne.1d0)) then
483  call cd_lab_radius(rad,s,hst,cd_radius)
484  aeff=aeff*cd_radius
485  else
486  cd_radius=1
487  endif
488 !
489 ! carry over factor (only for straight throught labyrinth)
490 !
491  if((n.gt.1).and.(hst.eq.0.d0)) then
492  cst=n/(n-1.d0)
493  szt=s/t
494  carry_over=cst/dsqrt(cst-szt/(szt+0.02))
495  aeff=aeff*carry_over
496  endif
497 !
498 ! calculation of the dynamic viscosity
499 !
500  if(dabs(dvi).lt.1e-30) then
501  write(*,*) '*ERROR in labyrinth: '
502  write(*,*) ' no dynamic viscosity defined'
503  write(*,*) ' dvi= ',dvi
504  call exit(201)
505  endif
506 !
507 ! calculation of the number of reynolds for a gap
508 !
509  reynolds=dabs(xflow)*2.d0*s/(dvi*a)
510 !**************************************
511 ! single fin labyrinth
512 ! the resolution procedure is the same as for the restrictor
513 !**************************************
514 !
515  if(n.eq.1)then
516 !
517 ! single fin labyrinth
518 !
519 ! incompressible basis cd , reynolds correction,and radius correction
520 !
521 ! "Flow Characteristics of long orifices with rotation and corner radiusing"
522 ! W.F. Mcgreehan and M.J. Schotsch
523 ! ASME 87-GT-162
524 !
525  dh=2*s
526  bdh=b/dh
527  rzdh=rad/dh
528 !
529  call cd_mcgreehan_schotsch(rzdh,bdh,reynolds,cdu)
530 !
531 ! compressibility correction factor
532 !
533 ! S.L.Bragg
534 ! "Effect of conpressibility on the discharge coefficient of orifices and convergent nozzles"
535 ! Journal of Mechanical engineering vol 2 No 1 1960
536 !
537  call cd_bragg(cdu,p2p1,cdbragg,itype)
538  cd=cdbragg
539  aeff=aeff*cd
540  endif
541 !
542 !****************************************
543 ! straight labyrinth & stepped labyrinth
544 ! method found in "Air system Correlations Part1 Labyrinth Seals"
545 ! H.Zimmermann and K.H. Wolff
546 ! ASME 98-GT-206
547 !****************************************
548 !
549  if(n.ge.2) then
550  num=(1.d0-p2p1**2)
551  denom=r*(n-log(p2p1)/log(e))
552 !
553 ! straight labyrinth
554 !
555  if((hst.eq.0.d0).and.(n.ne.1)) then
556  call cd_lab_straight(n,p2p1,s,b,reynolds,cd_lab)
557  aeff=aeff*cd_lab*cd_honeycomb*cd_radius
558 !
559 ! Stepped Labyrinth
560 !
561  else
562 ! corrective term for the first spike
563  p1p2=p1/p2
564  pt0zps1=(p1p2)**(1/prop(index+4))
565  call cd_lab_1spike (pt0zps1,s,b,cd_1spike)
566 !
567 ! corrective term for cd_lab_1spike
568 !
569  call cd_lab_correction (p1p2,s,b,cd_correction)
570 !
571 ! calculation of the discharge coefficient of the stepped labyrinth
572 !
573  cd=cd_1spike*cd_correction
574  cd_lab=cd
575 !
576  aeff=aeff*cd_lab*cd_radius*cd_honeycomb
577  endif
578 !
579  call lab_straight_ppkrit(n,ppkrit)
580 !
581  endif
582 !
583  xflow_oil=0
584 !
585  write(1,*) ''
586  write(1,55) ' from node',node1,
587  &' to node', node2,': air massflow rate= ',xflow,
588  &', oil massflow rate= ',xflow_oil
589  55 FORMAT(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a)
590 
591  if(inv.eq.1) then
592  write(1,56)' Inlet node ',node1,': Tt1=',t1,
593  & ', Ts1=',t1,', Pt1=',p1
594 
595  write(1,*)' Element ',nelem,lakon(nelem)
596  write(1,57)' dyn.visc.= ',dvi,', Re= ' ,
597  & reynolds,
598  &', Cd_radius= ',cd_radius,', Cd_honeycomb= ', 1+cd_honeycomb/100
599 !
600 ! straight labyrinth
601  if((hst.eq.0.d0).and.(n.ne.1)) then
602  write(1,58)' COF= ',carry_over,
603  & ', Cd_lab= ',cd_lab,', Cd= ',carry_over*cd_lab
604 
605 ! stepped labyrinth
606  elseif(hst.ne.0d0) then
607  write(1,59)' Cd_1_fin= ',
608  & cd_1spike, ', Cd= ',cd,', pt0/ps1= ',pt0zps1,
609  & ', p0/pn= ',p1/p2
610 
611 ! single fin labyrinth
612  elseif(n.eq.1) then
613  write(1,60) ' Cd_Mcgreehan= ',cdu,
614  & ', Cd= ',cdbragg
615  endif
616 
617  write(1,56)' Outlet node ',node2,': Tt2= ',t2,
618  & ', Ts2= ',t2,', Pt2= ',p2
619 
620 !
621  else if(inv.eq.-1) then
622  write(1,56)' Inlet node ',node2,': Tt1= ',t1,
623  & ', Ts1= ',t1,', Pt1= ',p1
624 
625  write(1,*)' element ',nelem,lakon(nelem)
626  write(1,57)' dyn.visc.=',dvi,', Re= '
627  & ,reynolds,
628  & ', Cd_radius= ',cd_radius,', Cd_honeycomb= ',1+cd_honeycomb/100
629 !
630 ! straight labyrinth
631  if((hst.eq.0.d0).and.(n.ne.1)) then
632  write(1,58)' COF = ',carry_over,
633  & ', Cd_lab= ',cd_lab,', Cd= ',carry_over*cd_lab
634 !
635 ! stepped labyrinth
636  elseif(hst.ne.0d0) then
637  write(1,59)' Cd_1_fin= ',
638  & cd_1spike,', Cd= ',cd,', pt0/ps1= ',pt0zps1,
639  & ', p0/pn= ',p1/p2
640 
641 ! single fin labyrinth
642  elseif(n.eq.1) then
643  write(1,60) ' Cd_Mcgreehan= ',
644  & cdu,' Cd= ',cdbragg
645  endif
646  write(1,56)' Outlet node ',node1,': Tt2= ',t2,
647  & ', Ts2= ',t2,', Pt2= ',p2
648 
649  endif
650 !
651  56 FORMAT(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a)
652  57 FORMAT(1x,a,e11.5,a,e11.4,a,e11.4,a,e11.4)
653  58 FORMAT(1x,a,e11.4,a,e11.4,a,e11.4)
654  59 FORMAT(1x,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
655  60 FORMAT(1x,a,e11.4,a,e11.4)
656  endif
657 !
658  xflow=xflow/iaxial
659  df(3)=df(3)*iaxial
660 !
661  return
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
static double * c1
Definition: mafillvcompmain.c:30
subroutine cd_bragg(cd, p2p1, cdbragg, itype)
Definition: cd_bragg.f:30
subroutine cd_lab_radius(rad, s, hst, cd_radius)
Definition: cd_lab_radius.f:30
subroutine cd_lab_honeycomb(s, lc, cd_honeycomb)
Definition: cd_lab_honeycomb.f:30
static double * e11
Definition: radflowload.c:42
subroutine cd_lab_1spike(pt0zps1, s, b, cd_1spike)
Definition: cd_lab_1spike.f:34
subroutine cd_mcgreehan_schotsch(rzdh, bdh, reynolds, cdu)
Definition: cd_Mcgreehan_Schotsch.f:30
subroutine cd_lab_correction(p1p2, s, b, cd_correction)
Definition: cd_lab_correction.f:34
subroutine lab_straight_ppkrit(n, ppkrit)
Definition: lab_straight_ppkrit.f:32
subroutine cd_lab_straight(n, p2p1, s, b, reynolds, cd_lab)
Definition: cd_lab_straight.f:31
Hosted by OpenAircraft.com, (Michigan UAV, LLC)