陳積瞻,蘇海東,祁勇峰

（長(zhǎng)江科學(xué)院a.水利部水工程安全與病害防治工程技術(shù)研究中心；b.材料與結(jié)構(gòu)研究所,武漢 430010)

?

基于CAD幾何的數(shù)值流形方法初步研究

陳積瞻a,b,蘇海東a,b,祁勇峰a,b

（長(zhǎng)江科學(xué)院a.水利部水工程安全與病害防治工程技術(shù)研究中心；b.材料與結(jié)構(gòu)研究所,武漢 430010)

摘 要：計(jì)算機(jī)輔助工程中的設(shè)計(jì)—分析—再設(shè)計(jì)的反復(fù)過(guò)程,蘊(yùn)含著計(jì)算機(jī)輔助設(shè)計(jì)（CAD)與計(jì)算機(jī)輔助工程分析（CAE)相融合的迫切需求。提出基于CAD幾何的數(shù)值流形方法:按照CAE真實(shí)物理場(chǎng)分布的復(fù)雜程度來(lái)布置數(shù)學(xué)網(wǎng)格和設(shè)置近似函數(shù)階次,比等幾何分析方法更合理；只需引入自動(dòng)且快速的切割操作,就能實(shí)現(xiàn)CAD模型進(jìn)入CAE后無(wú)需修改而直接建模；針對(duì)以往的流形法需要將曲線(xiàn)邊界離散成折線(xiàn)的問(wèn)題,給出了曲線(xiàn)邊界與網(wǎng)格直邊的切割算法,實(shí)現(xiàn)了幾何模型在CAE建模和網(wǎng)格細(xì)化中的保形性；針對(duì)流形法通常使用的多項(xiàng)式近似函數(shù),推導(dǎo)了曲線(xiàn)“近似”單純形的解析積分公式并應(yīng)用于帶有曲線(xiàn)邊界的流形元的精確積分運(yùn)算；最后通過(guò)平板內(nèi)的圓孔算例驗(yàn)證了方法的可行性。該方法對(duì)CAD和CAE的融合提出了全新的思路,為實(shí)現(xiàn)CAD進(jìn)入CAE后的自動(dòng)化分析打下了基礎(chǔ)。

關(guān)鍵詞：數(shù)值流形方法；等幾何分析方法；CAD幾何；曲線(xiàn)與直線(xiàn)的切割；單純形積分；CAD/ CAE協(xié)同

2016,33（02):137-143

1 研究背景

隨著計(jì)算機(jī)科學(xué)與數(shù)值分析技術(shù)的迅速發(fā)展,計(jì)算機(jī)輔助工程包括計(jì)算機(jī)輔助設(shè)計(jì)（CAD)與計(jì)算機(jī)輔助工程分析（CAE)等,已成為現(xiàn)代設(shè)計(jì)的主要工具和手段［1］。應(yīng)用CAE數(shù)學(xué)模型替代以往的物理模型對(duì)設(shè)計(jì)結(jié)果進(jìn)行分析和評(píng)價(jià),可為改進(jìn)、優(yōu)化設(shè)計(jì)提供重要依據(jù)。這種設(shè)計(jì)—分析—再設(shè)計(jì)的反復(fù)過(guò)程,蘊(yùn)含著CAD設(shè)計(jì)與CAE分析相融合（又稱(chēng)為CAD和CAE協(xié)同設(shè)計(jì)與分析)的迫切需求。我們認(rèn)為,最終目標(biāo)是實(shí)現(xiàn)CAD進(jìn)入CAE后的自動(dòng)化分析。

計(jì)算分析方法和軟件是CAE的關(guān)鍵因素［1］。目前,有限元方法是應(yīng)用最為廣泛的分析方法,但有限元網(wǎng)格需要適應(yīng)材料邊界并保持單元合適形狀的硬性要求帶來(lái)了2方面的影響:其一,CAD幾何模型在進(jìn)入CAE時(shí)一般要經(jīng)過(guò)修剪、整理與簡(jiǎn)化,如去掉小的圓形倒角并進(jìn)行細(xì)節(jié)修補(bǔ)、曲線(xiàn)或曲面（本文以下部分均以曲線(xiàn)為例)邊界的分段簡(jiǎn)化（如線(xiàn)性化)等；其二,對(duì)形狀復(fù)雜的結(jié)構(gòu)進(jìn)行有限元網(wǎng)格劃分往往非常困難。因此需要大量繁瑣且非常耗時(shí)的人工操作,而最令分析人員頭疼的是設(shè)計(jì)—分析反復(fù)過(guò)程導(dǎo)致的因設(shè)計(jì)變更而重復(fù)建模。雖然一些大型商用有限元軟件如ANSYS旗下的協(xié)同分析平臺(tái)以及MSC公司最新推出的APEX平臺(tái)等,提供了一定程度的“智能化”協(xié)同分析和網(wǎng)格劃分功能,但離分析人員的預(yù)期仍有不小的距離,這是由于有限元方法自身的局限性使之在CAD和CAE融合方面存在一定的困難。

