主题：[原创]杭州70码(欺实马)的验证

杭州70码(欺实马)的验证

2009年5月7日晚20时, 杭州市发生了一具飙车撞死行人惨案,  杭州师范大學体育系的学生胡斌驾驶一辆改裝三菱跑車在在路上行走的浙大毕业生谭卓撞死.
杭州交警召开新闻发布会，提及“当时车速在70码”，由此引发舆论不满.

5月14日，杭州飙車案事故鉴定完成，专家称“车速肯定不是70码”。杭州市公安局当日向媒体发布交通肇事案鉴定报告，认定事故车在事发路段的行车时速在84.1公里到101.2公里之间.
对这个数据, 很多人仍然持怀疑态度, 见到当时碰撞场面的目击者说, 车速至少在100码以上.  那么真相究竟如何? 我们一起来验证一下.

有一条重要的消息是, “有目击者声称，谭卓被撞出大约4-5米高后再重重摔在20米以外的地方，当场死亡”。这条信息非常重要, 目击者对这种基本情况的描述, 比较有可信性. 根据这一点数据, 已经可以对当时事发时的车速进行估算.

肇事者当时驾驶的是三菱的EVO 9型(Mitsubishi Lancer Evolution IX)改装跑车, 从网上查到这种车的自重是1490 kG. 加上驾驶员,及改装后的重量, 行驶时总重量应超过1600 KG, 该车的极速是249.4 kM/h, 最大输出功率213.3 kW, 车体宽1770mm, 高1450mm.

汽车与人体相撞时,发生的是非完全弹性碰撞, 如此被撞的人体才会飞出去, 因为人体不是一个刚性球,像乒乓球那样,一拍就跳起来,也不是一块烂泥, 那样的话汽车撞上以后人就粘车上了对不对.
如果用恢复系数e来表示,  e=( V2-V1)/(V10-V20), 其中,

V10 表示碰撞瞬间的车速,V20表示行人碰撞前的速度,  V1表示碰撞发生后的车速, V2表示碰撞以后行人获得的速度. 那么在完全弹性碰撞时, e=1, 在完全非弹性碰撞时, e=0.
我们知道人体是由骨头,肌肉,皮肤,组成,外面还有衣服, 所以并不是太容易受伤害的, 比如平时别人打你一拳, 肌肉的弹性可以消耗掉大多动能. 在非常硬的物体与人体发生碰撞时, 这种碰撞可以认为是完全非弹性碰撞, 可认为车撞人的时候, 人对车完全没有阻止作用, 则e=0, 应该是比较合理的数值.


谭卓在行走时,速度约为5 kM/h,正对车而行, 则V20=-1.3888 M/s  同时车与人碰撞时系统动量守恒, 我们有, M1V10+M2V20=M1V1+M2V2, 其中M1是车的质量, 上面已经看到约为1600 kg, M2是谭飞的质量,取70 kg 则可.
两个方程联立求解, 可以得到 V2=V20+M1/(M1+M2)*(1+e)*(V10-V20)    

前面根据目击者的报告, 谭卓被车撞后飞起4-5米高, 飞出20米远, 他在被车撞后,做的大致是一个自由抛体运动, 另外查到更加合理的消息: 警方告诉谭卓父亲, 实际被撞出37.3米远, 而不是20米远. 我们取中间值, 假设他被撞飞起4.5米高, 那么他在垂直方向的速度V2y=(2gh)^0.5=(2*9.81*4.5)^0.5=9.396 M/s. 被撞后到达最高点的时间约为 (2*4.5/g)^0.5=0.9578 s

谭卓在空中运动时, 不能忽略空气阻力的影响, 因为根据体育运动学, 100M赛跑如果运动员的成绩是10s左右,他受到的风阻大约为3-5 kGf, 谭卓被汽车撞飞,速度肯定远远超过10M/s,所以空气阻力是肯定不小的. 空气阻力f=0.5*ρ*Cd*S*V^2 ,与速度的二次方成正比, S为物体等效表面积, 相当于扬声器的有效振动面积, 对人体来说, 一般人体皮肤表面积S=71.84*G^0.425*h^0.725来估算, G为体重, h为身高(cm), 则谭卓体表面积S=71.84*70^0.425*170^0.725约为1.81 M^2, 在空中运动, 大概22%的面积受到气流正面作用, 则S=1.81*0.22=0.3981 M^2.  空气密度ρ=1.21 kG/M^3, Cd是空阻系数, 对人体来说, 我们取0.707比较适合. 则有f=0.5*1.21*0.707*0.3981*V^2=0.1703 V^2 [N]

实际上他在空中受到的真实阻力,根据牛顿定律, 有
MVx’=-kV^2*(Vx/V)=-KV*Vx
MVy’=-Mg-kV^2(Vy/V)=-Mg-kV*Vy

