北大力学老教授的“蓝天一梦”-资讯-知识分子

北大力学老教授的“蓝天一梦”

2022/11/02
导读

撰文|陈耀松(北京大学力学系教授、中国计算流体力学领域泰斗)

编辑|马修

2003年突然来了SARS(非典型肺炎),通知要求“断路”“封门”,各自留在家中“自我安排”。幸好已有计算机和互联网,依旧可以跟大家交流联系,不致“孤独难耐”。笔者就用这段时间写了点回忆,其中一篇关于航空技术,题为《重归蓝天》。时隔17年,来了“新冠”,又被“关”在家中,想起这篇旧文可以更新了。网络比以前快了,只是“回忆”的对象是历史,可以补充,不得更改。这个“补充”将包括事后的“成果”,由此可以检讨当初设想的是与非,算是笔者的“蓝天一梦”。
01
北大自制的计算机“6912”在国内首次实现CFD全机气动力计算
“文革”后期要求抓革命促生产,北京大学最“耀眼”的成绩,是于1973年与738厂合作制成“百万次电子计算机”。次年,由多家单位合作,生产了八台集成电路中型计算机,面向各行推广。当时我们在汉中分校,分到一台。“从无到有”,真是科技设备“大跃进”了。“文革”后期,初期的“革命热情”已经消退,学校除完成“工农兵”学生的教学以外,还鼓励各家自行寻找工作的“课题”。
自1958年国家仿制苏联BASM计算机(每秒万次,编号104机)以来已过近20年,我们只有伸着脖子”观望“的份儿。1974年,汉中分校有了北大自产的计算机。在山沟中最关心这架大型计算机的,是几个力学系教员。尤其对我而言,就像老天专为我准备的礼物。
解放前,我父亲曾任职“震旦铁工厂”,那里是国内唯一仿制IBM机械计算机的厂商。我考入清华后,发现航空系有一架损坏闲置的IBM原装机,容许我摆弄。于是,我有空就拿个改锥去捣鼓,希望碰巧能修好。最终我也没能修好,却激励了我去研究电子计算机。从二进制到冯·诺依曼计算机原理,再到以后的BD-200操作系统;从半导体、门电路、CPU原理,最后到“6912机电路图集”都通读完毕。
当时的电子计算机都是散装的,管理制度比较严格,各项都有专人负责,外人不得触动。而我在山沟中,因专职照相出了名,加之计算机的专业人员原本都是数力系的熟人,所以,除了主机以外的工作我都可以去帮忙(参加)。
科技计算都用“高级语言”,需要在计算过程中修改参数,这得知道它的内存地址。我专为此设计了一个搜索步序,照此操作,“门牌号码”随即显现。磁带读码有错,计算机会自动退至起点重新再试。这样非常花时间,我就改操作系统 ——退回五步就重读。操作系统要改坏了,无非重新装一遍。至于硬件,万一弄坏就得停机待修,责任非同小可,专家(都是朋友)修机时我陪着聊天,不敢去碰。要是辅机有问题,那就好办些。我们配套引进的六台纸带穿孔机,有一台开箱就不工作,机房答应我来修。我把说明书拿回家,每天一早起来就学习。吃完早饭拿个万用表就去机房检查、分析。最后判定是一只三极管坏了。我把分析结果告诉机房,由机房专家亲自动手来修理。其实,我的焊接技术并不比他差,只是万一碰到意外,这责任我负不起。
1965年,开始在“三线”布置重点院校和科研机构。北大汉中分校代号是“653”,离我们不远便是“29基地”。大家都从北京内迁,最初各自“闹革命”,互不来往。后期,“闲得发慌”,我们以促生产名义与29基地彼此交流。
我们两家是“同行”,自然先“交流”。29基地一位名叫安继光的老总,是航空空气动力学的著名教授,我们就请他来分校作“飞机空气动力学计算”的报告。当初,计算机资源有限,全机计算只限“位势流”,计算就是以飞机外形为“边界”,以数字形式求解“调和方程”。正好汉中邻县城固迁来大型飞机制造厂,厂校结合拿来“运八-全机气动力计算”的课题。
作者供图
全机,也就是机翼、机身、机尾、发动机舱都得考虑,我们的计算机要应付这样的大型课题有点难。主要的矛盾是工具小、对象大。组建与拆解大件都用到“分解”。但组建成大件必须解决小件之间的结合问题。我们面对的计算课题,就是设法算出分部之间相互影响的程序,以这些影响作为“粘合剂”将分部计算的结果结合成整体。
北大自制的计算机“6912”以磁芯作内存,外存仅有磁带。国产磁带质量不过关,一寸宽的磁带分12磁道,即便5道用作校验位,也未读成过一次。录入的数据回放总是出错,然而数字计算容不得一个错,以致我们没有外存可用,整个计算必须“一气呵成”。
真实飞机的外形十分复杂,将各部计算的程序凑成一体,计算中途的数据当然亦需保存。由此可知,计算机必须有非常大的内存。“6912”内存总共只有384k字节(byte)(注【1】)。计算机工作必须有操作系统,它占去200k。留给用户计算的内存不到200k。当年,一般计算机一个浮点数占用4字节(32bite)。我们穷,只得放低精度,用3字节。计算采用“面源法”,先用高斯消去法求解。N-面元对应的系数方程的系数就有N2个。仔细算来,“6912”只能解算 200个面元的机翼。要算全机——机翼、机身,再加两个发动机舱一共 800多个单元,根本不可能。
为解决这一困难,我们尝试:
1.改用雅可比求解线性方程组代替高斯消去法,以同样的内存,可解方程元素扩大一倍。
2、专门设计了“干扰提取算法”,以此计算分部之间的影响(注【2】)
对此作简单介绍:

