!> author: 左志华 !> date: 2022-07-02 !> !> Internal forces <br> !> 内力：粘性力、压力 !> @todo 找到解决初始时刻边界压力和排斥力扰动问题 module spx_internal_force use spx_kinds, only: rk use spx_particle, only: particle_type use spx_nnps_pairs, only: nnps_pairs_type use spx_eos, only: p_ideal_gas, eow_obj implicit none procedure(internal_force_fcn), pointer :: internal_force => internal_force_with_visc !! Internal force function pointer <br> !! 内力函数指针 logical :: viscosity !! Viscosity on/off <br> !! 是否考虑粘性力 abstract interface !> Internal force function <br> !> 内力函数 subroutine internal_force_fcn(particle, pairs, acc, n) import :: particle_type, nnps_pairs_type, rk type(particle_type), intent(inout) :: particle !! 区域 type(nnps_pairs_type), intent(in) :: pairs !! 粒子对 real(rk), intent(inout) :: acc(:, :) !! 加速度 integer, intent(in) :: n !! 实弹性粒子数, numbers(1) end subroutine internal_force_fcn end interface contains !> Internal force function without viscosity <br> !> 内力求解 !> @note 边界与浮体粒子都是光滑壁面，不参与内力计算 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 !> Internal force function with viscosity <br> !> 内力求解（含粘性力） !> @todo 粘性熵、温度变化 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 end module spx_internal_force