2005年,Hughes等［2］提出了等幾何分析方法,在幾何模型和分析模型中采用一套描述方式——通常是NURBS（Non-Uniform Rational B-Splines非均勻有理B樣條),解決了傳統(tǒng)數(shù)值計(jì)算的求解域描述與幾何設(shè)計(jì)之間不相兼容的問(wèn)題,實(shí)現(xiàn)了CAD和CAE的無(wú)縫連接。CAD模型可直接參與CAE分析而無(wú)需再劃分網(wǎng)格,避免了設(shè)計(jì)與分析之間過(guò)多的人工交互以及設(shè)計(jì)修改后的分析模型重建。同時(shí),該方法具有幾何模型的精確性,即保形性:在CAE建模時(shí)避免了因幾何模型簡(jiǎn)化而產(chǎn)生的誤差,且網(wǎng)格細(xì)分后仍保持原有邊界曲線(xiàn)的幾何特征［3］。另外還具有近似函數(shù)的高階連續(xù)性,計(jì)算精度相對(duì)較高。等幾何分析方法提出了CAD和CAE融合的新思路,成為國(guó)際上的研究熱點(diǎn)。

然而,等幾何分析方法也存在諸多不足:工業(yè)界通常采用以邊界數(shù)據(jù)表示的幾何模型,與等幾何分析要求的基于實(shí)體造型的幾何模型不匹配［4］,而要求工業(yè)界為了適應(yīng)某種計(jì)算方法而改變現(xiàn)狀是很難的,反之,將邊界模型轉(zhuǎn)化為實(shí)體模型,無(wú)異于做一次完整的網(wǎng)格剖分；該方法采用的樣條形函數(shù)具有多維基函數(shù)的張量積特征,形成的網(wǎng)格在拓?fù)渖鲜且?guī)整的［4］,對(duì)于復(fù)雜求解域的適應(yīng)性以及局部的網(wǎng)格細(xì)化不如有限元法方便；“等參”是等幾何分析的基礎(chǔ),而“等參”帶來(lái)的實(shí)際物理域映射到標(biāo)準(zhǔn)參數(shù)域的一些問(wèn)題,比如,網(wǎng)格扭曲（相對(duì)于規(guī)則網(wǎng)格)、曲線(xiàn)邊界（相對(duì)于直線(xiàn)邊界)所造成的多項(xiàng)式精度階次降低,數(shù)值積分精度難以把握等,也會(huì)被帶入等幾何分析方法［4］；筆者認(rèn)為該方法最主要的問(wèn)題是,“等幾何”將物理場(chǎng)的復(fù)雜性和幾何描述的復(fù)雜性等同起來(lái),然而實(shí)際情況是,幾何描述簡(jiǎn)單或復(fù)雜并不意味著物理場(chǎng)也同樣簡(jiǎn)單或復(fù)雜,因此在計(jì)算效率方面值得討論；同時(shí),等幾何分析中的NURBS基函數(shù)及其導(dǎo)數(shù)的計(jì)算遠(yuǎn)比有限元法復(fù)雜［4］。

本文基于數(shù)值流形方法,提出CAD與CAE融合的另一種新思路:根據(jù)該方法“構(gòu)造近似函數(shù)的數(shù)學(xué)網(wǎng)格與描述求解域的幾何模型相互獨(dú)立”的特點(diǎn),只需引入自動(dòng)且快速的切割操作,就能實(shí)現(xiàn)CAD幾何模型無(wú)需修改而直接進(jìn)行CAE建模,其中,采用曲線(xiàn)幾何邊界與直線(xiàn)網(wǎng)格邊界的切割,實(shí)現(xiàn)CAE建模和網(wǎng)格細(xì)化的保形性；針對(duì)多項(xiàng)式級(jí)數(shù)的近似函數(shù),采用單純形（包括具有曲線(xiàn)邊界的“近似”單純形)積分法進(jìn)行剛度矩陣和荷載向量各項(xiàng)中的被積函數(shù)的解析積分；最終目標(biāo)是采用自動(dòng)化的前處理并根據(jù)誤差估計(jì)來(lái)控制計(jì)算精度,實(shí)現(xiàn)CAD進(jìn)入CAE后的自動(dòng)化分析,即CAD和CAE的完全融合。

2 數(shù)值流形方法與CAD幾何

