Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-2019 Pablo Marchant & 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_tides
21 :
22 : use star_lib
23 : use star_def
24 : use const_def, only: dp, pi, pi4, secyer, rsun, msun, one_third, standard_cgrav, convective_mixing
25 : use utils_lib
26 : use math_lib
27 : use binary_def
28 :
29 : implicit none
30 :
31 :
32 : contains
33 :
34 :
35 0 : subroutine sync_spin_orbit_torque(id, ierr)
36 : integer, intent(in) :: id
37 : integer, intent(out) :: ierr
38 : type (star_info), pointer :: s
39 0 : real(dp) :: osep ! orbital separation (cm)
40 0 : real(dp) :: qratio ! mass_other_star/mass_this_star
41 0 : real(dp) :: rlr ! roche lobe radius (cm)
42 0 : real(dp) :: dt_next ! next timestep
43 0 : real(dp) :: Ftid ! efficiency of tidal synchronization. (time scale × FSYNC).
44 : character (len=strlen) :: sync_type
45 : character (len=strlen) :: sync_mode
46 : type (binary_info), pointer :: b
47 : ierr = 0
48 :
49 0 : call star_ptr(id, s, ierr)
50 0 : if (ierr /= 0) then
51 0 : write(*,*) 'failed in star_ptr'
52 0 : return
53 : end if
54 :
55 0 : call binary_ptr(s% binary_id, b, ierr)
56 0 : if (ierr /= 0) then
57 0 : write(*,*) 'failed in binary_ptr'
58 0 : return
59 : end if
60 :
61 0 : if (.not. b% do_tidal_sync) return
62 :
63 0 : call star_ptr(id, s, ierr)
64 :
65 0 : osep = b% separation
66 0 : qratio = b% m(b% a_i) / b% m(b% d_i)
67 0 : if (is_donor(b, s)) then
68 0 : rlr = b% rl(b% d_i)
69 : else
70 0 : qratio = 1d0/qratio
71 0 : rlr = b% rl(b% a_i)
72 : end if
73 0 : dt_next = s% dt
74 :
75 0 : if (b% point_mass_i /= 1 .and. id == b% s1% id) then
76 0 : Ftid = b% Ftid_1
77 0 : sync_type = b% sync_type_1
78 0 : sync_mode = b% sync_mode_1
79 : else
80 0 : Ftid = b% Ftid_2
81 0 : sync_type = b% sync_type_2
82 0 : sync_mode = b% sync_mode_2
83 : end if
84 :
85 0 : if (b% use_other_sync_spin_to_orbit) then
86 0 : call b% other_sync_spin_to_orbit(s% id, s% nz, osep, qratio, rlr, dt_next, Ftid, sync_type, sync_mode, ierr)
87 : else
88 0 : call sync_spin_to_orbit(s% id, s% nz, osep, qratio, rlr, dt_next, Ftid, sync_type, sync_mode, ierr)
89 : end if
90 :
91 : end subroutine sync_spin_orbit_torque
92 :
93 0 : subroutine sync_spin_to_orbit(id, nz, osep, qratio, rl, dt_next, Ftid, sync_type, sync_mode, ierr)
94 : ! initially based on spiba.f kindly provided by Norbert Langer and group.
95 : integer, intent(in) :: id
96 : integer, intent(in) :: nz
97 : real(dp), intent(in) :: osep ! orbital separation (cm)
98 : real(dp), intent(in) :: qratio ! mass_other_star/mass_this_star
99 : real(dp), intent(in) :: rl ! roche lobe radius (cm)
100 : real(dp), intent(in) :: dt_next ! next timestep
101 : real(dp), intent(in) :: Ftid ! efficiency of tidal synchronization. (time scale / Ftid ).
102 :
103 : character (len=strlen), intent(in) :: sync_type ! synchronization timescale
104 : character (len=strlen), intent(in) :: sync_mode ! where to put/take angular momentum
105 : integer, intent(out) :: ierr
106 :
107 : type (star_info), pointer :: s
108 0 : real(dp) :: G, m, t_sync, r_phot, delta_total_J, &
109 0 : sum_J_sync, sum_J_non_sync, tdyn, tkh, rho_face, cv_face, &
110 0 : T_face, csound_face, ff, omega_orb
111 :
112 0 : real(dp), dimension(nz) :: j_sync, delta_j, tdyn_div_tkh
113 0 : integer, dimension(nz) :: layers_in_sync
114 : integer :: k, num_sync_layers
115 : type (binary_info), pointer :: b
116 :
117 0 : real(dp) :: a1,a2
118 :
119 : include 'formats'
120 :
121 : ierr = 0
122 :
123 0 : call star_ptr(id, s, ierr)
124 0 : if (ierr /= 0) then
125 0 : write(*,*) 'failed in star_ptr'
126 0 : return
127 : end if
128 :
129 0 : if (osep <= 0) return
130 0 : if (qratio <= 0) return
131 0 : if (rl <= 0) return
132 0 : if (dt_next <= 0) return
133 0 : if (Ftid <= 0) return
134 0 : if (sync_type == "None") return
135 :
136 0 : call binary_ptr(s% binary_id, b, ierr)
137 0 : if (ierr /= 0) then
138 0 : write(*,*) 'failed in binary_ptr'
139 0 : return
140 : end if
141 :
142 0 : t_sync = 0
143 :
144 0 : G = standard_cgrav
145 :
146 0 : if (is_donor(b, s)) then
147 0 : m = b% m(b% d_i)
148 0 : r_phot = b% r(b% d_i)
149 : else
150 0 : m = b% m(b% a_i)
151 0 : r_phot = b% r(b% a_i)
152 : end if
153 :
154 0 : omega_orb = 2d0*pi/b% period
155 0 : do k=1,nz
156 0 : j_sync(k) = omega_orb*s% i_rot(k)% val
157 : end do
158 :
159 0 : if (sync_type == "Instantaneous") then ! instantaneous synchronisation
160 0 : do k=1,nz
161 0 : delta_j(k) = s% j_rot(k) - j_sync(k)
162 : end do
163 : else
164 : ! get synchronization timescale, i.e.
165 : ! 1/t_sync = \dot(omega_star)/(omega_orb-omega_star)
166 : ! to order e (e=eccentricity).
167 0 : if (.not. b% use_other_tsync) then
168 0 : call get_tsync(s% id, sync_type, Ftid, qratio, m, r_phot, osep, t_sync, ierr)
169 0 : if (ierr/=0) return
170 : else
171 0 : call b% other_tsync(s% id, sync_type, Ftid, qratio, m, r_phot, osep, t_sync, ierr)
172 0 : if (ierr/=0) return
173 : end if
174 0 : if (sync_mode == "Uniform") then
175 0 : a1 = f2(b% eccentricity)
176 0 : a2 = pow(1-pow2(b% eccentricity), 1.5d0)*f5(b% eccentricity)
177 0 : do k=1,nz
178 0 : delta_j(k) = (1d0 - exp(-a2*dt_next/t_sync))*(s% j_rot(k) - a1/a2*j_sync(k))
179 : end do
180 0 : else if (sync_mode == "tdyn_div_tkh") then
181 : ! set local timescale following Gautschy & Glatzel 1990 MNRAS 245, 597-613.
182 0 : do k=1,nz
183 : !tdyn_div_tkh := local dynamical time-scale / local thermal time-scale
184 0 : rho_face = star_interp_val_to_pt(s% rho,k,nz,s% dq,'binary_tides')
185 0 : cv_face = star_interp_val_to_pt(s% cv,k,nz,s% dq,"binary_tides")
186 0 : T_face = star_interp_val_to_pt(s% T,k,nz,s% dq,"binary_tides")
187 0 : csound_face = star_interp_val_to_pt(s% csound,k,nz,s% dq,"binary_tides")
188 0 : tkh = pi4*s% r(k)*s% r(k)*rho_face*cv_face*T_face/s% L(k) ! (4.4)
189 0 : tdyn = 1/csound_face
190 0 : tdyn_div_tkh(k) = tdyn/tkh
191 : end do
192 :
193 : ! j_sync = synchronized specific angular momentum
194 0 : delta_total_J = 0
195 0 : do k=1,nz
196 0 : delta_total_J = delta_total_J + abs(s% j_rot(k) - j_sync(k))*s% dm_bar(k)
197 0 : delta_j(k) = 0
198 0 : layers_in_sync(k) = 0
199 : end do
200 0 : delta_total_J = delta_total_J*(1d0 - exp(-dt_next/t_sync))
201 :
202 : ! Iteratively solve the scaling factor ff to add (or remove) delta_total_J.
203 : ! At each iteration, ff is solved such that each zone k has a change on
204 : ! its angular momentum J_k of the form:
205 : !
206 : ! Delta J_k = (J_k-J_{k,sync})*tdyn_div_tkh(k)**2*ff,
207 : !
208 : ! and the total change adds up to delta_total_J.
209 : ! Since tides can at most drive a layer into sync, ff cannot be
210 : ! solved right away. Then, each iteration solves ff ignoring layers
211 : ! going over sync, spreads the angular momentum to see which ones
212 : ! become synced, and then recalculates taking this into account
213 : ! until it converges.
214 0 : num_sync_layers = 0
215 : !write(*,*) "Entering tide iteration", delta_total_J
216 : do while (.true.)
217 0 : sum_J_sync = 0
218 0 : sum_J_non_sync = 0
219 0 : do k=1,nz
220 0 : delta_j(k) = s% j_rot(k) - j_sync(k)
221 0 : if (layers_in_sync(k) /= 1) then
222 : sum_J_non_sync = sum_J_non_sync + &
223 0 : s% dm_bar(k) * abs(delta_j(k)) * pow2(tdyn_div_tkh(k))
224 : else
225 0 : sum_J_sync = sum_J_sync + s% dm_bar(k) * abs(delta_j(k))
226 : end if
227 : end do
228 :
229 0 : ff = (delta_total_J - sum_J_sync) / sum_J_non_sync
230 :
231 0 : do k=1,nz
232 0 : if (layers_in_sync(k) == 1) cycle
233 0 : if (delta_j(k) >= 0) then
234 0 : delta_j(k) = min(delta_j(k),delta_j(k)*pow2(tdyn_div_tkh(k))*ff)
235 : else
236 0 : delta_j(k) = max(delta_j(k),delta_j(k)*pow2(tdyn_div_tkh(k))*ff)
237 : end if
238 0 : if (delta_j(k) == s% j_rot(k) - j_sync(k)) layers_in_sync(k) = 1
239 : end do
240 0 : if (num_sync_layers == sum(layers_in_sync)) exit
241 0 : num_sync_layers = sum(layers_in_sync)
242 :
243 : end do
244 : !write(*,*) "tide iteration", num_sync_layers, &
245 : ! delta_total_J, sum_J_sync, sum_J_non_sync, ff
246 : else
247 0 : write(*,*) "sync_mode = " , sync_mode , " not recognized"
248 0 : write(*,*) "not doing tides"
249 0 : do k=1,nz
250 0 : delta_j(k) = 0d0
251 : end do
252 : end if
253 :
254 : end if
255 :
256 0 : if (b% point_mass_i /= 1 .and. b% s1% id == s% id) then
257 0 : b% t_sync_1 = t_sync
258 0 : if (b% model_twins_flag) then
259 0 : b% t_sync_2 = t_sync
260 : end if
261 : else
262 0 : b% t_sync_2 = t_sync
263 : end if
264 :
265 0 : if (.not. b% doing_first_model_of_run) then
266 0 : do k=1,nz
267 0 : s% extra_jdot(k) = s% extra_jdot(k) - delta_j(k)/dt_next
268 : end do
269 : end if
270 :
271 : end subroutine sync_spin_to_orbit
272 :
273 :
274 0 : real(dp) function f2(e)
275 : real(dp), intent(in) :: e
276 :
277 0 : f2 = 1d0
278 :
279 : ! Hut 1981, A&A, 99, 126, definition of f2 after eq. 11
280 0 : if (e > 0d0) then
281 0 : f2 = 1d0 + 15d0/2d0*pow2(e) + 45d0/8d0*pow4(e) + 5d0/16d0*pow6(e)
282 : end if
283 :
284 0 : end function f2
285 :
286 0 : real(dp) function f3(e)
287 : real(dp), intent(in) :: e
288 :
289 0 : f3 = 1d0
290 :
291 : ! Hut 1981, A&A, 99, 126, definition of f3 after eq. 11
292 0 : if (e > 0d0) then
293 0 : f3 = 1d0 + 15d0/4d0*pow2(e) + 15d0/8d0*pow4(e) + 5d0/64d0*pow6(e)
294 : end if
295 :
296 0 : end function f3
297 :
298 :
299 0 : real(dp) function f4(e)
300 : real(dp), intent(in) :: e
301 :
302 0 : f4 = 1d0
303 :
304 : ! Hut 1981, A&A, 99, 126, definition of f4 after eq. 11
305 0 : if (e > 0d0) then
306 0 : f4 = 1d0 + 3d0/2d0*pow2(e) + 1d0/8d0*pow4(e)
307 : end if
308 :
309 0 : end function f4
310 :
311 :
312 0 : real(dp) function f5(e)
313 : real(dp), intent(in) :: e
314 :
315 0 : f5 = 1d0
316 :
317 : ! Hut 1981, A&A, 99, 126, definition of f5 after eq. 11
318 0 : if (e > 0d0) then
319 0 : f5 = 1d0 + 3d0*pow2(e) + 3d0/8d0*pow4(e)
320 : end if
321 :
322 0 : end function f5
323 :
324 0 : subroutine get_tsync(id, sync_type, Ftid, qratio, m, r_phot, osep, t_sync, ierr)
325 : integer, intent(in) :: id
326 : character (len=strlen), intent(in) :: sync_type ! synchronization timescale
327 : real(dp), intent(in) :: Ftid ! efficiency of tidal synchronization. (time scale / Ftid ).
328 : real(dp), intent(in) :: qratio ! mass_other_star/mass_this_star
329 : real(dp), intent(in) :: m
330 : real(dp), intent(in) :: r_phot
331 : real(dp), intent(in) :: osep ! orbital separation (cm)
332 : real(dp), intent(out) :: t_sync
333 : integer, intent(out) :: ierr
334 0 : real(dp) :: rGyr_squared, moment_of_inertia, kdivt
335 : type (binary_info), pointer :: b
336 : type (star_info), pointer :: s
337 : integer :: k
338 :
339 : include 'formats'
340 :
341 : ierr = 0
342 :
343 0 : call star_ptr(id, s, ierr)
344 0 : if (ierr /= 0) then
345 0 : write(*,*) 'failed in star_ptr'
346 0 : return
347 : end if
348 :
349 0 : call binary_ptr(s% binary_id, b, ierr)
350 0 : if (ierr /= 0) then
351 0 : write(*,*) 'failed in binary_ptr'
352 0 : return
353 : end if
354 : ! calculate the gyration radius squared
355 0 : moment_of_inertia = 0d0
356 0 : do k=1, s% nz
357 0 : moment_of_inertia = moment_of_inertia + s% i_rot(k)% val*s% dm_bar(k)
358 : end do
359 0 : rGyr_squared = (moment_of_inertia/(m*r_phot*r_phot))
360 0 : if (sync_type == "Hut_conv") then
361 : ! eq. (11) of Hut, P. 1981, A&A, 99, 126
362 :
363 0 : kdivt = k_div_T(b, s, .true., ierr)
364 0 : if(ierr/=0) return
365 :
366 0 : t_sync = 3d0*kdivt*(qratio*qratio/rGyr_squared)*pow6(r_phot/osep)
367 : ! invert it.
368 0 : t_sync = 1d0/t_sync
369 0 : else if (sync_type == "Hut_rad") then
370 : ! eq. (11) of Hut, P. 1981, A&A, 99, 126
371 :
372 0 : kdivt = k_div_T(b, s, .false., ierr)
373 0 : if(ierr/=0) return
374 :
375 0 : t_sync = 3d0*kdivt*(qratio*qratio/rGyr_squared)*pow6(r_phot/osep)
376 : ! invert it.
377 0 : t_sync = 1d0/t_sync
378 0 : else if (sync_type == "Orb_period") then ! sync on timescale of orbital period
379 0 : t_sync = b% period ! synchronize on timescale of orbital period
380 : else
381 0 : ierr = -1
382 0 : write(*,*) 'unrecognized sync_type', sync_type
383 0 : return
384 : end if
385 0 : t_sync = t_sync / Ftid
386 : end subroutine get_tsync
387 :
388 0 : real(dp) function k_div_T(b, s, has_convective_envelope, ierr)
389 : type(binary_info), pointer :: b
390 : type(star_info), pointer :: s
391 : logical, intent(in) :: has_convective_envelope
392 : integer, intent(out) :: ierr
393 :
394 : integer :: k
395 0 : real(dp) :: osep, qratio, m, r_phot,porb, m_env, r_env, tau_conv, P_tid, f_conv
396 0 : real(dp) :: e2
397 :
398 0 : ierr = 0
399 :
400 : ! k/T computed as in Hurley, J., Tout, C., Pols, O. 2002, MNRAS, 329, 897
401 : ! Kudos to Francesca Valsecchi for help implementing and testing this
402 :
403 0 : k_div_T = 0d0
404 :
405 0 : osep = b% separation
406 0 : qratio = b% m(b% a_i) / b% m(b% d_i)
407 0 : if (is_donor(b, s)) then
408 0 : m = b% m(b% d_i)
409 0 : r_phot = b% r(b% d_i)
410 : else
411 0 : qratio = 1d0/qratio
412 0 : m = b% m(b% a_i)
413 0 : r_phot = b% r(b% a_i)
414 : end if
415 0 : porb = b% period
416 :
417 0 : if (has_convective_envelope) then
418 0 : m_env = s% m(1) / Msun ! assume fully convective
419 0 : r_env = s% r(1) / Rsun
420 0 : do k=1, s% nz ! search if _not_ fully convective
421 0 : if (s% mixing_type(k) /= convective_mixing .and. &
422 0 : s% rho(k) > 1d5*s% rho(1)) then
423 0 : r_env = (r_phot - s% r(k))/Rsun
424 0 : m_env = (s% m(1) - s% m(k))/Msun
425 0 : exit
426 : end if
427 : end do
428 : tau_conv = 0.431d0*pow(m_env*r_env* &
429 0 : (r_phot/Rsun-r_env/2d0)/3d0/s% L_phot,one_third) * secyer
430 0 : if (1d0/porb /= s% omega_avg_surf/(2d0*pi)) then
431 0 : P_tid = 1d0/abs(1d0/porb-s% omega_avg_surf/(2d0*pi))
432 0 : f_conv = min(1.0d0, pow(P_tid/(2d0*tau_conv),b% tidal_reduction))
433 : else
434 : f_conv = 1d0
435 : end if
436 :
437 0 : k_div_T = 2d0/21d0*f_conv/tau_conv*m_env/(m/Msun)
438 : else
439 : !NOTE:There is a typo in eq. (42) of Hurley+ 2002,
440 : !correct expression is given in footnote 3 of
441 : !Sepinsky+ 2007
442 0 : k_div_T = 1.9782d4*sqrt(m*r_phot*r_phot/pow5(osep)/(Msun/pow3(Rsun)))
443 0 : k_div_T = k_div_T*pow(1d0+qratio,5d0/6d0)
444 :
445 0 : if (b% use_other_e2) then
446 0 : call b% other_e2(s% id, e2, ierr)
447 0 : if(ierr/=0) then
448 0 : write(*,*) "Error when computing other_e2 ",ierr
449 0 : return
450 : end if
451 : else
452 : ! E2 from Hurley 2002 eq 43 based on Zahn 1975
453 0 : e2 = 1.592d-9*pow(m/Msun,2.84d0)
454 : end if
455 0 : k_div_T = k_div_T*e2/secyer
456 : end if
457 :
458 : end function k_div_T
459 :
460 : end module binary_tides
461 :
|
