Demolition Dust Model

by Buddy McCloskey     4/02/2008     Contact

1. Purpose:

The goals are mentioned in the preceding article World Trade Center Dust Bowl Analysis. We want to extend the results discovered in the article, to apply to any building, if certain parameters are known, then extend to any major city seen in 9 below. As a summary, the WTC Dust Bowl Model was to predict quantity and region of possible dust settlement. Another more important function will be to extend this model to a multitude of buildings in metropolitan area to estimate loose dust generated that will be sucked into a vortex of a mushroom cloud during a nuclear blast so more accurate estimation of downwind radioactive fallout can be achieved in case of a nuclear attack on our major cities. City profiles of buildings heights and spacing density will be approximated and a total dust weight and volume from each city's sampled spectrum of buildings will be used. These tons of dust will be added to the ground crater's dust blown out by the mechanical energy from the nuclear blast as input into the radioactive fallout model RADFO as dust inventory.

2. Building Parameters:

For the nuclear, general conversion constants will be estimated, many of which will be the WTC findings. For the dust fall from a planned building demolition, more detailed and actual values should be used. In any case, the parameters derived or actual are listed below:

Densities g/cm^3
Accepted of Cement 2.4
Structural Light Weight Concrete 1.84 (WTC)
My Cement Block Rub 2.128
Accepted Average Cement 2.12
Reinforced Steel 7.85

% - Ratios
Cement/Building by Volume 3.461%
Steel/Cement by Volume 25.0021%
Actual/Gravity Rate of Collapse 1.6284
Max Fluff Factor 1.873

Building Dimensions Units
Height ft
Length ft
Width ft
Floor Slab Thickness ? 2.5" in WTC
Wall Thickness ? in
No of Floors ? #
Cement Weight * Tons
Cement Volume * ft^3

Special Derived Constants Value
Block Rub # 3 Pressure .5236 PSI
Block Rub # 3 Width 4.25 in
Block Rub # 3 Rub Length 98.4 ft
Block Rub # 3 Tons/Sec 1.767e-8
Block Rub # 4/# 3 Pressure Ratio 2.38861
Block Rub # 4/# 3 Grams Ratio 6.01928

For the building dimensions, * means one or the other, ? indicates sensitive but important data for accuracy. Also, if any others values are known in order to replace the densities or ratios listed.

3. Dust by Physical Processes: - B%D.c

In demolition, we assume the base of the building is destroyed and falls upon itself. It might be best to look at this general process here in trap3. If these displays are shifted and jumbled when you open them up with your web browser processor, open Word Pad or some other word processor, open them in MS Accessories - Notepad. If by default you can't, you might right click and save the display as a text file, then use Notepad.

The general idea is that the building walls are crushing the walls below them, cumulating rubble that falls equally on the floor slab within the walls and outside the walls - cumulating a large pile on the outside that rests against the accumulated slab heights - the "pancake" effect often mentioned. The shape is a rectangular prism with base radius Rc, height Hw in trap4. The R's for radius rubble and to represent another term later. Notice, the building and the next floor falling down.

Note the more detail in trap4, we note two processes. We see the chunks of cement and the airborne dust is divided equally to either side of the crushing wall. This is indicated by Dust # 1. Here, we have the bottom of the wall sliding in both directions < - >, causing dust by the rubbing surfaces under great pressures, but assuming here that is not enough to greatly exceed the compression limit of 6000 psi because the rubbing changes it to powder dust instantly. We calculate the dust by the wall edges rubbing, and the remaining chunks are divided equally between the outside and inside piles. The first results of compression is the crumbling of chunks caused by > 6000 psi on jagged edges.

In the lower diagram, we see the secondary processes at work indicated by Dust #2 & #3. Now Dust #2 is caused by rubbing against the outside cone pile or lateral surface area down side (hypotenuse L) of the triangle. We show spaces here but actually, the chunks are continuous extension of the wall with the total pressure still on the chunks as they traverse the outside wall, because the time between slabs falls is less than .15 sec or 62 mph.

To the right, we have chunks accumulating above each floor slab (tbs) until the floor slab above meets it at height Ht (height temporarily trapped). When the slab falls on each of it's respective pile, it's the crushing action by the crush power formula in use to produce dust #3 to combine with dust #1. Both of these are extruded with compressed air puffs explained below. If there were no smaller pile, the slabs would stack like pancakes, and land flat, producing no side rubbing. If they did slide, like the walls, the area is immensely larger but the psi is less than 50 initially so not as much dust produced as one might think.

4. Deriving Other Parameters from the Building Parameters:

       A. Dust Puffs

Many derived parameters were developed in a forerunner program before the final model here, not knowing if the data was needed but interesting to shed more light on the physical process. We will present most of them here and tell if was used in the final or not. Actually, we can view the data in the summary stat sheet ddmdat.

As a standard, we expand the WTC study again to this application. Referring again to the falling slabs, besides crushing the inside pile, the air volume can't be crushed into smaller chunks but is compressed. The inside office air volume is the slab area, length x width after subtracting out the wall thickness so the Inside Wall Slab Area (iwsa) = 43053.4'^2, then times the ceiling height which is the avg floor height - slab thickness. Thus the office volume is 5.253e+5'^3, neglecting the office furniture, wall dividers, stairs and elevators. This volume is compressed in just .134 sec or a flow volume of 3.918e+6'^3/sec. What happens, the windows are blown out and this dust blow out should be visible, which it was. Based on the windows alone, the surface area of the windows/floor is 3328'^2. So dividing this into the flow rate gives us the windows or broken glass puff distance of 157.8' occurring in .134 sec = 1177.2'/sec = 802.6 mph if no walls or ceilings collapsed.

Actually due to the crumbling walls and slabs, it is less than that. The falling rubble and dust umbrella were so wide, that only peaks of the puffs were visible. Screen freezing the videos when the puffs occurred and measuring them on the PC monitor, then ratio them to the building side, it appears that puffs near 45' were observed. Thus, this is 45/157.8 = 28.5% of the broken glass window openings account for the total holes, so 100% - 28.5% = 71.5% of the holes were through the cement openings. Apportioning this amount to the cement walls is wall area-glass area = Cement Wall Area = 6823'^2. The Inside Wall Slab Area = 43053.4'^2 from above so Slab/Total Cement areas ratio = 86.3% of all the cement areas (71.5%) = 61.7% of the holes were in the cement slabs, while the walls account for the rest or 100% - 61.7% = 38.3%.

By the WTC-911 videos, it appears that the walls crumbling rapidly so in reality, the walls should contribute the most and the slabs minor. However, the trapped dust puffs between the slabs are quite dusty and when any dust goes through the breaking up slabs, it escapes to the next floor level above to be extruded out the window and walls with some escaping to the next floor. Finally, as the last ceiling collapses, all the trapped dust escapes to the atmosphere. We calculated 53.8 tons @ a 60u dust layer total for the sum of the slabs. With steel beams falling in between, it becomes be more complicated.

       B. Separating Slabs and Walls

We treat the height of the rubble pile as the accumulated slab heights = 22.92'. The volume of each slab is 2.5" times length x width = 9.915e+5 ft^3. Subtracting this from the total cement volume leaves us the total volume of the walls. Pretend, we unravel the enclosed building area by stacking the 4 walls end to end (perimeter) at building height. Now to find the wall thickness, we divide the volume by the area perimeter x height to find the thickness of this new shaped volume. However, a minor adjustment, we must subtract off the heights of each slab or the stacked slab heights here for the new height - walls with no slabs. This gives us a wall thickness of 11.3". Multiplying the cement density by volume of each slab = 517.686 Tons and each floor's wall volume = 9.383 Tons.

        C. Rub Areas

Also needed is the rub areas of the wall. Since the slabs are accumulating vertically, we assume no horizontal rubbing of the slabs. The wall rubs on itself and the residual goes inside and outside the wall. Simply, the area is the wall thickness x perimeter = 784.18'^2 (wtcaw). The rub 3 block rub length is 98.4' x it's width 4.25" from the above table gives us an rub3a area of 34.85'^2 so the top of the wall area ratio is rarb=wtcaw/rub3a = 784.18/34.85 = 22.50 times as much dust as the rub 3 block area (rub3a) which by hand rubbing alone would be minute.

        D. Rub Time

We loop up the building height by the slab thickness. This gives us a better approximation for the wall dimensions but allows us to evaluate the whole slab when it's turn comes. So for the WTC, 1365'/(2.5/12) = 7872 evaluations. If it takes 15 sec for the building to fall, then we only evaluate for 15 sec/7872 = 0.00229 sec for the 2.5". This is multiplied by the Block Rub # 3 Tons/Sec (1.767e-8) to scale it down to the actual time the 2.5" of wall was rubbed. Thus now this appears very minor.

        E. PSI Application

