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