! Copyright (c) 2021-2022 Kannan Masilamani <kannan.masilamani@dlr.de> ! Copyright (c) 2021 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@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. ! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de> ! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.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. ! Copyright (c) 2014-2015, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.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. ! **************************************************************************** ! !> Boundary condition wall treatment routines !! !! This module contains higher order wall treatments !! A detailed description on the implementation details are given !! in [[tem_bc_module]]. !! module mus_bc_fluid_turbulent_module use iso_c_binding, only: c_f_pointer ! include treelm modules use env_module, only: rk, long_k use tem_param_module, only: cs2inv, cs2, rho0, rho0Inv, cs4inv, & & div1_3 use tem_time_module, only: tem_time_type use treelmesh_module, only: treelmesh_type use tem_varSys_module, only: tem_varSys_type use tem_debug_module, only: dbgUnit use tem_geometry_module, only: tem_ElemSizeLevel use tem_property_module, only: prp_solid use tem_construction_module, only: tem_levelDesc_type use tem_float_module, only: operator(.feq.), operator(.fne.) use tem_stencil_module, only: tem_stencilHeader_type use tem_aux_module, only: tem_abort use tem_logging_module, only: logUnit ! include musubi modules use mus_bc_header_module, only: boundary_type, glob_boundary_type use mus_scheme_layout_module, only: mus_scheme_layout_type use mus_field_prop_module, only: mus_field_prop_type use mus_derVarPos_module, only: mus_derVarPos_type use mus_param_module, only: mus_param_type use mus_physics_module, only: mus_physics_type use mus_mixture_module, only: mus_mixture_type use mus_varSys_module, only: mus_varSys_data_type use mus_relaxationParam_module, only: mus_viscosity_type use mus_turb_wallFunc_module, only: mus_turb_wallFunc_type use mus_turbulence_module, only: mus_turbulence_type implicit none private public :: turbulent_wall public :: turbulent_wall_libb public :: turbulent_wall_eq public :: turbulent_wall_eq_curved public :: turbulent_wall_noneq_expol public :: turbulent_wall_noneq_expol_curved public :: turbulent_wall_eq_incomp public :: turbulent_wall_eq_curved_incomp public :: turbulent_wall_noneq_expol_incomp public :: turbulent_wall_noneq_expol_curved_incomp !! Constant parameters for van-driest damping function real(kind=rk), parameter :: vd_Aplus = 26.0_rk contains ! ************************************************************************** ! !> BC routine for turbulent wall. !! It uses wall model to compute velocity on the boundary node. !! The implementation is based on the following paper: !! Haussmann, Marc; Ries, Florian; Jeppener-Haltenhoff, Jonathan B.; Li, !! Yongxiang; Schmidt, Marius; Welch, Cooper et al. (2020): Evaluation of a !! Near-Wall-Modeled Large Eddy Lattice Boltzmann Method for the Analysis of !! Complex Flows Relevant to IC Engines. In Computation 8 (2), p. 43. !! DOI: 10.3390/computation8020043. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall', !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos, normalInd_inv real(kind=rk) :: f_pre(layout%fStencil%QQ) real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: dens, dens_inv, vel(3), velSW(3) real(kind=rk) :: dens_neigh, dens_neigh_inv, vel_neigh(3) real(kind=rk) :: dens_bnd, vel_normal, normal(3) real(kind=rk) :: fEq, fEq_bnd, fEq_neigh real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state. ! computeNeighBuf uses FETCH which does implicit bounce back which ! is valid for qVal=0.5 f_pre = me%neigh(iLevel)%computeNeighBuf( (iElem-1)*QQ+1: iElem*QQ ) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens = sum(f_pre) dens_inv = 1.0_rk/dens dens_neigh = sum(f_neigh) dens_neigh_inv = 1.0_rk/dens_neigh ! velocity vel(1) = sum( f_pre * layout%fstencil%cxdirrk(1, :) ) * dens_inv vel(2) = sum( f_pre * layout%fstencil%cxdirrk(2, :) ) * dens_inv vel(3) = sum( f_pre * layout%fstencil%cxdirrk(3, :) ) * dens_inv vel_neigh(1) = sum( f_neigh * layout%fstencil%cxdirrk(1, :) ) * dens_neigh_inv vel_neigh(2) = sum( f_neigh * layout%fstencil%cxdirrk(2, :) ) * dens_neigh_inv vel_neigh(3) = sum( f_neigh * layout%fstencil%cxdirrk(3, :) ) * dens_neigh_inv ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) ! calculate density at boundary for straight wall using Zho-He approach normalInd_inv = layout%fStencil%cxDirInv( globBC%elemLvl(iLevel) & & %normalInd%val(iElem) ) normal = layout%fStencil%cxDirRK(:, normalInd_inv) vel_normal = dot_product(velSW, normal) dens_bnd = 0.0_rk f_pre = me%neigh(iLevel)%computeNeighBuf( (iElem-1)*QQ+1: iElem*QQ ) do iDir = 1, layout%fStencil%QQ if ( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then dens_bnd = dens_bnd + f_pre(layout%fStencil%cxDirInv(iDir)) else dens_bnd = dens_bnd + f_pre(iDir) end if end do dens_bnd = dens_bnd / (1.0_rk + vel_normal) do iDir = 1, layout%fStencil%QQN if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then ! Equilibrium on precollision density and velocity ! calculate equilibrium density feq = layout%weight( idir ) & & * dens & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) ! calculate equilibrium density feq_neigh = layout%weight( idir ) & & * dens_neigh & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_neigh(:)*vel_neigh(:)) * 0.5_rk * cs2inv ) ) ! Equilibrium on wall model velocity ! calculate equilibrium density feq_bnd = layout%weight( idir ) & & * dens_bnd & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:)*velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd + 0.5_rk*( f_pre(iDir) - fEq ) & & + 0.5_rk*( f_neigh(iDir) - fEq_neigh ) end if end do end do end subroutine turbulent_wall ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall. !! It uses wall model to compute velocity on the boundary node. !! The implementation is based on the following paper: !! Haussmann, Marc; Ries, Florian; Jeppener-Haltenhoff, Jonathan B.; Li, !! Yongxiang; Schmidt, Marius; Welch, Cooper et al. (2020): Evaluation of a !! Near-Wall-Modeled Large Eddy Lattice Boltzmann Method for the Analysis of !! Complex Flows Relevant to IC Engines. In Computation 8 (2), p. 43. !! DOI: 10.3390/computation8020043. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall', !! curved = true, !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_libb( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! real(kind=rk) :: fIn, fOut, fNgh real(kind=rk) :: cIn, cOut, cNgh integer :: iLink integer :: iElem, iDir, QQ, elemPos, invDir real(kind=rk) :: f_pre(layout%fStencil%QQ) real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: f_preBuffer(layout%fStencil%QQ * globBC%nElems(iLevel)) real(kind=rk) :: dens, dens_inv, vel(3), velSW(3) real(kind=rk) :: dens_neigh, dens_neigh_inv, vel_neigh(3), qVal real(kind=rk) :: fEq, fEq_bnd, fEq_neigh, fNEq real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! 1. First perform libb step i.e. curved boundary step !NEC$ ivdep !DIR$ ivdep !IBM* independent do iLink = 1, me%links(iLevel)%nVals cIn = me%bouzidi(iLevel)% cIn( iLink ) cOut = me%bouzidi(iLevel)%cOut( iLink ) cNgh = me%bouzidi(iLevel)%cNgh( iLink ) fIn = bcBuffer( me%bouzidi(iLevel)% inPos(iLink) ) fOut = bcBuffer( me%bouzidi(iLevel)%outPos(iLink) ) fNgh = me%neigh(iLevel)%computeNeighBuf(me%bouzidi(iLevel)%nghPos(iLink)) state( me%links(iLevel)%val(iLink) ) = cIn*fIn + cOut*fOut + cNgh*fNgh end do ! iLink ! 2. Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! 3. Velocity correction step ! Get pre-collision of boundary elementy in a buffer since it gets ! overwritten in the next element loop do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) do iDir = 1, layout%fStencil%QQ f_preBuffer((iElem-1)*QQ + iDir) = state( & & neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0 ) end do end do do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state after ! curved boundary step f_pre = f_preBuffer( (iElem-1)*QQ+1: iElem*QQ ) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens = sum(f_pre) dens_inv = 1.0_rk/dens dens_neigh = sum(f_neigh) dens_neigh_inv = 1.0_rk/dens_neigh ! velocity vel(1) = sum( f_pre * layout%fstencil%cxdirrk(1, :) ) * dens_inv vel(2) = sum( f_pre * layout%fstencil%cxdirrk(2, :) ) * dens_inv vel(3) = sum( f_pre * layout%fstencil%cxdirrk(3, :) ) * dens_inv vel_neigh(1) = sum( f_neigh * layout%fstencil%cxdirrk(1, :) ) * dens_neigh_inv vel_neigh(2) = sum( f_neigh * layout%fstencil%cxdirrk(2, :) ) * dens_neigh_inv vel_neigh(3) = sum( f_neigh * layout%fstencil%cxdirrk(3, :) ) * dens_neigh_inv ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) do iDir = 1, layout%fStencil%QQN if ( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then invDir = layout%fStencil%cxDirInv(iDir) ! qValues qVal = globBC%elemLvl( iLevel )%qVal%val( invDir, iElem ) ! Equilibrium on precollision density and velocity ! calculate equilibrium density feq = layout%weight( idir ) & & * dens & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel(:)*vel(:)) * 0.5_rk * cs2inv ) ) ! Equilibrium on wall model velocity ! calculate equilibrium density feq_bnd = layout%weight( idir ) & & * dens & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:)*velsw(:)) * 0.5_rk * cs2inv ) ) if (qVal >= 0.75_rk) then fNeq = f_pre(iDir) - fEq else ! calculate equilibrium density feq_neigh = layout%weight( idir ) & & * dens_neigh & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_neigh(:)*vel_neigh(:)) * 0.5_rk * cs2inv ) ) fNeq = qVal * (f_pre(iDir) - fEq) & & + (1.0_rk-qVal) * (f_neigh(iDir) - fEq_neigh) end if state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd + fNeq end if end do end do end subroutine turbulent_wall_libb ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on Guo's nonequilibrium extrapolation. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium plus non-equilibrium. !! Density is computed using Zho-He approach for straight walls. !! "On pressure and velocity boundary conditions for the lattice Boltzmann !! BGK model", Physics of Fluids 9, 1591-1598 (1997) !! https://doi.org/10.1063/1.869307 !! !! non-equilibrium are computed from PDF on neighbor and extrapolated to !! boundary. This routine is used for straight wall boundaries. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_noneq_expol', !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_noneq_expol( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos, normalInd_inv real(kind=rk) :: f_pre(layout%fStencil%QQ) real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: dens_neigh, dens_inv, velSW(3), vel_neigh(3) real(kind=rk) :: dens_bnd, vel_normal, normal(3) real(kind=rk) :: fEq, fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens_neigh = sum(f_neigh) dens_inv = 1.0_rk/dens_neigh ! velocity vel_neigh(1) = sum( f_neigh * layout%fstencil%cxdirrk(1, :) ) * dens_inv vel_neigh(2) = sum( f_neigh * layout%fstencil%cxdirrk(2, :) ) * dens_inv vel_neigh(3) = sum( f_neigh * layout%fstencil%cxdirrk(3, :) ) * dens_inv ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) ! calculate density at boundary for straight wall using Zho-He approach normalInd_inv = layout%fStencil%cxDirInv( globBC%elemLvl(iLevel) & & %normalInd%val(iElem) ) normal = layout%fStencil%cxDirRK(:, normalInd_inv) vel_normal = dot_product(velSW, normal) dens_bnd = 0.0_rk f_pre = me%neigh(iLevel)%computeNeighBuf( (iElem-1)*QQ+1: iElem*QQ ) do iDir = 1, layout%fStencil%QQ if ( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then dens_bnd = dens_bnd + f_pre(layout%fStencil%cxDirInv(iDir)) else dens_bnd = dens_bnd + f_pre(iDir) end if end do dens_bnd = dens_bnd / (1.0_rk + vel_normal) do iDir = 1, layout%fStencil%QQ ! Equilibrium on precollision density and velocity ! calculate equilibrium density feq = layout%weight( idir ) & & * dens_neigh & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_neigh(:)*vel_neigh(:)) * 0.5_rk * cs2inv ) ) ! Equilibrium on wall model velocity ! calculate equilibrium density feq_bnd = layout%weight( idir ) & & * dens_bnd & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:)*velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd + ( f_neigh(iDir) - fEq ) end do end do end subroutine turbulent_wall_noneq_expol ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on Guo's nonequilibrium extrapolation. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium plus non-equilibrium. Density and non-equilibrium are commputed !! from PDF on neighbor and extrapolated to boundary. !! This routine is used for curved boundaries. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_noneq_expol', !! curved = true, !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_noneq_expol_curved( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: dens_neigh, dens_inv, velSW(3), vel_neigh(3) real(kind=rk) :: fEq, fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens_neigh = sum(f_neigh) dens_inv = 1.0_rk/dens_neigh ! velocity vel_neigh(1) = sum( f_neigh * layout%fstencil%cxdirrk(1, :) ) * dens_inv vel_neigh(2) = sum( f_neigh * layout%fstencil%cxdirrk(2, :) ) * dens_inv vel_neigh(3) = sum( f_neigh * layout%fstencil%cxdirrk(3, :) ) * dens_inv ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) do iDir = 1, layout%fStencil%QQ ! Equilibrium on precollision density and velocity ! calculate equilibrium density feq = layout%weight( idir ) & & * dens_neigh & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * vel_neigh(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_neigh(:)*vel_neigh(:)) * 0.5_rk * cs2inv ) ) ! Equilibrium on wall model velocity ! calculate equilibrium density feq_bnd = layout%weight( idir ) & & * dens_neigh & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:)*velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd + ( f_neigh(iDir) - fEq ) end do end do end subroutine turbulent_wall_noneq_expol_curved ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on equilibrium BC. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium. !! Density is computed using Zho-He approach for straight walls. !! "On pressure and velocity boundary conditions for the lattice Boltzmann !! BGK model", Physics of Fluids 9, 1591-1598 (1997) !! https://doi.org/10.1063/1.869307 !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_eq', !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_eq( me, state, bcBuffer, globBC, levelDesc, tree, & & nSize, iLevel, sim_time, neigh, layout, fieldProp, varPos, & & nScalars, varSys, derVarPos, physics, iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos, normalInd_inv real(kind=rk) :: f_pre(layout%fStencil%QQ) real(kind=rk) :: velSW(3) real(kind=rk) :: dens_bnd, vel_normal, normal(3) real(kind=rk) :: fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) ! calculate density at boundary for straight wall using Zho-He approach normalInd_inv = layout%fStencil%cxDirInv( globBC%elemLvl(iLevel) & & %normalInd%val(iElem) ) normal = layout%fStencil%cxDirRK(:, normalInd_inv) vel_normal = dot_product(velSW, normal) dens_bnd = 0.0_rk f_pre = me%neigh(iLevel)%computeNeighBuf( (iElem-1)*QQ+1: iElem*QQ ) do iDir = 1, layout%fStencil%QQ if ( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then dens_bnd = dens_bnd + f_pre(layout%fStencil%cxDirInv(iDir)) else dens_bnd = dens_bnd + f_pre(iDir) end if end do dens_bnd = dens_bnd / (1.0_rk + vel_normal) do iDir = 1, layout%fStencil%QQ ! Equilibrium on wall model velocity ! calculate equilibrium density feq_bnd = layout%weight( idir ) & & * dens_bnd & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:)*velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd end do end do end subroutine turbulent_wall_eq ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on equilibrium BC. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium plus non-equilibrium. Density is commputed !! from PDF on neighbor and extrapolated to boundary. !! This routine is used for curved boundaries. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_eq', !! curved = true, !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_eq_curved( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: dens_neigh, velSW(3) real(kind=rk) :: fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens_neigh = sum(f_neigh) ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) do iDir = 1, layout%fStencil%QQ ! Equilibrium on wall model velocity ! calculate equilibrium density feq_bnd = layout%weight( idir ) & & * dens_neigh & & * ( 1._rk & & + ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & + (sum(layout%fstencil%cxdirrk(:,idir) * velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir) * velsw(:))) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:)*velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd end do end do end subroutine turbulent_wall_eq_curved ! ************************************************************************** ! ! ************************************************************************** ! !> Calculation stream-wise velocity compononent from wall function and !! friction velocity, stream-wise unit vector and turbulent viscosity with !! mixing length formulation. subroutine calcVelSW_unitSW_velTau_tVisc(velSW, unitSW, turbwallFunc, & & nElems, elemPos, normalInd, neighBufferPre, viscKine, & & turbulence, stencil, iLevel) ! -------------------------------------------------------------------------- !> Stream-wise velocity component from wall function real(kind=rk), intent(out) :: velSW(:) !> Stream-wise unit vector real(kind=rk), intent(out) :: unitSW(:,:) !> Turbulent wall model type contains viscosity, velTau, distToBnd and !! function pointers to compute velTau and velSW type(mus_turb_wallFunc_type), intent(inout) :: turbwallFunc !> Number of elements in current boundary integer, intent(in) :: nElems !> Current element position in state array. (Used to access viscosity) integer, intent(in) :: elemPos(:) !> Normal index for evey boundary element integer, intent(in) :: normalInd(:) !> Pre-collision from 1st and 2nd fluid neighbor real(kind=rk), intent(in) :: neighBufferPre(:,:) !> Kinematic viscosity type(mus_viscosity_type) :: viscKine !> turbulence model type type(mus_turbulence_type), intent(in) :: turbulence !> Fluid stencil type(tem_stencilHeader_type), intent(in) :: stencil !> Current level integer, intent(in) :: iLevel ! -------------------------------------------------------------------------- integer :: iElem, QQ, normalInd_inv real(kind=rk) :: f_neigh(stencil%QQ) real(kind=rk) :: vel_neigh(3), dens_neigh_inv, normal(3), vec(3), vec_mag real(kind=rk) :: velSW_neigh(nElems) real(kind=rk) :: yPlus, gradU, velSW_neigh_2 ! -------------------------------------------------------------------------- QQ = stencil%QQ ! Calculate stream-wise velocity component on first neighbor do iElem = 1, nElems ! Get density and velocity on second fluid element along the normal ! direction i.e. 1st neighbor. f_neigh = neighBufferPre(1, (iElem-1)*QQ+1:(iElem-1)*QQ+QQ) dens_neigh_inv = 1.0_rk / sum(f_neigh) ! Get velocity from pdf, inverse of density and stencil as input. vel_neigh(1) = sum( f_neigh * stencil%cxdirrk(1, :) ) * dens_neigh_inv vel_neigh(2) = sum( f_neigh * stencil%cxdirrk(2, :) ) * dens_neigh_inv vel_neigh(3) = sum( f_neigh * stencil%cxdirrk(3, :) ) * dens_neigh_inv ! Calculate local stream-wise unit vector ! e' = (u - (u . n) . n) / |(u - (u . n) . n| ! ! NormalInd in bc elem is pointing into the domain so we invert it normalInd_inv = stencil%cxDirInv( normalInd(iElem) ) normal = stencil%cxDirRK(:, normalInd_inv) vec = vel_neigh - dot_product( vel_neigh, normal ) * normal vec_mag = sqrt(dot_product(vec, vec)) ! Unit vector unitSW(:, iElem) = vec / vec_mag ! stream-wise velocity component velSW_neigh(iElem) = dot_product(vel_neigh, unitSW(:, iElem)) end do ! calculate friction velocity on neighbor element call turbwallFunc%calcFricVel( & & velTau = turbwallFunc%dataOnLvl(iLevel)%velTau, & & velSW = velSW_neigh, & & distToBnd = turbwallFunc%dataOnLvl(iLevel)%neighDistToBnd, & & viscKine = viscKine%dataOnLvl(iLevel)%val( elemPos(:) ), & & nElems = nElems ) ! calculate stream-wise velocity on boundary element call turbwallFunc%calcStreamWiseVel( & & velSW = velSW, & & velTau = turbwallFunc%dataOnLvl(iLevel)%velTau, & & distToBnd = turbwallFunc%dataOnLvl(iLevel)%distToBnd, & & viscKine = viscKine%dataOnLvl(iLevel)%val( elemPos(:) ), & & nElems = nElems ) ! Calculate Turbulent viscosity according to mixing length formulation ! with von karman constant. ! only if LES turbulence model is smagorinsky because turb viscosity ! from Vreman and WALE reduces towards the wall ! nu_t = (k*y)**2 * |dudy| if (turbulence%active .and. & & trim(turbulence%config%model) == 'smagorinsky') then do iElem = 1, nElems ! Get density and velocity on third fluid element along the normal ! direction i.e. 2nd neighbor. f_neigh = neighBufferPre(2, (iElem-1)*QQ+1:(iElem-1)*QQ+QQ) dens_neigh_inv = 1.0_rk / sum(f_neigh) ! Get velocity from pdf, inverse of density and stencil as input. vel_neigh(1) = sum( f_neigh * stencil%cxdirrk(1, :) ) * dens_neigh_inv vel_neigh(2) = sum( f_neigh * stencil%cxdirrk(2, :) ) * dens_neigh_inv vel_neigh(3) = sum( f_neigh * stencil%cxdirrk(3, :) ) * dens_neigh_inv velSW_neigh_2 = dot_product(vel_neigh, unitSW(:, iElem)) ! computed velocity gradient from second order forward difference gradU = abs( - 3.0_rk * velSW(iElem) + 4.0_rk * velSW_neigh(iElem) & & - velSW_neigh_2 ) / 2.0_rk turbwallFunc%dataOnLvl(iLevel)%tVisc(iElem) = ( turbWallFunc%vonKarman & & * turbwallFunc%dataOnLvl(iLevel)%distToBnd(iElem) )**2.0_rk & & * gradU end do ! Apply van Driest damping ! nu_t = nu_t * (1-exp(-yplus/Aplus))**2 if (turbWallFunc%useVanDriest) then do iElem = 1, nElems yPlus = turbwallFunc%dataOnLvl(iLevel)%distToBnd(iElem) & & * turbwallFunc%dataOnLvl(iLevel)%velTau(iElem) & & / viscKine%dataOnLvl(iLevel)%val(elemPos(iElem)) turbwallFunc%dataOnLvl(iLevel)%tVisc(iElem) & & = turbwallFunc%dataOnLvl(iLevel)%tVisc(iElem) & & * (1.0_rk-exp(-yPlus/vd_Aplus))**2.0_rk end do end if end if end subroutine calcVelSW_unitSW_velTau_tVisc ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on Guo's nonequilibrium extrapolation. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium plus non-equilibrium. !! Density is computed using Zho-He approach for straight walls. !! "On pressure and velocity boundary conditions for the lattice Boltzmann !! BGK model", Physics of Fluids 9, 1591-1598 (1997) !! https://doi.org/10.1063/1.869307 !! !! non-equilibrium are computed from PDF on neighbor and extrapolated to !! boundary. This routine is used for straight wall boundaries. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_noneq_expol', !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_noneq_expol_incomp( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos, normalInd_inv real(kind=rk) :: f_pre(layout%fStencil%QQ) real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: dens_neigh, velSW(3), vel_neigh(3) real(kind=rk) :: dens_bnd, vel_normal, normal(3) real(kind=rk) :: fEq, fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc_incomp( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens_neigh = sum(f_neigh) ! velocity vel_neigh(1) = sum( f_neigh * layout%fstencil%cxdirrk(1, :) ) * rho0inv vel_neigh(2) = sum( f_neigh * layout%fstencil%cxdirrk(2, :) ) * rho0inv vel_neigh(3) = sum( f_neigh * layout%fstencil%cxdirrk(3, :) ) * rho0inv ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) ! calculate density at boundary for straight wall using Zho-He approach normalInd_inv = layout%fStencil%cxDirInv( globBC%elemLvl(iLevel) & & %normalInd%val(iElem) ) normal = layout%fStencil%cxDirRK(:, normalInd_inv) vel_normal = dot_product(velSW, normal) dens_bnd = 0.0_rk f_pre = me%neigh(iLevel)%computeNeighBuf( (iElem-1)*QQ+1: iElem*QQ ) do iDir = 1, layout%fStencil%QQ if ( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then dens_bnd = dens_bnd + f_pre(layout%fStencil%cxDirInv(iDir)) else dens_bnd = dens_bnd + f_pre(iDir) end if end do dens_bnd = dens_bnd / (1.0_rk + vel_normal) do iDir = 1, layout%fStencil%QQ ! Equilibrium on precollision density and velocity feq = layout%weight( idir ) & & * ( dens_neigh + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_neigh(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_neigh(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_neigh(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_neigh(:) * vel_neigh(:)) * 0.5_rk * cs2inv ) ) ! Equilibrium on wall model velocity feq_bnd = layout%weight( idir ) & & * ( dens_bnd + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:) * velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd + ( f_neigh(iDir) - fEq ) end do end do end subroutine turbulent_wall_noneq_expol_incomp ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on Guo's nonequilibrium extrapolation. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium plus non-equilibrium. Density and non-equilibrium are commputed !! from PDF on neighbor and extrapolated to boundary. !! This routine is used for curved boundaries. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_noneq_expol', !! curved = true, !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_noneq_expol_curved_incomp( me, state, bcBuffer, & & globBC, levelDesc, tree, nSize, iLevel, sim_time, neigh, & & layout, fieldProp, varPos, nScalars, varSys, derVarPos, & & physics, iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: dens_neigh, dens_inv, velSW(3), vel_neigh(3) real(kind=rk) :: fEq, fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc_incomp( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens_neigh = sum(f_neigh) dens_inv = 1.0_rk/dens_neigh ! velocity vel_neigh(1) = sum( f_neigh * layout%fstencil%cxdirrk(1, :) ) * dens_inv vel_neigh(2) = sum( f_neigh * layout%fstencil%cxdirrk(2, :) ) * dens_inv vel_neigh(3) = sum( f_neigh * layout%fstencil%cxdirrk(3, :) ) * dens_inv ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) do iDir = 1, layout%fStencil%QQ ! Equilibrium on precollision density and velocity feq = layout%weight( idir ) & & * ( dens_neigh + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_neigh(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*vel_neigh(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*vel_neigh(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(vel_neigh(:) * vel_neigh(:)) * 0.5_rk * cs2inv ) ) ! Equilibrium on wall model velocity feq_bnd = layout%weight( idir ) & & * ( dens_neigh + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:) * velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd + ( f_neigh(iDir) - fEq ) end do end do end subroutine turbulent_wall_noneq_expol_curved_incomp ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on equilibrium BC. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium. !! Density is computed using Zho-He approach for straight walls. !! "On pressure and velocity boundary conditions for the lattice Boltzmann !! BGK model", Physics of Fluids 9, 1591-1598 (1997) !! https://doi.org/10.1063/1.869307 !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_eq', !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_eq_incomp( me, state, bcBuffer, globBC, levelDesc, & & tree, nSize, iLevel, sim_time, neigh, layout, fieldProp, varPos,& & nScalars, varSys, derVarPos, physics, iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos, normalInd_inv real(kind=rk) :: f_pre(layout%fStencil%QQ) real(kind=rk) :: velSW(3) real(kind=rk) :: dens_bnd, vel_normal, normal(3) real(kind=rk) :: fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc_incomp( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) ! calculate density at boundary for straight wall using Zho-He approach normalInd_inv = layout%fStencil%cxDirInv( globBC%elemLvl(iLevel) & & %normalInd%val(iElem) ) normal = layout%fStencil%cxDirRK(:, normalInd_inv) vel_normal = dot_product(velSW, normal) dens_bnd = 0.0_rk f_pre = me%neigh(iLevel)%computeNeighBuf( (iElem-1)*QQ+1: iElem*QQ ) do iDir = 1, layout%fStencil%QQ if ( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then dens_bnd = dens_bnd + f_pre(layout%fStencil%cxDirInv(iDir)) else dens_bnd = dens_bnd + f_pre(iDir) end if end do dens_bnd = dens_bnd / (1.0_rk + vel_normal) do iDir = 1, layout%fStencil%QQ ! Equilibrium on wall model velocity feq_bnd = layout%weight( idir ) & & * ( dens_bnd + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:) * velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd end do end do end subroutine turbulent_wall_eq_incomp ! ************************************************************************** ! ! ************************************************************************** ! !> BC routine for turbulent wall based on equilibrium BC. !! The implementation is based on the following paper: !! 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. !! !! It uses wall model to compute velocity on the boundary node. !! All directions of PDF in the boundary elements are updated with !! Equilibrium plus non-equilibrium. Density is commputed !! from PDF on neighbor and extrapolated to boundary. !! This routine is used for curved boundaries. !! !! Usage !! ----- !!```lua !!boundary_condition = { !! { !! label = 'wall', !! kind = 'turbulent_wall_eq', !! curved = true, !! wall_model = 'musker', !! nonlinear_solver = 'fixed_point' !! } !!} !! !! This subroutine's interface must match the abstract interface definition !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be !! callable via [[boundary_type:fnct]] function pointer. subroutine turbulent_wall_eq_curved_incomp( me, state, bcBuffer, globBC, & & levelDesc, tree, nSize, iLevel, sim_time, neigh, layout, & & fieldProp, varPos, nScalars, varSys, derVarPos, physics, & & iField, mixture ) ! ------------------------------------------------------------------------ ! !> global boundary type class(boundary_type) :: me !> Current state vector of iLevel real(kind=rk), intent(inout) :: state(:) !> size of state array ( in terms of elements ) integer, intent(in) :: nSize !> state values of boundary elements of all fields of iLevel real(kind=rk), intent(in) :: bcBuffer(:) !> iLevel descriptor type(tem_levelDesc_type), intent(in) :: levelDesc !> Treelm Mesh type(treelmesh_type), intent(in) :: tree !> fluid parameters and properties type(mus_field_prop_type), intent(in) :: fieldProp !> stencil layout information type(mus_scheme_layout_type), intent(in) :: layout !> the level On which this boundary was invoked integer, intent(in) :: iLevel !> connectivity array corresponding to state vector integer, intent(in) :: neigh(:) !> global time information type(tem_time_type), intent(in) :: sim_time !> pointer to field variable in the state vector integer, intent(in) :: varPos(:) !> number of Scalars in the scheme var system integer, intent(in) :: nScalars !> scheme variable system type(tem_varSys_type), intent(in) :: varSys !> position of derived quantities in varsys type(mus_derVarPos_type), intent(in) :: derVarPos !> scheme global boundary type type(glob_boundary_type), intent(in) :: globBC !> scheme global boundary type type(mus_physics_type), intent(in) :: physics !> current field integer, intent(in) :: iField !> mixture info type(mus_mixture_type), intent(in) :: mixture ! ------------------------------------------------------------------------ ! integer :: iElem, iDir, QQ, elemPos real(kind=rk) :: f_neigh(layout%fStencil%QQ) real(kind=rk) :: dens_neigh, velSW(3) real(kind=rk) :: fEq_bnd real(kind=rk) :: unitSW(3, globBC%nElems(iLevel)) real(kind=rk) :: velSW_bnd(globBC%nElems(iLevel)) ! ------------------------------------------------------------------------ ! ! Load stencil for velocity space (DdQq) with q=QQ QQ = layout%fStencil%QQ ! Calculate stream-wise velocity component from wall function call calcVelSW_unitSW_velTau_tVisc_incomp( & & velSW = velSW_bnd, & & unitSW = unitSW, & & turbwallFunc = me%turbwallFunc, & & nElems = globBC%nElems(iLevel), & & elemPos = globBC%elemLvl(iLevel)%elem%val(:), & & normalInd = globBC%elemLvl(iLevel)%normalInd%val(:), & & neighBufferPre = me%neigh(iLevel)%neighBufferPre_nNext, & & viscKine = fieldProp%fluid%viscKine, & & turbulence = fieldProp%fluid%turbulence, & & stencil = layout%fStencil, & & iLevel = iLevel ) ! Velocity correction on boundary element incoming direction do iElem = 1, globBC%nElems(iLevel) ! Current element position in state array. (Used for state in step 10.) elemPos = globBC%elemLvl(iLevel)%elem%val(iElem) ! Compute density and velocity from pre-collision state on neighbor f_neigh = me%neigh(iLevel)%neighBufferPre_nNext( 1, & & (iElem-1)*QQ+1:(iElem-1)*QQ+QQ ) dens_neigh = sum(f_neigh) ! stream-wise velocity on boundary element from wall model velSW = velSW_bnd(iElem) * unitSW(:, iElem) do iDir = 1, layout%fStencil%QQ ! Equilibrium on wall model velocity feq_bnd = layout%weight( idir ) & & * ( dens_neigh + rho0 & & * ( cs2inv & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & + ( sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) & & * sum(layout%fstencil%cxdirrk(:,idir)*velsw(:)) ) & & * cs2inv * cs2inv * 0.5_rk & & - sum(velsw(:) * velsw(:)) * 0.5_rk * cs2inv ) ) state( neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) & & = fEq_bnd end do end do end subroutine turbulent_wall_eq_curved_incomp ! ************************************************************************** ! ! ************************************************************************** ! !> Calculation stream-wise velocity compononent from wall function and !! friction velocity, stream-wise unit vector and turbulent viscosity with !! mixing length formulation. subroutine calcVelSW_unitSW_velTau_tVisc_incomp(velSW, unitSW, turbwallFunc, & & nElems, elemPos, normalInd, neighBufferPre, viscKine, & & turbulence, stencil, iLevel) ! -------------------------------------------------------------------------- !> Stream-wise velocity component from wall function real(kind=rk), intent(out) :: velSW(:) !> Stream-wise unit vector real(kind=rk), intent(out) :: unitSW(:,:) !> Turbulent wall model type contains viscosity, velTau, distToBnd and !! function pointers to compute velTau and velSW type(mus_turb_wallFunc_type), intent(inout) :: turbwallFunc !> Number of elements in current boundary integer, intent(in) :: nElems !> Current element position in state array. (Used to access viscosity) integer, intent(in) :: elemPos(:) !> Normal index for evey boundary element integer, intent(in) :: normalInd(:) !> Pre-collision from 1st and 2nd fluid neighbor real(kind=rk), intent(in) :: neighBufferPre(:,:) !> Kinematic viscosity type(mus_viscosity_type) :: viscKine !> turbulence model type type(mus_turbulence_type), intent(in) :: turbulence !> Fluid stencil type(tem_stencilHeader_type), intent(in) :: stencil !> Current level integer, intent(in) :: iLevel ! -------------------------------------------------------------------------- integer :: iElem, QQ, normalInd_inv real(kind=rk) :: f_neigh(stencil%QQ) real(kind=rk) :: vel_neigh(3), normal(3), vec(3), vec_mag real(kind=rk) :: velSW_neigh(nElems) real(kind=rk) :: yPlus, gradU, velSW_neigh_2 ! -------------------------------------------------------------------------- QQ = stencil%QQ ! Calculate stream-wise velocity component on first neighbor do iElem = 1, nElems ! Get density and velocity on second fluid element along the normal ! direction i.e. 1st neighbor. f_neigh = neighBufferPre(1, (iElem-1)*QQ+1:(iElem-1)*QQ+QQ) ! Get velocity from pdf, inverse of density and stencil as input. vel_neigh(1) = sum( f_neigh * stencil%cxdirrk(1, :) ) * rho0inv vel_neigh(2) = sum( f_neigh * stencil%cxdirrk(2, :) ) * rho0inv vel_neigh(3) = sum( f_neigh * stencil%cxdirrk(3, :) ) * rho0inv ! Calculate local stream-wise unit vector ! e' = (u - (u . n) . n) / |(u - (u . n) . n| ! ! NormalInd in bc elem is pointing into the domain so we invert it normalInd_inv = stencil%cxDirInv( normalInd(iElem) ) normal = stencil%cxDirRK(:, normalInd_inv) vec = vel_neigh - dot_product( vel_neigh, normal ) * normal vec_mag = sqrt(dot_product(vec, vec)) ! Unit vector unitSW(:, iElem) = vec / vec_mag ! stream-wise velocity component velSW_neigh(iElem) = dot_product(vel_neigh, unitSW(:, iElem)) end do ! calculate friction velocity on neighbor element call turbwallFunc%calcFricVel( & & velTau = turbwallFunc%dataOnLvl(iLevel)%velTau, & & velSW = velSW_neigh, & & distToBnd = turbwallFunc%dataOnLvl(iLevel)%neighDistToBnd, & & viscKine = viscKine%dataOnLvl(iLevel)%val( elemPos(:) ), & & nElems = nElems ) ! calculate stream-wise velocity on boundary element call turbwallFunc%calcStreamWiseVel( & & velSW = velSW, & & velTau = turbwallFunc%dataOnLvl(iLevel)%velTau, & & distToBnd = turbwallFunc%dataOnLvl(iLevel)%distToBnd, & & viscKine = viscKine%dataOnLvl(iLevel)%val( elemPos(:) ), & & nElems = nElems ) ! Calculate Turbulent viscosity according to mixing length formulation ! with von karman constant. ! only if LES turbulence model is smagorinsky because turb viscosity ! from Vreman and WALE reduces towards the wall ! nu_t = (k*y)**2 * |dudy| if (turbulence%active .and. & & trim(turbulence%config%model) == 'smagorinsky') then do iElem = 1, nElems ! Get density and velocity on third fluid element along the normal ! direction i.e. 2nd neighbor. f_neigh = neighBufferPre(2, (iElem-1)*QQ+1:(iElem-1)*QQ+QQ) ! Get velocity from pdf, inverse of density and stencil as input. vel_neigh(1) = sum( f_neigh * stencil%cxdirrk(1, :) ) * rho0inv vel_neigh(2) = sum( f_neigh * stencil%cxdirrk(2, :) ) * rho0inv vel_neigh(3) = sum( f_neigh * stencil%cxdirrk(3, :) ) * rho0inv velSW_neigh_2 = dot_product(vel_neigh, unitSW(:, iElem)) ! computed velocity gradient from second order forward difference gradU = abs( - 3.0_rk * velSW(iElem) + 4.0_rk * velSW_neigh(iElem) & & - velSW_neigh_2 ) / 2.0_rk turbwallFunc%dataOnLvl(iLevel)%tVisc(iElem) = ( turbWallFunc%vonKarman & & * turbwallFunc%dataOnLvl(iLevel)%distToBnd(iElem) )**2.0_rk & & * gradU end do ! Apply van Driest damping ! nu_t = nu_t * (1-exp(-yplus/Aplus))**2 if (turbWallFunc%useVanDriest) then do iElem = 1, nElems yPlus = turbwallFunc%dataOnLvl(iLevel)%distToBnd(iElem) & & * turbwallFunc%dataOnLvl(iLevel)%velTau(iElem) & & / viscKine%dataOnLvl(iLevel)%val(elemPos(iElem)) turbwallFunc%dataOnLvl(iLevel)%tVisc(iElem) & & = turbwallFunc%dataOnLvl(iLevel)%tVisc(iElem) & & * (1.0_rk-exp(-yPlus/vd_Aplus))**2.0_rk end do end if end if end subroutine calcVelSW_unitSW_velTau_tVisc_incomp ! ************************************************************************** ! end module mus_bc_fluid_turbulent_module ! *****************************************************************************!