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
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),*)
55 if(nactdog(2,node1).ne.0)
then 57 elseif(nactdog(2,node2).ne.0)
then 59 elseif(nactdog(1,nodem).ne.0)
then 63 elseif(iflag.eq.1)
then 69 if(lakon(nelem)(2:5).ne.
'LABF')
then 85 elseif(lakon(nelem)(2:5).eq.
'LABF')
then 86 nodea=nint(prop(index+1))
87 nodeb=nint(prop(index+2))
100 s=dsqrt((co(1,nodeb)+vold(1,nodeb)-
101 & co(1,nodea)-vold(1,nodea))**2)
109 t1=v(0,node1)-physcon(1)
114 t1=v(0,node2)-physcon(1)
135 xflow=inv*p1*aeff*dsqrt(2.d0*kdkm1*p2p1**(2.d0/kappa)
136 & *(1.d0-p2p1**(1.d0/kdkm1))/r)/dsqrt(t1)
141 xflow=inv*p1*aeff*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
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))))
166 xflow=inv*p1*aeff/dsqrt(t1)*dsqrt(2.d0/r)*ppkrit
170 elseif(iflag.eq.2)
then 178 xflow=v(1,nodem)*iaxial
179 t1=v(0,node1)-physcon(1)
180 t2=v(0,node2)-physcon(1)
189 xflow=-v(1,nodem)*iaxial
190 t1=v(0,node2)-physcon(1)
191 t2=v(0,node1)-physcon(1)
205 if(lakon(nelem)(2:5).ne.
'LABF')
then 210 n=nint(prop(index+5))
221 elseif(lakon(nelem)(2:5).eq.
'LABF')
then 222 nodea=nint(prop(index+1))
223 nodeb=nint(prop(index+2))
227 n=nint(prop(index+6))
236 s=dsqrt((co(1,nodeb)+vold(1,nodeb)-
237 & co(1,nodea)-vold(1,nodea))**2)
251 cd_honeycomb=1+cd_honeycomb/100
257 if((rad.ne.0.d0).and.(n.ne.1d0))
then 263 if((n.ge.2).and.(hst.eq.0.d0))
then 266 carry_over=cst/dsqrt(cst-szt/(szt+0.02))
272 if(dabs(dvi).lt.1e-30)
then 273 write(*,*)
'*ERROR in labyrinth: ' 274 write(*,*)
' no dynamic viscosity defined' 275 write(*,*)
' dvi= ',dvi
281 reynolds=dabs(xflow)*2.d0*s/(dvi*a*cd_honeycomb/cd_radius)
310 call cd_bragg(cdu,p2p1,cdbragg,itype)
321 c1=dsqrt(2.d0*kdkm1/r)*aeff
325 ca1=-
c1*x/(kappa*p1*y)
326 cb1=
c1*km1dk/(2.d0*p1)
327 ca2=-ca1*p2p1-xflow*dt1/(p1*p1)
329 f=xflow*dt1/p1-
c1*p2p1**(1.d0/kappa)*x
330 if(cb2.le.-(alambda+ca2)*x)
then 332 elseif(cb2.ge.(alambda-ca2)*x)
then 337 df(2)=xflow/(2.d0*p1*dt1)
339 if(cb1.le.-(alambda+ca1)*x)
then 341 elseif(cb1.ge.(alambda-ca1)*x)
then 347 c3=dsqrt(kappa/r)*(tdkp1)**(kp1/(2.d0*km1))*aeff
349 df(1)=-xflow*dt1/(p1)**2
350 df(2)=xflow/(2*p1*dt1)
365 denom=r*(n-log(p2p1)/log(e))
369 if((hst.eq.0.d0).and.(n.ne.1))
then 371 aeff=aeff*cd_lab*cd_honeycomb*cd_radius
378 pt0zps1=(p1p2)**(1/prop(index+4))
387 cd=cd_1spike*cd_correction
390 aeff=aeff*cd_lab*cd_radius*cd_honeycomb
397 if(p2p1.gt.ppkrit)
then 399 f=xflow*dt1/p1-dsqrt(num/denom)*aeff
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)
406 df(4)=-aeff/2.d0*dsqrt(denom/num)*(-2.d0*(p2/p1**2.d0)
407 & /denom)+num/denom**2.d0*r/p2
412 c2=dsqrt(2/r)*aeff*ppkrit
415 df(1)=-xflow*dt1/(p1**2)
416 df(2)=xflow/(2.d0*p1*dt1)
424 elseif(iflag.eq.3)
then 431 xflow=v(1,nodem)*iaxial
432 t1=v(0,node1)-physcon(1)
433 t2=v(0,node2)-physcon(1)
442 xflow=-v(1,nodem)*iaxial
443 t1=v(0,node2)-physcon(1)
444 t2=v(0,node2)-physcon(1)
455 n=nint(prop(index+4))
469 e=2.718281828459045d0
475 aeff=aeff*(1.d0+cd_honeycomb/100.d0)
482 if((rad.ne.0.d0).and.(n.ne.1d0))
then 491 if((n.gt.1).and.(hst.eq.0.d0))
then 494 carry_over=cst/dsqrt(cst-szt/(szt+0.02))
500 if(dabs(dvi).lt.1e-30)
then 501 write(*,*)
'*ERROR in labyrinth: ' 502 write(*,*)
' no dynamic viscosity defined' 503 write(*,*)
' dvi= ',dvi
509 reynolds=dabs(xflow)*2.d0*s/(dvi*a)
537 call cd_bragg(cdu,p2p1,cdbragg,itype)
551 denom=r*(n-log(p2p1)/log(e))
555 if((hst.eq.0.d0).and.(n.ne.1))
then 557 aeff=aeff*cd_lab*cd_honeycomb*cd_radius
564 pt0zps1=(p1p2)**(1/prop(index+4))
573 cd=cd_1spike*cd_correction
576 aeff=aeff*cd_lab*cd_radius*cd_honeycomb
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)
592 write(1,56)
' Inlet node ',node1,
': Tt1=',t1,
593 &
', Ts1=',t1,
', Pt1=',p1
595 write(1,*)
' Element ',nelem,lakon(nelem)
596 write(1,57)
' dyn.visc.= ',dvi,
', Re= ' ,
598 &
', Cd_radius= ',cd_radius,
', Cd_honeycomb= ', 1+cd_honeycomb/100
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
606 elseif(hst.ne.0d0)
then 607 write(1,59)
' Cd_1_fin= ',
608 & cd_1spike,
', Cd= ',cd,
', pt0/ps1= ',pt0zps1,
613 write(1,60)
' Cd_Mcgreehan= ',cdu,
617 write(1,56)
' Outlet node ',node2,
': Tt2= ',t2,
618 &
', Ts2= ',t2,
', Pt2= ',p2
621 else if(inv.eq.-1)
then 622 write(1,56)
' Inlet node ',node2,
': Tt1= ',t1,
623 &
', Ts1= ',t1,
', Pt1= ',p1
625 write(1,*)
' element ',nelem,lakon(nelem)
626 write(1,57)
' dyn.visc.=',dvi,
', Re= ' 628 &
', Cd_radius= ',cd_radius,
', Cd_honeycomb= ',1+cd_honeycomb/100
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
636 elseif(hst.ne.0d0)
then 637 write(1,59)
' Cd_1_fin= ',
638 & cd_1spike,
', Cd= ',cd,
', pt0/ps1= ',pt0zps1,
643 write(1,60)
' Cd_Mcgreehan= ',
644 & cdu,
' Cd= ',cdbragg
646 write(1,56)
' Outlet node ',node1,
': Tt2= ',t2,
647 &
', Ts2= ',t2,
', Pt2= ',p2
651 56
FORMAT(1x,a,i6,a,
e11.4,a,
e11.4,a,
e11.4,a)
655 60
FORMAT(1x,a,
e11.4,a,
e11.4)
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