As we loop up the building, actually building descending, it is better said that we are evaluating every 2.5" as it collapses, so we examine every 2.5" section increments up the building when it arrives at the ground base. As each 2.5" collapse, the weight of this 2.5" section is subtracted from the building weight for the next 2.5" evaluation. This reduces the PSI wall rub application. For the 1st 2.5", we calculate a PSI on the bottom of the wall perimeter as 2070 (press) psi to the Block Rub # 3 = .5236 psi from the table above = 3953.42. Now the Block Rub # 4/# 3 Pressure Ratio from above = 2.38861 (prat0) so for prat0^ex = prat, ex = log prat/log prat0 = 9.51215. From the above table, the Block Rub # 4/# 3 Grams Ratio is 6.01928, so we have 6.01928^9.51215 = 2.601e7 (doratb).

        F. Combining Factors C., D. and E.

Finally, we have 2.601e7 x Rub Area Ratios (22.50) x Block # 3 Tons/sec (1.767e-8) x Rub Time (0.00229 sec) = .0237 tons. For the next n (2.5") of the loop, press = 2069.84 psi. All the other constants are the same, only the new dorat = 6.01928^ex, finding ex as we did previously. In this first case, .0237 tons is just the maximum potential of dust that could be produced, denoted by the term tondb (tons of dust by building - wall rubbing on itself). Recall, in B. above, the tons/wall per slab thickness height (2.5") is 9.383 tons (tpw). So, we can't rub off more than what's available but we have .0237 << 9.383 so OK here. It doesn't appear that we would have a building where the potential tons rub is greater than what was available if the wall thickness is large enough influencing the psi's and tpw. In any case, we use the code check "if(tondb > tpw) tondb=tpw".

        G. Budgeting Dust and Residue

We accumulate that tons of dust from the building source by the term 'cumtd'. Now the solid chunks that weren't turned to dust is left as tonw = tpw-tondb. As we mention before, we split this to accumulate outside the wall and inside the wall so each receives tonw/2. Redefine tonw now as original tonw/2. The inside wall accumulating term is tons between slabs (tbs) incremented by tonw. It is reset to zero once a slab has fallen on the pile.

        H. Working With Volumes of Tonw

               (1) Outside Wall Rubble Volume

We need to convert tonw to volume (tonwv) '^3 by dividing by density. We dump these solid chunks to the pile outside the wall. Note in trap1, trap3 and trap4, that the rubble pile against the height of the accumulated slabs (Hw). Now note in trap1, a vertical cross sectional mid length cut shape of the outside pile while trap2 shows the top view. When we dump tonwv/2 on this pile, it is distributed all along the perimeter 2 Lr lengths and 2 Wr widths. However, since the perimeter wraps around at four 90 degree angles, what happens at the corners? Here we view four 1/4 cone bases indicated by @.

Thus we have the volume sum of 4 trapezoid prisms + 1 cone. For the cone, Vc=pi*Rc^2*Hw/3 and for the trapezoid Vt=4(1/2Rc*Hw)*2*(length+width). So total Vol = Vt+Vc = pi*Rc^2*Hw/3 + 4(1/2Rc*Hw)*2*(length+width) = pi*Rc^2*Hw/3 + 4*Rc*Hw*(length+width). Multiplying by 3, 3*vol=pi*Rc^2*Hw + 12*Rc*Hw*(length+width). We finally have a quadratic equation Rc^2*Hw*pi + Rc*12*Hw*(length+width) - 3Vol = 0. Solving for a quadratic ax^2 + bx + c= 0, x = -b +- sqrt(b^2-4ac)/2a. So here we have if x=Rc, a=Hw*pi, b=12*Hw*(length+width) and c=-3Vol, we solve for the positive root Rc = (-b + sqrt(b*b-4*a*c))/(2*a).

               (2) Inside Wall Rubble Volume

Before we can apply Hw in (1) above here, we need to know the height Hw of the rubble from inside the wall. We could define Hw as the accumulation of slab thicknesses if all the inside wall rubble (tbs) is turned to dust. This dust has to be calculated. For the loop, we start out (n=0), and accumulate the 1st floor slab. Also, we later make a cone pile out of tbs but initially, we have to set a tiny tbs=.1 ton so the math routines won't have to work with 0's and terminate the program. Assume, we have a few floors collapsed already and we have some outside wall height Hw. Now we accumulate the wall rubbing tonw to the tbs volume or tbsv '^3 by dividing by density. Referring to trap4, we want to find Ht and Rb of the triangle rubble against the collapsing inside wall. This triangle shape extends all around the inside perimeter of the building wall. So a cross section area cut (cac) = tbsv/2(length+width).

