在长达两千多年的科学发展史中,重大的科学成就与重要科学方法的应用是分不开的.牛顿发现了万有引力定律进而创立了科学的天文学.由于进行了光的分解创立了科学的光学,由于创立了二项式定律和无限理论而创立了科学的数学,由于认识了力的本性而创立了科学的力学[1].在牛顿力学中,如果已知初始条件,可以预测有序系统的未来的运动状态.牛顿力学对天文学问题的处理是最成功的,科学家可以准确预测宇宙飞船的着陆位置、行星运行的位置、日食与月食的发生时刻等.大学物理中最基本的概念(力、动量、能量等)和原理(牛顿运动定律、动量定理等)在教学中都会详细介绍,学生非常熟悉,经典教材中通常利用微积分方法分析解决一些力学问题.由于功能强大的计算机可以快速完成重复性计算工作,逐步预测方法正在成为现代科学和工程中的一项强有力的技术.近些年来,一些研究者在大学物理教学中引入Matlab计算机数值模拟,培养学生用计算机进行学习和研究问题的能力[2,3],使得学生更加注重研究物理现象本身,从而准确地构建物理模型.计算机数值模拟与实验研究、理论分析是理工科学生的必备素质.本文将介绍一种与传统的物理教学不同的新的方法,利用计算机建模对物体的动力学和运动学进行模拟.计算机建模不仅在物理理论和实验中起着重要作用,而且在几乎所有其他科学和工程领域中也起着重要作用.基于基本物理原理创建简单的计算模型可以更清楚地了解这些原理是如何控制物理系统的实时行为的,这些模型可以将物理原理应用于复杂的系统[4].

本文采用VPython软件对物体运动过程进行实时模拟与可视化研究,从基本物理原理开始,对自然现象建模,利用计算机仿真进行定量研究,用VPython编程语言描述具体的物理行为,通过对实际系统的简化、理想化、近似的估计进行建模过程.其核心是应用动量定理来预测物体或系统的行为;给定相互作用物体的初始条件,计算物体间的相互作用力,根据作用力和很短的作用时间,可以得出很短时间后物体的动量;然后利用新的动量计算物体的位置[5,6].只需编写一段程序进行迭代重复,不需进行复杂的数学运算,设定很小的时间步长即可得到足够高的计算精度.使用合适的程序语言,可以模拟物体在三维空间中的运动,实时观测物体的速度变化及其运动轨迹,动力学和运动学有机的结合,物体的运动源于相互作用力,同时深刻感受到牛顿力学随时间演变的特性.

本文以万有引力作用下天体(卫星、行星等)的运动为例,详细介绍如何应用VPython进行天体轨道运动模拟与可视化研究.

1、万有引力定律及程序语言VPython简介

1687年牛顿发现了万有引力定律,任何物体之间都有相互吸引力,万有引力的大小与相互作用的物体质量成正比,而与它们质心之间的距离的平方成反比,方向沿两物体之间的连心线,且与两物体的化学组成和其间介质种类无关.如果用m1、m2表示两个物体的质量,r表示两物体质心间的距离,则物体间相互吸引力大小为F=Gm1m2r2,G为万有引力常量.该定律对人类理解自然具有重要作用.在本文中只考虑天体间的万有引力作用,通过软件VPython来模拟天体的运动轨迹,有助于直观地理解现实世界中如何应用动量定理预测物体的运动行为.

Vpython是广泛使用的Python编程语言的扩展.VPython支持3D矢量运算,并允许创建实时动态3D动画作为物理计算的直观表现[4].用VPython模拟的两个物体间相互引力如图1所示.

图1两个物体相互引力模拟图

天体在万有引力的作用下,其运动轨迹与初速度有关,可能是匀速圆周运动、椭圆、抛物线、双曲线或其他运动形式.天体在万有引力作用下是如何运动的?初始条件和万有引力如何决定天体的运动?采用微元法的思想,当dt(时间间隔)足够小的时候,可以近似认为天体在dt(时间间隔)内做匀速直线运动.已知万有引力常量为G=6.678259×10-11N·m2/kg2.为了描述天体的运动,通常定义下列物理量:r,rˆ,F,v,s,其具体的物理意义如下:r=r1-r2表示一个天体相对于另一个天体的相对位置矢量;rˆ=r|r|表示两天体连线方向上的单位矢量;F=−GMmr2rˆ表示天体相互间的万有引力;P=P0+F×dt表示计算天体更新后的动量;v=Pm表示天体的瞬时速度;s=s0+v×dt表示天体更新后的新位置。

