| Home | Trees | Indices | Help |
|---|
|
|
The energy balance class.
Author: R. Lombaert, H. Olofsson & M. Maercker (Chalmers, Sweden)
For the use of the EnergyBalance module, please contact one of the three authors listed here.
Calculates energy balance and provides tools for reading and writing I/O, and setting input physics.
For now limited to the energy balance from the dust formation radius up to the outer radius.
An example for calculating the energy balance (see the respective functions for in-depth information): >>> from cc.modeling.physics import EnergyBalance as EB >>> eb = EB.EnergyBalance() >>> eb.iterT() >>> eb.plotT() >>> eb.plotRates() This calculates and plots the temperature profiles and heating/cooling rates across all iterations using the default settings given in cc.path.aux/inputEnergyBalance_standard.dat
Default in the inputEnergyBalance_standard file can be adapted through keywords passed to EnergyBalance. A few examples: >>> eb = EB.EnergyBalance(a=[0.005,0.25,300,1]) which creates a grain-size grid from 0.005 to 0.25 micron with 300 points distributed in log scale.
>>> eb = EB.EnergyBalance(hterms=['dg','dt']) adds dust-gas collisional heating and heat exchange to the adiabatic cooling as heat/cool terms.
>>> eb = EB.EnergyBalance(template='gastronoom') uses the input template to reproduce the GASTRoNOoM settings.
>>> eb = EB.EnergyBalance(gamma='h2') calculates the adiabatic coefficient from the temperature-dependent measured coefficient of molecular hydrogen.
>>> eb = EB.EnergyBalance(heatmode='gs2014') calculates the heating and cooling terms from dust-gas interaction based on the formalism of Gail & Sedlmayr 2014 as opposed to the classical implementations of Decin et al. 2006 and Schoïer et al. 2001.
Furthermore, the iterative procedure can be adapted to user-specific needs: >>> eb.iterT(conv=0.005,imax=200,step_size=0.05) sets the relative convergence criterion to 0.005 from default 0.01, the maximum iterations to 200, and the step_size of the gradual maximum allowed temperature change between iterations to 5%.
|
|||
|
|||
| list[str,dict] |
|
||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
| array |
|
||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
Inherited from |
|||
|
|||
|
Inherited from |
|||
|
|||
Creating an EnergyBalance object used to calculate the energy balance. Can work in conjunction with iteration with a RT code. Takes into account several contributors to the heating and cooling of an outflow:
Not yet implemented:
An instance of EnergyBalance contains three types of properties:
First set is given in aux/inputEnergyBalance.dat. Second set includes: radius, grain size, gas velocity, gas/dust mass-loss rates, dust opacities, dust temperature, stellar radiance. Third set includes: gas temperature, gas/dust density, drift velocity. Note that for now line cooling parameters are fixed. Hence after a new RT model was calculated, a new energy balance must be calculated. In the future, the level pops will be updateable.
|
Format an input keyword list for a variable. Can be given as a string, split by spaces, or as the list output from str.split(). Assumes the first element is the requested function or a constant, followed by key=value pairs that serve as the input keyword arguments for the function if the first element is not a constant. The first element thus does not contain a '=' sign. The key=value pairs are formatted into a dictionary. Example str: 'read_opacity species=AMCDHSPREI index=1 order=5' Example list: ['read_opacity','species=AMCDHSPREI','index=1','order=5']
|
Read the input parameters for the energy balance. If a filename is not given, the one set upon initialisation is used. If none was given upon initialisation the default one from aux/ is read. Resets any variables that depend on independent variables.
|
Initialise the coordinate grids: radius (cm), grain size (cm), wavelength (cm). The grids can be given as lists, which are used as input for Gridding.makeGrid. If not given, they are taken from the inputfile.
|
Set the stellar properties. If any of these are None, they are taken from the inputEnergyBalance.dat file. Defines the Radiance profile based on these properties.
|
Initialise the independent profiles:opacity, mass-loss rate, initial temperature, dust temperature The grids can be given as a 1- or 2-item list, which are used as input for the respective Profiler classes. If not given, they are taken from the inputfile. Format of the input lists:
|
Set the adiabatic coefficient. Three options:
|
Set the velocity profile, based on the requested function. Done separately from other profiles, because this can depend on the temperature at the inner radius (independent), and on the adiabatic coefficient (for the sound velocity). This is done before the line cooling term is set.
|
Set the line cooling parameters. This includes reading the collision rates and level populations. This is done per molecule, and upon initialisation of the EnergyBalance. A distinction is made between the source of the spectroscopy/level pops:
|
Set an initial set of level populations. Based on the Boltzmann distribution for a given excitation temperature. The kinetic temperature is taken for now. Should likely be scalable in the future. Solves the set of equations:
|
Set the abundance profile for a molecule. Two possibilities:
The second possibility requires a read method and a filename to be given in the molecule definition, e.g. molecule=12C16O np.loadtxt fname=waql.par usecols=[1,4] skiprows=9 unpack=1 Alternatively, a constant value can be given as well (not advised).
|
Calculate the excitation temperature for the level populations given in self.pop. This method is called when the pops are read for the first time (or set by setInitialPop), and when the level populations are updated. The excitation temperature is set as an array of dimensions (number of transition indices, number of impact parameter), hence for a given transition index, you can retrieve the excitation temperature as self.Texc[m][i-1,:] as a function of impact parameter. Alternatively, you can access Texc by calling getTexc and passing either the indices or the lup/llow. The impact parameter grid can be retrieved from self.pop[m].getP(). Note that the excitation temperature is calculated for all radiative transitions included in the molecular spectroscopy; not for all collisional transitions.
|
Return the excitation temperature profile as a function of impact parameter for given transition indices of a given molecule. Can also take upper and/or lower level indices to determine the transition indices. The excitation temperature is returned as an array of dimensions (number of transition indices, number of impact parameters).
|
Calculate the drift velocity profile. This assumes a stellar radiance profile is given, of which the wavelength grid is used to integrate Q_l * L_l, thereby interpolating the opacity. The two available modes are standard and beta:
Note that the vbeta mode follows the F mode for dust velocity in MCP. |
Set the dust velocity profile. The drift is averaged over grain size. |
Calculate the density profile for gas or dust. The dust density profile is not dependent on grain size, since the drift averaged over grain size.
|
Update the level populations. This reads the level populations from the same file as listed upon creation of the EnergyBalance and assumes those level populations are appropriate for the temperature profile of the iteration number self.i. This is done per molecule.
|
Set the current temperature profile for this iteration. When i == 0, this is the initial guess. This is done explicitly, because upon initialisation the state of the object depends on self.T being None or not. |
Iterate the temperature profile until convergence criterion is reached. Extra arguments are passed on to the spline1d interpolation of the total cooling and heating terms through the calcT() call.
|
Calculate the temperature profile based on the differential equation given by Goldreich & Scoville (1976). The function is defined in dTdr. The initial temperature is taken to be the second argument of the initial temperature profile given by T in inputEnergyBalance.dat. This is typically the condensation temperature. Additionals arguments are passed on to the spline1d interpolation of the total cooling and heating terms, e.g. k=3, ext=0 are defaults.
|
Calculate the adiabatic cooling rate. This is not used by the solver of the dTdr differential equation, and is only for plotting purposes. Uses the latest calculation of the T-profile. Equation derived from Decin et al. 2006 and Danilovich 2016. |
Calculate the heating rate by dust-gas collisions. Note that this term evaluates to zero for the inner wind at r<r0, if the dust density is 0 there. This can be enforced by choosing an mdot_step function. Sputtering is applied if a grain size distribution is used. For a constant grain size, this is not done, since it's not realistic to remove all dust if the drift becomes too large. The equation is derived from Goldreich & Scoville 1976, based on Schoier et al 2001, and Decin et al. 2006: General case: Hdg = n_d sigma_d w x 1/2 rho_g w^2 => Hdg = pi/2 n_H2 m_H (fH+2) (1+4fHe) (1-P)^(2/3) Int(n_d w^3 a^2 da) MCP/ALI case (average grain size a): Hdg = pi n_H2 m_H n_d w^3 a^2 GASTRoNOoM case (MRN): Hdg = pi/2 A n_H2^2 m_H (fH+2)^2 (1+4fHe) Int(w^3 a^-1.5 da) |
Calculate the heating rate by dust-gas heat exchange. The equation is derived from Burke & Hollenbach 1983, based on Groenewegen 1994, Schoier et al. 2001, and Decin et al. 2006: Note that this term evaluates to zero for the inner wind at r<r0, if the dust density is 0 there. This can be enforced by choosing an mdot_step function. Sputtering is applied if a grain size distribution is used. For a constant grain size, this is not done, since it's not realistic to remove all dust if the drift becomes too large. The older implementations used by MCP/ALI and GASTRoNOoM follow Burke & Hollenbach 1983 and Groenwegen 1994. The general-case implementation follows the more recent work of Gail & Sedlmayr, adding in proper vT and drift terms. The GS2014 is not yet fully implemented, but the Hdt term is already available here. General case (following Gail & Sedlmayr 2014, see Eq 15.19): Hdt = alpha pi k_b n_h2 (fH+2.)(1.-P)^(-2./3.) Int(n_d a^2 da) (T_d - T) sqrt(8*vT^2+drift^2) (1/(gamma-1)) MCP/ALI case (n_d: dust number density, average grain size a): Hdt = 2 alpha pi k_b n_h2 (n_d a^2) (T_d - T) vT GASTRoNOoM case (n_d: distribution following MRN): Hdt = 2 alpha pi k_b n_h2 (fH+2.) Int(n_d a^2 da) (T_d - T) vT In all cases, alpha is the accommodation coefficient (from Groenewegen 1994): alpha = 0.35 exp(-sqrt(0.002*(T_d + T)))+0.1 and vT is the thermal velocity: vT = sqrt(8 k_b t/(mu pi m_H)) |
Calculate the heating rate by the photoelectric effect. Two methods are available at present: 1) Following Draine 1978 and Huggins et al. 1988, as used by MCP (see also Crosas & Menten 1997) 2) Following Bakes & Tielens 1994., as implemented by Decin et al. 2006 in GASTRoNOoM. No options to tweak these methods has been implemented yet, but a lot of consistency checks should be done, and some of the assumptions can be improved on with a consistent calculation (e.g. Av in method 2). Photoelectric heating as derived by Bakes & Tielens is maybe a good basis for further development of the heating term, but see also Woitke 2015 (2015EPJWC.10200011W) for recent developments. The derivation in method 2 makes a lot of assumptions and are applied specifically to the case of GASTRoNOoM (e.g. fixed grain size distribution). Hence, the implementation is considered appropriate for the default settings in the inputEnergyBalance_gastronoom.dat file. Any deviations from that must be treated carefully. No sputtering is applied. Small grains are the most relevant for Hpe, and any sputtering would reduce the size of existing dust grains, thus increasing the amount of small grains. Since we cannot have Hpe directly depend on the grain size distribution (and instead work with a scaling factor amin_scale), we do not apply sputtering to the photoelectric heating. |
Calculate the heating rate by cosmic rays. The rate is taken from Goldsmith & Langer 1978, and is also used by Groenwegen 1994 and Decin et al 2006. This is a very approximate formula and will need updating in the future. Two modes are available: The standard one, using n(H2), and the one used by Groenewegen 1994 and Decin et al. 2006 (where all gas mass is placed into H2, even if fH or fHe are non-zero). Controlled through the keyword cr_method == 'groenewegen' or cr_method == 'standard' in the inputfile. |
Calculate the line cooling rate by vibrational excitation of H_2. Two modes are available (used by MCP/ALI and GASTRoNOoM, respectively): 1) Following Groenewegen 1994, based on Hartquist et al 1980. Based on fitting of tabulated data of H2 cooling under LTE conditions 2) Following Decin et al 2006, based on GS1976, Hollenbach & McKee 1979 and Hollenbach & McKee 1989. The keyword h2_method (groenewegen, or decin) determines which of the two is used. |
Calculate the radiative line cooling rate for a molecule. For GS2014, no correction term yet for the excitation temperature per level. Follows Sahai 1990. Cubic spline interpolation for the level populations. Linear interpolation and extrapolation for the collision rates (as for GASTRoNOoM). Currently not yet implemented to use a sqrt(T/T0) extrapolation at lower boundary as is done by MCP/ALI. Note that EnergyBalance will internally call this function with the molecule tag tacked onto the Clc function name.
|
Plot the heating or cooling rates for several iterations of a given mechanism and heat exchange sign (cooling or heating).
|
Plot the heating and cooling rates for a single iteration.
|
Plot the temperature profiles of the different iterations.
|
Plot the total heating and cooling term (excluding adiabatic cooling) calculated by calcT before the differential equation is solved. This includes the density factor that enters, and so is essentially (H - C)/rho. The velocity does not enter here, since that is calculated explicitly in dTdr. Hence the y-axis units K/s.
|
Plot the excitation temperature as a function of impact parameter for a selection of radiative transitions, given by either the transition indices or the upper and/or lower level indices. Retrieves the excitation temperature from self.Texc. They are appropriate for the level populations used by the current iteration, and are updated when the level populations are updated. At this time, plotting Texc for other iterations is not possible.
|
| Home | Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0.1 on Mon Nov 7 18:01:58 2016 | http://epydoc.sourceforge.net |