By examining crushed stone piles (1" diameter), we seem to have a consistent Rb/Ht, or the inside cone Ht, has the ratio near 9/10 call it tana, the tangent of it's angle of 42.0 degrees looking down from it's peak. Actually, these chunks are probably greater than 2-3" but the roughness of the cement will provide more friction with one another to produce a possibly steeper height. Actually, we tried to determine this angle with 3 types of sand. With a moist sand, we determined this angle to be 45 degrees. With a wet sand, it first flows like a liquid and makes a pancake type base. As the water keeps dripping, the dryer sand makes more consistency and the height builds up. Actually, this type of sand makes those high peaked dripping sand castles. We tried to keep a consistent slope angle so settled with a peak (smallest angle) of 35 degrees. For dry sand, the grains don't stick like being moist but flow like liquid wet sand in the beginning. Actually, this was controllable and stabilized with an angle of 59 degrees. Then, continued pouring the dry sand until the dimensions where increased by over 2.5 times, the angle was stabilized at 57 degrees. So averaging, we have 58 degrees dry, 45 degrees moist, and 35 degrees wet. So the stone pile angle of 42 degrees is very similar to the moist sand.

This is all due to the material trying to pour like a liquid but the increased horizontal friction caused by adhering of multiple sand grains when semi moist or the rough irregular shape of stone particles digging vertically, causing higher psi in small columns inhibiting horizontal movement of adjacent stones. Looking from the base, the slope upward is the 90 degree compliment or 32, 45 and 55 degrees for dry -> wet sand and 48 degrees for crushed stone. Visualize this base angle (Ht/Rb) instead of looking down from the angle of the peak then dividing into halves, indicates the material stacking is better. So generally, the slope is 1 near 45 degrees (Ht=Rb) for most material, depending on it's roughness for larger chunks and less than 45 degrees for smoother dry material. We want this cross section triangle to be in this ratio. So area of triangle (top view) = Rb/2 x Ht = cac. Substitute Rb = Ht*tana, we have cac=Ht*tana/2 x Ht = Ht^2*tana/2, Ht^2=2cac/tana, Ht=sqrt(2cac/tana).

                    (A) Crushing the Inside Pile Theory

We use the energy of the slab (sle) to crush this inside pile of rubble. The weight on the slab is the weight remaining from the building on this slab which travels a distance Ht for crushing. No work is done until the slab is in contact with the pile. Thus sle is in ft-lbs, converting to kwh by multiplying by 3.76616e-7 kwh/ft-lb. Now recall in WTC Dustbowl Model, we had a crushp=1.15 KWH/Ton converting all the particles to the particle distribution profile produced from the block rubs and considering 1/2 this value with cement on cement instead of metal on cement. So potential tons by crushing is then tonbcp=sle/crushp. Is this crushp in KWH/Ton based on collecting the amount of dust during some KWH of electricity by noting meter readings? If this is the case, is the machine efficiency taken into account. With just a wild guess, for now, we might use me=.25, a 75% loss of efficiency by friction.

A Mr. John Dochner from Rohrer's Quarry in Lititz, Pa. has forwarded me some limestone crushing parameters seen in rohrerq.dat. Using the data, I generated more parameters in QuarryCrush.dat. This allows redoing the crushp=1.15 KWH/Ton from above. In this dust created table, let's normalize the energy for the particle sizes back to the standard 60 micron (u) to compare to Jerry's Russell's 1.5 KWH/Ton, referenced at the bottom. The factors are under the sqrt(60/D) column, the resultant reduced tons under the tons column while keeping the same Kwh. We also could have kept the same tons but increased the corresponding Kwh. By the first method, the total tons were reduced by 49%, resulting in a raised KWH/Ton to 4.75 Kwh/Ton @ 60u for all 4 groups. By the way, the % dust = 6688/84811 = 7.89%.

Before we do anymore adjusting, the note by Jerry Russell that cement is softer, reduces 20 kwh/ton to 1.5 kwh/ton. Now the hardness factor by Rohrer's Quarry is Mohs Hardness 4.0. Jerry Russell's reference is In Table 2 there, we adjusted the long ton to short ton to increase the kw-hr/lt to kw-hr/st by multiplying by 1.12 for Blast = .482, Crush = 3.629 and Grind = 19.958 Kwh/Ton. Jerry must have chosen the 19.958 for his referenced but corrected to 20 kw-hr/st because of the 2cm --> 60u size. But this is Grind which we are treating independently. However to go from 20 kwh/ton to 1.5 kwh/ton based on hardness seems extreme plus the fact that grinding (rubbing) is so different than crushing. In grinding, the completely dislodged particle distribution is produced, the amount increases with more rubbing or more pressure but the same particle distribution. With crushing, the particle distribution produces large chunks, then with more crushing, the particle distribution shifts to smaller sizes.

Using Rohrer's Quarry 4.75 Kwh/Ton @ 60u to start here may be better. Now to reduce 4.75 because cement is softer than limestone, looking at hardness, we used a copper penny and fingernails to rub on our cement rub blocks. For the penny, the edge was worn off rapidly but most was cement. It appears as 4/5 was cement, 1/5 copper. For the fingernails, they went to stubs immediately and no cement rubbed off. So by the hardness table and note at the bottom, cement > 2.5 (fingernail) but < 3.5 (copper penny). Let's use that 1/5 of the interval below 3.5 or 3.3 for the hardness of cement. So limestone hardness = 4.0 to cement hardness = 3.3, we want to reduce the 4.75 Kwh/Ton to what? Simply linear, 3.3/4 = 82.5%. This is certainly a less decrease than from 20 to 1.5. Mostly all power and energy terms use a square somewhere. With the hardness, maybe the term squared would indicate the more crumbling, with less loss of energy crushing through the rock. So .825^2 = .681 x 4.75 = 3.235 Kwh/Ton. Now with an electric motor at 25% efficiency, we have .25 x 3.235 = .808 Kwh/Ton. Then finally, instead of cement to metal crushing surfaces, we have cement to cement, making it twice as efficient so we require half the energy and it becomes .404 Kwh/Ton. Before, we used the caf = 1.53 to upgrade to our observed rubbed particle distribution, but we can use this recent data because we will classify all dust < 600u. We will redo the next program using this adjustment to help in the building dust particle distribution again but for now, we will just call it all dust. Besides the electric motor at 20% efficiency, we don't know if the feed is at 100% volume in the crusher at all times. Spaces or gaps could reduce this further so even at 75% full, we would have .303 Kwh/Ton.

With all this, we looked at another approach using Rohrer's Vertical Shaft Crusher. Here we will use the more accurate HP output capacity, say 85% of this max operation. The rubble is coned to nearly 100% volume from 3/8" to 1 1/2" particle sizes, no estimate of halfing the value because it's stone upon stone. We still need the hardness factor. Also, 100% of the input load must go to dust. The HP rating is Motor = 250 hp x 745.7 Watts/hp = 186425 Watts = 186.425 KW and 186.425 KW/3600 sec/hr = .051785 KWH. We have 12 tph of #30 (600 micron) dust to crush so .051785 KWH/12 tph = .004315 KWH/Ton. If run at 85% of HP rating, .85 x .004315 = .003668 KWH/Ton. This is 1/83 of .303 above - a better estimate of both methods.

                    (B) Limits to Crushed Dust

We can't produce more dust than tonbcp and we have the pile (tbs) and the slab (tps) to work with. We set a variable ton = to the minimum of tbs and tps. With ton set to the minimum of slabs or pile, the crushing is effecting both pile and slab so as with the rub, dust is being generated from both so the maximum dust will be double the minimum between them or 2*ton. We may have to reduce tonbcp from it's potential to an actual by the following code. We have if(tonbcp > 2*ton) tonbcp=2*ton, bringing the max potential tonbcp down to the actual 2*ton. As an example, at the onset of the 2nd slab, tonbcp=945.7 with tbs=276.1 and tps=517.7. The min ton=tbs=276.1 which will usually be in all cases. Now 2*ton = 552.2 < potential (945.7) so we reduced the potential tonbcp down to 552.2 tons of dust produced by the slab crush power. Later, when there is less weight remaining from the building, the potential (tonbcp) is less than 2*ton. In this case, the actual tonbcp is limited to it's original potential (sle/crushp). This happens after 46 stories collapse so this procedure is crucial.

                    (C) Budgeting Dust and Residual

Again, we immediately accumulate this crushed dust to variable cumtd, where the wall dust accumulates also. Now with the residual chunks, we add it to the other pile inside the wall (cumti). How much debris do we add? Well we took 1/2 from each pile and slab source so what remains is (tps-tonbcp/2)+(tbs-tonbcp/2) or tps+tbs-tonbcp. The first expression explains what is really happening where the 2nd is mathematically equivalent and easier - or could say, add slab and pile then take away the dust they produced but it doesn't happen in that order - we really add what is left in the slab after this dust and also what is left in the pile after an equal amount of dust (276.1 tons) or 517.7-276.1=241.6 tons left in the slab and none in the pile between slabs in our example. We then reset the tons between slab (tbs) to 0 for the next slab.

                (3) Outside Wall Rubble Dust Process

Procedure and values are different than the rub wall application.

                    (A) Changing Rub Areas

Now that we have a wall of debris from (C) above, we have it's height Hw, where we apply the quadratic roots formula in H.(1) above. Here we have the base or cone radius Rc in trap4 by solving the quadratic. With Hw and Rc known, we have L = sqrt(Rc^2 + Hw^2), the length the residual tonw slides down the side of the outside pile (Airborne Dust # 2). The lateral surface area of the outside rubble trapezoid is L x perimeter or L*2*(length+width). Now L is also set to the lateral length of the four quarter cones whose total lateral surface area is the 1 cone pi x Rc x L. Both these areas are summed to represent the rub working area wca. Now just as we had in 4.C. Rub Areas above for the wall, we now have the debris trapezoid - cone area to the rub3 area as rarc=wca/rub3a.

                    (B) Changing Rub Time Factors

With Hw and Rc, we also have the 1/2 angle of the cone peak trap-cone as angc=Arc Tan(Rc/Hw). For doratc, the psi component normal towards the side of the trap-cone is sin(angc) so prat is reduced, so is ex and the resulting doratc. However, the velocity component may be slightly faster because of less friction from less pressure so we use a factor of doratb/doratc to approximate the speed up the velocity. However, the velocity component is now avs*cos(angc) where avs is the average speed of fall @ bh/sec from Actual/Gravity Rate of Collapse = 1.6284 from table above = 91.0'/sec or 62 mph done above in first paragraph of 3. Dust by Physical Processes:. So for the travel time secc along L, secc=L/(avs*cos(angc)*doratb/doratc).

                    (C) Combining All Factors

As with the wall rub above (4.F. Combining Factors), we still have the same physical terms but different values so different variable names using c's as in tondc=rubtps3*doratc*rarc*secc.

                    (D) Limits to Rubbed Dust

Again, we are careful to account for all the tons - dust and pile residual. This is similar to the dust caused by crushing in 4.H.(2)(C) above but important differences. Here the potential max dust is represented by tondc. Again, we set ton to the minimum of the two rub sources this time, tonw again - the residual from the wall dust rub and cumto - the cumulative tons of outside rubble against the wall. Recall, tonw also went to inside the wall to be crushed by the slab, this tonw will be rubbed by sliding down the trap-cone lateral side or rubbed the 2nd time. Again we have "if(tondc > 2*ton) tondc=2*ton" bringing the max potential tonbdc down to the actual 2*ton.

                    (E) Budgeting Dust and Residual

Now here, we want to add what residual is left to the outside wall cumto pile. This residual is tonw-tondc/2. The other half (tondc/2) of the dust was from the cumto pile itself. So we have another tondc/2 to subtract from cumto. Altogether, we add (tonw-tondc/2) - tondc/2 or tonw-tondc. Again, the simpler last expression conceals the actual process. We add tonw less half the dust that was rubbed off from it to the pile, which also contributed half the dust. The difference here is that the pile that we add the residual to is growing and provides a larger rubbing surface area itself rather than in the crushing, where we crush the whole pile between slabs each time. Also, we add the generated dust tondc to total dust term cumtd.

5. Results:

Reviewing the bottom of ddmdat, we see in Rubble:, we listed the tons of residual from within and outside the wall and total dust from all sources combined (cumtd). This should total the tons of the original building before it was demolished. Corresponding and underneath is the % of each contribution of tons. It is seen that for the WTC, this method produced 62.6% to dust compared 72.5% in the wtcdbm model. This lower value is interesting but acceptable. In the previous, we took a % of the total cone volumes for each foot of building fall. Since we used the resulting volume to help predict the total volume, there were some bias built in. However, there was some merit to quantify the physical process that we could use here. In this study, we used actual psi on walls when they occurred and treated the slabs separately as a crushing estimate. The pile and the rubble were accumulated but had no bias built in. We will be examining some variables more accurately as mentioned above in 4. Deriving Other Parameters from the Building Parameters:.

6. Adding Steel & Fluff to Pile for Final Dimensions: (Not Used in Final Model)

We will add steel (vols) to the inside of wall only as that is most likely. Now vols=psv*volc where psv is Steel/Cement by % Vol in table above. Thus, the new inside volume iwv=cumti*2.e3/den + vols. Now the new hw=iwv/iwsa where iwsa = inside wall slab area as used above. This means we have a new height for the total complex outside wall volume owv. Again, we apply the quadratic root equation in 4. H. (1) above for the radius/base Rb but substitute owv for Vol.

For the fluff factor, let's just boost owv directly because of the crushing slabs, more air spaces would be jammed up and filled. Also, since the WTC pile was maxed at 3 stories high, lets not increase that anymore. So we solved the quadratic root equation for Rb again and substitute owv for Vol again. Only the Rb extends outward but still short of the observed. The real wide umbrella cone in the WTC really brings up questions and a mystery when it should have just fallen vertically. Could it be that the added energy of the plane and fuel burn, primed the WTC for greater efficiency for when the gravity took over? It should appear that the fall velocity would be greater but the extreme lateral spreading indicates some other energy assistance was supplied horizontally internally?

7. Particle % Distribution from Crushing:

Above, we calculated % dust but assumed the particle distribution produced was from the rub blocks @ 12,800 dpi plot. This time, the dust produced from inside the wall by crushing will have a different particle distribution. In rubbing, the finest individual particles are removed and become airborne. The rubbing of larger particles breaking down to smaller is assumed minimal. In the crushing process, larger particles are breaking down to smaller particles, a different process than rubbing.

From the quarry, we refer to QuarryCrush.dat again. In the dust created, for the Vertical Shaft Crusher we have 4 Sieve Sizes - or screens with equivalent Diameters D. We see the % Passing, what was left = % Not Passing and the difference between these sizes is the % Trapped Particles. Taking the mid values of the classes down to 75u should represent some sort of distribution. The final dust (< 75u) was also sent to me from Rohrer's Quarry. So from 75 u to the smallest would complete this crushing particle distribution. We used the method previously as explained in World Trade Center Dust Bowl Analysis under 2. % Particle Distribution:. In order to scan the limestone sample, it appeared that if we were to drop pinches of the sample on the plate, then the individual particles would separate. Not so, they still were seen clumped together before the scanning. The only way was to make the sample airborne. So the method was to blow it out of the bag and onto the plate glass for scanning. A scan sample of this limestone is seen in # 30 dust @ 12,800 dpi plot. Notice, the 2-tone dust color this time has to be programmed to combined into one dust particle. The results of the scan are combined with above to a complete % distribution graph shown in sievecpd.gif using a Bezier smooth again. The total # of particles scanned was 3001, for the dust < 75u. Rotating and scaling again, reads the pixels along the graph to produce a particle distribution file we will use as simply crush.%d and rename the earlier rub profile as rub.%d. Recall our crush adjustment factor caf - the power needed to crush a particle distribution from a departure of all particles at 60u. With the quarry profile as seen in sievecpd.gif, caf=.593, because the bulk frequencies of the particles mean > 60u. We won't use this term but it is referenced for comparison to the rub profile part dist % frequency of caf=1.53.

        A. Split Rub and Crush Outputs in Program

We now modify the program DUSTBM.c (DUST Box Model), to output separate files of total dust from the outside wall rub and inside wall crush components. One output file is the total dust as before. Every time a slab crushes the inside wall residual pile, another output file contains variables rcaf and tonbcp. To refer to these terms, we view the diagram and explanation in crush.dia. Here we assume a slab travels 2 thicknesses of itself then breaks up. Before breaking up, the slab is attached to the building taking all the tons weight above it (ctr). So for this distance (2*slabt), x ctr*2e3 = ft-lb = sle (slab energy). For the rest of the distance, (hpi-2*slabt), we have the remaining tons/slab (tps) weight applied which is rat = 4*stl*(length+width)/iwsa where stl=cac/slabt, the cross sectional area cut of the pile (cac) divided by a new height of 1 slab thickness.

That is because the rubble triangle base is spread by being squashed like a pancake with height = 1 slab, and length stl (slab thickness length). The 4*(length+width) is the perimeter as used before so 4*stl*(length+width) is the area along the 4 sided building perimeter. This area to the inside wall area (iwsa)=(length-wallt)*(width-wallt), is the ratio rat. This detached energy contribution is minimal compared to building weight original sle but added to it by sle = sle + (hpi-2*slabt)*rat*tps. Now the total sle is converted to KWH by sle/ft-lb*3.76616e-7 kwh/ft-lb. Now then tonbcp=sle/crushp where crushp=.004315 KWH/Ton is Rohrer's Quarry Vertical Shaft Crusher from table QuarryCrush.dat and tonbcp = Tons of dust by crush power available. Actually, this is the potential (poton) for that production depending how much cement is also available as explained in 4.H.(2)(b)

In this update, we calculate the remaining tons per slab (tps) or rat*tps now instead of just tps because of the breakup of the slab crushing the residual pile. As before, we choose the minimum (ton) of the two (rat*tps or tbs= tons of rubble between slabs), then each (crusher and pile) supply ton or 2 ton total. Now before we used if(tonbcp > 2*ton) tonbcp = 2*ton which is still good but also, how many times can the actual tonbcp be crushed if available? Well that simply is rcaf=poton/tonbcp but of what use is this? Now we have how many rubbles this size could we crush to produce the quarry's dust profile as seen in the graph sievecpd.gif. Since there is only one rubble pile here, we use the remaining energy sle to keep crushing this pile to smaller particle sizes because we ran out available cement to crush.

        B. Keep on Crushing - CRUSH%D.c Update 10/22/08

In another program, we read the Quarry's Dust Profile values interpolated from the graph sievecpd.gif as we did before with a few programs in World Trade Center Dust Bowl Analysis under 2. % Particle Distribution:. We call this file "crush.%d" and will be used in all subsequent programs. Suppose we want to extend the crushing beyond this profile to produce a new particle % distribution profile reduced to smaller particles. Well we could determine how much energy was available and keep applying the sqrt(d/id) for each d until the crushing energy is matched. Rather than doing this over and over, lets determine this beforehand by crushing this original file crush.%d in increments to produce a different profile for each corresponding energy crush increment.

We determine how much crushing energy is produced from each slab and find it's corresponding calculated compressed energy - % profile increment. However, if crushed less than that required to produce crush.%d, we have larger particles, rather than smaller. We won't extrapolate this backwards here but remember that more larger dust particles could be available with less crushing. In the expression crush=sqrt(60/d), lets change the constant 60 to any diameter d and d to a new standard constant id. To crush from d microns (u) to id microns (u), requires crush=sqrt(d/id). Or in other words to increase the energy factor of sqrt(60/d), we divide by sqrt(60/id). That is sqrt(60/d)/sqrt(60/id) would also be sqrt[60/id)/(60/d)] = sqrt(60/60 * d/id) = sqrt(d/id) again.

