lctm2015flux Function

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

calculate the total flux through face for turbulent/transition flow (LCTM2015)

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:8):: inputs

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


Called by

proc~~lctm2015flux~2~~CalledByGraph proc~lctm2015flux~2 lctm2015flux proc~update_lctm2015~2 update_lctm2015 proc~update_lctm2015~2->proc~lctm2015flux~2 proc~update_with_lusgs update_with_lusgs proc~update_with_lusgs->proc~update_lctm2015~2 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 lctm2015flux(ql, qr, du, inputs)
      !< calculate the total flux through face for turbulent/transition flow (LCTM2015)
      !---------------------------------------
      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:8)    , intent(in) :: inputs
      real(wp), dimension(1:n_var)             :: Flux
      real(wp), dimension(1:n_var)             :: lctm2015flux
      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)    :: dtwdx
      real(wp)    :: dtwdy
      real(wp)    :: dtwdz
      real(wp)    :: dtgmdx
      real(wp)    :: dtgmdy
      real(wp)    :: dtgmdz
      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
      real(wp)    :: sigma_k
      real(wp)    :: sigma_w
      real(wp)    :: F1

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


      !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(8)   =   ql(1) * ql(8)

      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(8)   =   U(8) / U(1)
      W(6)   = W(6) + 0.5*(1.-sign(1.,W(6)))*(ql(6)-W(6))
      W(7)   = W(7) + 0.5*(1.-sign(1.,W(7)))*(ql(7)-W(7))
      W(8) = max(W(8), 0.0)
      !W(8) = min(W(8), 1.0)

      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) )   
      Flux(8) = ( W(8) * 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
      dtwdx  =  ( P(7) - W(7) ) * nx * Area / Volume
      dtwdy  =  ( P(7) - W(7) ) * ny * Area / Volume
      dtwdz  =  ( P(7) - W(7) ) * nz * Area / Volume
      dtgmdx  =  ( P(8) - W(8) ) * nx * Area / Volume
      dtgmdy  =  ( P(8) - W(8) ) * ny * Area / Volume
      dtgmdz  =  ( P(8) - W(8) ) * 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

      sigma_k = sigma_k1*F1 + sigma_k2*(1.0 - F1)
      sigma_w = sigma_w1*F1 + sigma_w2*(1.0 - F1)

      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_w*tmu)*(dtwdx*nx + dtwdy*ny + dtwdz*nz)
      Flux(8) = Flux(8) + (mmu + tmu)*(dtgmdx*nx + dtgmdy*ny + dtgmdz*nz)

      Flux    = Flux * Area
      lctm2015flux = Flux
    end function lctm2015flux