mus_wall_function_musker_module.f90 Source File

! Copyright (c) 2023 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.
! **************************************************************************** !
!> This module contains data types, function and routines for wall function
!! computations relative to Musker profile.
!! 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.
!!
!! author: Gregorio Gerardo Spinelli
module mus_wall_function_musker_module

  ! include treelm modules
  use env_module,                        only: rk
  use mus_wall_function_abstract_module, only: mus_wall_function_type

  implicit none

  private

  public :: mus_wall_function_musker_type

  !> extend the abstract subclass mus_wall_function_type
  type, extends(mus_wall_function_type) :: mus_wall_function_musker_type

  contains

    !> function to get uPlus
    procedure, nopass :: get_uPlus

    !> function to apply the newon method
    ! F = uPlus_n - uPlus_(n+1)
    ! here we compute the derivative of uPlus_(n+1) w.r.t. u_tau
    ! this function computes the derivative of uPlus with respect to uTau
    procedure, nopass :: get_d_uPlus_d_uTau

  end type mus_wall_function_musker_type


contains


  !> function to get uPlus
  pure function get_uPlus( yPlus ) result( uPlus )

    !> yPlus
    real(kind=rk), intent(in) :: yPlus
    !> output is uPlus
    real(kind=rk) :: uPlus
    !uPlus = ( 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   )

    ! above expression simplified with sympy in python
    uPlus = 0.434294481903252_rk * log( ( 0.0943396226415094_rk * yPlus + 1._rk )**9.6_rk /  &
      &     ( 0.0116279069767442_rk * yPlus**2 - 0.0947674418604651_rk * yPlus + 1_rk )**2 ) &
      &     + 5.424_rk * atan( 0.119760479041916_rk * yPlus - 0.488023952095808_rk )         &
      &     + 2.46666038465466_rk

  end function get_uPlus

  !> function to get the derivative of uPlus with respect to uTau
  pure function get_d_uPlus_d_uTau( y, uTau, nu ) result( d_uPlus_d_uTau )

    !> vertical distance from the wall
    real(kind=rk), intent(in) :: y
    !> uTau at iteration n
    real(kind=rk), intent(in) :: uTau
    !> dynamic viscosity
    real(kind=rk), intent(in) :: nu
    !> output is derivative of uPlus with respect to uTau
    real(kind=rk) :: d_uPlus_d_uTau
    ! ------------------------------------------------------------------------------
    real(kind=rk) :: uTau_y, yPlus, nu_sqr, uTau_y_sqr, nu_uTau_y, A, B, C, D, E
    uTau_y = uTau * y
    nu_uTau_y = uTau_y * nu
    yPlus = uTau_y / nu
    nu_sqr = nu**2
    uTau_y_sqr = uTau_y**2
    A = ( 1._rk + 0.0943396226415094_rk * yPlus )
    B = nu_sqr - 0.0947674418604651_rk * nu_uTau_y + 0.0116279069767442_rk * uTau_y_sqr
    C = ( nu - 0.245398773006135_rk * uTau_y )**2
    D = nu_sqr + 0.238167377819212_rk * C
    E = A**9.6_rk

    ! obtained as diff from sympy in python
    d_uPlus_d_uTau = 0.434294481903252_rk * y * ( 1.49571515501793_rk * nu_sqr        &
      & * E * B + D * ( nu * E * ( 0.18953488372093_rk * nu                           &
      & - 0.0465116279069767_rk * uTau_y ) + 0.905660377358491_rk * A**8.6_rk * B ) ) &
      & / ( nu * E * D * B )

  end function get_d_uPlus_d_uTau

end module mus_wall_function_musker_module