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! construction of the b matrix 20! 21 subroutine resultnet(itg,ieg,ntg, 22 & bc,nload,sideload,nelemload,xloadact,lakon,ntmat_, 23 & v,shcon,nshcon,ipkon,kon,co,nflow, 24 & iinc,istep,dtime,ttime,time, 25 & ikforc,ilforc,xforcact,nforc,ielmat,nteq,prop,ielprop, 26 & nactdog,nacteq,iin,physcon,camt,camf,camp,rhcon,nrhcon, 27 & ipobody,ibody,xbodyact,nbody,dtheta,vold,xloadold, 28 & reltime,nmethod,set,mi,ineighe,cama,vamt,vamf,vamp,vama, 29 & nmpc,nodempc,ipompc,coefmpc,labmpc,iaxial,qat,qaf,ramt, 30 & ramf,ramp,cocon,ncocon,iponoel,inoel,iplausi) 31! 32 implicit none 33! 34 logical identity 35 character*8 lakonl,lakon(*) 36 character*20 sideload(*),labmpc(*) 37 character*81 set(*) 38! 39 integer mi(*),itg(*),ieg(*),ntg,nteq,nflow,nload, 40 & ielmat(mi(3),*),iflag,ider,iaxial,nalt,nalf, 41 & nelemload(2,*),nope,nopes,mint2d,i,j,k,m,nrhcon(*), 42 & node,imat,ntmat_,id,ifaceq(8,6),ifacet(6,4),numf, 43 & ifacew(8,5),node1,node2,nshcon(*),nelem,ig,index,konl(20), 44 & ipkon(*),kon(*),idof,ineighe(*),idir,ncocon(2,*), 45 & iinc,istep,jltyp,nfield,ikforc(*),ipobody(2,*), 46 & ilforc(*),nforc,nodem,idirf(8),ieq,nactdog(0:3,*),nbody, 47 & nacteq(0:3,*),ielprop(*),nodef(8),iin,kflag,ibody(3,*),icase, 48 & inv, index2,nmethod,nelem0,nodem0,nelem1,nodem1,nelem2, 49 & nodem2,nelemswirl,nmpc,nodempc(3,*),ipompc(*), 50 & iponoel(*),inoel(2,*),iplausi,indexe 51! 52 real*8 bc(nteq),xloadact(2,*),cp,h(2),physcon(*),r,dvi,rho, 53 & xl2(3,8),coords(3),dxsj2,temp,xi,et,weight,xsj2(3), 54 & gastemp,v(0:mi(2),*),shcon(0:3,ntmat_,*),co(3,*),shp2(7,8), 55 & field,prop(*),tg1,tg2,dtime,ttime,time,g(3),eta, 56 & xforcact(*),areaj,xflow,tvar(2),f,df(8),camt(*),camf(*), 57 & camp(*),tl2(8),cama(*),vamt,vamf,vamp,vama,term, 58 & rhcon(0:1,ntmat_,*),xbodyact(7,*),sinktemp,kappa,A,T,Tt,pt, 59 & dtheta,ts1,ts2,xs2(3,7),pt1,pt2,Qred_crit,xflow_crit, 60 & xflow0,xflow1,reltime,coefmpc(*),qat,qaf,ramt,ramf,ramp, 61 & xflow2,heat,cocon(0:6,ntmat_,*), 62 & Cp_cor,U,Ct,vold(0:mi(2),*),xloadold(2,*),bnac, 63 & xflow360,Tt1,Tt2,T1,T2,heatnod,heatfac,A2,d,l,s 64! 65 include "gauss.f" 66! 67 data ifaceq /4,3,2,1,11,10,9,12, 68 & 5,6,7,8,13,14,15,16, 69 & 1,2,6,5,9,18,13,17, 70 & 2,3,7,6,10,19,14,18, 71 & 3,4,8,7,11,20,15,19, 72 & 4,1,5,8,12,17,16,20/ 73 data ifacet /1,3,2,7,6,5, 74 & 1,2,4,5,9,8, 75 & 2,3,4,6,10,9, 76 & 1,4,3,8,10,7/ 77 data ifacew /1,3,2,9,8,7,0,0, 78 & 4,5,6,10,11,12,0,0, 79 & 1,2,5,4,7,14,10,13, 80 & 2,3,6,5,8,15,11,14, 81 & 4,6,3,1,12,15,9,13/ 82 data iflag /2/ 83! 84 kflag=2 85 ider=0 86! 87 tvar(1)=time 88 tvar(2)=ttime+time 89! 90 heatnod=0.d0 91! 92! calculating the maximum correction to the solution in the 93! present network iteration 94! 95 camt(1)=0.d0 96 camf(1)=0.d0 97 camp(1)=0.d0 98 cama(1)=0.d0 99! 100 camt(2)=0.5d0 101 camf(2)=0.5d0 102 camp(2)=0.5d0 103 cama(2)=0.5d0 104! 105! maximum change in the solution since the start of the 106! network iterations (= entering radflowload) 107! 108 vamt=0.d0 109 vamf=0.d0 110 vamp=0.d0 111 vama=0.d0 112! 113 do i=1,ntg 114 node=itg(i) 115 do j=0,3 116 if(nactdog(j,node).eq.0) cycle 117 idof=nactdog(j,node) 118! 119 if(dabs(v(j,node)).gt.1.d-30) then 120 if(dabs(bc(idof)/v(j,node)).lt.1.d-13) bc(idof)=0.d0 121 endif 122 bnac=bc(idof) 123! 124 if(j.eq.0) then 125 if(dabs(bnac).gt.dabs(camt(1))) then 126 camt(1)=bnac 127 camt(2)=node+0.5d0 128 endif 129 elseif(j.eq.1) then 130 if(dabs(bnac).gt.dabs(camf(1))) then 131 camf(1)=bnac 132 camf(2)=node+0.5d0 133 endif 134 elseif(j.eq.2) then 135 if(dabs(bnac).gt.dabs(camp(1))) then 136 camp(1)=bnac 137 camp(2)=node+0.5d0 138 endif 139 else 140 if(dabs(bnac).gt.dabs(cama(1))) then 141 cama(1)=bnac 142 cama(2)=node+0.5d0 143 endif 144 endif 145 enddo 146 enddo 147! 148! updating v 149! vold is the solution at the entry of radflowload (= at 150! the start of the network iterations) 151! 152 do i=1,ntg 153 node=itg(i) 154 do j=0,2 155 if(nactdog(j,node).eq.0) cycle 156 v(j,node)=v(j,node)+bc(nactdog(j,node))*dtheta 157! 158 if((j.eq.0).and.(dabs(v(j,node)-vold(j,node)).gt.vamt)) then 159 vamt=dabs(v(j,node)-vold(j,node)) 160 elseif((j.eq.1).and. 161 & (dabs(v(j,node)-vold(j,node)).gt.vamf)) then 162 vamf=dabs(v(j,node)-vold(j,node)) 163 elseif((j.eq.2).and. 164 & (dabs(v(j,node)-vold(j,node)).gt.vamp)) then 165 vamp=dabs(v(j,node)-vold(j,node)) 166 endif 167! 168 enddo 169 enddo 170c write(*,*) 'resultnet' 171c write(*,*) 'mass flow in node 3: ',v(1,3) 172c write(*,*) 'pressure in node 4: ',v(2,4) 173c write(*,*) 'mass flow in node 5: ',v(1,5) 174! 175! update geometry changes 176! 177 do i=1,nflow 178 if(lakon(ieg(i))(6:7).eq.'GV') then 179 index=ipkon(ieg(i)) 180 node=kon(index+2) 181 if(nactdog(3,node).eq.0) cycle 182 index=ielprop(ieg(i)) 183 v(3,node)=v(3,node)+bc(nactdog(3,node))*dtheta 184 if(v(3,node).gt.0.99999d0) then 185 v(3,node)=0.99999d0 186 elseif(v(3,node).lt.0.12501d0) then 187 v(3,node)=0.12501d0 188 endif 189 if(dabs(v(3,node)).gt.vama) vama=dabs(v(3,node)) 190! 191! Geometry factor of an ACC tube 192! for all versions of the ACC tube 193! 194 elseif(lakon(ieg(i))(2:7).eq.'ACCTUB') then 195 index=ipkon(ieg(i)) 196 node=kon(index+2) 197 if(nactdog(3,node).eq.0) cycle 198 index=ielprop(ieg(i)) 199! Interval factor unknown 200 if(nint(prop(index+1)).eq.2) then 201! Using a smaller step gets better convergence(factor 0.5) 202 v(3,node)=v(3,node)+bc(nactdog(3,node))* 203 & dtheta 204 if(v(3,node).lt.0.1d0)then 205 v(3,node) = 0.1d0 206 endif 207! Hole diameter factor unknown 208 elseif(nint(prop(index+1)).eq.3) then 209 v(3,node)=v(3,node)+bc(nactdog(3,node))*dtheta 210 if(v(3,node).lt.0.1d0)then 211 v(3,node) = 0.1d0 212 endif 213 endif 214 if(dabs(v(3,node)).gt.vama) vama=dabs(v(3,node)) 215! 216! update location of hydraulic jump 217! 218 elseif(lakon(ieg(i))(2:5).eq.'LICH') then 219 if((lakon(ieg(i))(6:7).eq.'SG').or. 220 & (lakon(ieg(i))(6:7).eq.'WE').or. 221 & (lakon(ieg(i))(6:7).eq.'DS')) then 222 index=ipkon(ieg(i)) 223 node=kon(index+2) 224 if(nactdog(3,node).eq.0) cycle 225 index=ielprop(ieg(i)) 226 if(lakon(ieg(i))(6:7).eq.'SG') then 227 eta=prop(index+4)+bc(nactdog(3,node))*dtheta 228 prop(index+4)=eta 229 nelem=nint(prop(index+7)) 230 elseif(lakon(ieg(i))(6:7).eq.'WE') then 231 eta=prop(index+4)+bc(nactdog(3,node))*dtheta 232 prop(index+4)=eta 233 nelem=nint(prop(index+7)) 234 elseif(lakon(ieg(i))(6:7).eq.'DS') then 235 eta=prop(index+7)+bc(nactdog(3,node))*dtheta 236 prop(index+7)=eta 237 nelem=nint(prop(index+9)) 238 endif 239 v(3,node)=eta 240 vama=eta 241! 242! check whether 0<=eta<=1. If not, the hydraulic jump 243! does not take place in the element itself and has to 244! be forced out of the element by adjusting the 245! water depth of one of the end nodes 246! 247 if((eta.lt.0.d0).or.(eta.gt.1.d0)) then 248 index=ipkon(nelem) 249 node1=kon(index+1) 250 nodem=kon(index+2) 251 node2=kon(index+3) 252 xflow=v(1,nodem) 253! 254! determining the temperature for the 255! material properties 256! 257 if(xflow.gt.0) then 258 if(node1.eq.0) then 259 gastemp=v(0,node2) 260 else 261 gastemp=v(0,node1) 262 endif 263 else 264 if(node2.eq.0) then 265 gastemp=v(0,node1) 266 else 267 gastemp=v(0,node2) 268 endif 269 endif 270! 271 imat=ielmat(1,nelem) 272! 273 call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon, 274 & cp,r,dvi,rhcon,nrhcon,rho) 275! 276 do j=1,3 277 g(j)=0.d0 278 enddo 279 if(nbody.gt.0) then 280 index=nelem 281 do 282 j=ipobody(1,index) 283 if(j.eq.0) exit 284 if(ibody(1,j).eq.2) then 285 g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j) 286 g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j) 287 g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j) 288 endif 289 index=ipobody(2,index) 290 if(index.eq.0) exit 291 enddo 292 endif 293! 294 kflag=3 295 call flux(node1,node2,nodem,nelem,lakon,kon,ipkon, 296 & nactdog,identity, 297 & ielprop,prop,kflag,v,xflow,f,nodef,idirf,df, 298 & cp,r,rho,physcon,g,co,dvi,numf,vold,set,shcon, 299 & nshcon,rhcon,nrhcon,ntmat_,mi,ider,ttime,time, 300 & iaxial,iplausi) 301 kflag=2 302! 303 endif 304 endif 305 endif 306 enddo 307! 308c write(*,*) 'resultnet' 309c do i=1,ntg 310c node=itg(i) 311c write(*,'(i10,3(1x,e11.4))') node,(v(j,node),j=0,2) 312c enddo 313! 314! testing the validity of the pressures 315! 316 do i=1,ntg 317 node=itg(i) 318 if(v(2,node).lt.0) then 319 write(*,*) 'wrong pressure; node ',node,v(2,node) 320 iin=0 321 return 322 endif 323 enddo 324! 325! testing validity of temperatures 326! 327 do i=1,ntg 328 node=itg(i) 329 if(v(0,node).lt.0) then 330 write(*,*) 'wrong temperature; node ',node,v(0,node) 331 iin=0 332 return 333 endif 334 enddo 335! 336! testing the validity of the solution for branches elements 337! and restrictor. Since the element properties are dependent on 338! a predefined flow direction a change of this will lead to 339! wrong head losses 340! 341 do i=1,nflow 342 nelem=ieg(i) 343 if ((lakon(nelem)(4:5).eq.'ATR').or. 344 & (lakon(nelem)(4:5).eq.'RTA')) then 345 xflow=v(1,kon(ipkon(nelem)+2)) 346 if(xflow.lt.0d0)then 347 Write(*,*)'*WARNING in resultnet.f' 348 write(*,*)'Element',nelem,'of TYPE ABSOLUTE TO RELATIVE' 349 write(*,*)'The flow direction is no more conform ' 350 write(*,*)'to element definition' 351 write(*,*)'Check the pertinence of the results' 352 endif 353 elseif(lakon(nelem)(2:3).eq.'RE') then 354! 355 if(lakon(nelem)(4:5).ne.'BR') then 356 nodem=kon(ipkon(nelem)+2) 357 xflow=v(1,nodem) 358 if (xflow.lt.0) then 359 Write(*,*)'*WARNING in resultnet.f' 360 write(*,*)'Element',nelem,'of TYPE RESTRICTOR' 361 write(*,*)'The flow direction is no more conform ' 362 write(*,*)'to element definition' 363 write(*,*)'Check the pertinence of the results' 364 endif 365! 366 elseif(lakon(nelem)(4:5).eq.'BR') then 367 index=ielprop(nelem) 368! 369 nelem0=nint(prop(index+1)) 370 nodem0=kon(ipkon(nelem0)+2) 371 xflow0=v(1,nodem0) 372! 373 nelem1=nint(prop(index+2)) 374 nodem1=kon(ipkon(nelem1)+2) 375 xflow1=v(1,nodem1) 376! 377 nelem2=nint(prop(index+3)) 378 nodem2=kon(ipkon(nelem2)+2) 379 xflow2=v(1,nodem2) 380! 381 if((xflow0.lt.0).or.(xflow1.lt.0).or.(xflow2.lt.0)) then 382 Write(*,*)'*WARNING in resultnet.f' 383 write(*,*)'Element',nelem,'of TYPE BRANCH' 384 write(*,*)'The flow direction is no more conform ' 385 write(*,*)'to element definition' 386 write(*,*)'Check the pertinence of the results' 387 endif 388 endif 389 endif 390 enddo 391! 392! determining the static temperature and checking for 393! critical conditions (only for non-chambers, i.e. for 394! nodes purely belonging to gas pipes or restrictors except 395! restrictor wall orifice) 396! 397 iplausi=1 398 do i=1,ntg 399 node=itg(i) 400 nelem=ineighe(i) 401 if(nelem.eq.-1) then 402 v(3,node)=v(0,node) 403 elseif(nelem.gt.0) then 404! 405 nodem=kon(ipkon(nelem)+2) 406 T=v(3,node) 407 Tt=v(0,node) 408 pt=v(2,node) 409 xflow=v(1,nodem) 410! 411 icase=0 412 inv=1 413 imat=ielmat(1,nelem) 414 call materialdata_tg(imat,ntmat_,v(3,node), 415 & shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho) 416! 417 index=ielprop(nelem) 418 indexe=ipkon(nelem) 419! 420 kappa=(cp/(cp-R)) 421! 422 if(lakon(nelem)(2:5).eq.'GAPF') then 423! 424! distinguish between standard pipes and flexible 425! pipes 426! 427 if((lakon(nelem)(6:8).eq.'AFR').or. 428 & (lakon(nelem)(6:8).eq.'ARL').or. 429 & (lakon(nelem)(6:8).eq.'ARG').or. 430 & (lakon(nelem)(6:8).eq.'IFR').or. 431 & (lakon(nelem)(6:8).eq.'IRL')) then 432 call calcgeomelemnet(vold,co,prop,lakon,nelem, 433 & ttime,time,ielprop,mi,A,A2,d,l,s) 434 else 435 A=prop(index+1) 436 endif 437! 438 if(lakon(nelem)(2:6).eq.'GAPFA') then 439 icase=0 440 elseif(lakon(nelem)(2:6).eq.'GAPFI') then 441 icase=1 442 endif 443! 444 elseif(lakon(nelem)(2:5).eq.'GAPR') then 445! 446! rotating adiabatic gas pipe with variable cross section: 447! check whether section 1 or 2 448! 449 if(kon(ipkon(nelem)+1).eq.node) then 450 A=prop(index+1) 451 else 452 A=prop(index+2) 453 endif 454 icase=0 455! 456 elseif(lakon(nelem)(2:3).eq.'RE') then 457 node1=kon(indexe+1) 458 node2=kon(indexe+3) 459! 460 if(lakon(nelem)(4:5).eq.'EX') then 461 if(lakon(nint(prop(index+4)))(2:6).eq.'GAPFA') then 462 icase=0 463 elseif(lakon(nint(prop(index+4)))(2:6).eq.'GAPFI')then 464 icase=1 465 endif 466 else 467 icase=0 468 endif 469! 470! defining the sections 471! 472 if(lakon(nelem)(4:5).eq.'BE') then 473 A=prop(index+1) 474! 475 elseif(lakon(nelem)(4:5).eq.'BR') then 476 if(lakon(nelem)(4:6).eq.'BRJ') then 477 if(nelem.eq.nint(prop(index+2)))then 478 A=prop(index+5) 479 elseif(nelem.eq.nint(prop(index+3))) then 480 A=prop(index+6) 481 endif 482 elseif(lakon(nelem)(4:6).eq.'BRS') then 483 if(nelem.eq.nint(prop(index+2)))then 484 A=prop(index+5) 485 elseif(nelem.eq.nint(prop(index+3))) then 486 A=prop(index+6) 487 endif 488 endif 489! 490 else 491! 492 if((lakon(nelem)(2:5).eq.'RBSF')) then 493 call calcgeomelemnet(vold,co,prop,lakon,nelem, 494 & ttime,time,ielprop,mi,A,A2,d,l,s) 495 if(node.eq.node2) then 496 A=A2 497 endif 498 else 499 if(node.eq.node1) then 500 A=prop(index+1) 501 elseif(node.eq.node2) then 502 A=prop(index+2) 503 endif 504 endif 505 endif 506 elseif(lakon(nelem)(2:3).eq.'UP') then 507! 508! user elements whose names start with UP are assumed 509! to be Pipe-like (static and total temperatures differ) 510! 511 call calcgeomelemnet(vold,co,prop,lakon,nelem, 512 & ttime,time,ielprop,mi,A,A2,d,l,s) 513c A=prop(index+1) 514 icase=0 515 endif 516! 517 xflow360=dabs(xflow*iaxial) 518 call ts_calc(xflow360,Tt,pt,kappa,r,A,T,icase) 519! 520 v(3,node)=T 521! 522 if(lakon(nelem)(2:3).eq.'RE') then 523! 524! check whether the mass flow does not exceed 525! critical conditions (only for restrictor elements) 526! 527 if(xflow.lt.0d0) then 528 inv=-1 529 else 530 inv=1 531 endif 532! 533 if(icase.eq.0) then 534 Qred_crit=dsqrt(kappa/R)*(1.d0+0.5d0*(kappa-1.d0)) 535 & **(-0.5d0*(kappa+1.d0)/(kappa-1.d0)) 536 else 537 Qred_crit=dsqrt(1/R)*(1.d0+0.5d0*(kappa-1.d0)/kappa) 538 & **(-0.5d0*(kappa+1.d0)/(kappa-1.d0)) 539 endif 540 xflow_crit=Qred_crit*pt*A/dsqrt(Tt) 541! 542 if(xflow360.ge.xflow_crit) then 543! 544! next 3 lines: change on 20200929 (Thomas Petzke) 545! 546 if(dabs((xflow_crit-xflow360)/xflow_crit).gt.1.d-5) 547 & then 548 iplausi=0 549 endif 550! 551 v(1,nodem)=inv*xflow_crit/iaxial 552 if(icase.eq.1) then 553! 554 if(nactdog(0,node2).ne.0) then 555 node1=kon(indexe+1) 556 node2=kon(indexe+3) 557 v(3,node2)=v(3,node1) 558 v(0,node2)=v(3,node2) 559 & *(1+0.5d0*(kappa-1)/kappa) 560 561 endif 562 endif 563 endif 564 endif 565! 566 endif 567! 568 enddo 569! 570! testing if the flow direction is conform to the pressure 571! gradient 572! 573c iplausi=1 574 do i=1,nflow 575 nelem=ieg(i) 576 indexe=ipkon(nelem) 577 node1=kon(indexe+1) 578 nodem=kon(indexe+2) 579 node2=kon(indexe+3) 580! 581 if((node1.ne.0).and.(node2.ne.0)) then 582! 583! the mass flow rates must always be oriented in the direction 584! of the decreasing pressure gradient except for ACCTUBE, 585! CROSS SPLIT, MOEHRING, ROTATING CAVITY, RADIAL OUTFLOW and 586! INLET, RIMSEAL, S-PUMP and VORTEX 587! 588 if((lakon(nelem)(2:8).ne.'ACCTUBE').and. 589 & (lakon(nelem)(2:5).ne.'ATR').and. 590 & (lakon(nelem)(2:5).ne.'RTA').and. 591 & (lakon(nelem)(2:5).ne.'CROS').and. 592 & (lakon(nelem)(2:3).ne.'LI').and. 593 & (lakon(nelem)(2:3).ne.'LP').and. 594 & (lakon(nelem)(2:4).ne.'MRG').and. 595 & (lakon(nelem)(2:4).ne.'RCV').and. 596 & (lakon(nelem)(2:3).ne.'RO').and. 597 & (lakon(nelem)(2:5).ne.'RIMS').and. 598 & (lakon(nelem)(2:5).ne.'FDPF').and. 599 & (lakon(nelem)(2:5).ne.'FCVF').and. 600 & (lakon(nelem)(2:5).ne.'MFPC').and. 601 & (lakon(nelem)(2:6).ne.'SPUMP').and. 602 & (lakon(nelem)(2:5).ne.'GAPR').and. 603 & (lakon(nelem)(2:3).ne.'VO')) then 604 if(nactdog(1,nodem).ne.0) then 605 if(((v(1,nodem).gt.0.d0).and. 606 & (v(2,node1)/v(2,node2).lt.(1.d0-1.d-10))).or. 607 & ((v(1,nodem).lt.0.d0).and. 608 & (v(2,node2)/v(2,node1).lt.(1.d0-1.d-10)))) then 609c if(((v(1,nodem).gt.0.d0).and. 610c & (v(2,node1).lt.v(2,node2))).or. 611c & ((v(1,nodem).lt.0.d0).and. 612c & (v(2,node1).gt.v(2,node2)))) then 613 iplausi=0 614 write(*,*) '*WARNING in resultnet: the flow' 615 write(*,*) ' direction is not conform' 616 write(*,*) ' to the pressure gradient' 617 write(*,*) ' element',nelem 618 endif 619 endif 620 endif 621 endif 622 enddo 623! 624! reinitialisation of the Bc matrix 625! 626 do i=1,nteq 627 bc(i)=0.d0 628 enddo 629! 630 qat=0.d0 631 qaf=0.d0 632 nalt=0 633 nalf=0 634! 635! determining the residual 636! 637 do i=1,nflow 638 nelem=ieg(i) 639 index=ipkon(nelem) 640 node1=kon(index+1) 641 nodem=kon(index+2) 642 node2=kon(index+3) 643 xflow=v(1,nodem) 644! 645! gas: the property temperature is the static temperature 646! 647 if((lakon(nelem)(2:3).ne.'LP').and. 648 & (lakon(nelem)(2:3).ne.'LI')) then 649 if(node1.eq.0) then 650 tg1=v(0,node2) 651 tg2=tg1 652 ts1=v(3,node2) 653 ts2=ts1 654 elseif(node2.eq.0) then 655 tg1=v(0,node1) 656 tg2=tg1 657 ts1=v(3,node1) 658 ts2=ts1 659 else 660 tg1=v(0,node1) 661 tg2=v(0,node2) 662 ts1=v(3,node1) 663 ts2=v(3,node2) 664 endif 665 gastemp=(ts1+ts2)/2.d0 666 else 667! 668! liquid: only one temperature 669! 670 if(xflow.gt.0) then 671 if(node1.eq.0) then 672 gastemp=v(0,node2) 673 else 674 gastemp=v(0,node1) 675 endif 676 else 677 if(node2.eq.0) then 678 gastemp=v(0,node1) 679 else 680 gastemp=v(0,node2) 681 endif 682 endif 683! 684 if(node1.eq.0) then 685 tg2=v(0,node2) 686 tg1=tg2 687 elseif(node2.eq.0) then 688 tg1=v(0,node1) 689 tg2=tg1 690 else 691 tg1=v(0,node1) 692 tg2=v(0,node2) 693 endif 694 endif 695! 696 imat=ielmat(1,nelem) 697! 698 call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,cp,r,dvi, 699 & rhcon,nrhcon,rho) 700 kappa=Cp/(Cp-R) 701! 702! Definitions of the constant for isothermal flow elements 703! 704 if(lakon(nelem)(2:6).eq.'GAPFI') then 705 if((node1.ne.0).and.(node2.ne.0)) then 706! 707 if((lakon(nelem)(7:8).eq.'FR').or. 708 & (lakon(nelem)(7:8).eq.'RL')) then 709 call calcgeomelemnet(vold,co,prop,lakon,nelem, 710 & ttime,time,ielprop,mi,A,A2,d,l,s) 711 else 712 A=prop(ielprop(nelem)+1) 713 endif 714! 715 pt1=v(2,node1) 716 pt2=v(2,node2) 717! 718 xflow360=xflow*iaxial 719 icase=1 720! 721 if(pt1.ge.pt2)then 722 if(dabs(tg2/ts2-(1.d0+0.5d0*(kappa-1)/kappa)).lt.1d-5) 723 & then 724 pt2=dabs(xflow360)*dsqrt(Tg2*R)/A 725 & *(1.d0+0.5d0*(kappa-1.d0)/kappa) 726 & **(0.5d0*(kappa+1.d0)/(kappa-1.d0)) 727! 728 endif 729 Tt1=v(0,node1)-physcon(1) 730 Tt2=v(0,node2)-physcon(1) 731 T1=v(3,node1) 732 call ts_calc(xflow360,Tt1,pt1,kappa,r,A,T1,icase) 733 call ts_calc(xflow360,Tt2,pt2,kappa,r,A,T2,icase) 734 v(3,node1)=T1 735 v(3,node2)=T2 736 else 737 pt1=v(2,node2) 738 pt2=v(2,node1) 739 if(v(2,nodem).ge.(pt2/pt1))then 740 pt2=v(2,nodem)*pt1 741 endif 742! 743 Tt1=v(0,node2)-physcon(1) 744 call ts_calc(xflow360,Tt1,pt1,kappa,r,A,T1,icase) 745 Tt2=v(0,node1)-physcon(1) 746 call ts_calc(xflow360,Tt2,pt2,kappa,r,A,T2,icase) 747 endif 748! 749 endif 750 endif 751! 752 if(node1.ne.0) then 753! 754! energy equation contribution node1 755! 756 if (nacteq(0,node1).ne.0) then 757 ieq=nacteq(0,node1) 758! 759 if(nacteq(3,node1).eq.0) then 760 if (xflow.lt.0d0)then 761 term=cp*(tg1-tg2)*xflow 762 bc(ieq)=bc(ieq)+term 763 qat=qat+dabs(term) 764 nalt=nalt+1 765 endif 766! 767 elseif(lakon(nelem)(2:6).eq.'GAPFI') then 768 if((nacteq(3,node1).eq.node2)) then 769! 770 bc(ieq)=(T2-T1) 771! 772 endif 773 endif 774 endif 775! 776! mass equation contribution node1 777! 778 if (nacteq(1,node1).ne.0) then 779 ieq=nacteq(1,node1) 780 bc(ieq)=bc(ieq)-xflow 781 qaf=qaf+dabs(xflow) 782 nalf=nalf+1 783 endif 784 endif 785! 786 if(node2.ne.0) then 787! 788! energy equation contribution node2 789! 790 if (nacteq(0,node2).ne.0) then 791 ieq=nacteq(0,node2) 792! 793 if(nacteq(3,node2).eq.0) then 794 if (xflow.gt.0d0)then 795 term=cp*(tg2-tg1)*xflow 796 bc(ieq)=bc(ieq)-term 797 qat=qat+dabs(term) 798 nalt=nalt+1 799 endif 800! 801 elseif(lakon(nelem)(2:6).eq.'GAPFI') then 802 if(nacteq(3,node2).eq.node1) then 803! 804 bc(ieq)=(T2-T1) 805! 806 endif 807 endif 808 endif 809! 810! mass equation contribution node2 811! only in the case of an ACCTUBE but not in the case of an ACCTUBO 812! 813 if(lakon(nelem)(2:8).ne.'ACCTUBE') then 814 if (nacteq(1,node2).ne.0) then 815 ieq=nacteq(1,node2) 816 bc(ieq)=bc(ieq)+xflow 817 qaf=qaf+dabs(xflow) 818 nalf=nalf+1 819 endif 820 else 821 if (nacteq(1,node2).ne.0) then 822 if(nelem.ne.nint(prop(ielprop(nelem)+14))) then 823 ieq=nacteq(1,node2) 824 bc(ieq)=bc(ieq)+xflow-v(0,nodem) 825 qaf=qaf+dabs(xflow) 826 nalf=nalf+1 827 else 828 ieq=nacteq(1,node2) 829 bc(ieq)=bc(ieq)+xflow 830 qaf=qaf+dabs(xflow) 831 nalf=nalf+1 832 endif 833 endif 834 endif 835 endif 836! 837! element equation 838! 839 if (nacteq(2,nodem).ne.0) then 840 ieq=nacteq (2,nodem) 841! 842! for liquids: determine the gravity vector 843! 844 if((lakon(nelem)(2:3).eq.'LI').or. 845 & (lakon(nelem)(2:3).eq.'LP')) then 846 do j=1,3 847 g(j)=0.d0 848 enddo 849 if(nbody.gt.0) then 850 index=nelem 851 do 852 j=ipobody(1,index) 853 if(j.eq.0) exit 854 if(ibody(1,j).eq.2) then 855 g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j) 856 g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j) 857 g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j) 858 endif 859 index=ipobody(2,index) 860 if(index.eq.0) exit 861 enddo 862 endif 863 endif 864! 865 call flux(node1,node2,nodem,nelem,lakon,kon,ipkon, 866 & nactdog,identity, 867 & ielprop,prop,kflag,v,xflow,f,nodef,idirf,df, 868 & cp,r,rho,physcon,g,co,dvi,numf,vold,set,shcon, 869 & nshcon,rhcon,nrhcon,ntmat_,mi,ider,ttime,time, 870 & iaxial,iplausi) 871 bc(ieq)=-f 872 endif 873 enddo 874! 875! convection with the walls: contribution to the energy equations 876! 877 do i=1,nload 878 if(sideload(i)(3:4).eq.'FC') then 879 nelem=nelemload(1,i) 880 index=ipkon(nelem) 881 if(index.lt.0) cycle 882 lakonl=lakon(nelem) 883 node=nelemload(2,i) 884 ieq=nacteq(0,node) 885 if(ieq.eq.0) then 886 cycle 887 endif 888! 889 call nident(itg,node,ntg,id) 890! 891! calculate the area 892! 893 read(sideload(i)(2:2),'(i1)') ig 894! 895! number of nodes and integration points in the face 896! 897 if(lakonl(4:4).eq.'2') then 898 nope=20 899 nopes=8 900 elseif(lakonl(4:4).eq.'8') then 901 nope=8 902 nopes=4 903 elseif(lakonl(4:5).eq.'10') then 904 nope=10 905 nopes=6 906 elseif(lakonl(4:4).eq.'4') then 907 nope=4 908 nopes=3 909 elseif(lakonl(4:5).eq.'15') then 910 nope=15 911 else 912 nope=6 913 endif 914! 915 if(lakonl(4:5).eq.'8R') then 916 mint2d=1 917 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) 918 & then 919 if(lakonl(7:7).eq.'A') then 920 mint2d=2 921 else 922 mint2d=4 923 endif 924 elseif(lakonl(4:4).eq.'2') then 925 mint2d=9 926 elseif(lakonl(4:5).eq.'10') then 927 mint2d=3 928 elseif(lakonl(4:4).eq.'4') then 929 mint2d=1 930 endif 931! 932 if(lakonl(4:4).eq.'6') then 933 mint2d=1 934 if(ig.le.2) then 935 nopes=3 936 else 937 nopes=4 938 endif 939 endif 940 if(lakonl(4:5).eq.'15') then 941 if(ig.le.2) then 942 mint2d=3 943 nopes=6 944 else 945 mint2d=4 946 nopes=8 947 endif 948 endif 949! 950! connectivity of the element 951! 952 do k=1,nope 953 konl(k)=kon(index+k) 954 enddo 955! 956! coordinates of the nodes belonging to the face 957! 958 if((nope.eq.20).or.(nope.eq.8)) then 959 do k=1,nopes 960 tl2(k)=v(0,konl(ifaceq(k,ig))) 961 do j=1,3 962 xl2(j,k)=co(j,konl(ifaceq(k,ig)))+ 963 & v(j,konl(ifaceq(k,ig))) 964 enddo 965 enddo 966 elseif((nope.eq.10).or.(nope.eq.4)) then 967 do k=1,nopes 968 tl2(k)=v(0,konl(ifacet(k,ig))) 969 do j=1,3 970 xl2(j,k)=co(j,konl(ifacet(k,ig)))+ 971 & v(j,konl(ifacet(k,ig))) 972 enddo 973 enddo 974 else 975 do k=1,nopes 976 tl2(k)=v(0,konl(ifacew(k,ig))) 977 do j=1,3 978 xl2(j,k)=co(j,konl(ifacew(k,ig)))+ 979 & v(j,konl(ifacew(k,ig))) 980 enddo 981 enddo 982 endif 983! 984! integration to obtain the area and the mean 985! temperature 986! 987 do m=1,mint2d 988 if((lakonl(4:5).eq.'8R').or. 989 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 990 xi=gauss2d1(1,m) 991 et=gauss2d1(2,m) 992 weight=weight2d1(m) 993 elseif((lakonl(4:4).eq.'8').or. 994 & (lakonl(4:6).eq.'20R').or. 995 & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then 996 xi=gauss2d2(1,m) 997 et=gauss2d2(2,m) 998 weight=weight2d2(m) 999 elseif(lakonl(4:4).eq.'2') then 1000 xi=gauss2d3(1,m) 1001 et=gauss2d3(2,m) 1002 weight=weight2d3(m) 1003 elseif((lakonl(4:5).eq.'10').or. 1004 & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then 1005 xi=gauss2d5(1,m) 1006 et=gauss2d5(2,m) 1007 weight=weight2d5(m) 1008 elseif((lakonl(4:4).eq.'4').or. 1009 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 1010 xi=gauss2d4(1,m) 1011 et=gauss2d4(2,m) 1012 weight=weight2d4(m) 1013 endif 1014! 1015 if(nopes.eq.8) then 1016 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 1017 elseif(nopes.eq.4) then 1018 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 1019 elseif(nopes.eq.6) then 1020 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 1021 else 1022 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 1023 endif 1024! 1025 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 1026 & xsj2(3)*xsj2(3)) 1027 areaj=dxsj2*weight 1028! 1029 temp=0.d0 1030 do k=1,3 1031 coords(k)=0.d0 1032 enddo 1033 do j=1,nopes 1034 temp=temp+tl2(j)*shp2(4,j) 1035 do k=1,3 1036 coords(k)=coords(k)+xl2(k,j)*shp2(4,j) 1037 enddo 1038 enddo 1039! 1040 sinktemp=v(0,node) 1041 if(sideload(i)(5:6).ne.'NU') then 1042 h(1)=xloadact(1,i) 1043 else 1044 read(sideload(i)(2:2),'(i1)') jltyp 1045 jltyp=jltyp+10 1046 call film(h,sinktemp,temp,istep, 1047 & iinc,tvar,nelem,m,coords,jltyp,field,nfield, 1048 & sideload(i),node,areaj,v,mi,ipkon,kon,lakon, 1049 & iponoel,inoel,ielprop,prop,ielmat,shcon,nshcon, 1050 & rhcon,nrhcon,ntmat_,cocon,ncocon, 1051 & ipobody,xbodyact,ibody,heatnod,heatfac) 1052 endif 1053 if(lakonl(5:7).eq.'0RA') then 1054 term=2.d0*((temp-sinktemp)*h(1)+heatnod)*dxsj2*weight 1055 bc(ieq)=bc(ieq)+term 1056 qat=qat+dabs(term) 1057 nalt=nalt+1 1058 else 1059 term=((temp-sinktemp)*h(1)+heatnod)*dxsj2*weight 1060 bc(ieq)=bc(ieq)+term 1061 qat=qat+dabs(term) 1062 nalt=nalt+1 1063 endif 1064 enddo 1065 endif 1066 enddo 1067! 1068! prescribed heat generation: contribution to the energy equations 1069! 1070 do i=1,ntg 1071 node=itg(i) 1072 idof=8*(node-1) 1073 call nident(ikforc,idof,nforc,id) 1074 if(id.gt.0) then 1075 if(ikforc(id).eq.idof) then 1076 ieq=nacteq(0,node) 1077 if(ieq.ne.0) then 1078 term=xforcact(ilforc(id)) 1079 bc(ieq)=bc(ieq)+term 1080 qat=qat+dabs(term) 1081 nalt=nalt+1 1082 endif 1083 cycle 1084 endif 1085 endif 1086 enddo 1087! 1088! in the case of forced vortices, when temperature change 1089! is required, additional heat input is added in the energy 1090! equation for the downstream node 1091! 1092 do i=1,nflow 1093 nelem=ieg(i) 1094! 1095! free vortex and no temperature change 1096! 1097 if(lakon(nelem)(2:3).eq.'VO') then 1098 if(lakon(nelem)(4:5).eq.'FR') then 1099 if((nint(prop(ielprop(nelem)+8)).eq.0).or. 1100 & (nint(prop(ielprop(nelem)+8)).eq.1)) cycle 1101 elseif(lakon(nelem)(4:5).eq.'FO') then 1102 if(nint(prop(ielprop(nelem)+6)).eq.0) cycle 1103 endif 1104! 1105 call calcheatnet(nelem,lakon,ipkon,kon,v,ielprop,prop, 1106 & ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,ipobody,ibody, 1107 & xbodyact,mi,nacteq,bc,qat,nalt) 1108! 1109 elseif(lakon(nelem)(2:5).eq.'GAPR') then 1110! 1111 call calcheatnet(nelem,lakon,ipkon,kon,v,ielprop,prop, 1112 & ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,ipobody,ibody, 1113 & xbodyact,mi,nacteq,bc,qat,nalt) 1114! 1115 elseif(lakon(nelem)(2:2).eq.'U') then 1116 call calcheatnet(nelem,lakon,ipkon,kon,v,ielprop,prop, 1117 & ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,ipobody,ibody, 1118 & xbodyact,mi,nacteq,bc,qat,nalt) 1119 else 1120 cycle 1121 endif 1122 enddo 1123! 1124! transfer element ABSOLUTE TO RELATIVE / RELATIVE TO ABSOLUTE 1125! 1126 do i=1,nflow 1127 nelem=ieg(i) 1128! 1129 if((lakon(nelem)(2:4).eq.'ATR').or. 1130 & (lakon(nelem)(2:4).eq.'RTA')) then 1131! 1132 nodem=kon(ipkon(nelem)+2) 1133 xflow=v(1,nodem) 1134 if(xflow.gt.0d0) then 1135 node1=kon(ipkon(nelem)+1) 1136 node2=kon(ipkon(nelem)+3) 1137 else 1138 node1=kon(ipkon(nelem)+1) 1139 node2=kon(ipkon(nelem)+3) 1140 endif 1141! 1142! computing temperature corrected Cp=Cp(T) coefficient 1143! 1144 Tg1=v(0,node1) 1145 Tg2=v(0,node2) 1146 if((lakon(nelem)(2:3).ne.'LP').and. 1147 & (lakon(nelem)(2:3).ne.'LI')) then 1148 gastemp=(tg1+tg2)/2.d0 1149 else 1150 if(xflow.gt.0) then 1151 gastemp=tg1 1152 else 1153 gastemp=tg2 1154 endif 1155 endif 1156! 1157 imat=ielmat(1,nelem) 1158 call materialdata_tg(imat,ntmat_,gastemp, 1159 & shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho) 1160! 1161 call cp_corrected(cp,Tg1,Tg2,cp_cor) 1162! 1163 index=ielprop(nelem) 1164 U=prop(index+1) 1165 ct=prop(index+2) 1166! 1167! if a swirl element was given by the user, this takes 1168! precendence to a value of ct 1169! 1170 if(nint(prop(index+3)).ne.0) then 1171 nelemswirl=nint(prop(index+3)) 1172 index2=ielprop(nelemswirl) 1173! 1174! previous element is a preswirl nozzle 1175! 1176 if(lakon(nelemswirl)(2:5).eq.'ORPN') then 1177 ct=prop(index2+4) 1178! 1179! previous element is a forced vortex 1180! 1181 elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then 1182 ct=prop(index2+7) 1183! 1184! previous element is a free vortex 1185! 1186 elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then 1187 ct=prop(index2+9) 1188 endif 1189 endif 1190! 1191 if(lakon(nelem)(2:4).eq.'ATR') then 1192 heat=Cp/Cp_cor*(0.5d0*(U**2-2d0*U*Ct)*xflow) 1193! 1194 elseif(lakon(nelem)(2:4).eq.'RTA') then 1195 heat=Cp/Cp_cor*(-0.5d0*(U**2-2d0*U*Ct)*xflow) 1196 endif 1197! 1198! including the resulting additional heat flux in the energy equation 1199! 1200 if(xflow.gt.0d0)then 1201 ieq=nacteq(0,node2) 1202 if(ieq.ne.0) then 1203 bc(ieq)=bc(ieq)+heat 1204 qat=qat+dabs(heat) 1205 nalt=nalt+1 1206 endif 1207 else 1208 ieq=nacteq(0,node1) 1209 if(ieq.ne.0) then 1210 bc(ieq)=bc(ieq)+heat 1211 qat=qat+dabs(heat) 1212 nalt=nalt+1 1213 endif 1214 endif 1215 endif 1216 enddo 1217! 1218! additional multiple point constraints 1219! 1220 j=nteq+1 1221 do i=nmpc,1,-1 1222 if(labmpc(i)(1:7).ne.'NETWORK') cycle 1223 j=j-1 1224 if(labmpc(i)(8:8).eq.' ') then 1225! 1226! linear equation 1227! 1228 index=ipompc(i) 1229! 1230 do 1231 node=nodempc(1,index) 1232 idir=nodempc(2,index) 1233 bc(j)=bc(j)-v(idir,node)*coefmpc(index) 1234 index=nodempc(3,index) 1235 if(index.eq.0) exit 1236 enddo 1237 else 1238! 1239! user-defined network equation 1240! 1241 call networkmpc_rhs(i,ipompc,nodempc,coefmpc,labmpc, 1242 & v,bc,j,mi,ipkon,kon,lakon,iponoel, 1243 & inoel,ielprop,prop,ielmat, 1244 & shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon) 1245 endif 1246 enddo 1247! 1248! determining the typical energy flow 1249! 1250 if(nalt.gt.0) qat=qat/nalt 1251! 1252! determining the typical mass flow 1253! 1254 if(nalf.gt.0) qaf=qaf/nalf 1255! 1256! max. residual energy flow, residual mass flow 1257! and residuals from the element equations 1258! 1259 ramt=0.d0 1260 ramf=0.d0 1261 ramp=0.d0 1262 do i=1,ntg 1263 node=itg(i) 1264 if(nacteq(0,node).ne.0) then 1265 ramt=max(ramt,dabs(bc(nacteq(0,node)))) 1266 endif 1267 if(nacteq(1,node).ne.0) then 1268 ramf=max(ramf,dabs(bc(nacteq(1,node)))) 1269 endif 1270 if(nacteq(2,node).ne.0) then 1271 ramp=max(ramp,dabs(bc(nacteq(2,node)))) 1272 endif 1273 enddo 1274! 1275c write(*,*) 'bc in resultnet' 1276c do i=1,3 1277c write(*,'(1x,e11.4)') bc(i) 1278c enddo 1279! 1280 return 1281 end 1282