Internal force function without viscosity
内力求解
边界与浮体粒子都是光滑壁面,不参与内力计算
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(particle_type), | intent(inout) | :: | particle |
区域 |
||
type(nnps_pairs_type), | intent(in) | :: | pairs |
粒子对 |
||
real(kind=rk), | intent(inout) | :: | acc(:,:) |
加速度 |
||
integer, | intent(in) | :: | n |
实弹性粒子数, numbers(1) |
pure subroutine internal_force_without_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) :: h(2)
! ----------------------------- 根据 EoS 计算压强信息 -------------------------------- !
do concurrent(i=1:n)
select case (particle%itype(i)) ! TOCHECK: 虚粒子是否需要基于 EoS 计算压强
case (50:59); call eow_obj%p_water(particle%rho(i), particle%p(i), particle%c(i)) ! 淡水
case (110:119); call p_ideal_gas(particle%rho(i), particle%u(i), particle%p(i), particle%c(i)) ! 理想气体
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 ! 仅弹性实粒子参与内力计算
! 用到了虚粒子的压强对实粒子的影响,防止边界处实粒子有较大的、指向边界的初始压强合力
h(:) = -(particle%p(ip) + particle%p(jp))*pairs%wdwdx(2:3, i)/ &
&(particle%rho(ip)*particle%rho(jp))
acc(:, ip) = acc(:, ip) + h(:)*particle%mass(jp)
acc(:, jp) = acc(:, jp) - h(:)*particle%mass(ip)
end associate
end do
end subroutine internal_force_without_visc