1991年,石根華博士首次將現(xiàn)代數(shù)學(xué)“流形”思想引入工程計(jì)算,發(fā)明了數(shù)值流形方法［5-6］（以下簡(jiǎn)稱(chēng)流形法),其中很關(guān)鍵的一點(diǎn)是:與有限元法將求解域離散成網(wǎng)格的方式不同,流形法構(gòu)造物理場(chǎng)近似解的所謂“數(shù)學(xué)網(wǎng)格”與實(shí)際求解域分離（求解域只用于定義積分區(qū)域),要求數(shù)學(xué)網(wǎng)格在空間上完全覆蓋求解域。如圖1所示,圖中著色的橢圓形求解域被矩形數(shù)學(xué)網(wǎng)格完全覆蓋。

對(duì)應(yīng)于數(shù)學(xué)流形中的“覆蓋”,數(shù)學(xué)網(wǎng)格又被稱(chēng)為數(shù)學(xué)覆蓋。數(shù)值計(jì)算中,在各覆蓋上定義覆蓋函數(shù)Vi,通過(guò)單位分解函數(shù)φi聯(lián)系起來(lái)形成整體近似函數(shù)（n為覆蓋數(shù))。覆蓋函數(shù)Vi常用多項(xiàng)式級(jí)數(shù)。

圖1 矩形數(shù)學(xué)覆蓋與求解域示意圖Fig.1 Schematic diagrams of rectangular mathematical covers and solving domain

從數(shù)學(xué)網(wǎng)格（數(shù)學(xué)覆蓋)的形式上看,流形法的研究分為2大類(lèi):

（1)傳統(tǒng)流形法的數(shù)學(xué)網(wǎng)格大都采用有限元網(wǎng)格［7］,如圖1（a)所示的矩形網(wǎng)格,這種覆蓋形式稱(chēng)為完全重疊的數(shù)學(xué)覆蓋［8］。在有限元結(jié)點(diǎn)上定義級(jí)數(shù)作為覆蓋函數(shù),單位分解函數(shù)就是有限單元的形函數(shù)。其網(wǎng)格細(xì)化仍要采用有限元網(wǎng)格的加密方式。

（2)非有限元的數(shù)學(xué)網(wǎng)格,目前主要是文獻(xiàn)［8］提出的部分重疊覆蓋的數(shù)學(xué)網(wǎng)格,引入了“獨(dú)立覆蓋”,即單位分解函數(shù)φi＝1、近似函數(shù)V就是給定級(jí)數(shù)Vi的覆蓋區(qū)域。如圖1（b)所示,圖中的大矩形就是獨(dú)立覆蓋,各覆蓋僅在條形區(qū)域重疊。建議條形厚度取較小值,其單位分解函數(shù)φi取為有限單元形函數(shù)（如2個(gè)覆蓋之間的條形為一維線(xiàn)性函數(shù))實(shí)現(xiàn)覆蓋之間的線(xiàn)性過(guò)渡。該方法又被稱(chēng)為基于獨(dú)立覆蓋的數(shù)值流形方法［9］,其一大優(yōu)勢(shì)是覆蓋加密方便,如圖2所示,在大的矩形覆蓋中劃分小覆蓋,可保證近似函數(shù)的協(xié)調(diào)性,已被證明能有效提高精度［10］。

圖2 矩形獨(dú)立覆蓋的加密Fig.2 Refinement of rectangular independent covers

有別于等幾何分析法的CAE“等幾何”特性,流形法具有CAE數(shù)學(xué)網(wǎng)格與CAD幾何模型的相互獨(dú)立性,根據(jù)CAE真實(shí)物理場(chǎng)分布的復(fù)雜程度來(lái)布置數(shù)學(xué)網(wǎng)格和設(shè)置近似函數(shù)階次,這種分析方式更為合理。但等幾何分析法可以直接利用CAD實(shí)體幾何模型的網(wǎng)格作為CAE分析模型的網(wǎng)格,而流形法要單獨(dú)生成數(shù)學(xué)網(wǎng)格,通常采用如圖1所示的簡(jiǎn)單且可自動(dòng)生成的數(shù)學(xué)網(wǎng)格,需要付出數(shù)學(xué)網(wǎng)格與幾何邊界進(jìn)行切割的代價(jià)?？紤]到這種切割操作只涉及簡(jiǎn)單的幾何運(yùn)算,能夠自動(dòng)且快速地完成,因此流形法的前處理同樣比有限元法簡(jiǎn)單得多,可做到完全的自動(dòng)化。

