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