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

Go to the source code of this file.

Functions/Subroutines

subroutine gaspipe_fanno (node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, cp, r, physcon, dvi, numf, set, shcon, nshcon, rhcon, nrhcon, ntmat_, co, vold, mi, ttime, time, iaxial)
 

Function/Subroutine Documentation

◆ gaspipe_fanno()

subroutine gaspipe_fanno ( integer, intent(in)  node1,
integer, intent(in)  node2,
integer, intent(in)  nodem,
integer, intent(in)  nelem,
character*8, dimension(*), intent(in)  lakon,
integer, dimension(*), intent(in)  kon,
integer, dimension(*), intent(in)  ipkon,
integer, dimension(0:3,*), intent(in)  nactdog,
logical, intent(inout)  identity,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(in)  prop,
integer, intent(in)  iflag,
real*8, dimension(0:mi(2),*), intent(inout)  v,
real*8, intent(inout)  xflow,
real*8, intent(inout)  f,
integer, dimension(*), intent(inout)  nodef,
integer, dimension(*), intent(inout)  idirf,
real*8, dimension(*), intent(inout)  df,
real*8, intent(in)  cp,
real*8, intent(in)  r,
real*8, dimension(*), intent(in)  physcon,
real*8, intent(in)  dvi,
integer, intent(inout)  numf,
character*81, dimension(*), intent(in)  set,
real*8, dimension(0:3,ntmat_,*), intent(in)  shcon,
integer, dimension(*), intent(in)  nshcon,
real*8, dimension(0:1,ntmat_,*), intent(in)  rhcon,
integer, dimension(*), intent(in)  nrhcon,
integer, intent(in)  ntmat_,
real*8, dimension(3,*), intent(in)  co,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, dimension(*), intent(in)  mi,
real*8, intent(in)  ttime,
real*8, intent(in)  time,
integer, intent(in)  iaxial 
)
24 !
25 ! pipe with friction losses (Fanno Formulas) GAPF
26 !
27 ! author: Yannick Muller
28 !
29  implicit none
30 !
31  logical identity,crit,wrongdir
32  character*8 lakon(*)
33  character*81 set(*)
34 !
35  integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
36  & ielprop(*),nodef(*),idirf(*),index,iflag,
37  & inv,ipkon(*),kon(*),icase,k_oil
38  & ,nshcon(*),nrhcon(*),ntmat_,mi(*),nodea,nodeb,
39  & nodec,iaxial
40 !
41  real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,a,d,l,
42  & t1,t2,tt1,tt2,pt1,pt2,cp,physcon(*),p2p1,km1,dvi,
43  & kp1,kdkm1,reynolds,pi,lambda,lld,kdkp1,
44  & c2,tdkp1,ttime,time,pt2zpt1,ks,form_fact,xflow_oil,
45  & pt2zpt1_c,qred1_crit,qred,phi,m1,m2,qred1,co(3,*),
46  & shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),vold(0:mi(2),*),
47  & radius,bb,cc,ee1,ee2,dfdm1,dfdm2,m1_c,qred2
48 !
49  intent(in) node1,node2,nodem,nelem,lakon,kon,
50  & ipkon,nactdog,ielprop,prop,iflag,
51  & cp,r,physcon,dvi,set,
52  & shcon,nshcon,rhcon,nrhcon,ntmat_,co,vold,mi,ttime,time,
53  & iaxial
54 !
55  intent(inout) xflow,nodef,idirf,numf,f,df,v,identity
56 !
57  if(iflag.eq.0) then
58  identity=.true.
59 !
60  if(nactdog(2,node1).ne.0)then
61  identity=.false.
62  elseif(nactdog(2,node2).ne.0)then
63  identity=.false.
64  elseif(nactdog(1,nodem).ne.0)then
65  identity=.false.
66  endif
67 !
68  elseif(iflag.eq.1)then
69 !
70  pi=4.d0*datan(1.d0)
71 !
72  index=ielprop(nelem)
73  kappa=(cp/(cp-r))
74  a=prop(index+1)
75  d=prop(index+2)
76  l=prop(index+3)
77  ks=prop(index+4)
78  if(lakon(nelem)(2:6).eq.'GAPFA') then
79  icase=0
80  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
81  icase=1
82  endif
83  form_fact=prop(index+5)
84  xflow_oil=prop(index+6)
85  k_oil=nint(prop(index+7))
86 !
87  if(lakon(nelem)(7:8).eq.'FR') then
88 !
89 ! flexible radius
90 !
91  nodea=nint(prop(index+1))
92  nodeb=nint(prop(index+2))
93  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
94  & co(1,nodea)-vold(1,nodea))**2)
95 !
96  a=pi*radius**2
97  d=2*radius
98 !
99  elseif(lakon(nelem)(7:8).eq.'RL') then
100 !
101 ! flexible radius and length
102 !
103  nodea=nint(prop(index+1))
104  nodeb=nint(prop(index+2))
105  nodec=nint(prop(index+3))
106  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
107  & co(1,nodea)-vold(1,nodea))**2)
108  d=2*radius
109  a=pi*radius**2
110  l=dsqrt((co(2,nodec)+vold(2,nodec)-
111  & co(2,nodeb)-vold(2,nodeb))**2)
112  endif
113 !
114  pt1=v(2,node1)
115  pt2=v(2,node2)
116 c write(*,*) 'gaspipe_fanno a',nelem,pt1,pt2
117 !
118  if(pt1.ge.pt2) then
119  inv=1
120  tt1=v(0,node1)-physcon(1)
121  tt2=v(0,node2)-physcon(1)
122  else
123  inv=-1
124  pt1=v(2,node2)
125  pt2=v(2,node1)
126  tt1=v(0,node2)-physcon(1)
127  tt2=v(0,node1)-physcon(1)
128  endif
129 !
130  p2p1=pt2/pt1
131  km1=kappa-1.d0
132  kp1=kappa+1.d0
133  kdkm1=kappa/km1
134  tdkp1=2.d0/kp1
135  c2=tdkp1**kdkm1
136 !
137 ! estimate of the flow using the orifice relationships
138 ! the flow is needed for Reynolds, Reynolds is needed
139 ! for the friction coefficient
140 !
141  if(p2p1.gt.c2) then
142  xflow=inv*pt1*a*dsqrt(2.d0*kdkm1*p2p1**(2.d0/kappa)
143  & *(1.d0-p2p1**(1.d0/kdkm1))/r)/dsqrt(tt1)
144  else
145  xflow=inv*pt1*a*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
146  & dsqrt(tt1)
147  endif
148 !
149 ! calculation of the dynamic viscosity
150 !
151  if(dabs(dvi).lt.1e-30) then
152  write(*,*) '*ERROR in gaspipe_fanno: '
153  write(*,*) ' no dynamic viscosity defined'
154  write(*,*) ' dvi= ',dvi
155  call exit(201)
156  endif
157 !
158  reynolds=dabs(xflow)*d/(dvi*a)
159 !
160  call friction_coefficient(l,d,ks,reynolds,form_fact,
161  & lambda)
162 !
163 ! estimate of the flow using the incompressible relationships
164 ! for a gas pipe
165 !
166  xflow=inv*a*dsqrt(d/(lambda*l)*2*pt1/(r*tt1)*(pt1-pt2))
167 !
168  call pt2zpt1_crit(pt2,pt1,tt1,lambda,kappa,r,l,d,
169  & pt2zpt1_c,qred1_crit,crit,icase,m1_c)
170 !
171  qred=dabs(xflow)*dsqrt(tt1)/(a*pt1)
172 !
173  if(crit) then
174 !
175 ! the flow is set to half the critical value
176 !
177  xflow=0.5*inv*qred1_crit*pt1*a/dsqrt(tt1)
178  elseif(qred.gt.qred1_crit) then
179 !
180 ! the flow is set to half the critical value
181 !
182  xflow=0.5*inv*qred1_crit*pt1*a/dsqrt(tt1)
183  endif
184 c write(*,*) 'gaspipe_fanno a',nelem,xflow
185 !
186 ! isothermal case: correcting the temperatures
187 !
188  if(icase.eq.1) then
189  call ts_calc(xflow,tt1,pt1,kappa,r,a,t1,icase)
190  call ts_calc(xflow,tt2,pt2,kappa,r,a,t2,icase)
191  if(inv.eq.1) then
192  v(3,node1)=t1
193  v(3,node2)=t1
194  if(nactdog(0,node2).eq.1) then
195  v(0,node2)=t1*(tt2/t2)
196  endif
197  else
198  v(3,node2)=t1
199  v(3,node1)=t1
200  if(nactdog(0,node1).eq.1) then
201  v(0,node1)=t1*(tt2/t2)
202  endif
203  endif
204  endif
205 !
206  elseif(iflag.eq.2)then
207 !
208  numf=5
209 !
210  pi=4.d0*datan(1.d0)
211 !
212  kappa=(cp/(cp-r))
213  km1=kappa-1.d0
214  kp1=kappa+1.d0
215  kdkm1=kappa/km1
216  kdkp1=kappa/kp1
217 !
218  index=ielprop(nelem)
219  a=prop(index+1)
220  d=prop(index+2)
221 !
222  l=prop(index+3)
223  ks=prop(index+4)
224  if(lakon(nelem)(2:6).eq.'GAPFA') then
225  icase=0
226  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
227  icase=1
228  endif
229  form_fact=prop(index+5)
230  xflow_oil=prop(index+6)
231  k_oil=nint(prop(index+7))
232 !
233  if(lakon(nelem)(7:8).eq.'FR') then
234 !
235 ! flexible radius
236 !
237  nodea=nint(prop(index+1))
238  nodeb=nint(prop(index+2))
239  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
240  & co(1,nodea)-vold(1,nodea))**2)
241  a=pi*radius**2
242  d=2*radius
243 !
244  elseif(lakon(nelem)(7:8).eq.'RL') then
245 !
246 ! flexible radius and length
247 !
248  nodea=nint(prop(index+1))
249  nodeb=nint(prop(index+2))
250  nodec=nint(prop(index+3))
251  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
252  & co(1,nodea)-vold(1,nodea))**2)
253  d=2*radius
254  a=pi*radius**2
255  l=dsqrt((co(2,nodec)+vold(2,nodec)-
256  & co(2,nodeb)-vold(2,nodeb))**2)
257  endif
258 !
259  pt1=v(2,node1)
260  pt2=v(2,node2)
261  xflow=v(1,nodem)*iaxial
262 c write(*,*) 'gaspipe_fanno b',nelem,pt1,pt2,xflow
263 !
264 ! inv is the sign of the flow
265 ! xflow is replaced by its absolute value
266 ! wrongdir means that the flow goes from low
267 ! pressure to high pressure
268 !
269  inv=xflow/dabs(xflow)
270  xflow=dabs(xflow)
271  if((pt1-pt2)*inv.lt.0.d0) then
272  wrongdir=.true.
273  else
274  wrongdir=.false.
275  endif
276 !
277 ! the element is reoriented such that the mass flow
278 ! is directed from node 1 to node 2;
279 ! the pressure in node 1 may be more or less than
280 ! the pressure in node 2
281 !
282  if(pt1.gt.pt2) then
283 c inv=1
284 !
285  tt1=v(0,node1)-physcon(1)
286  call ts_calc(xflow,tt1,pt1,kappa,r,a,t1,icase)
287 !
288  nodef(1)=node1
289  nodef(2)=node1
290  nodef(3)=nodem
291  nodef(4)=node2
292  nodef(5)=node2
293  else
294 c inv=-1
295  pt1=v(2,node2)
296  pt2=v(2,node1)
297 c xflow=-v(1,nodem)*iaxial
298 !
299  tt1=v(0,node2)-physcon(1)
300  call ts_calc(xflow,tt1,pt1,kappa,r,a,t1,icase)
301 !
302  nodef(1)=node2
303  nodef(2)=node2
304  nodef(3)=nodem
305  nodef(4)=node1
306  nodef(5)=node1
307  endif
308 c write(*,*) 'gaspipe_fanno1 ',nelem,xflow,inv
309 !
310  idirf(1)=2
311  idirf(2)=0
312  idirf(3)=1
313  idirf(4)=2
314  idirf(5)=0
315 !
316  pt2zpt1=pt2/pt1
317 !
318 ! calculation of the dynamic viscosity
319 !
320  if(dabs(dvi).lt.1e-30) then
321  write(*,*) '*ERROR in gaspipe_fanno: '
322  write(*,*) ' no dynamic viscosity defined'
323  write(*,*) ' dvi= ',dvi
324  call exit(201)
325  endif
326 !
327  reynolds=xflow*d/(dvi*a)
328 cc reynolds=dabs(xflow)*d/(dvi*A)
329 c if(reynolds.lt.1) then
330 c reynolds=1.d0
331 c endif
332 !
333 ! if the flow goes from low to high pressure
334 ! large friction is applied (low Reynolds)
335 !
336 c if(wrongdir) then
337 c reynolds=1.d0
338 c else
339 c reynolds=xflow*d/(dvi*A)
340 c endif
341 !
342 ! calculation of the friction coefficient
343 !
344  if(xflow_oil.ne.0d0) then
345 !
346 ! two-phase-flow
347 !
348  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
349  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
350  & v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
351  & shcon,rhcon,ntmat_,mi)
352  lambda=lambda*phi
353  else
354 !
355 ! for pure air
356 !
357  phi=1.d0
358  call friction_coefficient(l,d,ks,reynolds,form_fact,
359  & lambda)
360  endif
361 !
362 ! calculating the critical conditions
363 !
364  call pt2zpt1_crit(pt2,pt1,tt1,lambda,kappa,r,l,d,
365  & pt2zpt1_c,qred1_crit,crit,icase,m1_c)
366 !
367  if(wrongdir) lambda=-lambda
368 !
369  qred1=xflow*dsqrt(tt1)/(a*pt1)
370 !
371 ! check whether flow is critical
372 ! assigning the physcical correct sign to xflow
373 !
374  if(crit) then
375  xflow=qred1_crit*a*pt1/dsqrt(tt1)
376 c xflow=inv*Qred1_crit*A*pt1/dsqrt(Tt1)
377  m1=dsqrt(2/km1*((tt1/t1)-1.d0))
378  if(icase.eq.0) then
379  m1=min(m1,0.999d0)
380  else
381  m1=min(m1,0.999d0/dsqrt(kappa))
382  endif
383  else
384  if(qred1.gt.qred1_crit) then
385 c xflow=inv*Qred1_crit*A*pt1/dsqrt(Tt1)
386  xflow=qred1_crit*a*pt1/dsqrt(tt1)
387  m1=m1_c
388  else
389 c xflow=inv*xflow
390  m1=dsqrt(2/km1*((tt1/t1)-1.d0))
391  endif
392 !
393 ! determining M2: Tt2 -> T2 -> Tt2/T2 -> M2
394 ! the actual value of Tt2 is not relevant
395 !
396  if(icase.eq.0) then
397  tt2=tt1
398  elseif(inv.eq.1) then
399  tt2=v(0,node2)-physcon(1)
400  else
401  tt2=v(0,node1)-physcon(1)
402  endif
403  call ts_calc(xflow,tt2,pt2,kappa,r,a,t2,icase)
404  m2=dsqrt(2/km1*((tt2/t2)-1.d0))
405 c Qred2=xflow*dsqrt(Tt1)/(A*pt2)
406 c write(*,*) 'gaspipe_fanno 1',M1,Qred1
407 c write(*,*) 'gaspipe_fanno 2',M1_c,Qred1_crit
408 c write(*,*) 'gaspipe_fanno 3',M2,Qred2
409 c write(*,*)
410  endif
411 c write(*,*) 'gaspipe_fanno2 ',xflow,inv,M1
412 c write(*,*) 'gaspipe_fanno3 ',M2,lambda,reynolds
413 !
414  bb=km1/2.d0
415  cc=-kp1/(2.d0*km1)
416  ee1=m1*(1.d0+bb*m1**2)/(1.d0+bb*m1**2*(1.d0+2.d0*cc))
417 !
418 ! definition of the coefficients
419 !
420  lld=lambda*l/d
421 !
422 ! adiabatic case
423 !
424  if(icase.eq.0) then
425 !
426  dfdm1=2.d0*(m1**2-1.d0)/(kappa*m1**3*(1.d0+bb*m1**2))
427 !
428  if(.not.crit) then
429 !
430 ! residual
431 !
432  ee2=m2*(1.d0+bb*m2**2)/(1.d0+bb*m2**2*(1.d0+2.d0*cc))
433  dfdm2=2.d0*(1.d0-m2**2)/(kappa*m2**3*(1.d0+bb*m2**2))
434 !
435  f=(1.d0/m1**2-1.d0/m2**2)/kappa+kp1/(2.d0*kappa)*
436  & dlog(((1.d0+bb*m2**2)*m1**2)/
437  & ((1.d0+bb*m1**2)*m2**2))-lld
438 !
439 ! pressure node1
440 !
441  df(1)=-dfdm1*ee1/pt1
442 !
443 ! temperature node1
444 !
445  df(2)=dfdm1*ee1/(2.d0*tt1)
446 !
447 ! mass flow
448 !
449  df(3)=(dfdm1*ee1+dfdm2*ee2)/(inv*xflow)
450 c df(3)=(dfdM1*ee1+dfdM2*ee2)/(xflow)
451 !
452 ! pressure node2
453 !
454  df(4)=-dfdm2*ee2/pt2
455 !
456 ! temperature node2
457 !
458  df(5)=dfdm2*ee2/(2.d0*tt2)
459 c write(*,*) 'gaspipe_fanno f',f
460 c write(*,*) 'gaspipe_fanno df(1)',df(1)
461 c write(*,*) 'gaspipe_fanno df(2)',df(2)
462 c write(*,*) 'gaspipe_fanno df(3)',df(3)
463 c write(*,*) 'gaspipe_fanno df(4)',df(4)
464 c write(*,*) 'gaspipe_fanno df(5)',df(5)
465 !
466  else
467 !
468  f=(1.d0/m1**2-1.d0)/kappa+kp1/(2.d0*kappa)*
469  & dlog(((1.d0+bb)*m1**2)/
470  & ((1.d0+bb*m1**2)))-lld
471 !
472 ! pressure node1
473 !
474  df(1)=-dfdm1*ee1/pt1
475 !
476 ! temperature node1
477 !
478  df(2)=dfdm1*ee1/(2.d0*tt1)
479 !
480 ! mass flow
481 !
482  df(3)=dfdm1*ee1/(inv*xflow)
483 c df(3)=dfdM1*ee1/(xflow)
484 !
485 ! pressure node2
486 !
487  df(4)=0.d0
488 !
489 ! temperature node2
490 !
491  df(5)=0.d0
492 !
493  endif
494  elseif(icase.eq.1) then
495 !
496 ! isothermal icase
497 !
498  dfdm1=2.d0*(kappa*m1**2-1.d0)/(kappa*m1**3)
499 c write(*,*) 'gaspipe dfdM1 ',dfdM1
500 !
501  if(.not.crit) then
502 !
503  ee2=m2*(1.d0+bb*m2**2)/(1.d0+bb*m2**2*(1.d0+2.d0*cc))
504  dfdm2=2.d0*(1.d0-kappa*m2**2)/(kappa*m2**3)
505 !
506 ! redidual
507 !
508  f=(1.d0/m1**2-1.d0/m2**2)/kappa+dlog((m1/m2)**2)-lld
509 !
510 ! pressure node1
511 !
512  df(1)=-dfdm1*ee1/pt1
513 !
514 ! temperature node1
515 !
516  df(2)=dfdm1*ee1/(2.d0*tt1)
517 !
518 ! mass flow
519 !
520  df(3)=(dfdm1*ee1+dfdm2*ee2)/(inv*xflow)
521 c df(3)=(dfdM1*ee1+dfdM2*ee2)/(xflow)
522 !
523 ! pressure node2
524 !
525  df(4)=-dfdm2*ee2/pt2
526 !
527 ! temperature node2
528 !
529  df(5)=dfdm2*ee2/(2.d0*tt2)
530 !
531  else
532 !
533 ! residual
534 !
535  f=(1.d0/m1**2-kappa)/kappa+dlog(kappa*m1**2)-lld
536 !
537 ! pressure node1
538 !
539  df(1)=-dfdm1*ee1/pt1
540 !
541 ! temperature node1
542 !
543  df(2)=dfdm1*ee1/(2.d0*tt1)
544 !
545 ! mass flow
546 !
547  df(3)=dfdm1*ee1/(inv*xflow)
548 c df(3)=dfdM1*ee1/(xflow)
549 c write(*,*) 'gaspipe f,df ',f,df(1),df(2),df(3)
550 !
551 ! pressure node2
552 !
553  df(4)=0.d0
554 !
555 ! temperature node2
556 !
557  df(5)=0.d0
558 !
559  endif
560  endif
561 !
562 ! output
563 !
564  elseif(iflag.eq.3) then
565 !
566  pi=4.d0*datan(1.d0)
567 !
568  kappa=(cp/(cp-r))
569  km1=kappa-1.d0
570  kp1=kappa+1.d0
571  kdkm1=kappa/km1
572  kdkp1=kappa/kp1
573 !
574  index=ielprop(nelem)
575  a=prop(index+1)
576  d=prop(index+2)
577  l=prop(index+3)
578 !
579  lambda=0.5
580 !
581  ks=prop(index+4)
582  if(lakon(nelem)(2:6).eq.'GAPFA') then
583  icase=0
584  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
585  icase=1
586  endif
587  form_fact=prop(index+5)
588  xflow_oil=prop(index+6)
589  k_oil=nint(prop(index+7))
590 !
591  pt1=v(2,node1)
592  pt2=v(2,node2)
593 !
594 c if(xflow.ge.0d0) then
595  if(pt1.gt.pt2) then
596  inv=1
597  xflow=v(1,nodem)*iaxial
598  tt1=v(0,node1)-physcon(1)
599  call ts_calc(xflow,tt1,pt1,kappa,r,a,t1,icase)
600  if(icase.eq.0) then
601  tt2=tt1
602  call ts_calc(xflow,tt2,pt2,kappa,r,a,t2,icase)
603  else
604  t2=t1
605  tt2=v(0,node2)-physcon(1)
606  call tt_calc(xflow,tt2,pt2,kappa,r,a,t2,icase)
607  endif
608 !
609  else
610  inv=-1
611  pt1=v(2,node2)
612  pt2=v(2,node1)
613  xflow=v(1,nodem)*iaxial
614  tt1=v(0,node2)-physcon(1)
615  call ts_calc(xflow,tt1,pt1,kappa,r,a,t1,icase)
616  if(icase.eq.0) then
617  tt2=tt1
618  call ts_calc(xflow,tt2,pt2,kappa,r,a,t2,icase)
619  else
620  t2=t1
621  tt2=v(0,node1)-physcon(1)
622  call tt_calc(xflow,tt2,pt2,kappa,r,a,t2,icase)
623  endif
624 !
625 c call ts_calc(xflow,Tt1,pt1,kappa,r,A,T1,icase)
626 c call ts_calc(xflow,Tt2,pt2,kappa,r,A,T2,icase)
627 !
628  endif
629 !
630  pt2zpt1=pt2/pt1
631 !
632 ! calculation of the dynamic viscosity
633 !
634  if(dabs(dvi).lt.1e-30) then
635  write(*,*) '*ERROR in gaspipe_fanno: '
636  write(*,*) ' no dynamic viscosity defined'
637  write(*,*) ' dvi= ',dvi
638  call exit(201)
639  endif
640 !
641  reynolds=dabs(xflow)*d/(dvi*a)
642 !
643  if(reynolds.lt.1.d0) then
644  reynolds= 1.d0
645  endif
646 !
647 ! definition of the friction coefficient for 2 phase flows and pure air
648 !
649  if(xflow_oil.ne.0d0) then
650  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
651  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
652  & v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
653  & shcon,rhcon,ntmat_,mi)
654  lambda=lambda*phi
655 !
656 ! for pure air
657 !
658  else
659  phi=1.d0
660  call friction_coefficient(l,d,ks,reynolds,form_fact,
661  & lambda)
662  endif
663 !
664  call pt2zpt1_crit(pt2,pt1,tt1,lambda,kappa,r,l,d,
665  & pt2zpt1_c,qred1_crit,crit,icase,m1_c)
666 !
667 ! definition of the coefficients
668 !
669  m1=dsqrt(2/km1*((tt1/t1)-1))
670  m2=dsqrt(2/km1*((tt2/t2)-1))
671 !
672  write(1,*) ''
673  write(1,55) ' from node ',node1,
674  & ' to node ', node2,': air massflow rate = ',xflow,
675  & ' , oil massflow rate = ',xflow_oil
676 !
677  if(inv.eq.1) then
678  write(1,53)' Inlet node ',node1,' : Tt1 = ',tt1,
679  & ' , Ts1 = ',t1,' , Pt1 = ',pt1,
680  & ' , M1 = ',m1
681  write(1,*)' Element ',nelem,lakon(nelem)
682  write(1,57)' dvi = ',dvi,' , Re = '
683  & ,reynolds
684  write(1,58)' PHI = ',phi,' , LAMBDA = ',
685  & lambda,
686  & ', LAMBDA*l/d = ',lambda*l/d,' , ZETA_PHI = ',
687  & phi*lambda*l/d
688  write(1,53)' Outlet node ',node2,' : Tt2 = ',
689  & tt2,
690  & ' , Ts2 = ',t2,' , Pt2 = ',pt2,
691  & ' , M2 = ',m2
692 !
693  else if(inv.eq.-1) then
694  write(1,53)' Inlet node ',node2,': Tt1= ',tt1,
695  & ' , Ts1= ',t1,' , Pt1= ',pt1,
696  & ' , M1= ',m1
697  write(1,*)' Element ',nelem,lakon(nelem)
698  write(1,57)' dvi = ',dvi,' , Re = '
699  & ,reynolds
700  write(1,58)' PHI = ',phi,' , LAMBDA = ',
701  & lambda,
702  & ', LAMBDA*l/d = ',lambda*l/d,' , ZETA_PHI = ',
703  & phi*lambda*l/d
704  write(1,53)' Outlet node ',node1,' : Tt2 = ',
705  & tt2,
706  & ' , Ts2 = ',t2,' , Pt2 =',pt2,
707  & ' , M2 = ',m2
708  endif
709  endif
710 !
711  55 format(1x,a,i6,a,i6,a,e11.4,a,e11.4)
712  53 format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
713  57 format(1x,a,e11.4,a,e11.4)
714  58 format(1x,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
715 !
716  xflow=xflow/iaxial
717  df(3)=df(3)*iaxial
718 !
719  return
pt2zpt1_crit
subroutine pt2zpt1_crit(pt2, pt1, Tt1, lambda, kappa, r, l, d, pt2zpt1_c, Qred1_crit, crit, icase, M1)
Definition: pt2zpt1_crit.f:35
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
min
#define min(a, b)
Definition: cascade.c:31
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
tt_calc
subroutine tt_calc(xflow, Tt, Pt, kappa, r, A, Ts, icase)
Definition: tt_calc.f:20
two_phase_flow
subroutine two_phase_flow(Tt1, pt1, T1, pt2, xflow_air, xflow_oil, nelem, lakon, kon, ipkon, ielprop, prop, v, dvi_air, cp, r, k_oil, phi, lambda, nshcon, nrhcon, shcon, rhcon, ntmat_, mi, iaxial)
Definition: two_phase_flow.f:23
e11
static double * e11
Definition: radflowload.c:42
Hosted by OpenAircraft.com, (Michigan UAV, LLC)