如圖1所示,經(jīng)過(guò)切割操作后,在數(shù)學(xué)網(wǎng)格內(nèi)有可能形成不能占滿(mǎn)整個(gè)網(wǎng)格的物理區(qū)域（數(shù)學(xué)網(wǎng)格內(nèi)的物理區(qū)域稱(chēng)為流形元),以圖1中的矩形網(wǎng)格c為例,如圖3所示,由頂點(diǎn)P1,P2,P3,P4,P5,P6,P7所組成的流形元具有不規(guī)則的形狀,其中P1至P5是由離散的折線(xiàn)來(lái)近似邊界曲線(xiàn)。對(duì)于任意形狀的流形元,當(dāng)近似函數(shù)為多項(xiàng)式級(jí)數(shù)時(shí),通常采用單純形積分公式［6］進(jìn)行解析積分。如圖3所示,將相鄰頂點(diǎn)構(gòu)成的線(xiàn)段依次與某一外點(diǎn)P0（通常取為坐標(biāo)原點(diǎn))相連形成一個(gè)有向單純形,可取逆時(shí)針頂點(diǎn)順序?yàn)檎?如P2,P1,P0等為正向,而P7,P6,P0和P1,P7,P0為負(fù)向。這樣,對(duì)于剛度矩陣和荷載向量各項(xiàng)中的被積函數(shù)經(jīng)整理后［11］形成的各單項(xiàng)式xm1ym2（m1, m2為正整數(shù)),在各單純形中可得到精確的積分,考慮積分區(qū)域的正、負(fù)向后,累加得到整個(gè)流形元上的精確積分。

圖3 單純形積分示意圖Fig.3 Schematic diagram of simplex integration

當(dāng)CAD幾何模型進(jìn)入CAE后,無(wú)需像有限元法那樣進(jìn)行“去倒角”等細(xì)節(jié)修補(bǔ)工作,因?yàn)榈菇鞘菐缀芜吔绲囊徊糠?它與CAE數(shù)學(xué)網(wǎng)格沒(méi)有關(guān)系,可以保留作為網(wǎng)格內(nèi)的流形元邊界的一部分。

圖4 CAD曲線(xiàn)邊界的線(xiàn)段化Fig.4 Process of changing CAD curve boundaries to line segments

與等幾何分析方法相比,現(xiàn)有的流形法仍然需要將曲線(xiàn)邊界線(xiàn)性化（線(xiàn)段化),如圖4所示。這會(huì)帶來(lái)幾個(gè)不利影響:首先,在CAE分析中會(huì)造成幾何模型誤差,且網(wǎng)格細(xì)分時(shí)面對(duì)的是已經(jīng)線(xiàn)性化的幾何邊界；其次,在CAE分析前就要預(yù)估線(xiàn)性化的程度（如將曲線(xiàn)劃分為多少直線(xiàn)段),這需要分析人員具有一定的經(jīng)驗(yàn)；同時(shí),對(duì)于每條線(xiàn)段都要進(jìn)行一次單純形積分,若為減小幾何誤差而劃分較多線(xiàn)段,可能會(huì)影響計(jì)算效率,特別是對(duì)于小的圓形倒角情況。以上問(wèn)題造成了流形法在CAD和CAE系統(tǒng)之間的不匹配。

因此,本文除了以上的論述——分析流形法在CAD和CAE融合中的特色外,下面的主要工作就是要做到保形性,即在分析模型構(gòu)造和網(wǎng)格細(xì)化以及計(jì)算分析過(guò)程中保持原有的幾何形狀。本文以平面計(jì)算為例,僅需兩項(xiàng)工作:曲線(xiàn)幾何邊界與直線(xiàn)網(wǎng)格邊界的切割；曲線(xiàn)邊界“近似”單純形的解析積分。

3 曲線(xiàn)幾何邊界與直線(xiàn)網(wǎng)格邊界的切割

數(shù)學(xué)網(wǎng)格和幾何模型邊界的切割操作,只涉及簡(jiǎn)單幾何運(yùn)算。目前的流形法基于直線(xiàn)段的切割,已有成熟算法和程序,主要過(guò)程如下［6］:①將數(shù)學(xué)網(wǎng)格和幾何邊界線(xiàn)段化,每個(gè)數(shù)學(xué)網(wǎng)格的邊和每條描述邊界的折線(xiàn)都各自作為一條線(xiàn)段；②所有線(xiàn)段之間相互切割；③根據(jù)每條線(xiàn)段上記錄的交點(diǎn)位置順序及每個(gè)交點(diǎn)對(duì)應(yīng)的線(xiàn)段序號(hào),按照一定規(guī)則確定環(huán)路,同時(shí)刪除不能形成環(huán)路的“樹(shù)枝”,在求解域內(nèi)的每個(gè)環(huán)路就是一個(gè)流形元；④根據(jù)流形元的形心判斷它屬于哪個(gè)數(shù)學(xué)網(wǎng)格。

圖5 坐標(biāo)轉(zhuǎn)換Fig.5 Coordinate transformation

