Add residual due to source terms of the k-kL 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_kkl_source(qp, residue, cells, Ifaces,Jfaces,Kfaces, dims)
!< Add residual due to source terms of the k-kL 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) :: Tau11
real(wp) :: Tau12
real(wp) :: Tau13
real(wp) :: Tau21
real(wp) :: Tau22
real(wp) :: Tau23
real(wp) :: Tau31
real(wp) :: Tau32
real(wp) :: Tau33
real(wp) :: S11
real(wp) :: S12
real(wp) :: S13
real(wp) :: S21
real(wp) :: S22
real(wp) :: S23
real(wp) :: S31
real(wp) :: S32
real(wp) :: S33
real(wp) :: delv
real(wp) :: d2udx2
real(wp) :: d2udy2
real(wp) :: d2udz2
real(wp) :: d2vdx2
real(wp) :: d2vdy2
real(wp) :: d2vdz2
real(wp) :: d2wdx2
real(wp) :: d2wdy2
real(wp) :: d2wdz2
real(wp) :: Lvk
real(wp) :: fp
real(wp) :: ud
real(wp) :: udd
real(wp) :: S_k
real(wp) :: S_kl
real(wp) :: D_k
real(wp) :: D_kl
real(wp) :: P_k
real(wp) :: P_kl
real(wp) :: density
real(wp) :: tk
real(wp) :: tkl
do k = 1,dims%kmx-1
do j = 1,dims%jmx-1
do i = 1,dims%imx-1
density = qp(i,j,k,1)
tk = qp(i,j,k,6)
tkl = qp(i,j,k,7)
S11 = 0.5*(gradu_x(i,j,k) + gradu_x(i,j,k))
S12 = 0.5*(gradu_y(i,j,k) + gradv_x(i,j,k))
S13 = 0.5*(gradu_z(i,j,k) + gradw_x(i,j,k))
S21 = 0.5*(gradv_x(i,j,k) + gradu_y(i,j,k))
S22 = 0.5*(gradv_y(i,j,k) + gradv_y(i,j,k))
S23 = 0.5*(gradv_z(i,j,k) + gradw_y(i,j,k))
S31 = 0.5*(gradw_x(i,j,k) + gradu_z(i,j,k))
S32 = 0.5*(gradw_y(i,j,k) + gradv_z(i,j,k))
S33 = 0.5*(gradw_z(i,j,k) + gradw_z(i,j,k))
delv = gradu_x(i,j,k) + gradv_y(i,j,k) + gradw_z(i,j,k)
Tau11 = mu_t(i,j,k)*(2*S11 - (2.0/3.0)*delv) - (2.0/3.0)*density*tk
Tau12 = mu_t(i,j,k)*(2*S12)
Tau13 = mu_t(i,j,k)*(2*S13)
Tau21 = mu_t(i,j,k)*(2*S21)
Tau22 = mu_t(i,j,k)*(2*S22 - (2.0/3.0)*delv) - (2.0/3.0)*density*tk
Tau23 = mu_t(i,j,k)*(2*S23)
Tau31 = mu_t(i,j,k)*(2*S31)
Tau32 = mu_t(i,j,k)*(2*S32)
Tau33 = mu_t(i,j,k)*(2*S33 - (2.0/3.0)*delv) - (2.0/3.0)*density*tk
P_k = 0.
P_k = P_k + Tau11*gradu_x(i,j,k) + Tau12*gradu_y(i,j,k) + Tau13*gradu_z(i,j,k)
P_k = P_k + Tau21*gradv_x(i,j,k) + Tau22*gradv_y(i,j,k) + Tau23*gradv_z(i,j,k)
P_k = P_k + Tau31*gradw_x(i,j,k) + Tau32*gradw_y(i,j,k) + Tau33*gradw_z(i,j,k)
D_k = (cmu**0.75)*density*(tk**2.5)/max(tkl,1.e-20)
P_k = min(P_k, 20*D_k)
! calculation of Lvk
! first get second order gradients
d2udx2 =(-(gradu_x(i-1,j ,k )+gradu_x(i,j,k))*Ifaces(i,j,k)%nx*Ifaces(i,j,k)%A &
-(gradu_x(i ,j-1,k )+gradu_x(i,j,k))*Jfaces(i,j,k)%nx*Jfaces(i,j,k)%A &
-(gradu_x(i ,j ,k-1)+gradu_x(i,j,k))*Kfaces(i,j,k)%nx*Kfaces(i,j,k)%A &
+(gradu_x(i+1,j ,k )+gradu_x(i,j,k))*Ifaces(i+1,j ,k )%nx*Ifaces(i+1,j ,k )%A &
+(gradu_x(i ,j+1,k )+gradu_x(i,j,k))*Jfaces(i ,j+1,k )%nx*Jfaces(i ,j+1,k )%A &
+(gradu_x(i ,j ,k+1)+gradu_x(i,j,k))*Kfaces(i ,j ,k+1)%nx*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
d2udy2 =(-(gradu_y(i-1,j ,k )+gradu_y(i,j,k))*Ifaces(i,j,k)%ny*Ifaces(i,j,k)%A &
-(gradu_y(i ,j-1,k )+gradu_y(i,j,k))*Jfaces(i,j,k)%ny*Jfaces(i,j,k)%A &
-(gradu_y(i ,j ,k-1)+gradu_y(i,j,k))*Kfaces(i,j,k)%ny*Kfaces(i,j,k)%A &
+(gradu_y(i+1,j ,k )+gradu_y(i,j,k))*Ifaces(i+1,j ,k )%ny*Ifaces(i+1,j ,k )%A &
+(gradu_y(i ,j+1,k )+gradu_y(i,j,k))*Jfaces(i ,j+1,k )%ny*Jfaces(i ,j+1,k )%A &
+(gradu_y(i ,j ,k+1)+gradu_y(i,j,k))*Kfaces(i ,j ,k+1)%ny*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
d2udz2 =(-(gradu_z(i-1,j ,k )+gradu_z(i,j,k))*Ifaces(i,j,k)%nz*Ifaces(i,j,k)%A &
-(gradu_z(i ,j-1,k )+gradu_z(i,j,k))*Jfaces(i,j,k)%nz*Jfaces(i,j,k)%A &
-(gradu_z(i ,j ,k-1)+gradu_z(i,j,k))*Kfaces(i,j,k)%nz*Kfaces(i,j,k)%A &
+(gradu_z(i+1,j ,k )+gradu_z(i,j,k))*Ifaces(i+1,j ,k )%nz*Ifaces(i+1,j ,k )%A &
+(gradu_z(i ,j+1,k )+gradu_z(i,j,k))*Jfaces(i ,j+1,k )%nz*Jfaces(i ,j+1,k )%A &
+(gradu_z(i ,j ,k+1)+gradu_z(i,j,k))*Kfaces(i ,j ,k+1)%nz*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
! gradient of v component
d2vdx2 =(-(gradv_x(i-1,j ,k )+gradv_x(i,j,k))*Ifaces(i,j,k)%nx*Ifaces(i,j,k)%A &
-(gradv_x(i ,j-1,k )+gradv_x(i,j,k))*Jfaces(i,j,k)%nx*Jfaces(i,j,k)%A &
-(gradv_x(i ,j ,k-1)+gradv_x(i,j,k))*Kfaces(i,j,k)%nx*Kfaces(i,j,k)%A &
+(gradv_x(i+1,j ,k )+gradv_x(i,j,k))*Ifaces(i+1,j ,k )%nx*Ifaces(i+1,j ,k )%A &
+(gradv_x(i ,j+1,k )+gradv_x(i,j,k))*Jfaces(i ,j+1,k )%nx*Jfaces(i ,j+1,k )%A &
+(gradv_x(i ,j ,k+1)+gradv_x(i,j,k))*Kfaces(i ,j ,k+1)%nx*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
d2vdy2 =(-(gradv_y(i-1,j ,k )+gradv_y(i,j,k))*Ifaces(i,j,k)%ny*Ifaces(i,j,k)%A &
-(gradv_y(i ,j-1,k )+gradv_y(i,j,k))*Jfaces(i,j,k)%ny*Jfaces(i,j,k)%A &
-(gradv_y(i ,j ,k-1)+gradv_y(i,j,k))*Kfaces(i,j,k)%ny*Kfaces(i,j,k)%A &
+(gradv_y(i+1,j ,k )+gradv_y(i,j,k))*Ifaces(i+1,j ,k )%ny*Ifaces(i+1,j ,k )%A &
+(gradv_y(i ,j+1,k )+gradv_y(i,j,k))*Jfaces(i ,j+1,k )%ny*Jfaces(i ,j+1,k )%A &
+(gradv_y(i ,j ,k+1)+gradv_y(i,j,k))*Kfaces(i ,j ,k+1)%ny*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
d2vdz2 =(-(gradv_z(i-1,j ,k )+gradv_z(i,j,k))*Ifaces(i,j,k)%nz*Ifaces(i,j,k)%A &
-(gradv_z(i ,j-1,k )+gradv_z(i,j,k))*Jfaces(i,j,k)%nz*Jfaces(i,j,k)%A &
-(gradv_z(i ,j ,k-1)+gradv_z(i,j,k))*Kfaces(i,j,k)%nz*Kfaces(i,j,k)%A &
+(gradv_z(i+1,j ,k )+gradv_z(i,j,k))*Ifaces(i+1,j ,k )%nz*Ifaces(i+1,j ,k )%A &
+(gradv_z(i ,j+1,k )+gradv_z(i,j,k))*Jfaces(i ,j+1,k )%nz*Jfaces(i ,j+1,k )%A &
+(gradv_z(i ,j ,k+1)+gradv_z(i,j,k))*Kfaces(i ,j ,k+1)%nz*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
!gradients of w components
d2wdx2 =(-(gradw_x(i-1,j ,k )+gradw_x(i,j,k))*Ifaces(i,j,k)%nx*Ifaces(i,j,k)%A &
-(gradw_x(i ,j-1,k )+gradw_x(i,j,k))*Jfaces(i,j,k)%nx*Jfaces(i,j,k)%A &
-(gradw_x(i ,j ,k-1)+gradw_x(i,j,k))*Kfaces(i,j,k)%nx*Kfaces(i,j,k)%A &
+(gradw_x(i+1,j ,k )+gradw_x(i,j,k))*Ifaces(i+1,j ,k )%nx*Ifaces(i+1,j ,k )%A &
+(gradw_x(i ,j+1,k )+gradw_x(i,j,k))*Jfaces(i ,j+1,k )%nx*Jfaces(i ,j+1,k )%A &
+(gradw_x(i ,j ,k+1)+gradw_x(i,j,k))*Kfaces(i ,j ,k+1)%nx*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
d2wdy2 =(-(gradw_y(i-1,j ,k )+gradw_y(i,j,k))*Ifaces(i,j,k)%ny*Ifaces(i,j,k)%A &
-(gradw_y(i ,j-1,k )+gradw_y(i,j,k))*Jfaces(i,j,k)%ny*Jfaces(i,j,k)%A &
-(gradw_y(i ,j ,k-1)+gradw_y(i,j,k))*Kfaces(i,j,k)%ny*Kfaces(i,j,k)%A &
+(gradw_y(i+1,j ,k )+gradw_y(i,j,k))*Ifaces(i+1,j ,k )%ny*Ifaces(i+1,j ,k )%A &
+(gradw_y(i ,j+1,k )+gradw_y(i,j,k))*Jfaces(i ,j+1,k )%ny*Jfaces(i ,j+1,k )%A &
+(gradw_y(i ,j ,k+1)+gradw_y(i,j,k))*Kfaces(i ,j ,k+1)%ny*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
d2wdz2 =(-(gradw_z(i-1,j ,k )+gradw_z(i,j,k))*Ifaces(i,j,k)%nz*Ifaces(i,j,k)%A &
-(gradw_z(i ,j-1,k )+gradw_z(i,j,k))*Jfaces(i,j,k)%nz*Jfaces(i,j,k)%A &
-(gradw_z(i ,j ,k-1)+gradw_z(i,j,k))*Kfaces(i,j,k)%nz*Kfaces(i,j,k)%A &
+(gradw_z(i+1,j ,k )+gradw_z(i,j,k))*Ifaces(i+1,j ,k )%nz*Ifaces(i+1,j ,k )%A &
+(gradw_z(i ,j+1,k )+gradw_z(i,j,k))*Jfaces(i ,j+1,k )%nz*Jfaces(i ,j+1,k )%A &
+(gradw_z(i ,j ,k+1)+gradw_z(i,j,k))*Kfaces(i ,j ,k+1)%nz*Kfaces(i ,j ,k+1)%A &
)/(2*cells(i,j,k)%volume)
udd = sqrt( (d2udx2+d2udy2+d2udz2)**2 &
+ (d2vdx2+d2vdy2+d2vdz2)**2 &
+ (d2wdx2+d2wdy2+d2wdz2)**2 )
ud = sqrt(2*(s11**2 + s12**2 + s13**2 &
+s21**2 + s22**2 + s23**2 &
+s31**2 + s32**2 + s33**2 ))
Lvk = kappa*abs(ud/max(udd,1.e-20))
fp = min(max(P_k/D_k, 0.5),1.0)
! Lvk limiter
Lvk = max(Lvk, tkl/max((tk*c11),1.e-20))
Lvk = min(Lvk, c12*kappa*dist(i,j,k)*fp)
eta = density*dist(i,j,k)*sqrt(0.3*tk)/(20*mu(i,j,k))
fphi = (1 + cd1*eta)/(1 + eta**4)
cphi2 = zeta3
cphi1 = (zeta1 - zeta2*((tkl/max(tk*Lvk,1.e-20))**2))
P_kl = cphi1*tkl*P_k/max(tk,1.e-20)
D_kl = cphi2*density*(tk**1.5)
S_k = P_k - D_k - 2*mu(i,j,k)*tk/(dist(i,j,k)**2) !Source term TKE
S_kl = P_kl - D_kl - 6*mu(i,j,k)*tkl*fphi/(dist(i,j,k)**2) !source term KL
S_k = S_k * cells(i, j, k)%volume
S_kl = S_kl * cells(i, j, k)%volume
residue(i, j, k, 6) = residue(i, j, k, 6) - S_k
residue(i, j, k, 7) = residue(i, j, k, 7) - S_kl
end do
end do
end do
end subroutine add_kkl_source