邊界積分方程再認(rèn)識及其軟件研發(fā)

張見明

序言：值杜慶華院士誕辰一百周年之際，感念與杜先生交往的點滴及杜先生對邊界元法研究的貢獻(xiàn)，對比有限差分法、有限元法和邊界元法各自適用問題的特點，呼吁學(xué)者們勇敢面對邊界元法及其自主CAD/CAE軟件研發(fā)勞動強度大、研究經(jīng)費難的現(xiàn)狀，直面邊界元法被日漸邊緣化的困境，遇山鉆洞、遇塹搭橋，繼承和發(fā)揚杜先生開創(chuàng)的邊界元法，并以此紀(jì)念杜先生！

值此杜慶華院士誕辰一百周年之際，我的導(dǎo)師姚振漢教授建議我寫一篇紀(jì)念短文。盡管研究杜先生倡導(dǎo)的邊界元法已有二十多年，可是我與杜先生的接觸不多，所以就寫一寫我與邊界元法的淵源以及在邊界元研究中的體會，以此來紀(jì)念杜先生。

我與邊界元法的緣分可追溯到20世紀(jì)80年代末，那時我在西安交通大學(xué)本科學(xué)習(xí)。由于杜先生的倡導(dǎo)和支持，西安交通大學(xué)的許多老師在研究邊界元法，包括樓志文教授、樂美峰教授以及當(dāng)時的系主任嵇醒教授等。當(dāng)時，在本科課程里有邊界元法的簡介，所以我從本科起就知道了杜先生的名字，印象中杜先生是一位高山仰止的大人物。記得那時我就被邊界積分方程的優(yōu)美所吸引，產(chǎn)生了幾點幼稚的想法，寫了一篇小短文，通過我的本科畢業(yè)設(shè)計指導(dǎo)老師李中華教授給樂美峰教授看，樂教授還以為是研究生寫的。由于種種原因，本科畢業(yè)后我沒有繼續(xù)讀研究生，而選擇去企業(yè)工作，這樣，我與邊界元法的緣分就中斷了。

到20世紀(jì)90年代末，我投身到姚振漢教授門下研究邊界元。親眼看到我仰慕的大學(xué)者時，杜先生已經(jīng)八十多歲了。雖然有點老態(tài)龍鐘，但大學(xué)者的風(fēng)范猶存。有一次，杜先生去上海交通大學(xué)訪問，姚老師指派我陪同，終于有機會與杜先生第一次親密接觸，發(fā)現(xiàn)杜先生是一位和藹可親的慈祥老人。后來，我去日本做博士后研究，杜先生為我寫了推薦信，還特意向我介紹了我在日本的合作者田中教授，提醒我在日本生活的注意事項。至今想來，仍心存感激。

在杜先生主導(dǎo)國內(nèi)邊界元法研究的一二十年間，可以說是邊界元法研究最輝煌的時期，我卻錯過了。拉格朗日遺憾沒有生在牛頓的時代，而我生在正確的年代，卻錯過了極好的機緣。當(dāng)十年后我返回高校學(xué)習(xí)、研究邊界元法時，邊界元法已經(jīng)開始衰落了。由于姚老師的倡導(dǎo)，研究快速多極算法在國內(nèi)曾掀起一陣熱潮。但是，整個大背景就要落幕，這一熱潮恐怕也只能算是最后的一點“余輝”罷。當(dāng)我從日本回國專門從事邊界元軟件研發(fā)時，邊界元法更是日薄西山，基本上被邊緣化了——獲得研究經(jīng)費愈加困難，文章不易發(fā)表。

我始終相信，邊界元法是一種科學(xué)的方法。有限元法是一種工程的方法，而科學(xué)的方法應(yīng)該是工程問題的最終解決方案。工程方法容易實現(xiàn)、效率高，因為其計算模型對真實結(jié)構(gòu)進(jìn)行大量抽象簡化。但是，一旦對結(jié)構(gòu)進(jìn)行幾何和拓?fù)涓淖?，網(wǎng)格自動劃分就變得極為困難。這就是幾十年來，盡管有實力雄厚的國際大公司的大量人力和財力投入，但CAE/CAD一體化研究至今沒有真正實現(xiàn)的原因。另外，有限元方程是原微分方程的一種近似（弱形式），有限元法得到的解不能滿足原微分方程。所以，其精度，特別是結(jié)構(gòu)危險部位的應(yīng)力精度，通常不能令人非常滿意。邊界積分方程可與原微分方程完全等價，而且對試函數(shù)的連續(xù)性要求比有限元法還要低。邊界積分方程的基本解其實就是物理定律在位勢理論中的表達(dá)式。工程師也許不欣賞基本解的奇異性，但奇異性是所有物理定律的共同特性，比如黑洞就是廣義相對論方程的一個奇異解。愛丁頓不相信黑洞存在，以致于錢德拉塞卡博士后出站時在英國找不到工作，然而30年后錢德拉塞卡獲得了諾貝爾獎。所以，正因為奇異基本解的存在，邊界元法才可以稱得上是科學(xué)的方法。

