From 33504209a65a9ed4402089ea7ff4e9e846950fa8 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 12:21:22 -0500 Subject: [PATCH] Set up simple disruption model --- src/collision/collision_generate.f90 | 24 +++++++++++++++++------- src/collision/collision_module.f90 | 1 + src/collision/collision_resolve.f90 | 3 ++- 3 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index a3ff0cbca..aa2e79761 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -328,28 +328,38 @@ end subroutine collision_generate_simple_pos_vec module subroutine collision_generate_simple_vel_vec(collider) !! Author: David A. Minton !! - !! Generates an initial "guess" for the velocitity vectors + !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be + !! 2x the escape velocity of the colliding pair. For hit and runs the velocity magnitude is set to be + !! 2x the escape velocity of the smallest of the two bodies. implicit none ! Arguments class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object ! Internals integer(I4B) :: i, j - logical :: lcat - real(DP), dimension(NDIM) :: vimp_unit + logical :: lhr + real(DP), dimension(NDIM) :: vimp_unit, vcom, rnorm real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise real(DP), parameter :: VNOISE_MAG = 0.10_DP real(DP) :: vmag associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) - lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + lhr = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! "Bounce" the first two bodies do i = 1,2 - fragments%vc(:,i) = impactors%vb(:,i) - impactors%vbcom(:) - 2 * impactors%vbimp(:) + vcom(:) = impactors%vb(:,i) - impactors%vbcom(:) + rnorm(:) = impactors%y_unit(:) + ! Do the reflection + vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:) + fragments%vc(:,i) = vcom(:) end do - vmag = .mag.impactors%vbcom(:) - vimp_unit(:) = .unit. impactors%vbimp(:) + ! Compute the escape velocity + if (lhr) then + vmag = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) + else + vmag = 2 * sqrt(2 * sum(impactors%Gmass(:)) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1)))) + end if call random_number(vnoise) vnoise = (2 * vnoise - 1.0_DP) * vmag do i = 3, nfrag diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index ef75b56e9..046341a93 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -58,6 +58,7 @@ module collision real(DP), dimension(NDIM,2) :: Lspin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: Lorbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision + real(DP), dimension(2) :: Gmass !! Two-body equivalent G*mass of the collider bodies prior to the collision real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision real(DP), dimension(2) :: radius !! Two-body equivalent radii of the collider bodies prior to the collision real(DP) :: Qloss !! Energy lost during the collision diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 419875a80..af4e4583c 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -57,7 +57,8 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa pl%kin(idx_parent(2))%parent = idx_parent(1) end if - impactors%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si + impactors%Gmass(:) = pl%Gmass(idx_parent(:)) + impactors%mass(:) = pl%mass(idx_parent(:)) impactors%radius(:) = pl%radius(idx_parent(:)) volume(:) = (4.0_DP / 3.0_DP) * PI * impactors%radius(:)**3