spx_artificial_force.f90 Source File

!> author: 左志华
!> date: 2022-07-05
!>
!> Artificial viscosity and heat <br>
!> 人工热量与人工粘性
!> @todo 人工热量，人工热量影响内能
module spx_artifical_force

    use spx_kinds, only: rk
    use spx_particle, only: particle_type
    use spx_nnps_pairs, only: nnps_pairs_type
    implicit none

    real(rk), parameter :: alpha = 1.0_rk   !! 剪切粘度
    real(rk), parameter :: beta = 1.0_rk    !! 容积（bulk）粘度
    real(rk), parameter :: etq = 0.1_rk     !! 避免奇点

contains

    !> Artificial viscosity <br>
    !> 人工粘度
    !> @todo 添加测试
    !> @note 边界与浮体粒子都是光滑壁面，不参与人工力计算
    pure subroutine artificial_viscosity(particle, pairs, acc, n)
        type(particle_type), intent(inout) :: particle     !! particle <br>
                                                    !! 计算域
        type(nnps_pairs_type), intent(in) :: pairs       !! Pair vector <br>
                                                    !! 粒子对
        real(rk), intent(inout) :: acc(:, :)        !! Visc acceleration <br>
                                                    !! 粘性力（加速度）
        integer, intent(in) :: n                    !! 实弹性粒子数, numbers(1)
        real(rk) :: vr, rr, dvel(2), h(2)
        real(rk) :: muv, mc, mrho, piv
        integer :: i

        do i = 1, pairs%len
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i))

                if (pairs%contact_type(i) /= 2) cycle  ! 仅弹性实粒子参与内力计算
                dvel(:) = particle%vel(:, ip) - particle%vel(:, jp)

                associate (dx => pairs%rdx(3*i - 1:3*i))
                    vr = sum(dvel(:)*dx(:))
                    rr = sum(dx(:)*dx(:))
                end associate

                if (vr < 0) then
                    muv = particle%hsml*vr/(rr + (particle%hsml*etq)**2)
                    mc = (particle%c(ip) + particle%c(jp))/2
                    mrho = (particle%rho(ip) + particle%rho(jp))/2
                    piv = (beta*muv - alpha*mc)*muv/mrho

                    h(:) = -piv*pairs%wdwdx(2:, i)
                    acc(:, ip) = acc(:, ip) + h(:)*particle%mass(jp)
                    acc(:, jp) = acc(:, jp) - h(:)*particle%mass(ip)
                end if

            end associate
        end do

    end subroutine artificial_viscosity

end module spx_artifical_force