[mesa-users] Temperature gradient from Canuto & Mazzitelli (1991)

[Warrick Ball]

Hi MESAlians,

What started off as a side interest following a conversation over a 
conference lunch 12 days ago has ended up with me possibly having 
implemented the temperature gradient from Canuto & Mazzitelli (1991) [1], 
or just CM1991.

Much of the original paper is committed to determining the convective flux 
according to a more detailed description of the spectrum of turbulent 
eddies than in mixing length theory.  These details are still a bit beyond 
me, but the results of those calculation are tabulated in the paper and an 
analytic fit given in their equation (32).  This introduces a new factor 
(eqn 34), involving the MLT convective flux (eqn 4), that must be 
incorporated in the traditional cubic equation of Cox & Giuli (eqn 14.82 
in the textbook; eqns 9 and 10 in CM1991).

The new form of the "cubic" (eqn 36) isn't really a cubic, so I used the 
safe_root routines from mesa/num to solve the equation, using the MLT 
value of the parameter Gamma as an initial guess.  I solved for the 
gradient of Gamma implicitly.  The parameter Zeta in the code becomes

Zeta = Phi_CM/(1 + Phi_CM)

(see eqn 63 in the paper) where Phi_CM is the new convective flux, given 
by the analytic fit (eqn 32).  Again, the derivative can be computed 
implicitly.

After some fumbling around, I seem to have succeeded in implementing the 
new temperature gradient in MESA.  I have attached a run_star_extras.f 
that contains a modified version of MESA's own MLT routine, in which the 
CM1991 temperature gradient is developed as a modification to the 'Cox' 
option.  Just set MLT_option = 'CM1991' and use_other_mlt = .true. in 
&controls.  (FYI, this was all developed in revision 7184.)

Attached is an HR track from a default work folder, modified to a 1 solar 
mass model.  The two tracks are for Cox MLT with alpha = 2 and CM1991 
MLT with alpha = 0.75, up to core hydrogen exhaustion.  They don't quite 
match, which isn't totally surprising because the conversion between the 
two is non-trivial.  I'm trying to calibrate a solar model now to see what 
value comes out.

Anyway, that's what I've done.  It started out as an exercise in teaching 
myself a bit more about convection in the models, and ended up with 
something I hope that someone, somewhere might find interesting or useful. 
I'd love to know if I've actually done the right thing.  I've iterated on 
the "cubic" equation (32), although the paper talks about iterating 
directly on the convective flux (equation 63).  I assume that these are 
equivalent...

The fact that I got somewhere with the code is testament to it's clarity: 
thanks again Bill and other contributors!  Finally, note that I haven't 
done anything regarding (a) the turbulent pressure, for which CM1991 also 
offers a different expression (I think?), or (b) the suggestion of 
limiting the mixing-length to the distance to the convective boundary.

Cheers,
Warrick

[1] http://esoads.eso.org/abs/1991ApJ...370..295C

------------
Warrick Ball
Postdoc, Institut für Astrophysik Göttingen
wball at astro.physik.uni-goettingen.de
+49 (0) 551 39 5069
-------------- next part --------------
A non-text attachment was scrubbed...
Name: run_star_extras.f
Type: text/x-fortran
Size: 68073 bytes
Desc: 
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141217/93e4d986/attachment.bin>
-------------- next part --------------
! inlist to evolve a 15 solar mass star

! For the sake of future readers of this file (yourself included),
! ONLY include the controls you are actually using.  DO NOT include
! all of the other controls that simply have their default values.

&star_job

  ! begin with a pre-main sequence model
    create_pre_main_sequence_model = .true.

  ! save a model at the end of the run
    save_model_when_terminate = .false.
    save_model_filename = '15M_at_TAMS.mod'

  ! display on-screen plots
    pgstar_flag = .false.

/ !end of star_job namelist


&controls
    use_other_mlt = .true.

    MLT_option = 'Cox' ! 'CM1991'
    mixing_length_alpha = 2 ! 0.75

  ! starting specifications
    initial_mass = 1 ! 15 ! in Msun units

  ! stop when the star nears ZAMS (Lnuc/L > 0.99)
    Lnuc_div_L_zams_limit = 0.99d0
    stop_near_zams = .false. ! .true.

  ! stop when the center mass fraction of h1 drops below this limit
    xa_central_lower_limit_species(1) = 'h1'
    xa_central_lower_limit(1) = 1d-3

/ ! end of controls namelist
-------------- next part --------------
A non-text attachment was scrubbed...
Name: HR_CM91.png
Type: image/png
Size: 73284 bytes
Desc: 
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141217/93e4d986/attachment.png>