Meeting Update --12/13/16
- Wire Turbulence.
- Merging Aluminum Cooling Table in the code
Table
Temperature(T) range 1-100K dT 0.1K Density range in the table - 1/cm3
Current Isothermal run without cooling
Wire Temperature 3.75E-12 K Wire Density 4.8E26 1/cm3 Wind Temperature 1.5E-8 K Wind Density 4.8E23 1/cm3
- 2D
- Bruce's OH231 runs
RGB star with sink particle
1) Understand how MESA stellar profiles translate to the grid.
It is important to make sure that the input from MESA is being gridded correctly, so compare with Ohlmann+16c https://arxiv.org/abs/1612.00008v1. See last pages of
this updated presentation. We see that the core gets cut out due to lack of resolution at the center, as expected. Though the units used are different from the Ohlman plots, I tried to lign them up. It can be seen that the initial pressure, density and specific internal energy profiles are in agreement with Ohlmann, as expected. This confirms the following:
(i) (Slide 5 & 6): The MESA density and pressure profiles that we put in are basically what we get out, except for the very center, which gets washed out due to lack of resolution (the resolution is slightly larger than 1 R_sun, which is where the profiles become flat). These profiles also agree with those of Ohlmann;
(ii) (Slide 7): A cutoff radius of 1 R_sun corresponds to a missing central core mass of about 0.4 Msun (which is the core mass used by Ohlmann);
(iii) (Slide 8): The profile for the specific internal energy looks very similar to Ohlmann's profile (though the units are different). This confirms that our code is using the same (ideal gas) equation of state as Ohlmann.
Now that this sanity check has been done, we can try introducing a sink particle at the center.
2) Introduce sink particle to replace the core. Here I've tried sink particles of the following masses:
a) 0 Msun,
b) 0.2 Msun,
c) 0.4 Msun (closest to the actual core mass),
d) 0.8 Msun,
e) 1.7 Msun,
f) 3.4 Msun,
g) 10 Msun.
Movies of density (both 2D slices and 1D profiles along x) are here. Simulations have resolution 1283 in a box length 140 R_sun.
Note that the star is not stable for any of the runs, but the 0.4 Msun sink particle does at least seem to help to stabilize the star compared with the case of no sink particle (see remarks below).
For the 0 Msun case (no sink particle) I've checked that the initial local sound-crossing time t_s(r) = r/c_s(r) is comparable to the time it takes for the outgoing pressure disturbance to propagate to radius r, for a few different values of r. Actually, t_s(r) is a few times larger than this expansion time, but within a factor of order unity, as expected. A plot of the initial sound-crossing time for the 0 Msun case is available here and can be compared with the expansion times in the movies showing density mentioned above. Movies of the 1D profile of the sound speed for the 0 Msun and 0.4 Msun cases are available here.
Another useful quantity is the fractional difference from hydrostatic equilibrium:
hef = | |grad P| - |rho g| |/max(|grad P|,|rho g|).
This has been plotted in the following fractional HE movies.
Remarks:
1) A sink particle with mass equal to the core mass helps to preserve the high density at the center, and slows the expansion due to hydrostatic imbalance, but also creates outgoing sound waves eminating from the center;
2) The presence of the sink particle somehow introduces unphysical 4-fold symmetry, followed by complete asymmetry. This seems to be a numerical problem that needs to be addressed;
3) Though it is not computationally optimal for this problem, it would make life simpler to change back to cgs units as the computational units here (rather than length units of R_sun, as in the current simulations). This would allow easier comparison with Ohlmann results and prevent unit conversion errors for these trial runs. I will implement this change in the next set of simulations.
Next step is to solve Ohlmann's modified Lane-Emden equation inside the truncation (cutoff) radius. The choice of truncation radius determines the core mass = m(r_truncation). Then match this solution onto the MESA profile at the truncation radius, as done by Ohlmann (density should be smooth across the transition at r_truncation). I will write a code to do this. Hopefully this will result in a more stable star.
Update 12/13
MHD
Anisotropic MHD simulation at a surface beta of 50 seems to have gone pretty well:
And top view:
The magnetic field goes a bit crazy at the end, but the flow still looks about what we'd expect, with outflow from the night side highly suppressed and material there concentrated near the equator by the field. Best, we don't get the strange fourfold outflow and fallback. Seems narrower than the suppression seen by Owen & Adams, but that may be because of a difference in temperature profiles (they use radiative transfer).
I'm currently running a simulation with a beta of 10, which may lead to more coherent field lines - I can next run a simulation with similar field strengths to Owen & Adams, if the anisotropic simulations appear to be good.
WASP
The planet has only existed for half an orbit, so there's not much to see for troubleshooting purposes yet. But so far it looks as expected, with lobes spiralling up-orbit and down-orbit. Side view shows strange wave patterns starting in stellar wind around frame 7 - perhaps I need to increase the size in z a bit more?
Putting RGB star on the grid: update
1) Movies —> all units are cgs
Simulation 1: Full star on the grid, resolution 2563
a) density, variable color bar b) density, fixed color bar c) density 1D profile d) density 1D profile, zoomed
Simulation 2: Central part of star on the grid, resolution 643 with 3 levels of amr, so up to 2563 (zoomed in by factor 50 from above movies, so most of star is "outside" the grid)
a) without mesh b) with mesh c) inner region close-up, with mesh
2) Original profiles from MESA slides
Notes on slides:
a) Slide 4 shows the scale height as a function of radius (upper plot) and pressure v. radius (lower plot). In the lower plot, the RGB star is orange. The scale height dips significantly below the local radius below r= 2e-2 Rsun. b) Slide 5 shows the density as a function of radius. The values can be compared with what is achieved in the simulations, e.g. at the center of the star.
3) Discussion
At yesterday's meeting it was decided that resolving the center is not possible, and that the most promising approach is to introduce a sink particle at the center and then modify the pressure profile accordingly to account for the modified potential, as done in Ohlmann 2016 (PhD thesis). Serendipitously, the paper on exactly this subject came out on arXiv today: https://arxiv.org/pdf/1612.00008.pdf