求解邊值問題有三大類數(shù)值方法：有限差分法、有限元法和邊界元法，三者進(jìn)行對比見表1。

由上面的對比可以看到，邊界元法是一個近似“圓滿”的方法。然而，如此“優(yōu)雅完美”的理論方法卻沒有產(chǎn)生一定規(guī)模的工程應(yīng)用，迅速地面臨被邊緣化的尷尬境地。每當(dāng)翻看杜先生編寫的《邊界積分方程方法——邊界元法》，我就想到前賢們做了那么多精深細(xì)致的工作，難道這些工作就要被淹沒、被淘汰、被遺忘？！每念及此，心沉鼻酸?！肮聵?biāo)傲世偕誰隱，一樣開花為底遲？”Heaviside曾說“Logic can be patient for it is eternal”。我相信，“邊界積分方程可以等待，因為它是科學(xué)的?！碑?dāng)今世界，技術(shù)雖然日新月異，但基礎(chǔ)理論幾乎不變。所以，基于科學(xué)方法開發(fā)的軟件，必定會經(jīng)久不衰。

邊界積分方程理論（即高等數(shù)學(xué)里的第二格林公式）的提出，最早可追溯到19世紀(jì)，比有限元理論早100多年。邊界元法的研究幾乎與有限元法同時起步，兩者“并駕齊驅(qū)”。有限元軟件越做越大越強，而邊界元軟件在初期雖有幾款，但不久都日漸式微、曇花一現(xiàn)。個中原因，我試想有以下幾種：

（1）早期計算機硬件不夠發(fā)達(dá)，計算機資源十分有限，算法的效率和內(nèi)存需求成為最突出的問題。

（2）邊界元法研究者沒有認(rèn)識到邊界積分方程的本質(zhì)優(yōu)點，錯誤地以為其最大優(yōu)點是降維（只劃分表面網(wǎng)格）。對于非均質(zhì)非線性問題，花費大量精力用于避免體積分，甚至在避免奇異積分方面也浪費了不少精力。雖然這些研究得到了一些數(shù)學(xué)形式上漂亮的結(jié)果，但最終都是徒勞。

（3）邊界元法的軟件實現(xiàn)套用有限元法的網(wǎng)格，對結(jié)構(gòu)的邊界變量描述不完備，也不能逾越有限元法在CAE/CAD一體化過程中面臨的鴻溝。

現(xiàn)今，邊界元軟件開發(fā)環(huán)境得到了極大改善。首先，計算機硬件已經(jīng)發(fā)展到相當(dāng)水平，足以適應(yīng)基于不加簡化的原始科學(xué)算法進(jìn)行一定規(guī)模的計算;其次，也特別重要的是，快速算法的出現(xiàn)將邊界元法的計算復(fù)雜度降低到線性或幾乎線性量級。在強調(diào)自動化、智能化的今天，針對有限元軟件在CAE/CAD一體化中難以克服的困難，利用邊界積分方程與生俱來的、與CAD完全無縫銜接的天然優(yōu)勢，開發(fā)一套CAD與CAE完全融合的計算力學(xué)軟件，應(yīng)該是振興邊界元法的一個重要契機。“蒼龍日暮還行雨，老樹春深更著花?！?/p>

十多年來，在“缺米愁柴”的境況中，我夜以繼日地工作，即使身體不好，也仍然不分節(jié)假日，只要有生存下去的條件，就潛心于程序研發(fā)中，心無旁騖。繞過九九八十一道彎，終于搭建成了一個完全融于CAD環(huán)境中的CAE分析軟件框架。該軟件框架全部用C++編寫，目前已有90個項目，近2 000個源文件，涉及4 714個類，共約80萬行代碼。在這個“舉世紛紛寶魚目，投人慎勿以明珠”的時代，僅憑一己之力，在短期內(nèi)從最基礎(chǔ)的數(shù)據(jù)結(jié)構(gòu)做起，開發(fā)一套與已經(jīng)發(fā)展了幾十年的商業(yè)有限元軟件相抗衡的分析軟件，談何容易！其艱難心酸，有幾首歪詩為證：

