Add residual due to source terms of SABC transition 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 | Store primitive variable at cell center |
|
| 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(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 add_saBC_source(qp, residue, cells, Ifaces,Jfaces,Kfaces, flow, dims)
!< Add residual due to source terms of SABC transition model
implicit none
type(extent), intent(in) :: dims
!< Extent of the domain:imx,jmx,kmx
type(flowtype), intent(in) :: flow
!< Information about fluid flow: freestream-speed, ref-viscosity,etc.
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
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) :: r
real(wp) :: S_v
real(wp) :: lamda
real(wp) :: dist_i
real(wp) :: dist_i_2
real(wp) :: Ji
real(wp) :: Ji_2
real(wp) :: Ji_3
real(wp) :: S
real(wp) :: Omega
real(wp) :: k2
real(wp) :: inv_k2_d2
real(wp) :: Shat
real(wp) :: inv_Shat
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
! transition modeling variables
real(wp) :: chi_1=0.002
real(wp) :: chi_2=5.0
real(wp) :: nu_BC
real(wp) :: nu_cr
real(wp) :: nu_t
real(wp) :: u,v,w
real(wp) :: glim
real(wp) :: g_6
real(wp) :: vmag
real(wp) :: Production
real(wp) :: Destruction
real(wp) :: re_v
real(wp) :: re_theta
real(wp) :: re_theta_t
real(wp) :: term1
real(wp) :: term2
real(wp) :: term_exponential
real(wp) :: gamma_BC
real(wp) :: tu
real(wp) :: tv
real(wp) :: density
tu = flow%tu_inf
do k = 1,dims%kmx-1
do j = 1,dims%jmx-1
do i = 1,dims%imx-1
!Local_vel_mag
density= qp(i,j,k,1)
u = qp(i,j,k,2)
v = qp(i,j,k,3)
w = qp(i,j,k,4)
tv = qp(i,j,k,6)
vmag = sqrt(u*u + v*v + w*w)
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*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*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*cells(i,j,k)%volume)
! __ vorticity __
Omega = 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))&
)
dist_i = dist(i,j,k)
dist_i_2 = dist_i*dist_i
k2 = kappa_sa*kappa_sa
nu = mu(i,j,k)/density
Ji = tv/nu
Ji_2 = Ji*Ji
Ji_3 = Ji_2*ji
! ___ functions ___
fv1 = (Ji_3)/((Ji_3) + (cv1_3))
fv2 = 1.0 - Ji/(1.0 + (Ji*fv1))
! ___ Shear stress for production ___
S = Omega
inv_k2_d2 = 1.0/(k2*dist_i_2)
Shat = S + tv*fv2*inv_k2_d2
Shat = max(Shat, 1.0e-10)
inv_Shat = 1.0/Shat
! ____ PRODUCTION term____
chi_1 = 0.002
chi_2 = 5.0
nu_t = tv*fv1
nu_cr = chi_2/flow%Reynolds_number
nu_bc = nu_t/(vmag*dist_i)
re_v = dist_i_2*Omega/nu
re_theta = re_v/2.193
re_theta_t = (803.73*((tu+0.6067)**(-1.027)))
!re_theta_t = 163.0 + exp(6.91-0.18)
term1 = sqrt(max(re_theta-re_theta_t,0.)/(chi_1*re_theta_t))
term2 = sqrt(max(nu_BC-nu_cr,0.0)/nu_cr)
term_exponential = (term1 + term2)
gamma_BC = 1.0 - exp(-term_exponential)
Production = gamma_BC*cb1*Shat*tv*cells(i,j,k)%volume
! ___ Destruction term___ !
r = min(tv*inv_Shat*inv_k2_d2, 10.0)
g = r + cw2*((r**6) - r)
g_6 = g**6
glim = ((1.0+cw3_6)/(g_6+cw3_6))**(1.0/6.0)
fw = g*glim
Destruction = (cw1*fw*tv*tv/dist_i_2)*(cells(i,j,k)%volume)
! ____ cross diffusion term ___
lamda = (density*CD1/sigma_sa - CD2*(nu+tv)/sigma_sa)*cells(i,j,k)%volume
S_v = (Production - Destruction + lamda)
residue(i, j, k, 6) = residue(i, j, k, 6) - S_v
end do
end do
end do
end subroutine add_saBC_source