若按图一的面元分布计算气动力,块数超过计算机内存容量;
现按图二重新分布面元。要求机身机翼接近、相互干扰大的区域面元块分布密度加大,远处密度较小,使总块数在计算机内存容许范围之内,进行气动计算。
保留机身所得面源强度分布(即保留机身对机翼的影响源)而对机翼细分(如图三)重新计算机翼气动力,得所需结果;
按同样的思路,可得所需的机身气动力。
就按这一方法,我们在1976年算得带四个发动机舱的“运八”全机气动力。这是国内第一份以真实飞机为对象的全机气动力计算。由此计算的机翼气动力分布与实测核对得到肯定。

最后计算是在西安总参同一型号的计算机中进行的,因为用到“雅可比求解方法”,耗时特别长,总共用了6小时。
02
对FTC2000新型飞机进行全机气动特性数值模拟
这是我们第一次为国内独立设计飞机作气动力计算服务。因为直接面向军工产品,我们组织三个计算组分头进行,做到结果一致后,再邀请603所气动计算负责人张锡金教授亲自审核。

由此有一“插曲“:

如今“山鹰号”大量生产,为培训航母起飞着陆以及出口作巨大贡献。
03
参加航空航天系统院所“CT-1”标模大攻角计算评比
2005年初,得知航空航天系统院所“气动计算基金”结题中有一项标模大攻角计算评比活动。
因为这是一项纯技术交流,我组申请参加,并得到容许,寄来计算目标。

这是一个超大攻角飞机计算的难题,攻角要求算到105度。
我们先组织两组,同时用 NS方程计算。攻角都算到45度就上不去了(按50度设置,不收敛)。于是改用LBM计算,收敛了。
powerFLOW是一份基于LBM的流体计算软件。它在计算有分离的流动和计算噪声方面非常成功。高档汽车的气动设计都用它,而没有在飞机气动计算上用过,我们这次是首次。
2005年6月底,我组派王一伟带了NS方程算到攻角45度和LBM算到攻角105度的两根曲线去参加评比。各家将自己计算的曲线交给秘书,开会那天,由秘书在开会时统一发布:

计算曲线发布后,秘书再将CT-1标模实验测得的曲线(黑色)加上。对比显示,我们的计算比较好。
其实,我们在出发前已在组内讨论,为何攻角大于 45度时 NS方程计算不收敛,而 LBM 能’收敛’?分析原因是所谓“NS方程计算”是用定常RANS方程迭代到“收敛”为解。而LBM的基础是“分子运动论”,天然是非定常的。它的所谓定常结果,无非是数个迭代步的平均值。当攻角较小时,尾涡“粒度”不大,“湍动”可以包含在RANS中的“湍流粘性项”中而按定常来计算。当攻角较大时,飞机后部左右双方交替出现大涡(卡曼涡街)。这时,以定常为前提的RANS方程计算就不可能收敛!除了我们一家以外,大家都用“NS方程计算”,攻角都算到105度。是否都有’特殊’的处理方法?

