flow_profile.f90

Go to the documentation of this file.

! Copyright (c) 2021, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * 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.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "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 THE
! COPYRIGHT OWNER 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.
!
module flow_profile
  use num_types, only : rp
  implicit none
  private
 
  abstract interface
     function blasius_profile(y, delta, u)
       import rp
       real(kind=rp), intent(in) :: y, delta, u
       real(kind=rp) :: blasius_profile
     end function blasius_profile
  end interface
 
  public :: blasius_profile, blasius_linear, blasius_cubic, blasius_quadratic, &
            blasius_quartic, blasius_sin, blasius_tanh
 
contains
 
  function blasius_linear(y, delta, u)
    real(kind=rp), intent(in) :: y, delta, u
    real(kind=rp) :: blasius_linear
    real(kind=rp) :: arg
 
    arg = y /delta
    if (arg .gt. 1.0_rp) then
       blasius_linear = u
    else
       blasius_linear = u * (y / delta)
    end if
 
  end function blasius_linear
 
  function blasius_quadratic(y, delta, u)
    real(kind=rp), intent(in) :: y, delta, u
    real(kind=rp) :: blasius_quadratic
    real(kind=rp) :: arg
 
    arg = ( 2.0_rp * (y / delta) - (y / delta)**2 )
 
    if (arg .gt. 1.0_rp) then
       blasius_quadratic = u
    else
       blasius_quadratic = u * arg
    end if
 
  end function blasius_quadratic
 
  function blasius_cubic(y, delta, u)
    real(kind=rp), intent(in) :: y, delta, u
    real(kind=rp) :: blasius_cubic
    real(kind=rp) :: arg
 
    arg = ( 3.0_rp / 2.0_rp * (y / delta) - 0.5_rp * (y / delta)**3 )
 
    if (arg .gt. 1.0_rp) then
       blasius_cubic = 1.0_rp
    else
       blasius_cubic = u * arg
    end if
 
  end function blasius_cubic
 
  function blasius_quartic(y, delta, u)
    real(kind=rp), intent(in) :: y, delta, u
    real(kind=rp) :: blasius_quartic
    real(kind=rp) :: arg
 
    arg = 2.0_rp * (y / delta) - 2.0_rp * (y / delta)**3 + (y / delta)**4
 
    if (arg .gt. 1.0_rp) then
       blasius_quartic = u
    else
       blasius_quartic = u * arg
    end if
 
  end function blasius_quartic
 
  function blasius_sin(y, delta, u)
    real(kind=rp), intent(in) :: y, delta, u
    real(kind=rp) :: blasius_sin
    real(kind=rp), parameter :: pi = 4.0_rp * atan(1.0_rp)
    real(kind=rp) :: arg
 
    arg = (pi / 2.0_rp) * (y/delta)
 
    if (arg .gt. 0.5_rp * pi) then
       blasius_sin = u
    else
       blasius_sin = u * sin(arg)
    end if
 
  end function blasius_sin
 
  function blasius_tanh(y, delta, u)
    real(kind=rp), intent(in) :: y, delta, u
    real(kind=rp) :: blasius_tanh
    real(kind=rp), parameter :: a = 0.33245_rp, &
                                n = 5.0_rp/3.0_rp, &
                                eta99 = 5.075
 
    blasius_tanh = u * tanh((a * eta99 * (y/delta))**n)**(1/n)
 
  end function blasius_tanh
 
end module flow_profile