[mesa-users] accretion

[Dean Townsley]

Hi Bill,

Sorry to drop out towards the end of the week...

A lot of stuff to reply to by this point  :)    But I think we are 
pretty much having the same thoughts.

I think concerns about advection and the cfl are good, but I'll say more 
about that in another reply.

On 10/04/2013 05:15 PM, Bill Paxton wrote:
> I'm extremely interested in understanding how the L > Cp*T*Mdot 
> constraint can be compatible with large dt*Mdot.

I haven't said a lot about the use of this constraint because I wanted 
to do some tests to try to understand it better.  Broxton and I have 
done some of those, so I think I have a better idea of how this 
constraint acts.

> My understanding is that we need to limit dt such that this inequality 
> holds in the outer dt*Mdot of mass.  So for time-constant Mdot and 
> space-constant L (near the surface), we need to make sure that for all 
> cells k where we want to use the Townsley&Bilsten formula for eps_grav 
> we have Cp(k)*T(k) < a value L/Mdot that is independent of k.  We 
> expect Cp(k)*T(k) to increase as we go inward and where it reaches 
> L/Mdot, we should switch to some other way to get eps_grav. You seem 
> to be saying that we can actually satisfy Cp(k)*T(k) < L/Mdot for very 
> thick layers at the surface, even reaching to the center in some 
> cases.  Excellent!  In that case, the timestep 
> limit L_div_CpTMdot_limit will not be setting the timesteps at all. 
>  Is that what you find when you run your tests with the previous 
> version of mesa?   Perhaps small Mdot allows nice large dt?

Yes.  We tried a case with Mdot=10^-11 Msun/yr and L_div_CpTMdot did not 
go below the default limit (=2) anywhere in the outer part of the star.  
(It goes to zero at the center because L goes to zero there.)

In terms of shape, L_div_CpTMdot falls pretty steeply in the outer 
layers, then has a minimum near the base of the accreted layer, goes 
back up for a bit and then goes to zero at the center.

For Mdot=10^-10, L_div_CpTMdot only goes down to about 1.  So the limit 
does do something there.



Two notes about the limiter that helped me understand it better:

It is phrased in a way that I had trouble with, but I think I've sorted 
out.  I think the thing we really want to be talking about is t_th/t_acc 
- this is what we want this limiter to keep small for the layer in which 
we use the quasi-steady approximation.  This is
   t_th/t_acc = Cp*T*dM/L / (dM/Mdot)  = 1/L_div_CpTMdot
where dM is basically a mass coordinate inward w.r.t. the surface. Note 
the default limit keeps L_div_CpTMdot > 2, corresponding to t_th/t_acc < 
0.5.  So I would say the default constraint is not very strict, which I 
think is neither good nor bad, it just is.  We usually set things like 
this as loosely as possible to still get reasonable results, and even 
this loose constraint, I think, does give "reasonable" results, even 
with the bad spot in the eps_grav profile shown in MESA paper 2.

I also want to just say how the L_div_CpTMdot limit is used in the code, 
as it was not obvious to me from the comments.  First how t_th/t_acc 
varies...  In most of the envelope L is constant.  Mdot is constant.  Cp 
is pretty much constant for an ideal gas, and T is increasing with 
depth.  Since the dM cancels we basically get
   t_th/t_acc \propto T
So it starts small for the outer part of the envelope and then 
increases.  The way the L_div_CpTMdot_limit works is that it sets a 
timstep constraint such that
   dt = dM( t_th/t_acc=limit) / Mdot
where dM( t_th/t_acc=limit) is the depth in the star at which t_th/t_acc 
is the set limit.

One issue with the way the limit works is that while all "newly 
accreted" material has t_th/t_acc < limit, the steady-state 
approximation is used up to dM=5*Mdot*dt by default.

We are currently in the process of trying out some other parameter 
combinations to see how they behave.  We'll report what we find.

Cheers,
Dean