Radioactive Fallout Model (Radfo) - Part II

by Buddy McCloskey     11/28/2008     Contact

1. Purpose:

The goals are mentioned in the first few paragraphs of World Trade Center Dust Bowl Analysis. We extended the results and theories from Radioactive Fallout Model (Radfo) - Part I and the previous studies (Demolition Dust Model and World Trade Center Dust Bowl Analysis ) to apply here with the building of components for the Radfo Model to become operational in real time. This involves the dust inventory that was the main core of Radfo Part I and the meteorology impact on these new modular components.

2. Bomb Yield Parameters:

The pic file Graph Plots shows all the plots on one graph. It is noted that they are almost perfect straight line in a log-log plot. So we use b=(ln y2/y1)/(ln x2/x1) and a=y/x^b as before. There are 2 distinct plots for the cap cloud base. Usually, two points that are on the line are chosen to calculate the line. For flash time the x=1 sec data point defines the a=2.5 kiloton in y=ax^b.

Origin Graph-Table         Parameter                 Formula             Color in Graph Plots
APT25a     Cap Radius in Statute Miles         10.86*mega^.4197         Green
APT25a     Cap Top in 1000's '                       70.38*mega^.1797         Bronze
APT25a     Cap Base in 1000's ' Mega > 20     46.1*mega^.1808         Blue
APT25a     Cap Base in 1000's ' Mega <= 20   43.75*mega^.1601       Blue
SRFM1     % Growth Cloud Top                     32*min^.6195               Purple
APT25b     Kilo from Flash Time Sec             2.5 sec^2.0323             Red or
APT25b     Flash Time Sec from Kilo             (kilo/2.5)^-2.0323         Red

3. Fireball Plume Rise Equations:

Adjustments for the cap top & bottom heights are needed. Looking at variety of graphs depicting cloud cap top vs. yield, there is a margin of error +,- 25,000'. Graph by FEMA shows an example of this. The cap & stem relative dimensions and % particle distribution is seen in pic SRFM2A.gif . Using 4 megatons in this example, the cap base goes from 54,000' to 90,000' about a mean of 72,000'. This range becomes larger, with larger bombs. At 1 megaton, it's about a 31,000' difference in range or +22%, -19%% error. FEMA's accompanying article states the atmospheric conditions and the winds are the major reasons for this variability. Let's examine these factors.

       A. Stability:

