Calculate the total flux through face for turbulent flow (SA)
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| 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 |
function SAFlux(ql, qr, du, inputs)
!< Calculate the total flux through face for turbulent flow (SA)
!--------------------------------------
! calculate the total flux through face
!---------------------------------------
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) :: SAFlux
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) :: dtvdx
real(wp) :: dtvdy
real(wp) :: dtvdz
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) :: muCap
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(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(6) = max(W(6), 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) )
! viscous terms
muCap = 0.25*(P(1)+W(1))*(P(6) + W(6))
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
dtvdx = ( P(6) - W(6) ) * nx * Area / Volume
dtvdy = ( P(6) - W(6) ) * ny * Area / Volume
dtvdz = ( P(6) - W(6) ) * 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
tmu = 0.5*(W(6) + P(6))
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 + muCap)*(dtvdx*nx + dtvdy*ny + dtvdz*nz)/sigma_sa
Flux = Flux * Area
SAFlux = Flux
end function SAFlux