巩伊鸿
(北京师范大学 物理学系,北京 100875)
蜡烛动力涡旋机常作为一种儿童手工玩具,其中蕴含丰富的物理知识. 追溯古代,有原理与蜡烛动力涡旋机类似的“走马灯”. 《西京杂记》中记载:“有青玉五枝灯,高七尺五寸. 作蟠螭,以口衔灯,灯燃,鳞甲皆动,焕炳若列星而盈室焉.”汉代时叫“蟠螭灯”,后从皇室流入民间,作为节日活动的灯具[1]. 走马灯底部放置蜡烛,热气流推动上方轴轮转动,四面与轴连接的剪纸随之转动,烛光将图案投影到灯罩上,连转不休,颇为精巧[1,2].
关于蜡烛何以提供动力,现有文献基本简言之[1,2]:热空气上升,四周冷空气补充,循环往复. 鲜少有探究扭矩的影响因素的,更未见定量计算蜡烛动力旋转物的扭矩的. 至于蜡烛产生的流场性质,有采用纹影进行可视化,观察到热气流上升与冷气流下降[3,4],暂无将纸螺旋旋转过程用纹影可视化或采用软件模拟的.
本文一方面研究蜡烛上方流场的性质,并采取合理假设以便理论推导扭矩的影响因素. 通过实验获取数据,带入理论公式计算出扭矩. 最后,理论与实验相结合,对3个主要因素如何影响扭矩得出半定量的结论.
在重力场中达到平衡的流体,其密度只与高度有关,放置燃烧的蜡烛后,力学平衡被破坏,产生自然对流. 本部分将就自然对流形成的流场性质进行说明.
首先假设流体不可压缩,那么压强沿流体流动方向变化要足够小,从而可忽略压强变化引起的密度变化. 文献给出近似不可压缩流体满足[5]
(1)
设研究的气柱高度h不太高;g为当地重力加速度;c为当地声速;α为定压体膨胀系数;Θ为特征温差.式(1)两边同乘初始密度ρ0后,左边表示这段高度气柱的压强导致的密度变化,右边则表示气柱温度变化导致的密度变化.要提醒一点,不可压缩的假设并不意味着要求温度变化导致的密度变化可忽略.取以下数据进行数量级估算,如表1.蜡烛上方气柱温度最高处为火焰外焰,取800 K,对应最小的α.室温取300 K,则特征温差相应取室温与外焰的差值500 K.计算得左边数值8.65×10-5,右边数值0.625,满足左边远小于右边.因此在本文实验环境下,对蜡烛上方气流不可压缩的假设是合理的.
表1 式(1)各量近似值
对表1中定压体膨胀系数α做一点补充说明.本文中
(2)
与一般热学教材表述[6]
(3)
(4)
然后假设蜡烛上方流场,即流体失稳导致的对流是定常的.以下证明这一假设.重力场中的纳维-斯托克斯方程、自由对流过程忽略了黏性项的传热方程以及连续性方程构成定常流自由对流的描述方程组
(5)
ν是动力黏度,χ是温导率.T′是所在高度的温度与某一恒定温度的差,不妨认为其参照点是点燃蜡烛前的室温.T′可表示为两项和,一不变差值与一小扰动量τ.同理有p′是压强与某恒定压强的差,参照点为点燃蜡烛前静力学平衡状态下的压强,对应的扰动量为ρε.ρ为此刻的流体密度.v是流速.设扰动对时间t的依赖关系为e-iωt.由式(5)可导出小扰动量满足的齐次线性微分方程组[6]
(6)
式(6)是无量纲化后的结果.设z竖直向上,n代表z方向单位矢量.Pr、Ra分别为普朗特数和瑞利数.近似认为蜡烛与空气接触面温度恒定,则接触面上边界条件为:v=0,τ=0.式(6)与边界条件组成定解问题,描述自然对流.将式(6)的前两式分别乘以v、τ的共轭v*、τ*,沿所研究的这段气柱积分
(7)
V表示流体体积.式(7)两式取共轭分别与对应原式相减得
(8)
(9)
将Ra乘以式(9)与(8)相加得
(10)
为流场定常这一概念进行必要的补充说明.定常指速度等参量只是坐标的函数,不随时间变化.与层流湍流无关.为使流体保持定常,在与流体接触的那些固体内必须有热源,使固体的温度保持不变[5].后续的实验中,蜡烛持续燃烧,一次实验持续短短几分钟,近似认为蜡烛外焰的温度不变,满足定常假设.
表2 Re各量取值
“圆管直径”取蜡烛直径;气流速度过小无法用风速仪测出,故流速一定小于风速仪可测最小值0.4 m/s; 气体运动黏度随温度升高而增大,故取室温25℃时的最小值[8].代入后计算出实验条件下最大雷诺数,Re≈1 282<2 000,故可认作稳定层流.
本部分基于文献,取实际数据验证,得出蜡烛产生的流场近似为不可压缩的定常层流的结论.
自然对流包括热空气上升和冷空气下沉,两者在蜡烛产生的流场中,哪个对纸螺旋的旋转提供扭矩?流场影响范围或气流截面有多大?风速太小无法明显感受,有无方式能直观地看到上述问题的答案?或许模拟和特殊的实验可以回答.
首先采用COMSOL Multiphysics软件对蜡烛产生的流场进行仿真模拟.选择“共轭传热-层流”物理场和“稳态”研究,设定相关参数,求解结果如图1所示.
图1 流场的仿真模拟
又用纹影法观察流场.温度不均使空气折射率分布不均,光线通过这部分流场后发生弯折,又经刀片遮挡,到达光屏时强度也不再均匀,屏幕上便出现了明暗图样.搭建反射平行光纹影装置,如图2所示.
图2 纹影装置示意图
“Z”形排列,使入射光与反射光分开,从而有效消除相差.刀片置于凹面镜焦点上,遮挡一半光线,成亮度减半的倒立实像。纹影法反映了折射率的梯度.通过纹影仪,观察了无其他外界干扰时燃烧的蜡烛上方的流场和将纸螺旋放在流场中旋转时的像,如图3所示.
图3 纹影法观察蜡烛流场
从图1可以看出,上升的热气流为主导,故在下一部分的理论推导中简化模型只考虑上升的热气流;从图3可以看出,气流截面较小且主要作用于纸螺旋的最下层,故实验部分对“高度”的定义为蜡烛距纸螺旋的最下层.
理论采用三维等速螺旋线,投影到平面即为阿基米德螺旋线,如图4所示.
图4 等速螺旋线
为方便推导,假设纸螺旋为大小、形状不变的刚体.应用于实验时,会进行修正,详见3.2.
设转轴为z轴,正方向竖直向上.无蜡烛时纸螺旋静止,由刚体平衡的充要条件[15]:
(11)
(12)
可知重力与支持力平衡,且对转轴的合力矩为0.因此纸螺旋重心在转轴上,不需要考虑重力矩.动力矩(即扭矩),是由气流施加的在水平方向的分力提供的.
取一小段宽度d,倾角θ的纸元进行分析.假设:1) 空气分子与纸面发生完全弹性碰撞;2) 只用一支蜡烛,提供动力的总力矩即为Mz.当密度为ρ,速度为v的气流冲撞到长l的一段纸螺旋上,Δt时间内气体分子在垂直纸元方向的动量改变量等于对纸元施加的冲量:
ρvΔtdl·2vcosθ=FΔt
(13)
整理得气流施加在垂直纸面的力为
F=2ρdlv2cosθ
(14)
分解到水平方向上有
Fxy=Fsinθ=2ρdlv2cosθsinθ
(15)
扭矩为(其中i的正方向与旋转的切向速度相同)
Mz=r×Fxy=iρdlrv2sin 2θ=
(16)
图5 气流施加冲量示意图
首先定义流管:在流体内作一微小闭合曲线,通过其上各点的流线所围成的细管[9].假设:1) 定常、不可压缩流场;2) 将纸螺旋面视为由实物构成管壁的流管的一小片;3) 流管沿螺旋面转弯,流管横截面积S不变.
图6 流管模型示意图
设螺旋面与水平方向夹角为α,由定常流场的连续性原理,一段流管流入的质量流量等于流出的质量
dQm1=dQm2
(17)
ρ1v1dS1=ρ2v2dS2
(18)
由假设“不可压缩流”知密度不变
ρ1=ρ2=ρ
(19)
结合假设3得v1=v2.得到各量几何关系如图7,且△ABC是等腰三角形.AD代表纸螺旋的表面,α是纸元面与水平方向的夹角,AB是气流冲撞纸元面后改变的方向,与原方向关于纸元面的法向AE对称.由几何关系得
(20)
∠EAB=π-α,∠CBA=α,则AE//BC,说明气流碰到纸螺旋后方向改变,且促使转向的力在它们的角平分线上,由纸螺旋提供.或者说气流冲撞纸螺旋,沿纸螺旋的法向产生冲力.
图7 方法二各量几何关系
在时间dt内流体受管壁作用力满足
Qm2dtv2-Qm1dtv1=Fdt
(21)
流体反作用于螺旋面的力为
F′=-F=Qm1v1-Qm2v2
(22)
由图7几何关系并结合式(17)—(19)得(22)中,力的大小为
(23)
(24)
Mz=r×F′xy=iρSrv2sin(π-2α)=
(25)
显然得到了与方法一相似的结论.两个方法角度不同相互映证,说明了纸螺旋扭矩与影响因素的关系,为后续实验提供研究方向.
三维等速螺旋线的参数方程为
每增加一层(二维投影角度转过2π),纸螺旋到转轴的距离Δr增加πa,三维高度Δz增加πc.纸螺旋的宽度d=πa.得比例关系:
(26)
图8 一小段三维螺旋线
设纸面密度为η,质量m、长度l、极角θ(ωt)等微分量的关系为
dm=ηaπ dl
(27)
(28)
(29)
dθ=ωdt
(30)
dz=cdt
(31)
将式(29)—(31)带入(27)、(28),得
(32)
则三维等速螺旋的转动惯量微分形式为
(33)
积分得
(34)
为与实验保持一致πa(d)取1.5 cm;ω取1 rad/s;Δz取范围[0,15]cm;t取范围[0,15π].用MATLAB绘制转动惯量随参数的变化图,如图9所示.
图9 转动惯量与层数、螺距的关系
可以看出转动惯量随层数、螺距的增加而增加.实验中采取的方法是:测量每一层的高度等,逐层代入积分式,再将各层的数值相加得到整个纸螺旋的转动惯量.
将纸螺旋置于一坐标纸前(实验装置详见4.1),对不同纸螺旋,总是从同一位置同一角度拍摄,并保持其末端大致在同一方向.拍摄后通过坐标纸测量每层的层距,最后按前述方法带入公式计算转动惯量.
图10 螺距测量
得到6种规格纸螺旋的转动惯量如表3.表中规格表示“克重-层数”,例如“150-6π”表示“150 g·m-2的牛皮纸制三层纸螺旋”.
表3 不同规格纸螺旋的转动惯量
通过理论分析和预实验,可供实验探究的变量有:纸螺旋的悬挂高度、层数和密度.理论中气流速度体现在高度上;受力面积和到转轴距离体现在层数上;而纸张密度则影响转动惯量.实验采用控制变量法探究三者对扭矩的影响.
不易形变的230/150克重牛皮纸;同规格蜡烛若干;MATLAB绘制t=6π;8π;10π的螺旋线(3层;4层;5层),均按比例调整宽度为1.5 cm,打印在纸张上裁剪出纸螺旋;顶部尖而摩擦小的铁签并搭建可升降支架;后面置一高筒,粘贴坐标纸,如图11.
图11 实验装置示意图
将手机固定在纸螺旋正上方,俯拍由静止开始转动、持续转动到撤走蜡烛后逐渐静止的全过程视频. 实景如图12所示.
图12 实验实景
将完整视频导入Tracker软件,定标纸螺旋宽度为1.5 cm,新建质点标为纸螺旋末端的红点,初始位置定为角度0°. 手动逐帧追踪质点全过程的角度变化,将θ与t的数据导入Excel.
图13 Tracker逐帧追踪
将数据导入MATLAB绘制θ-t散点图,观察到撤走蜡烛前的中间阶段疑似匀速转动. 示例如图14所示.
图14 角度随时间变化散点图
验证匀速:读出每条斜线上水平位置大致相同的点的时间坐标,相邻者相减求出每转一圈所用时间(周期). 用标准差除以平均周期求相对误差,最大为14%,最小4%. 故认为误差允许范围内,撤走蜡烛前很长一段时间内纸螺旋匀速转动. 因此这一阶段动力矩近似等于阻力矩.
接下来设法求减速阶段刚开始时的阻力矩.低速运动,空气阻力与速度一次方成正比,设另有一常阻力f,即
(35)
特征根方程为
(36)
θ(t)=C1er1t+C2er2t
(37)
用实际数据拟合后发现常阻力通常不占主导,当常阻力可忽略,用θ(t)=ae-bt+c就可以得到很好的结果.对拟合方程求两阶导,取t=0时的值即为阻力导致的角加速度k,则Ik近似为动力矩(扭矩).
1) 摘出撤走蜡烛后减速阶段的数据,θ、t各一列.将t列均减去第一个值,即起始时间设为0.
2) Tracker中,角度是由0到-180°,再跳到+180°,然后递减回到0.因此对θ列,每从负到正时就减360°,转化为单调有序数列.再统一减去第一个值使初始角度为0.最后化为弧度制.
3) 将上千个数据点导入MATLAB的curve fitting tool自定义函数拟合并计算扭矩(方法见4.2). 拟合图像如图15所示.
图15 θ(t)=22.64e-0.14t-22.31
图15拟合的均方根误差为0.167 7,实验中最大的均方根误差为0.385 9,表明拟合效果良好.
首先研究密度的影响. 由表3可知,层数较少的6π规格纸螺旋两者螺距基本一致,为控制变量,故用两个三层纸螺旋探究密度的影响. 选取了两个高度进行比较,第一列代表坐标纸上不同高度的标号,数字越大,高度越高. 计算得数据如表4. 150克重的纸螺旋,转动惯量与角加速度都更大,因此密度小的牛皮纸制的纸螺旋扭矩更大.
表4 密度对扭矩的影响
再看层数的影响. 由于150克重多层的纸螺旋随温度升高形变逐渐严重,使转动惯量计算值与实际有较大误差,故用230克重的进行本组实验. 保证纸螺旋下端与蜡烛的距离相同,得到数据如表5. 发现扭矩随层数的增加而增加,猜想二者间存在函数关系,用幂率拟合如图16所示,均方根误差为0.385 9. 图16表明扭矩与层数4次方成正比.
表5 层数对扭矩的影响
图16 函数关系y=0.034x3.907
最后看高度(蜡烛距离纸螺旋最下端)的影响.本组实验使用同一纸螺旋,具有相同的转动惯量,因此稳定时的转速ω一定与扭矩正相关,不妨比较ω,数据如表6.在很长一段高度内,转速变化并不大,仅随高度升高而略微变大.可推知扭矩随高度变化并不大.
表6 高度对转速的影响
蜡烛上方出现自然对流形成不可压缩的定常流场,处于流场中的纸螺旋的扭矩与纸张密度、层数、悬挂高度等因素有关. 密度越小、层数越多,扭矩越大;而扭矩随高度变化不大,并有随高度增加而略增加的趋势.
创新点主要有:用纹影观察纸螺旋在流场中的转动. 纸螺旋太过轻小,难用平衡法测扭矩,故采用数学方法计算扭矩的近似值.
实验数据点较少是因为,逐帧追点处理整个视频耗时过长. 最好获取更多数据去检验结论.