本文中,過(guò)程①的曲線(xiàn)幾何邊界無(wú)需變成折線(xiàn),直接進(jìn)入過(guò)程②與直線(xiàn)網(wǎng)格邊界進(jìn)行切割。切割算法如下:如圖5所示,設(shè)曲線(xiàn)n可表示為

求其函數(shù)與直線(xiàn)段A （xa,ya)—B（xb,yb)的交點(diǎn)。

將整體坐標(biāo)系轉(zhuǎn)換到以A（xa,ya)為原點(diǎn)、AB方向?yàn)閤t軸的局部坐標(biāo)系（xt,yt)

式中θ為AB線(xiàn)段與整體坐標(biāo)x軸的夾角。則式（1)為

設(shè)局部坐標(biāo)下的交點(diǎn)為（xt,0) ,令xm＝xtcosθ＋xa,則式（3)為

利用二分法求得xm的值,即可求出xt,再利用式（2)得到整體坐標(biāo)系下的交點(diǎn)坐標(biāo)（x,y)。在多個(gè)交點(diǎn)情況下可參考文獻(xiàn)［12］進(jìn)行求解。

對(duì)于圓、橢圓曲線(xiàn),雖然整體上不能簡(jiǎn)單表示成式（1)的函數(shù)形式,但可以用分段函數(shù)表示?；蛘卟捎闷渌绞?如圖6所示的以O(shè)為圓點(diǎn)的圓弧與直線(xiàn)MN相交,記交點(diǎn)為J（x,y) ,則OM大于半徑,ON小于半徑,按二分法迭代求解,當(dāng)OJ等于半徑時(shí)即為交點(diǎn)。

可見(jiàn),曲線(xiàn)與直線(xiàn)求交方式多樣化,是由曲線(xiàn)種類(lèi)的多樣性決定的。在CAD中,可將一些復(fù)雜曲線(xiàn)統(tǒng)一由NURBS表示,而NURBS與直線(xiàn)求交的算法已很成熟,限于篇幅不再詳述。

圖6 圓弧與直線(xiàn)段相交Fig.6 A circle arc intersected by a line segment

當(dāng)然,曲線(xiàn)與直線(xiàn)求交的運(yùn)算量要遠(yuǎn)大于直線(xiàn)之間求交的運(yùn)算量,因此為提高效率,事先要排除無(wú)關(guān)的直線(xiàn)段,如圖7（a)所示的曲線(xiàn)n ,采用包圍盒子法［13］排除盒子（圖中虛線(xiàn))外的無(wú)關(guān)線(xiàn)段,再區(qū)分圖7（b)的幾種直線(xiàn)情況:首先按函數(shù)值排除法,如直線(xiàn)段5和6的首尾端點(diǎn)的函數(shù)值均在曲線(xiàn)上方或下方；然后,直線(xiàn)段2和3的某一端點(diǎn)在曲線(xiàn)上,則交點(diǎn)為其端點(diǎn)；豎直的線(xiàn)段4,將x坐標(biāo)代入曲線(xiàn)方程即可求得y坐標(biāo)；最后才用圖5所示的坐標(biāo)轉(zhuǎn)換及二分法處理直線(xiàn)段1。

圖7 包圍盒子法及交點(diǎn)的多種情況Fig.7 Bounding box method and some circumstances of intersection points

4 曲線(xiàn)邊界“近似”單純形的解析積分

仍以圖1中的矩形網(wǎng)格c為例,與圖3中將流形元曲線(xiàn)邊界轉(zhuǎn)化為多條折線(xiàn)的方式不同,如圖8（a)所示,P1到P5仍用原始曲線(xiàn)表示,將其與P0相連,構(gòu)成一條邊為曲線(xiàn)的“近似”單純形,如圖8 （b)所示,其解析積分分為2部分:一部分是直線(xiàn)邊的單純形,按常規(guī)單純形公式［14］計(jì)算；另一部分是P1PkP5,其積分下限為P1到P5的直線(xiàn),積分上限為P1到P5的曲線(xiàn)。以下討論后一部分。

圖8 含一條曲線(xiàn)邊的“近似”單純形Fig.8 An approximate simplex containing a curved edge

設(shè)作為積分上限的曲線(xiàn)表述為多項(xiàng)式函數(shù)

m為多項(xiàng)式的最高階次,不能寫(xiě)成函數(shù)形式時(shí)可分段表示為函數(shù)。作為積分下限的直線(xiàn)為

其中, ai, bi為已知系數(shù)。曲線(xiàn)兩端點(diǎn)為（x1,y1)、（x2,y2) ,被積函數(shù)f（x,y)為xm1ym2,則

記:L＝1/（m2＋1)；Am＝[a0a1…am]T；Xm＝[x0x1…xm]T；B＝[b0b1]T；m＝m＋1,化簡(jiǎn)

102為

