J/A+A/633/A25 meteor_spect.f90 code package (Park, 2020)
Radiation phenomenon for large meteoroids.
Park C.
<Astron. Astrophys. 633, A25 (2020)>
=2020A&A...633A..25P 2020A&A...633A..25P (SIMBAD/NED BibCode)
ADC_Keywords: Models
Keywords: meteorites, meteors, meteoroids - radiation mechanisms: thermal -
shock waves
Abstract:
The relationship between the intrinsic properties of a large
meteoroid, that is, mass, density, chemical composition, and flight
velocity, and its spectrum are identified.
Assuming thermochemical equilibrium and inviscid flow and using a
quasi-onedimensional approach, the flow-field and radiative transfer
problems are solved from the ablating wall of a meteoroid to the
observer on the ground.
The study discovers that the precursor region immediately ahead of the
shock layer absorbs nitrogen and oxygen lines, and thereby modifies
the spectrum received by the ground observer. The study leads to a
computer code with which the observed spectra and light-curves for
Benesov and Sumava meteoroids are numerically and approximately
reproduced with the help of a meteoroid fragmentation theory. Luminous
efficiency value is obtained as a byproduct.
Description:
README for meteor_spect.f90 code package
A. Overview
In the meteor_spect.f90 code, first the absorption coefficients are
calculated in the three regions: in the air region, in the
ablation-product vapor region, and in the precursor region. With nwave
set to 400000, the absorption coefficient calculation requires 20
minutes on a personal computer. This calculation does not involve any
iteration, and therefore never fails. The calculated absorption
coefficients are recorded on three different tables, for air shock
layer, ablation-product vapor layer and the precursor region.
Next, the flow solution and radiative transfer calculation is carried
out for the given condition. This calculation uses the absorption
coefficient tables created above. The calculation involves iteration:
The calculated radiative transfer rate determines the local energy
input (divergence of radiative flux is the right-hand side of the
energy conservation equation). When this local energy input is added
to the local enthalpy value, the local enthalpy value changes. This in
turn changes the radiative flux values. In order to prevent from
overshooting, the corrections on enthalpy changes are made only a
small fraction of the calculated value. This action occurs separately
for air shock layer and the ablation-product layer. Therefore there
are two sets of these fractional change values to be specified in the
input file. These allowed fractional changes are: about 0.5 for the
air shock layer and about 0.05 or smaller for the ablation-product
vapor layer. If these fractional changes are too large, the solution
diverges.
The time required by this calculation varies with the number of
iteration specified in the input file. With nwave set to 400000, and
the minimum acceptable number of iterations, the radiative transfer
calculation requires about 40 minutes on a personal computer. If the
solution diverges, you restart with different values of allowed
fractional changes. But in this case, you can use the absorption
coefficient tables already created: there is no need to repeat that
calculation.
The meteor_spect.f90 code package consists of:
(1) meteor spect[2].f90 Fortran f90 source code. 17201 lines,
(2) atom.dat input file 51667 lines, giving atomic spectroscopic data,
(3) diatom.dat input file 5736 lines, giving diatomic spectroscopic data,
(4) triatom.dat input file 1284 lines, giving triatomic spectroscopic data,
(5) eqtabairvis.dat input file 3837 lines, giving equilibrium air data,
(6) meteor_spect.inp input file 10717 lines, giving instruction on run.
A successful run produces the following files:
(1) atom.out file, repeating the atom.dat file, to verify that the
atom.dat file is correctly read in. After a successful run (which means
that atom.dat file is correctly read in), this file can be discarded.
(2) diatom.out file, repeating the diatom.dat file, to verify that the
diatom.dat file is correctly read in. After a successful run (which means
that diatom.dat file is correctly read in), this file can be discarded.
(3) triatom.out file, repeating the triatom.dat file, to verify that the
triatom.dat file is correctly read in. After a successful run (which means
that triatom.dat file is correctly read in), this file can be discarded.
(4) absb_air.out file, giving a table of absorption coefficients for air
for use in the air shock layer. Keep this file: you may need this file
in the next run.
(5) absb.vap.out file, giving a table of absorption coefficients for
ablation-product vapor for use in the ablation-product vapor layer.
Keep this file: you may need this file in the next run.
(6) absb.low.out file, giving a table of absorption coefficients for air at
low temperatures for use in the precursor region. Keep this file:
you may need this file in the next run.
(7) amdot.tally file, giving a short summary of the fluid-mechanical
aspects of output.
(8) flowfield.out file, giving distribution within the air shock layer and
the ablation-product layer of temperature, mass flow rate, density,
radiative heat fluxes in two (outward and inward) directions.
(9) inner_conv.out file, giving the convergence history of the solution in
the ablation-product layer.
(10) outer_conv.out file, giving the convergence history of the solution in
the air shock layer.
(11) meteor_spect.out file, giving the main results of the computation.
(12) short_file giving a coarse table of spectrum at four locations, i.e.
at the interface, at the shock wave, at the edge of precursor region,
and ground.
(13) spect_air.plt file, giving the spectrum of the air shock layer.
(14) spect_vap.plt file, giving the spectrum of the ablation-product
vapor file.
(15) temp_display file, giving the progression of temperature distribution
within the ablation-product vapor layer. After a successful run,
this file can be discarded.
(16) trouble_shooting file, giving what happens at the portion of the code
calculating the distribution of temperature in the ablation-product vapor
layer. After a successful run, this file can be discarded.
B. meteor_spect.f90 source code
The source code is assembled as a single Fortran f90 code. Subroutines
are arranged in alphabetical order. All quantities are in double
precision.
There are two parts in this source code: (1) absorption coefficient
calculation, and (2) flow-field and radiative transfer calculation.
Absorption coefficients are calculated at five (for air) or eleven
(for ablation-product vapor) temperatures. In order to calculate
absorption coefficients, first one must know the number densities of
all component species. The number densities are calculated from the
partition functions assuming equilibrium, using the rate-equation
method. For this purpose, energy levels and their multiplicity values
must be known. These are given in the input files atom.dat,
diatom.dat, and triatom.dat files. There are three parts here: air
shock layer, ablation-product vapor layer, and the precursor region.
After the equilibrium condition calculation, absorption coefficient
calculation is carried out in the subroutine radipac which calls
subroutine emis_absb. The intensity parameters, the transition
probability for an atom, electronic transition moment for diatom, or
absorption cross section for triatom, are in atom.dat, diatom.dat, and
triatom.dat files.
The calculation of flow-field and radiative transfer is carried out in
subroutine comet_lbl. comet_lbl calls subroutine calc and subroutine
calc1, and actual calculation is made in these two subroutines. calc
and calc1 are nearly the same. In these two subroutines, radiative
transfer and flow-field, i.e. enthalpy field and density-velocity
field, calculation are carried out. Radiative transfer calculation is
made in subroutines flux_outer and flux_inner. flux_outer is for the
air shock layer and flux_inner is for the ablation-product vapor
layer. The radiative transfer calculation is made in two directions,
outward and inward. The divergence of radiative heat flux is the
source in the energy equation. Energy equation is solved in
subroutines calc and calc1. The density-velocity product, which is the
solution of mass and momentum equations, is made in subroutine ffp5.
Both these solutions, enthalpy and density-velocity, are made partly
numerically and partly analytically. This is because there is a zero
singularity at the interface. Numerical integration is made to 80% of
the way. The solution so obtained is fitted with third-order
polynomial. The polynomial is extended to the singular point. These
calculations produce ablation rate and radiation spectrum at
everywhere.
Emissivity of the wall is given as an input. The value of 0.5 is
chosen in most of the calculations. This value is a good value because
the surface is melting and therefore smooth. A smooth surface of the
melting rock should have emissivity of about 0.5. This also means the
surface has a reflectance of 0.5. The surface, therefore, should
radiate with intensity equal to black-body value times emissivity plus
the reflection of the incidence radiation times reflectance. This is
how the wall value of radiation intensity is determined.
Inside subroutine ffp5, the solution to mass and momentum conservation
equation is obtained through a fifth-order polynomial. Again, the
governing equation, a nonlinear second-order ordinary differential
equation has a zero singularity at the interface. Several different
order of polynomial was tried. Fifth order was found to be the best.
Higher order polynomial solution was found to diverge.
C. atom.dat, diatom.dat files
atom.dat file gives spectroscopic parameters of atoms. The file is
self-explanatory, except the starkg, starkn, ngi, and ngk. Stark
half-width at half-height is assumed to obey the expression w =
starkg(Ne/10e16)**starkn in A, Ne is electron density in cm-3. These
starkg and starkn are given in atom.dat file in columns 97 to 112,
starkg from column 97 to 106 and starkn from 107 to 112. When they are
left unwritten or put in the values of 0, the code internally
determines starkg and starkn by the default formula. starkg =
0.042(???2/(Ip - E))2 A, and starkn = 033. Where wavel is wavelength
in A, Ip is ionization potential in cm-1, and E is the upper state
energy level in cm-1. The quantities ngi and ngk are used for
calculating non-Boltzmann distribution of states, and are not used in
the present work. The Stark width values in Griem's Plasma
Spectroscopy are in the file for N and O. Other than that, 77 selected
lines in total are given a finite value with the intention of
numerically reproducing the spectrum for Benesov meteoroid.
There are two different ways in specifying the bound-free continuum in
atom.dat, by cross sections or by Gaunt factors. The twenty-seventh
column after the species in the bound-free continuum section, there
will be either character c (for cross section) or G (for Gaunt
factor). By reading this character, the code understands whether the
data give cross sections or Gaunt factors. The cross-section values
are entirely from Peach's table, G. Peach, Memorandum of Royal
Astronomical Soceity, Vol 73,(1970) 1-123. Gaunt factors are
calculated assuming that the net calculated Gaunt factors will be
unity.
The diatom.dat file will also be self-explanatory. The electronic
transition moment values include Franck-Condon effects. They go
directly into the calculation of intensity of diatomic bands. The
electronic transition moment values are obtained from the published
data. In order to put the published data into diatom.dat, one needs to
carry out computation involving band origin, r-centroid, and
Franck-Condon factors. First, one must produce these parameters using
a Shroedinger equation solver. This task was carried out to produce
the values in diatom.dat. In the future, anyone wishing to expand the
diatom.dat file must repeat the same procedure.
Column 43 to 46 in the line after the species line gives four
characters that specify splitting of branches. The present code
distinguishes all branches in doublets and triplets faithfully
following Herzberg's book and the book by Kovacs, I., "Rotational
Structure in the Spectra of Diatomic Molecules," Adam Hilger, 1969.
The details are given by Hyun Hyun, S. Y., "Radiation Code
SPRADIAN07 and Its Applications," Ph.D Thesis, Department of
Aerospace Engineering, Korea Advanced Institute of Science and
Technology, 2009.
D. meteor_spect.inp file (the main input file)
The parameters appearing in meteor_spect.inp file are explained in the
file. There are two most important parts in this file: the first 23
lines and from line 1219 to line 1235.
nwave in line 2 is the number of wavelength points. For an accurate
run, it should be set to 400000. But smaller values are allowed,
though sacrificing accuracy.
nair in line 3 and nabl in line 4 can be either 0 or 1. When set to 0,
calculation is made of the absorption coefficient tables. When set to
1, the absorption coefficient calculation is skipped, and use is made
of the absorption coefficient tables made in the earlier run. Of
course, these tables must be available to read in.
rhoinf, vinf, rnose in lines 8, 9, 10 are as indicated in the file.
ndm in line 11 is the number of node points in the air + vapor layer.
nout in line 12 is the number of iterations on the air layer at one
time. There are four times of this nout iterations.
nim in line 13 is the number of iterations on the vapor layer at one
time. There are four times of this nim iterations plus 10 times.
facout1 in line 14 is the fraction by which the old enthalpy profile
is modified by the new enthalpy profile for air shock layer in the
first two sweeps of iteration.
facin1 in line 15 is the fraction by which the old enthalpy profile is
modified by the new enthalpy profile for vapor layer in the first two
sweeps of iteration.
facout2 in line 16 is the fraction by which the old enthalpy profile
is modified by the new enthalpy profile for air shock layer in the
second two sweeps of iteration.
facin2 in line 17 is the fraction by which the old enthalpy profile is
modified by the new enthalpy profile for vapor layer in the second two
sweeps of iteration.
famdot in line 18 is the fraction by which the old ablation rate is
modified by the new ablation rate during iteration.
enfac in line 19 is (shock layer enthalpy)/(flow enthalpy) accounting
for radiative cooling. User must choose this parameter so that, in the
last statement of the output file, engy_sum2 equals 0.5*vinf*vinf
(freestream enthalpy). By choosing enfac appropriately, Eq. (32) in
"Radiation phenomena in large meteoroids" is satisfied. enfac is
the ratio of the second term on the right-hand side of Eq. (32) to the
left-hand side.
elev_angle in line 20 is the elevation angle of the meteoroid when
seen by the ground observer.
hwhh in line 21 is the half-width at half-height of the Gaussian
profile that describes the slit function of the spectrograph used by
the ground observer.
sigma in line 22 is the emissivity of the ablating wall of the
meteoroid.
spallf in line 23 is the spallation factor, which is the ratio of the
effective ablation rate to the thermo-chemically determined ablation
rate. The effective ablation rate is larger than the thermo-chemically
determined ablation rate because of particle emission by spallation.
From line 24 to line 1218, spectroscopic data on air species are
given. These data must not be changed.
In line 1221, number of elements and species in the vapor layer are
given. These numbers must not be changed.
From line 1222 to line 1235, mass fractions of 14 chemical components
of the meteoroids are given. These mass fractions must be chosen by
the user. "Radiation phenomena for large meteoroids" gives several
important sets of these mass fractions.
From line 1236 to line 10714, spectroscopic data on vapor are given.
These data must not be changed.
E. meteor_spect.out file (the main output file)
The output file starts copying the input file. Then it prints out the
progress of the calculation. Only toward the end are the output
parameters.
Toward the end of the file, you find:
frozen inviscid shock standoff distance, chemical relaxation distance,
viscous mixing/cooling distance, and equilibrium inviscid shock
standoff distance. These quantities are some of important output
parameters. Chemical relaxation distance must be much smaller than
equilibrium inviscid shock standoff distance, or better still, the
thickness of the air shock layer or vapor layer found in the
flowfield.out file, in order for the equilibrium flow assumption to be
satisfied. Equilibrium flow assumption is the fundamental assumption
made in 'Radiation phenomenon in large meteoroids'. Also, viscous
mixing/cooling distance must be much smaller than the thickness of the
air shock layer or vapor layer in order for the assumption of inviscid
flow to be valid. This assumption is also one fundamental assumption
made in 'Radiation phenomenon in large meteoroids'.
They are followed by qps(radiative flux at shock wave), and a table
giving spectrum named 'normal intensity entering precursor
region'. In that table, iw is the index of wavelength values, int_s
is the intensity of the normal ray leaving the shock wave in the units
of W/(cm2-µm-sr), and cros_sec is the absorption cross sections of
air in the free-stream.
They are follows by 'enthalpy in precursor region raised by
absorption', and a table giving igrid, r(m), q(J/m3), enth(J/kg).
igrid is the index of positions, r(m) is the distance from the shock
wave in meters, q(j/m3), is the divergence of radiative heat flux in
J/kg, enth(J/kg) is enthalpy in J/kg.
They are followed by 'flux in precursor region above 1750 A'. The
variables, igrid r(m), temp, flux(W/m2) are now self-explanatory.
Next, the brightness magnitude is printed out. This is, of course, for
this one meteoroid only.
Next, Boltzmann plot is carried out. Four Mg I lines and twenty-seven
Fe I lines are tried. First, four Mg I lines are plotted. E is the
energy level in cm-1, wavel is the wavelength in A, y is I/(gA), and
log(y) is the natural logarithm of y. Then a least-square
straight-line fit is made on the four points. bx(1,1) and bx(2,1) are
the coefficients of the straight line. Next, the straight line
coordinates are given. Then the Boltzmann temperature obtained from
the straight line is given.
Next, the Boltzmann plot with 27 Fe I lines is given. For each data
point, the details of the readings are given. Because the intensity
value is an integration under the curve, one must know the starting
point wavelength and the intensity value, the peak point wavelength
and intensity value, and the last point wavelength and intensity value
are given. In this way, the quality of the reading of each line is
assessed. wavelc is the tabulated wavelength of the line. ier is an
error message. If it is 1, there is no error. After going through 27
lines in this way, the Boltzmann plot is made. area is the area under
the intensity curve; y is I/(gA), and log(y) is natural logarithm of
y. Least-square straight-line fit is made with these 27 points, and
the temperature so determined is printed out. Frequently, one of these
lines results in a negative area. When this happens, temperature is
naturally indeterminate.
Next comes the energy tally. flux_s is the radiative heat flux at
shock, flux_prec is the radiative heat flux at the edge of the
precursor region, and flux_grd is the radiative heat flux at the
ground. The most important quantity here is the engy_sum2. It is enths
+ flux_prec/(rhoinf*vinf). This quantity must equal 0.5*vinf*vinf
(freestream enthalpy), to satisfy Eq. (32) in 'Radiation phenomena
in large meteoroids'. The error in accomplishing this is given as
error fraction = engy_sum2/(0.5*vinf*vinf) - 1.0). The absolute
magnitude of this error should be below 0.01 (1%). In order to bring
error to such a small value, the enfac in input file (19th line in
meteor_spect.inp) must be adjusted. Brightness magnitude and luminous
efficiency values are given here.
Next Borovicha parameters are printed out. These parameters are needed
in the equations of meteor defined by Borovicka in "Radiation Study
of the Two Very Bright Terrestrial Bolides and an Application to the
Comet S-L 9 Collision with Jupiter", ICARUS, Vol 121, 1996, pp.
484-510. F. spect_air.plt, spect_vap.plt, and short_file files These
two files give the calculated spectra. In spect_air.plt file, int_s is
the intensity of the normal ray at the shock wave, int_prec is the
intensity of the normal ray at the edge of the precursor region, and
int_grd is the intensity of the normal ray at ground, all in the units
of W/(cm2-µm-sr). camera is the intensity of the normal ray captured
by the spectrograph. Unlike other intensity values, this intensity has
been processed by the slit function. The slit function is assumed to
be Gaussian. The half-width at half-height of the Gaussian slit
function is given by the line 21 in the input meteor_spect.inp file.
The spect_vap.plt file gives the intensity of normal ray given by the
ablation-product vapor. int_e is the intensity of normal ray at the
edge of the vapor layer. wall(+) is the intensity of the normal ray
emitted by the wall. This intensity is the sum of the black-body
radiation times emissivity + partial reflection by the wall.
short-file gives the intensity of (outward) normal ray at four places,
edge of vapor layer, shock wave, edge of precursor region, and ground
roughly. That is, at 1000 wavelength points: 1000 wavelength points
are selected out of the 400000 wavelength points. This file is
provided to give a quick overview of the radiation intensity.
G. absb_air.out, absb_vap.out and absb_low.out files
These files give the tables of absorption coefficients at several
temperatures. absb_air.out is for the air layer; absb_vap.out is for
the vapor layer, and absb_low.out file is for the precursor region.
For absb_air.out, absorption coefficients are given at five
temperatures; absb_vap.out and absb_low.out, absorption coefficients
are given at eleven temperatures. These coefficient values are for a
given altitude and given flight velocity. In a fresh case, these
absorption coefficient values are calculated and stored in tables
first, by setting nair and nabl in lines 3 and 4 in meteor_spect.inp
file as 0. Next time absorption coefficients are needed for the same
altitude and same flight velocity, these table files can be read in by
setting nair and nabl both as 1.
H. flowfield.out file
flowfield.out file presents the variation in various thermodynamic and
fluid mechanical variables in the air shock layer and the vapor layer.
y is the distance from wall in meters; rho is the density in kg/m3;
rhov is density-normal velocity product in kg/(m2-s); temp is
temperature in K; qp is the radiative heat flux in the plus (outward)
direction in W/m2; qm is the radiative heat flux in the minus (inward)
direction in W/m2; dqdy is the divergence of radiative heat flux in
W/m3; and enth is enthalpy in J/kg.
I. inner_conv.out and outer_conv.out files
These two files show the convergence pattern of the solutions.
inner_conv.out is for the vapor layer and outer_conv.out is for the
air layer.
J. amdot.tally file
amdot.tally file lists twelve output parameters. amdot is the ablation
rate, which will be needed in determining luminous efficiency. delH is
the ablation energy. fw is the so-called blowing parameter. If it is
much larger than unity, then the flow-field over the meteoroid is said
to be massively blowing. In a massively blowing situation, convective
heat transfer rate is negligible. The list includes shock stand-off
distance and mixing thickness dmix. If dmix is much smaller than the
shock stand-off distance, the flow-field can be considered inviscid.
The list also gives four important temperatures, wall temperature,
temperature at the edge of ablation-product vapor layer, interface
temperature which is the temperature of the air layer at the
interface, and the temperature immediately behind the shock wave.
K. temp_display and trouble_shooting
temp_display file displays the distribution of temperature in the
ablation-product vapor layer at each step of iteration. One can judge
how iteration is progressing. trouble_shooting file gives temperature,
enthalpy, and mass flow rate at the time of starting the iteration on
the ablation-product layer. These variables can tell whether the
choice of the correction factors in the iteration for the vapor layer
facin1, facin2, and famdot.
File Summary:
--------------------------------------------------------------------------------
FileName Lrecl Records Explanations
--------------------------------------------------------------------------------
ReadMe 80 . This file
diatom.dat 289 5736 diatom.dat file
meteor_spectX2X.f90 181 17201 Fortran code
triatom.dat 186 1307 triatom.dat file
absb_air.out1 60 200002 Example for absorption coefficients
absb_air.out2 60 200001 Example for absorption coefficients
--------------------------------------------------------------------------------
Acknowledgements:
Chul Park, cpark216(at)yahoo.com
(End) Patricia Vannier [CDS] 04-Dec-2019