Increase in MinDensity seems to be a fix & status of wire turbulence jobs
Here is the page. Increasing the minimum density helps with the issue I brought up in the previous post. See for yourself:
Here is a side-by-side of the same parameters (beta = 10, Ma = 3, Thickness=0.01) but with different min densities at frame 40. On the left you have a smaller minimum density of 1d-10, where on the left we have 1d-2.
Here are two side-by-side gifs. Sorry if it slows down the webpage.
MinDensity = 1d-10 | MinDensity = 1d-2 |
![]() | ![]() |
Now after doing the restarts for mindensity = 1d-3, the first one to complete with that value was the beta=10, ma=20, thickness=0.06 case.
Here is a GIF of the new frames with that mindensity. With the old mindensity of 1d-2, there would be this wave that would propagate out immediately from the restart.
Job Update
Hydro/MHD | Ma | Thickness | Frame (Restart) |
Hydro | 20 | 0.01 | 200, completed |
" | " | 0.03 | 200, completed |
" | " | 0.06 | 200, completed |
Beta=0.1 | 20 | 0.01 | 200, completed |
" | " | 0.03 | 45* |
Beta=10 | 20 | 0.01 | 200, completed |
" | " | 0.03 | 87* |
" | " | 0.06 | 121* |
" | 3 | 0.01 | 76 |
" | " | 0.03 | Restart from 0 |
" | " | 0.06 | Restart from 0 |
- Beta=0.1, Ma = 20, Thickness = 0.03 is queued up.
- Beta=10, Ma = 20, Thickness = 0.03 is queued up.
- Beta=10, Ma = 20, Thickness = 0.06 is queued up.
- Beta=10, Ma = 3, Thickness = 0.01 is running currently (7 hrs at this post). Expected to end in 11.9 hrs at frame 76.1 !
- Beta=10, Ma = 3, Thickness = 0.03 is queued up.
- Beta=10, Ma = 3, Thickness = 0.06 is queued up.
- All Ma =3 were done with 1d-2 to go faster. As you can see from previous blog post, there were some errors. So I decided to restart them with 1d-10 (original mindensity).
- * indicates those changed with mindensity of 1d-3 which I had 1d-2 a few weeks prior. However the visualization, as we saw, looks bad, so I need to restart them from last good frame.
- 5 jobs are completed.
Job Update (Week of 05-04-2015)
Hydro/MHD | Ma | Thickness | Frame (Restart) |
Hydro | 20 | 0.01 | 200, completed |
" | " | 0.03 | 200, completed |
" | " | 0.06 | 200, completed |
Beta=0.1 | 20 | 0.01 | 200, completed |
" | " | 0.03 | 50 |
Beta=10 | 20 | 0.01 | 200, completed |
" | " | 0.03 | 100 |
" | " | 0.06 | 135 |
" | 3 | 0.01 | 200, completed |
" | " | 0.03 | 49 |
" | " | 0.06 | 71 |
Meeting Update 04/27/2015 -- Baowei
- Planetary Wind
- wiki:u/bliu/PlanetaryWind
- Back flow for gamma=1.01, mach=0.8
- Stand off distance Vs. ns
- Eddie's runs
- status of runs
d2.5_M15 | 22/50 |
d4.5_M10 | 46/50 |
d4.5_M15 | 23/50 |
d6.5_M10 | 48/50 |
d6.5_M15 | 23/50 |
- transferring from Gordon to BH. ~5TB total.
- current space usage on local machines — blog:bliu04242015
- XSEDE Allocation
- Current SUs — wiki:ProjectRuns
- New SDSC Machine — wiki:Czarships/ExternalResources/Machines (~47000 computing cores, standard queue: 1700 cores for 48 hours, fast, currently free)
A reduced grid space for Planetary Wind Runs
The stand-off distance for the stellar wind/planetary wind R_so can be computed simply by balancing the ram pressure
Thus we have
Thus dropping the stellar wind density by 10 should increase the stand-off distance by ~ 3.
Also lets keep keep v_s the same and change the T_s to modify the Mach number
Here are the runs I think we should focus on
# | lambda | Solar Wind Mach # | |
1 | 1.01 | 0.8 | 10-5 |
2 | 1.01 | 5.0 | 10-5 |
3 | 1.3 | 0.8 | 10-5 |
4 | 1.3 | 5.0 | 10-5 |
Using AstroBEAR to exam the theory
I have calculated the density variation in equatorial plane of the giant in a binary system. The parameters are:
Temperature: 3000 K
Giant mass: 1.01 SM
Companion mass: 0.5 SM
Separation: 4 AU
This plot shows the density variation of 0.2 AU, 0.4 AU and 1 AU in equatorial plane of the giant compared to the highest density at that radius.
Absolute density profile at 1 AU.
I plotted the isodensity contour of equatorial plane. It is not a circle indeed. The 3d contour has a hole because of the singularity on z-axis.
If we use AstroBEAR to exam the theory, I suggest:
Simulation zone: 10 AU * 10 AU * 5 AU
Giant boundary at 0.2 AU
Base resolution at 100*100*50
AMR level at 5 or 6
the smallest cell will be 1/32 AU = 0.03125 AU or 1/64 AU = 0.015625 AU.
Updated Outflow Module
The outflow module has been updated with a few new features, so that I can use it for my 3D, pulsed, magnetized, cooling jets. I basically added 3 new profile objects to the outflow object:
- SpeedProfile: Before, you could only give the outflow object one outflow speed, and you could make it oscillate sinusoidally if needed. Now, the outflow object will first check to see if there is a SpeedProfile. This is a profile of the outflow speed as a function of time. I put this in here so that I can eventually set up a profile of random velocities to use.
- TempProfile: The original outflow object just takes one temperature for the entire outflow. Now, with the TempProfile, you can give it a profile of temperature with respect to radius. This was necessary to use the hydromagnetic pressure equilibrium profile that I have used in my jets.
- MagProfile: The original outflow object took one value for the peak toroidal magnetic field. Now, you can pass it a MagProfile which a profile of toroidal magnetic field strengths as a function of radius. The functionality for other field geometries is not there yet, but I, and many other projects, would only need toroidal fields anyways. We can implement helical fields later.
I built a new PulsedJet problem module which uses this new outflow object set up. In my module, I define an appropriate temperature and magnetic field profile. I have not done the random velocity profile yet, but that should be pretty straight forward. I plan on putting all of these profiles into profile_functions.f90 so that anyone can easily use them.
Below is a very low resolution example of a 3D, magnetized, pulsed, cooling jet. It is a 2D slice of density.
The bow shock seems pretty wide to me, and I have not yet figured out why. I have checked my temperature profile several times, and it is correct. Perhaps there is some confusion between thermal pressure and internal energy within the outflow object. I'm going to investigate this further just to make sure everything is working accurately. The magnetic field part looks like it's working perfectly, and I haven't seen any restarts yet.
This simulation is pretty slow since I'm just using 16 procs on alfalfa. I'll let it run further in order to see if there are any more bugs that need to be worked out.
UPDATE
I fixed a bunch of stuff related to the new outflow module:
- Got the velocity pulses working.
- Edited the new pulsed jet problem module to handle hydro (non-MHD).
- Edited profiles.f90 to make it more user-friendly.
- Fixed some issues with the magnetic field part of the outflow object. This involved a time dependent issue, and also replacing use of MagProfile with !EMFProfile.
Number 4 above has been the trickiest to solve, and I'm not convinced that we have it right yet. One issue is that you cannot step on magnetic fields because that could create divergence. So what the code does is add a certain amount of field to the outflow object. However, you can't just keep adding your initial magnetic field profile because then the field will build up and you'll have too much magnetic energy/pressure.
The amount of field you need to add is equal to the amount of field that has advected out of the outflow object. Everything in the outflow object is done in terms of the EMF, so the EMF source that needs to be in the outflow object is as follows:
This should, in theory, take care of time-dependent issues. One other issue that remains unresolved is that the initial Bphi profile in the chombo does not match the Bphi profile that I supply.
Many of these points look ok, but I'm not sure why the profile has a discrepancy around r = 0.5. This might be solved by increasing resolution or subcycling.
Below are low resolution simulations with the current implementation. Much better than before, and maybe good enough, but I want to figure out the Bphi issue first.
Current Disk Usage
/clover: 11 T johannjc 3.0 T shuleli 2.5 T /alfalfa: 5.3 T shuleli 1.6 T martinhe 1.1 T johannjc 422 G bliu 262 G ehansen 17 G ckrau 13 G /bamboo: 13 T madams 3.2 T shuleli 2.4 T erica 2.2 T johannjc 1.1 T bliu 973 G /grass: 5.5 T erica 3.4 T shuleli 762 G johannjc 119 G
Update on Wire Turbulence Study (04-22-2015)
I have been updated the Wire Turbulence Studies page I have made daily.
These runs have been taking a while. The dt in the time step is very small. Jonathan asked me to try to increase the mindensity in the physics.data file.
MinDensity = .01 !1d-10 ! Minimum computational density before protection is triggered [1d-10]
I noticed there is something weird with the runs for which I did that, either from the beginning of the simulation, or at the point where I introduced the new mindensity during the restart.
Run | B10-T.03-M20 | B10-T.03-M3 | B10-T.01-M3 | B.1-T.03-M20 |
GIF | GIF | GIF | GIF | GIF |
It looks like a boundary conditions problem, however I have not changed the Gmthbc for any of these runs:
Gmthbc = 1,2,2,1,2,2
Update 04/20/2015 - Baowei
- Planet Wind 1 high resolution results for different lambda and different ns/np — smaller ns/np will make the planet wind relax faster. https://astrobear.pas.rochester.edu/trac/wiki/u/bliu/PlanetaryWind
Update 04/20/2015 - Eddie
The 3D runs for the Mach stem paper are going on several machines. Here is the status of all the runs (d = separation distance in rclump, a = angle in degrees):
Run | # Clumps | Defining parameters | Machine | Frames (out of 50) |
---|---|---|---|---|
A | 2 | d = 2.5, M = 7, 10 | Stampede | 43 |
B | 2 | d = 4.5, M = 7, 10 | Gordon | 16 |
C | 2 | d = 6.5, M = 7, 10 | Gordon | 24 |
D | 2 | d = 2.5, M = 7, 15 | Gordon | 15 |
E | 2 | d = 4.5, M = 7, 15 | Gordon | 11 |
F | 2 | d = 6.5, M = 7, 15 | Gordon | 11 |
G | 3 | d = 2.5, a = 30 | BlueHive | 0 |
H | 3 | d = 2.5, a = 60 | BlueHive | 32 |
I | 3 | d = 2.5, a = 90 | BlueHive | 0 |
J | 3 | d = 4.5, a = 30 | BlueStreak | 36 |
K | 3 | d = 4.5, a = 60 | BlueStreak | 34 |
L | 3 | d = 4.5, a = 90 | BlueStreak | 15 |
M | 10 | BlueStreak | 0 |
Notes:
- The Mach numbers are now all integer numbers: 7, 10, and 15.
- The "slow" clump remains at M = 7 for all runs. This means that the runs that contain an M = 15 clump now have larger domains to accommodate the faster clump with the same slow reference frame.
- The 10 clump run is coded up now. I will test it and run this one after all the others are finished.
- Data transfer is holding me up. I need to get the chombos/bovs from Stampede to Clover so I can visualize but this is taking a very long time.
- I have reservations on BlueHive and BlueStreak starting on Wednesday at 12:01 AM. I could only get 10 nodes on BlueHive and 512 nodes on BlueStreak for one week. This is not much compared to what I have access to already, but I won't have to wait in the queue or do so many restarts which is nice.
Below are column density maps and emission maps from run A.
Other Stuff
- I will do some more paper edits today and tomorrow and post them.
- I submitted an abstract for the CIRC poster session, and I have my poster so that's all set.
- I typed up a quick abstract for AstroNUM and attached it to this blog post. Haven't submitted it yet, please let me know what you think.
draft
draft
Grid of Runs for Planetary Wind Problem
Lets begin with all runs having
Hot Jupiter
Stellar Wind
km/s
Other parameters.
Resolution: 60 zone/planet radius
Allow planetary wind to fill grid for 2 of its own crossing times ( L_grid/V_pw(max) ) then turn on stellar wind and let it run for 10 crossing times ( L_grid/V_s )
# | lambda | Solar Wind Mach # | |
1 | 1.01 | 0.8 | 10-4 |
2 | 1.01 | 0.8 | 10-4 |
3 | 1.3 | 0.8 | 10-4 |
4 | 1.3 | 0.8 | 10-4 |
5 | 1.01 | 5.0 | 10-4 |
6 | 1.01 | 5.0 | 10-4 |
7 | 1.01 | 5.0 | 10-6 |
8 | 1.01 | 5.0 | 10-3 |
9 | 1.3 | 5.0 | 10-4 |
10 | 1.01 | 0.8 | 10-3 |
Parameters of binary simulations
Stellar rotation speed is important in this simulation. Research in the rotation period and gyrochronology has been down in M35 (S., Meibom, R., D., Mathieu and K., G., Stassun, 2009). The picture in the paper shows the period distribution of K and M type of stars.
The sun's rotation period is 24 days and the orbital period at sun's surface is 0.12 days. Therefore the picture above shows that:
\[\frac{\Omega}{\omega_{orb}}=\frac{P_{orb}}{P}\in(0.001,1)\]
and Skumanich's law is the relation between angular speed and the age of the star:
\[\Omega\sim t{-0.5}\]
The older the star, the slower it rotate. K and M type of stars may be low-mass orange dwarfs or K type giants, but the giants can not have such short rotation period, it must be main sequence low-mass orange stars. That's the tool the paper use to determine the age of the star but it also provide us some data.
Basically, the period of AGB or giants cannot be measured at this time because the Doppler effect would be too small at its photon sphere and they indeed rotation very slow. Thus its rotation speed is negligible compared to the break up speed. Here we assume there is a third star that is so close to the primary and it spin up the envelop of the primary's.
So we need some other mechanism to give the envelop angular momentum. In N., Soker and A., Harpaz, 2000, they give a formula of the ratio of the final angular frequency and break up angular frequency for horizontal branch stars:
\[\frac{\Omega}{\omega_{Kep}}=0.1(\frac{M_p}{0.01M_{env}})(\frac{a_0}{R}){0.5}(\frac{\alpha}{0.1}){-1}\]
where
, is the coefficient that determined by the envelop's density profile.is the mass of the planet (Jupiter like or BD), is the mass of the envelop's. is the initial separation of the planet and the star. is the radius of the star.
An example can be:
, , and
then the envelop will reach half of the break up speed.
I just realize that
\[\frac{\Omega}{\omega_{orb}}\]
should incorporate the effect of
for gas. (Or should I? Because observation can not determine the )# | separation | equatorial open angle | mass ratio | temperature | links to movie | ||||
1 | 4AU | 2pi/9 | 0.5 | 100 K | 0.67 | 0.005 | 1.mov | ||
2 | 4AU | pi/3 | 0.5 | 100 K | 0.67 | ||||
3 | 4AU | pi/3 | 0.5 | 100 K | 0.05 | 0.8 | contour | 4AUtopview | 4AUsideview |
4 | 4AU | pi/3 | 0.5 | 200 K | 0.2 | ||||
5 | 3AU | pi/3 | 0.5 | 100 K | 0.05 | 1.08 | 3AUsideview | 3AUtopview | |
6 | 3AU | pi/3 | 0.8 | 100 K | 0.2 | 1.08 | |||
7 | 4AU | pi/3 | 0.5 | 100 K | 0.05 | 1.08 | 4AUsideview | ||
8 | 6AU | pi/3 | 0.5 | 100 K | 0.05 | 1.08 | 6AUtopview | 6AUsideview | |
9 | 4AU | pi/2 | 0.5 | 100 K | 0.67 | 0.0005 | contour | 9_sideview.gif | |
10 | 4AU | pi/3 | 0.5 | 400 K | 0.2 | 0.8 | https://www.youtube.com/watch?v=tefDvaO4o9k | ||
11 | 4AU | 7pi/9 | 0.5 | 100 K | 0 | 1.14 | |||
12 | 3 AU | pi | 0.5 | 100 K | 0 | 1.08 |
(1, 2, 9) can investigate the effect of open angle.
(2, 3, 7) can investigate the effect of angular momentum of the equatorial outflow.
(3, 4, 10) can investigate the effect of temperature.
(3, 5, 8) can investigate the effect of separation.
(5, 6) can investigate the effect of mass ratio.
Meeting Update 04/13/2015 - Eddie
Not much to report today, because I'm still waiting on a bunch of simulations.
- Had issues with the Stampede file system over the weekend, so I had to restart the 2-clump, Mfast = 10 runs. They are still waiting in the queue, but they are almost 50% done.
- On bluehive, I'm running some low resolution tests before moving to production runs:
- 2-clump, Mfast = 15, Mslow ~ 11 which will replace the Mfast = 5, Mslow = 3.66 runs
- 3-clump, Mfast = 15, Mmed = 10, Mslow = 7
- Will continue writing the sections of the paper that I can without results
Meeting Update 04/13/2015
- OutflowWind rotating frame
- gamma=1.01; Planet is ~0.1 AU away from mass center. Star orbital angular velocity 0.5
- Resutls: density; temperature
- 2.5D ambient as stellar wind double-check
- gamma=1.01; set ambient density,velocity and temperature same as the stellar wind
- low res results — no turbulence and falling back at the tail: density and temperature; zoomed in
- High res results:3AMR
Angular Momentum and Sink Particles
Looked over Federrath Outflow Model and noticed significant issues with how we handle angular momentum in sink particles.
Both of us update the mass, momentum, and center of mass the same
But our treatment of angular momentum has glossed over some subtleties
A better treatment is to conserve total angular momentum.
where
and S is the spin.
which if we approximate to first order in
which implies our method is only accurate if
andFall Back Disk - Circumbinary Disk
Binary separation = 4 AU
Binary orbital period = 6 years
Simulation temp = 100 K
Primary:
Mass = 1 SM
Quiet phase outflow = 1E-7 SM/yr
Equatorial outflow open angle = 20 degree
Equatorial outflow = 2E-5 SM/yr
It is calculated in this way:
\[\dot{M}=\rho S v\]
\[S=2 \pi r2 \times 2\int_{7\pi/18}{\pi/2} \sin\theta d\theta\]
\[v=\sqrt{v_r2+(\Omega r)2}\]
Equatorial outflow density is increased by 1000 times compared to quite phase, varies with time as the picture below shows, it last for 10 years. The total mass ejected from the equatorial outflow is:
\[M=\delta t \times \dot{M}=2\times10{-4}M_{\odot}\]
Equatorial outflow has angular velocity of:
\[\Omega=1.36294\times10-7 s{-1}\]
This is equivalent to a period of 1.46 years. Since the equatorial outflow is launched at 1 AU, this means
.Secondary:
Mass = 0.5 SM
The simulation:
This is a contour plot of density. The more transparent the color, the low the density. Therefore, red, which is on the top of the color bar, is the most opaque color and represent the highest density, the next highest the grey and so on.
We can see there is a primitive disk like structure around two stars. Spiral shocks compress the disk, push the disk outwards and give it angular momentum.
Since the angular frequency of the equatorial outflow is below the orbital angular frequency and we still see a circumbinary disk, it means that, it the outflow is induced by the in fall on a planet.
Future works:
- Include mass, momentum and angular momentum conservation of the sink particles so that we can study the orbital evolution.
- Add cooling. NK cooling preferred (D., Neufeld and M., Kaufman (1993)).
Page for wire turbulence study
Find it here: WireTurbulenceStudies
Mach stem papers
There are only 2 relevant papers on ADS with "Mach stem" in the title:
- My HEDLA conference paper: http://adsabs.harvard.edu/abs/2014arXiv1412.6495H
- Kris' hysteresis paper: http://adsabs.harvard.edu/abs/2013HEDP....9..251Y
More results are found by searching for "Mach reflection":
- Kemm 2014, interesting paper on double Mach stem problem: http://adsabs.harvard.edu/abs/2014arXiv1404.6510K
- Ivanov 2002, nice reference on regular vs. Mach reflection, experimental and numerical: http://adsabs.harvard.edu/abs/2002ESASP.487..341I
Other references refer to optics, water waves, detonation waves, etc. I don't think these are very applicable to our scenario. Papers 1-4 might have enough references within for me to use. Paper 4 has led me to some purely fluid mechanics papers:
- Chpoun 1995, wind tunnel experiments regular vs. Mach reflection, http://adsabs.harvard.edu/abs/1995JFM...301...19C
- Hornung 1982, more on transition from regular to Mach reflection: http://adsabs.harvard.edu/abs/1982JFM...123..155H
- Ivanov 1995, more on hysteresis effects: http://adsabs.harvard.edu/abs/1995PhFl....7..685I
There are a few books that give good descriptions of shock reflection phenomena which I have used as references in my previous paper as well:
- Ben-Dor: Shock Wave Reflection Phenomena
- Courant and Friedrichs: Supersonic Flows and Shock Waves in Applied Mathematical Sciences
- Landau and Lifshitz: Fluid Mechanics
There is also a derivation of the critical angle in this paper:
- de Rosa 1992: http://adsabs.harvard.edu/abs/1992PhRvA..45.6130D
Ben-Dor is definitely seems to be a key player in shock-wave reflection. He has a lengthy review paper just on shock wave reflection in different regimes:
- Ben-Dor 1988: http://adsabs.harvard.edu/abs/1988PrAeS..25..329B
I've skimmed through these papers, but I will be reading them more closely over the next few days as I start to write my paper. Some key things that I have learned:
Historically, there have been two minimum critical angles for Mach stem formation based on two different criterion. The one that I have been using is known as the "detachment" criterion (alpha_d), and it is commonly cited in text books. The other criterion is the "von Neumann" or "mechanical equilibrium" criterion (alpha_n), and it is a pressure argument which leads to an angle less than the detachment criterion. Many papers talk about a dual-solution region where the angle is alpha_n < alpha < alpha_d. In this region, both regular reflection and Mach reflection are theoretically possible. The figure below illustrates this point:
Ben-Dor's review paper suggests that for steady flows, the mechanical equilibrium criterion is the correct one to use, and that regular reflection only occurs above alpha_n in what he calls psuedo-steady flows. One complication that I see is that alpha_n is much more dependent on M than alpha_d. Thus, strong shock approximations might not be appropriate, and I would have to use a critical angle formula that is also dependent on M. I have yet to find a good formula for that like I did for alpha_d.
Ben-Dor also has an region showing no reflection in his figures, but there doesn't seem to be a good explanation of the transition criterion or a useful formula anywhere.
Experimental results vary. Some agree with alpha_n, and some agree with alpha_d, while others show hysteresis effects and support the dual-solution region argument.
In conclusion, I might be quoting a minimum critical angle of ~ 40 deg, when in actuality it could be more like 30 deg. The issue could be avoided for Mach numbers ~ 2.2, since alpha_n = alpha_d at this M. However, then I would have to use the exact formula for alpha_d that contains M since a strong shock approximation would no longer be valid for M=2.2.
Meeting Update 04/06/2015
- OutflowWind module
- gamma=5/3 no wind
- Rotate:
1)reproduce corotate binary;
2)Roate_no_wind; Rotate works for 2.5D ?
3)Weak_wind_along x?;
4) set the wind(stellar) direction along y
Proposed fiducial parameters for Wire Turbulence study
Proposed fiducial parameters for Wire Turbulence | |
---|---|
Wire diameter (in units of wire spacing) | .25 |
Mach Number | 20 |
Beta | 1 |
Gamma | 1.1 |
Field orientation | y (flow is in x and mesh is aligned with yz unit vectors) |
wire density (in units of wind density) | 1000 |
ambient density (in units of wind density) | .01 |
distance from wind boundary to wire center (in units of wire spacing) | 1 |
Number of zones per wire spacing | 1024 |
Some Planetary Wind Tests
Results of some \gamma = 5/3 runs.
Test Run A: standard. 200x400. 1 level of refinement. Tmax = 1.0
Test Run B: standard. 100x200. 0 level of refinement. Tmax = 1.0
- seeing weird behavior?
Test Run C: standard. 100x200. 0 level of refinement. Tmax = 4.9; 50 frames
- seeing weird behavior.
- Wind does not expand uniformly
attachment:C_0404015_rho.gif
Test Run D: standard. 100x200. 1 level of refinement. Tmax = 4.9; 50 frames
- no change from C
attachment:D_0404015_rho.gif
Test Run E: standard. 200x400. 1 level of refinement. Tmax = 4.9; 50 frames
- no change from C, just better resolution
- is it the \gamma=5/3 or something about timing? Almost looks like the wind is turning off.
attachment:E_040415_Full.gif
Test Run F: standard. 200x400. 1 level of refinement. Tmax = 4.9; 50 frames
- \gamma=1.01
Sensitivity of H-alpha Routine
The H-alpha routine in the code determines the amount of H-alpha emission (neutral hydrogen de-excitation line from level 3 to level 2). It is dependent on electron number density, temperature, ionization fraction, and hydrogen density. There are 2 main components to this emission: a recombination term and a collisional excitation term.
Collisional excitation is what most people usually think about when talking about H-alpha emission: electrons collide with neutral hydrogen and the excited hydrogen atom goes through the H-alpha transition. The recombination term comes into play when you have ionized hydrogen that wants to recombine. When it recombines, there is a nonzero probability that it goes through the H-alpha transition. So, if you have a lot of ionized hydrogen at a low temperature, there should be a lot of recombination and the recombination term can dominate the H-alpha emission.
Below are two plots showing the H-alpha emission as a function of temperature for two different ionization fractions. The hydrogen number density is the same in both plots, and for simplicity the gas was assumed to be all hydrogen. nH = 5560 was chosen because it is the average of the post-shock densities of two simulations which I will talk about.
Notice how in the x=0.99 plot, the recombination term is much stronger and dominates the H-alpha at low temperatures longer. This is because 99% ionization is unrealistically high for these temperatures, so the gas should rapidly recombine and thus go through many more H-alpha transitions.
3D Mach stem simulations
For reference, see ehansen02242015.
For the 3D Mach stem simulations, we noticed a significant decrease in H-alpha emission from the M = 10 runs to the M = 5 runs. The best way to explain the difference is by comparing the temperatures in a couple of these runs. Left is M = 5, right is M = 10, both at separation distances of 3 rclump.
I measured the temperatures in the bow shock intersection regions and got about 6900 K for M = 5, and 19,700 K for M = 10. This is approximately a factor of 3 which seems reasonable given that postshock temperatures are approximately proportional to vshock2
These simulations have relatively low ionization fractions throughout, and at the bow shock intersection points, the ionization remains near 0.01 for both runs. So we should look at the x = 0.01 H-alpha plot.
It would appear that these simulations have hit the sweet spot. There is a sharp increase in the H-alpha going from 6900 K to 19,700 K. This is where the collisional excitation term takes over and dominates. I estimate a jump of about 2 or 3 orders of magnitude which is consistent with the emission maps.
Also, look at the maximum temperatures in the simulations, and then look at the x = 0.01 H-alpha plot. This explains why the M = 5 model does not show a wide range of H-alpha intensities like the M = 10 model does.
In conclusion, the H-alpha routine is working just fine. The parameters I chose for the simulations just happen to illustrate the temperature sensitivity in our code.
Binary Star Fall Back Disk Variable Set
Here are variables which go into the HAFBAD problem (Highly Abstracted Fall Back Disk)
Stars
Primary Mass (solar masses):
Binary Mass ratio:
Orbital Separation:
Outflow in the form of a brief Shell
Velocity of Shell (< escape velocity but >> sound speed:
Density of Shell:
Temperature of Shell (with sound speed constraint:
Rotation of Shell:
Duration of shell ejection
Simulation
Resolution defined in terms of injection radius R_0:
fall back disk
Assume a third star with mass between 0.001 SM to 0.05 SM (hot Jupiter to Brown Dwarf) around the primary fall into the primary during common envelop phase evolution from an orbit of 1 AU to 0. So it will induce an equatorial outflow J. Nordhaus and E. Blackman 2006. The angular momentum, the potential energy and the kinetic energy of the planet will partially be transferred to the outflow. A detailed calculation will not be presented here.
Let's assume the equatorial outflow is
, lasting for 1 - 10 yrs, so the total mass of the outflow is .Below are some parameters:
: Binary separation, can vary from 4 AU to 20 AU. A close binary will transfer more angular momentum to emitted gas but will also challenge the small planet around the primary. Because they are basically three body system. (Boldly assume the small planet is 0.8 - 1 AU from the primary such that the secondary will impose so much effect.)
: Primary mass, let's fix it to be 1 SM and low metallicity. Metallicity is important in AGB evolution since it will provide dust thus drive the wind away from the primary.
: Ratio of secondary's mass to primary's mass. It can be 0.2 to 1, varies in correspond to the binary separation. The greater the ratio the less the orbital period thus the more angular momentum will be transferred to the gas emitted from the primary. Note here when the gas is emitted from the primary, it receives some angular momentum from the primary because the primary has orbital motion.
: Open angle of the equatorial outflow. Let it vary from . The more mass the planet (or Brown Dwarf), the wider the open angle. This open angle is closely associated with turbulence driven by the in falling planet rather than viscosity. companion driven dynamo happens in this process.
: Spinning angular velocity of the primary. The envelop of the primary is spun up by the in fall of the third planet and may induce a dynamo. This may also come from its own spinning (though very small, data missing here). A value to be noticed is the orbital angular velocity at the radii which the outflow will be emitted. .
Given that:
\[\omega=\sqrt{G M/r3}\]
\[l_z=r\times v=r2 \omega = \sqrt{G M r}\]
This means, any object with angular momentum larger than 0 can find its orbit, mathematically, if the radial velocity is 0.
: Actual temperature can be around 2000 K, like the temperature in the envelop. In the test simulations, we can choose 200 K simply because it reduce the pressure and the sound speed.
: Radial velocity of the outflow. The best scenario is we can depict its angular dependence. Here we can use step function due to low angular resolution in test runs. This speed should below escaping velocity. Let's say, ¾ of the escaping velocity.
: The radius which the outflow is emitted. It is set to be 1 AU.
A recent simulation:
Primary star has a constant alpha of 0.1339 throughout space. Simulation temp is 100 K.
Proposed Runs for 3D Mach Stems
2-clump runs
We vary the separation and Mach numbers.
Run | Separation (rclump) | Mfast, Mslow |
---|---|---|
A | 2.5 | 10, 7.32 |
B | 4.5 | 10, 7.32 |
C | 6.5 | 10, 7.32 |
D | 2.5 | 5, 3.66 |
E | 4.5 | 5, 3.66 |
F | 6.5 | 5, 3.66 |
On bluehive, I used about 3800 SUs per run for 40 cells/rclump resolution. Adding a level of refinement to get 80 cells/rclump will increase SUs by approximately a factor of 16 (a little less due to filling fractions). Accounting for the fact that Stampede cpus are slightly faster than bluehive (2.7 GHz vs. 2.4 GHz) leads to approximately 52,000 SUs per run. For 6 runs, this is a total of 312,000 SUs.
3-clump runs
Clumps stay at M = 2.5, 5, 10 for all runs. The M = 2.5 and M = 5 clumps remain in same plane at a fixed separation distance, and the M = 10 clump is positioned at some angle and distance from the M = 2.5 clump.
Looking down the barrel, the 3-clump setup is like this:
Run | Separation (rclump) | 3rd Clump Position (deg) |
---|---|---|
G | 2.5 | 30 |
H | 2.5 | 60 |
I | 2.5 | 90 |
J | 5 | 30 |
K | 5 | 60 |
L | 5 | 90 |
We would get longer run times because filling fractions are higher with 3 clumps. My best guess is that these runs could take about 77,000 SUs each. A total of 462,000 SUs for 6 runs.
10-clump run
Clumps at random positions with varying separations and velocities, but all moving in same general direction.
With many clumps on a large grid, this could be a very expensive run. I'm guessing it would be at least 100,000 SUs if we still want 80 cells/rclump.
In total, my crude estimates give on the order of 900,000 SUs for the entire study. Again, we're playing with a factor of ~16 here if we want the increased resolution. So if we keep the resolution at 40 cells/rclump, then we can do all these runs with less than 60,000 SUs. Maybe we only do the 2-clump runs at the increased resolution and the rest at the original, lower resolution.
Equatorial outflow
After an equatorial outflow from the primary star.
Bounding box is 192 AU * 192 AU * 48 AU
Primary mass = 1 SM
Seconary = 0.5 SM
alpha of primary = 0.1339 throughout
Isothermal temp = 1500 K
post wind has pushed the dense shell outwards