1 : C BCOND2 SOURCE DURA 06/07/11 21:16:48 5481 2 : SUBROUTINE BCOND2 (X,DX,XL,RUG,XW,XN,TN,EN,BN,KIMP,PSLIM,REL, 3 : & P1,PS1,T1,Y1,QAE,QEE1,PHI1, 4 : & P2,T2,Y2,U2,QEE2,PHI2,QW2,RE,H,PSQ,RINDEX, 5 : & NPP,ITP,PF,PP,DPF,DPP,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4) 6 : 7 : IMPLICIT INTEGER(I-N) 8 : IMPLICIT REAL*8 (A-H,O-Z) 9 : C operateur FUITE option gaz reel 10 : C>>> air + vapeur + liquide 11 : C QAE : QA*E*B invariant air flowrate 12 : C QEEi : QEEi*E*B invariant water flowrate 13 : C XW=1.(0.) with (without) wall condensation 14 : C indice 1,2 : entree et sortie de troncon 15 : C RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4 : coef lois de frot utilisateur 16 : C H correspond a l'energie totale cedee a la paroi 17 : C H=HM+RAPH*HW transmis 18 : C H=HM+HW pour le calcul de T2 19 : C 20 : DIMENSION XN(NPP),TN(NPP),EN(NPP),BN(NPP) 21 : 22 : IF(KIMP.GE.1) THEN 23 : write(6,*) 'entree bcond2 X T1 Y1 PHI1= ',X,T1,Y1,PHI1 24 : ENDIF 25 : RINDEX=1.0 26 : 27 : ERPS1 = ABS (1.-(PS1/PSATT0(T1))) 28 : IF (ERPS1.GT.1D-4) THEN 29 : WRITE(6,*) ' bcond :erreur sur PS1 ',ERPS1 30 : ENDIF 31 : 32 : C IMPLEMENTATION DES PROPRIETES PHYSIQUES 32 : CALL BPHYS(T0,P0,RA,RV,CA,XCV,XCL,XLAT0,XROL,XKL,XKT,XREL) 33 : 34 : C TEMPERATURE PAROI A ENTREE (KELVIN) 34 : CALL BTPAR(XN,TN,X,TP,NPP,ITP) 35 : TP=TP+T0 36 : 37 : C EPAISSEUR DE FISSURE EN ENTREE 37 : CALL BTPAR(XN,EN,X,E1,NPP,ITP) 38 : C etendue a l'entree 38 : CALL BTPAR(XN,BN,X,B1,NPP,ITP) 39 : 40 : C------------------------------------- 41 : C DEFINITION DES DIFFERENTS DEBITS (homogenes a ro.u) 42 : C------------------------------------- 43 : C debit air 44 : QA=QAE/E1/B1 45 : 46 : C debit liquide 47 : QE1=QEE1/E1/B1 48 : 49 : C debit vapeur 50 : QV1=QE1*Y1 51 : 52 : C debit liquide 53 : QL1= (1-Y1)*QE1 54 : 55 : C debit total 56 : Q = QA + QE1 57 : 58 : C------------------------------------- 59 : C DEFINITION DES MASSES VOLUMIQUES 60 : C CONSTITUANT: 61 : C : air 62 : C : vapeur 63 : C : eau 64 : C : eau + vapeur 65 : C : air + eau + vapeur 66 : C-------------------------------------- 67 : C PS1=PHI1*P1 68 : 69 : C masse volumique de air 70 : ROA1=(P1-PS1)/RA/T1 71 : C masse volumique vapeur 72 : ROV1=ROVAP0(PS1,T1) 73 : 74 : C masse volumique eau liquide 75 : ROL1=ROLIQ0(P1,T1) 76 : 77 : C calcul la fraction volumique de gaz 78 : AL =BALPHA(Y1,ROL1,ROV1) 79 : 80 : C calcul la masse volumique air + eau + vapeur 81 : RO1=AL *(ROA1+ROV1)+(1.-AL )*ROL1 82 : C------------------------------------------------- 83 : C CALCUL DE LA VITESSE DE ECOULEMENT EN ENTREE 84 : C------------------------------------------------- 85 : 86 : U1=QA/AL /ROA1 87 : 88 : C------------------------------------------------- 89 : C CALCUL DU COEFFICIENT DE COMPRESSIBILITE Z 90 : C CONSTITUANT: 91 : C : vapeur 92 : C------------------------------------------------- 93 : ZV = ZVAP0(ROV1,T1) 94 : 95 : C------------------------------------------------- 96 : C CALCUL DES VISCOSITES DYNAMIQUES 97 : C CONSTITUANT: 98 : C : air + vapeur 99 : C : air + eau + vapeur 100 : C------------------------------------------------- 101 : 102 : C viscosite dynamique air+vapeur 103 : XMU=BMUG(T1,PS1/P1,T0) 104 : 105 : C calcul de la viscosite dynamique air+vapeur+eau_liquide 106 : XMU=AL*XMU+BMUL(T1,T0)*(1.-AL) 107 : C WRITE(6,*) 'XMU',XMU 108 : 109 : C-------------------------------------------------- 110 : C CALCUL DES ENTHALPIES SPECIFIQUES 111 : C CONSTITUANT: 112 : C : eau liquide 113 : C : vapeur 114 : C-------------------------------------------------------- 115 : 116 : C calcul enthalpie specifique eau liquide 117 : HLS01 = HLS0(P1,T1) 118 : C WRITE(6,*) 'HL',HLS01 119 : C calcul enthalpie specifique de la vapeur 120 : HVS01 = HVS0(PS1,T1) 121 : C WRITE(6,*) 'HV',HVS01 122 : 123 : C-------------------------------------------------------- 124 : C CALCUL CHALEUR LATENTE DU MELANGE 125 : C CONSTITUANT: 126 : C : eau+vapeur 127 : C-------------------------------------------------------- 128 : 129 : C calcul la chaleur latente du melange 130 : XLAT = (HVS01-HLS01) 131 : C WRITE(6,*) 'XLAT',XLAT 132 : C-------------------------------------------------------- 133 : C CALCUL DES CONDUCTIVITES THERMIQUES 134 : C CONSTITUANT: 135 : C : air 136 : C-------------------------------------------------------- 137 : 138 : XLA = BLA(T1,T0) 139 : 140 : C-------------------------------------------------------- 141 : C CALCUL DES CHALEURS SPECIFIQUES 142 : C CONSTITUANT: 143 : C : vapeur 144 : C : eau liquide 145 : C : air + vapeur + eau 146 : C-------------------------------------------------------- 147 : 148 : C calcul chaleur specifique de la vapeur 149 : CV1 = DHVDT0(PS1,T1) 150 : 151 : C calcul chaleur specifique eau liquide 152 : CL1 = CPLIQ0(P1,HLS01) 153 : 154 : 155 : C calcul chaleur specifique du melange air vapeur liquide 156 : CP1 = (QA*CA+QE1*(Y1*CV1+(1-Y1)*CL1))/(QA+QE1) 157 : 158 : 159 : C WRITE(6,*) 'cv1,cl1,cp1',CV1,CL1,CP1 160 : C-------------------------------------------------------- 161 : C ACTUALISATION DES VARIABLES 162 : C CONSTITUANT: 163 : C : vitesse U 164 : C : titre Y 165 : C : densite du melange RO 166 : C : pression vapeur PS 167 : C : pression totale P 168 : C : temperature T 169 : C : debit eau (vapeu+liquide) QE 170 : C : debit total QX 171 : C------------------------------------------------------ 172 : 173 : E=E1 174 : B=B1 175 : U=U1 176 : Y=Y1 177 : RO=RO1 178 : PS=PS1 179 : P=P1 180 : T=T1 181 : QE=QE1 182 : QX=QA+QE 183 : 184 : 185 : 186 : IF ((X+DX).GT.1.D0) DX=1.D0-X 187 : 188 : X=X+DX 189 : DU=0 190 : T2=0 191 : P2=0 192 : 192 : CALL BTPAR(XN,EN,X,E2,NPP,ITP) 192 : CALL BTPAR(XN,BN,X,B2,NPP,ITP) 193 : 194 : NITER = 0 195 : 10 CONTINUE 196 : 197 : NITER = NITER + 1 198 : C write(6,*) 199 : C write(6,*) ' niter P2 T2 ',NITER,P2,T2 200 : TT2=T2 201 : PP2=P2 202 : YY2=Y2 203 : 204 : 205 : C-------------------------------------------------------- 206 : C CALCUL DU NOMBRE DE REYNOLDS DE ECOULEMENT 207 : C CONSTITUANT: 208 : C : air + eau + vapeur 209 : C-------------------------------------------------------- 210 : 211 : C calcul du nombre de REYNOLDS 212 : RE=QX*2*E/XMU 213 : C-------------------------------------------------------- 214 : C CALCUL DE LA PERTE CHARGE DU AU FORCE DE FROTTEMENT 215 : C CONSTITUANT: 216 : C : air + eau + vapeur 217 : C-------------------------------------------------------- 218 : 219 : DPP = XL*QX*QX /4/E 220 : 221 : C calcul de la perte de charge DP 222 : BK=BKFRO(RE,REL,XKL,XKT,2*E,RUG,RECU,XKUL,XKUT1,XKUT2,XKUT3, 223 : * XKUT4) 224 : DP=DX*BK*DPP/RO 225 : IF(KIMP.GE.1) THEN 226 : WRITE(6,*) 'bcond2 ap BKFRO : RE,QX,DPP= ',RE,QX,DPP 227 : WRITE(6,*) 'bcond2 ap BKFRO : DX,BK,RO= ',DX,BK,RO 228 : ENDIF 229 : P2=P1-DP 230 : IF (P2.LT.PSLIM) THEN 231 : PSQ = -1. 232 : RETURN 233 : ENDIF 234 : 235 : 236 : C-------------------------------------------------------- 237 : C CALCUL DU NOMBRE DE PRANDT DE ECOULEMENT 238 : C CONSTITUANT: 239 : C : air + eau + vapeur 240 : C-------------------------------------------------------- 241 : 242 : C calcul du nombre de Prandt estime en amont 243 : PR = XMU*CP1/BLA(T1,T0) 244 : 245 : 246 : 247 : C-------------------------------------------------------- 248 : C CALCUL DU COEFFICIENT ECHANGE THERMIQUE 249 : C CONSTITUANT: 250 : C : air + eau + vapeur 251 : C-------------------------------------------------------- 252 : 253 : 254 : C calcul du coefficient echange thermique 255 : HM = BHECH(RE,REL,PR,XLA,E) 256 : H = HM 257 : HW=0.D0 258 : 259 : C si condensation en film 260 : IF(T1.GE.TP) THEN 261 : HW=XW*BWALL(XLAT,CV1,PS/P,H,(DPSAT0(T1)/PSATT0(T1))) 262 : H = HM+HW 263 : ENDIF 264 : 265 : 266 : C calcul de la fraction molaire de vapeur 267 : PHI=PS/P 268 : 269 : C-------------------------------------------------------- 270 : C CALCUL DES COEFFICIENTS ALPHAi 271 : C CONSTITUANT: 272 : C : on se sert des fonctions de VARI 273 : C-------------------------------------------------------- 274 : 275 : 276 : C definition des alphai pour resolution systeme 277 : 278 : alpha1=DHVDT0(PS ,T ) 279 : alpha2=DHVDP0(PS ,T ) 280 : alpha3=DHLDP0(P ,T ) 281 : HLIQ=HLS0(P ,T ) 282 : alpha4=CPLIQ0(P ,HLIQ ) 283 : alpha5=DPSAT0(T ) 284 : alpha6=DZVDP0(PS ,T ) 285 : alpha7=DZVDT0(PS ,T ) 286 : 287 : 288 : C-------------------------------------------------------- 289 : C INTEGRATION DES PARAMETRES DE EQUATION THERMIQUE 290 : C Qy dT/dx= (H)(Tw-T)+S resolution BT3 291 : C-------------------------------------------------------- 292 : 293 : 294 : C definition de Qc 295 : Qc=QAE*CA+Y *QEE1*(alpha1+(alpha2*alpha5))+(1-Y )*alpha4*QEE1 296 : C WRITE(6,*) 'QA QE1 Qc ',QA,QE1,Qc 297 : 298 : 299 : C definition de a 300 : a=(RA/RV)*QAE 301 : C write(6,*) 'a ',a 302 : 303 : 304 : C DEFINITION DE ALPHA9 305 : x1=(PHI *alpha5)/(ZV*PSATT0(T )*(1-PHI )*(1-PHI )) 306 : x2=(PHI *(alpha5*alpha6+alpha7))/(1-PHI ) 307 : x2=x2/ZV/ZV 308 : alpha9=x1-x2 309 : 310 : C WRITE(6,*) 'ALPHA9 x1 x2 ',alpha9 ,x1,x2 311 : 312 : C DEFINITION DE ALPHA10 313 : x3=(PHI *XLAT*a)/(ZV*P *(1-PHI )*(1-PHI )) 314 : x4=alpha3*(1-Y )*QEE1 315 : alpha10=x3-x4 316 : C WRITE(6,*) 'alpha10 x3 x4 ',alpha10,x3,x4 317 : 318 : 319 : C calcul du terme source de la thermique S 320 : C rappel DP=P1-P2 et non l'inverse 321 : S=-(alpha10)*(DP/(XL*DX)) 322 : 323 : Qy=Qc+a*XLAT*alpha9 324 : 325 : 326 : C---------------------------------------------------- 327 : C Calcul la temperature T2 de sortie du troncon 328 : C---------------------------------------------------- 329 : 330 : HEFF = 2 * H 331 : C bt3 traite l echange de chaleur par unite d etendue 332 : Qy = Qy/B 333 : S = S/B 334 : T2 = BT3(T1,TP,S,HEFF,Qy,DX,XL) 335 : 336 : IF(KIMP.GE.1) THEN 337 : WRITE(*,*) 'bcond2 ap BT3: niter Qy B HEFF S',niter,Qy,B,HEFF,S 338 : WRITE(6,*) 'bcond2 ap BT3: T1 T2 TP= ',T1,T2,TP 339 : ENDIF 340 : 341 : C calcul de increment de temperature DT 342 : DT = T2-T1 343 : 344 : C---------------------------------------------------- 345 : C FIN DU TRONCON CALCULS DES NOUVELLES VALEURS 346 : C---------------------------------------------------- 347 : 348 : C calcul de nouvelle pression vapeur saturante 349 : PS2=PSATT0(T2) 350 : 351 : C-------------------------------------------------------- 352 : C DEFINITION DES NOUVELLES MASSES VOLUMIQUES 353 : C CONSTITUANT: 354 : C : air 355 : C : vapeur 356 : C : eau 357 : C : eau + vapeur 358 : C : air + eau + vapeur 359 : C-------------------------------------------------------- 360 : 361 : C masse volumique de air 362 : ROA2=(P2-PS2)/RA/T2 363 : 364 : C masse volumique vapeur 365 : ROV2=ROVAP0(PS2,T2) 366 : 367 : C-------------------------------------------------------- 368 : C CALCUL DU COEFFICIENT DE COMPRESSIBILITE Z de vapeur 369 : C CONSTITUANT: 370 : C : vapeur 371 : C-------------------------------------------------------- 372 : Z2 = ZVAP0(ROV2,T2) 373 : 374 : C masse volumique eau liquide 375 : ROL2=ROLIQ0(P2,T2) 376 : 377 : 378 : C-------------------------------------------------------- 379 : C Calcul du dQe 380 : C-------------------------------------------------------- 381 : 382 : C calcul increment de debit eau 383 : DQEE= -2.*HW*DX*XL*B*(T-TP)/XLAT 384 : 385 : C calcul nouveau debit eau (liq+vap) 386 : QEE2=QEE1+DQEE 387 : 388 : C calcul debit eau 389 : QE2=QEE2/E2/B2 390 : 391 : C calcul du titre de vapeur sortie troncon 392 : Y2=(ROV2/ROA2)*(QA/QE2) 393 : 394 : IF(KIMP.GE.3) THEN 395 : WRITE(6,*) ' ROA2 ROV2 Y2',ROA2,ROV2,Y2 396 : write(6,*) ' QE1 QE2 ',QE1,QE2 397 : ENDIF 398 : 399 : 400 : C-------------------------------------------------------- 401 : C On raffine le maillage dans 2 cas 402 : C-------------------------------------------------------- 403 : 404 : IF ( ((QEE2/QEE1).LT.1.E-5).OR.(Y2.GE.1.)) THEN 405 : C WRITE(6,*) ' WALL FLUX ',QEE2/QEE1 406 : RINDEX=0.5 407 : RETURN 408 : ENDIF 409 : 410 : 411 : QV2 = QE2 * Y2 412 : QL2 = QE2 * (1.-Y2) 413 : 414 : 415 : C calcul la fraction volumique de gaz 416 : AL2=BALPHA(Y2,ROL2,ROV2) 417 : 418 : 419 : C calcul la masse volumique air + eau + vapeur 420 : RO2=AL2*(ROA2+ROV2)+(1.-AL2)*ROL2 421 : 422 : 423 : 424 : 425 : 426 : C-------------------------------------------------------- 427 : C CALCUL DE LA NOUVELLE VITESSE DE ECOULEMENT 428 : C-------------------------------------------------------- 429 : 430 : U2= QA/AL2/ROA2 431 : 432 : C-------------------------------------------------------- 433 : C DEFINITION DES NOUVELLES VALEURS MOYENNES 434 : C 435 : C : Pression 436 : C : Temperature 437 : C : pression vapeur 438 : C-------------------------------------------------------- 439 : 440 : P=(P1+P2)/2 441 : T=(T1+T2)/2 442 : PS=(PS1+PS2)/2 443 : 444 : PHI = PS/P 445 : 446 : C-------------------------------------------------------- 447 : C CALCUL DES VITESSES MOYENNES SUR TRONCON 448 : C CONSTITUANT: melange 449 : C : masse volumique 450 : C : vitesse 451 : C : debit eau 452 : C : debit total 453 : C : titre 454 : C-------------------------------------------------------- 455 : 456 : 457 : C calcul de la masse volumique moyenne melange 458 : RO=(RO1+RO2)/2 459 : 460 : C calcul de la vitesse moyenne du melange 461 : U=Q/RO 462 : 463 : 464 : C calcul du debit eau moyen du melange 465 : QE=(QE1+QE2)/2. 466 : 467 : C calcul du debit total moyen du melange 468 : QX=QA+QE 469 : 470 : C calcul du titre moyen sur le troncon 471 : Y=(Y1+Y2)/2 472 : 473 : C----------------------------------------------------------------- 474 : C CALCUL DU TAUX DE CONDENSATION EN FILM 475 : C----------------------------------------------------------------- 476 : 477 : QW2=(QEE2-QEE1)/DX/XL/B 478 : 479 : 480 : 481 : IF(KIMP.GE.3) write(6,*) ' P2 T2 Y2 ',P2,T2,Y2 482 : 483 : C definition des conditions fin de boucle 484 : 485 : 486 : C critere sur P,T,Y 487 : 488 : ERT=DABS((T2-TT2)/T2) 489 : ERP=DABS((P2-PP2)/P2) 490 : ERY=DABS((Y2-YY2)/Y2) 491 : 492 : IF (((ERT.GT.1E-4).OR.(ERP.GT.1E-4).OR.(ERY.GT.1E-4)) 493 : $ .AND.(NITER.LE.10)) GOTO 10 494 : 495 : IF(KIMP.GE.2) THEN 496 : write(6,2300) X, P2/P0, T2-T0, Y2, RO2, QE1, QE2 497 : 2300 FORMAT((1X,'bcond2: X P2 T2 Y2 RO2 QE1 QE2 ',7E12.5)) 498 : ENDIF 499 : IF(KIMP.GE.3) THEN 500 : WRITE(6,*) ' HM HW HW/HM ',HM,HW,HW/HM 501 : ENDIF 502 : PHI2=PHI1 503 : 504 : 505 : 506 : C-------------------------------------------------------- 507 : C BILAN ENERGIE SUR TRONCON 508 : C-------------------------------------------------------- 509 : 510 : C write(6,*) ' E1 E2 ',E1,E2 511 : 512 : HH1=(QA*CA*T1)+(QV1*HVS0(PS1,T1))+(QL1*HLS0(P1,T1)) 513 : HH1=HH1*E1*B1 514 : HH2=(QA*CA*T2)+(QV2*HVS0(PS2,T2))+(QL2*HLS0(P2,T2)) 515 : HH2=HH2*E2*B2 516 : C puissance fluide 517 : C ! PF = HH2-HH1 518 : PF = (HH2-HH1) - (DQEE*HVS0(PS1,T1)) 519 : 520 : C puissance paroi 521 : PP = H *(T-TP)*2*DX*XL*B 522 : 523 : 524 : IF(KIMP.GE.2) THEN 525 : write(6,2200) X,HH1,HH2,PF,PP 526 : 2200 format(1X,'bcond2 X HH1 HH2 DH HDT ',5E12.5) 527 : ENDIF 528 : IF(KIMP.NE.0.AND.X.EQ.1.) THEN 529 : write(6,2110) X,RE,BK 530 : 2110 format(1X,'bcond2 X RE BKRO ',3E12.5) 531 : ENDIF 532 : DPF=PF/DX/XL/B 533 : DPP=PP/DX/XL/B 534 : 535 : RETURN 536 : END 537 : 538 : 539 : 540 : 541 : 542 : 543 :
© Cast3M 2003 - All rights reserved.
Disclaimer