十年辛苦不尋常，此身難逃名利坊。

不尋天路攀玉桂，單向荒蠻踏冰霜。

癡僧西天求圣法，頑猴東海乞?qū)毑亍?/p>

今日握得金箍棒，莫復(fù)冷灶愁米糧。

心事總在程序中，奮不顧身十年同。

算法亦有蓬萊景，編程豈無鴉片功。

狷狂本是書生色，疏愚不合世流叢。

早知前路非坦蕩，九轉(zhuǎn)千回夢未空。

十年編程似遠(yuǎn)征，風(fēng)雨無阻費勞神。

難忍壓力幾度棄，已無退路續(xù)前程。

麓山腳下見孤影，湘江岸邊遺淚痕。

敝履柴車一瘦客，不是黃粱夢里人。

與現(xiàn)有的商業(yè)有限元軟件和邊界元軟件相比，我們研發(fā)的這個軟件框架是一個全新的計算力學(xué)軟件平臺。在研發(fā)過程中遇山鉆洞、遇塹搭橋，面對每一個難于越過的坎，沒有采用變通的辦法繞過去，而是進(jìn)行硬生生的攻堅戰(zhàn)，將他們徹底征服。算法上盡量做到自然而科學(xué)，避免投機取巧和拼湊。新的主要算法有：

（1）邊界面法。將邊界積分方程和計算機圖形學(xué)結(jié)合起來，繼承邊界元法的所有優(yōu)點，同時盡量避免結(jié)構(gòu)的幾何簡化，直接在結(jié)構(gòu)的三維實體模型上進(jìn)行應(yīng)力分析，是一種真正的等幾何分析方法。

（2）雙層插值法。將單元插值與無網(wǎng)格法相結(jié)合，統(tǒng)一連續(xù)和不連續(xù)單元，同時還盡量完備描述邊界積分方程的邊界變量。雙層插值法既可以插值非連續(xù)函數(shù)，也可以用不連續(xù)網(wǎng)格插值連續(xù)函數(shù)，使得變網(wǎng)格計算簡單易行，因而適用于非線性問題。利用簡單的二叉樹網(wǎng)格生成算法產(chǎn)生非連續(xù)網(wǎng)格，可實現(xiàn)任意復(fù)雜結(jié)構(gòu)的全自動CAE分析，徹底將CAE從網(wǎng)格劃分的牢籠中解放出來。

（3）自適應(yīng)快速多極算法和自適應(yīng)幾何交叉近似算法。實現(xiàn)分級矩陣的所有代數(shù)運算及其Out Of Core算法（包括直接求逆和LU分解），使得邊界積分方程在大規(guī)模問題上的應(yīng)用不再是難題。

（4）單元球面細(xì)分法。實現(xiàn)對任意基本解類型、任意單元形狀、任意源點位置的近奇異、奇異和超奇異積分的任意精度計算。

該軟件從原始CAD模型到CAE分析再到后處理全過程全自動，稱之為5aCAE（詳見http：//www.5aCAE.com）。用戶不需要懂得計算力學(xué)專門知識，不需要選擇單元類型，只需根據(jù)實際情況定義材料、載荷和邊界約束即可。施加約束和載荷也在CAD環(huán)境中完成，且均具體直觀地施加在幾何體上。

雖然，我以為邊界元法比有限元法更科學(xué)，但并不認(rèn)為他可以取代有限元法，正像量子力學(xué)計算中從頭算法比密度泛函更科學(xué)但不能取代密度泛函、分子動力學(xué)不能取代連續(xù)介質(zhì)力學(xué)一樣。不管是工程問題還是科學(xué)問題，在計算機資源條件允許的情況下，應(yīng)優(yōu)先選取更科學(xué)的算法。對于如機械這樣一類精細(xì)結(jié)構(gòu)的應(yīng)力分析，邊界面法可能是更好的選擇。

我以為，繼承和發(fā)揚杜先生開創(chuàng)的事業(yè)，是紀(jì)念先生最好的方式。下面用一首七律紀(jì)念杜先生，并以之結(jié)束本文：

仿真原本萬年計，廿載毋論百年功。

振臂一聲開新域，搖旗幾度聚舊朋。

多極運算遠(yuǎn)近別，雙層插值斷連同。

先生有靈應(yīng)喜甚，邊元振興紹宗風(fēng)。