Fix wind accretion (#477)

* Fix to the orbital velocity in binary_wind.f90.
also add support for different BH_beta according to hurley+2002

* modify beta in the HMXB testcase

* fix assignement of wind_BH_beta controls

* adjust inlist of wind_fed_bhhmxb to use old value of wind_BH_beta_1

---------

Co-authored-by: Matthias Fabry <[email protected]>
Co-authored-by: Pablo Marchant <[email protected]>

Author: 3 people

binary/defaults/binary_controls.defaults

@@ -813,15 +813,19 @@
 ! wind_BH_alpha_*
 ! ~~~~~~~~~~~~~~~
 
-! Bondi-Hoyle accretion parameter for each star. The default is 3/2 taken
-! from Hurley et al. 2002, MNRAS, 329, 897, in agreement with
-! Boffin & Jorissen 1988, A&A, 205, 155
+! Bondi-Hoyle accretion parameter for each star. The default for alpha is 3/2
+! taken from Hurley et al. 2002, MNRAS, 329, 897, in agreement with
+! Boffin & Jorissen 1988, A&A, 205, 155.
+! The default for beta is 1/8=0.125 in accordance for results of cool
+! supergiants from Kucinskas A., 1999, Ap&SS, 262, 127
 ! "\_1" refers to first star, "\_2" to the second one.
 
 ! ::
 
     wind_BH_alpha_1 = 1.5d0
     wind_BH_alpha_2 = 1.5d0
+    wind_BH_beta_1 = 1.25d-1
+    wind_BH_beta_2 = 1.25d-1
 
 
 ! max_wind_transfer_fraction_*

binary/private/binary_ctrls_io.f90

@@ -125,6 +125,8 @@ module binary_ctrls_io
          do_wind_mass_transfer_2, &
          wind_BH_alpha_1, &
          wind_BH_alpha_2, &
+         wind_BH_beta_1, &
+         wind_BH_beta_2, &
          max_wind_transfer_fraction_1, &
          max_wind_transfer_fraction_2, &
 
@@ -428,6 +430,8 @@ subroutine store_binary_controls(b, ierr)
          b% do_wind_mass_transfer_2 = do_wind_mass_transfer_2
          b% wind_BH_alpha_1 = wind_BH_alpha_1
          b% wind_BH_alpha_2 = wind_BH_alpha_2
+         b% wind_BH_beta_1 = wind_BH_beta_1
+         b% wind_BH_beta_2 = wind_BH_beta_2
          b% max_wind_transfer_fraction_1 = max_wind_transfer_fraction_1
          b% max_wind_transfer_fraction_2 = max_wind_transfer_fraction_2
 
@@ -618,6 +622,8 @@ subroutine set_binary_controls_for_writing(b, ierr)
          do_wind_mass_transfer_2 = b% do_wind_mass_transfer_2
          wind_BH_alpha_1 = b% wind_BH_alpha_1
          wind_BH_alpha_2 = b% wind_BH_alpha_2
+         wind_BH_beta_1 = b% wind_BH_beta_1
+         wind_BH_beta_2 = b% wind_BH_beta_2
          max_wind_transfer_fraction_1 = b% max_wind_transfer_fraction_1
          max_wind_transfer_fraction_2 = b% max_wind_transfer_fraction_2
 

binary/private/binary_wind.f90

@@ -94,8 +94,8 @@ subroutine Bondi_Hoyle_wind_transfer(binary_id, s_i, ierr)
 
       type(binary_info), pointer :: b
       type (star_info), pointer :: s
-      real(dp) :: v_orb, v_wind, b_BH
-      real(dp) :: alpha  ! Bondi-Hoyle alpha for that star
+      real(dp) :: v_orb, v_wind
+      real(dp) :: alpha, beta  ! Bondi-Hoyle alpha, beta for that star
       real(dp) :: max_xfer  ! Maximum transfer fraction
 
       call binary_ptr(binary_id, b, ierr)
@@ -107,18 +107,20 @@ subroutine Bondi_Hoyle_wind_transfer(binary_id, s_i, ierr)
       if (s_i == 1) then
          s => b% s1
          alpha = b% wind_BH_alpha_1
+         beta = b% wind_BH_beta_1
          max_xfer = b% max_wind_transfer_fraction_1
       else
          s => b% s2
          alpha = b% wind_BH_alpha_2
+         beta = b% wind_BH_beta_2
          max_xfer = b% max_wind_transfer_fraction_2
       end if
 
       ! orbital speed Hurley et al 2002 eq. 8
-      v_orb = sqrt(standard_cgrav * b% m(s_i) / b% separation) !cm/s
+      v_orb = sqrt(standard_cgrav * (b% m(1) + b% m(2)) / b% separation) !cm/s
 
       ! windspeed from Hurley et al 2002 eq. 9
-      v_wind = sqrt( 2d0 / 8d0 *  standard_cgrav * b% m(s_i) / b% r(s_i) )
+      v_wind = sqrt(2d0 * beta *  standard_cgrav * b% m(s_i) / b% r(s_i))
 
       ! Bondi-Hoyle transfer fraction Hurley et al. 2002 eq. 6
       b% wind_xfer_fraction(s_i) = alpha / pow2(b% separation) /&

binary/public/binary_controls.inc

@@ -84,6 +84,7 @@ logical :: do_enhance_wind_1, do_enhance_wind_2
 real(dp) :: tout_B_wind_1, tout_B_wind_2
 logical :: do_wind_mass_transfer_1, do_wind_mass_transfer_2
 real(dp) :: wind_BH_alpha_1, wind_BH_alpha_2
+real(dp) :: wind_BH_beta_1, wind_BH_beta_2
 real(dp) :: max_wind_transfer_fraction_1, max_wind_transfer_fraction_2
 
 ! orbital jdot controls

binary/test_suite/wind_fed_bhhmxb/inlist_project

@@ -22,6 +22,8 @@
 
    limit_retention_by_mdot_edd = .true.
    do_wind_mass_transfer_1 = .true.
+   wind_BH_beta_1 = 1.25d-1 ! wind_BH_beta_1 = 7d0  ! 7 since donor is O star, using this value would make more sense, but we keep the default of 1/8
+                           ! which results in wind accretion below and above the Eddington limit (which we want to test).
    do_jdot_mb = .false.
 
    report_rlo_solver_progress = .true.