In this program, we read the quarry's crush.%d for d and % particle distribution pd. We start and reduce all d's to 1u by crush=sqrt(d/id), then times pd (weighted) then cumulated by ccrush = ccrush + pd*crush. Since the sum of pd's = 100%, then the sum of crush's (ccrush) is the final weighted energy required to crush all particles in crush.%d down to id=1u. We won't take all the particles to that extreme and may not have that required energy. Also, the WTC photos of 911 show dust particles settling, where 1u particles would essentially stay afloat in the observed downwind distance.

Starting from 2"-3" cement chunks, we have crushed to a profile as seen in sievecpd.gif again, it appears as a double normal curve. Other normal curves are seen in NormDis.png. Here u is the mean=median=mode at u=0, and u=-2. The standard deviation is sigma but shows the plots of various sig^2's. The formula for a normal curve is given by pdf.png. A sig represents one standard deviation of the data. Now the area of the peak to 1 sig on either side of any normal curve represents 68.26% of the data. If the area under the whole curve = 1.0, then -1 sig to + 1 sig = .6826. The exp( ) part formula above can be rewritten as exp[-(x-u)/sig]^2/2 where (x-u)/sig is the number of sig units we call (ri). So if x=10, u=4 and sig=2, then ri=(10-4)/2 = 3 sigs. The area of the curve between + - 3 sig is 99.74%. Actually, most use 2.5 sigs to represent all the data (98.76%). For our case, we require not to lose any dust so we choose + - 5 sigs = 100.00%, 100.000% or more decimal placed zeros.

What do we plan to do with this normal curve? We want to break particles in half, or try to, that even though most won't be exactly half, determined by what precision you judge by, they will distribute about the mean=mode=median of 1/2 the size to crush. If we consider a particle size of 504u, which we won't use that large in our model, then 1/2 that or 252u will be the mean=mode=median as indicated by pb2.gif. Here we don't have the sig but estimating 68%, we might eyeball guess sig = 65u to get the feel of a normal distribution. If we want ri = 5 sigs to represent all the data, then sig must be 252u/5 = 50.4 u, less than our eyeball estimate.

