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! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine orifice(node1,node2,nodem,nelem,lakon,kon,ipkon, 20 & nactdog,identity,ielprop,prop,iflag,v,xflow,f, 21 & nodef,idirf,df,cp,R,physcon,dvi,numf,set,co,vold,mi,ttime, 22 & time,iaxial,iplausi) 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,iplausi 38! 39 real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,R,a,d,dl, 40 & p1,p2,T1,Aeff,C1,C2,C3,cd,cp,physcon(*),p2p1,km1,dvi, 41 & kp1,kdkm1,tdkp1,km1dk,x,y,ca1,cb1,ca2,cb2,dT1,alambda, 42 & rad,beta,reynolds,theta,k_phi,c2u_new,u,pi,xflow_oil, 43 & ps1pt1,uref,cd_chamf,angle,vid,cdcrit,T2,radius, 44 & initial_radius,co(3,*),vold(0:mi(2),*),offset,ttime,time, 45 & x_tab(100), y_tab(100),x_tab2(100),y_tab2(100),curve,xmach 46! 47! 48! 49 pi=4.d0*datan(1.d0) 50 if(iflag.eq.0) then 51 identity=.true. 52! 53 if(nactdog(2,node1).ne.0)then 54 identity=.false. 55 elseif(nactdog(2,node2).ne.0)then 56 identity=.false. 57 elseif(nactdog(1,nodem).ne.0)then 58 identity=.false. 59 endif 60! 61 elseif(iflag.eq.1)then 62 if(v(1,nodem).ne.0.d0) then 63 xflow=v(1,nodem) 64 return 65 endif 66! 67 index=ielprop(nelem) 68 kappa=(cp/(cp-R)) 69 a=prop(index+1) 70 d=prop(index+2) 71 dl=prop(index+3) 72! 73 if(lakon(nelem)(2:5).eq.'ORFL') then 74 nodea=nint(prop(index+1)) 75 nodeb=nint(prop(index+2)) 76 offset=prop(index+4) 77 radius=dsqrt((co(1,nodeb)+vold(1,nodeb)- 78 & co(1,nodea)-vold(1,nodea))**2)-offset 79 initial_radius=dsqrt((co(1,nodeb)-co(1,nodea))**2)-offset 80 A=pi*radius**2 81 d=2*radius 82 endif 83! 84 p1=v(2,node1) 85 p2=v(2,node2) 86 if(p1.ge.p2) then 87 inv=1 88 T1=v(0,node1)-physcon(1) 89 else 90 inv=-1 91 p1=v(2,node2) 92 p2=v(2,node1) 93 T1=v(0,node2)-physcon(1) 94 endif 95! 96 cd=1.d0 97! 98 p2p1=p2/p1 99 km1=kappa-1.d0 100 kp1=kappa+1.d0 101 kdkm1=kappa/km1 102 tdkp1=2.d0/kp1 103 C2=tdkp1**kdkm1 104 Aeff=A*cd 105 if(p2p1.gt.C2) then 106 xflow=inv*p1*Aeff*dsqrt(2.d0*kdkm1*p2p1**(2.d0/kappa) 107 & *(1.d0-p2p1**(1.d0/kdkm1))/r)/dsqrt(T1) 108 else 109 xflow=inv*p1*Aeff*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/ 110 & dsqrt(T1) 111 endif 112! 113 elseif(iflag.eq.2)then 114! 115 numf=4 116 alambda=10000.d0 117 index=ielprop(nelem) 118 kappa=(cp/(cp-R)) 119 a=prop(index+1) 120! 121 p1=v(2,node1) 122 p2=v(2,node2) 123 if(p1.ge.p2) then 124 inv=1 125 xflow=v(1,nodem)*iaxial 126 T1=v(0,node1)-physcon(1) 127 nodef(1)=node1 128 nodef(2)=node1 129 nodef(3)=nodem 130 nodef(4)=node2 131 else 132 inv=-1 133 p1=v(2,node2) 134 p2=v(2,node1) 135 xflow=-v(1,nodem)*iaxial 136 T1=v(0,node2)-physcon(1) 137 nodef(1)=node2 138 nodef(2)=node2 139 nodef(3)=nodem 140 nodef(4)=node1 141 endif 142! 143 idirf(1)=2 144 idirf(2)=0 145 idirf(3)=1 146 idirf(4)=2 147! 148! calculation of the dynamic viscosity 149! 150! 151 if(dabs(dvi).lt.1d-30) then 152 write(*,*) '*ERROR in orifice: ' 153 write(*,*) ' no dynamic viscosity defined' 154 write(*,*) ' dvi= ',dvi 155 call exit(201) 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)) 384c 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! 391c if(iaxial.ne.0) then 392c A=pi*radius**2/iaxial 393c else 394 A=pi*radius**2 395c 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.1d-30) then 478 write(*,*) '*ERROR in orifice: ' 479 write(*,*) ' no dynamic viscosity defined' 480 write(*,*) ' dvi= ',dvi 481 call exit(201) 482 endif 483! 484 index=ielprop(nelem) 485 kappa=(cp/(cp-R)) 486 a=prop(index+1) 487! 488 if((lakon(nelem)(4:5).ne.'BT').and. 489 & (lakon(nelem)(4:5).ne.'PN').and. 490 & (lakon(nelem)(4:5).ne.'C1')) then 491 d=prop(index+2) 492 dl=prop(index+3) 493 u=prop(index+7) 494 nelemswirl=nint(prop(index+8)) 495 if(nelemswirl.eq.0) then 496 uref=0.d0 497 else 498! swirl generating element 499! 500! preswirl nozzle 501 if(lakon(nelemswirl)(2:5).eq.'ORPN') then 502 uref=prop(ielprop(nelemswirl)+5) 503! rotating orifices 504 else if((lakon(nelemswirl)(2:5).eq.'ORMM').or. 505 & (lakon(nelemswirl)(2:5).eq.'ORMA').or. 506 & (lakon(nelemswirl)(2:5).eq.'ORPM').or. 507 & (lakon(nelemswirl)(2:5).eq.'ORPA')) then 508 uref=prop(ielprop(nelemswirl)+7) 509! forced vortex 510 elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then 511 uref=prop(ielprop(nelemswirl)+7) 512! 513! free vortex 514 elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then 515 uref=prop(ielprop(nelemswirl)+9) 516! Moehring 517 elseif(lakon(nelemswirl)(2:4).eq.'MRG') then 518 uref=prop(ielprop(nelemswirl)+10) 519! RCAVO 520 elseif((lakon(nelemswirl)(2:4).eq.'ROR').or. 521 & (lakon(nelemswirl)(2:4).eq.'ROA'))then 522 uref=prop(ielprop(nelemswirl)+6) 523! RCAVI 524 elseif(lakon(nelemswirl)(2:4).eq.'RCV') then 525 uref=prop(ielprop(nelemswirl)+5) 526 else 527 write(*,*) '*ERROR in orifice:' 528 write(*,*) ' element',nelemswirl 529 write(*,*) 'refered by element',nelem 530 write(*,*) 'is not a swirl generating element' 531 endif 532 endif 533! write(*,*) 'nelem',nelem, u, uref 534 u=u-uref 535 angle=prop(index+5) 536! 537 endif 538! 539! calculate the discharge coefficient using Bragg's Method 540! "Effect of Compressibility on the discharge coefficient 541! of orifices and convergent nozzles" 542! journal of mechanical Engineering 543! vol2 No 1 1960 544! 545 if((lakon(nelem)(2:5).eq.'ORBG')) then 546! 547 p2p1=p2/p1 548 d=dsqrt(a*4/Pi) 549 reynolds=dabs(xflow)*d/(dvi*a) 550 cdcrit=prop(index+2) 551! 552 itype=2 553 call cd_bragg(cdcrit,p2p1,cd,itype) 554! 555 elseif(lakon(nelem)(2:5).eq.'ORMA') then 556! 557! calculate the discharge coefficient using own table data and 558! using Dr.Albers method for rotating cavities 559! 560 reynolds=dabs(xflow)*d/(dvi*a) 561! 562 call cd_own_albers(p1,p2,dl,d,cd,u,T1,R,kappa) 563! 564! chamfer correction 565! 566 if(angle.gt.0.d0)then 567 call cd_chamfer(dl,d,p1,p2,angle,cd_chamf) 568 cd=cd*cd_chamf 569 endif 570! 571 elseif(lakon(nelem)(2:5).eq.'ORMM') then 572! 573! calculate the discharge coefficient using McGreehan and Schotsch method 574! 575 rad=prop(index+4) 576! 577 reynolds=dabs(xflow)*d/(dvi*a) 578! 579 call cd_ms_ms(p1,p2,T1,rad,d,dl,kappa,r,reynolds,u,vid,cd) 580! 581 if(cd.ge.1) then 582 write(*,*) '' 583 write(*,*) '**WARNING**' 584 write(*,*) 'in RESTRICTOR ',nelem 585 write(*,*) 'Calculation using' 586 write(*,*) ' McGreehan and Schotsch method:' 587 write(*,*) ' Cd=',Cd,'>1 !' 588 write(*,*) 'Calcultion will proceed will Cd=1' 589 write(*,*) 'l/d=',dl/d,'r/d=',rad/d,'u/vid=',u/vid 590 cd=1.d0 591 endif 592! 593! chamfer correction 594! 595 if(angle.gt.0.d0) then 596 call cd_chamfer(dl,d,p1,p2,angle,cd_chamf) 597 cd=cd*cd_chamf 598 endif 599! 600 elseif (lakon(nelem)(2:5).eq.'ORPA') then 601! 602! calculate the discharge coefficient using Parker and Kercher method 603! and using Dr. Albers method for rotating cavities 604! 605 rad=prop(index+4) 606! 607 beta=prop(index+6) 608! 609 reynolds=dabs(xflow)*d/(dvi*a) 610! 611 call cd_pk_albers(rad,d,dl,reynolds,p2,p1,beta,kappa, 612 & cd,u,T1,R) 613! 614! chamfer correction 615! 616 if(angle.gt.0.d0) then 617 call cd_chamfer(dl,d,p1,p2,angle,cd_chamf) 618 cd=cd*cd_chamf 619 endif 620! 621 elseif(lakon(nelem)(2:5).eq.'ORPM') then 622! 623! calculate the discharge coefficient using Parker and Kercher method 624! and using Mac Grehan and Schotsch method for rotating cavities 625! 626 rad=prop(index+4) 627! 628 beta=prop(index+6) 629 reynolds=dabs(xflow)*d/(dvi*a) 630! 631 call cd_pk_ms(rad,d,dl,reynolds,p2,p1,beta,kappa,cd, 632 & u,T1,R) 633! 634! chamfer correction 635! 636 if(angle.gt.0.d0) then 637 call cd_chamfer(dl,d,p1,p2,angle,cd_chamf) 638 cd=cd*cd_chamf 639 endif 640! 641 elseif(lakon(nelem)(2:5).eq.'ORC1') then 642! 643 d=dsqrt(a*4/Pi) 644 reynolds=dabs(xflow)*d/(dvi*a) 645 cd=1.d0 646! 647 elseif(lakon(nelem)(2:5).eq.'ORBT') then 648! 649! calculate the discharge coefficient of bleed tappings (OWN tables) 650! 651 d=dsqrt(A*Pi/4) 652 reynolds=dabs(xflow)*d/(dvi*a) 653 ps1pt1=prop(index+2) 654 curve=nint(prop(index+3)) 655 number=nint(prop(index+4)) 656 reynolds=dabs(xflow)*d/(dvi*a) 657 if(number.ne.0.d0)then 658 do i=1,number 659 x_tab(i)=prop(index+2*i+3) 660 y_tab(i)=prop(index+2*i+4) 661 enddo 662 endif 663! 664 call cd_bleedtapping(p2,p1,ps1pt1,number,curve,x_tab,y_tab, 665 & cd) 666! 667 elseif(lakon(nelem)(2:5).eq.'ORPN') then 668! 669! calculate the discharge coefficient of preswirl nozzle (OWN tables) 670! 671 d=dsqrt(4*A/pi) 672 reynolds=dabs(xflow)*d/(dvi*a) 673 curve=nint(prop(index+4)) 674 number=nint(prop(index+6)) 675! 676 if(number.ne.0.d0)then 677 do i=1,number 678 x_tab2(i)=prop(index+2*i+5) 679 y_tab2(i)=prop(index+2*i+6) 680 enddo 681 endif 682! 683 call cd_preswirlnozzle(p2,p1,number,curve,x_tab2,y_tab2,cd) 684! 685 theta=prop(index+2) 686 k_phi=prop(index+3) 687! 688 if(p2/p1.gt.(2/(kappa+1.d0))**(kappa/(kappa-1.d0))) then 689 c2u_new=k_phi*cd*sin(theta*Pi/180.d0)*r* 690 & dsqrt(2.d0*kappa/(r*(kappa-1)))* 691 & dsqrt(T1*(1.d0-(p2/p1)**((kappa-1)/kappa))) 692! 693 else 694 c2u_new=k_phi*cd*sin(theta*Pi/180.d0)*r* 695 & dsqrt(2.d0*kappa/(r*(kappa-1)))* 696 & dsqrt(T1*(1.d0-2/(kappa+1))) 697 endif 698 prop(index+5)=c2u_new 699 endif 700! 701 if(cd.gt.1.d0) then 702 Write(*,*) '*WARNING:' 703 Write(*,*) 'In RESTRICTOR',nelem 704 write(*,*) 'Cd greater than 1' 705 write(*,*) 'Calculation will proceed using Cd=1' 706 cd=1.d0 707 endif 708 xflow_oil=0 709! 710 if(iflag.eq.3)then 711! 712 xmach=dsqrt(((p1/p2)**((kappa-1.d0)/kappa)-1.d0)*2.d0/ 713 & (kappa-1.d0)) 714 write(1,*) '' 715 write(1,55) ' from node ',node1, 716 & ' to node ', node2,' : air massflow rate = ' 717 &,inv*xflow,' ', 718 & ', oil massflow rate = ',xflow_oil,' ' 719! 720 if(inv.eq.1) then 721 write(1,56)' Inlet node ',node1,' : Tt1 = ',T1, 722 & ' , Ts1 = ',T1,' , Pt1 = ',p1, ' ' 723! 724 write(1,*)' Element ',nelem,lakon(nelem) 725 write(1,57)' dyn.visc = ',dvi,' , Re = ' 726 & ,reynolds,', M = ',xmach 727 if(lakon(nelem)(2:5).eq.'ORC1') then 728 write(1,58)' CD = ',cd 729 else if((lakon(nelem)(2:5).eq.'ORMA').or. 730 & (lakon(nelem)(2:5).eq.'ORMM').or. 731 & (lakon(nelem)(2:5).eq.'ORPM').or. 732 & (lakon(nelem)(2:5).eq.'ORPA'))then 733 write(1,59)' CD = ',cd,' , C1u = ',u, 734 & ' , C2u = ', prop(index+7), ' ' 735 endif 736! special for bleed tappings 737 if(lakon(nelem)(2:5).eq.'ORBT') then 738 write(1,63) ' P2/P1 = ',p2/p1, 739 &' , ps1pt1 = ', ps1pt1, ' , DAB = ',(1-p2/p1)/(1-ps1pt1), 740 &' , curve N� = ', curve,' , cd = ',cd 741! special for preswirlnozzles 742 elseif(lakon(nelem)(2:5).eq.'ORPN') then 743 write(1,62) ' cd = ', cd, 744 &' , C2u = ',c2u_new,' ' 745! special for recievers 746 endif 747! 748 write(1,56)' Outlet node ',node2,' : Tt2 = ',T2, 749 & ' , Ts2 = ',T2,' , Pt2 = ',p2,' ' 750! 751 else if(inv.eq.-1) then 752 write(1,56)' Inlet node ',node2,': Tt1 = ',T1, 753 & ' , Ts1 = ',T1,' , Pt1 = ',p1, ' ' 754 & 755 write(1,*)' element R ',set(numf)(1:30) 756 write(1,57)' dyn.visc. = ',dvi,' , Re =' 757 & ,reynolds,', M = ',xmach 758 if(lakon(nelem)(2:5).eq.'ORC1') then 759 write(1,58)' CD = ',cd 760 else if((lakon(nelem)(2:5).eq.'ORMA').or. 761 & (lakon(nelem)(2:5).eq.'ORMM').or. 762 & (lakon(nelem)(2:5).eq.'ORPM').or. 763 & (lakon(nelem)(2:5).eq.'ORPA'))then 764 write(1,59)' CD = ',cd,' , C1u = ',u, 765 & ' , C2u = ', prop(index+7), ' ' 766 endif 767! special for bleed tappings 768 if(lakon(nelem)(2:5).eq.'ORBT') then 769 write(1,63) ' P2/P1 = ',p2/p1, 770 &' , ps1pt1 = ', ps1pt1, ' , DAB = ',(1-p2/p1)/(1-ps1pt1), 771 &' , curve N� = ', curve,' , cd = ',cd 772! special for preswirlnozzles 773 elseif(lakon(nelem)(2:5).eq.'ORPN') then 774 write(1,*) ' cd = ', cd,' , C2u = ' 775 &,c2u_new,' ' 776 endif 777! 778 write(1,56)' Outlet node ',node1,': Tt2 = ',T2, 779 & ' , Ts2 = ',T2,' , Pt2 = ',p2, ' ' 780 endif 781 endif 782! 783! 784 endif 785! 786 55 format(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a) 787 56 format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a) 788 57 format(1x,a,e11.4,a,e11.4,a,e11.4) 789 58 format(1x,a,e11.4) 790 59 format(1x,a,e11.4,a,e11.4,a,e11.4,a) 791 62 format(1x,a,e11.4,a,e11.4,a,e11.4,a) 792 63 format(1x,a,e11.4,a,e11.4,a,e11.4,a,i2,a,e11.4) 793! 794 xflow=xflow/iaxial 795 df(3)=df(3)*iaxial 796! 797 return 798 end 799