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 gaspipe_rot(node1,node2,nodem,nelem,lakon,kon, 20 & ipkon,nactdog,identity,ielprop,prop,iflag,v,xflow,f, 21 & nodef,idirf,df,cp,r,physcon,dvi,numf,set, 22 & shcon,nshcon,rhcon,nrhcon,ntmat_,co,vold,mi,ttime,time, 23 & iaxial,iplausi) 24! 25! rotating pipe with friction and variable cross section 26! 27 implicit none 28! 29 logical identity,crit 30 character*8 lakon(*) 31 character*81 set(*) 32! 33 integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,ielprop(*), 34 & idirf(*),index,iflag,nodef(*),inv,ipkon(*),kon(*),icase,ith, 35 & nshcon(*),nrhcon(*),ntmat_,mi(*),iaxial,iplausi 36! 37 real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,A,d,xl, 38 & T1,T2,Tt1,Tt2,pt1,pt2,cp,physcon(*),km1,dvi,kp1,kdkm1, 39 & reynolds,pi,lambda,kdkp1,rho,km1d2,ttime,time,pt2zpt1,ks, 40 & form_fact,pt2zpt1_c,Qred1_crit,Qred,phi,M1,M2,Qred1,co(3,*), 41 & shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),vold(0:mi(2),*), 42 & bb,cc,ee1,ee2,dfdM1,dfdM2,M1_c,Qred2,za,zb,zc,Z1,Z2,r1,r2, 43 & term,omega,om2,d1,d2,A1,A2,alpha,beta,coef,M2_c,Qred2_crit 44! 45 ith=0 46! 47 if(iflag.eq.0) then 48 identity=.true. 49! 50 if(nactdog(2,node1).ne.0)then 51 identity=.false. 52 elseif(nactdog(2,node2).ne.0)then 53 identity=.false. 54 elseif(nactdog(1,nodem).ne.0)then 55 identity=.false. 56 endif 57! 58 elseif(iflag.eq.1)then 59 if(v(1,nodem).ne.0.d0) then 60 xflow=v(1,nodem) 61 return 62 endif 63! 64 pi=4.d0*datan(1.d0) 65! 66 index=ielprop(nelem) 67 kappa=(cp/(cp-r)) 68! 69 A1=prop(index+1) 70 A2=prop(index+2) 71 xl=prop(index+3) 72 ks=prop(index+4) 73 form_fact=prop(index+5) 74 d1=prop(index+6) 75 if(form_fact.eq.1.d0) then 76 d1=2.d0*dsqrt(A1/pi) 77 endif 78 d2=prop(index+7) 79 if(form_fact.eq.1.d0) then 80 d2=2.d0*dsqrt(A2/pi) 81 endif 82 r1=prop(index+8) 83 r2=prop(index+9) 84 omega=prop(index+10) 85! 86 pt1=v(2,node1) 87 pt2=v(2,node2) 88! 89 Tt1=v(0,node1)-physcon(1) 90 Tt2=v(0,node2)-physcon(1) 91! 92 A=(A1+A2)/2.d0 93 d=(d1+d2)/2.d0 94 if(r2.ge.r1) then 95 om2=omega**2 96 else 97 om2=-omega**2 98 endif 99! 100! calculation of the dynamic viscosity 101! 102 if(dabs(dvi).lt.1d-30) then 103 write(*,*) '*ERROR in gaspipe_fanno: ' 104 write(*,*) ' no dynamic viscosity defined' 105 write(*,*) ' dvi= ',dvi 106 call exit(201) 107 endif 108! 109! assumed value for the reynolds number 110! 111 reynolds=3.e3 112! 113 call friction_coefficient(xl,d,ks,reynolds,form_fact, 114 & lambda) 115! 116! estimate of the flow using the incompressible relationships 117! for a gas pipe 118! 119! mean density 120! 121 rho=(pt1/(r*Tt1)+pt2/(r*Tt2))/2.d0 122 term=(rho*(pt1-pt2)+rho**2*om2*(r2**2-r1**2)/2.d0)* 123 & 2.d0*d/(lambda*xl) 124! 125 if(term.ge.0.d0) then 126 inv=1.d0 127 if(v(1,nodem).le.0.d0) then 128 xflow=A*dsqrt(term) 129 else 130 xflow=v(1,nodem)*iaxial 131 endif 132 else 133! 134! if the term underneath the square root is negative, 135! lambda must have a negative sign, which means that the flow 136! direction has to be reversed 137! 138 inv=-1.d0 139 lambda=-lambda 140 if(v(1,nodem).ge.0.d0) then 141 xflow=-A*dsqrt(-term) 142 else 143 xflow=v(1,nodem)*iaxial 144 endif 145 endif 146! 147 kp1=kappa+1.d0 148 km1=kappa-1.d0 149 km1d2=km1/2.d0 150! 151! check whether the flow does not exceed the critical one 152! 153 alpha=-4.d0*(d2-d1)/(xl*d)-(r1+r2)*om2*kp1/ 154 & (cp*(Tt1+Tt2)*km1) 155 if(alpha.eq.0.d0) then 156 write(*,*) '*ERROR in gaspipe_rot: looks like the' 157 write(*,*) ' cross section is constant and' 158 write(*,*) ' the rotational speed is zero;' 159 write(*,*) ' please use the gaspipe_fanno' 160 write(*,*) ' element instead' 161 call exit(201) 162 endif 163 beta=lambda*kappa/d 164! 165! for subsonic flow: 166! coef>0.d0 means that the Mach number increases from 1 to 2 167! (i.e. sonic conditions can only occur at 2) 168! coef<0.d0 means that the Mach number decreases from 1 to 2 169! (i.e. sonic conditions can only occur at 1) 170! 171 call ts_calc(xflow,Tt1,pt1,kappa,r,A1,T1,ith) 172 M1=dsqrt(((Tt1/T1)-1.d0)/km1d2) 173 call ts_calc(xflow,Tt2,pt2,kappa,r,A2,T2,ith) 174 M2=dsqrt(((Tt2/T2)-1.d0)/km1d2) 175 coef=alpha+beta*((M1+M2)/2.d0)**2 176! 177! icase tells where sonic conditions can occur, if any. 178! 179 if(coef.ge.0.d0) then 180 icase=2 181 else 182 icase=1 183 endif 184! 185 za=1.d0/alpha 186 zb=(alpha+beta)*beta/(alpha*(alpha*km1d2-beta)) 187 zc=kp1*km1d2/(2.d0*beta-alpha*km1) 188! 189 if(omega.eq.0.d0) then 190 call pt2zpt1_rot(pt2,pt1,kappa,r,xl,pt2zpt1_c,crit,icase, 191 & M1_c,M2_c,za,zb,zc,alpha,beta,Qred1_crit,Qred2_crit,A1,A2) 192 else 193 crit=.false. 194 endif 195! 196! check for critical flow only in the absence of rotational 197! speed 198! 199 if((icase.eq.1).and.(omega.eq.0.d0)) then 200c if(icase.eq.1) then 201! 202! decreasing Mach number from 1 to 2 203! 204 Qred2=dabs(xflow)*dsqrt(Tt2)/(A2*pt2) 205! 206! check whether flow is critical 207! assigning the physcical correct sign to xflow 208! 209 if(crit) then 210! 211! the flow is set to half the critical value or 212! one of the pressures is adapted (depending on 213! which variable is unknown) 214! 215c if(nactdog(1,nodem).ne.0) then 216 xflow=0.5d0*inv*Qred2_crit*pt2*A2/dsqrt(Tt2) 217c elseif(nactdog(2,node1).ne.0) then 218c if(beta.ge.0.d0) then 219c v(2,node1)=pt2/pt2zpt1_c*0.99d0 220c else 221c v(2,node1)=pt2/pt2zpt1_c*1.01d0 222c endif 223c elseif(nactdog(2,node2).ne.0) then 224c if(beta.ge.0.d0) then 225c v(2,node2)=pt2zpt1_c*pt1*1.01d0 226c else 227c v(2,node2)=pt2zpt1_c*pt1*0.99d0 228c endif 229c endif 230 elseif(Qred.gt.Qred2_crit) then 231! 232! the flow is set to half the critical value 233! 234 xflow=0.5d0*inv*Qred2_crit*pt2*A2/dsqrt(Tt2) 235 endif 236 else 237! 238! increasing Mach number from 1 to 2 239! 240 Qred=dabs(xflow)*dsqrt(Tt1)/(A1*pt1) 241 if(crit) then 242! 243! the flow is set to half the critical value or 244! one of the pressures is adapted (depending on 245! which variable is unknown) 246! 247c if(nactdog(1,nodem).ne.0) then 248 xflow=0.5d0*inv*Qred1_crit*pt1*A1/dsqrt(Tt1) 249c elseif(nactdog(2,node1).ne.0) then 250c if(beta.ge.0.d0) then 251c v(2,node1)=pt2/pt2zpt1_c*0.99d0 252c else 253c v(2,node1)=pt2/pt2zpt1_c*1.01d0 254c endif 255c elseif(nactdog(2,node2).ne.0) then 256c if(beta.ge.0.d0) then 257c v(2,node2)=pt2zpt1_c*pt1*1.01d0 258c else 259c v(2,node2)=pt2zpt1_c*pt1*0.99d0 260c endif 261c endif 262 elseif(Qred.gt.Qred1_crit) then 263! 264! the flow is set to half the critical value 265! 266 xflow=0.5d0*inv*Qred1_crit*pt1*A1/dsqrt(Tt1) 267 endif 268 endif 269! 270 elseif(iflag.eq.2)then 271! 272 numf=5 273! 274 pi=4.d0*datan(1.d0) 275! 276 index=ielprop(nelem) 277 kappa=(cp/(cp-r)) 278! 279 A1=prop(index+1) 280 A2=prop(index+2) 281 xl=prop(index+3) 282 ks=prop(index+4) 283 form_fact=prop(index+5) 284 d1=prop(index+6) 285 if(form_fact.eq.1.d0) then 286 d1=2.d0*dsqrt(A1/pi) 287 endif 288 d2=prop(index+7) 289 if(form_fact.eq.1.d0) then 290 d2=2.d0*dsqrt(A2/pi) 291 endif 292 r1=prop(index+8) 293 r2=prop(index+9) 294 omega=prop(index+10) 295! 296 pt1=v(2,node1) 297 pt2=v(2,node2) 298! 299 Tt1=v(0,node1)-physcon(1) 300 Tt2=v(0,node2)-physcon(1) 301! 302 xflow=v(1,nodem)*iaxial 303! 304 write(*,*) 'gaspipe_rot: ',pt1,pt2,xflow,Tt1,Tt2 305! 306 idirf(1)=2 307 idirf(2)=0 308 idirf(3)=1 309 idirf(4)=2 310 idirf(5)=0 311! 312 nodef(1)=node1 313 nodef(2)=node1 314 nodef(3)=nodem 315 nodef(4)=node2 316 nodef(5)=node2 317! 318 pt2zpt1=pt2/pt1 319! 320 A=(A1+A2)/2.d0 321 d=(d1+d2)/2.d0 322 if(r2.ge.r1) then 323 om2=omega**2 324 else 325 om2=-omega**2 326 endif 327! 328! calculation of the dynamic viscosity 329! 330 if(dabs(dvi).lt.1d-30) then 331 write(*,*) '*ERROR in gaspipe_fanno: ' 332 write(*,*) ' no dynamic viscosity defined' 333 write(*,*) ' dvi= ',dvi 334 call exit(201) 335 endif 336! 337 reynolds=dabs(xflow)*d/(dvi*A) 338! 339! calculation of the friction coefficient 340! 341 call friction_coefficient(xl,d,ks,reynolds,form_fact, 342 & lambda) 343! 344 if(xflow.lt.0.d0) then 345 lambda=-lambda 346 inv=-1 347 else 348 inv=1 349 endif 350! 351 kp1=kappa+1.d0 352 km1=kappa-1.d0 353 km1d2=km1/2.d0 354! 355 alpha=-4.d0*(d2-d1)/(xl*d)-(r1+r2)*om2*kp1/ 356 & (cp*(Tt1+Tt2)*km1) 357 if(alpha.eq.0.d0) then 358 write(*,*) '*ERROR in gaspipe_rot: looks like the' 359 write(*,*) ' cross section is constant and' 360 write(*,*) ' the rotational speed is zero;' 361 write(*,*) ' please use the gaspipe_fanno' 362 write(*,*) ' element instead' 363 call exit(201) 364 endif 365 beta=lambda*kappa/d 366! 367! for subsonic flow: 368! coef>0.d0 means that the Mach number increases from 1 to 2 369! (i.e. sonic conditions can only occur at 2) 370! coef<0.d0 means that the Mach number decreases from 1 to 2 371! (i.e. sonic conditions can only occur at 1) 372! 373 call ts_calc(xflow,Tt1,pt1,kappa,r,A1,T1,ith) 374 M1=dsqrt(((Tt1/T1)-1.d0)/km1d2) 375 call ts_calc(xflow,Tt2,pt2,kappa,r,A2,T2,ith) 376 M2=dsqrt(((Tt2/T2)-1.d0)/km1d2) 377 coef=alpha+beta*((M1+M2)/2.d0)**2 378 write(*,*) 'gaspipe_rot M1,M2',M1,M2 379! 380! icase tells where sonic conditions can occur, if any. 381! 382 if(coef.ge.0.d0) then 383 icase=2 384 else 385 icase=1 386 endif 387! 388 za=1.d0/alpha 389 zb=(alpha+beta)*beta/(alpha*(alpha*km1d2-beta)) 390 zc=kp1*km1d2/(2.d0*beta-alpha*km1) 391 write(*,*) 'gaspipe_rot alpha beta',alpha,beta 392 write(*,*) 'za zb zc',za,zb,zc 393! 394 if(omega.eq.0.d0) then 395 call pt2zpt1_rot(pt2,pt1,kappa,r,xl,pt2zpt1_c,crit,icase, 396 & M1_c,M2_c,za,zb,zc,alpha,beta,Qred1_crit,Qred2_crit,A1,A2) 397 else 398 crit=.false. 399 endif 400 write(*,*) 'gaspipe_rot crit ',crit 401! 402! check for critical flow only in the absence of rotational 403! speed 404! 405 if((icase.eq.1).and.(omega.eq.0.d0)) then 406! 407! decreasing Mach number from 1 to 2 408! 409 Qred2=dabs(xflow)*dsqrt(Tt2)/(A2*pt2) 410! 411! check whether flow is critical 412! assigning the physcical correct sign to xflow 413! 414 if(crit) then 415 write(*,*) '*WARNING in gaspipe_rot' 416 write(*,*) ' critical conditions detected' 417 write(*,*) ' in element ',nelem 418 xflow=inv*Qred2_crit*A2*pt2/dsqrt(Tt2) 419! 420! check whether flow has changed; if so, update v 421! for consistency 422! 423 if(dabs((xflow-iaxial*v(1,nodem))/xflow).gt.1.d-5) then 424 iplausi=0 425 if(nactdog(1,nodem).ne.0) v(1,nodem)=xflow/iaxial 426 endif 427 M2=dsqrt(((Tt2/T2)-1.d0)/km1d2) 428 M2=min(M2,0.999d0) 429 if((alpha+beta*M2*M2)/(alpha+beta).lt.0.d0) then 430 M2=M2_c 431 endif 432! 433! recalculate M1 434! 435 call ts_calc(xflow,Tt1,pt1,kappa,r,A1,T1,ith) 436 M1=dsqrt(((Tt1/T1)-1.d0)/km1d2) 437 else 438 if(Qred2.gt.Qred2_crit) then 439 xflow=inv*Qred2_crit*A2*pt2/dsqrt(Tt2) 440! 441! check whether flow has changed; if so, update v 442! for consistency 443! 444 if(dabs((xflow-iaxial*v(1,nodem))/xflow).gt.1.d-5) 445 & then 446 iplausi=0 447 if(nactdog(1,nodem).ne.0) v(1,nodem)=xflow/iaxial 448 endif 449! 450 M2=M2_c 451! 452! recalculate M1 453! 454 call ts_calc(xflow,Tt1,pt1,kappa,r,A1,T1,ith) 455 M1=dsqrt(((Tt1/T1)-1.d0)/km1d2) 456c else 457c call ts_calc(xflow,Tt2,pt2,kappa,r,A2,T2,ith) 458c M2=dsqrt(((Tt2/T2)-1.d0)/km1d2) 459 endif 460! 461c Tt1=Tt2 462c call ts_calc(xflow,Tt1,pt1,kappa,r,A1,T1,ith) 463c M1=dsqrt(((Tt1/T1)-1.d0)/km1d2) 464 endif 465 elseif((icase.eq.2).and.(omega.eq.0.d0)) then 466! 467! increasing Mach number from 1 to 2 468! 469 Qred1=dabs(xflow)*dsqrt(Tt1)/(A1*pt1) 470! 471! check whether flow is critical 472! assigning the physcical correct sign to xflow 473! 474 if(crit) then 475 write(*,*) '*WARNING in gaspipe_rot' 476 write(*,*) ' critical conditions detected' 477 write(*,*) ' in element ',nelem 478 xflow=inv*Qred1_crit*A1*pt1/dsqrt(Tt1) 479! 480! check whether flow has changed; if so, update v 481! for consistency 482! 483 if(dabs((xflow-iaxial*v(1,nodem))/xflow).gt.1.d-5) then 484 iplausi=0 485 if(nactdog(1,nodem).ne.0) v(1,nodem)=xflow/iaxial 486 endif 487! 488 M1=dsqrt(((Tt1/T1)-1.d0)/km1d2) 489 M1=min(M1,0.999d0) 490 if((alpha+beta)/(alpha+beta*M1*M1).lt.0.d0) then 491 M1=M1_c 492 endif 493! 494! recalculate M2 495! 496 call ts_calc(xflow,Tt2,pt2,kappa,r,A2,T2,ith) 497 M2=dsqrt(((Tt2/T2)-1.d0)/km1d2) 498 else 499 if(Qred1.gt.Qred1_crit) then 500 xflow=inv*Qred1_crit*A1*pt1/dsqrt(Tt1) 501! 502! check whether flow has changed; if so, update v 503! for consistency 504! 505 if(dabs((xflow-iaxial*v(1,nodem))/xflow).gt.1.d-5) 506 & then 507 iplausi=0 508 if(nactdog(1,nodem).ne.0) v(1,nodem)=xflow/iaxial 509 endif 510! 511 M1=M1_c 512! 513! recalculate M2 514! 515 call ts_calc(xflow,Tt2,pt2,kappa,r,A2,T2,ith) 516 M2=dsqrt(((Tt2/T2)-1.d0)/km1d2) 517c else 518c call ts_calc(xflow,Tt1,pt1,kappa,r,A1,T1,ith) 519c M1=dsqrt(((Tt1/T1)-1.d0)/km1d2) 520 endif 521! 522c Tt2=Tt1 523 call ts_calc(xflow,Tt2,pt2,kappa,r,A2,T2,ith) 524 M2=dsqrt(((Tt2/T2)-1.d0)/km1d2) 525 endif 526c elseif(icase.eq.1) then 527c else 528 endif 529! 530 bb=km1d2 531 cc=-kp1/(2.d0*km1) 532! 533! definition of the coefficients 534! 535 if(icase.eq.1) then 536! 537! decreasing Mach number from 1 to 2 538! 539 Z2=M2**2 540 ee2=M2*(1.d0+bb*Z2)/(1.d0+bb*Z2*(1.d0+2.d0*cc)) 541 dfdM2=2.d0*M2*(za/Z2+zb/(alpha+beta*Z2) 542 & +zc/(1.d0+km1d2*Z2)) 543! 544 if(.not.crit) then 545! 546! residual 547! 548 Z1=M1**2 549 ee1=M1*(1.d0+bb*Z1)/(1.d0+bb*Z1*(1.d0+2.d0*cc)) 550 dfdM1=-2.d0*M1*(za/Z1+zb/(alpha+beta*Z1) 551 & +zc/(1.d0+km1d2*Z1)) 552! 553 f=dlog((Z2/Z1)**za 554 & *((alpha+beta*Z2)/(alpha+beta*Z1))**(zb/beta) 555 & *((1.d0+km1d2*Z2)/(1.d0+km1d2*Z1))**(zc/km1d2)) 556 & -xl 557! 558! pressure node1 559! 560 df(1)=-dfdM1*ee1/pt1 561! 562! temperature node1 563! 564 df(2)=dfdM1*ee1/(2.d0*Tt1) 565! 566! mass flow 567! 568 df(3)=(dfdM1*ee1+dfdM2*ee2)/xflow 569! 570! pressure node2 571! 572 df(4)=-dfdM2*ee2/pt2 573! 574! temperature node2 575! 576 df(5)=dfdM2*ee2/(2.d0*Tt2) 577! 578 else 579! 580 f=dlog(Z2**za 581 & *((alpha+beta*Z2)/(alpha+beta))**(zb/beta) 582 & *((1.d0+km1d2*Z2)/(1.d0+km1d2))**(zc/km1d2)) 583 & -xl 584! 585! pressure node1 586! 587 df(1)=0.d0 588! 589! temperature node1 590! 591 df(2)=0.d0 592! 593! mass flow 594! 595 df(3)=(dfdM2*ee2)/xflow 596! 597! pressure node2 598! 599 df(4)=-dfdM2*ee2/pt2 600! 601! temperature node2 602! 603 df(5)=dfdM2*ee2/(2.d0*Tt2) 604! 605 endif 606 elseif(icase.eq.2) then 607! 608! increasing Mach number from 1 to 2 609! 610 Z1=M1**2 611 ee1=M1*(1.d0+bb*Z1)/(1.d0+bb*Z1*(1.d0+2.d0*cc)) 612 dfdM1=-2.d0*M1*(za/Z1+zb/(alpha+beta*Z1) 613 & +zc/(1.d0+km1d2*Z1)) 614! 615 if(.not.crit) then 616! 617! residual 618! 619 Z2=M2**2 620 ee2=M2*(1.d0+bb*Z2)/(1.d0+bb*Z2*(1.d0+2.d0*cc)) 621 dfdM2=2.d0*M2*(za/Z2+zb/(alpha+beta*Z2) 622 & +zc/(1.d0+km1d2*Z2)) 623! 624 f=dlog((Z2/Z1)**za 625 & *((alpha+beta*Z2)/(alpha+beta*Z1))**(zb/beta) 626 & *((1.d0+km1d2*Z2)/(1.d0+km1d2*Z1))**(zc/km1d2)) 627 & -xl 628! 629! pressure node1 630! 631 df(1)=-dfdM1*ee1/pt1 632! 633! temperature node1 634! 635 df(2)=dfdM1*ee1/(2.d0*Tt1) 636! 637! mass flow 638! 639 df(3)=(dfdM1*ee1+dfdM2*ee2)/xflow 640! 641! pressure node2 642! 643 df(4)=-dfdM2*ee2/pt2 644! 645! temperature node2 646! 647 df(5)=dfdM2*ee2/(2.d0*Tt2) 648! 649 else 650! 651 f=dlog((1.d0/Z1)**za 652 & *((alpha+beta)/(alpha+beta*Z1))**(zb/beta) 653 & *((1.d0+km1d2)/(1.d0+km1d2*Z1))**(zc/km1d2)) 654 & -xl 655! 656! pressure node1 657! 658 df(1)=-dfdM1*ee1/pt1 659! 660! temperature node1 661! 662 df(2)=dfdM1*ee1/(2.d0*Tt1) 663! 664! mass flow 665! 666 df(3)=dfdM1*ee1/xflow 667! 668! pressure node2 669! 670 df(4)=0.d0 671! 672! temperature node2 673! 674 df(5)=0.d0 675! 676 endif 677 endif 678! 679! output 680! 681 elseif(iflag.eq.3) then 682! 683 pi=4.d0*datan(1.d0) 684! 685 kappa=(cp/(cp-r)) 686 km1=kappa-1.d0 687 km1d2=km1/2.d0 688 kp1=kappa+1.d0 689 kdkm1=kappa/km1 690 kdkp1=kappa/kp1 691! 692 index=ielprop(nelem) 693 A1=prop(index+1) 694 A2=prop(index+2) 695 xl=prop(index+3) 696! 697 lambda=0.5d0 698! 699 ks=prop(index+4) 700 form_fact=prop(index+5) 701 d1=prop(index+6) 702 if(form_fact.eq.1.d0) then 703 d1=2.d0*dsqrt(A1/pi) 704 endif 705 d2=prop(index+7) 706 if(form_fact.eq.1.d0) then 707 d2=2.d0*dsqrt(A2/pi) 708 endif 709 r1=prop(index+8) 710 r2=prop(index+9) 711 omega=prop(index+10) 712! 713 pt1=v(2,node1) 714 pt2=v(2,node2) 715! 716 if(xflow.ge.0.d0) then 717 inv=1 718 xflow=v(1,nodem)*iaxial 719 Tt1=v(0,node1)-physcon(1) 720 Tt2=v(0,node2)-physcon(1) 721 else 722 inv=-1 723 pt1=v(2,node2) 724 pt2=v(2,node1) 725 xflow=v(1,nodem)*iaxial 726 Tt1=v(0,node2)-physcon(1) 727 Tt2=v(0,node1)-physcon(1) 728 endif 729 call ts_calc(xflow,Tt1,pt1,kappa,r,A1,T1,ith) 730c Tt2=Tt1 731 call ts_calc(xflow,Tt2,pt2,kappa,r,A2,T2,ith) 732! 733 pt2zpt1=pt2/pt1 734! 735 A=(A1+A2)/2.d0 736 d=(d1+d2)/2.d0 737! 738! calculation of the dynamic viscosity 739! 740 if(dabs(dvi).lt.1d-30) then 741 write(*,*) '*ERROR in gaspipe_fanno: ' 742 write(*,*) ' no dynamic viscosity defined' 743 write(*,*) ' dvi= ',dvi 744 call exit(201) 745 endif 746! 747 reynolds=dabs(xflow)*d/(dvi*A) 748! 749 if(reynolds.lt.1.d0) then 750 reynolds= 1.d0 751 endif 752! 753! definition of the friction coefficient for pure air 754! 755 phi=1.d0 756 call friction_coefficient(xl,d,ks,reynolds,form_fact, 757 & lambda) 758! 759! definition of the coefficients 760! 761 M1=dsqrt(((Tt1/T1)-1)/km1d2) 762 M2=dsqrt(((Tt2/T2)-1)/km1d2) 763! 764 write(1,*) '' 765 write(1,55) ' from node ',node1, 766 & ' to node ', node2,': air massflow rate = ',xflow 767! 768 if(inv.eq.1) then 769 write(1,53)' Inlet node ',node1,' : Tt1 = ',Tt1, 770 & ' , Ts1 = ',T1,' , Pt1 = ',pt1, 771 & ' , M1 = ',M1 772 write(1,*)' Element ',nelem,lakon(nelem) 773 write(1,57)' dvi = ',dvi,' , Re = ' 774 & ,reynolds 775 write(1,58)' PHI = ',phi,' , LAMBDA = ', 776 & lambda, 777 & ', LAMBDA*l/d = ',lambda*xl/d,' , ZETA_PHI = ', 778 & phi*lambda*xl/d 779 write(1,53)' Outlet node ',node2,' : Tt2 = ', 780 & Tt2, 781 & ' , Ts2 = ',T2,' , Pt2 = ',pt2, 782 & ' , M2 = ',M2 783! 784 else if(inv.eq.-1) then 785 write(1,53)' Inlet node ',node2,': Tt1= ',Tt1, 786 & ' , Ts1= ',T1,' , Pt1= ',pt1, 787 & ' , M1= ',M1 788 write(1,*)' Element ',nelem,lakon(nelem) 789 write(1,57)' dvi = ',dvi,' , Re = ' 790 & ,reynolds 791 write(1,58)' PHI = ',phi,' , LAMBDA = ', 792 & lambda, 793 & ', LAMBDA*l/d = ',lambda*xl/d,' , ZETA_PHI = ', 794 & phi*lambda*xl/d 795 write(1,53)' Outlet node ',node1,' : Tt2 = ', 796 & Tt2, 797 & ' , Ts2 = ',T2,' , Pt2 =',pt2, 798 & ' , M2 = ',M2 799 endif 800 endif 801! 802 55 format(1X,a,i6,a,i6,a,e11.4,a,e11.4) 803 53 format(1X,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4) 804 57 format(1X,a,e11.4,a,e11.4) 805 58 format(1X,a,e11.4,a,e11.4,a,e11.4,a,e11.4) 806! 807 xflow=xflow/iaxial 808 df(3)=df(3)*iaxial 809! 810 return 811 end 812