边界积分方程再认识及其软件研发

张见明

序言：值杜庆华院士诞辰一百周年之际，感念与杜先生交往的点滴及杜先生对边界元法研究的贡献，对比有限差分法、有限元法和边界元法各自适用问题的特点，呼吁学者们勇敢面对边界元法及其自主CAD/CAE软件研发劳动强度大、研究经费难的现状，直面边界元法被日渐边缘化的困境，遇山钻洞、遇堑搭桥，继承和发扬杜先生开创的边界元法，并以此纪念杜先生！

值此杜庆华院士诞辰一百周年之际，我的导师姚振汉教授建议我写一篇纪念短文。尽管研究杜先生倡导的边界元法已有二十多年，可是我与杜先生的接触不多，所以就写一写我与边界元法的渊源以及在边界元研究中的体会，以此来纪念杜先生。

我与边界元法的缘分可追溯到20世纪80年代末，那时我在西安交通大学本科学习。由于杜先生的倡导和支持，西安交通大学的许多老师在研究边界元法，包括楼志文教授、乐美峰教授以及当时的系主任嵇醒教授等。当时，在本科课程里有边界元法的简介，所以我从本科起就知道了杜先生的名字，印象中杜先生是一位高山仰止的大人物。记得那时我就被边界积分方程的优美所吸引，产生了几点幼稚的想法，写了一篇小短文，通过我的本科毕业设计指导老师李中华教授给乐美峰教授看，乐教授还以为是研究生写的。由于种种原因，本科毕业后我没有继续读研究生，而选择去企业工作，这样，我与边界元法的缘分就中断了。

到20世纪90年代末，我投身到姚振汉教授门下研究边界元。亲眼看到我仰慕的大学者时，杜先生已经八十多岁了。虽然有点老态龙钟，但大学者的风范犹存。有一次，杜先生去上海交通大学访问，姚老师指派我陪同，终于有机会与杜先生第一次亲密接触，发现杜先生是一位和蔼可亲的慈祥老人。后来，我去日本做博士后研究，杜先生为我写了推荐信，还特意向我介绍了我在日本的合作者田中教授，提醒我在日本生活的注意事项。至今想来，仍心存感激。

在杜先生主导国内边界元法研究的一二十年间，可以说是边界元法研究最辉煌的时期，我却错过了。拉格朗日遗憾没有生在牛顿的时代，而我生在正确的年代，却错过了极好的机缘。当十年后我返回高校学习、研究边界元法时，边界元法已经开始衰落了。由于姚老师的倡导，研究快速多极算法在国内曾掀起一阵热潮。但是，整个大背景就要落幕，这一热潮恐怕也只能算是最后的一点“余辉”罢。当我从日本回国专门从事边界元软件研发时，边界元法更是日薄西山，基本上被边缘化了——获得研究经费愈加困难，文章不易发表。

我始终相信，边界元法是一种科学的方法。有限元法是一种工程的方法，而科学的方法应该是工程问题的最终解决方案。工程方法容易实现、效率高，因为其计算模型对真实结构进行大量抽象简化。但是，一旦对结构进行几何和拓扑改变，网格自动划分就变得极为困难。这就是几十年来，尽管有实力雄厚的国际大公司的大量人力和财力投入，但CAE/CAD一体化研究至今没有真正实现的原因。另外，有限元方程是原微分方程的一种近似（弱形式），有限元法得到的解不能满足原微分方程。所以，其精度，特别是结构危险部位的应力精度，通常不能令人非常满意。边界积分方程可与原微分方程完全等价，而且对试函数的连续性要求比有限元法还要低。边界积分方程的基本解其实就是物理定律在位势理论中的表达式。工程师也许不欣赏基本解的奇异性，但奇异性是所有物理定律的共同特性，比如黑洞就是广义相对论方程的一个奇异解。爱丁顿不相信黑洞存在，以致于钱德拉塞卡博士后出站时在英国找不到工作，然而30年后钱德拉塞卡获得了诺贝尔奖。所以，正因为奇异基本解的存在，边界元法才可以称得上是科学的方法。

求解边值问题有三大类数值方法：有限差分法、有限元法和边界元法，三者进行对比见表1。

由上面的对比可以看到，边界元法是一个近似“圆满”的方法。然而，如此“优雅完美”的理论方法却没有产生一定规模的工程应用，迅速地面临被边缘化的尴尬境地。每当翻看杜先生编写的《边界积分方程方法——边界元法》，我就想到前贤们做了那么多精深细致的工作，难道这些工作就要被淹没、被淘汰、被遗忘？！每念及此，心沉鼻酸。“孤标傲世偕谁隐，一样开花为底迟？”Heaviside曾说“Logic can be patient for it is eternal”。我相信，“边界积分方程可以等待，因为它是科学的。”当今世界，技术虽然日新月异，但基础理论几乎不变。所以，基于科学方法开发的软件，必定会经久不衰。

边界积分方程理论（即高等数学里的第二格林公式）的提出，最早可追溯到19世纪，比有限元理论早100多年。边界元法的研究几乎与有限元法同时起步，两者“并驾齐驱”。有限元软件越做越大越强，而边界元软件在初期虽有几款，但不久都日渐式微、昙花一现。个中原因，我试想有以下几种：