In our program, let's use 1u intervals, determine the size to crush (max), and ri= -5 sig -> 0 sig -> +5 sig or 0 u -> max/2 u -> max u or 0 u -> 252 u -> 504u. So we have sr2pi=sqrt(2*pi) and sig=(max/2)/5, and ri= (max/2 - i)/sig where i=particle diameter index [i] and u's. We loop i=0 -> max-1 and ri2=ri^2. We have y or erf[i]=exp(-ri2/2)/(sig*sr2pi). Two examples are seen in erf.dat. The simple 1st case is sig=1u centered around diameter = 6u. Notice the symmetry of the erf x 100 values and the cumulative area of 100.0%. In the 2nd example, we still have the cumulative area of 100.0%, symmetry around 26 u but many more values with sig= 5.00 making a wider bell curve similar to the blue plot in NormDis.png.

We decide to use diameter (max). What decides the size max? Well, we could just start with the max diameter size, crush a bit, then the next largest size down to 10u or less. The most likely of sequence particles to be crushed will be all the volumes occupied by particle diameter i. Simply, this is the % particle distribution pd x volume (vol) where vol=4/3*pi*r^3. However, vol x npar (%) is also termed particle radio activity so pact[i]=pd*vol for each particle i to np. For our 1st max, we read diameter or index i & pd from "crush.%d" and calculate vol x np, for pact[i]. For our 1st case, max pact[i] (or pd*vol) was np=346u.

Here mean=mode=median max/2 = 173u. So the 1st crushing will be on the np max = 346u particles. How much of this volume should we crush? If we crush all of it, that is too much and we might show a rough % particle distribution. Too little, and the program and data output files will become too large. We can take a percentage, so we did and the best choice is 20%. Why? Well with the 6000 loops, the np=346u repeated about 5 times as the max? This is what probably happens in the real world. Just 20% of the volume is taken off the max particle volume and then after 2 or 3 times, there was another size to replace 346u as max. Then 346u comes back a few more times as max again until most of the 346u particles have been depleted by being crushed smaller. So even though vol is small with smaller particles, the np has accumulated immensely and hardly any more larger particles are to be crushed.

We find the total volume occupied by 346u's as vols=tvol*pd, correcting pd(%)/100 where tvol is total crushed dust volume. Lets use delta volume (dvol) of vols. For dvol, we use a factor 20% = .2. So we have dvol=.2*vols. We take this dvol value from vols. What we do now is to reduce the % particle distribution of 346u by pd[346]=(vols-dvol)/tvol then * 100 for %. Next we distribute dvol by the normal curve's erf values as just seen by the two examples in erf.dat. To see this probability or % distribution normal (pdn) applied to all particles < 346u, see crush.346. Note, at the bottom, it is subtracted from 346u, while for the other particles it is added. The added effect on pd x 100 appears nil sometimes showing no change within the 4 decimal % near the tails. However, as you approach the max/2 = 173u, it's added effect (->) by the erf's pdn is the greatest and diminishes towards the beginning (1u) and end (345u) of the profile tails.

We call this the crushing loops l. This is repeated 6000 times so there are few larger particles to be crushed into smaller particles. If we go more than 6000 crush increments, they jam into the 1-10 u range and becomes a blocky distribution. Also, 1-10u particles are considered to stay aloft for many miles so no need to distinguish between these small sizes and neither did the earlier models mentioned. For 6000 crushes, we end up with the profile as seen in loop.6e3. Here we see the jamming tendency already in the 1 - 10u range (cpd=90%) and we have our cumulated % distribution of 99.9% within 1 - 100u range.

For each of the 6000 crushes, we calculate the crush energy required for each new change of % particle distribution. So we have crush=sqrt(max/id) for each diameter max crushed to diameter id as explained above where diameter max replaces 60u. We have id in the denominator because we crush to the small sizes being in the denominator < numerator so crush > 1. Now we accumulate crush by ccrush. However, crush is really weighted by the number of particles np or pd, really pd/100 since pd is given by %'s. So the cumulated crush would be ccrush = ccrush + crush*(pd/100), but we're applying the % of the normal erf curve (pdn) representing how the crushed particles are distributed. Here, pdn @ max/2 would be the highest. So now, ccrush = ccrush + crush*(pd/100)*pdn when pdn is in decimal form and pd in %. For each of the 6000 cumulative crush intervals (ccrush), we have a unique particle % distribution. Both the ccrush and associated new % (pd) are written as a 3 character binary file with ccrush for each of the 449 % pd. This file (cumcaf.pd) is over 8 mb in size.

        C. Matching Crush Energy Between Files

So we read the total dust from one file, then the WTC rcaf,tonbcp - dust produced by each slab tonbcp (cton) from another file. Now real caf (rcaf) above is defined as poton/tonbcp, or a factor of how many times tonbcp could be crushed to the quarry profile % dist (crush.%d). Now in review from above, poton (potential tons) = sle/crushp = tonbcp from ton (lesser of tbs=tons between slabs, tps=tons/slab). if(tonbcp > 2*ton) tonbcp=2*ton. It's 2*ton because tons from slab & pile are crushing each other equally. In topic above, we just calculated the cumulated crushing adjustment factor (ccrush). So now we match these energy factors, after we rewrite the WTC rcaf,tonbcp file in ascending order ( to match the ascending order of the 8 mb cumcaf.pd file. We read the file first, then ccrush in the cumcaf.pd - %pd file. As soon as we match the real rcaf by the code if(ccrush > rcaf), we continue the id loop to read the np particles % distribution (pd). An example of this is seen in slabs.dat. The search for the ccrush > slab's rcaf is given for the 1st and last 11 slabs. Note, we start the loop 1 back from the previous matched slab in the loop.

Now, we matched ccrush with slab's rcaf so we have a particles % distribution stored for each slab in an array of particles dimension pact[id] = pact[id] + pd*tonbcp. We used pact[] again, which originally represented particle activity but now used here as tons of dust for each particle size diameter accumulated for each slab. Actually, tons - volumes can be used interchangeably by a constant of density so it actually is still particle activity pact[] again but expressed in terms of tons this time but really means ton[i]. As a check, we now totaled the final pact[]'s, from all collapsed slabs, and the crushed tons cton = cton + tonbcp from the file = Crush Tons = 59455.6. The sum of pact[i]'s was 59525.0 or .12% too high. With rounding off and accuracy, this is excellent agreement, however all pact[i]'s were scaled by 1/1.0012 = .9988 correction factor anyway.

        D. Recombining Rub and Crush Dust % Dist

Finally, we have the total crushed tonbcp's (cton) as dust from slabs. We also read the total tons of dust (tton) from a file earlier, then rton=tton-cton is the total rubbed tons of dust. From our previous model runs, we have Total Tons 73788.063 = Crush Tons 59455.637 + Rub Tons 14332.426. In the rubbed % part dist file (rub.%d), the pd read after each particle diameter i produces rubt=pd/100*rton of dust. This added to the crushed dust ton[i] or denoted by pact[i], the new % of dust for each particle i is dustp=(pact[i]+rubt)/tton*100. No need to work with an array [i] any more, so i and dustp are dumped to the final dust output file dust.%d. However, we did not have any rubbed particles > 399u, so the rest of the tons/particle i diameter is the crushed tons only from 400 -> 449u. To see a comparison of rubbed and crushed tons tabulated for each particle i = id diameter, see crush.rub. Note, the crushed tons < rub tons from 1-13u then crushed tons > rub tons from >= 14u.

Using dust.%d in the DUSTBM.c again, and removing the fluff factor because of rethinking the packing of tiny particles layered to a few inches, compared to pulverized limestone, we have the updated output file DUSTBM.dat. The rubble pile has not changed nor the amount of dust but because we have overall larger particles in crushing than rubbing. As seen in crush.dia, we have a best fit with the observed WTC street level dust, so we repeat here, experimenting, we found that 3.6 x slab thickness rather than 2, and the occurrences of no lateral movement is at 1/1.4 times the slab travel distance of 3.6 times it's thickness. So 1/3 --> 1/1.4 and 2 --> 3.6 correlates better with the observed WTC street dust photos from the web. So finally for the WTC's, The max dust was 2.57" @ .147 mi and .45" @ 2 mi with 52.5% of dust deposited within the 2 mile range. 62.9% of WTC buildings went to dust. - End Update 10/22/08

8. Kinetic Energy Effects of a Nuclear Blast:

        A. Pressure Wave

After viewing 1 Megaton PSI Effects, we chose within a radius of 5 psi from a any nuclear blast, to collapse all buildings, at least the cement. All the calculated dust will be combined with the ground crater dust from a surface blast. Each component has it's own particle distribution to be inhaled by the updraft velocity of the stem in the nuclear mushroom cloud. We plotted each of their radial distance x vs. psi. For quick reference, it was as follows:

R mi PSI
1.7 12
2.7 5
4.7 2
7.4 1