上述各个物理量(位置、万有引力、速度、动量)在计算时均采用矢量的形式,这样便于在VPython中模拟星体运动.

2、用VPython模拟不同情况下的天体运动

2.1模拟火星绕太阳的运动

利用天文学中真实的数据进行分析与模拟,可以与真正的天体运动联系起来,例如在探究火星绕太阳运动轨迹时,忽略其他星体,只单独对该双星系统进行研究,太阳的半径为Rsun=6.955×108m,其质量为Msun=1.989×1030kg;设程序运行初始时刻,太阳中心距离火星中心的距离为L=2.28×1011m,火星的半径大小为RMars=3.3895×106m,其质量为MMars=6.42×1023kg[4](部分数据来自于百度搜索).利用VPython创建一个引力系统模型来模拟火星绕太阳的运动.

首先对太阳、火星形状和位置进行定义,火星和太阳都是球状物,根据质量和半径以及相互间距,利用VPython显示出太阳和火星的大小、颜色以及相互位置.然后将球体的轨道可视化,在VPython模拟时可以直观地模拟并观察它的运动过程及轨迹,程序中使用的物理公式及对应的程序循环语句如下:

whileTrue:

rate(10000)

r=rmars-rsun

r=mars.pos-sun.pos

rˆ=r|r|rhat=r/mag(r)

Fon_mars=−GMsunmmarsr2rˆ

Fon_mars=(G*M*m/mag(r)**2)*(-rhat)

dPmars=Fon_mars×dt

Pmars=Pmars+Fon_m*dt

dSmars=vmars×dt,P=mmars×vmars

mars.pos=mars.pos+(Pmars/m)*dt

Fon_sun=-Fon_mars

Fon_sun=-Fon_mars

dPsun=Fon_sun×dt

Psun=Psun+Fon_sun*dt

dSsun=vsun×dt,Psun=Msun×vsun

sun.pos=sun.pos+(Psun/M)*dt

其理论原理利用了微积分的思想,将运动时间分成许多很短的时间间隔,在很短的时间间隔里(足够小),物体运动的速度可以近似看成是不变的,因此在dt时间内,物体近似做匀速运动,利用运动学公式s=s0+vdt,(v,s均为矢量),可以计算出物体各个时刻的位置.利用动量定理,即物体的动量改变量等于所受冲量的大小Pf-Pi=Fon_mdt,v=P/m可以计算出更新后的速度v,物体更新后的位置s,故经过程序循环可获得其各个时刻物体的位置以及运动轨迹.利用天文学真实数据以及Vpython模拟,火星绕太阳运动轨迹如图2中黄线所示,红色圆点表示太阳(为显示更清楚,将太阳的半径增大了4倍),由图2可知火星绕太阳运动的轨迹是一个椭圆,并且太阳的位置近似处于椭圆的一个焦点位置处,其模拟结果与开普勒第一定律及教科书上理论推导结果相一致.

图2火星绕太阳运动的轨迹

2.2卫星、月球绕地球的运动(三体问题)

在上述火星绕太阳运动模拟过程中,只需求出两个星球之间的万有引力,然后通过VPython来创建计算模型模拟运动轨迹,在研究三体运动过程中,则需考虑各个物体之间的万有引力作用,因此单个星体所受的力的个数会变成两个;用解析的方法研究三体运动问题求解过程十分复杂.本文利用天文学真实数据,模拟卫星、月球绕地球的运动,现已知地球的半径REarth=6.4×106m,质量为MEarth=6×1024kg,月球距地球的起始距离L=3.84×108m,月球半径为RMoon=1.7×106m,月球质量为MMoon=7.3×1022kg,卫星在离地球中心十倍于地球半径的地方运动,卫星的半径为RSatellite=6.7m,质量为MSatellite=1.5×104kg,分别给定vSatellite=<0,2×103,0>(m/s),vEarth=<0,0,0>(m/s)(忽略地球的自转的影响)[4].vmoon=<0,1×103,0>(m/s),根据3个物体间的相互引力作用,在设置了初始条件后,通过VPython来模拟.程序中使用的物理公式及对应的程序部分语句如下:

