ic_2dcrvpPressure_for Function

public function ic_2dcrvpPressure_for(me, coord, n) result(pressure)

This function defines the density of the spinning (= co-rotating) vortex pair See the matlab file where the pressure is plot in the ase-testcases/ repo in musubi/crvp/matlab/crvp_velPress_plot.m

As a reference, see [1] R. Fortenbach, 'Mehrskalenmodellierung von aeroakustischen Quellen in schwach kompressiblen Stroemungen,' pp. 1-151, Jul. 2006.

Source: complex velocity potential of both vortices complex coordinates: z = x+i*y Gamma ... circulation

 b = r0*exp(i*omega*t) ... rotation orbit
 w(z,t) = Gamma/(2Pi*i)*ln(z^2-b^2) ... potential function
 dw/dz = Gamma/(2Pi*i)*z/(z^2-b^2)  ... derivative of potential
 u =  Re( dw/dz( z, t=0 )           ... x -velocity
 v = -Im( dw/dz( z, t=0 )           ... y -velocity
 u0 = Gamma/(4Pi*r0)                ... rotation velocity at center of each
                                        vortice
 omega = 2Pi/T0 = u0/r0 = Gamma/(4Pi*ro^2) ... rotation angular frequency
 Ma = u0/cs        ... rotation Mach number
 rho = rho0 - Ma^2*rho0/cs^2*( Re{ -omega*Gamma/Pi * b^2/(z^2-b^2)}
                                                   + (u^2+v^2)/2 )

Unit of the result is in kg/m^3, as the coordinates are given in physical coordinates and hence all other parameters also have to be physical ones

Arguments

Type IntentOptional Attributes Name
type(ic_2dcrvp_type), intent(in) :: me

global gauss pulse data
real(kind=rk), intent(in) :: coord(n,3)

coordinate of an element
integer, intent(in) :: n

number of return values

Return Value real(kind=rk), (n)

return value which is sent to state variable


Calls

proc~~ic_2dcrvppressure_for~~CallsGraph proc~ic_2dcrvppressure_for ic_2dcrvpPressure_for proc~ic_2dcrvpx_for ic_2dcrvpX_for proc~ic_2dcrvppressure_for->proc~ic_2dcrvpx_for proc~cutoff_factor cutoff_factor proc~ic_2dcrvppressure_for->proc~cutoff_factor proc~ic_2dcrvpy_for ic_2dcrvpY_for proc~ic_2dcrvppressure_for->proc~ic_2dcrvpy_for proc~ic_2dcrvpx_for->proc~cutoff_factor proc~ic_2dcrvpy_for->proc~cutoff_factor

Called by

proc~~ic_2dcrvppressure_for~~CalledByGraph proc~ic_2dcrvppressure_for ic_2dcrvpPressure_for proc~tem_spatial_for_coord tem_spatial_for_coord proc~tem_spatial_for_coord->proc~ic_2dcrvppressure_for proc~tem_spatial_scalar_for_index tem_spatial_scalar_for_index proc~tem_spatial_scalar_for_index->proc~tem_spatial_for_coord interface~tem_spatial_for tem_spatial_for interface~tem_spatial_for->proc~tem_spatial_for_coord interface~tem_spatial_for->proc~tem_spatial_scalar_for_index proc~tem_spacetime_for_treeids tem_spacetime_for_treeIDs proc~tem_spacetime_for_treeids->interface~tem_spatial_for proc~tem_spacetime_vector_for_index tem_spacetime_vector_for_index proc~tem_spacetime_vector_for_index->interface~tem_spatial_for proc~tem_spacetime_vector_for_coord tem_spacetime_vector_for_coord proc~tem_spacetime_vector_for_index->proc~tem_spacetime_vector_for_coord proc~tem_spatial_scalar_storeval tem_spatial_scalar_storeVal proc~tem_spatial_scalar_storeval->interface~tem_spatial_for proc~tem_spacetime_for_coord tem_spacetime_for_coord proc~tem_spacetime_for_coord->interface~tem_spatial_for proc~tem_spacetime_scalar_for_index tem_spacetime_scalar_for_index proc~tem_spacetime_scalar_for_index->interface~tem_spatial_for proc~tem_spacetime_scalar_for_index->proc~tem_spacetime_for_coord proc~tem_spacetime_vector_for_treeids tem_spacetime_vector_for_treeIDs proc~tem_spacetime_vector_for_treeids->interface~tem_spatial_for proc~tem_spacetime_vector_for_coord->interface~tem_spatial_for proc~tem_spatial_vector_storeval tem_spatial_vector_storeVal proc~tem_spatial_vector_storeval->interface~tem_spatial_for interface~tem_spacetime_for tem_spacetime_for interface~tem_spacetime_for->proc~tem_spacetime_for_treeids interface~tem_spacetime_for->proc~tem_spacetime_vector_for_index interface~tem_spacetime_for->proc~tem_spacetime_for_coord interface~tem_spacetime_for->proc~tem_spacetime_scalar_for_index interface~tem_spacetime_for->proc~tem_spacetime_vector_for_treeids interface~tem_spacetime_for->proc~tem_spacetime_vector_for_coord proc~tem_spacetime_for_stcoord tem_spacetime_for_stcoord proc~tem_spacetime_for_stcoord->proc~tem_spacetime_for_coord interface~tem_spatial_storeval tem_spatial_storeVal interface~tem_spatial_storeval->proc~tem_spatial_scalar_storeval interface~tem_spatial_storeval->proc~tem_spatial_vector_storeval

Contents

Source Code

ic_2dcrvpPressure_for

Source Code

  function ic_2dcrvpPressure_for( me, coord, n ) result(pressure)
    ! ---------------------------------------------------------------------------
    !> number of return values
    integer, intent(in) :: n
    !> global gauss pulse data
    type(ic_2dcrvp_type), intent(in) :: me
    !> coordinate of an element
    real(kind=rk), intent(in) :: coord(n, 3)
    !> return value which is sent to state variable
    real(kind=rk) :: pressure(n)
    ! ---------------------------------------------------------------------------
    complex(kind=rk) :: z_coord, b_coord, z_compl, b_coordChg, z_Core
    real(kind=rk) :: omega
    real(kind=rk) :: u(1),v(1)
    complex(kind=rk),parameter :: i_complex = (0._rk,1._rk)
    real(kind=rk) :: r, rC
    real(kind=rk) :: dist2Left, dist2Right, pMin, phiTC, uxC, uyC
    integer :: iElem
    ! ---------------------------------------------------------------------------

    do iElem = 1, n
      ! complex coordinate
      z_coord    = ( coord(iElem,1) - me%center(1))                            &
        &        + i_complex*(coord(iElem,2) - me%center(2))
      z_compl    = ( coord(iElem,1) - me%center(1))                            &
        &        - i_complex*(coord(iElem,2) - me%center(2))

      omega = me%circulation/(Pi*me%radius_rot*me%radius_rot*4._rk)
      ! Current radial coordinate
      r = abs( z_coord )
      rC = me%radius_C
      b_coord    = me%radius_rot*( cos( omega*me%t )                           &
        &        + i_complex*sin( omega*me%t))
      b_coordChg = b_coord ! *3._rk
      z_Core = b_coord + rC

      ! Calculate velocity
      u(:)   = ic_2dcrvpX_for( me, coord( iElem, :), 1 )
      v(:)   = ic_2dcrvpY_for( me, coord( iElem, :), 1 )

      ! Gaussian pulse approimxation inside core radius
      ! Compute the distance^2
      dist2Left  = real((z_coord + b_coord)*(z_compl + b_coordChg), kind=rk)
      dist2Right = real(((z_coord) - b_coord)*(z_compl-b_coordChg), kind=rk)
      if( dist2Right .lt. rC*rC .or. dist2Left .lt. rC*rC ) then
        uxC = real( me%circulation/(2._rk*PI*i_complex*(z_Core-b_coord)),      &
          &                                                            kind=rk)
        uyC = real(-aimag( me%circulation/(2._rk*PI*i_complex*                 &
          &                                        (z_Core-b_coord))), kind=rk)
        phiTC = real( -omega*me%circulation/PI*b_coord**2/                     &
          &                                   (z_Core**2-b_coord**2), kind=rk )
        pMin = me%rho0*( phiTC + (uxC**2+uyC**2)*0.5_rk);
        pMin = pMin * me%matchFactor
      end if

      if( dist2Right .lt. rC*rC ) then
        dist2Right = -0.5_rk/(rC*rC)*dist2Right
        if( me%pressGaussModel ) then
          ! velocity at the core radius: real(b_coord) + rC
          ! right vortex
          pressure( iElem ) = me%p0 - cutoff_factor( me%cutoff, r )            &
            &               * pMin *exp( dist2Right )
        else
          pressure( iElem ) = me%p0
        end if

      elseif( dist2Left .lt. rC*rC ) then
        dist2Left  = -0.5_rk/(rC*rC)*dist2Left
        if( me%pressGaussModel ) then
          ! velocity at the core radius: real(b_coord) + rC
          pressure( iElem ) = me%p0 - cutoff_factor( me%cutoff, r )            &
            &               * pMin*exp( dist2Left )
        else
          pressure( iElem ) = me%p0
        end if
      else
        pressure( iElem ) = ( - cutoff_factor(me%cutoff, r)                    &
          &               * ( real( -omega*me%circulation/PI*b_coord*b_coord/  &
          &                                 (z_coord*z_coord-b_coord*b_coord)) &
          &               + (u(1)*u(1) + v(1)*v(1))*0.5_rk ))*me%rho0 + me%p0
      end if
    end do

  end function ic_2dcrvpPressure_for