KKLFlux Function

public function KKLFlux(ql, qr, du, inputs)

calculate the total flux through face for turbulent flow (k-kL)

Arguments

Type IntentOptional AttributesName
real(kind=wp), intent(in), dimension(1:n_var):: ql
real(kind=wp), intent(in), dimension(1:n_var):: qr
real(kind=wp), intent(in), dimension(1:n_var):: du
real(kind=wp), intent(in), dimension(1:7):: inputs

Return Value real(kind=wp), dimension(1:n_var)


Called by

proc~~kklflux~~CalledByGraph proc~kklflux KKLFlux proc~update_kkl_variables update_KKL_variables proc~update_kkl_variables->proc~kklflux proc~update_with_lusgs update_with_lusgs proc~update_with_lusgs->proc~update_kkl_variables proc~get_next_solution get_next_solution proc~get_next_solution->proc~update_with_lusgs proc~iterate_one_more_time_step iterate_one_more_time_step proc~iterate_one_more_time_step->proc~get_next_solution program~main main program~main->proc~iterate_one_more_time_step

Contents

Source Code


Source Code

    function KKLFlux(ql, qr, du, inputs)
      !< calculate the total flux through face for turbulent flow (k-kL)
      implicit none
      real(wp), dimension(1:n_var), intent(in) :: ql !left state
      real(wp), dimension(1:n_var), intent(in) :: qr !right state
      !conservative form of updated neighbour
      real(wp), dimension(1:n_var), intent(in) :: du
      real(wp), dimension(1:7)    , intent(in) :: inputs
      real(wp), dimension(1:n_var)             :: Flux
      real(wp), dimension(1:n_var)             :: KKLFlux
      real(wp), dimension(1:n_var)             :: U ! conservative variables
      real(wp), dimension(1:n_var)             :: W ! new primitive variables
      real(wp), dimension(1:n_var)             :: P ! primitive variables of right cell

      !for extraction of the inputs
      real(wp) :: area
      real(wp) :: nx
      real(wp) :: ny
      real(wp) :: nz
      real(wp) :: volume
      real(wp) :: mmu
      real(wp) :: tmu


      real(wp)    :: dudx
      real(wp)    :: dudy
      real(wp)    :: dudz
      real(wp)    :: dvdx
      real(wp)    :: dvdy
      real(wp)    :: dvdz
      real(wp)    :: dwdx
      real(wp)    :: dwdy
      real(wp)    :: dwdz
      real(wp)    :: dTdx
      real(wp)    :: dTdy
      real(wp)    :: dTdz
      real(wp)    :: dtkdx
      real(wp)    :: dtkdy
      real(wp)    :: dtkdz
      real(wp)    :: dtkldx
      real(wp)    :: dtkldy
      real(wp)    :: dtkldz
      real(wp)    :: T1, T2
      real(wp)    :: uface
      real(wp)    :: vface
      real(wp)    :: wface
      real(wp)    :: trace
      real(wp)    :: Tauxx
      real(wp)    :: Tauyy
      real(wp)    :: Tauzz
      real(wp)    :: Tauxy
      real(wp)    :: Tauxz
      real(wp)    :: Tauyz
      real(wp)    :: Qx
      real(wp)    :: Qy
      real(wp)    :: Qz
      real(wp)    :: HalfRhoUsquare
      real(wp)    :: RhoHt
      real(wp)    :: K_heat
      real(wp)    :: FaceNormalVelocity
      real(wp)    :: mu

      area   = inputs(1)
      nx     = inputs(2)
      ny     = inputs(3)
      nz     = inputs(4)
      volume = inputs(5)
      mmu    = inputs(6)
      tmu    = inputs(7)


      !save the old stat in P
      P = qr

      ! find conservative variable
      U(1)   =   ql(1)
      U(2)   =   ql(1) * ql(2)
      U(3)   =   ql(1) * ql(3)
      U(4)   =   ql(1) * ql(4)
      U(5)   = ( ql(5) / (gm-1.0) ) + ( 0.5 * ql(1) * sum(ql(2:4)**2) )
      U(6)   =   ql(1) * ql(6)
      U(7)   =   ql(1) * ql(7)

      U(1:n_var) = U(1:n_var) + du(1:n_var)
      

      W(1)   =   U(1)
      W(2)   =   U(2) / U(1)
      W(3)   =   U(3) / U(1)
      W(4)   =   U(4) / U(1)
      W(5)   = (gm-1.0) * ( U(5) - ( 0.5 * SUM(U(2:4)**2) / U(1) ) )
      W(6)   =   U(6) / U(1)
      W(7)   =   U(7) / U(1)
      W(6) = max(W(6), 1e-8)
      W(7) = max(W(7), 1e-8)

      FaceNormalVelocity = (W(2) * nx) + (W(3) * ny) + (W(4) * nz)
      uface = 0.5 * ( W(2) + P(2) )
      vface = 0.5 * ( W(3) + P(3) )
      wface = 0.5 * ( W(4) + P(4) )


      Flux(1) =   W(1) * FaceNormalVelocity
      Flux(2) = ( W(2) * Flux(1) ) + ( W(5) * nx )
      Flux(3) = ( W(3) * Flux(1) ) + ( W(5) * ny )
      Flux(4) = ( W(4) * Flux(1) ) + ( W(5) * nz )

      HalfRhoUsquare = 0.5 * W(1) * ( W(2)*W(2) + W(3)*W(3) + W(4)*W(4) )
      RhoHt          = ( (gm/(gm-1.0)) * W(5) ) + HalfRhoUsquare
      Flux(5)        = RhoHt * FaceNormalVelocity
      Flux(6) = ( W(6) * Flux(1) )   
      Flux(7) = ( W(7) * Flux(1) )   


      ! viscous terms
      mu = mmu + tmu
      T1     =    W(5) / ( W(1) * R_gas )
      T2     =    P(5) / ( P(1) * R_gas )
      dTdx   =  ( T2   - T1   ) * nx * Area / Volume
      dTdy   =  ( T2   - T1   ) * ny * Area / Volume
      dTdz   =  ( T2   - T1   ) * nz * Area / Volume
      dudx   =  ( P(2) - W(2) ) * nx * Area / Volume
      dudy   =  ( P(2) - W(2) ) * ny * Area / Volume
      dudz   =  ( P(2) - W(2) ) * nz * Area / Volume
      dvdx   =  ( P(3) - W(3) ) * nx * Area / Volume
      dvdy   =  ( P(3) - W(3) ) * ny * Area / Volume
      dvdz   =  ( P(3) - W(3) ) * nz * Area / Volume
      dwdx   =  ( P(4) - W(4) ) * nx * Area / Volume
      dwdy   =  ( P(4) - W(4) ) * ny * Area / Volume
      dwdz   =  ( P(4) - W(4) ) * nz * Area / Volume
      dtkdx  =  ( P(6) - W(6) ) * nx * Area / Volume
      dtkdy  =  ( P(6) - W(6) ) * ny * Area / Volume
      dtkdz  =  ( P(6) - W(6) ) * nz * Area / Volume
      dtkldx  =  ( P(7) - W(7) ) * nx * Area / Volume
      dtkldy  =  ( P(7) - W(7) ) * ny * Area / Volume
      dtkldz  =  ( P(7) - W(7) ) * nz * Area / Volume

      trace = dudx + dvdy + dwdz
      Tauxx =  2. * mu * (dudx - trace/3.0)
      Tauyy =  2. * mu * (dvdy - trace/3.0)
      Tauzz =  2. * mu * (dwdz - trace/3.0)
      Tauxy = mu * (dvdx + dudy)
      Tauxz = mu * (dwdx + dudz)
      Tauyz = mu * (dwdy + dvdz)

      K_heat = ( mmu / Pr  + tmu/tpr) * gm * R_gas / ( gm - 1.0 )
      Qx = K_heat*dTdx
      Qy = K_heat*dTdy
      Qz = K_heat*dTdz

      Flux(2) = Flux(2) - ( Tauxx * nx + Tauxy * ny + Tauxz * nz )
      Flux(3) = Flux(3) - ( Tauxy * nx + Tauyy * ny + Tauyz * nz )
      Flux(4) = Flux(4) - ( Tauxz * nx + Tauyz * ny + Tauzz * nz )
      Flux(5) = Flux(5) - ( Tauxx * uface + Tauxy * vface + Tauxz * wface + Qx ) * nx
      Flux(5) = Flux(5) - ( Tauxy * uface + Tauyy * vface + Tauyz * wface + Qy ) * ny
      Flux(5) = Flux(5) - ( Tauxz * uface + Tauyz * vface + Tauzz * wface + Qz ) * nz
      Flux(6) = Flux(6) + (mmu + sigma_k*tmu)*(dtkdx*nx + dtkdy*ny + dtkdz*nz)
      Flux(7) = Flux(7) + (mmu + sigma_phi*tmu)*(dtkldx*nx + dtkldy*ny + dtkldz*nz)

      Flux    = Flux * Area
      KKLFlux = Flux

    end function KKLFlux