whileTrue:

rate(1000)

r1=rcraft-rearth

r1=craft.pos-earth.pos

rˆ1=r1|r1|

rhat1=r1/mag(r1)

r2=rcraft-rmoon

r2=craft.pos-moon.pos

rˆ2=r2|r2|

rhat2=r2/mag(r2)

r3=rmoon-rearth

r3=moon.pos-earth.pos

rˆ3=r3|r3|

rhat3=r3/mag(r3)

Fon_craft1=−GmcraftMearthr21rˆ1F_oncraft1=G*m*M*(−rhat1)/mag(r1)**2

Fon_craft2=−GM1moonmcraftr22rˆ2F_oncraft2=G*M1*m*(−rhat2)/mag(r2)**2

Fon_moon1=−GM1moonMearthr23r3F_onmoon1=G*M1*M*(−rhat3)/mag(r3)**2

地球、月球、卫星的运动轨迹如图3所示.

图3卫星、月亮绕地球运动的运动轨迹

在上述模拟中,作用在每个星体上的引力都是两个,其所受合力将会改变星体的动量,在很小的一段时间内,星体的运动可以近似看成匀速直线运动,通过动量求出星体各个时刻的速度,然后根据运动学公式计算出星体在各个时刻的位置,通过转化成VPython程序来模拟,从图3中可以看出,月球的轨迹近似为圆,卫星的轨迹近似为以地球为焦点的椭圆.

2.3五颗卫星绕地球的运动(多体问题)

在谈论了三体问题后,如果所研究的天体数目继续增加,那么通过VPython模拟出的运动过程,又会有什么不同呢?现假设有五颗卫星绕地球运动,由于卫星的质量远远小于地球的质量,在编辑程序时,忽略各卫星之间的引力影响,并通过VPython建立计算模型来模拟卫星运动,现在5颗卫星相对于地球的初始位置(单位为m)分别为:

r1=<−13×107,6.5×107,0>,r2=<−6.5×107,6.5×107,0>,r3=<6.5×107,6.5×107,0>,r4=<0,6.5×107,0>,r5=<13×107,6.5×107,0>

人造卫星的半径大小为RX=6.7m,质量为MX=1.5×104kg,给定每个卫星的初始动量为(单位kg·m/s)[6]:

<1.65×107,0,0>,<3.1×107,0,0>,<3.75×107,0,0>,<3.1×107,0,0><2.5×107,0,0>

程序中使用的物理公式及相应的程序循环部分语句如下:

whileTrue:

rate(500)

r1=rs1-rplanet

r1=s1.pos-planet.pos

rˆ1=r1|r1|

rhat1=r1/mag(r1)

Fon_s1=−GMplanetms1r21rˆ1F_ons1=(G*M*m/mag

(r1)**2)*(-rhat1)

dPs1=Fon_s1×dt

Ps1=Ps1+F_on1*dt

dSs1=vs1×dt,Ps1=ms1×vs1

s1.pos=s1.pos+(Ps1/m)*dt

5颗卫星围绕地球运动轨迹如图4所示.由于地球半径远远大于人造卫星的半径,故人造卫星在图中不可视(与轨迹线重合).

图45颗卫星绕地球运动轨迹

其中有的卫星的轨迹为椭圆,有的卫星的轨迹为圆,这是由于各个卫星的初始位置以及初动量的不同而造成的,故初始位置和初速度会影响星体运动轨迹的形状,下面通过将改变卫星的初速度来观察其运动轨迹的变化,当把第一颗卫星的初速度由v=<1100,0,0>(m/s)改为v=<1667,0,0>(m/s)时,第一颗卫星的运动轨迹如图5所示.由图5可知此卫星的轨迹由斜椭圆变成了正椭圆.

图55颗卫星绕地球运动轨迹

2.4卫星绕固定不动的双星运动

