artificial_viscosity Subroutine

    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