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 preinitialnet(ieg,lakon,v,ipkon,kon,nflow,prop,ielprop, 20 & ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,mi,iponoel,inoel, 21 & itg,ntg) 22! 23! this routine only applies to compressible networks 24! 25! determination of initial values based on the boundary conditions 26! and the initial values given by the user by propagating these 27! through the network (using information on the mass flow direction 28! derived from unidirectional network elements or mass flow given 29! by the user (boundary conditions or initial conditions)) 30! 31! a mass flow with size 1.d-30 is used to propagate the sign 32! of the mass flow (in case only the sign and not the size is known) 33! 34! it is assumed that pressures, temperatures and mass flows cannot 35! be identically zero (a zero pressure does not make sense, 36! temperature has to be the absolute temperature, 37! a zero mass flow leads to convergence problems). 38! 39 implicit none 40! 41 character*8 lakon(*) 42! 43 integer mi(*),ieg(*),nflow,i,ielmat(mi(3),*),ntmat_,node1,node2, 44 & nelem,index,nshcon(*),ipkon(*),kon(*),nodem,imat,ielprop(*), 45 & nrhcon(*),neighbor,ichange,iponoel(*),inoel(2,*),indexe, 46 & itg(*),ntg,node,imin,imax,iel,nodemnei,ierror,nelemnei, 47 & nodenei,ibranch,numel,noderef,nelemref,ierr,j 48! 49 real*8 prop(*),shcon(0:3,ntmat_,*),xflow,v(0:mi(2),*),cp,r, 50 & dvi,rho,rhcon(0:1,ntmat_,*),kappa,cti,Ti,ri,ro,p1zp2,omega, 51 & p2zp1,xmin,xmax,fluxtot,ratio,pref,prefnew,r1,r2,xl,om2, 52 & Tt2mTt1,a1,c1,c2,d1,d2,disc 53! 54! the user should assign an initial pressure to any 55! - node which is connected to an inlet or an outlet 56! - node belonging to more than 2 network elements 57! this is checked in the next lines 58! 59 ierror=0 60! 61 do i=1,ntg 62 node=itg(i) 63 index=iponoel(node) 64! 65! node not connected to any element 66! 67 if(index.eq.0) cycle 68 nelem=inoel(1,index) 69! 70! midside node 71! 72 if(kon(ipkon(nelem)+2).eq.node) cycle 73! 74! initial pressure assigned 75! 76 if(v(2,node).gt.0.d0) cycle 77! 78! check whether any node in the element has a zero label 79! 80c if((kon(ipkon(nelem)+1).eq.0).or. 81c & (kon(ipkon(nelem)+3).eq.0)) then 82c write(*,*) '*ERROR in preinitialnet:' 83c write(*,*) ' node',node, 84c & ' is connected to an inlet or outlet, yet' 85c write(*,*) ' no initial pressure was assigned' 86c ierror=1 87c cycle 88c endif 89! 90! check whether node belongs to more than 1 element 91! 92 index=inoel(2,index) 93 if(index.eq.0) then 94 write(*,*) '*ERROR in preinitialnet:' 95 write(*,*) ' node',node, 96 & ' is connected to an inlet or outlet, yet' 97 write(*,*) ' no initial pressure was assigned' 98 ierror=1 99 cycle 100 endif 101! 102 call networkneighbor(nelem,node,nelemnei,nodenei,ibranch, 103 & iponoel,inoel,ipkon,kon) 104 if(nodenei.eq.0) then 105 write(*,*) '*ERROR in preinitialnet:' 106 write(*,*) ' node',node, 107 & ' is connected to an inlet or outlet, yet' 108 write(*,*) ' no initial pressure was assigned' 109 ierror=1 110 cycle 111 endif 112! 113 index=inoel(2,index) 114 if(index.ne.0) then 115 write(*,*) '*ERROR in preinitialnet:' 116 write(*,*) ' node',node, 117 & ' belongs to more than 2 network elements, yet' 118 write(*,*) ' no initial pressure was assigned' 119 ierror=1 120 endif 121 enddo 122 if(ierror.eq.1) call exit(201) 123! 124! for directional elements: small mass flow if none specified 125! 126 do i=1,nflow 127 nelem=ieg(i) 128 indexe=ipkon(nelem) 129! 130 nodem=kon(indexe+2) 131! 132 if(lakon(nelem)(2:3).eq.'OR') then 133! 134! orifice 135! 136 if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30 137 elseif(lakon(nelem)(2:4).eq.'LAB') then 138! 139! stepped labyrinth 140! 141 if((lakon(nelem)(5:6).eq.'SP').or. 142 & (lakon(nelem)(6:7).eq.'SP')) then 143 if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30 144 endif 145 elseif(lakon(nelem)(2:3).eq.'RE') then 146! 147! enlargement, contraction, wall orifice, entrance or exit 148! 149 if((lakon(nelem)(4:5).eq.'EL').or. 150 & (lakon(nelem)(4:5).eq.'CO').or. 151 & (lakon(nelem)(4:7).eq.'WAOR').or. 152 & (lakon(nelem)(4:5).eq.'EN').or. 153 & (lakon(nelem)(4:5).eq.'EX')) then 154 if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30 155 endif 156 elseif(lakon(nelem)(4:5).eq.'BR') then 157! 158! joint or split 159! 160 if((lakon(nelem)(6:6).eq.'J').or. 161 & (lakon(nelem)(6:6).eq.'S')) then 162 if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30 163 endif 164 elseif(lakon(nelem)(2:7).eq.'CROSPL') then 165! 166! cross split 167! 168 if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30 169 elseif(lakon(nelem)(2:3).eq.'MR') then 170! 171! Moehring 172! 173 if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30 174 endif 175 enddo 176! 177 do 178 ichange=0 179! 180! propagation of the pressure through the network 181! 182 do i=1,nflow 183 nelem=ieg(i) 184 indexe=ipkon(nelem) 185! 186 node1=kon(indexe+1) 187 if(node1.eq.0) cycle 188 nodem=kon(indexe+2) 189 node2=kon(indexe+3) 190 if(node2.eq.0) cycle 191! 192! exactly one total pressure value unknown in the element 193! (only this case is considered!) 194! 195 if(((v(2,node1).ne.0.d0).and.(v(2,node2).eq.0.d0)).or. 196 & ((v(2,node1).eq.0.d0).and.(v(2,node2).ne.0.d0))) then 197! 198 if(lakon(nelem)(2:3).eq.'VO') then 199! 200! vortex: pressure ratio can be determined 201! from geometry 202! 203 index=ielprop(nelem) 204! 205 if(prop(index+1).lt.prop(index+2)) then 206! 207! r2 < r1 208! 209 if(v(0,node2).ne.0.d0) then 210 ri=prop(index+1) 211 ro=prop(index+2) 212 Ti=v(0,node2) 213 if(lakon(nelem)(4:5).eq.'FO') then 214 omega=prop(index+5) 215 else 216 omega=prop(index+7) 217 endif 218! 219 imat=ielmat(1,nelem) 220 call materialdata_tg(imat,ntmat_,Ti,shcon, 221 & nshcon,cp,r,dvi,rhcon,nrhcon,rho) 222! 223 kappa=cp/(cp-r) 224 cti=omega*ri 225! 226 if(lakon(nelem)(4:5).eq.'FO') then 227! 228! forced vortex 229! 230 p1zp2=(1.d0+cti**2*((ro/ri)**2-1.d0)/ 231 & (2.d0*cp*Ti))**(kappa/(kappa-1.d0)) 232 else 233! 234! free vortex 235! 236 p1zp2=(1.d0+cti**2*(1.d0-(ri/ro)**2)/ 237 & (2.d0*cp*Ti))**(kappa/(kappa-1.d0)) 238 endif 239! 240 if(v(2,node1).eq.0.d0) then 241 v(2,node1)=v(2,node2)*p1zp2 242 else 243 v(2,node2)=v(2,node1)/p1zp2 244 endif 245 ichange=1 246 endif 247 else 248! 249! r1 <= r2 250! 251 if(v(0,node1).ne.0.d0) then 252 ri=prop(index+2) 253 ro=prop(index+1) 254 Ti=v(0,node1) 255 if(lakon(nelem)(4:5).eq.'FO') then 256 omega=prop(index+5) 257 else 258 omega=prop(index+7) 259 endif 260! 261 imat=ielmat(1,nelem) 262 call materialdata_tg(imat,ntmat_,Ti,shcon, 263 & nshcon,cp,r,dvi,rhcon,nrhcon,rho) 264! 265 kappa=cp/(cp-r) 266 cti=omega*ri 267! 268 if(lakon(nelem)(4:5).eq.'FO') then 269! 270! forced vortex 271! 272 p2zp1=(1.d0+cti**2*((ro/ri)**2-1.d0)/ 273 & (2.d0*cp*Ti))**(kappa/(kappa-1.d0)) 274 else 275! 276! free vortex 277! 278 p2zp1=(1.d0+cti**2*(1.d0-(ri/ro)**2)/ 279 & (2.d0*cp*Ti))**(kappa/(kappa-1.d0)) 280 endif 281! 282 if(v(2,node1).eq.0.d0) then 283 v(2,node1)=v(2,node2)/p2zp1 284 else 285 v(2,node2)=v(2,node1)*p2zp1 286 endif 287 ichange=1 288 endif 289 endif 290 elseif((lakon(nelem)(2:5).eq.'GAPR').and. 291 & (prop(ielprop(nelem)+10).gt.0.d0)) then 292! 293! truly rotating pipe 294! 295 index=ielprop(nelem) 296 r1=prop(index+8) 297 r2=prop(index+9) 298 if(((r2.lt.r1).and.(v(0,node2).gt.0.d0)).or. 299 & ((r2.ge.r1).and.(v(0,node1).gt.0.d0))) then 300! 301! estimating the pressure ratio across the pipe 302! 303 if(r2.lt.r1) then 304 Ti=v(0,node2) 305 om2=-prop(index+10)**2 306 else 307 Ti=v(0,node1) 308 om2=prop(index+10)**2 309 endif 310! 311 imat=ielmat(1,nelem) 312 call materialdata_tg(imat,ntmat_,Ti,shcon, 313 & nshcon,cp,r,dvi,rhcon,nrhcon,rho) 314 kappa=cp/(cp-r) 315! 316 xl=prop(index+3) 317! 318c p1zp2=dexp(kappa*om2*(r1+r2)*xl/ 319c & ((1.d0-kappa)*cp*Ti*2.d0)) 320! 321! improved formula 322! 323 p1zp2=(1.d0+om2*(r1+r2)*xl/(cp*v(0,node1)*2.d0)) 324 & **(kappa/(1.d0-kappa)) 325! 326c if(v(0,node1).gt.0.d0) then 327c a1=xl*r1/(r2-r1) 328c disc=a1**2-2.d0*xl*cp*v(0,node1)/ 329c & (om2*(r2-r1)) 330c if(disc.lt.0.d0) then 331c write(*,*) '*ERROR in preinitialnet:' 332c write(*,*) ' negative discriminant' 333c stop 334c endif 335c d1=-a1+dsqrt(disc) 336c d2=-a1-dsqrt(disc) 337c c1=(d1+a1)/(d1-d2) 338c c2=(d2+a1)/(d2-d1) 339c p2zp1=((d1-xl)/d1)**(2.d0*kappa*c1/(kappa-1))* 340c & ((d2-xl)/d2)**(2.d0*kappa*c2/(kappa-1)) 341cc p1zp2=1.d0/p2zp1 342c write(*,*) 'preinitialnet p1zp2 ', 343c & p1zp2,1.d0/p2zp1 344c else 345c a1=xl*r2/(r1-r2) 346c disc=a1**2-2.d0*xl*cp*v(0,node2)/ 347c & (om2*(r2-r1)) 348c if(disc.lt.0.d0) then 349c write(*,*) '*ERROR in preinitialnet:' 350c write(*,*) ' negative discriminant' 351c stop 352c endif 353c d1=-a1+dsqrt(disc) 354c d2=-a1-dsqrt(disc) 355c c1=(d1+a1)/(d1-d2) 356c c2=(d2+a1)/(d2-d1) 357c p2zp1=1.d0/( 358c & ((d1-xl)/d1)**(2.d0*kappa*c1/(kappa-1))* 359c & ((d2-xl)/d2)**(2.d0*kappa*c2/(kappa-1))) 360cc p1zp2=1.d0/p2zp1 361c write(*,*) 'preinitialnet p1zp2 ', 362c & p1zp2,1.d0/p2zp1 363c endif 364! 365 if(v(2,node1).eq.0.d0) then 366 v(2,node1)=v(2,node2)*p1zp2 367 else 368 v(2,node2)=v(2,node1)/p1zp2 369 endif 370 ichange=1 371 endif 372 elseif(v(1,nodem).ne.0.d0) then 373! 374! mass flow is given (either size or just the 375! direction): small slope 376! 377 ierror=0 378 if(v(1,nodem).gt.0.d0) then 379 if(v(2,node1).eq.0.d0) then 380 v(2,node1)=v(2,node2)*1.01d0 381 call networkneighbor(nelem,node1,nelemnei, 382 & nodenei,ibranch,iponoel,inoel,ipkon,kon) 383 if(v(2,nodenei).le.v(2,node1)) ierror=1 384 else 385 v(2,node2)=v(2,node1)*0.99d0 386 call networkneighbor(nelem,node2,nelemnei, 387 & nodenei,ibranch,iponoel,inoel,ipkon,kon) 388 if(v(2,nodenei).ge.v(2,node2)) ierror=2 389 endif 390 else 391 if(v(2,node1).eq.0.d0) then 392 v(2,node1)=v(2,node2)*0.99d0 393 call networkneighbor(nelem,node1,nelemnei, 394 & nodenei,ibranch,iponoel,inoel,ipkon,kon) 395 if(v(2,nodenei).ge.v(2,node1)) ierror=1 396 else 397 v(2,node2)=v(2,node1)*1.01d0 398 call networkneighbor(nelem,node2,nelemnei, 399 & nodenei,ibranch,iponoel,inoel,ipkon,kon) 400 if(v(2,nodenei).le.v(2,node2)) ierror=2 401 endif 402 endif 403! 404! if a discontinuity is detected in the pressure 405! the complete branch (i.e. the elements between 406! branch points and/or inlets and/or outlets) is 407! reanalyzed 408! 409 if(ierror.ne.0) then 410! 411! looking in the direction of node 1 until a 412! branch point/inlet/outlet 413! 414 node=node1 415 do 416 call networkneighbor(nelem,node,nelemnei, 417 & nodenei,ibranch,iponoel,inoel,ipkon,kon) 418 if((ibranch.eq.1).or.(nodenei.eq.0)) exit 419 node=nodenei 420 nelem=nelemnei 421 enddo 422! 423! noderef is branch point/inlet/outlet 424! the adjacent element into the branch is nelemref 425! 426 noderef=node 427 nelemref=nelem 428! 429 indexe=ipkon(nelemref) 430 if(kon(indexe+1).eq.noderef) then 431 node=kon(indexe+3) 432 else 433 node=kon(indexe+1) 434 endif 435! 436! looking in the other direction until 437! branch point/inlet/outlet 438! counting the non-vortex elements 439! storing the pressure ratio over the vortices 440! 441 ratio=1.d0 442 numel=0 443! 444 if(lakon(nelem)(2:3).eq.'VO') then 445 ratio=ratio*v(2,node)/v(2,noderef) 446 else 447 numel=numel+1 448 endif 449! 450 ierr=0 451 do 452 call networkneighbor(nelem,node,nelemnei, 453 & nodenei,ibranch,iponoel,inoel,ipkon,kon) 454 if((ibranch.eq.1).or.(nodenei.eq.0)) exit 455 if(lakon(nelemnei)(2:3).eq.'VO') then 456 if((v(2,node).eq.0.d0).or. 457 & (v(2,nodenei).eq.0.d0)) then 458 ierr=1 459 exit 460 endif 461 ratio=ratio*v(2,nodenei)/v(2,node) 462 else 463 numel=numel+1 464 endif 465 node=nodenei 466 nelem=nelemnei 467 enddo 468! 469 if(ierr.eq.0) then 470! 471! determining the required pressure ratio over the 472! non-vortex elements from the pressure at the ends of 473! the branch and the pressure ratio over the vortices 474! 475 ratio=(v(2,node)/(v(2,noderef)*ratio)) 476 & **(1.d0/numel) 477! 478! going through the branch again; determining the 479! initial values 480! 481 indexe=ipkon(nelemref) 482 if(kon(indexe+1).eq.noderef) then 483 node=kon(indexe+3) 484 else 485 node=kon(indexe+1) 486 endif 487! 488 nelem=nelemref 489 pref=v(2,node) 490 if(lakon(nelem)(2:3).ne.'VO') then 491 v(2,node)=v(2,noderef)*ratio 492 endif 493! 494 do 495 call networkneighbor(nelem,node,nelemnei, 496 & nodenei,ibranch,iponoel,inoel,ipkon,kon) 497 if((ibranch.eq.1).or.(nodenei.eq.0)) exit 498 if(lakon(nelemnei)(2:3).eq.'VO') then 499 prefnew=v(2,nodenei) 500 v(2,nodenei)=v(2,node)*v(2,nodenei)/pref 501 pref=prefnew 502 else 503 pref=v(2,nodenei) 504 v(2,nodenei)=v(2,node)*ratio 505c write(*,*) nodenei,v(2,nodenei) 506 endif 507 node=nodenei 508 nelem=nelemnei 509 enddo 510 else 511! 512! the pressure in the nodes adjacent to at least one 513! vortex element were not available 514! 515! reset values; no change; 516! 517 if(ierror.eq.1) then 518 v(2,node1)=0.d0 519 else 520 v(2,node2)=0.d0 521 endif 522 cycle 523 endif 524 endif 525! 526 ichange=1 527 endif 528 endif 529 enddo 530! 531! propagation of the mass flow through the network 532! 533 loop1: do i=1,nflow 534 nelem=ieg(i) 535 indexe=ipkon(nelem) 536 nodem=kon(indexe+2) 537! 538 if(dabs(v(1,nodem)).le.1.d-30) then 539! 540! no initial mass flow given yet 541! check neighbors for mass flow (only if not 542! branch nor joint) 543! 544! first end node 545! 546 node1=kon(indexe+1) 547! 548 if(node1.ne.0) then 549 index=iponoel(node1) 550! 551 if(inoel(2,inoel(2,index)).eq.0) then 552! 553! no branch nor joint; determine neighboring element 554! 555 if(inoel(1,index).eq.nelem) then 556 neighbor=inoel(1,inoel(2,index)) 557 else 558 neighbor=inoel(1,index) 559 endif 560! 561! initial mass flow in neighboring element 562! 563 xflow=v(1,kon(ipkon(neighbor)+2)) 564! 565 if(dabs(v(1,nodem)).gt.0.d0) then 566! 567! propagate initial mass flow 568! 569 if(dabs(xflow).gt.1.d-30) then 570 if(kon(ipkon(neighbor)+1).eq.node1) then 571 v(1,nodem)=-xflow 572 else 573 v(1,nodem)=xflow 574 endif 575 ichange=1 576 cycle 577 endif 578 else 579! 580! propagate only the sign of the mass flow 581! 582 if(dabs(xflow).gt.0.d0) then 583 if(kon(ipkon(neighbor)+1).eq.node1) then 584 v(1,nodem)=-xflow 585 else 586 v(1,nodem)=xflow 587 endif 588 ichange=1 589 cycle 590 endif 591 endif 592! 593! if more than 2 elements meet: check whether 594! the flux in all but the element at stake is 595! known. If so, apply the mass balance 596! 597 fluxtot=0.d0 598 do 599 if(inoel(1,index).ne.nelem) then 600 iel=inoel(1,index) 601 nodemnei=kon(ipkon(iel)+2) 602 if(dabs(v(1,nodemnei)).le.1.d-30) exit 603! 604! convention: inflow = positive 605! 606 if(kon(ipkon(iel)+1).eq.node1) then 607 fluxtot=fluxtot-v(1,nodemnei) 608 else 609 fluxtot=fluxtot+v(1,nodemnei) 610 endif 611 endif 612 if(inoel(2,index).eq.0) then 613 v(1,nodem)=fluxtot 614 ichange=1 615 cycle loop1 616 else 617 index=inoel(2,index) 618 endif 619 enddo 620! 621 endif 622 endif 623! 624! second end node 625! 626 node2=kon(indexe+3) 627! 628 if(node2.ne.0) then 629 index=iponoel(node2) 630! 631 if(inoel(2,inoel(2,index)).eq.0) then 632! 633! no branch nor joint; determine neighboring element 634! 635 if(inoel(1,index).eq.nelem) then 636 neighbor=inoel(1,inoel(2,index)) 637 else 638 neighbor=inoel(1,index) 639 endif 640! 641! initial mass flow in neighboring element 642! 643 xflow=v(1,kon(ipkon(neighbor)+2)) 644! 645 if(dabs(v(1,nodem)).gt.0.d0) then 646! 647! propagate initial mass flow 648! 649 if(dabs(xflow).gt.1.d-30) then 650 if(kon(ipkon(neighbor)+3).eq.node2) then 651 v(1,nodem)=-xflow 652 else 653 v(1,nodem)=xflow 654 endif 655 ichange=1 656 cycle 657 endif 658 else 659! 660! propagate only the sign of the mass flow 661! 662 if(dabs(xflow).gt.0.d0) then 663 if(kon(ipkon(neighbor)+3).eq.node2) then 664 v(1,nodem)=-xflow 665 else 666 v(1,nodem)=xflow 667 endif 668 ichange=1 669 cycle 670 endif 671 endif 672! 673! if more than 2 elements meet: check whether 674! the flux in all but the element at stake is 675! known. If so, apply the mass balance 676! 677 fluxtot=0.d0 678 do 679 if(inoel(1,index).ne.nelem) then 680 iel=inoel(1,index) 681 nodemnei=kon(ipkon(iel)+2) 682 if(dabs(v(1,nodemnei)).le.1.d-30) exit 683! 684! convention: outflow = positive 685! 686 if(kon(ipkon(iel)+3).eq.node2) then 687 fluxtot=fluxtot-v(1,nodemnei) 688 else 689 fluxtot=fluxtot+v(1,nodemnei) 690 endif 691 endif 692 if(inoel(2,index).eq.0) then 693 v(1,nodem)=fluxtot 694 ichange=1 695 cycle loop1 696 else 697 index=inoel(2,index) 698 endif 699 enddo 700! 701 endif 702 endif 703 endif 704 enddo loop1 705! 706! propagation of the temperature 707! 708 do i=1,nflow 709 nelem=ieg(i) 710 indexe=ipkon(nelem) 711 node1=kon(indexe+1) 712 if(node1.eq.0) cycle 713 node2=kon(indexe+3) 714 if(node2.eq.0) cycle 715! 716! only case in which exactly 1 temperature is unknown 717! is considered 718! 719 if(((v(0,node1).ne.0.d0).and.(v(0,node2).ne.0.d0)).or. 720 & ((v(0,node1).eq.0.d0).and.(v(0,node2).eq.0.d0))) cycle 721! 722! If the element is an adiabatic gas pipe the 723! total temperature at both ends is equal 724! 725 if(lakon(nelem)(2:6).eq.'GAPFA') then 726 if(v(0,node1).eq.0.d0) then 727 v(0,node1)=v(0,node2) 728 else 729 v(0,node2)=v(0,node1) 730 endif 731 ichange=1 732 cycle 733 elseif(lakon(nelem)(2:5).eq.'GAPR') then 734! 735! total temperature change due to the rotation 736! 737 index=ielprop(nelem) 738 xl=prop(index+3) 739 r1=prop(index+8) 740 r2=prop(index+9) 741! 742 if(v(0,node1).eq.0.d0) then 743 Ti=v(0,node2) 744 else 745 Ti=v(0,node1) 746 endif 747! 748! if r2 > r1 then the centrifugal force points in 749! the direction from node1 to node2 750! 751 if(r2.lt.r1) then 752 om2=-prop(index+10)**2 753 else 754 om2=prop(index+10)**2 755 endif 756! 757 imat=ielmat(1,nelem) 758 call materialdata_tg(imat,ntmat_,Ti,shcon, 759 & nshcon,cp,r,dvi,rhcon,nrhcon,rho) 760! 761 Tt2mTt1=om2*(r1+r2)*xl/(2.d0*cp) 762! 763 if(v(0,node1).eq.0.d0) then 764 v(0,node1)=v(0,node2)-Tt2mTt1 765 else 766 v(0,node2)=v(0,node1)+Tt2mTt1 767 endif 768! 769 ichange=1 770 cycle 771 endif 772! 773 nodem=kon(indexe+2) 774! 775 if(v(1,nodem).eq.0.d0) then 776! 777! direction of mass flow unknown in the element 778! 779 cycle 780 elseif(v(1,nodem).gt.0.d0) then 781! 782! positive mass flow (i.e. going from node1 to node2) 783! 784 if(v(0,node1).eq.0.d0) cycle 785! 786! propagating the temperature to node2 787! 788 v(0,node2)=v(0,node1) 789 ichange=1 790 cycle 791 else 792! 793! negative mass flow (i.e. going from node2 to node1) 794! 795 if(v(0,node2).eq.0.d0) cycle 796! 797! propagating the temperature to node1 798! 799 v(0,node1)=v(0,node2) 800 ichange=1 801 cycle 802 endif 803 enddo 804c write(*,*) 'preinitialnet ' 805c do i=1,ntg 806c write(*,'(i10,3(1x,e11.4))') itg(i),(v(j,itg(i)),j=0,2) 807c enddo 808 if(ichange.eq.0) exit 809 enddo 810! 811! set of mass flow of +-1.d-30 to zero 812! 813 do i=1,nflow 814 nelem=ieg(i) 815 indexe=ipkon(nelem) 816 nodem=kon(indexe+2) 817 if(dabs(v(1,nodem)).eq.1.d-30) v(1,nodem)=0.d0 818 enddo 819! 820! check the pressures: set pressures to zero (i.e. no initial condtion) 821! which lie in between the neighboring pressures (i.e. at least one 822! neighbor has a lower pressure and at least one neighbor has a higher 823! pressure) => Laplace method is applied in initialnet 824! 825 loop: do i=1,ntg 826 node=itg(i) 827! 828! neighboring elements (excluding nodes which do not belong to 829! any network element or just to one network element (middle nodes) 830! 831 index=iponoel(node) 832 if((index.eq.0).or.(inoel(2,index).eq.0)) cycle 833! 834 imin=node 835 imax=node 836 xmin=v(2,node) 837 xmax=v(2,node) 838! 839 do 840 nelem=inoel(1,index) 841 if((lakon(nelem)(2:3).eq.'VO').or. 842 & (lakon(nelem)(2:5).eq.'GAPR')) cycle loop 843 indexe=ipkon(nelem) 844! 845! neighboring vertex node 846! 847 if(kon(indexe+1).ne.node) then 848 neighbor=kon(indexe+1) 849 else 850 neighbor=kon(indexe+3) 851 endif 852 if(neighbor.eq.0) cycle loop 853! 854! check its value 855! 856 if(dabs(v(2,neighbor)).lt.xmin) then 857 xmin=dabs(v(2,neighbor)) 858 imin=neighbor 859 elseif(dabs(v(2,neighbor)).gt.xmax) then 860 xmax=dabs(v(2,neighbor)) 861 imax=neighbor 862 endif 863! 864 index=inoel(2,index) 865 if(index.eq.0) exit 866 enddo 867! 868! if value lies in between the neighboring values => assign a 869! negative sign as marker 870! 871 if((imin.ne.node).and.(imax.ne.node)) then 872 v(2,node)=-v(2,node) 873 endif 874 enddo loop 875! 876! set marked values to zero => Laplace equation will be used 877! in initialnet.f 878! 879 do i=1,ntg 880 node=itg(i) 881 index=iponoel(node) 882 if((index.eq.0).or.(inoel(2,index).eq.0)) cycle 883c if(v(2,node).lt.0.d0) v(2,node)=0.d0 884 if(v(2,node).lt.0.d0) v(2,node)=-v(2,node) 885 enddo 886c write(*,*) 'preinitialnet end ' 887c do i=1,ntg 888c write(*,'(i10,3(1x,e11.4))') itg(i),(v(j,itg(i)),j=0,2) 889c enddo 890! 891! same for temperatures? 892! 893 return 894 end 895 896 897