在多体问题中,卫星的初始速度的改变会导致其运动轨迹的改变,本文进一步研究了随着速度的变化,其运动轨迹变化是否有所规律.为了易于观察,假定地球的位置为posEarth=<0,0,0>(单位为m),月球的位置为posMoon=<-3.84×108,0,0>(单位为m),都固定不动,人造卫星的初始位置为posSatellite=<6.4×107,0,0>(单位为m),人造卫星的半径为RSatellite=6.7m,质量为MSatellite=1.5×104kg[4](部分数据来自于百度搜索).设定人造卫星的任意初速度v=<0,3.4×103,0>(m/s),用VPython模拟其运动轨迹.程序中使用的物理公式及对应的程序部分循环语句如下:

WhileTrue:

rate(10000)

r1=rcraft1-rearth

r1=craft1.pos-earth.pos

rˆ1=r1|r1|

rhat1=r1/mag(r1)

r2=rcraft2-rmoon

r2=craft2.pos-moon.pos

rˆ2=r2|r2|

rhat2=r2/mag(r2)

Fon_craft1=−Gmcraft1Mearthr21rˆ1F_oncraft1=G*m*M*(−rhat1)/mag(r1)**2

Fon_craft2=−GM1moonmcraftr22rˆ2F_oncraft2=G*M1*m*(−rhat2)/mag(r2)**2

本文模拟了具有不同初速度的卫星的运动情况:

v0<v<v1,v1<v<v2,v>v2

v0=〈0,1.9×103,0〉m/s

v1=<0,3.26×103,0>m/s,v2=<0,3.26×103,0>m/s

经过模拟实验,当v0<v<v1时,卫星会绕地球做近似圆周运动(如图6(a)所示);当v1<v<v2时,会绕固定双星系统做闭合曲线运动,轨迹变为椭圆,如图6(b)所示.继续增大初速度,v>v2时,卫星的运动状态由束缚态变为非束缚态,轨迹不再是闭合曲线,卫星开始脱离地球,运动轨迹为抛物线(如图6(c))或者双曲线(如图6(d)).

图6卫星绕地球、月球系统运动图

有趣的是,在系统受力相同的情况下,即所设情景和上例相同,地球和月球固定不动,改变卫星不同的初速度数值,用VPython模拟得到一些特殊的运动轨迹,由模拟实验可以得到当设定初速度为v=<0,3.2634×103,0>(m/s),模拟出的运动轨迹呈八字形(如图7所示).

图7卫星绕地球、月球系统运动,轨迹为8字形

对于三体或多体相互作用及运动问题,用解析的方法分析物体的运动是非常复杂的,因为涉及复杂的积分运算,但是本文通过迭代程序方法,可以模拟多体运动问题,根据设定相互作用力及初始条件即可预测物体未来的行为,更重要的是,可以观测其运动过程,非常直观形象且与物理原理紧密联系,应用基本的物理原理即可解决复杂的多体运动问题.

3、结语

本文采用VPython软件来模拟天体的运动,利用天文学真实的数据以及相应的理论知识,可以模拟得出天体的运动变化轨迹,将物体三维运动可视化,描述了天体运动实时动态变化过程,更加清楚、直观地反映了星体运动的轨迹,从而准确地构建物理模型.有趣的是,可以实时预测复杂体系的运动,如多个万有引力相互作用物体的运动轨迹(行星,宇宙飞船,双恒星),使得对复杂体系的动态研究变得形象直观.这些建模研究是经典的大学物理教学内容的有益补充,通过建模培养学生使用计算机进行学习和研究问题的能力,使得学生更加注重研究物理现象本身.

参考文献：

[1]查有梁.牛顿力学的方法——纪念牛顿《自然哲学之数学原理》发表300周年(1687～1987)[J].自然辩证法研究,1987,5:1.

[2]何钰,吴平,张晓,等.Matlab在物理学中的应用[J].大学物理,2013,12:39-42.

[3]彭芳麟,梁颖,忻蓓.计算软件在计算物理课程中的地位和作用[J].大学物理,2013,8:6-11.

吴永昊,刘玉颖,宋敏.基于VPython的天体轨道运动模拟与可视化[J].大学物理,2020,39(02):69-73+77.

基金:2019年高等学校教学研究立项项目(DJZW201907hb);2017年高等学校教学研究立项项目(DJZW201703hb)资助.