關(guān)鍵是求出下式的精確解,

展開(kāi)式中仍有同類(lèi)項(xiàng)需要合并,若各項(xiàng)的冪次0·N0＋1·N1＋2·N2…＋m·Nm相等,則標(biāo)記為同類(lèi)項(xiàng),合并后一共有m·m0＋1項(xiàng),記各項(xiàng)系數(shù)為Kj（j是從0至m·m0＋1遞增的整數(shù)),則式（9)可寫(xiě)為

可見(jiàn),隨j增大的各項(xiàng)按照x的冪差為1的次序遞增,為編程提供了良好的遞推規(guī)律,即后項(xiàng)比前項(xiàng)多乘一個(gè)x1和x2。

上述方法理論上對(duì)于積分上限為任意多項(xiàng)式的曲線(xiàn)均可以計(jì)算,但如果描述曲線(xiàn)的多項(xiàng)式階數(shù)過(guò)高,則計(jì)算量較大,且容易引入數(shù)值誤差。考慮到一般計(jì)算中的數(shù)學(xué)網(wǎng)格不至于很大,網(wǎng)格內(nèi)的曲線(xiàn)也不會(huì)過(guò)于復(fù)雜,2階或3階多項(xiàng)式就可以很好地描述,因此,建議首先用2階或3階多項(xiàng)式逼近網(wǎng)格內(nèi)的復(fù)雜曲線(xiàn)（可以是任意表達(dá)式,如NURBS,還包括高于3階的多項(xiàng)式),再利用上述公式進(jìn)行解析積分。

對(duì)于曲線(xiàn)邊界l上的分布力荷載,設(shè)分布力p＝[Px,Py]T,其荷載向量中的被積函數(shù)也分解成單項(xiàng)式xm1ym2的形式［11］?？紤]l為二次曲線(xiàn)y＝∑i＝20aixi,導(dǎo)數(shù)為y′＝a1＋2a2x。若計(jì)算的正壓力荷載為線(xiàn)性分布P＝c0＋c1x＋c2y ,則Px＝Pcosθ,Py＝Psinθ（θ為曲線(xiàn)上某點(diǎn)的法線(xiàn)與坐標(biāo)x軸的夾角),

其最終仍可歸結(jié)為求式（9)中R的解析積分。

5 算 例

中部開(kāi)小圓孔（半徑為1 m)的矩形板在上、下兩端受均布拉力p＝10 kN/ m2,基于對(duì)稱(chēng)性,采用如圖9（a)所示的1/4計(jì)算模型,左邊和底部用積分形式的罰函數(shù)法施加法向約束。計(jì)算圓孔中部邊界點(diǎn)的應(yīng)力集中系數(shù),理論值是3。

采用如圖1（b)所示的基于矩形獨(dú)立覆蓋的新型流形法,充分利用其覆蓋網(wǎng)格自動(dòng)生成、h型覆蓋加密以及p型升階方便的優(yōu)勢(shì),在誤差指標(biāo)控制下嘗試h-p型混合自適應(yīng),初步實(shí)現(xiàn)了CAE自動(dòng)化分析:結(jié)構(gòu)外形（含曲線(xiàn)邊界)由CAD輸入,人工只需輸入材料參數(shù)和邊界條件,其它所有工作完全交由計(jì)算機(jī)完成（另文介紹)。

最終的流形元網(wǎng)格如圖9（a)所示,圓孔周邊局部在圖9（b)中放大,1/4圓弧被矩形網(wǎng)格分為5段,分別位于3個(gè)獨(dú)立覆蓋及2個(gè)條形內(nèi)。在包含上述曲線(xiàn)邊界的流形元的積分運(yùn)算中,取每段圓弧中點(diǎn)及2個(gè)端點(diǎn)的坐標(biāo)為插值點(diǎn),采用2階多項(xiàng)式逼近圓弧邊界。計(jì)算程序自動(dòng)設(shè)置獨(dú)立覆蓋的多項(xiàng)式階次:大多數(shù)為2階,只在圓孔周邊的少數(shù)幾個(gè)覆蓋用到3至4階多項(xiàng)式。自由度總數(shù)為462。

圖9 圓孔算例的1/4計(jì)算模型Fig.9 1/4 model of a rectangular plate with a small circular hole

圖10 垂直向應(yīng)力與理論值的對(duì)比Fig.10 Comparison between calculated data and analytical data of vertical stress

計(jì)算得到的底部一條線(xiàn)上的垂直向應(yīng)力σy見(jiàn)圖10（標(biāo)注“曲線(xiàn)邊界”),可見(jiàn)與理論解［16］很好地符合。其中圓孔邊界點(diǎn)的σy為30.16 kN/ m2,應(yīng)力集中系數(shù)約為3.02,與理論值非常接近,自由表面應(yīng)力σx和τxy分別為0.09 kN/ m2和0.05 kN/ m2,與理論上的零應(yīng)力也非常接近。