Most power laws have energy decreasing by 1/R^2. A plot of this is seen in log log plot rpsi.gif. Here the eyeball straight line graph, rather than the best fit least squares method with such a few data points, takes on the form, y=ax^b where y=psi and x = distance in miles. The real guess is to interpolate the value a at x=1, where y=a. Best guess was about 28.5 PSI, when plotted on more detailed log-log paper. Now for slope b, we could do a dy/dx on a print out here and that would be ok. However, using a simple graphic arts program (MS Paint), and noting that y-pixels increases from the top, we have b=dyp/dxp = (141-455)/(642-72)= -314/570 = -.5509. However note, that they expanded the x scale to show more detail. How much, 10v = 219 pixs, 10h = 655 pixs. So 655/219 = 2.99 so they must have used 3. Now -314/(570/3) = -1.653.

From the graph, the mid points 2 & 3 seem to really parallel the line, so we use them for slope b. We have b = dy/dx, but with logs, b = ln dy/ln dx. Then dx=x2-x1 and dy=y2-y1 with b=ln(y2-y1)/ln(x2-x1) = ln(y2/y1)/ln(x2/x1). Now subbing, b=ln(5/2)/ln(2.7/4.7)=ln(2.5)/ln(.574) = .9163/-.5543 = -1.653 ditto matching the visual slope by using dx & dy pixels. Finally Y=28.5 X^-1.653 or psi=28.5 mi^-1.653. Note, the ideal power law would have b = -2 or x-^2 = 1/x^2, the inverse power law. Since 1.653 < 2, the rate of psi decreases with less than the expected 1/x^2. Thus the psi wave reflections from the ground overcompensates the loss by absorption through the ground.

However, we want psi=5 and we need the equation for R mi. We know from the table that we have 2.7 mi @ 5 psi, but by the formula, y=a x^b, x^b=y/a. Take 1/b root of both sides, x^(b/b) = (y/a)^1/b, x^1 = x = (y/a)^1/b or mi = (psi/a)^1/b. Subbing, mi=(5/28.5)^(1/-1.653) =.1754^-.605 = 2.87 mi > 2.7 mi table value or a 6% error here. For other size nuclear blasts, we still have the same decay rate of energy (b), or slope of the line. So all other bomb sizes will have parallel lines but their intercept a @ x=1 is different. If the nuclear explosion is doubled, the intercept a is doubled. That is a=2*28.5 = 57. So a 2 megaton bomb has a psi of 57 at 1 mi, a 1/2 megaton bomb will have a psi of 14.25 at 1 mi or expressed as a=mega*28.5. Then back to mi=(psi/a)^1/b = [5/(mega*28.5)]^(1/-1.653) = (.1754/mega)^(-.605). Then corresponding to y=a x^b is psi=mega*28.5 mi^-1.653. In rpsi.gif, we show parallel plots of mega=3.5 and mega=1.

        B. Crater Dimensions

For a surface burst of 1 megaton, the PBS website indicates a hole of 200' deep, 1000' wide. This is seen but inverted in crater.dia. The crater is the top segment of the sphere (hole in the ground) with a radius of curvature R. The width of 1000' is broken into the smaller radii r, which is the radius of the hole at the ground surface - remember, this is inverted. This dirt segment volume of the sphere is also called the cap of a sphere. We can compute this dirt volume by either of the 2 formulas shown v=pi/3*h^2*(3*R-h) that requires the radius of curvature R=(r^2+h^2)/2h or v=pi/6*h*(3*r^2 + h^2) using r and h directly. The radius of curvature here has it's center 725' above the bottom of the hole or 525' above the ground. In both cases, v = 8.273e7'^3 and the volume V of the whole sphere, of which the segment is a part is, V = 1.596e9'^3. This means the dirt blown out the hole is only 5.2% of the energy used as if the bomb were implanted below the surface deep enough before the explosion, but near the surface enough to just blow the hole out in the ground above it. This 5.2% will hold for all bomb sizes so we simply multiply v = 8.273e7'^3/mega x megatons for any other surface burst bombs. However, we don't have the particle size % distribution yet. Now we could take dirt, shake it up and blow it onto the plate glass for scanning as we did before with the limestone from the quarry, however, this is certainly what a nuclear blast does. We may not have to do this as we explain later.

9. Extrapolate WTC Dust Generation to Major Cities: - CITY%D.c

        A. City Building Profiles

We now have the amount of dust for each particle size from 1-449u from the combination of rubbing and crushing. We will apply this method to different cities as we estimate each city's building density distribution from photos. If any reader has done this with better information, please contribute to this model development as we will acknowledge all those that contribute.

All the buildings within the 5 psi ring will use some of the other building parameters from the table above rather than the WTC's data. As a first case, since we did WTC, we might as well start with New York City (NYC). Looking at a map, we classify Manhattan. We will have approximate city profile by building estimates on each city in a file called city.dat. Such a line of data is as follows: NYC cityl= 9 cityw= 2 cityd= 65 floors1= 65 floors2= 25 length= 240 width= 225. Now cityl=city length in mi, cityw=city width in mi, cityd = city building density in %, floors1 = highest average # of floors (stories) in center of city, floors2 = average # of floors (stories) at end of city (cityl), length = average building length in ft throughout city, and width = average building width in ft throughout city. Past the length of the city, we have floors=4 stories.

However, using this city's building height profile and size, we need to estimate the height for every .1 mile from ground zero to the city's outer limit shown by the comma symbol with city's equivalent area radius Rc in citydust.dia. So for each dr=.1 mi, we round off to the nearest number of buildings floors, rather than the building height because the distance between floors is assumed to be divisible by the slab thickness. For the WTC, it was 12.4 ft, and rule of thumb is 10' floor-ceiling. We will center between 10' & 12.4' to near 11' including the slab. We can run the model and observe the effects of different slab thicknesses. For whole # slabs, nsc=11/slabt+.5. Then the actual ceiling height ch=slabt*nsc. We loop up the building by the slab thicknesses, so the number of times to loop n is ln=nsc*floors. With every nsc, mod = remainder n/nsc = 0, we evaluate the slab crushing routine. We linearly interpolate and round off the # of floors between floors1 -> floors2 to distance x5psi. Any distance past Rc, floors=4 stories. We have lx=Rc/dr+1 buildings to evaluate along the radial, since we start at center of city ix=0, to ix<=lx at x5psi.

        B. How Many Building Fall?

