internal_force_with_visc Subroutine

    pure subroutine internal_force_with_visc(particle, pairs, acc, n)
        type(particle_type), intent(inout) :: particle  !! 区域
        type(nnps_pairs_type), intent(in) :: pairs   !! 粒子对
        real(rk), intent(inout) :: acc(:, :)    !! 更新加速度
        integer, intent(in) :: n                !! 实弹性粒子数, numbers(1)
        integer :: i
        real(rk), parameter :: r23 = 2.0_rk/3.0_rk  !! 2/3
        real(rk) :: h(2), dvel(2), hxx, hxy, hyy
        real(rk), dimension(n) :: txx, txy, tyy, eta  ! @todo 分配内存消耗少量时间?

        txx(:) = 0.0_rk; txy(:) = 0.0_rk; tyy(:) = 0.0_rk    ! 初始化本地数组

        ! ----------------------------- 根据 EoS 计算压强等信息 -------------------------------- !
        do concurrent(i=1:n)

            select case (particle%itype(i))
            case (50:59)    ! 淡水
                call eow_obj%p_water(particle%rho(i), particle%p(i), particle%c(i))
                eta(i) = 1.0e-3_rk  ! 运动学粘性系数
            case (110:119)    ! 理想气体
                call p_ideal_gas(particle%rho(i), particle%u(i), particle%p(i), particle%c(i))
                eta(i) = 0.0_rk
            end select

        end do

        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(:, jp) - particle%vel(:, ip)
                hxx = r23*(2*dvel(1)*pairs%wdwdx(2, i) - dvel(2)*pairs%wdwdx(3, i))
                hxy = r23*(dvel(1)*pairs%wdwdx(3, i) + dvel(2)*pairs%wdwdx(2, i))
                hyy = r23*(2*dvel(2)*pairs%wdwdx(3, i) - dvel(1)*pairs%wdwdx(2, i))

                txx(ip) = txx(ip) + particle%mass(jp)*hxx/particle%rho(jp)
                txx(jp) = txx(jp) + particle%mass(ip)*hxx/particle%rho(ip)
                txy(ip) = txy(ip) + particle%mass(jp)*hxy/particle%rho(jp)
                txy(jp) = txy(jp) + particle%mass(ip)*hxy/particle%rho(ip)
                tyy(ip) = tyy(ip) + particle%mass(jp)*hyy/particle%rho(jp)
                tyy(jp) = tyy(jp) + particle%mass(ip)*hyy/particle%rho(ip)

                ! ----------------------------- 计算压力、剪切力 ------------------------------ !
                ! 用到了虚粒子的压强对实粒子的影响，防止边界处实粒子有较大的、指向边界的初始压强合力
                ! 采用 \( (p_i + p_j)/(rho_i*rho_j) \) 方式
                h(:) = -(particle%p(ip) + particle%p(jp))*pairs%wdwdx(2:, i)
                ! he = sum((particle%vel(:, jp) - particle%vel(:, ip))*h(:))
                h(1) = h(1) + &
                   &(eta(ip)*txx(ip) + eta(jp)*txx(jp))*pairs%wdwdx(2, i) + &
                    (eta(ip)*txy(ip) + eta(jp)*txy(jp))*pairs%wdwdx(3, i)
                h(2) = h(2) + &
                   &(eta(ip)*txy(ip) + eta(jp)*txy(jp))*pairs%wdwdx(2, i) + &
                    (eta(ip)*tyy(ip) + eta(jp)*tyy(jp))*pairs%wdwdx(3, i)

                associate (rho_ij => particle%rho(ip)*particle%rho(jp))
                    acc(:, ip) = acc(:, ip) + h(:)*particle%mass(jp)/rho_ij
                    acc(:, jp) = acc(:, jp) - h(:)*particle%mass(ip)/rho_ij
                end associate

            end associate
        end do

    end subroutine internal_force_with_visc