【思想与探索|武际可】计算力学的形成与发展
风流知音 【思想与探索|武际可】 谈谈工匠精神 (2020年7月28日) CFDST(2020)1028
武际可
北京大学
1 计算机的发明对力学的推动
人类研究计算工具有很久的历史，从算筹、算盘、手摇计算机、电动计算机，已经有几千年了。1945年在美国诞生的电子计算机既是计算工具的革命，又是牵动整个科学技术的大革命。
最早的电子计算机的设计方案是莫希利（J.W.Mauchly,1907－1980）提出的。第二次世界大战中，莫希利当时在宾夕法尼亚大学电工系同在附近的一所美国陆军弹道研究所共同负责一项任务，这项任务要求每天为海军计算数以千计的弹道。弹道计算是一个典型的力学课题，每条弹道以飞行60秒计，用当时的台式计算机，一个熟练的计算员要花20小时。用当时的大型模拟机也需要15分钟。即使雇用200人计算，两个月也不一定能完成大炮一天打炮的任务。莫希利为了解决这一困难，于1942年8月提出了题为《高速电子管计算装置的使用》的报告。这个报告1943年得到批准，估算投资为15万美元，所造的机器名为“电子数值积分计算机”，简称ENIAC(Electronic NumericalIntegrator and Computer)，于1943年6月开工。
研制小组的总工程师埃克特（J.P.Eckert,1919－）是一位硕士，年仅24岁。莫希利本人也只有36岁。1945年底，ENIAC宣告竣工，共耗资48万美元。1946年2月15日举行投产典礼，起初是专门用于计算弹道，后来又经过多次改进而成为通用计算机。
世界第一台通用计算机 ENIAC
计算机一旦来到世界上，便受到人们的热情关注与不断改进。先后经过了四次换代：从1945年到1958年以电子管来制造计算机为第一代，从1959年到1963年以晶体管来制造计算机为第二代，从1964年到20世纪70年代初以集成电路来制造计算机为第三代，从70年代以后用大规模集成电路制造巨型计算机为第四代。特别是，从70年代中开始的微处理机，使计算机的性能大大提高，并且由于价格便宜使计算机大为普及。据有人统计，从1945年第一台计算机诞生起，每18个月计算机的性能提高一倍、价格降低一半。
第二代晶体管计算机
自有史以来人类发明的各种工具，都是延长人的器官，如望远镜、显微镜是延长人的眼。而计算机则是延长人的脑。所以人们又把计算机称为“电脑”。计算机从它的原理、设计、制造与应用已经形成了一个新的庞大的学科群，这就是计算机科学。
在20世纪初，英国著名的力学家乐甫在他的名著《数学弹性理论》一开头，总结力学发展的规律时说：定理越来越少，计算越来越繁。意思是说一些有比较狭义意义的定理被一些更广泛意义的定理所包含，而计算公式越来越复杂。因而力学研究极大的困难在于计算太慢。计算工具太慢就成为力学研究与发展的瓶颈。美国人发明电子计算机的初衷也正是为了解决计算弹道这个典型复杂的力学问题开始的。计算机的产生使力学学科发生了巨大变化。结构分析、弹道计算、空气动力计算、数值天气预报、渗流与地下水的运动规律、天体力学中的轨道计算等等越来越多、越复杂的问题都可以交给计算机计算了。
计算机产生后，力学学科的研究手段，从只有理论、实验，增加为理论、实验与计算三种手段。计算机的强大威力淘汰了一些不适应计算机的过时方法，适应计算机的特点发展了新的计算方法，在计算机的帮助下发现了许多新现象，如奇怪吸引子与混沌就是在计算机的帮助下发现的。
计算力学这一名词的出现大约是20世纪50年代末的事情。它是研究借助计算机求解力学问题、探索力学规律、处理力学数据的新学科。计算力学又是力学、数学、计算机科学的交叉学科。
2 有限元法的提出与普及
在计算机发明后的早期，用计算机求解力学问题或别的问题仅仅利用了计算机快这一优点。紧接着而来的问题是程序工作量不能适应计算机的高速度。一台计算机需要数以百计的工作人员编程序才能喂饱。于是编写程序又成了合理使用计算机的瓶颈。人们想出了许多方法去解决这一困难。从20世纪50年代先后出现的符号汇编语言、FORTRAN语言、ALGOL语言等以及随之而迅速发展起来的软件产业，就是为解决这一问题应运而生的。
在适应于计算机求解力学问题节约程序人力方面，最成功的就是有限元方法的产生与发展。它的产生也是计算力学作为力学一个独立的分支学科形成的标志。
有限元法的思想尽管可以追溯得更早，如有人说有限元的思想是20世纪40年代美国人库朗（R. Courant）在1943年提出来的，有人说有限元是加拿大人辛格（J.L. Synge）在40年代提出来的，更有人说有限元是欧拉的折线法就包含的，还有人说在东汉刘徽的割圆术就是有限元法，不一而足。当然这些说法也不是完全没有道理。因为有限元法的思想的确是有一部分同上述人的工作有点联系。但是要知道，有限元法是同计算机紧紧相联系的。
事实是，在20世纪50年代中期世界各国都有一批人在思考用计算机求解结构力学与连续介质力学问题。如曾经在英、德工作过的希腊人阿吉里斯（J.H. Argyris）1956年、美国的特纳（M.J. Turner）、克拉夫（R.W. Clough）与马丁（H.C. Mardin）在1956年、苏联的符拉索夫（В.З.Власов）在20世纪50年代、中国的冯康在60年代初都提出了帽子函数插值或单元刚度的矩阵表示。所以很难说有限元的思想是哪一个人的发明，它是一种世界性思潮的产物。
不过在有限元法的发展历史上的重要事件是，20世纪50年代末加里福尼亚大学伯克利分校的威耳孙（E.L. Wilson,1930－）在克劳夫指导下的博士论文《二维结构的有限元分析》，该论文于1963年完成了世界上第一个解决平面弹性力学问题的通用程序。这个程序的主旨是借助于它解算任何平面弹性力学问题不需再编程序了，只要按说明输入必要的描述问题的几何、材料、荷载数据，机器就可以进行计算，并且按照要求输出计算结果。有限元法的程序一经投产，立刻显出它的无比优越性，原来在弹性力学领域内对付平面问题，只有复变函数方法与平面光弹性方法两手，这两种方法在有限元法的对比下便渐渐退出了历史舞台。威耳孙在有限元程序系统方面后来还进行过许多有意义的研究，他编写了有限元的多种单元的程序SAP（StructuralAnalysis Program），在他的指导下，他的研究生编写了非线性结构分析程序NONSAP，1981年他还最早编写了适应微处理机的程序SAP81。SAP程序经曲圣年、邓成光、吴良芝等移植与修正、SAP81程序经袁明武扩充改造形成独立的版本SAP84，这两个程序在我国工程建设中发挥了重大作用。NONSAP经过美国巴特（Bathe）的改进形成有世界影响的非线性分析程序ADINA。
随后，结构分析的有限元软件迅速发展。包含二维元、三维元、梁单元、杆单元、板单元、壳单元、流体单元等多种单元、能解决弹性、塑性、流变、流体以及温度场、电磁场各种复杂耦合问题的软件以及软件系统不断出现。在10多年内生产与销售有限元软件形成了有相当规模的社会新产业，而且使用有限元法解决实际问题迅速在工程技术部门普及。
1960年克劳夫在匹兹堡举行的美国土木学会电子计算会议上的《平面应力分析中的有限元法》是最早提到有限元的论文。之后有限元的论文、文集、专著大量涌现，专题学术会议不断召开。新的单元、新的求解器不断提出，先后有等参元、高次元、不协调元、拟协调元、杂交元、样条元、边界元、罚单元等不同的单元，有带宽与变带宽消去法、超矩阵法、波前法、子结构法、子空间迭代法等求解方法，还有网格自动剖分等前后处理的研究，这些工作大大加强了有限元法的解题能力，使有限元方法逐渐趋于成熟。1988年出版的《有限元法手册》是有限元法发展的一个阶段总结。[1]
3 非线性计算力学
计算力学的迅速发展，以及为他所取得的成功所鼓舞，使得一些学者对于计算力学的成就产生了过分乐观的估计。例如在20年前美国就有人说，再过10年风洞就要被计算机代替，20年过去了，计算机还不能取代任何风洞。计算力学所取得的成就，大体上说，对于可以用线性理论来近似的那些问题，靠计算机大部可以较好地解决了，可是对于实质上是非线性的那种力学问题，目前计算机几乎还是无能为力的。
钱学森先生说，力学“是一门用计算机计算去回答一切宏观的实际科学技术问题，计算方法非常重要；另一个辅助手段是巧妙的实验。”[1]如果说，目前在宏观力学问题中线性问题有百分之九十的可以依靠计算机来求解，百分之十靠实验求解，那么在非线性领域内，情形正好反过来。所以自从计算力学这个学科产生以来，它努力的方向就有两个方面。一方面对于线性问题，主要是扩大它求解问题的规模；另一方面，对于非线性问题来说，在努力寻求计算方法。近年来非线性问题的求解已经成为计算力学学科发展的主攻方向了。现在看来钱学森先生的看法对于宏观问题中的线性问题，已经是一种现实，而对于宏观问题中的非线性问题，这只能当作计算力学这一学科的努力方向，我们还必须准备走很长的路。
钱学森像
从20世纪60年代开始，在结构分析的有限元程序中，逐渐计入非线性项。例如讨论结构材料的塑性性质的，称为物理非线性问题，讨论结构的大变形引起的修正，称为几何非线性问题。最初的计算方案都是采用荷载增量法，即逐步给荷载一个小的增量，求相应的变形增量。大约从20世纪60年代末，人们在实际解题中发现有的问题在荷载达到极大值时计算机总是溢出而停机。这个问题困惑了人们许多年，直到20世纪70年代末80年代初才解决。1971年美国学者温泊纳（G.A.Wempner）、1978年荷兰学者瑞克斯（E. Riks）分别从理论上提出解决这个问题的方法，20世纪80年代初人们在程序上实现了这个方法。这个方法后来被称为弧长法。
结构的优化设计是计算力学中一个重要的非线性研究领域，它的主要目的是在满足一系列条件下（这些条件也被称为约束）寻求结构最优参数。通常这类问题是非线性的，而且计算量非常大，只有靠计算机的帮助才能解决。在钱令希（1916－）教授的大力提倡、组织与推动下，大连理工大学的程耿东、钟万勰得到了一些重要结果，结构优化的研究在我国有很好的发展。
求解非线性问题紧接着而来的是遇到分叉的问题。在有限元的通用程序中，对于结构稳定性的问题，通常是将问题化归于一个特征值问题，它的基础还是线性理论。在用非线性程序来求解时，往往由于遇到分叉而不能前进。这是因为在分叉点结构的总体刚度矩阵退化问题无法继续求解。为了克服这一困难，对于高维系统中的平衡解的静分叉以及霍普夫分叉，人们又发展了一系列的方法，但是在实践上还不能说已经彻底解决了。这方面的总结可参阅武际可与苏先樾著的《弹性系统的稳定性》一书（科学出版社，1994年）。关于高维系统的同宿轨道与异宿轨道的计算，以及高维系统向混沌转化的计算，迄今仍是难题。
[2]钱学森，《我对今日力学的认识》，载中国力学学会主编，《力学，迎接21世纪新的挑战》，北京理工大学出版社，1997年，第41页.
附记：国内有限单元法应用与推广的一些早期重要事件
1964年夏，中国科学院计算研究所三室承担的刘家峡水坝应力分析得到满意的结果，计算中采用了差分变分法。
1965年中国科学院计算研究所冯康在《应用数学与计算数学》上发表《基于变分原理的差分格式》，最早给出有限单元法求解平面二阶椭圆型方程问题收敛性的证明。
1772年北京大学曲圣年在国产108机上用汇编程序编制了求解弹性力学平面问题的有限元通用程序，命名为BD程序。
1973年武际可王大钧等编制了用有限单元法求解弹性旋转薄壳的静力分析程序BS和动力分析程序BSD。这两个程序解决了当时电力系统大型冷却塔应力分析问题。为此获1988年国家教委科技进步二等奖和国家科技进步三等奖。
1973年-1975年，北京大学的曲圣年用汇编语言编制的含20节点的三维求解弹性力学问题的YD程序。1975年曲圣年与殷有泉用YD程序与BD程序指导大学生对碧口水坝进行了计算。
1974年河海大学（原华东水利学院）徐芝纶编写并由水利电力出版社出版中国第一部系统讲述有限单元法的专著《弹性力学问题的有限单元法》。
1975年大连工学院的钟万勰编制有限元法分析程序JIGFIX投入使用，以后多次补充改进形成新的版本，以适应海洋平台、高层建筑、结构优化、抗震等多方面的需要。
1979年开始，北京大学曲圣年、邓成光、吴良芝等，对美国加利福尼亚州立大学伯克利分校威尔逊开发的结构分析SAP系列进行了消化、移植和推广工作。先后经过近20年的工作，举办培训班多期，为数以百计的重要工程解决了计算与应力分析问题。该项目获1985年国家科技进步二等奖。
1979年，飞机结构强度研究所的冯钟越主持编制的航空结构分析系统（HAJIFI）通过鉴定它是航空静力结构分析应用软件系统。后来1981年结构分析II型（HAJIFII）是用于动力分析的软件系统，投入应用。1985年结构分析系统III型（HAJIFIII）用于非线性分析，通过鉴定投入使用。这些软件获1985年国家科技进步二等奖。
1984年北京大学袁明武将加利福尼亚州立大学伯克利分校威尔逊开发的在微机上运行的结构分析通用程序SAP81扩充改进为SAP84，为国内大量工程计算应用。它是国内结构分析首个微机版程序。
1988年，中国建筑科学院推出专门用于建筑结构分析的软件PKPM，这个软件是继承该院从1970年代开始编制结构分析程序，不断积累与改进的结果。
