module asbm_m
! This module incorporates into the main code the ASBM algorithm by calling the
! calc_struct_tensors_at_node subroutine. Additionally, a smoothing expression is applied
! to the Reynolds stress and mean velocity tensors respectively, in order to increase the
! stability of the algebraic model. In the specific code, the smoothing expression assumes
! uniform grids, however more sophisticated expressions can be used.
! The following expressions are expressed as given in the following reference paper
!=========================================================================================
! The ASBM-SA turbulence closure: Taking advantage of structure-based
! modeling in current engineering CFD codes
!
! Panagiotou,C. and Kassinos,S.
! 2015
!
! Journal= International Journal of Heat Fluid Flow,
volume=52,
pages=111-128
!=========================================================================================
!define variables
================================================================================
!set number of times that smoothing is applied. For the 2D cases, 5 iterations
!are believed to be sufficient enough to provide smooth, fully converged solutions.
! integers
max_filters_Gij ! number of iterations where smoothing is applied on the mean velocity gradient tensor Gij
max_filters_Rij ! number of iterations where smoothing is applied on the Reynolds stress tensor Rij
ktrmax ! maximum iterations of newton-raphson method, used as input in the
! ASBM subroutine
!logicals
bound ! used for refinement reasons, always set to FALSE
fix ! used for refinement reasons, always set to FALSE
FILTERING ! used to determine if smoothing is needed
!Scalars
etaR ! mean rotation over strain parameter, used as output in the ASBM subroutine,
! defined in eq.11a
etaF ! frame rotation over strain parameter, used as output in the ASBM subroutine,
! defined in eq.11b
a2 ! invariant of the eddy-axis tensor a_ij x a_ji used as output in the ASBM subroutine,
! defined in eq.11c
bphi ! wall blockage scalar
gradB2 ! squared magnitude of the wall blockage gradient, appeared as denominator in eq.13
nu ! kinematic viscosity
k_SPA ! kinetic energy, needed to dimensionalize the reynolds streeses (Rij=2k x rij)
ts ! turbulent time scale
omega ! frame rotation vector
! Vectors
AxisEddyS_xj ! irrotational part of a_xj, defined in eq.5
AxisEddyS_yj ! irrotational part of a_yj, defined in eq.5
AxisEddyS_zj ! irrotational part of a_zj, defined in eq.5
AxisEddyH_xj ! homogeneous part of a_xj, defined in eq.6
AxisEddyH_yj ! homogeneous part of a_yj, defined in eq.6
AxisEddyH_zj ! homogeneous part of a_zj, defined in eq.6
AxisEddyB_xj ! inhomogeneous (final) eddy axis tensor a_xj, defined in eq.14
AxisEddyB_yj ! inhomogeneous (final) eddy axis tensor a_yj, defined in eq.14
AxisEddyB_zj ! inhomogeneous (final) eddy axis tensor a_zj, defined in eq.14
Hblock_xj ! partial projection operator P_xj, defined in eq.15
Hblock_yj ! partial projection operator P_yj, defined in eq.15
Hblock_zj ! partial projection operator P_zj, defined in eq.15
Bblock_xj ! wall blocking tensor B_xj, defined in eq.13
Bblock_yj ! wall blocking tensor B_yj, defined in eq.13
Bblock_zj ! wall blocking tensor B_zj, defined in eq.13
Rx !Reynolds stress components involving -x direction R_xj ( Rxy, Rxz,Rzx, Rxy,Rxx)
Ry !Reynolds stress components involving -y direction R_yj ( Ryx, Ryz,Rzy, Rxy,Ryy)
Rz !Reynolds stress components involving -z direction R_zj ( Rzx, Rzy,Rxz, Ryz,Rzz)
Dx !same as before, but for the Dimensionality components D_xj
Dy !same as before, but for the Dimensionality components D_yj
Dz !same as before, but for the Dimensionality components D_zj
Fx !same as before, but for the circulicity components F_xj
Fy !same as before, but for the circulicity components F_yj
Fz !same as before, but for the circulicity components F_zj
gradB ! gradient of wall blockage scalar (phi), defined as d_phi/d x_j
omega ! frame rotation vector
!Tensors
st ! strain rate tensor S_ij
wt ! Rotation rate tensor W_ij
lb ! wall blocking tensor B_ij, used as input in the ASBM subroutine
las ! irrotational part of a_ij, used as input/out in the ASBM subroutine
lar ! homogeneous part of a_ij, used as input/output in the ASBM subroutine
la ! final eddy-axis tensor a_ij, used as output in the ASBM subroutine
lwt ! nondimensional mean rotation rate tensor, used as input in the ASBM subroutine
lst ! nondimensional strain rate tensor, used as input in the ASBM subroutine
lwft ! non-dimensional frame rotation tensor, used as input in the ASBM subroutine
lHblock ! partial projection tensor, used as output in the ASBM subroutine
lR, lD, lF ! normalized structure tensors (r_ij, d_ij, f_ij), used as outputs
! in the ASBM subroutine
!end of variable assignment
!======================================================================================
Next, for simplicity we consider the whole process at one node, since the generalization
to the whole domain is straight forward.
!In order to improve the stability of the algorithm,
!ASBM takes the previous time step values for the Newton-Rhapson iteration during the
!entire computation.
!At the first time that asbm is called, we consider nearly isotropic turbulence
AxisEddyS_xj(1) = 0.3
AxisEddyS_xj(2) = 0.0
AxisEddyS_xj(3) = 0.0
AxisEddyS_yj(1) = 0.0
AxisEddyS_yj(2) = 0.3
AxisEddyS_yj(3) = 0.0
AxisEddyS_zj(1) = 0.0
AxisEddyS_zj(2) = 0.0
AxisEddyS_zj(3) = 0.4
AxisEddyH_xj(1) = 0.3
AxisEddyH_xj(2) = 0.0
AxisEddyH_xj(3) = 0.0
AxisEddyH_yj(1) = 0.0
AxisEddyH_yj(2) = 0.3
AxisEddyH_yj(3) = 0.0
AxisEddyH_zj(1) = 0.0
AxisEddyH_zj(2) = 0.0
AxisEddyH_zj(3) = 0.4
!If a previous ASBM solution is used as initial guess for the current simulation, then the
! previous solution is chosen as initial guess.
'Initial guess for a_ij = old values from restart file'
!Accordingly, SMOOTHING is applied on the dimensionless Strain rate and
!rotation rate tensors.
do i = 1,max_filters_Gij
SMOOTHING FUNCTION IS CALLED AND APPLIED ON S_ij and W_ij
end do
Set blockage scalar bphi bounds from 0 to 1.
Define the Blockage tensor B_ij, as given in eq.13.
if gradB2 less than a very small value then
Bblock_xj(1:3) = 0 ! for the cases considered in the paper,
Bblock_yj(1:3) = 0 ! setting the right hand side equal to zero did not
Bblock_zj(1:3) = 0 ! change the results.
else
Bblock_xj(1:3) = bphi(ino)*gradB(1)*gradB(1:3)/gradB2
Bblock_yj(1:3) = bphi(ino)*gradB(2)*gradB(1:3)/gradB2
Bblock_zj(1:3) = bphi(ino)*gradB(3)*gradB(1:3)/gradB2
end if
!For convenience, we set the wall blockage tensor lb, which will be used as input
! in the asbm subroutine
lb(1,1:3) = Bblock_xj(1:3)
lb(2,1:3) = Bblock_yj(1:3)
lb(3,1:3) = Bblock_zj(1:3)
!account for the frame rotation effects.
!define the non-dimensional frame rotation tensor
lwft(1,1) = 0.0
lwft(2,2) = 0.0
lwft(3,3) = 0.0
lwft(1,2) = -omega(3)*ts !Use the expression W_qp= e_pqi x omega_i x ts
lwft(2,1) = omega(3)*ts
lwft(2,3) = -omega(1)*ts
lwft(3,2) = omega(1)*ts
lwft(3,1) = -omega(2)*ts
lwft(1,3) = omega(2)*ts
! non-dimensionalize the strain rate and rotation rate tensors
! in respect to the time scale
lwt = wt(:,:,ino)*ts
lst = st(:,:,ino)*ts
! The relevant eddy-axis tensors are chosen to take the values of the previous time
!step solution in order to provide a better initial guess and accelerate the convergence
! of the iterative methods.
las(1,1:3) = AxisEddyS_xj(:)
las(2,1:3) = AxisEddyS_yj(:)
las(3,1:3) = AxisEddyS_zj(:)
lar(1,1:3) = AxisEddyH_xj(:)
lar(2,1:3) = AxisEddyH_yj(:)
lar(3,1:3) = AxisEddyH_zj(:)
!Set zero the structure tensors, which are used as inputs/outputs in the asbm
! subroutines.
lr = 0.0; lf = 0.0; ld = 0.0
ierr = 0
! call the ASBM subroutine, which calculates the relevant structure tensors
! at the corresponding node.
call calc_struct_tensors_at_node(lst, lwt, lwft, lb, las, lar, la, lHblock, &
lr, ld, lf, lktrmax, ierr, etaF(ino), etaR(ino), a2(ino), bound, fix)
!update the turbulent quantities, needed for convergence purposes at the next time step
AxisEddyS_xj(:) = las(1,:)
AxisEddyS_yj(:) = las(2,:)
AxisEddyS_zj(:) = las(3,:)
AxisEddyH_xj(:) = lar(1,:)
AxisEddyH_yj(:) = lar(2,:)
AxisEddyH_zj(:) = lar(3,:)
AxisEddyB_xj(:) = la(1,:)
AxisEddyB_yj(:) = la(2,:)
AxisEddyB_zj(:) = la(3,:)
Hblock_xj(:) = lHblock(1,:)
Hblock_yj(:) = lHblock(2,:)
Hblock_zj(:) = lHblock(3,:)
!dimensionalize the stresses
Rx(:,ino) = lr(1,:)*2.0*k_SPA
Ry(:,ino) = lr(2,:)*2.0*k_SPA
Rz(:,ino) = lr(3,:)*2.0*k_SPA
!for two dimensional flows (with the axis of independence parallel to -z direction),
!the relevant stress components are set zero in order to improve the stability
! of the algorithm.
if (two dimensional flow) then
Rx(3) = 0.0 !Rxz
Ry(3) = 0.0 !Ryz
Rz(1) = 0.0 !Rzx
Rz(2) = 0.0 !Rzy
end if
!dimensionalize the remaining structure tensors.
Dx(:,ino) = ld(1,:)*2.0*k_SPA
Dy(:,ino) = ld(2,:)*2.0*k_SPA
Dz(:,ino) = ld(3,:)*2.0*k_SPA
Fx(:,ino) = lf(1,:)*2.0*k_SPA
Fy(:,ino) = lf(2,:)*2.0*k_SPA
Fz(:,ino) = lf(3,:)*2.0*k_SPA
end do
!!!!!!!!!!!!!!!!!!!!!!!1 end of fotos routine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!applying smoothing method on the updated Reynolds stresses which are used in the momentum
!equations.
!Accordingly, SMOOTHING is applied on the Reynolds stress tensor.
do i = 1,max_filters_Gij
SMOOTHING FUNCTION IS CALLED AND APPLIED ON R_ij
end do
end module