Compute viscous fluxes for first five Navier-Stokes equation
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(inout), | dimension(:, :, :, :) | :: | F | Flux array |
|
| real(kind=wp), | intent(in), | dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var) | :: | qp | Store primitive variable at cell center |
|
| type(celltype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2) | :: | cells | Input cell quantities: volume |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2+flags(1),-2:dims%jmx+2+flags(2),-2:dims%kmx+2+flags(3)) | :: | faces | Face quantities: area and unit normal |
|
| integer, | intent(in), | dimension(3) | :: | flags | flags for direction switch |
|
| type(schemetype), | intent(in) | :: | scheme | finite-volume Schemes |
||
| type(flowtype), | intent(in) | :: | flow | Information about fluid flow: freestream-speed, ref-viscosity,etc. |
||
| type(extent), | intent(in) | :: | dims | Extent of the domain:imx,jmx,kmx |
subroutine compute_viscous_fluxes_laminar(F, qp, cells, faces, flags, scheme, flow, dims)
!< Compute viscous fluxes for first five Navier-Stokes equation
implicit none
real(wp), dimension(:, :, :, :), intent(inout) :: F
!< Flux array
type(schemetype), intent(in) :: scheme
!< finite-volume Schemes
type(flowtype), intent(in) :: flow
!< Information about fluid flow: freestream-speed, ref-viscosity,etc.
type(extent), intent(in) :: dims
!< Extent of the domain:imx,jmx,kmx
integer, dimension(3), intent(in) :: flags
!< flags for direction switch
real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var), intent(in) :: qp
!< Store primitive variable at cell center
type(celltype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: cells
!< Input cell quantities: volume
type(facetype), dimension(-2:dims%imx+2+flags(1),-2:dims%jmx+2+flags(2),-2:dims%kmx+2+flags(3)), intent(in) :: faces
!< Face quantities: area and unit normal
! local variables
real(wp) :: dudx, dudy, dudz
real(wp) :: dvdx, dvdy, dvdz
real(wp) :: dwdx, dwdy, dwdz
real(wp) :: dTdx, dTdy, dTdz
real(wp) :: normal_comp
real(wp) :: d_LR
real(wp) :: T_RE
real(wp) :: T_LE
real(wp) :: K_heat, Qx, Qy, Qz
real(wp) :: mu_f
real(wp) :: mut_f
real(wp) :: total_mu
real(wp) :: delx
real(wp) :: dely
real(wp) :: delz
real(wp) :: delu
real(wp) :: delv
real(wp) :: delw
real(wp) :: delT
real(wp) :: Tau_xx
real(wp) :: Tau_xy
real(wp) :: Tau_xz
real(wp) :: Tau_yx
real(wp) :: Tau_yy
real(wp) :: Tau_yz
real(wp) :: Tau_zx
real(wp) :: Tau_zy
real(wp) :: Tau_zz
real(wp) :: nx
real(wp) :: ny
real(wp) :: nz
real(wp) :: area
real(wp) :: uface
real(wp) :: vface
real(wp) :: wface
integer :: i, j, k
integer :: ii, jj, kk
ii = flags(1)
jj = flags(2)
kk = flags(3)
!---------------------------------------------------------------------
! Calculating the fluxes at the faces
!--------------------------------------------------------------------
do k = 1, dims%kmx - 1 + kk
do j = 1, dims%jmx - 1 + jj
do i = 1, dims%imx - 1 + ii
!--- FACE Gradients ---!
! Gradients at face as average of gradients at cell centres
dudx = 0.5 * (gradu_x(i-ii, j-jj, k-kk) + gradu_x(i, j, k))
dudy = 0.5 * (gradu_y(i-ii, j-jj, k-kk) + gradu_y(i, j, k))
dudz = 0.5 * (gradu_z(i-ii, j-jj, k-kk) + gradu_z(i, j, k))
dvdx = 0.5 * (gradv_x(i-ii, j-jj, k-kk) + gradv_x(i, j, k))
dvdy = 0.5 * (gradv_y(i-ii, j-jj, k-kk) + gradv_y(i, j, k))
dvdz = 0.5 * (gradv_z(i-ii, j-jj, k-kk) + gradv_z(i, j, k))
dwdx = 0.5 * (gradw_x(i-ii, j-jj, k-kk) + gradw_x(i, j, k))
dwdy = 0.5 * (gradw_y(i-ii, j-jj, k-kk) + gradw_y(i, j, k))
dwdz = 0.5 * (gradw_z(i-ii, j-jj, k-kk) + gradw_z(i, j, k))
dTdx = 0.5 * (gradT_x(i-ii, j-jj, k-kk) + gradT_x(i, j, k))
dTdy = 0.5 * (gradT_y(i-ii, j-jj, k-kk) + gradT_y(i, j, k))
dTdz = 0.5 * (gradT_z(i-ii, j-jj, k-kk) + gradT_z(i, j, k))
!--- For ODD-EVEN coupling error ---!
! distance between cell center of adjacent cell for the i,j,k face
delx = cells(i, j, k)%centerx - cells(i-ii, j-jj, k-kk)%centerx
dely = cells(i, j, k)%centery - cells(i-ii, j-jj, k-kk)%centery
delz = cells(i, j, k)%centerz - cells(i-ii, j-jj, k-kk)%centerz
d_LR = sqrt(delx*delx + dely*dely + delz*delz)
! Finding the temperature of left and right element to the face i,j,k
T_LE = qp(i-ii, j-jj, k-kk, 5) / (qp(i-ii, j-jj, k-kk, 1) * flow%R_gas)
T_RE = qp(i, j, k, 5) / (qp(i, j, k, 1) * flow%R_gas)
! difference in state across face
delu = qp(i, j, k, 2) - qp(i-ii, j-jj, k-kk, 2) !x_speed
delv = qp(i, j, k, 3) - qp(i-ii, j-jj, k-kk, 3) !y_speed
delw = qp(i, j, k, 4) - qp(i-ii, j-jj, k-kk, 4) !z_speed
delT = T_RE - T_LE
!normal_comp = ( delta(phi) - (grad(phi).dot.delR) )/magnitudeR
!new grad(phi) = grad(phi) + correction(normal_comp.dot.delR/magnitudeR)
normal_comp = (delu - (dudx*delx + dudy*dely + dudz*delz))/d_LR
dudx = dudx + (normal_comp * delx / d_LR)
dudy = dudy + (normal_comp * dely / d_LR)
dudz = dudz + (normal_comp * delz / d_LR)
normal_comp = (delv - (dvdx*delx + dvdy*dely + dvdz*delz))/d_LR
dvdx = dvdx + (normal_comp * delx / d_LR)
dvdy = dvdy + (normal_comp * dely / d_LR)
dvdz = dvdz + (normal_comp * delz / d_LR)
normal_comp = (delw - (dwdx*delx + dwdy*dely + dwdz*delz))/d_LR
dwdx = dwdx + (normal_comp * delx / d_LR)
dwdy = dwdy + (normal_comp * dely / d_LR)
dwdz = dwdz + (normal_comp * delz / d_LR)
normal_comp = (delT - (dTdx*delx + dTdy*dely + dTdz*delz))/d_LR
dTdx = dTdx + (normal_comp * delx / d_LR)
dTdy = dTdy + (normal_comp * dely / d_LR)
dTdz = dTdz + (normal_comp * delz / d_LR)
!--- end of ODD-EVEN coupling correction ---!
mu_f = 0.5*(mu(i-ii, j-jj, k-kk) + mu(i,j,k))
if(trim(scheme%turbulence)/='none') then
mut_f = 0.5*(mu_t(i-ii, j-jj, k-kk) + mu_t(i, j, k))
else
mut_f = 0.0
end if
! effective viscosity
total_mu = mu_f + mut_f
! Using lambda = -2 * mu / 3
! diagonal terms of stress tensor
Tau_xx = 2. * total_mu * (dudx - ((dudx + dvdy + dwdz) / 3.))
Tau_yy = 2. * total_mu * (dvdy - ((dudx + dvdy + dwdz) / 3.))
Tau_zz = 2. * total_mu * (dwdz - ((dudx + dvdy + dwdz) / 3.))
! off diagonal symmetrical part of stress tensor
Tau_xy = total_mu * (dvdx + dudy)
Tau_xz = total_mu * (dwdx + dudz)
Tau_yz = total_mu * (dwdy + dvdz)
Tau_yx = Tau_xy
Tau_zx = Tau_xz
Tau_zy = Tau_yz
! Pr: Prandtl Number and tPr: Turbulent Prandtl number
! Qx, Qy, Qz: Conduction fluxes
K_heat = (mu_f/flow%Pr + mut_f/flow%tPr)* flow%gm * flow%R_gas / (flow%gm - 1)
Qx = K_heat*dTdx
Qy = K_heat*dTdy
Qz = K_heat*dTdz
! calling some element from memory and keep them handy for calculation
nx = faces(i,j,k)%nx
ny = faces(i,j,k)%ny
nz = faces(i,j,k)%nz
area = faces(i,j,k)%A
uface = 0.5 * (qp(i-ii, j-jj, k-kk, 2) + qp(i, j, k, 2))
vface = 0.5 * (qp(i-ii, j-jj, k-kk, 3) + qp(i, j, k, 3))
wface = 0.5 * (qp(i-ii, j-jj, k-kk, 4) + qp(i, j, k, 4))
! adding viscous fluxes to stored convective flux
F(i, j, k, 2) = F(i, j, k, 2) - ((Tau_xx*nx + Tau_xy*ny + Tau_xz*nz)*area)
F(i, j, k, 3) = F(i, j, k, 3) - ((Tau_yx*nx + Tau_yy*ny + Tau_yz*nz)*area)
F(i, j, k, 4) = F(i, j, k, 4) - ((Tau_zx*nx + Tau_zy*ny + Tau_zz*nz)*area)
! Energy flux
! (TijVj + Qi)ni
F(i, j, k, 5) = F(i, j, k, 5) - (area * ( &
((Tau_xx*uface + Tau_xy*vface + Tau_xz*wface + Qx)*nx) + &
((Tau_yx*uface + Tau_yy*vface + Tau_yz*wface + Qy)*ny) + &
((Tau_zx*uface + Tau_zy*vface + Tau_zz*wface + Qz)*nz) ) )
end do
end do
end do
end subroutine compute_viscous_fluxes_laminar