1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23!=============================================================================== 24! Function: 25! --------- 26 27!> \file cs_gas_mix_physical_properties.f90 28!> 29!> \brief This subroutine fills physical properties which are variable in time 30!> for the gas mixtures modelling with or without steam inside the fluid 31!> domain. In presence of steam, this one is deduced from the 32!> noncondensable gases transported as scalars 33!> (by means of the mass fraction of each species). 34!> 35!------------------------------------------------------------------------------- 36 37!------------------------------------------------------------------------------- 38! Arguments 39!______________________________________________________________________________. 40! mode name role ! 41!______________________________________________________________________________! 42!_______________________________________________________________________________ 43 44 45subroutine cs_gas_mix_physical_properties 46 47!=============================================================================== 48 49!=============================================================================== 50! Module files 51!=============================================================================== 52 53use paramx 54use numvar 55use optcal 56use cstphy 57use cstnum 58use entsor 59use pointe 60use albase 61use parall 62use period 63use ppppar 64use ppthch 65use ppincl 66use mesh 67use field 68use cavitation 69use cs_c_bindings 70use cs_f_interfaces 71 72!=============================================================================== 73 74implicit none 75 76! Arguments 77 78! Local variables 79 80integer iel , iscal, ifcvsl, iesp, jesp, ierror, f_id 81 82character(len=80) :: name_i, name_j, name_d 83 84double precision xsum_mu, xsum_lambda, phi_mu, phi_lambda, x_k 85double precision mu_i, mu_j, lambda_i, lambda_j 86 87type(gas_mix_species_prop), pointer :: s_j, s_i 88type(gas_mix_species_prop), target :: s_d 89type(gas_mix_species_prop), dimension(:), allocatable, target :: s_k 90 91double precision, allocatable, dimension(:), target :: lam_loc 92double precision, allocatable, dimension(:), target :: tk_loc 93 94double precision, dimension(:), pointer :: cpro_rho 95double precision, dimension(:), pointer :: cpro_viscl, cpro_cp 96double precision, dimension(:), pointer :: cpro_venth, cpro_vyk 97double precision, dimension(:), pointer :: cvar_enth , cvar_yk, tempk 98double precision, dimension(:), pointer :: cvar_yi, cvar_yj 99double precision, dimension(:), pointer :: y_d, ya_d 100double precision, dimension(:), pointer :: mix_mol_mas 101double precision, dimension(:), pointer :: lambda 102 103!=============================================================================== 104 105!=============================================================================== 106! 1. Initializations 107!=============================================================================== 108 109ierror = 0 110if (irovar.le.0) ierror = 1 111if (ivivar.le.0) ierror = 1 112if (icp.lt.0) ierror = 1 113 114if (ierror.gt.0) then 115 call csexit(1) 116endif 117 118! In compressible, the density is updated after the pressure step (cfmspr) 119if (ippmod(icompf).lt.0) then 120 call field_get_val_s(ivarfl(isca(iscalt)), cvar_enth) 121 ! Density value 122 call field_get_val_s(icrom, cpro_rho) 123 allocate(tk_loc(ncel)) 124 tempk => tk_loc 125else 126 call field_get_val_s(ivarfl(isca(itempk)), tempk) 127endif 128 129! Molecular dynamic viscosity value 130call field_get_val_s(iviscl, cpro_viscl) 131! Specific heat value 132if (icp.ge.0) then 133 call field_get_val_s(icp, cpro_cp) 134 135! Stop if Cp is not variable 136else 137 write(nfecra,1000) icp 138 call csexit (1) 139endif 140 141! Lambda/Cp value 142call field_get_key_int(ivarfl(isca(iscalt)), kivisl, ifcvsl) 143call field_get_val_s(ifcvsl, cpro_venth) 144 145call field_get_val_s(igmxml, mix_mol_mas) 146 147! Deduce mass fraction (y_d) which is 148! y_h2o_g in presence of steam or 149! y_he/y_h2 with noncondensable gases 150if (ippmod(igmix).eq.0) then 151 name_d = "y_he" 152elseif (ippmod(igmix).eq.1) then 153 name_d = "y_h2" 154elseif (ippmod(igmix).ge.2.and.ippmod(igmix).lt.5) then 155 name_d = "y_h2o_g" 156else ! ippmod(igmix).eq.5 157 name_d = "y_o2" 158endif 159call field_get_val_s_by_name(name_d, y_d) 160call field_get_val_prev_s_by_name(name_d, ya_d) 161call field_get_id(name_d, f_id) 162call field_get_key_struct_gas_mix_species_prop(f_id, s_d) 163 164if (ippmod(icompf).lt.0) then 165 allocate(lam_loc(ncelet)) 166 lambda => lam_loc 167else 168 call field_get_key_int(ivarfl(isca(itempk)), kivisl, ifcvsl) 169 call field_get_val_s(ifcvsl, lambda) 170endif 171 172allocate(s_k(nscasp+1)) 173 174!=============================================================================== 175!2. Define the physical properties for the gas mixture with: 176! - the density (rho_m) and specific heat (cp_m) of the gas mixture function 177! temperature and species scalar (yk), 178! - the dynamic viscosity (mu_m) and conductivity coefficient (lbd_m) of 179! the gas mixture function ot the enthalpy and species scalars, 180! - the diffusivity coefficients of the scalars (Dk, D_enh). 181!=============================================================================== 182 183!Storage the previous value of the deduced mass fraction ya_d 184do iel=1, ncelet 185 ya_d(iel) = y_d(iel) 186enddo 187 188!----------------------------------------- 189! Compute the mass fraction (y_d) deduced 190! from the mass fraction (yk) transported 191!----------------------------------------- 192 193! Initialization 194do iel = 1, ncel 195 y_d(iel) = 1.d0 196 197 ! Mixture specific heat 198 cpro_cp(iel) = 0.d0 199 200 ! Mixture molar mass 201 mix_mol_mas(iel) = 0.d0 202 203 ! Mixture molecular diffusivity 204 cpro_viscl(iel) = 0.d0 205 206 ! Thermal conductivity 207 lambda(iel) = 0.d0 208enddo 209 210do iesp = 1, nscasp 211 ! Mass fraction array of the different species 212 call field_get_val_s(ivarfl(isca(iscasp(iesp))), cvar_yk) 213 214 call field_get_key_struct_gas_mix_species_prop( & 215 ivarfl(isca(iscasp(iesp))), s_k(iesp)) 216 217 do iel = 1, ncel 218 y_d(iel) = y_d(iel)-cvar_yk(iel) 219 mix_mol_mas(iel) = mix_mol_mas(iel) + cvar_yk(iel)/s_k(iesp)%mol_mas 220 enddo 221enddo 222 223! Clipping 224do iel = 1, ncel 225 y_d(iel) = max(y_d(iel), 0.d0) 226enddo 227 228!Finalize the computation of the Mixture molar mass 229s_k(nscasp+1) = s_d 230do iel = 1, ncel 231 mix_mol_mas(iel) = mix_mol_mas(iel) + y_d(iel)/s_d%mol_mas 232 mix_mol_mas(iel) = 1.d0/mix_mol_mas(iel) 233enddo 234 235!============================================================== 236! Mixture specific heat function of species specific heat (cpk) 237! and mass fraction of each gas species (yk), as below: 238! ----------------------------- 239! - noncondensable gases and the mass fraction deduced: 240! cp_m(iel) = Sum( yk.cpk)_k[0, nscasp] 241! + y_d.cp_d 242! ----------------------------- 243! remark: the mass fraction deduced depending of the 244! modelling chosen by the user. 245! with: 246! - igmix = 0 or 1, a noncondensable gas 247! - igmix > 2 , a condensable gas (steam) 248!============================================================== 249do iesp = 1, nscasp 250 call field_get_val_s(ivarfl(isca(iscasp(iesp))), cvar_yk) 251 252 do iel = 1, ncel 253 cpro_cp(iel) = cpro_cp(iel) + cvar_yk(iel)*s_k(iesp)%cp 254 enddo 255enddo 256 257! Finalization 258do iel = 1, ncel 259 cpro_cp(iel) = cpro_cp(iel) + y_d(iel)*s_d%cp 260enddo 261 262!=========================================================== 263! gas mixture density function of the temperature, pressure 264! and the species scalars with taking into account the 265! dilatable effects, as below: 266! ---------------------- 267! - with inlet/outlet conditions: 268! [idilat=2]: rho= p0/(R. temp(1/sum[y_i/M_i])) 269! - with only wall conditions: 270! [idilat=3]: rho= pther/(R. temp(1/sum[y_i/M_i])) 271! ---------------------- 272! i ={1, .. ,N} : species scalar number 273! y_i : mass fraction of each species 274! M_i : molar fraction [kg/mole] 275! R : ideal gas constant [J/mole/K] 276! p0 : atmos. pressure (Pa) 277! pther : pressure (Pa) integrated on the 278! fluid domain 279!=========================================================== 280 281if (ippmod(icompf).lt.0) then 282 do iel = 1, ncel 283 ! Evaluate the temperature thanks to the enthalpy 284 tempk(iel) = cvar_enth(iel)/ cpro_cp(iel) 285 if (idilat.eq.3) then 286 cpro_rho(iel) = pther*mix_mol_mas(iel)/(cs_physical_constants_r*tempk(iel)) 287 else 288 cpro_rho(iel) = p0*mix_mol_mas(iel)/(cs_physical_constants_r*tempk(iel)) 289 endif 290 enddo 291endif 292 293!================================================== 294! Dynamic viscosity and conductivity coefficient 295! the physical properties associated to the gas 296! mixture with or without condensable gas. 297!================================================== 298 299! Loop over ALL the species 300do iesp = 1, nscasp+1 301 302 s_i => s_k(iesp) 303 304 ! Mass fraction deduced 305 ! (as steam or noncondensable gas) 306 if (iesp.eq.nscasp+1) then 307 cvar_yi => y_d 308 name_i = name_d 309 ! Noncondensable species 310 else 311 call field_get_val_s(ivarfl(isca(iscasp(iesp))), cvar_yi) 312 call field_get_name (ivarfl(isca(iscasp(iesp))), name_i) 313 endif 314 315 do iel = 1, ncel 316 317 ! Viscosity and conductivity laws 318 ! for each mass fraction species 319 320 if (ivsuth.eq.0) then 321 ! With a linear law 322 call cs_local_physical_properties & 323 !================================ 324 ( mu_i, lambda_i, tempk(iel), tkelvi, s_i, name_i) 325 else 326 ! Or : with a Sutherland law 327 call cs_local_physical_properties_suth & 328 !================================ 329 ( mu_i, lambda_i, tempk(iel), s_i, name_i) 330 endif 331 332 xsum_mu = 0.d0 333 xsum_lambda = 0.d0 334 335 ! Loop over ALL the species 336 do jesp = 1, nscasp+1 337 338 s_j => s_k(jesp) 339 340 ! Mass fraction deduced 341 ! (as steam or noncondensable gas) 342 if (jesp.eq.nscasp+1) then 343 cvar_yj => y_d 344 name_j = name_d 345 ! Noncondensable species 346 else 347 call field_get_val_s(ivarfl(isca(iscasp(jesp))), cvar_yj) 348 call field_get_name (ivarfl(isca(iscasp(jesp))), name_j) 349 endif 350 351 if (ivsuth.eq.0) then 352 ! With a linear law 353 call cs_local_physical_properties & 354 !================================ 355 ( mu_j, lambda_j, tempk(iel), tkelvi, s_j, name_j) 356 else 357 ! Or : with a Sutherland law 358 call cs_local_physical_properties_suth & 359 !================================ 360 ( mu_j, lambda_j, tempk(iel), s_j, name_j) 361 endif 362 363 phi_mu = (1.d0/sqrt(8.d0)) & 364 *(1.d0 + s_i%mol_mas / s_j%mol_mas)**(-0.5d0) & 365 *(1.d0 + (mu_i / mu_j )**(+0.5d0) & 366 * (s_j%mol_mas / s_i%mol_mas)**(+0.25d0))**2 367 368 369 phi_lambda = (1.d0/sqrt(8.d0)) & 370 *(1.d0 + s_i%mol_mas / s_j%mol_mas)**(-0.5d0) & 371 *(1.d0 + (lambda_i / lambda_j )**(+0.5d0) & 372 * (s_j%mol_mas / s_i%mol_mas)**(+0.25d0))**2 373 374 x_k = cvar_yj(iel)*mix_mol_mas(iel)/s_j%mol_mas 375 xsum_mu = xsum_mu + x_k * phi_mu 376 377 xsum_lambda = xsum_lambda + x_k * phi_lambda 378 379 enddo 380 381 ! Mixture viscosity defined as function of the scalars 382 !----------------------------------------------------- 383 x_k = cvar_yi(iel)*mix_mol_mas(iel)/s_i%mol_mas 384 cpro_viscl(iel) = cpro_viscl(iel) + x_k * mu_i / xsum_mu 385 386 lambda(iel) = lambda(iel) + x_k * lambda_i / xsum_lambda 387 388 enddo 389enddo 390 391!===================================================== 392! Dynamic viscosity and conductivity coefficient 393! the physical properties filled for the gas mixture 394!===================================================== 395! Same diffusivity for all the scalars except the enthalpy 396do iesp = 1, nscasp 397 398 iscal = iscasp(iesp) 399 400 call field_get_key_int (ivarfl(isca(iscal)), kivisl, ifcvsl) 401 call field_get_val_s(ifcvsl, cpro_vyk) 402 403 do iel = 1, ncel 404 cpro_vyk(iel)= cpro_viscl(iel) 405 enddo 406 407enddo 408 409if(ippmod(icompf).lt.0) then 410 ! --- Lambda/Cp of the thermal scalar 411 do iel = 1, ncel 412 cpro_venth(iel) = lambda(iel)/cpro_cp(iel) 413 enddo 414 415 ! deallocate local arrays if not compressible 416 deallocate(tk_loc) 417 deallocate(lam_loc) 418 tempk => null() 419 lambda => null() 420endif 421 422deallocate(s_k) 423 424!=============================================================================== 425! 3. Checking of the user values 426!=============================================================================== 427 428!-------- 429! Formats 430!-------- 431 432 1000 format( & 433'@',/, & 434'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 435'@',/, & 436'@ @@ WARNING: stop when computing physical quantities',/, & 437'@ =======',/, & 438'@ Inconsistent calculation data',/, & 439'@',/, & 440'@ usipsu specifies that the specific heat is uniform',/, & 441'@ icp = ',i10 ,' while',/, & 442'@ cs_user_physical_properties prescribes a variable specific heat.',/,& 443'@',/, & 444'@ The calculation will not be run.',/, & 445'@',/, & 446'@ Modify usipsu or cs_user_physical_properties.',/, & 447'@ ',/,& 448'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 449'@',/) 450 451!---- 452! End 453!---- 454 455return 456end subroutine cs_gas_mix_physical_properties 457 458!=============================================================================== 459! Purpose: 460! ------- 461 462!> \file cs_gas_mix_physical_properties.f90 463!> 464!> \brief This user subroutine is used to compute the dynamic viscosity and 465!> conductivity coefficient associated to each gas species. 466!------------------------------------------------------------------------------- 467 468!------------------------------------------------------------------------------- 469! Arguments 470!______________________________________________________________________________. 471! mode name role ! 472!______________________________________________________________________________! 473!> \param[out] mu dynamic viscosity associated to the gas species 474!> \param[out] lambda conductivity coefficient of the gas species 475!> \param[in] tk temperature variable in kelvin 476!> \param[in] tkelvin reference temperature value 477!> \param[in] spro constants used for the physcial laws 478!> \param[in] name name of the field associated to the gas species 479!_______________________________________________________________________________ 480 481subroutine cs_local_physical_properties(mu, lambda, tk, tkelvin, spro, name) 482!=============================================================================== 483 484use field 485use cs_c_bindings 486 487!=============================================================================== 488 489implicit none 490 491! Arguments 492 493double precision mu, lambda 494double precision tk, tkelvin 495 496character(len=80) :: name 497 498type(gas_mix_species_prop) spro 499 500!=============================================================================== 501! The viscosity law for each species is defined 502! as below: 503! ---------------------------------- 504! 1/. for Steam species: 505! mu = mu_a.(tk - tkelvi), with a law in (°C) unit 506! 2/. for Helium species: 507! mu = mu_a .(tk/tkelvi)**c + mu_b, with nodim. law 508! 3/. for Hydrogen species: 509! mu = mu_a .(tk-tkelvi) + mu_b, with t (°C) 510! 4/. for Oxygen and Nitrogen species: 511! mu = mu_a .(tk) + mu_b, with t in (°K) 512! 513! The conductivity expression for each species is 514! defined as: 515! ---------------------------------- 516! 1/. for Steam species: 517! lambda = lambda_a .(tk-tkelvi) + lambda_b, with a law in (°C) unit 518! 2/. for Helium species: 519! lambda = lambda_a .(tk/tkelvi)**c + lambda_b, with nodim. law 520! 3/. for Hydrogen species: 521! lambda = lambda_a .tk + lambda_b, with tk (°K) 522! 4/. for Oxygen and Nitrogen species : 523! lambda = lambda_a .(tk) + lambda_b, with t (°K) 524!=============================================================================== 525 526if (name.eq.'y_h2o_g') then 527 mu = spro%mu_a *(tk-tkelvin) + spro%mu_b 528 lambda = spro%lambda_a*(tk-tkelvin) + spro%lambda_b 529elseif(name.eq.'y_he') then 530 mu = spro%mu_a * (tk/tkelvin)**0.7d0 531 lambda = spro%lambda_a * (tk/tkelvin)**0.7d0 532elseif(name.eq.'y_h2') then 533 mu = spro%mu_a * (tk-tkelvin) + spro%mu_b 534 lambda = spro%lambda_a * tk + spro%lambda_b 535elseif (name.eq.'y_o2'.or.name.eq.'y_n2') then 536 mu = spro%mu_a * tk + spro%mu_b 537 lambda = spro%lambda_a * tk + spro%lambda_b 538else 539 call csexit(1) 540endif 541 542!---- 543! End 544!---- 545return 546 547end subroutine cs_local_physical_properties 548 549subroutine cs_local_physical_properties_suth(mu, lambda, tk,spro,name) 550!=============================================================================== 551use field 552use cs_c_bindings 553use cstphy 554use ppthch 555!=============================================================================== 556 557implicit none 558 559! Arguments 560 561double precision mu, lambda 562double precision tk 563 564character(len=80) :: name 565type(gas_mix_species_prop) spro 566 567! Local variables 568 569double precision muref, lamref 570double precision trefmu, treflam, smu, slam 571 572!=============================================================================== 573! Sutherland law for viscosity and thermal conductivity 574! The viscosity law for each specie is defined 575! as below: 576! ---------------------------------- 577! mu = muref*(T/Tref)**(3/2)*(Tref+S1)/(T+S1) 578 579! The conductivity expression for each specie is 580! defined as: 581! ---------------------------------- 582! lambda = lambdaref*(T/Tref)**(3/2)*(Tref+S2)/(T+S2) 583! ------------------------------------ 584! S1 and S2 are respectively Sutherland temperature for conductivity and 585! Sutherland temperature for viscosity of the considered specie 586! Tref is a reference temperature, equal to 273K for a perfect gas. 587! For steam (H20), Tref has not the same value in the two formulae. 588! Available species : O2, N2, H2, H20 and He 589! The values for the parameters come from F.M. White's book "Viscous Fluid Flow" 590!================================================================================ 591 592if (name.ne.'y_h2o_g' .and. name.ne.'y_he' .and. name.ne.'y_o2' & 593 .and. name.ne.'y_n2' .and. name.ne.'y_h2') then 594 call csexit(1) 595endif 596 597muref = spro%muref 598lamref = spro%lamref 599trefmu = spro%trefmu 600treflam = spro%treflam 601smu = spro%smu 602slam = spro%slam 603 604mu = muref * (tk / trefmu)**1.5d0 & 605 * ((trefmu+smu) / (tk+smu)) 606lambda = lamref * (tk / treflam)**1.5d0 & 607 * ((treflam+slam) / (tk+slam)) 608 609!---- 610! End 611!---- 612return 613 614end subroutine cs_local_physical_properties_suth 615