!-------------------------------------------------------------------------------
! Code_Saturne version 4.0.6
! --------------------------
! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2016 EDF S.A.
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 2 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!-------------------------------------------------------------------------------
!===============================================================================
! Purpose:
! -------
!> \file cs_user_physical_properties.f90
!> \brief Definition of physical variable laws.
!>
!> usphyv
!> \brief Definition of physical variable laws.
!>
!> \section Warning
!>
!> It is \b forbidden to modify turbulent viscosity \c visct here
!> (a specific subroutine is dedicated to that: \ref usvist)
!>
!> - icp = 1 must have been specified
!> in \ref usipsu if we wish to define a variable specific heat
!> cpro_cp (otherwise: memory overwrite).
!>
!> - the kivisl field integer key (scalar_diffusivity_id)
!> must have been specified
!> in \ref usipsu if we wish to define a variable viscosity
!> \c viscls.
!>
!>
!> \remarks
!> - This routine is called at the beginning of each time step
!> Thus, AT THE FIRST TIME STEP (non-restart case), the only
!> values initialized before this call are those defined
!> - in the GUI or \ref usipsu (cs_user_parameters.f90)
!> - density (initialized at \c ro0)
!> - viscosity (initialized at \c viscl0)
!> - in the GUI or \ref cs_user_initialization
!> - calculation variables (initialized at 0 by defaut
!> or to the value given in the GUI or in \ref cs_user_initialization)
!>
!> - We may define here variation laws for cell properties, for:
!> - density: rom kg/m3
!> - density at boundary faces: romb kg/m3)
!> - molecular viscosity: cpro_viscl kg/(m s)
!> - specific heat: cpro_cp J/(kg degrees)
!> - diffusivities associated with scalars: cpro_vscalt kg/(m s)
!>
!> \b Warning: if the scalar is the temperature, cpro_vscalt corresponds
!> to its conductivity (Lambda) in W/(m K)
!>
!>
!> The types of boundary faces at the previous time step are available
!> (except at the first time step, where arrays \c itypfb and \c itrifb have
!> not been initialized yet)
!>
!> It is recommended to keep only the minimum necessary in this file
!> (i.e. remove all unused example code)
!>
!>
!> \section usphyv_cell_id Cells identification
!>
!> Cells may be identified using the \ref getcel subroutine.
!> The syntax of this subroutine is described in the
!> \ref cs_user_boundary_conditions subroutine,
!> but a more thorough description can be found in the user guide.
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
! mode name role !
!______________________________________________________________________________!
!> \param[in] nvar total number of variables
!> \param[in] nscal total number of scalars
!> \param[in] mbrom indicator of filling of romb array
!> \param[in] dt time step (per cell)
!_______________________________________________________________________________
subroutine usphyv &
( nvar , nscal , &
mbrom , &
dt )
!===============================================================================
!===============================================================================
! Module files
!===============================================================================
use paramx
use pointe
use numvar
use optcal
use cstphy
use entsor
use parall
use period
use ppppar
use ppthch
use ppincl
use field
use mesh
!===============================================================================
implicit none
! Arguments
integer nvar , nscal
integer mbrom
double precision dt(ncelet)
! Local variables
integer ivart, iel, ifac
integer ith, iscal, ii, ifcvsl
integer nlfac1, nlfac2, n, jj, j, ivar, i, ilelt, Hnum
integer, allocatable, dimension(:) :: lstfac1, lstfac2
!double precision vara, varb, varc, varam, varbm, varcm, vardm
!double precision varal, varbl, varcl, vardl
!double precision varac, varbc
!double precision xvart
double precision theta, thcond, temp, riel_det_out, riel_det_inn, ratio
double precision loc_Rout, loc_Rinn, en, elt0, elt, dens, del, angiel
double precision, dimension(:), pointer :: bfpro_rom, cpro_rom
double precision, dimension(:), pointer :: cpro_viscl, cpro_vscalt, cpro_cp, cpro_beta
double precision, dimension(:), pointer :: cvar_scalt
double precision, dimension(:), pointer :: cpro_crom
double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp
double precision, dimension(:,:), pointer :: cvar_vel
double precision relax, pi, PT, x1,x2,x3,x4, a1,a2,a3,a4, c1,c2,c3,c4
double precision ngtan, temp1, temp2, xm, ym, zm, taniel
double precision, allocatable, dimension(:) :: splth, splttan, R1D, vertical_z,treco, sch_ang
double precision, allocatable, dimension(:,:) :: Rdet_inn, Rdet_out, tbulk, schR_inn, schR_out
double precision, allocatable, dimension(:,:) :: tfac_inn, flux_inn, tfac_out, flux_out
!allocate arrays
allocate(splth(1:6),splttan(1:5))
allocate(R1D(1:108),vertical_z(1:161))
allocate(Rdet_inn(1:9,1:5),Rdet_out(1:9,1:5))
allocate(tbulk(1:4,1:160),tfac_inn(1:4,1:160),flux_inn(1:4,1:160))
allocate(schR_inn(1:4,1:160),schR_out(1:4,1:160))
allocate(tfac_out(1:4,1:160),flux_out(1:4,1:160))
allocate(treco(nfabor))
allocate(sch_ang(9))
relax=0.2
pi=3.141592653
!read R-value data from result of detailed model Rdet_inn(9,5), Rdet_out(9,5)
OPEN(3,FILE='R-values-from-detmod.txt',status='old')
do i=1,108
read(3,*)R1D(i)
enddo
close(3)
n=0
do j=1,5
do i=1,9
n=n+1
Rdet_inn(i,j)=R1D(n)
Rdet_out(i,j)=R1D(n+54)
enddo
enddo
!read vertical corrdinates of mesh nodes
OPEN(4,FILE='vertical_z.txt',status='old')
do i=1,161
read(4,*)vertical_z(i)
enddo
close(4)
!split height H1-5
splth(1)=0.d0
splth(2)=1.092
splth(3)=3.198
splth(4)=5.304
splth(5)=7.228
splth(6)=8.33
!split angle (4 cell at narrow gap of orig cfd model)
splttan(1)=1.0241399
splttan(2)=1.53282695
splttan(3)=2.472112956
splttan(4)=5.2104926632
splttan(5)=10**5.d0
PT=40
x1=1.79297
x2=406.813
x3=0.839844
x4=40160
a1=0.42163
a2=1.17385*(10**(-3.d0))
a3=-1.77174*(10**(-6.d0))
a4=1.29644*(10**(-9.d0))
c1=1.5896
c2=-0.22458
c3=38.0755*(10**(-3.d0))
c4=-2.5235*(10**(-3.d0))
ngtan=1.0241399
sch_ang=(/ 95, 105, 115, 125, 135, 145, 155, 165, 175 /)
!calculate schRth(4,160,2), sub-channel thermal resistance
iscal = iscalt
ivar = isca(iscal)
call field_get_val_s(icrom, cpro_crom)
call field_get_val_s(ivarfl(isca(iscalt)), cvar_scalt)
call field_get_val_v(ivarfl(iu), cvar_vel)
call field_get_coefa_s(ivarfl(ivar), coefap)
call field_get_coefb_s(ivarfl(ivar), coefbp) !for face temperature
call field_get_coefaf_s(ivarfl(ivar), cofafp)
call field_get_coefbf_s(ivarfl(ivar), cofbfp) !for face heat flux
do ifac = 1, nfabor
iel = ifabor(ifac)
treco(ifac) = cvar_scalt(iel)
enddo
allocate(lstfac1(max(ncel,nfac,nfabor)))
call getfbr('annulus_inn_face', nlfac1, lstfac1)
allocate(lstfac2(max(ncel,nfac,nfabor)))
call getfbr('annulus_out_face', nlfac2, lstfac2)
do i=1,4
do j=1,160
!calculate tbulk
temp1=0
temp2=0
do iel=1,ncel
xm=xyzcen(1,iel)
ym=xyzcen(2,iel)
zm=xyzcen(3,iel)
taniel=ym/xm
if(taniel .gt. splttan(i) .and. taniel .lt. splttan(i+1)) then
if(zm .gt. vertical_z(j) .and. zm .lt. vertical_z(j+1)) then
temp1=temp1+cvar_vel(3,iel)*cpro_crom(iel)*volume(iel)*cvar_scalt(iel)
temp2=temp2+cvar_vel(3,iel)*cpro_crom(iel)*volume(iel)
endif
endif
enddo
call parsom(temp1)
call parsom(temp2)
tbulk(i,j)=temp1/temp2
!calculate Tinn and Flux_inn
do ilelt=1, nlfac1
ifac = lstfac1(ilelt)
iel = ifabor(ifac)
xm=xyzcen(1,iel)
ym=xyzcen(2,iel)
zm=xyzcen(3,iel)
taniel=ym/xm
if(taniel .gt. splttan(i) .and. taniel .lt. splttan(i+1)) then
if(zm .gt. vertical_z(j) .and. zm .lt. vertical_z(j+1)) then
tfac_inn(i,j) = coefap(ifac)+treco(ifac)*coefbp(ifac) !face temperature
flux_inn(i,j) = abs(cofafp(ifac) + cofbfp(ifac)*treco(ifac)) !face heat flux
endif
endif
enddo
!calculate Tout and Flux_out
do ilelt=1, nlfac2
ifac = lstfac2(ilelt)
iel = ifabor(ifac)
xm=xyzcen(1,iel)
ym=xyzcen(2,iel)
zm=xyzcen(3,iel)
taniel=ym/xm
if(taniel .gt. splttan(i) .and. taniel .lt. splttan(i+1)) then
if(zm .gt. vertical_z(j) .and. zm .lt. vertical_z(j+1)) then
tfac_out(i,j) = coefap(ifac)+treco(ifac)*coefbp(ifac) !face temperature
flux_out(i,j) = abs(cofafp(ifac) + cofbfp(ifac)*treco(ifac)) !face heat flux
endif
endif
enddo
enddo
enddo
do i=1,4
do j=1,160
schR_inn(i,j)=abs(flux_inn(i,j)/(tfac_inn(i,j)-tbulk(i,j)))
schR_out(i,j)=abs(flux_out(i,j)/(tfac_out(i,j)-tbulk(i,j)))
enddo
enddo
do iel=1,ncel
thcond=0.5
cpro_vscalt(iel)=thcond
enddo
!do iel=1,ncel
! xm=xyzcen(1,iel)
!enddo
deallocate(splth,splttan)
deallocate(R1D,vertical_z)
deallocate(Rdet_inn,Rdet_out)
deallocate(tbulk,tfac_inn,flux_inn)
deallocate(tfac_out,flux_out)
deallocate(treco)
deallocate(sch_ang)
return
end subroutine usphyv