另外按傳統(tǒng)方法,1/4圓弧用10段折線(xiàn)離散,其他條件不變。計(jì)算得到圓孔邊界點(diǎn)的σy為29.66 kN/ m2,應(yīng)力集中系數(shù)約為2.97。底部一條線(xiàn)上的垂直向應(yīng)力σy如圖10所示（標(biāo)注“折線(xiàn)邊界”),雖然也與理論解較吻合,但不如曲線(xiàn)邊界的計(jì)算精度高,不難看出幾何誤差帶來(lái)的影響。

6 結(jié) 語(yǔ)

本文提出基于CAD幾何的數(shù)值流形方法,與等幾何分析方法相比具有以下特點(diǎn):采用工業(yè)界通用的以邊界數(shù)據(jù)表示的CAD幾何模型；在CAE建模和網(wǎng)格細(xì)化中同樣具有幾何模型的保形性,避免了幾何誤差；針對(duì)流形法通常使用的多項(xiàng)式近似函數(shù)（其計(jì)算過(guò)程比等幾何分析中的NURBS基函數(shù)要簡(jiǎn)單),采用單純形精確積分法實(shí)現(xiàn)任意形狀流形元的解析積分,沒(méi)有等幾何分析方法的“等參”方式所帶來(lái)的諸多問(wèn)題。

在等幾何分析方法中,CAD實(shí)體模型與CAE分析模型相同?？紤]到幾何描述的復(fù)雜性并不等同于物理場(chǎng)的復(fù)雜性,本文方法根據(jù)真實(shí)物理場(chǎng)分布的復(fù)雜程度來(lái)布置數(shù)學(xué)網(wǎng)格和設(shè)置近似函數(shù)階次,似更為合理。而且只需采用自動(dòng)且快速的切割操作,也能實(shí)現(xiàn)CAD模型無(wú)需修改而直接在CAE中建模,甚至完全自動(dòng)化的前處理。配合使用基于獨(dú)立覆蓋的新型流形法,網(wǎng)格加密很方便,對(duì)復(fù)雜求解域的適應(yīng)性更強(qiáng)。只是本文方法的近似函數(shù)不具備等幾何分析方法的高階連續(xù)性。

因此,本文方法為CAD與CAE融合提供了全新的思路,最終能夠?qū)崿F(xiàn)CAD進(jìn)入CAE后的自動(dòng)化分析（另文介紹),即CAD和CAE的完全融合。當(dāng)然,本文僅對(duì)二維問(wèn)題開(kāi)展了初步研究,提出的切割算法、“近似”單純形積分算法只為說(shuō)明方法的基本學(xué)術(shù)思想,算法的適用性及效率等問(wèn)題還有待進(jìn)一步研究。在算例中雖然采用了基于獨(dú)立覆蓋的新型流形法,但本文方法也適用于采用有限元數(shù)學(xué)網(wǎng)格的傳統(tǒng)流形法。將來(lái)擴(kuò)展到三維研究,涉及曲面幾何邊界與平面網(wǎng)格邊界的切割和曲面“近似”單純形的解析積分,本文的思路同樣適用。

致謝：部分重疊覆蓋的思想來(lái)自于石根華博士，筆者的研究一直受到石根華博士的指導(dǎo)，文中有關(guān)常規(guī)單純形積分的部分內(nèi)容采用了長(zhǎng)江科學(xué)院林紹忠教授的研究成果，在此表示由衷的感謝。

參考文獻(xiàn)：

［1］ 百度百科.CAE［EB/ OL］.（2015-06-11)［2015-07-10］. http:/ / baike.baidu.com/ view/112169.htm.

［2］ HUGHES T J R,COTTRELL J A,BAZILEVS Y. Isogeometric Analysis: CAD, Finite Elements, NURBS, Exact Geometry and Mesh Refinement［J］. Computer Methods in Applied Mechanics and Engineering,2005,194:39-41.

［3］ 吳紫俊,黃正東,左兵權(quán),等.等幾何分析研究概述［J］.機(jī)械工程學(xué)報(bào),2015,51（5):114-129.

［4］ 夏 陽(yáng).假設(shè)位移協(xié)調(diào)有限元及其在精確幾何分析中的應(yīng)用［D］.大連:大連理工大學(xué),2013.

［5］ SHI G H. Manifold Method of Material Analysis［C］/ / U. S. Army Research Office. Transactions of the Ninth Army Conference on Applied Mathematics and Computing, Minneapolis, Minnesota, U S A, June 18-21,1991:57-76.