Since cap top will usually be above the mixing depth boundary layer (say 8000' max), the upper air lapse rates will be > -1 deg C/100m. Lets use a range of -.5 -> 1 for this upper atmosphere and this will be an input parameter with most likely an 'E' stability class in most all cases. However, due to these greater heights, the available temperature profiles will dictate the stability parameters more accurately.

       B. Air Temperature:

This determines the magnitude of the buoyancy and the kinetic momentum term and should also be considered as we will have temperature profile as an input parameter. Air temperature is also needed for density calculations for the particle terminal fall velocity through each layer of the atmosphere.

       C. Wind Vectors:

Wind velocity is likely the main important atmospheric variable for final plume heights. We also need wind vectors to calculate particle cylinder translation and shear in a downwind direction.

       D. Stack Plume Rise Formulas:

Many authors have developed formulas for plume rise of stacks. What is given is a variety of stack conditions and observed plume photos vs wind velocities and stabilities. With all the data needed available near the surface, the exponent wind profile law can approximate the wind velocity at stack heights. Many papers have been published and it's possible we have now reached a point of no return in prediction improvement. This is probably due to mainly the unstable case where thermals below the boundary layer mixing depth cause looping of the plumes and the accompanying wind gustiness causes fast variation of the plume rise. This thermal reaches the mixing depth which could be over 6000' in summer. The dimension of a thermal diameter near 1/2 mile can have quite an influence on the plume rise for 5-10 min as the plume is rising with the thermal to it's final plume rise height and then visa versa with the downdraft. Also the intense super adiabatic variations are so extreme that this plays havoc with the buoyancy term. Recall, that class A stability is open ended ( <- 1.8 deg C/100m). The stacks can be considered steady state eddy (in one direction) within a thermal but with a nuclear blast, this air thermal is considered as one big turbulence eddy within the stem of a nuclear blast. Also, the nuclear blast thermal rises to 10-20 times the height of the boundary layer mixing height. Since this is a region of rarely D but mostly E & F class stabilities, the final plume rise should be more consistent in this steady state stable environment and easier to predict. As for nuclear blasts, we have a much smaller sampling of accurate met data available at the time of fewer test blasts than stack tests.

            (1) Basic Stack Plume Rise Formulas:

The best source of plume rise data is taken from TVA. You would think that these large stack diameters, high exhaust temperatures and large stack velocities would better approximate a nuclear blast. To best analyze this, most basic plume formulas are usually of the form dh = aX^b F^c/U^d. Now, a is a constant that varies with atmospheric stability, X is the downwind distance (m) and needed where the plume levels (final rise when dh=0), and F is the buoyancy & momentum flux terms, either/or depending on the greater or both, and U is the mean wind velocity at stack height. The exponents b and c are positive but < 1 and we have d as negative. Usually X is 2000m or less and b=.2 -> .6 again may be sometimes linked to stability, and c is usually known as the 1/3 power law. However, for U, d could be -1 -> -1/2, which depends on the author of the selected formula. U can also be combined with other stability factors. Higher exponents of U may also be used for entrainment.

            (2) Stack or Thunderstorm?:

The problem we also have is what to considered as the physical stack height for a fireball? Most bombs producing their maximum effect from a blast would be 2000'- 5000' above ground level and is usually considered this level to produce the particle distribution. Actually, a nuclear blast could be considered a large singled cell thunderstorm. However, since the synoptic situation is such that multiple cells won't develop like regular thunderstorms do to develop a squall line. Cloud tops in thunderstorms can be calculated by comparing positive and negative energy areas, but can this be applied here? How big is the positive buoyancy area? We will consider it a surface burst to simplify our dh calculations. Our idea is to use the formulas derived from APT25a & the pic Graph Plots above for cap cloud height and scale that height in proportion by ratio of the calculated plume height rather than predict the absolute height. So, we first try to convert the nuclear blast to a super size stack to model in hope to have at least a best estimate for scaling.

The first problem is to determine a fictitious stack height (Hs). Comparing this to a real stack, the walls of a stack inhibits entrainment of the outside air by no wind blowing horizontally through the stack into the vertical escape jet. In a stack, we have continuous heat and momentum fluxes. In our phony nuke stack, we have only buoyancy at first as the fireball initiates the process. Once the buoyancy starts, we have a momentum term with the inrush of adjacent ambient air. However, during the flash time, all this inrush air is also super heated, just as the initial air volume is, in and surrounding the fireball. However, later we have strong surface winds being sucked into the fireball as seen by pic Fig. 1.3-1 .

            (3) Review of Nuclear Cloud Growth:

Based on the bomb size, we do have some empirical parameters generated. We have the dimensions of the cloud and stem. That is, we take 1/3 of the cap diameter as seen in pic SRFM2A.gif again, to indicate that the stem feeds the cloud cap top. This area of violent updraft best represents the stack but how high? Well the 'pg' formula represents the % growth of the cloud top, independent of bomb size, as depicted in SRFM1 again. Notice, the slowing down of cap cloud height with time. We might take the near surface initial growth to approximate the fastest cloud growth. Since the cloud height and size stabilize in 9 minutes, then let's take the first minute which is 1/3 the final cap height. However, notice the flash time in Mega . For 10 megatons, the flash time is near 1 minute, for 1 mega, it's 19 sec. Therefore, if we choose the 1st minute, the steepest slope on the graph for the fastest cloud rise, the fireball has stopped heating the gases for 41 seconds so the violent updrafts should be less by that time so the smaller average is not representative.

            (4) Constant Updrafts?:

To normalize all bomb blasts, let's consider only the time that the fireball exists. That is, the illumination time, or the older term of smaller bombs - flash time. This best simulates a stack as we consider we have a near constant stem velocity, as long as the fireball exists, as well as constant stem temperature. That is: Hs=pg/100 Cap*1000' and Vs=Hs/flash where Hs is a fictitious stack height in ft. and flash is the flash time in sec. Then Vs is stack (stem) velocity in ft/sec. We are assuming that the growing cap top is caused only by the vertical motion of the hot gases in the stem. Measured temperature from test blasts put the stem near 3000 - 4000 deg F. For our 1st estimate calculations, we use 3000F. The estimated ambient air temp at Hs for 1 megaton, say 32F and a wind velocity of 18 knots for Hs=11,000'. If we considered the final predicted cap height of 70,380' by the bronze line in pic Graph Plots , for 1 meg, pg = 15.7% x 70,380'= 11,050' after 19.1 sec flash time. Therefore, the nuclear cloud is probably mostly stem at that point before the cap develops. So, updraft Vs = 11,050'/19.1 sec = 582'/sec = 177 mps = 396 mph. Actually, this procedure was repeated for .1,.5,1,2,3 and 5 megatons. They all produced a Vs that ranged from 180 mps @ .1 mega -> 175 mps @ 5 mega so the 2 independent graphs for % growth and cap top for bomb size are in excellent agreement @ 1.5% error. This must mean that the gas expansion from the hot fireball is constant, but just a bigger stems with larger bombs.

            (5) Testing Stack Plume Rise Formulas:

Now 3 Delta Height (dh) formulas are based on TVA's power plant stacks. One produced plume heights over a factor of 5 times the empirical cap top y=ax^b formula. This appears so wildly high but 2 out of 3 equations were approximately the same. Another formula was tested and it appears to fall within the graph by FEMA range. This was such a surprise, that it might be best to work with. It appears that the other 2 equations were too sensitive to diameter for this application.

This equation used is explained below, but the important thing is, why was it so much more accurate? It is by Dr. Gary Briggs, the main plume height author used by EPA. He must have adapted the plume rise to fit all types of stacks for standardizing the model for regulation purposes. Since he uses smaller stacks, the parameters that he introduce were entrainment parameters. To think that a large diameter behaves like a stack, shielded from the ambient air due to the mere size may be erroneous. His added feature of the entrainment parameters must have made that difference but was only 2/3 rd's of the observed cap top. Since it was closer and crossed unity to become under predicted, it is worth examining. The equation by Briggs is an older EPA model (1980's) but the main theory is there and will be adjusted by fudge factors anyhow. However, the most important feature of this equation is the somewhat separate momentum and buoyancy terms for plume height impact within theory to be explained below. While working for the City of Phila., as their air pollution meteorologist, I tried that formula on a forest fire and I can only recall that it worked well as confirmed by aircraft. What could be another test case is an active volcano during it's steady state emission rate.

            (6) Gary Briggs Stack Plume Rise Parameters:

