1. Running model 3 now. Model 2 has finished. I will post model 2 and model 3 result on Friday.

2. 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.

Using the Bondi routine with the particle's position fixed:

Using the Bondi routine with the particle's position fixed:

Using the uniform collapse module:

Case 1. Box radius= bondi radius

movie​

movie​

movie​

Case 2. Box radius= 2xbondi radius

Movies:

Wire Turbulence

1. ​hydro & mhd post processing
2. ​two terms Gaussian fit for mhd

• 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:

Jacobian = \begin{bmatrix} 0 & 2c^3 \\ \frac{GM}{2c^3} & 0 \end{bmatrix}

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).

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.

• 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 r^2 \rho u = C, 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 r_c = \frac{G M m_p}{4 T} so that u_c^2 = \frac{2 T}{m_p} (i.e., sonic radius); substitute and rearrange to find

Integrate 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.

For lambda=5, may not reach quasi-steady state.

T_p = 5x10⁴ K:

M_p = 0.25 M_J:

Sonic radii:

• Continuing to work my way through Zel'dovich.

• Schedule for the Visitor
  1. office: 476?
  2. No hotel shuttle on weekends. Transporting (Zhuo?)
  3. will send out schedule draft soon

• New external hard disk for archiving (~160USD for 4TB, 250USD for 6TB)

• Wire Turbulence
  1. Mixture distribution: two types of materials (wire & wind)? — scatter plot: ​https://astrobear.pas.rochester.edu/trac/wiki/u/bliu/mixtureDistribution

• Problem for Bruce's JetClump module: ticket #446

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:

\frac{\partial \vec{B}}{\partial t} = \vec{\nabla} \times (\vec{v} \times \vec{B}).

In resistive MHD, the induction equation becomes:

\frac{\partial \vec{B}}{\partial t} = \vec{\nabla} \times (\vec{v} \times \vec{B}) + \frac{\eta}{4\pi} \nabla^2 \vec{B}

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:

\frac{\partial \vec{B}}{\partial t} = \vec{\nabla} \times (\vec{v} \times \vec{B}) + \frac{\eta}{4\pi} \nabla^2 \vec{B} + \frac{1}{e} \vec{\nabla} \times (\frac{\vec{\nabla} P}{n})

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. n = n_e).

Nernst Effect

Lastly, we have the Nernst term. Here's how the induction equation becomes even more complicated:

\frac{\partial \vec{B}}{\partial t} = \vec{\nabla} \times ((\vec{v} + \vec{v_N}) \times \vec{B}) + \frac{\eta}{4\pi}\nabla^2 \vec{B} + \frac{1}{e} \vec{\nabla} \times (\frac{\vec{\nabla} P}{n})

This \vec{v_N} is the Nernst velocity, and it can be written as:

\vec{v_N} = \frac{\vec{q}}{2.5 P} = \frac{-\kappa \cdot \vec{\nabla} T}{2.5 P}

where \vec{q} is the heat flux and \kappa is the thermal conductivity. \kappa is technically a tensor since it is different in different directions (anisotropic). We can simplify the tensor and write \vec{q} as follows:

\vec{q} = - \kappa_\parallel \vec{b}(\vec{\nabla}T \cdot \vec{b}) - \kappa_\perp \vec{b} \times (\vec{\nabla}T \times \vec{b}) -\kappa_\wedge (\vec{b} \times \vec{\nabla}T)

where \vec{b} is a unit vector in the direction of \vec{B}. 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 v_N in the induction equation and treat it as a source term.

I will run model 3, model 4 and model 5.

q is the mass ratio of the two stars, defined by M_{com}/M_{AGB}

z is the tidal force coefficient, defined by (\frac{M_{com}}{(d-r_{AGB})^2}-\frac{M_{com}}{(d+r_{AGB})^2})/\frac{M_{AGB}}{d^2} where d is the binary separation.

​linktosummary

1. Jonathan helped in debugging. The code can restart on bluehive without producing error.

2. 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)

3. 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.

1. Plan for runs

1. 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.

2. 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.

3. The 4 AU separation run is still going on.

sad eddie is gone

other than that, finished the newest draft and am next turning towards testing the outflow module jonathan installed.

• Wire Turbulence
  1. 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 \rho V^{2} 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
  1. hydro spectra data: running
  2. linear fit: will do

• Solving parker wind equation,

psi − ln(psi) = 4 ln(xi) + 4/xi − 3, where psi = (v/v_s)² 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:

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.

1. About the simulation zone. The physical size of the current simulation has 48au^3. 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 96au^3 for 3au and 4au binary separation. The simulation zone may still need to be larger for 6au and 10au separation. The physical size of the finest cell should be 0.05au or less. Therefore, a 4 level AMR with 120^3 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.

2. 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.

3. Run a 4au simulation with high resolution. 96au*96au*48au

1. 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 3\times10^{-6}M_{\odot}/yr. 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 10^{-7}M_{\odot}/yr. 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 2343L_{\odot} and mass is 1M_{\odot}. So the mass loss rate and size is between RGB and super AGB.

2. 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 (2\times10^{-4}cm^{2}/g). The dust formation radii is calculated by

R_d=\sqrt{\frac{\alpha L \exp{(-\tau)}}{4 \pi \sigma T_{sub}^{4}}}

Where \tau=\int_{rphotosphere}^{rmax} \rho\kappa_{gas}dr is first calculated by assuming that there is only gas in the simulation. When the R_d is determined, I update the opacity where there is dust. The dust opacity is 8cm^{2}/g.

A questionable number is \alpha=0.3.

3. 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 48au^3. The previous resolution is 120^3 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.

4. I calculate the strength of tidal force by:

(\frac{M_{com}}{(d-r_{AGB})^2}-\frac{M_{com}}{(d+r_{AGB})^2})/\frac{M_{AGB}}{r_{AGB}^2}

where d is the binary separation. Substitute d=3au, r=1.5au, M_{AGB}=M_{\odot} and M_{com}=0.5M_{\odot} you will get the value is close to 0.5 which is 50%.

5. Currently I have some 3au and 4au simulations.

1. 3au simulation:

I use T_{eff}=2900K or 2343L_{\odot} 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 T_{eff}=2600K or 1513L_{\odot}. 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:

densitytop3au.gif​

2. 4au simulation:

Interestingly, the 4au separation binary star will have non-escaping spiral shock when I put the normal AGB star (T_{eff}=2900K) in the simulation. Also, the simulation box is too small.

The movie is here:

densitytop4au.gif​

It seems that 3au simulation has Roche lobe overflow while 4au simulation has wind Roche lobe overflow?