Add residual due to source terms of SA turbulence model
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in), | dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var) | :: | qp | ||
| real(kind=wp), | intent(inout), | dimension(:, :, :, :) | :: | residue | Store residue at each 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+3,-2:dims%jmx+2,-2:dims%kmx+2) | :: | Ifaces | Input varaible which stores I faces' area and unit normal |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2) | :: | Jfaces | Input varaible which stores J faces' area and unit normal |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3) | :: | Kfaces | Input varaible which stores K faces' area and unit normal |
|
| type(extent), | intent(in) | :: | dims | Extent of the domain:imx,jmx,kmx |
subroutine add_sa_source(qp, residue, cells, Ifaces,Jfaces,Kfaces, dims)
!< Add residual due to source terms of SA turbulence model
implicit none
type(extent), intent(in) :: dims
!< Extent of the domain:imx,jmx,kmx
real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var), intent(in) :: qp
real(wp), dimension(:, :, :, :), intent(inout) :: residue
!< Store residue at each 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+3,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: Ifaces
!< Input varaible which stores I faces' area and unit normal
type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2), intent(in) :: Jfaces
!< Input varaible which stores J faces' area and unit normal
type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3), intent(in) :: Kfaces
!< Input varaible which stores K faces' area and unit normal
integer :: i,j,k
real(wp) :: CD1
real(wp) :: CD2
real(wp) :: fv1
real(wp) :: fv2
real(wp) :: fw
real(wp) :: g
real(wp) :: Scap
real(wp) :: r
real(wp) :: vort
real(wp) :: S_v
real(wp) :: D_v
real(wp) :: P_v
real(wp) :: lamda
real(wp) :: kd2
real(wp) :: xi
real(wp) :: nu
real(wp) :: gradrho_x
real(wp) :: gradrho_y
real(wp) :: gradrho_z
real(wp), dimension(6) :: RhoFace
real(wp), dimension(6) :: Area
real(wp), dimension(6,3) :: Normal
real(wp) :: density
real(wp) :: tv
do k = 1,dims%kmx-1
do j = 1,dims%jmx-1
do i = 1,dims%imx-1
density = qp(i,j,k,1)
tv = qp(i,j,k,6)
RhoFace(1) = qp(i-1,j ,k ,1)+density
RhoFace(2) = qp(i ,j-1,k ,1)+density
RhoFace(3) = qp(i ,j ,k-1,1)+density
RhoFace(4) = qp(i+1,j ,k ,1)+density
RhoFace(5) = qp(i ,j+1,k ,1)+density
RhoFace(6) = qp(i ,j ,k+1,1)+density
Area(1) = Ifaces(i,j,k)%A
Area(2) = Jfaces(i,j,k)%A
Area(3) = Kfaces(i,j,k)%A
Area(4) = Ifaces(i+1,j ,k )%A
Area(5) = Jfaces(i ,j+1,k )%A
Area(6) = Kfaces(i ,j ,k+1)%A
Normal(1,1:3) = (/Ifaces(i,j,k)%nx,Ifaces(i,j,k)%ny,Ifaces(i,j,k)%nz/)
Normal(2,1:3) = (/Jfaces(i,j,k)%nx,Jfaces(i,j,k)%ny,Jfaces(i,j,k)%nz/)
Normal(3,1:3) = (/Kfaces(i,j,k)%nx,Kfaces(i,j,k)%nx,Kfaces(i,j,k)%nx/)
Normal(4,1:3) = (/Ifaces(i+1,j ,k)%nx,Ifaces(i+1,j ,k)%ny,Ifaces(i+1,j ,k)%nz/)
Normal(5,1:3) = (/Jfaces(i ,j+1,k)%nx,Jfaces(i ,j+1,k)%ny,Jfaces(i ,j+1,k)%nz/)
Normal(6,1:3) = (/Kfaces(i ,j,k+1)%nx,Kfaces(i ,j,k+1)%ny,Kfaces(i ,j,k+1)%nz/)
gradrho_x = (-(RhoFace(1))*Normal(1,1)*Area(1) &
-(RhoFace(2))*Normal(2,1)*Area(2) &
-(RhoFace(3))*Normal(3,1)*Area(3) &
+(RhoFace(4))*Normal(4,1)*Area(4) &
+(RhoFace(5))*Normal(5,1)*Area(5) &
+(RhoFace(6))*Normal(6,1)*Area(6) &
)/(2.0*cells(i,j,k)%volume)
gradrho_y = (-(RhoFace(1))*Normal(1,2)*Area(1) &
-(RhoFace(2))*Normal(2,2)*Area(2) &
-(RhoFace(3))*Normal(3,2)*Area(3) &
+(RhoFace(4))*Normal(4,2)*Area(4) &
+(RhoFace(5))*Normal(5,2)*Area(5) &
+(RhoFace(6))*Normal(6,2)*Area(6) &
)/(2.0*cells(i,j,k)%volume)
gradrho_z = (-(RhoFace(1))*Normal(1,3)*Area(1) &
-(RhoFace(2))*Normal(2,3)*Area(2) &
-(RhoFace(3))*Normal(3,3)*Area(3) &
+(RhoFace(4))*Normal(4,3)*Area(4) &
+(RhoFace(5))*Normal(5,3)*Area(5) &
+(RhoFace(6))*Normal(6,3)*Area(6) &
)/(2.0*cells(i,j,k)%volume)
! __ vorticity __
vort = sqrt( ((gradw_y(i,j,k)- gradv_z(i,j,k))**2 &
+ (gradu_z(i,j,k)- gradw_x(i,j,k))**2 &
+ (gradv_x(i,j,k)- gradu_y(i,j,k))**2 &
)&
)
! ___ cross diffusion ___
CD1 = cb2*((gradtv_x(i,j,k)*gradtv_x(i,j,k))&
+ (gradtv_y(i,j,k)*gradtv_y(i,j,k))&
+ (gradtv_z(i,j,k)*gradtv_z(i,j,k))&
)
! ___ addition cross diffusion result conservative form of tv ___
CD2 = ((gradrho_x*gradtv_x(i,j,k))&
+ (gradrho_y*gradtv_y(i,j,k))&
+ (gradrho_z*gradtv_z(i,j,k))&
)
kd2 = (kappa_sa*dist(i,j,k))**2
nu = mu(i,j,k)/density
xi = tv/nu
! ___ functions ___
fv1 = (xi**3)/((xi**3) + (cv1**3))
fv2 = 1.0 - xi/(1.0 + (xi*fv1))
! ___ Shear stress for production ___
scap = max(vort + (tv*fv2/(kd2)), 0.3*vort)
! ___ wall function ___
r = min(tv/(Scap*kd2), 10.0)
g = r + cw2*((r**6) - r)
fw = g*( (1.0+(cw3**6))/((g**6)+(cw3**6)) )**(1.0/6.0)
! ____ Dissipation term ___
D_v = density*cw1*fw*((tv/dist(i,j,k))**2)
! ____ PRODUCTION term____
P_v = density*cb1*Scap*tv
! ____ cross diffusion term ___
lamda = density*CD1/sigma_sa - CD2*(nu+tv)/sigma_sa
S_v = (P_v - D_v + lamda)*cells(i,j,k)%volume
residue(i, j, k, 6) = residue(i, j, k, 6) - S_v
end do
end do
end do
end subroutine add_sa_source