能量估计

能量估计
一、问题重述
1945年7月16日，美国科学家在新墨西哥州Los Alamos 沙漠试爆了世界上第一颗原子弹，这一事件令全球震惊. 但在当时有关原子弹爆炸的任何资料都是保密的，而很多其他国家的科学家非常想知道这次爆炸的威力有多大.
两年之后，美国政府首次公开了这次爆炸的录像带，而其他数据和资料仍然不被外界所知. 英国物理学家G. I. Taylor (1886 ~ 1975)通过研究原子弹爆炸的录像带，建立数学模型对爆炸所释放出的能量进行了估计，得到估计值与若干年后正式公布的爆炸能量21 kt 相当接近。

Taylor 是如何根据爆炸录像估计的呢？主要是通过测量爆炸形成的“蘑菇云”半径来进行估计的. 因为爆炸产生的冲击波从中心点向外传播，爆炸的能量越大，在相同时间内冲击波传播得越远、蘑菇云的半径就越大. Taylor 通过研究录像带，测量了从爆炸开始的不同时刻t 所对应的蘑菇云半径r(t),如下表所示：
表1 时刻t(ms)所对应的“蘑菇云”半径r(m)
t r(t) t r(t) t r(t) t r(t) t r(t) 0.10 11.1 0.80 34.2 1.50 44.4 3.53 61.1 15.0 106.5 0.24 19.9 0.94 36.3 1.65 46.0 3.80 62.9 25.0 130.0 0.38 25.4 1.08 38.9 1.79 46.9 4.07 64.3 34.0 145.0 0.52 28.8 1.22 41.0 1.93 48.7 4.34 65.6 53.0 175.0 0.66 31.9
1.36 4
2.8
3.26 59.0
4.61 67.3
62.0 185.0
二、预备知识：
原子弹爆炸所释放的能量估计将涉及很多因素，在这种复杂的情况下采用量纲分析法将极大的简化问题，其中涉及线性代数中矩阵的计算，以及用最小二乘法对数据进行拟合求相关的系数。

对于MATLAB 等数学软件的掌握也是必不可少的，能起到事半功倍的效果。

三、问题假设
1 同一时间只有一点发生爆炸，传播的空间没有大型障碍物阻止。

2 爆炸开始的时间定为在t=0,爆炸的能量完全释放。

四、模型建立及求解
1 考虑到原子弹爆炸在极短的时间内释放出巨大的能量，蘑菇云半径r 主要与时间t 、爆炸能量E 、以及空气密度ρ等几个参数有关，要寻求的关系式记作：0),,,(=ψt E r ρ。

2 这是一个力学问题，基本量纲选为L ，T ，M 。

上述各物理量的量纲表为：
M
L MT L E T t L r 32
2
1][][][][---====ρ
其中n=3<m=4。

3 由以上各物理量的量纲表可写出量纲矩阵：
ρ
(((()()()(02-1-011
003-2014
3E t r
T M L ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=A ⨯
并且计算
RANK （A ）=3
4 解齐次方程 Ay =0
方程有m-r=4-3=1个基本解，可取为：
[]T
y 11251-=
5 给出这一个独立的无量纲量
ρπ1251-=E t r
6 为了更为准确的计算爆炸能量，将蘑菇云半径公式改写为：
t r E
ln 2ln 5)ln()ln(1+=+πρ
在实验室通过做实验，测定E ，ρ,r,t 后可计算得524.10ln 1-=π。

此时可根据测量数据得到t r ln 2ln 5+对应的数据的平均值c= 21.5628,由1
ln -πρc e E ≈，由
此可计算得J E 13
106123.8⨯≈，单位换算后约等于20.65kt 。

五、解及模型分析与应用推广
首先改写蘑菇云半径的公式为b
at r =的形式，而要作线性最小二乘法拟合，进一步改写公式为：t b a r ln ln ln +=。

根据测量数据我们得到lnr 和lnt 的数据，将它们的函数关系拟合为一次多项式，得到系数b=0.4058,其值与前面分析的结果2/5非常接近，从而验证了量纲分析得到的公式的正确性。

此模型在和平年代可以用于军事演习以及房屋桥梁等的爆破上，便于领导在不懂得专业知识的时候依然能估算炸药用量，以防下属虚报数字，贪污国家或公司财产。

MATLAB源程序：
% 读入数据t和对应的蘑菇云半径r
[t,r]=textread('c:/energy.txt','%f %f',-1); % 以下将R,T,E,p用基本量纲L,M,T表示
clc
R=[1 0 0]';
T=[0 0 -1]';
E=[2 1 -2]';
p=[-3 1 0]';
% 创建A矩阵
A=[R T E p]
% 确定方程基本解
r0=rank(A)
y=null(A,'r0');
y=y./y(4,1) %将基本解的系数化为整数
% 拟合求解E
x1=5*log(r);
y1=2*log(t);
c=sum(x1+y1)/length(t)
E=exp(c+10.524)
% 解及模型分析
y2=log(r);
x2=log(t);
polyfit(x2,y2,1)
程序结果：
A =
1 0
2 -3
0 0 1 1
0 -1 -2 0
r0 =
3
y =
5.0000
2.0000
-1.0000
1.0000
c =
21.5628
E =
8.6123e+13
ans =
0.4058 3.5903。