Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 Joris Vos & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module binary_wind
21 :
22 : use const_def, only: dp, standard_cgrav
23 : use star_lib
24 : use star_def
25 : use math_lib
26 : use binary_def
27 :
28 : implicit none
29 :
30 : contains
31 :
32 0 : subroutine eval_wind_xfer_fractions(binary_id, ierr)
33 : integer, intent(in) :: binary_id
34 : integer, intent(out) :: ierr
35 : type (binary_info), pointer :: b
36 : ierr = 0
37 0 : call binary_ptr(binary_id, b, ierr)
38 0 : if (ierr /= 0) then
39 0 : write(*,*) 'failed in binary_ptr'
40 0 : return
41 : end if
42 :
43 : ! for the primary
44 0 : if (b% point_mass_i /= 1) then
45 0 : if (.not. b% do_wind_mass_transfer_1 .or. b% model_twins_flag) then
46 0 : b% wind_xfer_fraction(1) = 0d0
47 0 : else if(.not. b% use_other_binary_wind_transfer) then
48 0 : call Bondi_Hoyle_wind_transfer(b% binary_id, 1, ierr)
49 0 : if (ierr /=0) then
50 0 : write(*,*) "Error in Bondi_Hoyle_wind_transfer(b% binary_id, 1, ierr)"
51 0 : return
52 : end if
53 : else
54 0 : call b% other_binary_wind_transfer(b% binary_id, 1, ierr)
55 0 : if (ierr /=0) then
56 0 : write(*,*) "Error in other_binary_wind_transfer(b% binary_id, 1, ierr)"
57 0 : return
58 : end if
59 : end if
60 : end if
61 :
62 : ! check if secondary needs wind transfer
63 0 : if (b% point_mass_i /= 2) then
64 0 : if (.not. b% do_wind_mass_transfer_2) then
65 0 : b% wind_xfer_fraction(2) = 0d0
66 0 : else if(.not. b% use_other_binary_wind_transfer) then
67 0 : call Bondi_Hoyle_wind_transfer(b% binary_id, 2, ierr)
68 0 : if (ierr /=0) then
69 0 : write(*,*) "Error in Bondi_Hoyle_wind_transfer(b% binary_id, 2, ierr)"
70 0 : return
71 : end if
72 : else
73 0 : call b% other_binary_wind_transfer(b% binary_id, 2, ierr)
74 0 : if (ierr /=0) then
75 0 : write(*,*) "Error in other_binary_wind_transfer(b% binary_id, 2, ierr)"
76 0 : return
77 : end if
78 : end if
79 : end if
80 :
81 : end subroutine eval_wind_xfer_fractions
82 :
83 0 : subroutine Bondi_Hoyle_wind_transfer(binary_id, s_i, ierr)
84 : integer, intent(in) :: binary_id, s_i ! s_i is index of the wind mass losing star
85 : integer, intent(out) :: ierr
86 :
87 : ! wind transfer fraction based on Bondi-Hoyle mechanism as described in
88 : ! Hurley et al. 2002, MNRAS, 329, 897-928
89 :
90 : type(binary_info), pointer :: b
91 : type (star_info), pointer :: s
92 0 : real(dp) :: v_orb, v_wind
93 0 : real(dp) :: alpha, beta ! Bondi-Hoyle alpha, beta for that star
94 0 : real(dp) :: max_xfer ! Maximum transfer fraction
95 :
96 0 : call binary_ptr(binary_id, b, ierr)
97 0 : if (ierr /= 0) then
98 0 : write(*,*) 'failed in binary_ptr'
99 0 : return
100 : end if
101 :
102 0 : if (s_i == 1) then
103 0 : s => b% s1
104 0 : alpha = b% wind_BH_alpha_1
105 0 : beta = b% wind_BH_beta_1
106 0 : max_xfer = b% max_wind_transfer_fraction_1
107 : else
108 0 : s => b% s2
109 0 : alpha = b% wind_BH_alpha_2
110 0 : beta = b% wind_BH_beta_2
111 0 : max_xfer = b% max_wind_transfer_fraction_2
112 : end if
113 :
114 : ! orbital speed Hurley et al 2002 eq. 8
115 0 : v_orb = sqrt(standard_cgrav * (b% m(1) + b% m(2)) / b% separation) ! cm/s
116 :
117 : ! windspeed from Hurley et al 2002 eq. 9
118 0 : v_wind = sqrt(2d0 * beta * standard_cgrav * b% m(s_i) / b% r(s_i))
119 :
120 : ! Bondi-Hoyle transfer fraction Hurley et al. 2002 eq. 6
121 : b% wind_xfer_fraction(s_i) = alpha / pow2(b% separation) / &
122 : (2d0 * sqrt(1d0 - pow2(b% eccentricity))) * &
123 : pow2(standard_cgrav * b% m(3-s_i) / pow2(v_wind)) * &
124 0 : pow(1d0 + pow2(v_orb/v_wind),-1.5d0)
125 :
126 : ! limit to provided maximum
127 0 : b% wind_xfer_fraction(s_i) = min(max_xfer, b% wind_xfer_fraction(s_i))
128 :
129 : end subroutine Bondi_Hoyle_wind_transfer
130 :
131 0 : subroutine Tout_enhance_wind(b, s)
132 : type (binary_info), pointer :: b
133 : type (star_info), pointer :: s
134 :
135 : ! Tidaly enhance wind mass loss as described by
136 : ! Tout & Eggleton 1988,MNRAS,231,823 (eq. 2)
137 0 : real(dp) :: B_wind ! enhancement parameter, B in eq. 2
138 : integer :: i, s_i
139 0 : real(dp) :: dm
140 0 : real(dp), DIMENSION(b% anomaly_steps):: rl_d, r_rl, mdot
141 :
142 0 : if (s% id == b% s1% id) then
143 0 : if (.not. b% do_enhance_wind_1) return
144 0 : B_wind = b% tout_B_wind_1
145 0 : s_i = 1
146 : else
147 0 : if (.not. b% do_enhance_wind_2) return
148 0 : B_wind = b% tout_B_wind_2
149 0 : s_i = 2
150 : end if
151 :
152 0 : do i = 1,b% anomaly_steps !limit radius / roche lobe
153 : ! phase dependent roche lobe radius
154 0 : rl_d(i) = (1d0 - pow2(b%eccentricity)) / (1d0 + b%eccentricity*cos(b% theta_co(i))) * b% rl(s_i)
155 0 : r_rl(i) = min(pow6(b% r(s_i) / rl_d(i)), pow6(0.5d0))
156 : end do
157 :
158 : ! actual enhancement
159 0 : mdot = s% mstar_dot * (1d0 + B_wind * r_rl)
160 :
161 : dm = 0d0
162 0 : do i = 2,b% anomaly_steps ! trapezoidal integration
163 0 : dm = dm + 0.5d0 * (mdot(i-1) + mdot(i)) * (b% time_co(i) - b% time_co(i-1))
164 : end do
165 :
166 : ! remember mass-loss is negative!
167 : !b% mdot_wind_theta = b% mdot_wind_theta + mdot ! store theta dependance for edot
168 0 : s% mstar_dot = dm ! return enhanced wind mass loss
169 :
170 0 : end subroutine Tout_enhance_wind
171 :
172 : end module binary_wind
173 :
|