（1）早期计算机硬件不够发达，计算机资源十分有限，算法的效率和内存需求成为最突出的问题。

（2）边界元法研究者没有认识到边界积分方程的本质优点，错误地以为其最大优点是降维（只划分表面网格）。对于非均质非线性问题，花费大量精力用于避免体积分，甚至在避免奇异积分方面也浪费了不少精力。虽然这些研究得到了一些数学形式上漂亮的结果，但最终都是徒劳。

（3）边界元法的软件实现套用有限元法的网格，对结构的边界变量描述不完备，也不能逾越有限元法在CAE/CAD一体化过程中面临的鸿沟。

现今，边界元软件开发环境得到了极大改善。首先，计算机硬件已经发展到相当水平，足以适应基于不加简化的原始科学算法进行一定规模的计算;其次，也特别重要的是，快速算法的出现将边界元法的计算复杂度降低到线性或几乎线性量级。在强调自动化、智能化的今天，针对有限元软件在CAE/CAD一体化中难以克服的困难，利用边界积分方程与生俱来的、与CAD完全无缝衔接的天然优势，开发一套CAD与CAE完全融合的计算力学软件，应该是振兴边界元法的一个重要契机。“苍龙日暮还行雨，老树春深更著花。”

十多年来，在“缺米愁柴”的境况中，我夜以继日地工作，即使身体不好，也仍然不分节假日，只要有生存下去的条件，就潜心于程序研发中，心无旁骛。绕过九九八十一道弯，终于搭建成了一个完全融于CAD环境中的CAE分析软件框架。该软件框架全部用C++编写，目前已有90个项目，近2 000個源文件，涉及4 714个类，共约80万行代码。在这个“举世纷纷宝鱼目，投人慎勿以明珠”的时代，仅凭一己之力，在短期内从最基础的数据结构做起，开发一套与已经发展了几十年的商业有限元软件相抗衡的分析软件，谈何容易！其艰难心酸，有几首歪诗为证：

十年辛苦不寻常，此身难逃名利坊。

不寻天路攀玉桂，单向荒蛮踏冰霜。

痴僧西天求圣法，顽猴东海乞宝藏。

今日握得金箍棒，莫复冷灶愁米粮。

心事总在程序中，奋不顾身十年同。

算法亦有蓬莱景，编程岂无鸦片功。

狷狂本是书生色，疏愚不合世流丛。

早知前路非坦荡，九转千回梦未空。

十年编程似远征，风雨无阻费劳神。

难忍压力几度弃，已无退路续前程。

麓山脚下见孤影，湘江岸边遗泪痕。

敝履柴车一瘦客，不是黄粱梦里人。

与现有的商业有限元软件和边界元软件相比，我们研发的这个软件框架是一个全新的计算力学软件平台。在研发过程中遇山钻洞、遇堑搭桥，面对每一个难于越过的坎，没有采用变通的办法绕过去，而是进行硬生生的攻坚战，将他们彻底征服。算法上尽量做到自然而科学，避免投机取巧和拼湊。新的主要算法有：

（1）边界面法。将边界积分方程和计算机图形学结合起来，继承边界元法的所有优点，同时尽量避免结构的几何简化，直接在结构的三维实体模型上进行应力分析，是一种真正的等几何分析方法。

（2）双层插值法。将单元插值与无网格法相结合，统一连续和不连续单元，同时还尽量完备描述边界积分方程的边界变量。双层插值法既可以插值非连续函数，也可以用不连续网格插值连续函数，使得变网格计算简单易行，因而适用于非线性问题。利用简单的二叉树网格生成算法产生非连续网格，可实现任意复杂结构的全自动CAE分析，彻底将CAE从网格划分的牢笼中解放出来。

（3）自适应快速多极算法和自适应几何交叉近似算法。实现分级矩阵的所有代数运算及其Out Of Core算法（包括直接求逆和LU分解），使得边界积分方程在大规模问题上的应用不再是难题。

（4）单元球面细分法。实现对任意基本解类型、任意单元形状、任意源点位置的近奇异、奇异和超奇异积分的任意精度计算。

该软件从原始CAD模型到CAE分析再到后处理全过程全自动，称之为5aCAE（详见http：//www.5aCAE.com）。用户不需要懂得计算力学专门知识，不需要选择单元类型，只需根据实际情况定义材料、载荷和边界约束即可。施加约束和载荷也在CAD环境中完成，且均具体直观地施加在几何体上。

虽然，我以为边界元法比有限元法更科学，但并不认为他可以取代有限元法，正像量子力学计算中从头算法比密度泛函更科学但不能取代密度泛函、分子动力学不能取代连续介质力学一样。不管是工程问题还是科学问题，在计算机资源条件允许的情况下，应优先选取更科学的算法。对于如机械这样一类精细结构的应力分析，边界面法可能是更好的选择。

我以为，继承和发扬杜先生开创的事业，是纪念先生最好的方式。下面用一首七律纪念杜先生，并以之结束本文：

仿真原本万年计，廿载毋论百年功。

振臂一声开新域，摇旗几度聚旧朋。

多极运算远近别，双层插值断连同。

先生有灵应喜甚，边元振兴绍宗风。