spx_float_object.f90 Source File

!> 浮体粒子集模块
!> @todo 核实浮体的加速度求解、速度、位移更新公式
module spx_float_object

    use spx_kinds, only: rk
    use spx_particle, only: particle_type
    use spx_math, only: dim
    implicit none

    private
    public :: float_object_type

    !> 浮体刚体粒子集数据结构
    type float_object_type
        integer :: istart, iend !! 浮体粒子集的索引
        real(rk) :: mass(2)     !! 浮体质量、转动惯量Izz
        real(rk) :: loc(2)      !! 浮体质心
        real(rk) :: vel(2)      !! 浮体速度
        real(rk) :: acc(2)      !! 浮体加速度
        real(rk) :: ang_vel     !! 浮体角速度
        real(rk) :: ang_acc     !! 浮体角加速度
    contains
        procedure :: init, get_acc, update_particle_init, update_particle_step
    end type float_object_type

contains

    !> 初始化浮体的质量、质心、转动惯量、速度、角速度
    pure subroutine init(self, particle, istart, iend)
        class(float_object_type), intent(inout) :: self  !! 浮体粒子集
        type(particle_type), intent(in) :: particle  !! 区域
        integer, intent(in) :: istart, iend  !! 浮体粒子集的索引
        integer :: i

        self%istart = istart
        self%iend = iend
        associate (mass => particle%mass(self%istart:self%iend), &
                   loc => particle%loc(:, self%istart:self%iend))
            self%mass(1) = sum(mass)
            do i = 1, dim
                self%loc(i) = sum(mass*loc(i, :))/self%mass(1)
            end do
            self%mass(2) = sum(mass*((loc(1, :) - self%loc(1))**2 + &
                                     (loc(2, :) - self%loc(2))**2))
        end associate
        self%vel = 0.0_rk
        self%ang_vel = 0.0_rk

    end subroutine init

    !> 获取浮体受力，更新浮体的加速度
    !> @todo 指向文档公式
    pure subroutine get_acc(self, particles)
        class(float_object_type), intent(inout) :: self  !! 浮体粒子集
        type(particle_type), intent(in) :: particles  !! 区域
        real(rk) :: force(2), torque
        integer :: i

        force = 0.0_rk
        torque = 0.0_rk
        do i = self%istart, self%iend
            force = force + particles%acc(:, i)*particles%mass(i)
            torque = torque + cross_product_local(particles%loc(:, i) - self%loc(:), &
                                                  particles%acc(:, i)*particles%mass(i))
        end do
        self%acc = force/self%mass(1)
        self%ang_acc = torque/self%mass(2)

    end subroutine get_acc

    !> 简单叉乘, 2 维叉乘
    pure function cross_product_local(a, b) result(c)
        real(rk), intent(in), dimension(2) :: a, b
        real(rk) :: c

        c = a(1)*b(2) - a(2)*b(1)

    end function cross_product_local

    !> 初始化leapfrog算法，基于加速度，推进速度半步长
    !> @todo 修复更新运动公式
    pure subroutine update_particle_init(self, particles, dt)
        class(float_object_type), intent(inout) :: self  !! 浮体粒子集
        type(particle_type), intent(inout) :: particles  !! 区域
        real(rk), intent(in) :: dt  !! 时间步长
        integer :: i

        associate (halfdt => dt/2)
            self%vel = self%vel + self%acc*halfdt
            self%ang_vel = self%ang_vel + self%ang_acc*halfdt
        end associate

        ! 更新粒子
        do i = self%istart, self%iend
            particles%vel(:, i) = self%vel(:) + cross_product_local2( &
                                  self%ang_vel, &
                                  particles%loc(:, i) - self%loc(:))
            particles%loc(:, i) = particles%loc(:, i) + particles%vel(:, i)*dt
        end do

        self%loc = self%loc + self%vel*dt

    end subroutine update_particle_init

    !> 更新粒子, 速度比位移、加速度快半步长
    pure subroutine update_particle_step(self, particles, dt)
        class(float_object_type), intent(inout) :: self  !! 浮体粒子集
        type(particle_type), intent(inout) :: particles  !! 区域
        real(rk), intent(in) :: dt                       !! 时间步长
        integer :: i

        ! 浮体表面光滑，加速度由位移决定
        self%vel = self%vel + self%acc*dt
        self%ang_vel = self%ang_vel + self%ang_acc*dt

        ! 更新粒子
        do i = self%istart, self%iend
            particles%vel(:, i) = self%vel(:) + cross_product_local2( &
                                  self%ang_vel, &
                                  particles%loc(:, i) - self%loc(:))
            particles%loc(:, i) = particles%loc(:, i) + particles%vel(:, i)*dt
        end do

        self%loc = self%loc + self%vel*dt

    end subroutine update_particle_step

    !> 简单叉乘 2, @tocheck \\( \omega \cdot r \\)
    pure function cross_product_local2(a, b) result(c)
        real(rk), intent(in) :: a
        real(rk), intent(in), dimension(2) :: b
        real(rk), dimension(2) :: c

        c(1) = -a*b(2)
        c(2) = a*b(1)

    end function cross_product_local2

end module spx_float_object