There are 3 possible cases where we must examine the 5 psi ring for city and metro-suburbs. These are presented in citydust.dia. In the top diagram, this is the most likely case, the 5 psi ring is greater than the city limits and encloses it. All the city buildings (#) are gone, as well as some of the suburbs (^) buildings out to the 5 psi ring. In the 2nd case, we reverse the 1st case, the city limits are bigger than the 5 psi ring. The 5 psi ring never reaches to the end of the city buildings or any of the suburbs. Many city buildings and all the suburbs buildings are still standing. In both of these cases, we can use the city building density cityd without modification because all the city lies within the 5 psi radius or the 5 psi radius is a circle that includes only a part of the city.

Our 3rd case is complex in that the 5 psi ring criss crosses the city's boundary. If we use an equivalent city radius Rc, this won't criss cross but will be either in or out of the circle boundary as either of the 2 cases above. This situation is likened to NYC's Manhattan's Island, where we have a 2 mi x 9 mi dimension. It will be best with such an elongated dimension to perform an alternative method. By the diagram, we see in reality, the length of the city is not reached by the 5 psi ring so that forces the buildings in the city's length into it's width which over estimates because for NYC, much of it is water, and should not be included. If the city's dimensions are almost a square, use of the equivalent radius Rc is fine.

We need to know, the city area inside the 5 psi ring. So by integration, we first use the equation of a circle and solve for y in x^2 + y^2 = a^2, y = (a^2 - x^2)^1/2 where a = radius 5 psi ring (x5psi). Now the x limits are from x1=0, to x2=cityw/2 for right upper quad the of circle. So INT y dx = [x*sqrt(a^2 - x^2) + a^2*asin(x/a)]/2, remove /2 and we have right semicircle since % in area of quad = % in area of semicircle and circle. Substituting limits x1 and x2 for x, we have NYC's 5psi area mi^2 cpsi5a = x2*sqrt(a^2 - x2^2) + a^2*asin(x2/a) - (x1*sqrt(a^2 - x1^2) + a^2*asin(x1/a)). Subbing with x1=0, we are left with cpsi5a = x2*sqrt(a^2 - x2^2) + a^2*asin(x2/a). Now the area of the 5 psi ring is psi5a=pi*x5psi^2 so the ratio of the integral area to it is rat=cpsi5a/(psi5a/2), using /2 for a semicircle. Now we can reduce the cityd to rat*cityd.

Back to the 5 psi ring (psi5a), we convert mi^2 to feet^2 and divide by the area of the building (width*length), we would cram all the buildings into the 5 psi ring or 100% and have the max # of buildings (nbm). When we multiply this, nba = nbm*cityd/100, we now have the actual # of buildings (nba) within the 5 psi ring. In the each building model run loop, we have the tton of dust but that was for lx = 28 buildings along the radial. So the average for 1 building is tton/lx. Now times nba actual buildings, we have total tons of dust for the city. However, this is only a fair estimate.

            (1). Weighting by Concentric Rings Areas - Update 11/20/08

We use dr=.1 mi out to the 5 psi ring, we have 28 buildings dr's. There is a building dust calculation for each dr. For any 2 radii, the area of a ring is the difference between their areas or pi (r2^2-r1^2) or pi (r2+r1) (r2-r1). If dr is the thickness of a ring, and r is the inner radius, then r1=r, r2=r+dr. So ring area = pi (r+dr+r) (r+dr-r) = pi (2r+dr) dr. Or if r2=r and r1=r-dr, area = pi (r2+r1) (r2-r1) = pi (r+r-dr) (r-(r-dr) = pi (2r-dr) dr. In the 1st case, r can start the loop at r=0, but not in the 2nd case so we would start at r2 = dr, then r1=0 for the start of the loop. We actually took this case as we didn't want to start with r=0 since it may have caused problems later in using r later in the program loops. For the 28 dr ring buildings, we see the results in Here we see the Dust Ring Area and it's accumulated columns. The decimal or Fraction column is the cumulated areas ratio to the 5psi area (psi5a) up to 1.000.

            (2). Another Look at PSI's

After some subsequent runs of the later programs, to cut off the city's building dust at the 5 psi ring seemed too short and abrupt. In the real world it, would rather be a diminishing of dust . Reviewing 1 Megaton PSI Effects again, less buildings will fall at 2psi so let's cut off any buildings falling at 1 psi instead. Using our formula above to work with in the programs, instead of the table of points, we use x5psi=2.76, x2psi=4.69 and x1psi=7.00 mi. Now for our new radial distance r, how do we reduce the # of buildings falling between x5psi and x1psi? We could linear interpolate 100% at x5psi and 0% at x1psi. However above, we saw it as a power reduction formula as seen in rpsi.gif as psi=mega*28.5 mi^-1.653 or a slightly better fit as psi=mega*29 mi^-1.73. Here mega=1 and we just solved for mi for x5psi,x2psi and x1psi above. We introduce a new ratio, pr as the Power Ratio or PSI Ratio as pr=(x2psi/x5psi)^-1.73 = 0.4000 as an example. So the ratio of these distances is still that inverse power law, but not quite the inverse power law squared (-2). Subbing for the rest for mi, x1psi/x2psi = 0.5000 and x1psi/x5psi = 0.2000. If just the ratios of distances, we have 2.76/4.69=.588 instead of .4, 4.69/7=.67 instead of .5, and 2.76/7= .394 instead of .2 indicating the power law as a better and faster reduction of the pressure wave with distance. Notice that the simpler form is just the ratios of the psi's, psi2/psi5 = 2/5 = .4 yields same results. However, we need the psi with each dr so pr = psi/5 = mega*29/5 mi^-1.73 where mi=r. We will use the pr to represent the reduction factor for the # of buildings to fall, but for x=r>x1psi, we have 0.0. We rather drop pr from .2 to .0 past x1psi rather than 1.0 to .0 past x5psi. For r <= x5psi, we keep pr=1.

            (3). A 70 Ring Circus

Where do we apply pr? Rather than put it off and apply it's correction in some of the subsequent programs and maybe make a mistake or lose it's meaning in the later complexity, we thought it best to apply it directly as a correction factor now to the Dust Ring Area that is stored in an array called perimeter area pa[] for use in the following major nested loop. Also, since we now use x1psi instead of x5psi, with dr=.1 mi out to the 1 psi ring (7.00), we now have 70 dr rings or 70 buildings. Now instead of working with the max # of buildings (nbm) for the whole 1psi area, we compute nbm and adjust to nba (actual) by the city's density for each ring area pa[ix] and now further adjust pa[ix] by the factor pr. The dust components and new ring areas with pr applied is seen in city.%pd. Notice here, that the Dust Ring Area = pa[ix] until pa[ix] is corrected by the pr < 1 factor. Also notice the dust components and the crater dust volume comparison. Not only do we have a 70 concentric ring circus, but we have more smaller performers in the larger outer rings.

            C. Super Duper Looper I - CITY%PD.c

The 1st file to read is city.dat - the city's building profile. The 2nd file to read is rub1.NYC which is the Rubbed Tons (rubt) for each of the 70 ix (dr=.1mi) ring buildings from CITY%D.c . Since we have assumed, but not exactly, to distinguish between the two dust generation methods, we have the same particle distributions from all rubbed slabs. That is, more dust, not more smaller particles as crushing does. We open the 3rd file rub.%d which is the Rubbed % Particle Distribution (pd) up to 399u. With these last two files, we initialize the tons of dust for each particle size id as ton[id]=rubt*pd/100*nba*pa[ix]. The pd is the % distribution of the particle size id and the nba*pa[ix] is the real # of buildings per ring [ix] as explained in the previous paragraph. Still reading the 1st ix ring data, we now open the 4th file crush1r.NYC, which is the dust potential from the crushing slabs program CITY%D.c reading each slab's rcaf (column 2) and tonbcp (column 3) from all 70 ix ring buildings (column 1). In crush1r.NYC, we have an indicator record (1.000 0.001) when we are reading a new stack of slabs in the next ix building ring. We now open the 5th file, cumcaf.pd , the binary file from CRUSH%D.c and run another nested loop to match 'ccrush' @ (id=0) with rcaf as explained above. Since there is no id=0, it's a perfect place to identify 'ccrush'. With pd = particle size (id) decimal distribution for each slab in building ix, we now add this to the rubbed tons id array ton[id] = ton[id] + tonbcp*pd/100*nba*pa[ix], which is the same procedure as for rubbed tons except tonbcp replaces rubt.

As the first ring of building's (ix=0 -> .1mi) is complete, we write it out to ctydst1.NYC, as the last ix (lix) and the ton[id], and reset ton[id]=0 for the next ix. We then repeat all these loops by continuing to read the next ix ring and must rewind files rub.%d and cumcaf.pd again for particle size distributions. Finally, when all the buildings are read, we reread all the ix and tons from ctydst1.NYC for a check to see if it matches the totals before all the pd's and pa[ix]'s sub divided the total against the 6th file totald1.NYC . This difference was < .05% error so no errors in the loops or the summated pd factors were evident.

            D. Super Duper Looper II CDRCRD.c & PD%R.c Update 03/08/12

Particle Distribution % by Regression replaces a few loops in CITY%PD.c but is now further complicated by the list of discrete buildings with their random radii (ir) from the LAS format instead of generated consecutive looped building previously. Since we have discrete buildings, we would like to store the % pd x tons 'tons' by radial distance 'ir' & particle diameter 'id' depicting a 3D surface, for a final output file but is too large for our CPU ram here. We output to a file then read it by file positions generated by id,ir. Then we have to do the same over again for the crushed dust after the rubbed dust. Most all large files in these RADFO programs use this procedure.

We finally approach this problem because Super Duper Looper I above was complex and was developed to combine all these files in CITY%PD.c rather than all the effort in CITY%D.c instead. A two-program approach was to make this easier. However, we are forced to a new method in CDRCRD.c because we now evaluate different building dimensions in the same psi ring (ir) and the number of buildings create their own actual density spacing and height instead of estimating it.

The matching of all the files in CITY%PD.c was done by sequencing through until a match was found. Matching 'rcaf' with 'crush' could apply to 6000 records for every building vertically for every slab thickness. In the estimated building height profile, linear decreasing with distance (ir), we could estimate where to start searching in the files since it is a few feet lower than the last building. We could also reverse the file order and start with the smaller buildings and the next larger building's matches should be the next slab group into the file. Here, any building height appears in any 'ir' ring. Therefore, we may have excess CPU times if we continue this method here.

Once the match is found, we extract all 1->449 particle id's pd's (% distribution) in the binary file's (displayed here as ASCII) seen to the right of id=0 column in cumcaf.pd which was calculated from CRUSH%D.c . PD%R.c is just a modified CITY%D.c from above. Instead of incrementing each ir psi ring and approximate a common building height vs. distance ir with fixed dimension for each building, then estimating the number of these buildings in this ir psi ring by a density guess, we calculated each ir psi ring by vector range from each building to the city's ground zero building as well as their real dimensions and even their orientation for psi impacts. All this will be developed later but now seen in 920515dr.lwh as a later input file.

                (1) New Files & Order

Now the new type of city data (1st file) to read is Miami.dat - the Ground Zero's building name & location by 3 different coordinate systems (Lat/Lon, UTM & State Plane). The 2nd file to read is rub1.MIA which is the Rubbed Tons (rubt) for each of the 70-ir (dr=.1mi) ring buildings from CDRCRD.c . Since we have assumed, but not exactly, to distinguish between the two dust generation methods, we have the same particle distributions from all wall-pile rubs. That is, more dust, not more smaller particles as crushing does. We open the 3rd file rub.%d which is the Rubbed % Particle Distribution (pd) With these last two rub files (tons/ir & % dist/id), we have the % tons of dust for each particle size id up to 399u for each given radial increment 'ir' to complete our 3D surface area.

Still reading the 1st 'ir' ring data, we now open the revised 4th file loop crush1.MIA which is the dust potential from crushing slabs program CDRCRD.c reading each slab's rcaf (column 2) and tonbcp (column 3) from all 70-ir ring buildings (column 1). Notice each building's 'ir' is in random order but remains the same until vertically profiled by slab thicknesses.

We now open up that same 5th file cumcaf.pd , the binary file from, CRUSH%D.c but run in our new program PD%R.c that replaces the end portion of CITY%PD.c . Now we develop a regression formula to match it's 'crush' with 'rcaf' directly.

        E. Quick Match Reference?

As seen by the first id=0, which is a diameter of 0 in CRUSH%D.c , is the ccrush term. Ignoring the rest of the file record, which is the 449 particles % distribution, we converted the 3-char ccrush term to ASCII as seen in rcumcaf.b2 . Note that the interval for 'ccrush' changes and sometimes erratic near the bottom. The test for this is to plot ccrush for each 'n' and try a linear regression analysis. Note that we are trying to match the 6000 record #'s (y) from the match of ccrush <-> rcaf (x). We ignore the first record and handle this later. Our first plot is seen in RCAF # 1.

                (1) Regression Analysis

The 6000 rcaf's (x) plot vs. the record # (y) is seen in blue and the best fit linear line is seen in yellow. Actually, it's non-linear but appears linear. Since the raw data's blue line is non-linear, we modify the regression equation's to handle the non-linear form y=ax^b instead of the linear y=bx+a. We need a working formula, that is, we don't have to go back and read each x,y again after we have the mean of x & y. This is developed in reg.dia.

                (2) Regression Analysis Results

The results are b=0.914 & a=144.421 with correlation coefficient r=.976. Actually, the term r^2 = .953 is more important meaning that 95.3% of the observed y variation can be predicted by the regression equation, or 4.7% is random noise? Well this seems pretty good to use, however, suppose we try to tweak and fudge it for better results. Actually, we could have written a program to loop a & b until the 'r' becomes better.

                (3) A Bad Surprise

Before this, we want to test the fit for other known equations so we generate x,y by a formula and then try to see if the non-linear regression discovers a,b exactly. Generating 150 points for y=10 x^b, where we have b from 1.0 -> .2, we see in axeb.txt the diminishing results of r with b. Not bad but thought we would have an exact match instead. So the conclusion is that with a set of related data, we may have a best fit line but not perfect on exact (r=1.000) correlated data.

                (4) Play Graphics

What happens if we force the a,b and plot the curves to see if the equation matches the raw data. It may be better to see what's happening rather than just watch the end result of r. Actually, it produced stunning results. This is time consuming but we only need a,b for one time. Now with the new a,b, we have the best visual fit in RCAF # 2. This is almost a perfect fit but we didn't calculate r but b=.678 & a=325. However, another statistical term called the standard of estimate (S.E.E) went from 538.8 down to 118.3.

                (5) Play Calculus

Well the black gaps between the yellow and blue lines in RCAF # 2 indicate where the data isn't an exact fit. Can we improve on this? Well we have to match rcaf with ccrush. Using y=a*x^b with b=0.678 & a=325, we find y (rec #) to open to read ccrush. If it isn't a match, we find another ccrush, say 10 records further by file position filp2=filp1+10*rl where we just read ccrush1 with filp1 using record length 'rl'. So the delta ccrush is dc=(ccrush2-ccrush1)/10 so then the number of records to match rcaf is delta records dr=(rcaf-ccrush1)/dc. These results are seen in rcaf2.b1 , rcaf2.b2 and rcaf2.b3. The last column tells what % error the final ccrush from rcaf occurs. The real absolute average difference is given by 'dif avg' on the last line. The last ( rcaf2.b3) is probably the most important range of rcaf's used which average less than a 1% error.

Could we still improve on this? We have calculated a dy/dx slope near the ccrush by using 10 records above to find a closer approximation to rcaf. Since the data has some noise and we may have jumpy slopes along the raw data, suppose we take the slope of the non-linear curve y=a*x^b. This is a classic basic derivative theorem in Calculus. So with y' = dy/dx = dy(ax^b)/dx = a*dy(x^b)/dx = ab*x^(b-1), we have the slope of the best fitted curve at ccrush1. Note this is displayed in rcaf.dia. The choice of this method also smoothes some erratic data. So we have y'= 325*.678*x^(.678-1). Applying this techniques, we have rcaf4.b1 , rcaf4.b2 and rcaf4.b3. Even though the 1st 2 appear worse than the other method above, the last is the important range used since the higher values with small rcaf values calculate to small tons of dust. We are satisfied by using this method which should produce < .4% error in the most weighted used rcaf range.

With pd = particle size (id) % by decimal distribution for each slab thickness in each building's ir, we now add this to the rubbed tons' diameter ir array ton[id]. This means we will initialize the final output file 'ctydst1.NYC' for NYC's example with 'ir' range vertically 0 -> 70 (dr=.1)'s and tons/id diameter's horizontally 0-449u. Here, we don't know what building will be the first and last ring 'ir' so we just keep summating the file array with it's % pr until complete. We have no more loops by continuing to read the next 'ir' ring and rewind the files rub.%d and cumcaf.pd for matching then write out particle size distributions. The final tons/id/ir/ is seen in ctydst1.NYC. Notice, the variable tons value length from Super Dooper Looper I. However, in our new file PD%R.c , we need a fixed value length to produce a constant record length to use file positions 'filp' to read, summate and write again - such as ctydst1.MIA. We also changed a few variable names and defined them in a glossary with a better pic in Trap4 Modified.dia . - End Update 03/08/12

        E. PSI Wave => Buildings: - 02/16/11

        F. Crater vs. Building Dust Results

We convert the tons of city's building to volume and compare it to the volume of the crater. Usually the density of dirt is considered at 2.5 g/cm^3, and we are using cement as 2.4, since we have a mixture of road and sidewalk, they are so close already, let's consider them equal, which we should do in the radioactive fallout model to really save CPU time. Running the model for NYC using 1,5 & 10 megaton bombs, we have the tally below:

                1 Mega
City Buildings: 1.197e+8 Crushed Tons of Dust = 1.597e+9'^3 = 19.31 x Crater Vol (8.273e+7'^3)
City Buildings: 1.415e+8 Crushed + Rubbed Tons of Dust = 1.889e+9'^3 = 22.84 x Crater Vol (8.273e+7'^3)

                5 Mega
