Add residual due to source terms of the LCTM2015 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 | ||
| 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(schemetype), | intent(in) | :: | scheme | finite-volume Schemes |
||
| type(extent), | intent(in) | :: | dims | Extent of the domain:imx,jmx,kmx |
subroutine add_sst_source_lctm2015(qp, residue, cells, Ifaces, Jfaces, Kfaces, scheme, dims)
!< Add residual due to source terms of the LCTM2015 transition model
implicit none
type(schemetype), intent(in) :: scheme
!< finite-volume Schemes
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) :: CD
!< Cross-diffustion term
real(wp) :: F1
!< single cell blending function
real(wp) :: vort
!< vorticity magnitude
real(wp) :: S_K
!< Total source term in TKE equation
real(wp) :: S_w
!< Total source term in Omega equation
real(wp) :: S_gm
!< Total source term in Gamma equation
real(wp) :: D_k
!< Destruction term in TKE equation
real(wp) :: D_w
!< Destruction term in Omega equation
real(wp) :: D_gm
!< Destruction term in Gamma equation
real(wp) :: P_k
!< production term in TKE equation
real(wp) :: P_w
!< production term in Omega equation
real(wp) :: P_gm
!< production term in Gamma equation
real(wp) :: lamda
!< additional source term in Omega equation
real(wp) :: Fonset1,Fonset2, Fonset3,Fonset
!< Transition onset term
real(wp) :: Rev
!< Reynodlds number based on vorticity
real(wp) :: RT
!< Turbulent reynolds number
real(wp) :: Fturb
real(wp) :: Re_theta
!< Cutt-off reynolds number based on momentum thickness
real(wp) :: TuL
!< local turbulence intensity
real(wp) :: strain
!< Strain rate magnitude
real(wp) :: intermittency
!< intermittency
real(wp) :: Pk_lim
!< production lim term
real(wp) :: Fon_lim
real(wp) :: lamd
!< pressure gradient
real(wp) :: Fpg
!< pressure gradient functin
real(wp) :: divergence
!< del.V
real(wp) :: dvdy
!< pressure gradient sensor
integer :: limiter
!< production limiter
real(wp) :: density
!< single cell Density
real(wp) :: tk
!< single cell TKE
real(wp) :: tw
!< single cell Omega
if(trim(scheme%turbulence) == 'sst2003')then
limiter = 10
gama1 = 5.0/9.0
gama2 = 0.44
else
limiter = 20
end if
!for pressure gradient calculation
call find_DCCVn(qp, cells, Ifaces, Jfaces, Kfaces, dims)
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)
tw = qp(i,j,k,7)
intermittency = qp(i,j,k,8)
! __ 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 &
)&
)
strain = 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) &
+ 2*(gradu_x(i,j,k)**2) &
+ 2*(gradv_y(i,j,k)**2) &
+ 2*(gradw_z(i,j,k)**2) &
)&
)
CD = 2*density*sigma_w2*(gradtk_x(i,j,k)*gradtw_x(i,j,k)&
+ gradtk_y(i,j,k)*gradtw_y(i,j,k)&
+ gradtk_z(i,j,k)*gradtw_z(i,j,k)&
)/tw
CD = max(CD, 10.0**(-limiter))
F1 = sst_F1(i,j,k)
sigma_k = sigma_k1*F1 + sigma_k2*(1. - F1)
sigma_w = sigma_w1*F1 + sigma_w1*(1. - F1)
gama = gama1*F1 + gama2*(1. - F1)
beta = beta1*F1 + beta2*(1. - F1)
! ____ Dissipation term ___
D_k = bstar*density*tw*tk
D_w = beta*density*tw**2
! ____ PRODUCTION term____
divergence = gradu_x(i,j,k) + gradv_y(i,j,k) + gradw_z(i,j,k)
P_k = mu_t(i,j,k)*(vort*strain) - ((2.0/3.0)*density*tk*divergence)
P_k = min(P_k, limiter*D_k)
P_w = (density*gama/mu_t(i,j,k))*P_k
! ____ cross diffusion term ___
lamda = (1. - F1)*CD
! ____Transition modeling ____
! --pressure gradient
dvdy = DCCVnX(i,j,k)*CCnormalX(i,j,k) &
+ DCCVnY(i,j,k)*CCnormalY(i,j,k) &
+ DCCVnZ(i,j,k)*CCnormalZ(i,j,k)
lamd =(-7.57e-3)*(dvdy*dist(i,j,k)*dist(i,j,k)*density/mu(i,j,k)) + 0.0128
lamd = min(max(lamd, -1.0), 1.0)
if(lamd>=0.0)then
Fpg = min(1.0 + 14.68*lamd, 1.5)
else
Fpg = min(1.0 - 7.34*lamd, 3.0)
end if
Fpg = max(Fpg, 0.0)
! --gradient
TuL = min(100.0*sqrt(2.0*tk/3.0)/(tw*dist(i,j,k)),100.0)
Re_theta = 100.0 + 1000.0*exp(-TuL*Fpg)
!Re_theta = 100.0 + 1000.0*exp(-TuL)
Rev = density*dist(i,j,k)*dist(i,j,k)*strain/mu(i,j,k)
RT = density*tk/(mu(i,j,k)*tw)
Fturb = exp(-(0.5*Rt)**4)
Fonset1 = Rev/(2.2*Re_theta)
Fonset2 = min(Fonset1, 2.0)
Fonset3 = max(1.0 - (RT/3.5)**3, 0.0)
Fonset = max(Fonset2 - Fonset3, 0.0)
P_gm = 100*density*strain*intermittency*(1.0 - intermittency)*Fonset
D_gm = 0.06*density*vort*intermittency*Fturb*((50.0*intermittency) - 1.0)
Fon_lim = min(max((Rev/(2.2*1100.0))-1.0, 0.0), 3.0)
Pk_lim = 5*max(intermittency - 0.2,0.0)*(1.0 - intermittency)*Fon_lim*max(3*mu(i,j,k) - mu_t(i,j,k), 0.0)*strain*vort
S_k = intermittency*P_k - max(intermittency,0.1)*D_k + Pk_lim !Source term gm
S_W = P_w - D_w + lamda !Source term gm
S_gm = P_gm - D_gm !Source term gm
S_k = S_k * cells(i, j, k)%volume
S_w = S_w * cells(i, j, k)%volume
S_gm= S_gm* 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_w
residue(i, j, k, 8) = residue(i, j, k, 8) - S_gm
end do
end do
end do
end subroutine add_sst_source_lctm2015