得到一组非线性微分方程组,该方程组不存在解析解,无法用常规方法解出, 可以作一些简化处理,

人体在垂直方向获得的速度较小,主要受到重力加速度影响,垂直方向速度迅速变为零, 但垂直方向有阻力,所以实际的到到最高点时间略小于0.9578 s, 可认为是0.9秒, 而后,人体向下运动, 垂直方向受到的空气阻力仍然与运动方向相反, 起到减速作用, 所以从最高点落地时间略长于0.9秒,我们假定它为1.1秒,则得出被撞人体在空中总共飞行时间约为2 s左右.

人体在水平方向,受到阻力较大, 则将方程简化为mVx’=-kVx^2, k前面已经算出, k=0.1703. 解这个微分方程,得到V=0.1703t+M/Vo, Vo是人体水平方向的初速度, 是我们需要求的, 但水平方向飞行路程S=37.3 M, 则37.3=∫(0.1703t+m/Vo)dt, t∈[0,2] , 解该方程,可得到Vo=70(e^(37.3*0.1702/70)-1)/(2*0.1702)=19.5219 M/s  
V2=(Vo^2+Vy^2)^0.5=(19.5219^2+9.396^2)^0.5=21.6654 M/s,
将此值代入前面解出的碰撞方程, 21.6654=-1.3888+1600/1670*(V10+1.3888) 解得, V10=22.6740 M/s
这就我们要的结果, 即车与谭卓碰撞瞬间, 车速约为22.674 M/s, 合81.626 kM/h.

这与警方最后公布的结果的速度下限84.1kM/h很接近. 

误差考虑: 没有考虑到当时的风速, 这是主要的问题, 另外微分方程的算法会存在一点偏差.实际上, 人有一定的反映时间, 胡斌在撞到谭卓之前, 肯定是踩了刹车的, 所以他在撞人之前的车速肯定大于84.1kM/h, 在刹车制动后,速度降了下来,但仍然撞到了谭卓, 所以之前他的车速肯定在100kM/h以上.如果再多一点数据, 我们可以算出胡斌在之前的真实车速是多少.
所以七十码肯定是不对的, 至少是八十多码,但是再多也不可能了. 民众不相信警方,主要是因为官方缺少公信力,实际上官方也有专家,并没有愚弄民众.

结论: 警方没撒谎. 另外从这件事也可以得出一个推论, 浙江大学虽然名气不小,但与南大相比差距还是明显的.  主要是因为, 浙大在杭州只是一支独秀,而南京名校林立, 所以浙大没办法比.

周末讲课--提高班

昨天, 我们看到杭州飙车案的简单计算, 其中很重要的是计算人体在空中的实际轨迹, 一般来说实际轨迹的计算非常困难,
物体在空中受到的空气阻力大小与它的速度模值的二次方成正比,方向与运动速度方向相反,它由一组微分方程组控制,
这组方程组无法用常规数学方法解出.一般在低速运动时这个空气阻力可以化简, 用f=-kV来表示, 当然这里的k与上面二次
型的系数是不同的. 化简单以后的方程比较简单,很容易解出来. 但我们感兴趣的还是如何解决实际的方程, 即描述真实物理
世界发生的事情.

我们来看看能做些什么, 变换一下昨天那组微分方程组的形势,得到它的微分形式:

dVx=-k/m*V*Vx*dt,
dVy=-g-k/m*V*Vy*dt,

约定向右与向上方向分别为x轴和y轴的正方向, 易发现在物体上升过程中与下降过程中垂直方向的阻力加速度是不同的,
上升时阻力为负向(与重力加速度g相同), 下降时阻力为正向(与g相反),
所以下降时的y方向的微分方程还要改一改, 变为dVy=-g+k/m*V*Vy*dt

这里在进一步修正阻力系数, 被撞后的人一般是屈着身子飞的, 迎风面积会减小, 但要考虑人穿着衣服, 所以有效阻力面积要加大,
应该定为40-50%体表面积,取45%,代入前面的阻力方程里算出k大约为0.34839左右, 比0.1703要大一些,这应该更符合实际情况.
m还是70kG, g取9.81 m/s^2.

我们把这个物理过程化分为无限小的时间片断, 假设每一小片段为0.01秒, 那么用dt=0.01代入上述方程组中,逐级累加,就可知道
最后真实运动的情况, 借助EXCEL强大的计算列表功能, 就可得出最后的物体运动轨迹.

经过实验以后,我使用万分之一秒为小时间片断, 计算2秒多的物体运动过程,共计算两万多点的数据, 得出的结果精度是相当高的.
(一般用百分之一秒作时间元,计算两百个数据就可得到满意的精度).