［6］ 石根華.數(shù)值流形方法與非連續(xù)變形分析［M］.裴覺(jué)民譯.北京:清華大學(xué)出版社,1997.

［7］ 張湘?zhèn)?章?tīng)?zhēng)榮,呂文閣,等.數(shù)值流形方法研究及應(yīng)用進(jìn)展［J］.力學(xué)進(jìn)展, 2010, 40（1): 1-12 .

［8］ 祁勇峰,蘇海東,崔建華.部分重疊覆蓋的數(shù)值流形方法初步研究［J］.長(zhǎng)江科學(xué)院院報(bào), 2013,30（1):65-70 .

［9］ 蘇海東,頡志強(qiáng),龔亞琦,等.基于獨(dú)立覆蓋的流形法的收斂性及覆蓋網(wǎng)格特性［J］.長(zhǎng)江科學(xué)院院報(bào), 2016,33（2):131-136.

［10］蘇海東,祁勇峰.部分重疊覆蓋流形法的覆蓋加密方法［J］.長(zhǎng)江科學(xué)院院報(bào), 2013, 30（7):95-100.

［11］林紹忠,祁勇峰,蘇海東.基于矩陣特殊運(yùn)算的高階流形單元分析［J］.長(zhǎng)江科學(xué)院院報(bào), 2006,23（3):36-39.

［12］林 永,陳 浩.用二分法求解一元實(shí)系數(shù)多項(xiàng)式方程的全部實(shí)根［J］.大學(xué)數(shù)學(xué),2008,24（4):89-91 .

［13］孫家廣,楊長(zhǎng)貴.計(jì)算機(jī)圖形學(xué)［M］.北京:清華大學(xué)出版社,1995.

［14］林紹忠.單純形積分的遞推公式［J］.長(zhǎng)江科學(xué)院院報(bào), 2005,22（3):32-34.

［15］蘇方方.動(dòng)態(tài)規(guī)劃算法計(jì)算組合數(shù)Cmn［J］.現(xiàn)代計(jì)算機(jī),2011,（10):35-37.

［16］徐芝綸.彈性力學(xué)（第2版)［M］.北京:人民教育出版社,1982.

（編輯:劉運(yùn)飛)

Preliminary Study of Numerical Manifold Method Based on CAD Geometry

CHEN Ji-zhan1,2,SU Hai-dong1,2,QI Yong-feng1,2
（1. Research Center on Water Engineering Safety and Disease Control Engineering Technology under Ministry of Water Resources, Yangtze River Scientific Research Institute, Wuhan 430010, China；2. Material and Engineering Structure Department, Yangtze River Scientific Research Institute, Wuhan 430010, China)

Abstract:The iterative process of design-analysis-redesign implies urgent requests of the integration of computer aided design（CAD) and computer aided engineering（CAE). In this paper, we present a numerical manifold method（NMM) based on CAD geometry: we can arrange mathematical meshes and set order of approximation functions according to the complexity degree of physical field distribution in CAE, which is more reasonable than isogeometric analysis（IGA) method. Through introducing automatic and fast cutting operations, we can realize direct modeling from CAD model to CAE model without modifications. Moreover, to solve the problem that curves on geometric boundaries are required to be discretized into line segments in present NMM, we put forward algorithms to cut the curves of geometric boundary with the lines of mesh boundary, hence preserving the shape of the geometric model in the procedures of CAE modeling and mesh refinement. Furthermore, as for the polynomial approximation functions usually used in NMM, we deduce analytical integral formula of“approximate”simplex with a curved edge and use it to obtain precise integral calculations of manifold elements with curved boundaries. Finally, we verify the feasibility of the method through an example of a circular hole in a plate. The research offers new thinking for the integration of CAD and CAE, and lays foundation for the automatic analysis from CAD models to CAE.

Key words:numerical manifold method（NMM)；isogeometric analysis（IGA)；CAD geometry；cutting of curves and lines；simplex integration；CAD/ CAE cooperativity

通訊作者：蘇海東（1968-),男,湖北武漢人,教授級(jí)高級(jí)工程師,博士,從事水工結(jié)構(gòu)數(shù)值分析工作和計(jì)算方法研究,（電話(huà))027-82927167（電子信箱)suhd@ mail.crsri.cn。

作者簡(jiǎn)介：陳積瞻（1990-),男,海南樂(lè)東人,碩士研究生,從事水工結(jié)構(gòu)數(shù)值分析工作和計(jì)算方法研究,（電話(huà))18827620840（電子信箱)xincunyihaoqiao@163.com。

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目（51409012)；中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)項(xiàng)目（CKSF2014054/ CL)

收稿日期：2015-07-16;修回日期：2015-09-11

doi:10.11988/ ckyyb.20150599

中圖分類(lèi)號(hào)：TB115;TV311

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1001-5485(2016)02-0137-07