Data setting for drift scalars

Introduction

This page provides an example of code blocks that may be used to perform a calculation with drift scalars.

Physical properties

Local variables to be added

The following local variables need to be defined for the examples in this section:

integer          ivart, iel, ifac
integer          iscal, iflid, iscdri, ifcvsl
integer          f_id, keydri, nfld, keysca
double precision rho, viscl
double precision diamp, rhop, cuning
double precision xvart, xk, xeps, beta1
double precision turb_schmidt
 
character*80     fname
 
double precision, dimension(:), pointer :: cpro_taup
double precision, dimension(:), pointer :: cpro_taufpt
double precision, dimension(:), pointer :: cpro_rom
double precision, dimension(:), pointer :: cvar_k, cvar_ep, cvar_omg
double precision, dimension(:), pointer :: cvar_r11, cvar_r22, cvar_r33
double precision, dimension(:), pointer :: cpro_viscl, cpro_vscal
double precision, dimension(:), pointer :: cvar_scalt
 

Initialization and finalization

The following initialization block needs to be added for the following examples:

call field_get_val_s(iviscl, cpro_viscl)
call field_get_val_s(icrom, cpro_rom)
 
! Key id for drift scalar
call field_get_key_id("drift_scalar_model", keydri)
 
! Key id for scalar id
call field_get_key_id("scalar_id", keysca)
 
! Number of fields
call field_get_n_fields(nfld)

In theory Fortran 95 deallocates locally-allocated arrays automatically, but deallocating arrays in a symmetric manner to their allocation is good practice, and it avoids using a different logic for C and Fortran.

Body

This example set the scalar laminar diffusivity (for Brownian motion) to take thermophoresis into account.

Here is the corresponding code:

 
! Loop over fields which are scalar with a drift
do iflid = 0, nfld-1
 
  call field_get_key_int(iflid, keydri, iscdri)
 
  ! We only handle here scalar with a drift
  if (btest(iscdri, drift_scalar_add_drift_flux)) then
 
    ! Position of variables, coefficients
    ! -----------------------------------
 
    ! Index of the scalar
    call field_get_key_int(iflid, keysca, iscal)
 
    ! --- Number of the thermal variable
    if (iscalt.gt.0) then
      ivart = isca(iscalt)
      call field_get_val_s(ivarfl(ivart), cvar_scalt)
    else
      write(nfecra,9010) iscalt
      call csexit (1)
    endif
 
    ! --- Scalar's diffusivity (Brownian motion)
 
    call field_get_key_int(ivarfl(isca(iscal)), kivisl, ifcvsl)
    if (ifcvsl.ge.0) then
      call field_get_val_s(ifcvsl, cpro_vscal)
    else
      cpro_vscal => null()
    endif
 
    ! --- Coefficients of drift scalar CHOSEN BY THE USER
    !       Values given here are fictitious
 
    ! diamp: is the diameter of the particle class
    ! cuning: is the Cuningham correction factor
    ! rhop: particle density
    diamp = 1.d-4
    cuning = 1.d0
    rhop = 1.d4
 
    ! Name of the scalar with a drift
    call field_get_name(iflid, fname)
 
    ! Index of the corresponding relaxation time (cpro_taup)
    call field_get_id('drift_tau_'//trim(fname), f_id)
    call field_get_val_s(f_id, cpro_taup)
 
    ! Index of the corresponding interaction time particle--eddies (cpro_taufpt)
    if (btest(iscdri, drift_scalar_turbophoresis)) then
      call field_get_id('drift_turb_tau_'//trim(fname), f_id)
      call field_get_val_s(f_id, cpro_taufpt)
    endif
 
    ! Computation of the relaxation time of the particles
    !----------------------------------------------------
 
    if (diamp.le.1.d-6) then
      ! Cuningham's correction for submicronic particules
      do iel = 1, ncel
        cpro_taup(iel) = cuning*diamp**2*rhop/(18.d0*cpro_viscl(iel))
      enddo
    else
      do iel = 1, ncel
        cpro_taup(iel) = diamp**2*rhop/(18.d0*cpro_viscl(iel))
      enddo
    endif
 
    ! Compute the interaction time particle--eddies (tau_fpt)
    !--------------------------------------------------------
 
    if (btest(iscdri, drift_scalar_turbophoresis)) then
 
      ! k-epsilon or v2-f models
      if (itytur.eq.2 .or. itytur.eq.5) then
        call field_get_val_s(ivarfl(ik), cvar_k)
        call field_get_val_s(ivarfl(iep), cvar_ep)
        call field_get_key_double(ivarfl(isca(iscal)), ksigmas, turb_schmidt)
        do iel = 1, ncel
          xk = cvar_k(iel)
          xeps = cvar_ep(iel)
          cpro_taufpt(iel) = (3.d0/2.d0)*(cmu/turb_schmidt)*xk/xeps
        enddo
 
      ! Rij-epsilon models
      else if (itytur.eq.3) then
        call field_get_val_s(ivarfl(ir11), cvar_r11)
        call field_get_val_s(ivarfl(ir22), cvar_r22)
        call field_get_val_s(ivarfl(ir33), cvar_r33)
        call field_get_val_s(ivarfl(iep), cvar_ep)
        beta1  = 0.5d0+3.d0/(4.d0*xkappa)
        do iel = 1, ncel
          xk = 0.5d0*( cvar_r11(iel)                    &
                      +cvar_r22(iel)                    &
                      +cvar_r33(iel) )
          xeps = cvar_ep(iel)
          cpro_taufpt(iel) = xk/xeps/beta1
        enddo
 
      ! k-omega models
      else if (iturb.eq.60) then
        call field_get_val_s(ivarfl(ik), cvar_k)
        call field_get_val_s(ivarfl(iomg), cvar_omg)
        call field_get_key_double(ivarfl(isca(iscal)), ksigmas, turb_schmidt)
        do iel = 1, ncel
          xk = cvar_k(iel)
          xeps = cmu*xk*cvar_omg(iel)
          cpro_taufpt(iel) = (3.d0/2.d0)*(cmu/turb_schmidt)*xk/xeps
        enddo
      endif
 
    endif
 
    ! Brownian diffusion at cell centers
    !-----------------------------------
 
    ! --- Stop if the diffusivity is not variable
    call field_get_key_int(ivarfl(isca(iscal)), kivisl, ifcvsl)
    if (ifcvsl.lt.0) then
      write(nfecra,1010) iscal
      call csexit (1)
    endif
 
    ! Homogeneous to a dynamic viscosity
    do iel = 1, ncel
      xvart = cvar_scalt(iel)
      rho = cpro_rom(iel)
      viscl = cpro_viscl(iel)
      cpro_vscal(iel) = rho*kboltz*xvart*cuning/(3.d0*pi*diamp*viscl)
    enddo
 
  endif ! --- Tests on drift scalar
enddo