pi = 3.14159 Pi
At = absolute temp 273.15 [deg C]
Ga = 9.806 Acceleration of gravity [m sec-2]
U = Ambient air velocity at stack height Hs [m sec-1]
Tg = Temperature Gradient [deg C/100m]
PTg = Tg+.976 Potential Temp Gradient [deg C/100m]
Hs = stack height [m]
Ta = Temp of ambient air [deg K]
Ts = Temp of stack gas [deg K]
Ds = Stack Diameter (1/3 of the cap diameter) [m]
Vs = Stack Velocity [m sec-1]

Fb = Ga*Vs*(Ds/2)^2*(Ts-Ta)/Ts Buoyancy Flux Term [m^4 sec^-3]
Fm = Ta/Ts*(Vs*Ds/2)^2 Momentum Flux Term [m^3 sec-2]
s = Ga*PTg/100/Ta Stability parameter term [sec -2]
Es = U/Vs+.333 Stack jet entrainment parameter for momentum term
Ea = .6+(Tg+1.7)*(-.06364) air entrainment parameter for buoyancy term

B = 3*Fb/(Ea^2*U*s) Buoyancy Term
M = 3*Fm/(Es^2*U*s^1/2) Momentum Term
Xf = pi*Us^-1/2 Distance to final plume rise [m]
B3 = B*(1-cos(Xf*s^1/2/U) Final Buoyancy Term Cubed m^3
M3 = M*sin(Xf*s^1/2/U) Final Momentum Term Cubed m^3
dh = (B3+M3)^1/3 Final Plume Rise [m]
H = Hs+dh Final Plume Height [m]

            (7) Adaption to Nuclear Clouds:

The only terms that need to be adjusted for our large scale stack environment is the jet and ambient air entrainment parameters because of the relatively larger Vs than in power plant stacks and because of no protective barrier between the ambient air and the stem air to such high heights. The influx of the air into the cloud base also means that the entrainment is the source of the stem velocity already and is heated - probably not uniformly so Es seems to hold one of the keys here. Since we calculated a relatively high updraft velocity Vs (177 m/sec = 396 mph) above for the stem, as Vs increases, U/Vs decreases and should probably decrease more but U/Vs+.333 is never < .333. In a stack, there is no jet entrainment from the sides until after the jet escapes the enclosed stack, here we have no stack walls to prevent jet entrainment and we are also bringing entrained air in from the surface. We adjust Es to match cap height with plume height dh. We find if we replace the .333 in Es=U/Vs+.333 with .2, it raises dh by increasing the momentum term.

            (8) Buoyancy - Momentum - Entrainment:

Why should we be decreasing the entrainment with such higher Vs than most stacks? Because in a stack, the entrainment begins at the stack tip when high Vs meets the air approaching at right angles. The vertical component of the ambient air is relatively close to 0. In heat thermals from the ground, it may be 1-2 mps. So entrainment due to turbulence eddies it very effective here. The barrier offered by the stack disappears at stack height. In a nuclear blast, the stem has no such barrier so there is frictional drag on the adjacent air and what we observed is the resulting Vs with this air, not the sudden turbulence at the top with inert air a few inches away. The ambient air has been dragged along upward for 2 miles. Also, a nuclear blast stem may have a diameter 100 times larger than these large stacks, actually then, 2500 times as much area. So, the horizontal ambient wind effect across the entire stem diameter take longer to impact so the relative entrainment is reduced.

However, if we examine the B3 & M3 term, we note that M3 is larger than B3. This indicates that the momentum term is larger than the buoyancy term that created it. The reaction can't be larger than the action. Also, due to the fact that there is frictional surface drag on the ground surface, building rubble and particles from the violent inflow of wind, the momentum term should be somewhat less. Let momentum to buoyancy be near 85% as a good first estimate. Now, we take the cube root of M3 & B3 first for the M/B ratio. However, the B or M to the total dh term is not so apparent but we are only interested in the M/B ratio. However, we just increased the momentum term, by decreasing Es, so the approximate M/B ratio is now even higher than before. Therefore, we still need to raise the dh term more by increasing the buoyancy term by decreasing Ea. The best results were obtained for Ea=.15 with tg=-.5. The corresponding M/B ratio became 87%. Therefore, the new constant replacing (-.06364) will be as follows:
Let Ea=.15. With .6+(-.5+1.7)*(-.06364), let (-.06364) = x, then .15=.6+1.2*x, .15-.6=1.2*x, x= -.45/1.2 = -.375 so now the new Ea=.6+(Tg+1.7)*(-.375).

Decreasing both entrainment terms Ea and Es increases buoyancy and momentum terms respectively. Since they are both in the denominator, the reciprocal equivalent of this term, will cause the terms to behave in their sense of meaning. Thus we could redefine the entrainment terms as 1/Ea and 1/Es where the results of these terms are now in the numerator where an increase in these causes an increase in dh. However, B & M are more directly affected by Vs/U and 1/Tg. Thus, we increase dh by increasing momentum by increased Vs, but this causes more horizontal wind shear entrainment. We increase buoyancy by decreasing Tg for more horizontal heat difference causing thermal turbulence as seen in the lower PTg stability parameter 's' in the denominator. It's obvious, higher vertical stack velocities cause higher protruding heights and decreased stability causes less negative (downward) buoyancy caused when rising air partials cool by expansion that is, the air partial - ambient air difference is now less that it is expanding into. This is just another fudge factor correction to the main buoyancy term. The cap top by the bronze line in Graph Plots from APT25a was 70380' where we now calculated H=Hs+dh = 3375m+18070m = 70356'. That is close enough since we are within an error less than .1%.

            (9) Result Runs:

Running this program again with a lower yield, the estimated height is too low and for higher yield, too high. This should be true since we have held the wind speed constant when actually, larger bomb yields produce higher stem heights so usually higher wind velocity will lower our over estimate. Also the reverse case with lower yields. The use of real wind profiles should help. Since we have a dh formula that matched the graphed predicted cap height, we may now use this to calculate the absolute cap heights now rather than scale them to these dh formulas. The fine tuning of this equation for other bomb yields using real weather data profiles may be required later.

            (10) Weighted Parameters: Update 03/11/10

Even though the formulas were applied to large power plant stacks, the plume rise of a nuclear blast will probably be at least 10 times as high to the final rise. Instead of the calculated temp, lapse rate and wind velocity nearest the stack height at 11,000', we should consider more of the air that it passes through but still with more emphasis near 11,000'. If we have a slight inversion or wind speed lull at this stack height, the final plume height may be too sensitive to a possible anomaly. Also, this anomaly could be averaged out within a larger layer. What would be a representative height? If we go from 11,000' @ 1 mega to 45,000', this should include most of the dh. So this is 45-11 = 34 e3' to average. Now better than an average, we weight the air more heavily by the weight itself or pressure. We accumulate the pressures (cp) from the layer below 11,000' to 150 mb. Each of the temps, lapse rates and wind velocities are multiplied (weighted) by the decreasing pressure aloft and summated. Then the final summated parameters is divided by the accumulate the pressures (cp) for the mean weighted parameters. Since there may not be a level exactly at 11,000', we choose a standard 700 mb height. End Update 03/11/10

            (11) Modified Stem Diameter: Update 05/23/10

With all the trial runs and the modifications used in (10) Weighted Parameters: above, it appears that the final dh height calculations were still too high. After reviewing all the nuclear cloud pics again, it appears that even though the final stem was 1/3 the cap diameter consisting of the dust inflow, initially the flash diameter may be less. Looking at Nuclear Picture # 5 and Nuclear Picture # 3 again, the long stem does not have that funnel shaped base. Therefore, for the heat flash diameter for plume rise only, we will now use stem = cap/4 instead of cap/3 for the new diameter. This brings most all cases into the graph range from FEMA , but still subject to met conditions that still may cause dh to exceed the graph range a little. This is probable because the limited number of test blasts didn't cover the complete range of winds aloft. End Update 05/23/10

4. Particle Fall Velocity:

       A. Theory:

For particles falling through air, the gravitational force is given by mGa and the frictional drag force by the fluid is denoted by the term Fd. The gravity force begins to fall until the frictional drag force increases to equal it. When this equilibrium occurs, we have a final terminal velocity u called the particle's settling velocity. That means that when a particle begins to fall, Fd force must be increasing while force mGa is constant. What increases the drag force? It is the air flowing along the surface area of the sphere. The more the surface, the more it's drag. Also, the greater the air velocity around the sphere creates more drag on the particle.

The volume of a spherical particle is 4/3 pi r^3 and it's surface area exposed to air friction is 4 pi r^2. It's mass is proportional to it's volume by it's particle density so it's mass (or volume) to it's surface area is 4/3 pi r^3 / 4 pi r^2 = 1/3 r = r/3. This means, as it's mass or r increases, the frictional surface area becomes less important compared to it's relative mass. So larger particles will fall through air faster than smaller ones with the same density. When Galileo dropped the 2 canon balls from the leaning tower of Pisa, they were very dense, hence more mass, and with the relatively short distance of 183', the drag effect wasn't apparent and it's terminal fall velocity was never reached. Therefore, an eyeball estimate couldn't detect the difference speeds, indicated by different landing times. If he had done this with very large and very small cotton balls, it might have been evident - calculate??

       B. Parameters:

Force x distance falling = work-energy = (m*cm/sec^2)*cm = m*cm^2/sec^2 = mV^2. Now average V = (Vo + Vf)/2 where Vo is initial velocity and Vf is the final velocity. Here Vo = 0, so V = Vf/2. In a vacuum then, the Kinetic Energy (K.E.) would be m(Vf/2)^2, so with V0=0, let V be renamed the terminal velocity, rather than Vf for simplicity, remembering that there is no initial velocity. Below, we will use u to replace V. CEH1, CEH2 and CEH3 are Page (5-59 & 60 & 61) from a Chemical Engineering Handbook - 19??

A = cross section area of particle sphere (O=pi r^2) through fluid cm^2, ft^2
Cd,C = drag coefficient
D = particle diameter cm, ft = 2r
r = particle radius cm, ft = D/2
Fd = drag or resistance force to motion of particle in fluid g, lb force
Ga,g = gravity acceleration 1-lb force constant = 32.17 ft/sec^2 or 980.6 cm/sec^2
a = buoyancy acceleration Re = Reynolds number = Dud/v
u,V = relative terminal velocity between particle and fluid cm/sec, ft/sec
v = fluid viscosity = lb/ft/sec or (4.339e-5*Tc+.01691)*1.e-2 g cm-1 sec-1
d = fluid density = m/vol = g/cm^3, lb/ft^3
Dp,@ = particle density = m/vol = g/cm^3, lb/ft^3
m = mass of particle g, lb = 4/3 pi r^3 Dp
Tc = temp in deg Celsius
cm = centimeters
ft = feet
vol = cm^3 or ft^3 = 4/3 pi r^3
g = grams

       C. Equations:

Thus, drag force Fd balances the fall velocity force mGa or mu^2/2 K = mGa. That is Fd = mu^2/2 x some constant K. Above, we mentioned the frictional surface area of a sphere above that is 4 pi r^2. The cross sectional area of a particle sphere (A) passes through the fluid resistance. Both include the term r^2 so are proportional to each other by a factor of 4, or better yet, 4 pi. This increased area friction balances the gravity acceleration term Ga to a steady terminal velocity u, and a stabilized K, depending upon the size of the particle and it's relative velocity. For this constant K, we use the dimensionless C,Cd, a drag coefficient, coupled with area A. So with K = Cd A, Fd = mu^2/2 (CdA). Examining the mass of air in mu^2/2, d=m/vol, m=d vol. Fd = d vol u^2/2 (A Cd). For vol = 1 cm^3, Fd = du^2/2 (A Cd). This balances the weight of the particle mGa so Fd = mGa. Examining m in mGa, Dp=m/vol, m=Dp vol and vol = 4/3 pi r^3 so:

    Fd              =              mGa
du^2/2 (A Cd) = Dp 4/3 pi r^3 Ga, Eq. (1) Journal of Applied Meteorology - 2 Nov 1959, pp 463.
                         Subbing for A = pi r^2,
du^2/2 (pi r^2 Cd) = Dp 4/3 pi r^3 Ga, divide by pi & r^2
du^2/2 Cd = Dp4/3rGa, clear denominator 2
du^2 Cd = Dp8/3rGa, solve for u^2
u^2 = Dp8/3r Ga/dCd, rearrange terms
u^2 = 8/3r Dp Ga/dCd, derived Eq. (2) Journal of Meteorology (JAM0) 2 Nov 1959, pp 463. , sub r=D/2
u^2 = 8/3 D/2 Dp Ga/dCd, GCD 8/2
u^2 = 4/3 D Dp Ga/dCd, square root
u     = (4/3 D Dp Ga/dCd)^1/2, combine denominator terms and rearrange numerator terms
u     = (4Ga D Dp/3dCd)^1/2 = Eq. (5-206) in CEH1.

Actually, it's really Dp-d instead of Dp. For this case, Dp-d ~ Dp, air d=.075 lb/ft^3 at sea level, Dp=2.5 for soil with density = 1 for water so Dp=2.5*62.3 lb/ft^3 = 155.75 lb/ft^3. With d/Dp = .075/155.75 or .05%, means the Dp-d term can be considered as Dp with d < .05% of Dp and 1/2 of that for heights > 18,000'. Here, we have linked and cross checked the JAM article Eq. (1) JAM0 with Eq. (5-206) in CEH1. Below, the Cd term is denoted by C, Dp by @ and Ga by g, not g for grams.

For spheres:
u = (4gD@/3dC)^1/2 = Eq. (5-206) in CEH1. Sub D=2r
u = (4g2r@/3dC)^1/2, multiply constant 2 & 4, square
u^2 = 8gr@/3dC, rearrange terms
u^2 = 8/3r@ g/dC, sub terms Cd=C, V=u and rearrange terms
V^2 = 8/3r@ g/(Cd d) = Eq. (2) from JAM0.

Here, we have linked and cross checked Eq. (5-206) in CEH1 with the JAM article Eq. (2) JAM0.

       D. Plotting Table Values:

Table 5-25 values of CEH3 was plotted for Fig. 5-70 in CEH2 . However, the y = Drag Coefficient C value and the x = Reynolds Number Re value both contain the variable u which we are trying to find. We can eliminate the u term by using Eq. 5-216 in CEH3 , to find C*Re. Solving for C in Eq. (5-206) in CEH1 again:

u    = (4gD@/3dCd)^1/2 = Eq. (5-206) in CEH1, square
u^2 = 4gD@/3dC, clear denominator
3dCu^2 = 4gD@, solve for C
C=4gD@/3du^2, rearrange terms

            (1) Dimensionless X & Y New Parameters:

Proving that C is dimensionless, a dimensional analysis on 4/3gD@/du^2 =
cm/sec^2 cm g/cm^3 / (g/cm^3 cm^2/sec^2), cancel densities g/cm^3
cm/sec^2 cm / (cm^2/sec^2), combine or cancel cm & sec terms
cm^2/sec^2 / (cm^2/sec^2) = 1

In Fig 5-71 of CEH2 , Re=Dud/v
Proving that Re is dimensionless, a dimensional analysis on Dud/v =
cm cm/sec g/cm^3 / [g/(cm sec)], combine numerator terms and rearrange
cm^2/cm^3 g/sec / [g/(cm sec)], combine cm terms and denominator
g/(cm sec) / [g/(cm sec)] = 1. So any C and Re combos are also dimensionless.

C Re^2 = 4/3gD@/du^2   (Dud/v)^2, multiply & combine like terms
C Re^2 = 4/3gD^3@(u^2/u^2)(d^2/d)/v^2 combine exponents of like terms
C Re^2 = 4/3gD^3@d/v^2, rearrange terms
C Re^2 = 4/3gdD^3@/v^2 is Eq. (5-216) in CEH3 .

            (2) X & Y Parameter Values:

Now for Eq. (5-217) in CEH3 , C=4/3gD@/du^2, Re=Dud/v and @-d ~ @

C/Re = (4/3gD@/du^2)/(Dud/v), multiply, cancel, combine exponents
C/Re = (4/3g@/d^2u^3)/(/v), invert 1/v
C/Re = 4/3gv@/(d^2u^3), combine denominators
C/Re = 4gv@/(3d^2u^3) = Eq. (5-217) in CEH3

            (3) A Better Y Parameter:

After deriving the C/Re term, there was no need to do this since it is really not needed. Even though the article suggests above Eq. (5-216) & Eq. (5-217) in CEH3-Application. , we suggest that in Eq. (5-216), we calculate C Re^2, look on a plot of C/Re vs. C Re^2 instead. By using the values in Table 5-25, find the corresponding C/Re value using the calculated C*Re^2 values from the line graph then solve for u in Eq. (5-217). This new graph is almost (log-log) scale so multiple line segments here will suffice. Another similar idea is shown in JAM1 with C Re^2 vs. Re instead of C/Re. We still calculate C Re^2 by Eq. (5-216), since we don't know u but have all the other parameters. We go to Table 5-25 and find the two tabular values that C Re^2 falls between.

            (4) Equations for Log-Log Plots:

In a linear graph, y= bx + a, so y = a at intercept x = 0 and slope b = (del y)/(del x) = (y2-y1)/(x2-x1). On a log-log scale,

y=ax^b take ln's
ln y = ln ax^b, apply log laws
ln y = ln a + ln x^b, apply log laws
ln y = ln a + b ln x, solve for b
b=(ln y-ln a)/ln x, apply log laws
b = ln (y/a)/ln (x/1)

In Fig. 5-70 CEH2 , we have a log-log scale which is the form of y=ax^b. Here, we already have b=-1 by the slope of -45 deg for Re <=.3. Extending any line segment for any b slopes to x=1, y=a. To find the other b slopes of the line segments, b=ln(y2-y1)/ln(x2-x1) = (ln y2/y1)/(ln x2/x1) and y=a @ intercept x=1.

       E. Solving for the Unknown u:

            (1) Method # 1:

Note the tabular values for C,Re in Table 5-25 CEH3 are slightly different from Cd,R in Table 1 in JAM1 . We will use the newer JAM1 table values. If we cover the range of the graph, we want these different segments of the curve. With y=ax^b, here y=C, the drag coefficient and x=Re, the Reynolds number defined between Eq. (5-205) & Eq. (5-206) in CEH1 and the x-axis for Figs. 5-72 & 5-73 in CEH3. Then, y=ax^b, C = a Re^b = a(Dud/v)^b by substitution. Substituting this into form for C into Eq. (5-206) in CEH1, we have the following below. Actually, we decided this wasn't the best method so you may skip it.

u = {4gD(@-d)/[3d( C )]}^1/2, d/@ = .05% so d=0 here
{4gD@ /[3d( C )]}^1/2, sub for C = a(Dud/v)^b
{4gD@ /[3d(a(Dud/v)^b)]}^1/2, combine constants
{4/3da gD@/(Dud/v)^b}^1/2, take b exponent of each variable
{4/3da gD@/(D^b u^b d^b /v^b)}^1/2, convert /^b to -b & //v^b v^b
{4/3da gD@D^-b u^-b d^-b v^b}^1/2, combine like bases
{4/3 d^1-b ag D^1-b @ u^-b v^b}^1/2, combine b term exponents
{4/3 (dD)^1-b ag@u^-b v^b}^1/2, combine b term exponents
{4/3 ag@(dD)^1-b u^-b v^b}^1/2, factor out u^-b
{4/3 ag@(dD)^1-b v^b}^1/2 (u^-b)^1/2, now (u^-b)^1/2 = u^(-b/2) to divide by
u/u^(-b/2) = u^[1-(-b/2)]
u^(1+b/2) = u^[(2+b)/2], now square
u^(2+b) = 4/3ag@(dD)^1-b v^b, take 1/(2+b) root
u = [4/3ag@(dD)^1-b v^b]^1/(2+b)

In this way, we find the corresponding Re vs C, actually in Eq. (5-206) in CEH1, C is in the equation for u, not by appearance and reading off a value from a line which is difficult on a log log scale graph but letting the formula of the line graph (y=ax^b) corresponding to it directly and exactly.

            (2) Method # 2:

The better method would be to have the plot y=Re vs x=C Re^2 in Fig.1 JAM1 . Here, R=Re, Cd=C are the terms. The appearance of this slope is approximately b=-1 on the log-log scale. See Fig.1 JAM1. Here we have the Stokes straight line b=-1 as before. This time, we will find the n-points that Cd*R^2 Table 1 values fall between by CR2n=C[n]*Ren[n]^2 and CR2n1=C[n-1]*Ren[n-1]^2 for x as before but, now we will find the slope b then intercept a of this line segment defining y=R by x=CdR^2. Here in Table 1 JAM1, we add R=.3. For R<.3, slope b=-1, we will let the program calculate the CR2n values rather than use any Table 1 values. This save time and gives us more accuracy and eliminates any between zones approximation of the slope. On this wise, b=(ln y2/y1)/(ln x2/x1) where now we have y2=Ren[n] & y1=Ren[n-1] and x2=CR2n and x1=CR2n1. So we have b=ln (Ren[n]/Ren[n1])/ln (CR2n/CR2n1) and then a=y/x^b = R/CdR^2)^b or Ren[n]/CR2n^b which is a little more complicated for b & a than before. However, the simpler form of u by R is worth it because now we can calculate the corresponding intermediate and noteworthy parameter R for y by R=a(Cd R^2)^b. Then we simply have u=Rv/Dd from the R=Dud/v referenced above.

        F. Results:

The formulas calculated values were checked against Fig. 1 in JAM2 and Fig. 2 in JAM3. The results are very close as seen by


The number stands for the micron diameter. Fig. 2 in JAM3 and our FVTTxxx have the particles falling through layers of the US standard atmosphere so each new air density is calculated and any errors would change the total travel time to the surface. Actually, our table values may be more accurate because their graphs were pre-computer time and they may not have done increments for every 1000' layer and may not be using this derived formula. Please note that Figs in JAM2 and JAM3, the particle size is by radii where our FVTT's has particle size by diameter.

5. Go To Radiosonde Soundings: