! Copyright (c) 2022 Kannan Masilamani <kannan.masilaman@dlr.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! **************************************************************************** ! !> This module contains turbulent wall function type and routines to calculate !! friction velocity and stream-wise velocity component. module mus_turb_wallFunc_module ! include treelm modules use env_module, only: rk, long_k, labelLen use tem_aux_module, only: tem_abort use tem_logging_module, only: logUnit use tem_float_module, only: operator(.fle.) ! include aotus modules use aotus_module, only: flu_State, aoterr_Fatal, aoterr_NonExistent, & & aoterr_WrongType, aot_get_val implicit none private public :: mus_turb_wallFunc_type public :: mus_load_turb_wallFunc !> Contains friction velocity and turbulent viscosity on boundary elements on !! each level type mus_turb_wallFunc_data_type !> Turbulent viscosity on boundary element computed using mixing length !! formulation in lattice unit !! nu_t = (vonKarman*distToBnd)**2 * |du/dy| real(kind=rk), allocatable :: tVisc(:) !> Friction velocity on first neighbor in normal direction in lattice unit !! computed from wall model real(kind=rk), allocatable :: velTau(:) !> Distance to boundary from first fluid in normal direction !! in lattice unit. real(kind=rk), allocatable :: distToBnd(:) !> Distance to boundary from first fluid neighbor in normal direction !! in lattice unit. real(kind=rk), allocatable :: neighDistToBnd(:) end type mus_turb_wallFunc_data_type !> Contains function pointers to compute fuction velocity and stream-wise !! velocity component type mus_turb_wallFunc_type !> is true if wall function is active logical :: isActive = .false. !> Wall model function character(len=labelLen) :: wall_func !> Nonlinear solver type character(len=labelLen) :: nonlinear_solver !> Von-Karman constant. Default = 0.4_rk real(kind=rk) :: vonKarman = 0.4_rk !> Use vanDriest damping function to damp turbulent viscosity logical :: useVanDriest = .true. !> Contains data computed in turbulent wall bc routine on each level type(mus_turb_wallFunc_data_type), allocatable :: dataOnLvl(:) !> Function pointer to compute friction velocity procedure(mus_proc_calcFricVel), pointer, nopass :: calcFricVel => null() !> Function pointer to compute strean-wise velocity component procedure(mus_proc_calcStreamWiseVel), pointer, nopass :: & & calcStreamWiseVel => null() end type mus_turb_wallFunc_type !> Interface definition for the turbulent wall bc routines abstract interface !> This abstract interface defines the interface to calculate turbulent wall !! friction velocity from given velocity and distance to boundary. !! All inputs and output are in lattice units. subroutine mus_proc_calcFricVel(velTau, velSW, distToBnd, & & viscKine, nElems) import :: rk !> Friction velocity computed from wall model. !! it is inout to provide velTau from previous timestep as initial velTau !! for fixed-point or Newton iteration solver real(kind=rk), intent(inout) :: velTau(:) !> Stream-wise velocity component from which friction velocity is computed real(kind=rk), intent(in) :: velSW(:) !> Distance to the boundary in the discrete normal direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems end subroutine mus_proc_calcFricVel !> This abstract interface defines the interface to calculate stream-wise !! velocity component from friction velocity and distance to boundary. !! All inputs and output are in lattice units. subroutine mus_proc_calcStreamWiseVel(velSW, velTau, distToBnd, & & viscKine, nElems) import :: rk !> Stream-wise velocity component from wall model real(kind=rk), intent(out) :: velSW(:) !> Friction velocity computd from wall model real(kind=rk), intent(in) :: velTau(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems end subroutine mus_proc_calcStreamWiseVel end interface !! Constant parameters for Implicit equation solver integer, parameter :: imEq_nIter = 1000 real(kind=rk), parameter :: imEq_tol = 1e-10 !! Constant parameters for Reichardt's law real(kind=rk), parameter :: reivonKA = 0.4_rk real(kind=rk), parameter :: oneOvervonKA = 1.0_rk/reivonKA real(kind=rk), parameter :: reiC = 7.8_rk real(kind=rk), parameter :: reiB1 = 11.0_rk real(kind=rk), parameter :: reiB2 = 3.0_rk !! constants for Werner & Wengle model !! http://trace-portal.de/userguide/trace/page_wallFunctionWernerWengle.html real(kind=rk), parameter :: ww_C_m = 8.3_rk real(kind=rk), parameter :: ww_m = 1.0_rk/7.0_rk real(kind=rk), parameter :: ww_uLmt_yPlus = ww_C_m**(1.0_rk/(1.0_rk-ww_m)) real(kind=rk), parameter :: ww_uLmt_vel = WW_uLmt_yPlus**2.0_rk real(kind=rk), parameter :: ww_onePlusM = 1.0_rk+ww_m real(kind=rk), parameter :: ww_onePlusMInv = 1.0_rk/ww_onePlusM real(kind=rk), parameter :: ww_mDivOnePlusM = ww_m/ww_onePlusM !! Werner & Wengle two layer model does not cover buffer layer so !! three different layer model from Schmidt is used to compute stream-wise !! velocity at first element. !! Viscous sublayer < 5 !! buffer layer > 5 and < 30 !! log layer > 30 real(kind=rk), parameter :: sc_lLmt = 5.0_rk ! lower limit of buffer layer real(kind=rk), parameter :: sc_uLmt = 30.0_rk ! upper limit for buffer layer real(kind=rk), parameter :: sc_vonKarman = 0.4_rk real(kind=rk), parameter :: sc_a = ((1.0_rk/sc_vonKarman)*log(30.0_rk) & & + 0.2_rk)/log(6.0_rk) real(kind=rk), parameter :: sc_b = 5.0_rk-sc_a*log(5.0_rk) contains ! -------------------------------------------------------------------------- ! !> This routine loads wall model and nonlinear solver type for nonlinear !! equation subroutine mus_load_turb_wallFunc(me, conf, parent) ! -------------------------------------------------------------------- ! !> Turbulent wall model type to fill assign wall model type(mus_turb_wallFunc_type), intent(inout) :: me !> lua flu state type( flu_state ) :: conf !> bc parent handle integer, intent(in) :: parent ! -------------------------------------------------------------------- ! integer :: iError ! -------------------------------------------------------------------- ! write(logUnit(10), "(A)") ' Loading turbulent wall_function ...' me%isActive = .true. ! wall model function name call aot_get_val(L = conf, & & thandle = parent, & & key = 'wall_function', & & val = me%wall_func, & & default = 'musker', & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*)'FATAL Error occured, while retrieving wall_function:' if (btest(iError, aoterr_WrongType)) then write(logUnit(1),*)'Variable has wrong type!' write(logUnit(1),*)'STOPPING' call tem_abort() end if end if write(logUnit(1),*) ' wall function: ', trim(me%wall_func) ! Get nonlinear solver for wall model other then power_law if (trim(me%wall_func) /= 'power_law') then call aot_get_val(L = conf, & & thandle = parent, & & key = 'nonlinear_solver', & & val = me%nonlinear_solver, & & default = 'fixed_point', & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*)'FATAL Error occured, while retrieving ' & & //'nonlinear_solver:' if (btest(iError, aoterr_WrongType)) then write(logUnit(1),*)'Variable has wrong type!' write(logUnit(1),*)'STOPPING' call tem_abort() end if end if write(logUnit(1),*) ' nonlinear solver: ', trim(me%nonlinear_solver) end if ! Von-Karman constant call aot_get_val(L = conf, & & thandle = parent, & & key = 'von_karman', & & val = me%vonKarman, & & default = 0.4_rk, & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*)'FATAL Error occured, while retrieving ' & & //'von_karman:' if (btest(iError, aoterr_WrongType)) then write(logUnit(1),*)'Variable has wrong type!' write(logUnit(1),*)'STOPPING' call tem_abort() end if end if write(logUnit(1),*) ' von Karman: ', me%vonKarman ! Use van Driest damping function call aot_get_val(L = conf, & & thandle = parent, & & key = 'use_van_driest', & & val = me%useVanDriest, & & default = .true., & & ErrCode = iError ) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*)'FATAL Error occured, while retrieving ' & & //'use_van_driest:' if (btest(iError, aoterr_WrongType)) then write(logUnit(1),*)'Variable has wrong type!' write(logUnit(1),*)'STOPPING' call tem_abort() end if end if write(logUnit(1),*) ' use van Driest: ', me%useVanDriest ! Assign turbulent wall model function pointer for friction and streamwise ! velocity calculation call assign_turb_wallFunc(me) end subroutine mus_load_turb_wallFunc ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> Assign wall models pointers to calculate friction velocity and stream-wise !! velocity component subroutine assign_turb_wallFunc(me) ! -------------------------------------------------------------------- ! !> Turbulent wall model type to fill assign wall model type(mus_turb_wallFunc_type), intent(inout) :: me ! -------------------------------------------------------------------- ! ! -------------------------------------------------------------------- ! select case (trim(me%wall_func)) case ('musker') me%calcStreamWiseVel => streamWiseVel_musker select case (trim(me%nonlinear_solver)) case ('fixed_point') me%calcFricVel => fricVel_musker_fixedPoint case ('newton') me%calcFricVel => fricVel_musker_newton case default call tem_abort('Error: Unsupported nonlinear solver for wall function') end select case ('reichardt') me%calcStreamWiseVel => streamWiseVel_reichardt select case (trim(me%nonlinear_solver)) case ('fixed_point') me%calcFricVel => fricVel_reichardt_fixedPoint case ('newton') me%calcFricVel => fricVel_reichardt_newton case default call tem_abort('Error: Unsupported nonlinear solver for wall function') end select case ('power_law') me%calcStreamWiseVel => streamWiseVel_powerLaw me%calcFricVel => fricVel_powerLaw case default call tem_abort("Error: Unknown turbulent wall function") end select end subroutine assign_turb_wallFunc ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine computes friction velocity from Musker wall model profile !! using fixed-point iteration method subroutine fricVel_musker_fixedPoint(velTau, velSW, distToBnd, viscKine, & & nElems) ! -------------------------------------------------------------------- ! !> Friction velocity computed from wall model. !! it is inout to provide velTau from previous timestep as initial velTau !! for fixed-point or Newton iteration solver real(kind=rk), intent(inout) :: velTau(:) !> Stream-wise velocity component from which friction velocity is computed real(kind=rk), intent(in) :: velSW(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! real(kind=rk) :: error integer :: iElem, iter real(kind=rk) :: velTau_old, velTau_new, yPlus, yPlus_fac ! -------------------------------------------------------------------- ! do iElem = 1, nElems ! Fixed point iteration iter = 0 velTau_old = velTau(iElem) yPlus_fac = distToBnd(iElem) / viscKine(iElem) fixed_point: do yPlus = velTau_old * yPlus_fac velTau_new = velSW(iElem) / ( 5.424_rk & & * atan((2.0_rk*yPlus - 8.15_rk)/16.7_rk) & & + log10( (yPlus + 10.6_rk)**9.6_rk & & / ( yPlus**2.0_rk - 8.15_rk*yPlus & & + 86.0_rk)**2.0_rk ) - 3.5072790194_rk ) error = abs(velTau_new - velTau_old) velTau_old = velTau_new iter = iter+1 if ((error .fle. imEq_tol) & & .or. (iter >= imEq_nIter)) exit fixed_point end do fixed_point velTau(iElem) = velTau_new end do end subroutine fricVel_musker_fixedPoint ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine computes friction velocity from Musker wall model profile !! using Newton iteration method subroutine fricVel_musker_newton(velTau, velSW, distToBnd, viscKine, & & nElems) ! -------------------------------------------------------------------- ! !> Friction velocity computed from wall model. !! it is inout to provide velTau from previous timestep as initial velTau !! for fixed-point or Newton iteration solver real(kind=rk), intent(inout) :: velTau(:) !> Stream-wise velocity component from which friction velocity is computed real(kind=rk), intent(in) :: velSW(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! real(kind=rk) :: error integer :: iElem, iter real(kind=rk) :: velTau_old, velTau_new, yPlus, yPlus_fac real(kind=rk) :: log_num, log_den, tan_num, fx, dfdx ! -------------------------------------------------------------------- ! do iElem = 1, nElems ! Fixed point iteration iter = 0 velTau_old = velTau(iElem) yPlus_fac = distToBnd(iElem) / viscKine(iElem) fixed_point: do yPlus = velTau_old * yPlus_fac velTau_new = velSW(iElem) / ( 5.424_rk & & * atan((2.0_rk*yPlus - 8.15_rk)/16.7_rk) & & + log10( (yPlus + 10.6_rk)**9.6_rk & & / ( yPlus**2.0_rk - 8.15_rk*yPlus & & + 86.0_rk)**2.0_rk ) - 3.5072790194_rk ) log_num = (yPlus + 10.6_rk)**9.6_rk log_den = (yPlus**2.0_rk - 8.15_rk * yPlus + 86.0_rk)**2.0_rk tan_num = 2.0_rk*yPlus - 8.15_rk ! Function of u_tau fx = velTau_old*(5.424_rk*atan(tan_num/16.7_rk) & & + log10(log_num/log_den) - 3.5072790194_rk)-velSW(iElem) ! derivative of function w.r.t u_tau dfdx = velTau_old * ( ( 9.6_rk*yPlus_fac*(yPlus+10.6_rk)**8.6_rk & & / log_den - (2_rk*log_num)*(2_rk*yPlus-8.15_rk*yPlus_Fac) & & / log_den**1.5_rk ) * ( log_den / (log(10.0_rk)*log_num) ) & & + 0.649581*yPlus_fac/(0.00358564_rk*(2.0_rk*yPlus-8.15)**2 + 1) ) & & + log(log_num/log_den)/log(10._rk) & & + 5.424_rk*atan(0.0598802_rk*(tan_num))-3.5072790194_rk velTau_new = velTau_old - (fx/dfdx) error = abs(velTau_new - velTau_old) velTau_old = velTau_new iter = iter+1 if ((error .fle. imEq_tol) & & .or. (iter >= imEq_nIter)) exit fixed_point end do fixed_point velTau(iElem) = velTau_new end do end subroutine fricVel_musker_newton ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routines computes streamWise velocity component from friction velocity !! and distance to boundary using Musker profile. subroutine streamWiseVel_musker(velSW, velTau, distToBnd, & & viscKine, nElems) ! -------------------------------------------------------------------- ! !> Stream-wise velocity component from wall model real(kind=rk), intent(out) :: velSW(:) !> Friction velocity computd from wall model real(kind=rk), intent(in) :: velTau(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! integer :: iElem real(kind=rk) :: yPlus ! -------------------------------------------------------------------- ! do iElem = 1, nElems yPlus = distToBnd(iElem) * velTau(iElem) / viscKine(iElem) velSW(iElem) = velTau(iElem) * ( 5.424_rk & & * atan((2.0_rk*yPlus - 8.15_rk)/16.7_rk) & & + log10( (yPlus + 10.6_rk)**9.6_rk & & / ( yPlus**2.0_rk - 8.15_rk*yPlus & & + 86.0_rk)**2.0_rk ) - 3.5072790194_rk ) end do end subroutine streamWiseVel_musker ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine computes friction velocity from reichardt wall model profile !! using fixed-point iteration method subroutine fricVel_reichardt_fixedPoint(velTau, velSW, distToBnd, viscKine, & & nElems) ! -------------------------------------------------------------------- ! !> Friction velocity computed from wall model. !! it is inout to provide velTau from previous timestep as initial velTau !! for fixed-point or Newton iteration solver real(kind=rk), intent(inout) :: velTau(:) !> Stream-wise velocity component from which friction velocity is computed real(kind=rk), intent(in) :: velSW(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! real(kind=rk) :: error integer :: iElem, iter real(kind=rk) :: velTau_old, velTau_new, yPlus, yPlus_fac ! -------------------------------------------------------------------- ! do iElem = 1, nElems ! Fixed point iteration iter = 0 velTau_old = velTau(iElem) yPlus_fac = distToBnd(iElem) / viscKine(iElem) fixed_point: do yPlus = velTau_old * yPlus_fac velTau_new = velSW(iElem) / ( oneOvervonKA & & * log( 1.0_rk + reivonKA*yPlus ) & & + reiC * ( 1.0_rk-exp(-yPlus/reiB1) & & - (yPlus/reiB1) * exp(-yPlus/reiB2) ) ) error = abs(velTau_new - velTau_old) velTau_old = velTau_new iter = iter+1 if ((error .fle. imEq_tol) & & .or. (iter >= imEq_nIter)) exit fixed_point end do fixed_point velTau(iElem) = velTau_new end do end subroutine fricVel_reichardt_fixedPoint ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine computes friction velocity from reichardt wall model profile !! using Newton iteration method subroutine fricVel_reichardt_newton(velTau, velSW, distToBnd, viscKine, & & nElems) ! -------------------------------------------------------------------- ! !> Friction velocity computed from wall model. !! it is inout to provide velTau from previous timestep as initial velTau !! for fixed-point or Newton iteration solver real(kind=rk), intent(inout) :: velTau(:) !> Stream-wise velocity component from which friction velocity is computed real(kind=rk), intent(in) :: velSW(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! ! -------------------------------------------------------------------- ! call tem_abort("Error: fricVel_reichardt_newton not implemented yet") end subroutine fricVel_reichardt_newton ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routines computes streamWise velocity component from friction velocity !! and distance to boundary using reichardt profile. subroutine streamWiseVel_reichardt(velSW, velTau, distToBnd, & & viscKine, nElems) ! -------------------------------------------------------------------- ! !> Stream-wise velocity component from wall model real(kind=rk), intent(out) :: velSW(:) !> Friction velocity computd from wall model real(kind=rk), intent(in) :: velTau(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! integer :: iElem real(kind=rk) :: yPlus ! -------------------------------------------------------------------- ! do iElem = 1, nElems yPlus = distToBnd(iElem) * velTau(iElem) / viscKine(iElem) velSW(iElem) = velTau(iElem) * ( oneOvervonKA & & * log( 1.0_rk + reivonKA*yPlus ) & & + reiC * ( 1.0_rk-exp(-yPlus/reiB1) & & - (yPlus/reiB1) * exp(-yPlus/reiB2) ) ) end do end subroutine streamWiseVel_reichardt ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routine computes friction velocity from explicit two equation !! Werner and Wengle wall model. !! Haussmann, Marc; BARRETO, Alejandro CLARO; KOUYI, Gislain LIPEME; !! Rivière, Nicolas; Nirschl, Hermann; Krause, Mathias J. (2019): !! Large-eddy simulation coupled with wall models for turbulent channel !! flows at high Reynolds numbers with a lattice Boltzmann method !! — Application to Coriolis mass flowmeter. In Computers & Mathematics !! with Applications 78 (10), pp. 3285–3302. DOI: 10.1016/j.camwa.2019.04.033. !! !! The explicit power-law in terms of friction velocity given in Eq. 33 in !! S. Wilhelm, J. Jacob, and P. Sagaut, "An explicit power-law-based wall !! model for lattice Boltzmann method–Reynolds-averaged numerical simulations !! of the flow around airfoils", Physics of Fluids 30, 065111 (2018) !! https://doi.org/10.1063/1.5031764 !! !! The model constants are chosen according to Werner and Wengle: !! Wengle, H. and Werner, H. (1993) ‘Large-eddy Simulation of Turbulent Flow !! Over Sharp-edged Obstacles in a Plate Channel’, (1985), pp. 192–199. subroutine fricVel_powerLaw(velTau, velSW, distToBnd, viscKine, & & nElems) ! -------------------------------------------------------------------- ! !> Friction velocity computed from wall model. !! it is inout to provide velTau from previous timestep as initial velTau !! for fixed-point or Newton iteration solver real(kind=rk), intent(inout) :: velTau(:) !> Stream-wise velocity component from which friction velocity is computed real(kind=rk), intent(in) :: velSW(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! integer :: iElem real(kind=rk) :: visc_div_dist, fac ! -------------------------------------------------------------------- ! do iElem = 1, nElems visc_div_dist = viscKine(iElem) / distToBnd(iElem) if (velSW(iElem) <= ww_uLmt_vel * visc_div_dist) then ! y_plus below 11.81 velTau(iElem) = sqrt(velSW(iElem) * visc_div_dist) else ! power-law fac = (visc_div_dist)**ww_mDivOnePlusM velTau(iElem) = fac * (velSW(iElem) / ww_C_m)**ww_onePlusMInv end if end do end subroutine fricVel_powerLaw ! -------------------------------------------------------------------------- ! ! -------------------------------------------------------------------------- ! !> This routines computes streamWise velocity component from friction velocity !! and distance to boundary using two equation Werner and Wengle profile. !! Ref to following paper for Schmitt three layer equations. !! Haussmann, M. et al. (2019) ‘Large-eddy simulation coupled with wall models !! for turbulent channel flows at high Reynolds numbers with a lattice !! Boltzmann method — Application to Coriolis mass flowmeter’, Computers & !! Mathematics with Applications. Elsevier Ltd, 78(10), pp. 3285–3302. subroutine streamWiseVel_powerLaw(velSW, velTau, distToBnd, & & viscKine, nElems) ! -------------------------------------------------------------------- ! !> Stream-wise velocity component from wall model real(kind=rk), intent(out) :: velSW(:) !> Friction velocity computd from wall model real(kind=rk), intent(in) :: velTau(:) !> Distance to the boundary in the discrete normai direction real(kind=rk), intent(in) :: distToBnd(:) !> Kinematic viscosity real(kind=rk), intent(in) :: viscKine(:) !> Number of elements in input and output arrays integer, intent(in) :: nElems ! -------------------------------------------------------------------- ! integer :: iElem real(kind=rk) :: yPlus ! -------------------------------------------------------------------- ! do iElem = 1, nElems yPlus = distToBnd(iElem) * velTau(iElem) / viscKine(iElem) if (yPlus >= sc_uLmt) then ! Interial layer, use powerlaw profile from werner and wengle velSW(iElem) = (ww_C_m * ( yPlus )**ww_m) * velTau(iElem) else if ( yPlus >= sc_lLmt .and. yPlus < sc_uLmt) then ! Buffer layer, use logarithmic profile velSW(iElem) = (sc_a*log(yPlus) + sc_b) * velTau(iElem) else ! viscous sublayer, use linear profile. velSW(iElem) = yPlus * velTau(iElem) end if end do end subroutine streamWiseVel_powerLaw ! -------------------------------------------------------------------------- ! end module mus_turb_wallFunc_module