COMMON ENVELOPE SIMULATIONS
New Work
- Update on Yisheng's work
- Update on Paper
Update on Yisheng's work
- Detailed study of initial conditions. Here are the notes.
- Improving graphs
- Sampling rate for inter-particle separation
- Mark locations of periastrons and apastrons on accretion vs time graphs
Update on Paper
Here is the pdf so far.
- Explanation for inter-particle separation graph.
- Other comments I inserted in the text (about issues to address).
- Decide the way forward.
- Next paper?
Update 3/26
Submitted WASP-12b revision.
Parameter space paper coming along. What still needs work:
- Intro: Have a sketch, but still needs citations and probably needs to be filled out more.
- Methods: Regime (Energy vs cooling limited)
- Discussion: Have an outline. Anything else that stands out from the figures?
Run6 for paper stalled, Run5 almost done.
Radiation Pressure
The higher-flux run looks good. Switched from plotting X to plotting 1-X, so log plots now work well (also changed this in the figures for the paper).
2d14, top
2d14, side
2d14, force comparisons
2d15, top
2d15, side
2d15, force comparisons
Update 3/26: Wind Tunnel Run 007
Finished Runs
WT Run 007
Same setup with Run006B: softening length resolved by 16 pixels. Region of highest refinement et to 30 CUs. Box size is set to be 20483 CUs, 23 as large as those of Run 006s, 43 as large as Run 005s. This enables a increase of run time to 2CU.
Log density & velocity vectors movie:
As shown in the movie, although the density disturbance around the particle is almost always contained within the box, near the end of the simulation gas velocities at the boundaries are altered due to gas self-gravity. In other words the gas across the box begins to collapse onto itself. This is to be expected since the box is larger than the Jeans scale of the gas, and self-gravity introduces positive feedback to "accretion" around the point particle.
The dynamic friction also confirms this positive feedback. In the shorter timescale used before (¼ of total runtime), Run 007 does not stand out -
But the drag blows up later -
Which can be fitted by an exponential function - F = f0*Exp(t/tau) + c
Here f0 = 0.157, tau = 7.442, and c = 2.324.
Constant offset c here can represent the drag force as predicted by BHL, ignoring self-gravity.
To avoid using incorrect DF measurements due to limited box size, I am running currently Run 007B, which is the same setup as Run 007 but with a factor of 23 box.
Next Steps - parameter space? Non-dimensionalization?
I think the simulation module has become mature enough that exploring a larger portion of the parameter space becomes possible. MacLeod suggested me to non-dimensionalize the problem, to better aim at the physics. I am studying his (their) approach - there might be some difficulties in, for instance, the scalability of density. But even in cgs we can still artificially alter Mach number, density and the particle mass, etc.
Radiation Pressure Testing
Comparisons of force from radiation pressure and magnitude of gradient of thermal pressure, and magnitude of gradient of ram pressure (not actually a force, doesn't seem useful for comparison):
Comparisons of ambient pressures in frame 150 and frame 151. Here is where the problem lies - ambient pressure is increasing due to OUTFLOW_WIND boundary conditions being set in a different way than before. This disrupts the steady state and masks the effect of radiation pressure on the stellar side of the wind.
Neutral fraction, temperature, density, and
, which is dependent on these. At edge of wind, between ~4 and 5, the temperature drops and the neutral fraction increases enough that most of the Lyman-alpha radiation is absorbed near the edge.COMMON ENVELOPE SIMULATIONS
New Work
- Corrected the main figures from last post
Plots relating to the force
- The quantity being plotted is
where (repeated from last post):
- and are the potentials due to the primary (RG core) and secondary (companion), respectively. Actually, inside the softening radius, we have also corrected for the spline potential by including extra terms, not written out here;
- is the potential of the gas;
- is the centrifugal potential in the frame corotating with the particles' orbit. Here is the angular velocity of this frame with respect to the lab frame and is the distance from the secondary in the x-y plane;
- is the potential due to the shift of the rotation axis of the rotating frame from the CM to the secondary;
- Last post, there was an error in which had the incorrect sign.
- is the Coriolis term;
- is the angular velocity of gas about the secondary in the frame corotating about the secondary with the instantaneous orbital angular velocity of the particles, and is the centrifugal force due to the motion of the gas IN THE COROTATING FRAME (i.e. we have already accounted for the rotation of the reference frame, but here we account for the rotation of the gas within the rotating reference frame);
- We normalize with respect to , that is, the magnitude of g due to the secondary alone at the location of the primary.
- Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right.
- Comments (itmes repeated from last post):
- Vectors have been ommited from inside the softening circle of the secondary for clarity as the magnitudes were large in some cases.
- In both cases, some of the gas around the secondary is accelerating inward while some is accelerating outward.
- So it is probably incorrect to conclude that gas is bound to the secondary.
Plots relating to the energy
- The quantity being plotted in pseudocolor is the Bernouilli parameter , where the first term is the specific kinetic energy, the second term is the specific enthalpy ( ), and the last term is the external potential.
- Since and , we have . Here is just the total energy density minus the bulk KE density of the gas.
- For plots in the frame rotating around the secondary, we set (so exclude the gas).
- For plots in the frame rotating around the particle center of mass, we set ,
where
.- Contours show . Contour levels are the same in each panel.
- Both pseudocolor and contours are normalized by .
In the frame rotating around the particle CM:
In the frame rotating around the secondary:
Update 3/15
Exoplanets II abstract:
Using AstroBEAR, a 3D MHD fluid dynamics simulation code, we have developed a new platform on which to study the dynamics of planetary winds, including a variety of microphysics such as ionization heating, radiation pressure, and charge exchange. We will present the results of a study of the effects of planet mass and stellar flux on the structure of hot Jupiter winds, as well as changes to this structure that result from introducing radiation pressure.
WASP-12b paper sent out to authors for review of referee response.
Parameter space paper coming along. First cut of methods done, with a start on results. How much repetition is appropriate for figure captions/commentary? Intro is next on list of things to write.
Related, running out of space on BlueHive again. My personal quotas are full on BlueStreak (could request more - only have 200 GB here), BlueHive, and Stampede, and afrank_lab is full on BlueHive. We could back up and remove all but the final frames of the four low-flux cases (do we need the low-flux, high-density, low-pressure test for anything?).
On the simulations, the low-pressure ambient test for Run5 is proceeding nicely. The low-pressure ambient tests are running significantly faster than the higher-pressure ones - something to keep in mind when designing future planet simulations.
LyAlphaFlux = 0
This seems to indicate a problem. Running a couple of tests indicates that there's no radiation pressure in the first dozen or so steps, and the flux is zero throughout the grid, so it may just be that we haven't reached a steady state yet in Run2 (on further consideration, though, this doesn't seem likely to be the case - Run2 looks very stable between frames ~120 and 150, at least).
Two possible tests spring to mind - let test run go further with debug output in the line transfer routines, or run Run2 further with iLineTransfer=1, see if the same thing occurs. Some minor parts (that shouldn't matter on a restart) were changed between the original run and the radiation pressure run, primarily the number of zones between the outer boundary and the zero radius.
Radiation Pressure, Run2, top
Radiation Pressure, Run2, side
Planet Updates 3/13
Planet density comparisons
Go here for gif/stills; test2_side and test2_top are the final frames acquired (frame 97)(the GIFs don't make it there).
2c max ram pressure @ 1.5 Rp, frame 65 = 1.74x10-4
2c max ram pressure @ 1.5 Rp, frame 20 = 1.04x10-4
2c initial ambient pressure = 3.9x10-8
2 max ram pressure @ 1.5 Rp, tfinal = 8.375x10-4
2 max ram pressure @ 1.5 Rp, frame 98 = 7.77x10-4
2 max ram pressure @ 1.5 Rp, frame 65 = 8.47x10-5
2 max ram pressure @ 1.5 Rp, frame 20 = 2.7x10-5
2 initial ambient pressure = 1.7x10-8
Note that the original Run2 didn't have any buffer zones (pressure matching was essentially right at the outer boundary), which is why the ambient pressure is lower.
COMMON ENVELOPE SIMULATIONS
New Work
- Figure of radial force/unit mass (=instantaneous acceleration) pseudocolor with force/unit mass vectors overplotted.
- Figure of Bernouilli parameter pseudocolor with contours of potential overplotted.
Plots relating to the force
- The quantity being plotted is
where
- and are the potentials due to the primary (RG core) and secondary (companion), respectively. Actually, inside the softening radius, we have also corrected for the spline potential by including extra terms, not written out here;
- is the potential of the gas;
- is the centrifugal potential in the frame corotating with the particles' orbit. Here is the angular velocity of this frame with respect to the lab frame and is the distance from the secondary in the x-y plane;
- is the potential due to the shift of the rotation axis of the rotating frame from the CM to the secondary;
- is the Coriolis term;
- is the angular velocity of gas about the secondary in the frame corotating about the secondary with the instantaneous orbital angular velocity of the particles, and is the centrifugal force due to the motion of the gas IN THE COROTATING FRAME (i.e. we have already accounted for the rotation of the reference frame, but here we account for the rotation of the gas within the rotating reference frame);
- We normalize with respect to , that is, the magnitude of g due to the secondary alone at the location of the primary.
- Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right.
- Comments:
- Vectors have been ommited from inside the softening circle of the secondary for clarity as the magnitudes were large in some cases.
- It would be much easier to simply plot the instantaneous acceleration instead of computing it, but we do not have this quantity available in the chombo.
- In both cases, some of the gas around the secondary is accelerating inward while some is accelerating outward.
- So it is probably incorrect to conclude that gas is bound to the secondary.
- What we need are streamlines (or actually, pathlines) https://en.wikipedia.org/wiki/Streamlines,_streaklines,_and_pathlines to really get a sense for what kinds of orbits the fluid is making in the vicinity of the secondary.
- Below are the plots of the various contributions to the force per unit mass, first for run 143 (no sub-grid accretion). In order from left to right, we have
- Pressure term ;
- Gravity terms ;
- Gravity terms (excluding gas);
- Gravity terms (lab frame);
- Centrifugal term .
- Coriolis term ;
- Now for Run 132:
Plots relating to the energy
- The quantity being plotted in pseudocolor is the Bernouilli parameter , where the first term is the specific kinetic energy, the second term is the specific enthalpy ( ), and the last term is the external potential.
- Since and , we have . Here is just the total energy density minus the bulk KE density of the gas.
- We set (so exclude the gas).
- Contours show .
- Both pseudocolor and contours are normalized by .
- Comments:
- Contours show the values -2.5, -2.25, -2, -1.75, -1.5, -1.25. The contour -1 is plotted but does not appear, which means that the values are < -1 everywhere on the plot.
- The divide gives a sense for which gas is bound to the particles (red) or unbound (blue).
- The gas in the vicinity of the particles is mostly bound.
- To get a sense of whether the gas is bound to the secondary, in particular, we can refer to the Roche potential contours.
- Note that alternatively we could have omitted from the calculation, to try to get a sense of whether material is bound to the secondary (rather than to both the secondary and primary together). However, this would not really be correct because it would imply that any "red" material near the primary was bound to the secondary when in reality it would be more tightly bound to the primary than to the secondary.
Discussion
- I think it might be worth including three plots in the paper (for each of the two runs). All are in the frame corotating with the particles' orbit and rotating around the secondary:
- Tangential velocity (blue/red pseudocolor), normalized with Keplerian speed for a point potential, with velocity vectors overplotted (see Feb 15 post).
- Force/unit mass (ie. instantaneous acceleration) along radial direction (blue/red pseudocolor, normalized) with acceleration vectors overplotted (this post).
- Bernouilli parameter (blue/red pseudocolor, normalized) with external potential overplotted (this post).
- We should also have the pathlines, ideally, but if this is not possible, then streamlines, to get a better sense of the trajectory of the gas near the secondary.
Next steps
- Finalize the remaining plots and insert into paper.
- ADAF papers; read and discuss.
- Text of paper.
- Explore energy conservation in the simulations.
- Is there enough power to drive a jet through the envelope? Estimates.
COMMON ENVELOPE SIMULATIONS
New Work
- Generated skeleton of paper including most of the "polished" figures and put on sharelatex.
- "Settled" the "enthalpy vs. thermal energy" issue and also the "gas bound to itself" issue with Eric's help.
- Wrote up method to calculate velocities in the corotating frame.
- Made polished plots of edge-on (to the orbital plane) views of gas density at two different zoom levels and inserted into paper.
- Made zoomed in plots of face-on gas density similar to the edge-on ones.
- Plots of the potential in the frame corotating with the particle orbit, centered on the companion, with and without the gas potential.
- Plots of , which relates to the force per unit mass along the radial direction with respect to the companion.
- Plot of square of tangential velocity with respect to the companion, in the corotating frame, but normalized with respect to .
Skeleton of paper
Here is the pdf so far.
"Enthalpy vs. thermal energy" and "gas bound to itself" issues
- Enthalpy vs thermal energy
- When assessing whether the gas is bound we evaluate whether a given quantity relating to the gas energy density is >0 (unbound) or <0 (bound).
- The question arose whether it should be the internal energy density that enters the equation (more specifically, the internal translational kinetic energy density of the gas particles related to the gas kinetic temperature), or the enthalpy per unit volume (the work per unit volume required to place a fluid element at that position within its environment, equal to the thermal energy density + the pressure).
- Since the pressure can also contribute to the outward motion of the gas (work against gravity), it makes sense to include the enthalpy per unit volume in the calculation, not just the thermal energy per unit volume.
- Gas bound to itself
- Similarly, the question arose whether to include the negative contribution of potential energy arising from gas self-gravity.
- On one hand, one would think yes because this gravity contributes to the binding of the gas-particle system.
- However, some of the gas may be unbound and escape, and the binding energy associated with this unbound gas should not contribute. For example, for a gas parcel that is bound to itself but not bound to the rest of the system, the potential energy associated with its self-gravity should not be included.
- To resolve this conundrum it helps to be precise about what we really want to call "bound". What we are probably most interested in is whether gas is bound to the particles (the RG core and companion). We are not as interested in whether the gas is bound to itself. Therefore, it is best NOT to include the potential energy due to the gas self-gravity.
The conclusion then is that 1) we should consider the enthalpy per unit volume rather than just internal translational kinetic energy per unit volume, and 2) we should consider the potential energy associated with the gas-particle interaction only, not including the gas-gas interaction. Both of these choices work in the direction of making the gas less "bound".
Method to calculate velocities in corotating frame
Here are some notes.
Edge-on plots of density
- Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right. Units are grams/cm3. Color bar is the same for both plots.
- Zoomed-in with different color scheme. Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right. Units are grams/cm3. Color bar is the same for both plots.
New face-on plots of density
- As above but now face-on. Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right. Units are grams/cm3. Color bar is the same for both plots.
New face-on plots of potential
- Total potential in frame rotating with orbital angular velocity and centered on the companion in cgs units. Includes gas potential. Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right.
- Same as above but NOT including gas potential in cgs units. Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right.
- Ratio of gas potential to total potential not including gas potential in cgs units. Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right.
- Comments: the gas potential contributes at about the 20% level at most.
Plots relating to the force
- The quantity in cgs units. Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right.
- Comments: Here includes the particle and gas gravitational potentials in the frame corotating with the orbit and centered on the companion. The force therefore includes the centrifugal force due to the rotating frame. HOWEVER, it does not include the Coriolis force/unit mass .
- Positive values indicate that the radial force with respect to the companion is directed outward (in the corotating frame), while negative values indicate that it is directed inward (in the corotating frame).
Plot of tangential velocity squared (with respect to the companion) normalized against the quantity plotted above (with a minus sign inserted)
- The quantity in cgs units. Run 143 (no sub-grid accretion) on the left and Run 132 (Krumholz sub-grid accretion) on the right.
- Comments: Values are plotted from 0 to 1. Below 0, the radial force is outward (white). The colored region is for inward radial force. (But note that the colored region does have parts that are almost white, which is a drawback of the color scheme used).
- For circular motion we would expect this quantity to be equal to unity because the centrifugal force in the frame of a moving gas parcel would balance the inward radial force. Therefore < 1 suggests accretion while > 1 suggests outflow.
- However, we have not considered the radial or vertical gas motions, so it is not clear from the above quantity alone whether gas will end up moving toward or away from the companion.
- Moreover, we have neglected the Coriolis force , which should be included in the analysis. Positive (counter-clockwise) values of would produce an outward radial acceleration.
Next steps
- Finalize the remaining plots and insert into paper.
- ADAF papers; read and discuss.
- Text of paper.
- Explore energy conservation in the simulations.
- Is there enough power to drive a jet through the envelope? Estimates.
Update 3/6
Updated WASP-12b paper to address referee concerns. Should new figure be included in the paper, or just the referee response?
Updated summary post with run statuses.
Started parameter space paper.
Radiation pressure:
Meeting update -- 03/06/2018
- JetClump module with 3D MHD compiles and runs
- Tested the post-processing python code for polarization map on the 3D MHD JetClump results.
- Matlab code for plotting polarization map.
- Details can be found on this page
Update 3/5: Wind Tunnel Runs 006AB&D, as well as reviewing Run 005B
Work in the past two weeks has been focused on testing an AMR scheme that only forces maximal refinement within a radius of the particle. Currently the module is working fine, except that it cannot have the box be relaxed first before introducing the particle.
Finished Runs
WT Run 006A
Same setup with Run 004A. No accretion, highest resolution 0.5CU, softening length (spline) is resolved by 8 pixels. Radius of highest refinement is set to 40CU, which is more than twice as large as the Bondi radius 18.7CU. Significant run time improvement (order of 8) & low storage requirement (~300M per frame vs. ~4G per frame w/ normal AMR).
WT Run 006B
Similar setup with Run 006A, except resolution is improved by a factor of two to 0.25CU, and softening length is resolved by 16 pixels. To compensate, radius of highest refinement is set to 20CU, still larger than Bondi radius. Similar run time as the previous run.
Log density movie with mesh grids:
WT Run 006C
Similar to 006B, but with radius of refinement set at 40CU. The runtime becomes too long (longer than the normal AMR), and was left unfinished.
WT Run 006D
Similar to 006A, low resolution, but larger radius of highest refinement - 80CU. The run time is less than 2 days, similar to the AMR case. But the result looks unreliable.
Log density movie with mesh grids and velocity vectors:
Comparison with Run 004A
More one Run 005B
Zoomed-in log density movie: Zoomed-in Mach number and velocity vector movie;
CE to PN parameters
There are two simulations distinguished by the momentum in flux: higher momentum flux and lower momentum flux.
: density of the outflow
: velocity of the outflow, spherical
ID | higher | lower |
300 | 300 | |
quiet time(day) | 0-500 | 0-6000 |
outburst time(day) | 500-1000 | 6000-indefinite |