以上这一计算至今已快20年了。在这期间,LBM又有大发展,前景看好。虽然RANS和LBM都要用到“模式”,但前者是根据应用对象建模,而后者则是面向流体建模,通用性更强。
04
发展大型螺旋桨飞机的滑流计算的设计方法
2005年,国家决定独立发展新的水陆两栖大型运输机,主要为特种民用需要,以螺旋桨驱动。从引进到仿制再到独立设计大型运输机,我国已积累了气动力设计的经验,但螺旋桨对飞机的气动力影响,特别它所产生的滑流产生“偏航”尚无有效的算法,以致实用上需要驾驶员在起飞时依靠经验操纵来补偿。这对“自动驾驶”来说是一个尚待解决的难题。
当今计算机性能和容量的发展,对整个螺旋桨以及飞机同时用CFD整体计算已无“困难”,但计算工作量太大,很难为工程设计所采纳(有文献估计,有待计算机继续发展到2045年才有可能)。

现在已用的简化办法,是采用具有压力差的“盘面”代替。

螺旋桨-机翼模型、等效盘-机翼模型的几何及网格划分
为了反应气流经螺旋桨后的扭转,再在盘面加上沿盘面旋转的切向速度差——涡盘。这就是所谓“等效扩展盘”的模型,如今它已经被实际应用,但仅适用于( 0~5度)机翼小攻角。攻角大了,机翼对螺旋桨周围气流的影响就不能忽略。为了解决这一问题,我们将“盘面”上的均匀压力差和涡度换成“非均匀”但有一定的分布规律。规律中的参数由真实机翼表面上的流场分布来确定,即数据驱动。

这个流场分布可由实验测得,亦可以通过非简化的螺旋桨—机翼组合计算所得。我们采取后者,以此确定分布规律。然后用它选定“分区等效模型”对整机进行计算,在较大攻角时得到与真实螺旋桨飞机计算几乎相同的结果。由此,工作效果有较大的改进,而计算工作量与以前的简化方法比并无增加。

四发螺旋桨飞机带螺旋桨几何及表面网络


四发螺旋桨飞机带等效盘几何及表面网络


15度攻角时三个模型计算的升阻力系数对比
对于复杂的全机构型,为降低网格生成难度,节省网格生成时间,我们使用非结构网格进行模拟,两种模型的全机表面网格见图1和图2。分别对机翼前缘、后缘、垂尾、平尾等容易发生流动分离的位置处进行网格加密,两个模型的网格量分别为 1300万和800万。在采用32核并行计算的情况下,两个模型计算每个工况的时间分别为24小时和2.5小时。由此,得到的流场结构比较精细,可用于改善尾翼的设计。

关于这一课题的研究,我组姜哲博士的论文已在《中国科学·技术科学》(2016年,8月刊)发表。
以上是北京大学参与我国航空工业飞机气动设计计算的部分成果。其中,第二项是北大与贵航飞机设计所双方自愿并有实际成效的科技合作。其余,只能说北大是具有完成飞机空气动力学计算的实力,而不能参与国家航空工业建设的。究其原因,是国家按纵向“管理”的结果。改革开放后国家提倡“寓军于民”,但长期形成的人脉关系却不是一纸命令就能改变的。

注【1】“6912”机台集成电路中型计算机(当时命名为DJS-18),运算速度为每秒钟15万次,全部内存 384K字节。

注【2】陈耀松《机翼-机身-短舱组合体的气动力计算》,《力学与实践》1982,4(3):39~43

制版|马修


欢迎关注我们,投稿、授权等请联系

saixiansheng@zhishifenzi.com





参与讨论
0 条评论
评论
暂无评论内容
《赛先生》微信公众号创刊于2014年7月,创始人为饶毅、鲁白、谢宇三位学者,成为国内首个由知名科学家创办并担任主编的科学传播新媒体平台,共同致力于让科学文化在中国本土扎根。
订阅Newsletter

我们会定期将电子期刊发送到您的邮箱

GO