但这里的解题方向需要转换, 假如我们知道1.初速度大小, 2. 起飞角度, 两个已知条件, 用上述方法就可以很轻松地求出物体运动轨迹.
但现在的条件是, 只知道物体大概飞的高度,以及飞行距离, 要反求初速度和起飞夹角,怎么办? 只好一个个去试,最后将方程的解试出来.

虽然空气阻力复杂,但是基于基本常识可想象到: 如果初速度越大, 飞得越高, 飞得越远. 而起飞角度, 在一定范围内, 也可认为与飞行高度
和距离成正向关系. 并且由于实际物理过程的复杂化, 而且受到柯西初值条件的约束(即t=0时, V=Vo, h=0, S=0),
这个方程的解是唯一的. 也就是说, 如果确定了最大飞行高度与飞行距离, 物体的初速度与起飞角度都被唯一固定.

现在目击者说, 谭卓飞起了大约4-5米高, 但到底是4米还是5米,没个准儿.因为当时的情况不能重新去复演一遍.  而最大飞行距离倒是固定的,
因为人落地死了以后就不会动了,警方去测量他的飞行距离比较可靠, 这个距离已知是37.3米.

所以题目变成: Hmax=4-5 M,  S=37.3 M, 求Vo与起飞角度?

经过仔细寻找试解,发现如果Hmax=4 M, Vo=24.88 M/s, 起飞角θ=21.64°, 空中总飞行时间1.7774秒.

将轨迹图形画出来, 见下图, 以红色线表示.
作为对比, 同时画出了在同样出始条件下,理想状态下(即不计空气阻力)的运动轨迹,以蓝色线表示.
蓝色线是一条标准抛物线, 而真实的轨迹红色线是一条无法用任何函数表示的曲线,只是一条光滑曲线.
可以看出, 如不计空气阻力, 最远飞行距离可达43米多, 与真实情况相差达6米,还是比较大的, 所以不能不考虑空气阻力.

 
此主题相关图片如下：truetrace.jpg


用同样的方法, 找出Hmax=5 M, S=37.3 M的根, 得到Vo=23.12 M/s, 起飞角θ=26.36°

将飞4米高的轨迹同5米高的轨迹画在一起, 对比一下,如下图


此主题相关图片如下：trcompar.jpg


我们得到了一个令人惊讶的结果, 即在真实物理环境下, 如果飞行距离固定了, 飞行高度主要受到起飞角度的影响, 而受初速度的影响不大.
所以,不论谭卓飞4米高还是5米高,他的初速度都在24 M/s 左右.

如果是4米, Vx=23.12644 M/s  用水平方向动量守恒, 算出最后车速是 23.12644*1670/1600=24.138 M/s, 合86.898 kM/h
如果是5米, Vx=20.71603 M/s 车速为20.71603*1670/1600=21.622 M/s, 合77.840 kM/h

现在观察肇事车辆三菱EVO9, 发现它的前脸并不是很扁,车窗的角度与发动机盖的角度不同, 而撞人后车窗裂了,说明人是与前车窗有接触,
显然人是被高速行驶的车"铲"到了发动机盖上, 然后又与车窗撞击, 所以起飞角度不可能太大, 21度左右比较靠谱. 所以相撞时86.9km/h的车速
比较可信,应该就是真实撞击瞬间的车速(实际上事故前该车应该以110-120km/h的速度在飙,驾驶员胡斌发现谭卓后急踩刹车,但是已经来不及减速了),


此主题相关图片如下：mitsubishi_lancerevolutionix8.jpg



扬声器的实际阻抗曲线问题与此类似, 很多电声软件都模拟过喇叭的阻抗曲线,但你会发现无论是LEAP的扬声器先进模型也好其它的使用半电感理论也好,它们都与
喇叭的实际阻抗曲线有差异,具体表现在低频时与真实较吻合,而高频部分偏差较大.而真实抛物曲线在前半部分与模拟曲线(标准抛物线)吻合较好,而后半部分偏差较大,
这是随着时间的推移, 误差越来越大. 与喇叭阻抗曲线的随着频率的升高,误差越来越大有同理之处.

这是因为所有的电声软件都意图使用一种算法或函数去表示一个喇叭的阻抗,
简单的使用几个T-S参数, 复杂的如LEAP 5的先进模型使用更多的参数,但仍然不能完美模拟出喇叭的曲线. 这是因为: 任何一个扬声器的阻抗曲线都不能被表示为
一个函数,不存在解析表达式.只能去逐点扫频去测试, 类似于我使用EXCEL逐点去计算函数值.(一般扫频的话扫几百个点就可保证测试精度)

这种解决问题的方法叫做数值法, 得益于计算机的发展,借助于计算机强大的计算能力,使再现真实物理世界成为可能.