City Buildings: 5.461e+8 Crushed Tons of Dust = 7.289e+9'^3 = 17.62 x Crater Vol (4.136e+8'^3)
City Buildings: 1.275e+9 Crushed + Rubbed Tons of Dust = 1.702e+10'^3 = 41.16 x Crater Vol (4.136e+8'^3)

                10 Mega
City Buildings: 9.469e+8 Crushed Tons of Dust = 1.264e+10'^3 = 15.28 x Crater Vol (8.273e+8'^3)
City Buildings: 3.606e+9 Crushed + Rubbed Tons of Dust = 4.814e+10'^3 = 58.19 x Crater Vol (8.273e+8'^3)

From 1-5 megatons, if we don't include all the city buildings, the ratio is smaller. If we include all the city buildings, the ratio is higher but if we include many of the city's outskirt smaller buildings, it is lower again because we weighted more smaller buildings while the dust from the crater hole is proportional to the size of the megaton blast within a large city's linear building profile. We noticed that the larger the megaton blast, the building dust increases faster than the crater hole dust increasing the ratio to 50+ times as much for 10 megatons. - End Update 11/20/08

Experimenting with the different slab and wall thicknesses produces somewhat complex offsetting effects. It is difficult to pinpoint what effects each have on the total dust. Finally, we decided to go with the 4" slab thickness, as many have quoted for older buildings combining with even thicker slabs with pillar supports. We averaged this with the lightweight strong cement of 2.5" The wall was considered a constant of 12" thick. These are the ground rules we will apply to all the cities. The % to dust for a 65 story building was 38.0%, and for a 20 story building was 35.6%.

10. Combining Crater and Building Dust into Mushroom Cloud:

For the crater dust, it is obliterated and immediately sucked into the nuclear cloud plasma and vaporized into the near smallest observed dust particle of 1u so the % particle distribution of the crater dust doesn't matter as it will all be 1u when passing through the plasma. Now, also for the building dust, much of it will also be sucked into the stem of the nuclear mushroom cloud and pass through the plasma reducing it's size to 1u also. The question is, how much of the building dust will be swept up?

In order to answer this, let's assume a minimum horizontal wind speed for the largest of particles (450u). The derived fall velocity formula is complex and will be referenced more when developing the RADFO model components. For reference, the calculated terminal velocity of 1.3564e+1 fps or 9.25 mph for <= 450u to 2.5082e-2 fps or .02 mph for <= 10u. How far out from ground zero is the inrush of air?

11. Radioactive Fallout Model: