Plan
- Running model 3 now. Model 2 has finished. I will post model 2 and model 3 result on Friday.
model number | q | z | burst | separation | link |
1 | 0.5 | 0.031 | no | 4AU | model1 |
2 | 0.5 | 0.031 | yes | 4AU | model2 |
3 | 0.1 | 0.016 | no | 3AU | model3 |
4 | 0.5 | 0.008 | no | 6AU | model4 |
5 | 0.5 | 0.1 | no | 3AU | model5 |
- I am thinking about how to analyze the results. Also, three dimensional single AGB star is not so symmetric (see the wind velocity profile summary). Implying that there should be something not right in the code.
Tests of the feedback routine
Using the Bondi routine with the particle's position fixed:
Mass conserving? | yes |
Can it run in AMR? | |
Can it run on multiple processors? |
Using the Bondi routine with the particle's position fixed:
Mass conserving? | yes |
Momentum conserving? | |
Angular momentum conserving? |
Using the uniform collapse module:
Meeting Update --7/27/16
Wire Turbulence
Update 7/27
- Ran simulation with lambda=5 for longer time, to get to steady state:
- Numerical integration technique for isothermal Parker wind - starting from differential equation in terms of u and r, parameterize in terms of s, so that
We can calculate the Jacobian for this dynamical system, and evaluate it at the critical point:
This gives the linear approximation to the system near the critical point; calculating the eigenvalues and eigenvectors, travel a small distance along the eigenvector with positive eigenvalue (corresponding to tangent to the stable manifold of the original system) and use numerical integrator of choice outward from the critical point.
- Also started determining B fields for planet and star - looking for sigma of both to be 0.1 initially (based on Jonathan's example parameter space).
Meeting update
Fixed a few bugs with the outflow feedback routine and test module. Outflow looks to be qualitatively reasonable. Next have a suite of tests to run to try and break the new routines.
For now, am starting with the bondi module which prescribes a bondi flow onto the origin given a central mass. The central sink particle is .4 solar masses. It is accreting via bondi accretion routines and is launching an outflow at .50 efficiency (that is, 50% of the accreted mass is launched in the outflow). The opening angle is 30 degrees.
Update 7/20
- Looked again at the Parker wind differential equation:
If we assume spherically symmetric outflow in an isothermal corona, then by conservation of momentum, continuity, and the ideal gas law, we have
Continuity shows that
, where C is a constant. If we substitute for p and rho into the momentum equation, we get a differential equation in u and r:Define
so that (i.e., sonic radius); substitute and rearrange to findIntegrate both sides to find the Parker wind equation,
- Modified the domain on a couple of simulations to better resolve the planet/capture sonic radius. Also plotted contour of sonic radius.
Value | Mdot (g/s) | ||
lambda=5 | 5.6x10^{11} | (unchanged) | |
M_{p} = 0.25M_{J} | > 1.67x10^{9} | (not well resolved) | 3.7x10^{8} |
T_{p} = 5x10^{4} K | Negative throughout | (not resolved) | 3.4x10^{8} |
For lambda=5, may not reach quasi-steady state.
T_{p} = 5x10^{4} K:
M_{p} = 0.25 M_{J}:
Sonic radii:
- Continuing to work my way through Zel'dovich.
Meeting Update --7/20/16
- Schedule for the Visitor
- office: 476?
- No hotel shuttle on weekends. Transporting (Zhuo?)
- will send out schedule draft soon
- New external hard disk for archiving (~160USD for 4TB, 250USD for 6TB)
- Wire Turbulence
- Mixture distribution: two types of materials (wire & wind)? — scatter plot: https://astrobear.pas.rochester.edu/trac/wiki/u/bliu/mixtureDistribution
Using AstroBEAR for Lab Experiments
AstroBEAR needs more physics in order to make it more applicable to lab experiments. I've been toying with this idea, and I don't think I'm going to get very far because it is quite complicated. However, I have been somewhat successful at hard-coding some 0th order approximations into a separate problem module (not usable by the rest of the code). Even the act of attempting to implement these things has forced me to learn about new physics, yay!
Laser Deposition
A simple dump of some fraction of laser energy onto a target surface (resonance absorption). Target heats up and the thermal pressure gradient launches a plume of plasma.
One of the tricky things here was the directional nature of the laser heating. Once the laser hits some critical density, it is only absorbed to a certain depth. Thus far, my code only works with single proc and a fixed grid.
A more clever implementation could make use of multiple processors and AMR. And to do laser deposition more accurately, you'd need ray tracing.
Now, the fun MHD stuff…
Biermann Battery
In ideal MHD, we have the induction equation:
.
In resistive MHD, the induction equation becomes:
This is in some version of the code thanks to Shule. There are other MHD terms that can be important depending on the context of the problem. Here is the induction equation with the Biermann term:
The Biermann term makes self-generated B-fields possible. In other words, you don't need an initial B to generate more B. This is also in my new problem module, and could someday be implemented into AstroBEAR for general use. This formulation assumes the plasma is fully ionized (i.e.
).
Nernst Effect
Lastly, we have the Nernst term. Here's how the induction equation becomes even more complicated:
This
is the Nernst velocity, and it can be written as:
where
is the heat flux and is the thermal conductivity. is technically a tensor since it is different in different directions (anisotropic). We can simplify the tensor and write as follows:
where
is a unit vector in the direction of . The first two terms are implemented somewhere in the code for conduction, but the third term is missing. Furthermore, I believe the conduction implementation in the code is strictly for transporting heat; in other words, there is no magnetic field generation.This is something that I haven't implemented yet. There should be a way to add the Nernst velocity to the fluxes within the code, but I think I will first try to separate
in the induction equation and treat it as a source term.Plan
I will run model 3, model 4 and model 5.
model number | q | z | burst | separation | link |
1 | 0.5 | 0.031 | no | 4AU | model1 |
2 | 0.5 | 0.031 | yes | 4AU | model2 |
3 | 0.1 | 0.016 | no | 3AU | model3 |
4 | 0.5 | 0.008 | no | 6AU | model4 |
5 | 0.5 | 0.1 | no | 3AU | model5 |
q is the mass ratio of the two stars, defined by
z is the tidal force coefficient, defined by
where is the binary separation.Result
- Jonathan helped in debugging. The code can restart on bluehive without producing error.
- I ran a 4 AU simulation this week. This simulation has a pulsating AGB star which described before. The difference is: I injected an outflow during 16.83 year and 26.83 year. So the total outflow time window is 10 years. Then the AGB star is set back to normal pulsating state. The result is very exciting. Below are three density plot and contour plot from the simulation. The burst of the outflow has ended and the AGB star is back to pulsating state. The speed of the outflow is 20km/s (remember the pulsating amplitude is only 5km/s), and the escape velocity (with respect to the AGB star only) where the outflow is ejected is about 40km/s. However, we can not treat the gas as ballistic matter. The pressure and temperature will push much of the gas out.
A large amount of the gas is gravitationally bound. I mistakenly overwrite the result of later time. A new simulation is running. I saw the late time evolution, a smooth circumbinary disk will form around the binary. (click the image to enlarge)
A disk forms around the secondary. (click the image to enlarge)
- Baowei helped me to set up the code on Bluestreak. I can use 8192 cores on bluestreak. The 3 AU, q=0.1 simulation is schedules to run on bluestreak.
Plan
- Plan for runs
mass ratio | 3AU | 4AU | 6AU | 10AU(maybe) |
q=0.1 | yes | unknown | no | no |
q=0.5 | no | yes | yes | maybe |
Result
- I was improving the code's efficiency and debugging. I restrict the highest level only to the AGB star and the secondary. The radii of dust formation zone is less than 2.4 AU. The finest cells consist of two parts, one is the cells that are within the 2.4 AU sphere of the AGB star and the other part is the cells that are within the cylinder which has 2.4 AU height and 3 AU radius. AMR is set to capture the spiral shock, too.
- The current simulation is: 96AU*96AU*48AU physical size with 5 level of AMR. The base resolution is 60*60*30. The coarsest cell is 1.6 AU and the finest cell is 0.05 AU. The AGB star has 18 cells across its radius and 21 cells across its photosphere. The simulation will take 8 days on 120 cores, correspond to 12 orbits, or 72 years.
- The 4 AU separation run is still going on.
meeting update
sad eddie is gone
other than that, finished the newest draft and am next turning towards testing the outflow module jonathan installed.
Meeting Update --07/07/16
- Wire Turbulence
- Updates of figures
- Updates of figures
1) square-root (standard?) velocity variance: different variance for x & yz comes from the x-direction flow?
2) redid all PDF figures with the new weighted-with-area histogram data from visit: still cannot do the time average due to the different values of x-axis/frequency values
3) redid the energy with total pressure instead of thermal energy, although the total pressure is so small due to gamma=1.001. Pressure rather than the thermal energy matters here?
4) magnetic energy with
plot: will do
5) physical meaning of Gaussian 2 fit: Mixture of two types of turbulence or Turbulence with two components
- Post processing Spectra Results
- hydro spectra data: running
- linear fit: will do
Update 7/7
- Solving parker wind equation,
psi − ln(psi) = 4 ln(xi) + 4/xi − 3, where psi = (v/v_{s})^{2} and xi = r/r_{s}.
There are 4 types of solutions - combinations of sub/supersonic at the surface and v → 0/v_{f} as xi → infinity. The only physical solution is subsonic at the surface and has a nonzero velocity infinitely far away, and passes through xi = 1, psi = 1.
Attempts at solutions using approximations and a single solve:
Alternate implementation of approximation from Jonathan's code; works for r > r_{s}, but not r < r_{s}.
Pertinent line is: yy(i)=vpasolve(y - log(y) == -3 + 4*log(xx(i))+4/xx(i),y,(xx(i) > 1) * xx(i) + (xx(i) ⇐ 1)*exp(3-4*log(xx(i))-4/xx(i)))
Interpreted as:
if xx(i) ≥ 1
yy(i)=vpasolve(y - log(y) == -3 + 4*log(xx(i))+4/xx(i));
else
yy(i)=vpasolve(y - log(y) == exp(3-4*log(xx(i))-4/xx(i)));
end
Correct & incorrect plots:
- Calculated mass loss rate in VisIt. Used a scalar expression for the radial mass (momentum) flux and summed over isosurfaces (can also use spherical slices) to get mass loss rate. Plots were unfortunately uniform in display; however, in experimenting, discovered the flux operator makes plots that appear qualitatively correct for the flux:
Mass loss rates for the various simulations from last time:
Value | Mdot (g/s) | |
Aniso | 6.1x10^{9} | |
lambda=5 | 5.6x10^{11} | |
lambda=15 | 3x10^{7} | |
M_{p} = M_{J} | 9.9x10^{9} | |
M_{p} = 0.5M_{J} | 1.5x10^{9} | |
M_{p} = 0.25M_{J} | > 1.67x10^{9} | (not well resolved) |
T_{p} = 5x10^{4} K | Negative throughout | (not resolved) |
T_{p} = 5x10^{3} K | 1.4x10^{9} | |
T_{amb} = 3 K | 4.5x10^{9} | |
T_{amb} = 25 K | 4.5x10^{9} | |
T_{amb} = 50 K | 4.5x10^{9} |
Aniso appears correct, based on plot from planet paper.
Working on calculating the sound speed to get the sonic radius.
- Started working through occasional prelim problems, as well.
Plan
- About the simulation zone. The physical size of the current simulation has . However, the hydrodynamic behavior suggest that the gas will go out further and will take a long time to fall back. I estimate that a minimum physical size would be for and binary separation. The simulation zone may still need to be larger for and separation. The physical size of the finest cell should be or less. Therefore, a 4 level AMR with base grid can handle the simulation. The time should take (estimated) more than 15 days on 120 cores (3 level AMR is 5 days). If I can run with 2*120 cores on bluehive. I think the simulation can be done by the end of August.
- The pulsation period is 1 yr. This is the highest macroscopic frequency in this simulation (compared to orbital period and/or any other). To resolve the mass loss rate, I need at least 8 outputs per pulsation period. Therefore if the simulation is 300 years. The output will have 2400. Each output will take 200 MB at least. So the total disk space needed for one simulation is at least 500 GB. Currently I have 1200 GB on bluehive.
- Run a 4au simulation with high resolution. 96au*96au*48au
Result
- High mass loss rate AGB star need a large atmosphere. Its photosphere will be around 1.5au. Liljegren et al. (2016) use 1.5au AGB and achieved . However, 1.5au will be too large for us. The tidal force will be 50% of the gravitational force from the AGB star. The highest mass loss rate I can do is . Below is the mass loss rate in 3-D simulation. I sum the mass flux through a 10 AU spherical shell. The luminosity of this AGB star is and mass is . So the mass loss rate and size is between RGB and super AGB.
- A detailed description of this AGB star model is:
AGB star model
There are two most important aspects in binary star simulation. One is the radiative transfer and the other is the AGB star model. The first one determine the large scale motion and morphology while the second one determine the robustness of inner boundary.
A sketch of the AGB star's structure is shown in the picture.
The hard pulsating sphere has velocity variation. Its amplitude is 5km/s and its period is 1 year.
The radius of the photosphere is also fixed. In reality, it must change but it is difficult to calculate. Within the photosphere, reduced gravity (80% of the original gravity) is being used. The reduced gravity has its physical stand - higher opacity inside the star but it is also because the resolution is poor. The simulation is resolving 2~3 order of magnitude density drop within 3 cells.
Between the photosphere and the dust formation zone, it is the gas opacity
. The dust formation radii is calculated by
Where
is first calculated by assuming that there is only gas in the simulation. When the is determined, I update the opacity where there is dust. The dust opacity is .A questionable number is
.- A 3 AU separation simulation now will take about 5 days on 120 cores. I tested the resolution on single AGB star and it has a difference. The size of the simulation zone is . The previous resolution is base grid with 2 level of AMR. The current AMR is 3 while the base grid resolution is the same. The finest physical resolution is 0.05au thus the AGB star has 18 cells across its radii. The photosphere has 21 cells.
- I calculate the strength of tidal force by:
where
is the binary separation. Substitute , , and you will get the value is close to 0.5 which is 50%.- Currently I have some 3au and 4au simulations.
- 3au simulation:
I use
or at the first 23 years. The AGB star is the same as in the single star case. After the simulation has come to a pattern state - the spiral shock escapes from the system. Below shows the pattern. Left figure shows the density plot and the right figure shows the temperature plot.After 23 years, I make the AGB star dimmer by setting
or . The spiral shock is not escaping in this situation but the simulation box is too small so we can not see its fall back.The movie is here:
- 4au simulation:
Interestingly, the 4au separation binary star will have non-escaping spiral shock when I put the normal AGB star (
) in the simulation. Also, the simulation box is too small.The movie is here:
It seems that 3au simulation has Roche lobe overflow while 4au simulation has wind Roche lobe overflow?