奇異性彈性場分析的新型超級奇異單元法綜述

王從曼，平學(xué)成，王醒醒，陳夢成

（1.天津科技大學(xué)機(jī)械工程學(xué)院，天津 300222； 2. 天津市輕工與食品工程機(jī)械裝備集成設(shè)計(jì)與在線監(jiān)控重點(diǎn)實(shí)驗(yàn)室，天津 300222； 3. 華東交通大學(xué)土木建筑學(xué)院，江西 南昌 330013）

應(yīng)力奇異性經(jīng)常發(fā)生在幾何形狀不連續(xù)的結(jié)構(gòu)中，典型的例子有裂紋、孔洞和夾雜等。 分析奇異性應(yīng)力場時，只有在簡單構(gòu)型情況下才能獲得解析解，對于復(fù)雜的情況，通常采用數(shù)值方法。 有限元法由于其理論的完備性與強(qiáng)大的通用性成為應(yīng)用最廣泛的數(shù)值方法之一。 常規(guī)有限元法求解奇異性應(yīng)力場時精度普遍不高，即使可以通過加密應(yīng)力奇異點(diǎn)領(lǐng)域網(wǎng)格的方法進(jìn)一步減少誤差，但卻又會導(dǎo)致求解效率大大降低。后來研究人員采用應(yīng)力/位移外插法求解奇異性應(yīng)力場， 但是越接近裂紋尖端，應(yīng)力和位移的計(jì)算值與實(shí)際理論值誤差越大，導(dǎo)致應(yīng)力/位移外插法無法得到高精度的結(jié)果。為了提高精度，Nisitani 等[1-2]和Ping 等[3]采用奇異點(diǎn)應(yīng)力法計(jì)算缺陷尖端奇異性應(yīng)力場，這種方法雖然計(jì)算精度較高， 但在求解時需要先給定一個參考問題的精確解，然后利用待求問題和參考問題的應(yīng)力比值來求解奇異性應(yīng)力場，限制了此方法的使用。 在對具有雙材料界面的界面邊進(jìn)行建模時，Chen 等[4]開發(fā)了一種富集元，可以解釋僅受機(jī)械載荷影響的雙材料界面裂紋尖端處的應(yīng)力奇異性。 與Chen 的工作類似，Gadi 等[5]，Pageau 等[6-7]也采用富集有限元來求解奇異性應(yīng)力場，解釋了不同材料連接處的奇異性應(yīng)力行為，但這種方法仍依賴于單元尺寸。 為了克服常規(guī)有限元的缺點(diǎn)，研究人員提出了一系列專門用來處理缺陷尖端應(yīng)力奇異性的特殊單元，這類單元通常被稱為奇異單元。 通過在缺陷尖端使用這類奇異單元，可以對那些含有幾何不連續(xù)部位的應(yīng)力場進(jìn)行有效的數(shù)值分析。

在分析二維平面奇異性應(yīng)力場時，最早使用的裂紋尖端奇異單元是Benzley 型單元[8]和Barsoum型單元[9]，其中Barsoum 型單元就是如今使用廣泛的四分之一點(diǎn)奇異單元。 Abdelaziz 等[10]提出了一種改進(jìn)的四分之一點(diǎn)奇異單元，用于對裂紋尖端的奇異點(diǎn)進(jìn)行建模。 Liu[11]提出了基于邊緣的平滑有限元，Jiang 等[12]在他們的基礎(chǔ)上開發(fā)了一個7 節(jié)點(diǎn)奇異單元， 用來模擬裂紋尖端周圍的奇異性應(yīng)力行為。這些單元通常是通過簡單地改變標(biāo)準(zhǔn)單元的節(jié)點(diǎn)位置來開發(fā)的，不需要重新配制，方便實(shí)際使用。 其后，Tong 等[13]在雜交元概念的基礎(chǔ)上提出了超級混合裂紋尖端奇異單元，將此超級奇異單元與常規(guī)有限元結(jié)合起來，可以分析平面裂紋尖端奇異性應(yīng)力場。Lin 等[14]開發(fā)了一種混合裂紋尖端奇異單元用于雙材料界面分析, 其中混合單元的假定應(yīng)力和位移場是基于復(fù)勢技術(shù)推導(dǎo)出來的。 Tan 等[15]開發(fā)了一種奇異單元來分析在面內(nèi)熱和耦合載荷下雙材料楔形體頂點(diǎn)處的奇異性應(yīng)力場。 擴(kuò)展Mote[16]和Bradford 等[17]的工作，Madenci 等[18]開發(fā)了一 種混合單元，并將其與常規(guī)有限元耦合，可以用于分析在機(jī)械和均勻熱載荷下由材料和幾何不連續(xù)性引起的奇異性應(yīng)力行為。Barut 等[19]基于機(jī)械和熱載荷下的特征函數(shù)展開方法，利用應(yīng)力和位移場的精確解開發(fā)了一種奇異單元，用于分析多相材料的奇異性應(yīng)力場。Zhang 等[20]開發(fā)了一種帶有彈性夾雜物的n邊形奇異單元，用于對具有隨機(jī)分散夾雜物的異質(zhì)材料進(jìn)行力學(xué)分析。Cai 等[21]開發(fā)了一種用于動態(tài)載荷下雙材料界面裂紋的新奇異單元。 Li 等[22]基于混合應(yīng)力函數(shù)開發(fā)了一種新奇異單元，可以用于各向異性材料中裂紋尖端奇異性應(yīng)力場分析，且此奇異單元形狀和節(jié)點(diǎn)的數(shù)量可靈活調(diào)整。

現(xiàn)有的特殊奇異元大多是針對二維問題，對于三維問題的探究還處于一個不太成熟的階段。 從二維到三維，除了模型自由度的增多，相當(dāng)一部分二維假設(shè)不再適用，需要重新對其定義。 與平面斷裂問題相比， 三維斷裂問題的數(shù)值分析難度更大，需要解決的問題更多， 但與實(shí)際工程應(yīng)用更接近，具有很大的研究價值。 在三維裂紋的研究方面，目前也有了一定的研究成果。

Kuna 等[23]在混合應(yīng)力模型的基礎(chǔ)上，開發(fā)了一種形函數(shù)包含裂紋尖端應(yīng)力的已知奇異行為的特殊裂紋尖端奇異單元， 用來分析三維彈性斷裂問題。李翠華[24]構(gòu)造了一種新的三維奇異單元，可以用于受遠(yuǎn)場載荷的圓環(huán)形和半橢圓形表面裂紋尖端應(yīng)力場分析，擴(kuò)展了Ingraffea 等[25]和Banks-skill[26]的工作。 Grummitt 等[27]考慮了用于三維斷裂問題的拉格朗日奇異單元的細(xì)節(jié)，并證明了它的準(zhǔn)確性和效率。 Tracey[28]提出了三維奇異單元分析裂紋的幾何參數(shù)，隨后Joao 等[29]使用此奇異單元對一系列斷裂力學(xué)問題進(jìn)行建模。 結(jié)果表明，在使用正確設(shè)計(jì)的網(wǎng)格時，分析具有一定的準(zhǔn)確性。 Ariza 等[30]提出了一種具有平面幾何形狀的9 節(jié)點(diǎn)二次奇異單元用于3D 斷裂力學(xué)分析。 陳偉華[31]在三維裂紋的尖端區(qū)域采用二維平面和反平面裂紋問題的辛解析本征解作為插值函數(shù)，構(gòu)造出一類三維辛解析奇異單元，將此奇異單元與三維常規(guī)有限元相結(jié)合可對三維斷裂問題進(jìn)行數(shù)值分析。Hu 等[32]提出了一種新辛解析奇異單元用于三維裂紋奇異性應(yīng)力場分析。

目前為止， 奇異單元已取得廣泛發(fā)展與應(yīng)用，但是隨著制造和材料技術(shù)的不斷進(jìn)步，奇異單元仍存在需要完善之處，例如數(shù)值結(jié)果仍然依賴于奇異單元的大??；除混合奇異單元外，由于節(jié)點(diǎn)處自由度的不同，奇異單元和常規(guī)單元之間的單元相容性無法得到滿足， 這些奇異單元很難收斂到精確解；多數(shù)奇異單元依賴于解析解，通用性不強(qiáng)。 針對這些不足，建立了一系列用于奇異性彈性場分析的新型超級奇異單元。 新型超級奇異單元是基于奇異性位移場和應(yīng)力場的數(shù)值特征解和H-R 變分原理構(gòu)造的包含缺陷角部的特殊單元，能夠用于分析奇異性彈性場。 與其它奇異單元相比，新型超級奇異單元在保證精度和穩(wěn)定性的基礎(chǔ)上，不依賴于單元尺寸，無需任何過渡單元即可與常規(guī)單元連接。由于其處理過程是基于數(shù)值解建立的，整個問題的實(shí)施過程具有很強(qiáng)的通用性，可以應(yīng)用于任意楔形角和任意材料組合的奇異性彈性場分析。 主要綜述了新型超級奇異單元在二維或三維中的裂紋尖端、夾雜角尖端、V 型缺口角尖端、孔洞角尖端、多相材料界面角尖端的應(yīng)用情況，對新型超級奇異單元的理論推導(dǎo)、單元建立和組裝以及分析過程進(jìn)行了介紹。

1 平面楔形角尖端奇異單元

工程結(jié)構(gòu)中， 如動力系統(tǒng)中轉(zhuǎn)向架焊接接頭，由于剛度失配導(dǎo)致應(yīng)力奇點(diǎn)，是產(chǎn)生裂紋的重要原因之一，奇異場分析對結(jié)構(gòu)組件安全要求有重要作用。工程材料中經(jīng)常存在V 型孔洞和缺口等非橢圓型微缺陷。 在小應(yīng)變彈性理論的背景下，這些微缺陷中V 型角尖端引起了奇異應(yīng)力，繼而引起裂紋萌生。 要評價含此類缺陷的疲勞強(qiáng)度，必須準(zhǔn)確有效地確定V 型角尖端奇異場。為獲得高精度和高效率的奇異彈性場數(shù)值解，基于Yamada 等[33]有限元特征法，建立一種新型有限元模型，通過引入陳夢成和Sze 的非協(xié)調(diào)元[34]，得到楔形角部彈性場數(shù)值解。

1.1 V 型角尖端奇異單元

如圖1 所示為V 型角尖端域。為了得到V 型角尖端奇異彈性場變化函數(shù)， 需要依據(jù)虛功原理，由特征問題方程來確定特征值和特征向量。 位移形函數(shù)采用非協(xié)調(diào)元， 而材料界面間則滿足協(xié)調(diào)條件。通過位移與應(yīng)變關(guān)系式、Hooke 定理、 材料彈性矩陣，將特征問題方程離散為有限元格式，建立標(biāo)準(zhǔn)特征方程，求解特征值和特征向量[35]。如果求解結(jié)果為多特征值，依據(jù)Teocaris[36]理論，通過曲線擬合和形函數(shù)插值得到的奇異彈性場為

圖1 V 型角尖端域及局部坐標(biāo)系（n，t）Fig.1 V-shaped corner tip domain and the definition of local coordinate system（n，t）

為了準(zhǔn)確獲得V 型角尖端奇異性應(yīng)力場，開發(fā)了超級V 型角尖端單元[37-38]。 建立新型有限元方法的思想是將原問題圖2（a）分解為兩個邊值問題，如圖2（b）和圖2（c）所示：①在域Ωc內(nèi)的混合邊界值問題，Γc和Cc為圓周邊界條件，該域內(nèi)應(yīng)用Pian 和Sumihara[39]開發(fā)的4 節(jié)點(diǎn)四邊形混合應(yīng)力單元；②在楔形域Ωw內(nèi)的混合邊界值問題，Γw為圓周邊界條件。

圖2 二維角尖端網(wǎng)格劃分Fig.2 Mesh division in a 2D corner-shaped domain

用H-R 原理建立超級V 型角尖端單元， 在Ωw域上的混合泛函形式簡化如下

式中：帶～的參數(shù)為邊界上的定義；Sw為單元柔度矩陣；t 為邊界應(yīng)力向量；D 為材料彈性矩陣；σw為V型角尖端應(yīng)力向量；uw為V 型角尖端位移向量；

通過構(gòu)造位移場和應(yīng)力場來簡化H-R 泛函，使用散度定理進(jìn)行降維運(yùn)算。 對泛函πw進(jìn)行變分，求解其穩(wěn)態(tài)值，最終得到超級V 型角尖端單元的單元剛度矩陣。

式中：σt，r→0，σnt，r→0分別為V 型角部尖端基于局部坐標(biāo)系原點(diǎn)（n，t）附近的t 向應(yīng)力和剪應(yīng)力;λ1和λ2分別為張開型（Ⅰ模型）和滑開型（Ⅱ模型）對應(yīng)的應(yīng)力奇異階數(shù)。

如圖3 所示，在拉伸載荷下，兩個沿x 方向排列的類鉆石型孔洞。 由于對稱，對右半平面進(jìn)行了建模。超級V 型角尖端單元可用于計(jì)算孔洞的局部奇異性應(yīng)力。 需要指出的是，角尖端單元的節(jié)點(diǎn)總數(shù)可能隨其大小和位置而變化。 V 型角尖端o 周圍網(wǎng)格的細(xì)化如圖4 所示。 然而，在傳統(tǒng)的有限元分析中，需要采用精細(xì)網(wǎng)格來獲得V 型角尖端附近的奇異性彈性場。 在V 型角尖端處，細(xì)化后的單元邊長最小尺寸約為10-5l。 新型超級奇異角單元法在網(wǎng)格劃分方面具有簡捷性[41]。

圖3 拉伸載荷下兩個水平排列的類鉆石型孔洞Fig.3 Two horizontally arranged diamond-like holes under tensile load

圖4 類鉆石型孔洞角尖端o 附近的傳統(tǒng)有限元網(wǎng)格Fig.4 Traditional finite element mesh near the corner tip o of the diamond-like hole

各向異性材料的奇異性彈性場研究與各向同性材料大體相同， 不同點(diǎn)在于彈性矩陣不再為定值，會隨著高斯點(diǎn)的變化而變化[42]。

各向異性雙材料楔形問題的廣義應(yīng)力強(qiáng)度因子表達(dá)式為[43]

1.2 雙材料楔形角尖端奇異單元

一般雙材料楔形體如圖5 所示，由兩個楔形體組成。 兩個楔形沿公共邊結(jié)合在一起，該邊構(gòu)成一個界面。 雙材料楔形體尖端的奇異性應(yīng)力場取決于兩方面問題：①是雙材料楔形體的結(jié)構(gòu)；②是雙材料之間的屬性匹配關(guān)系。 在求解奇異性應(yīng)力場過程中，雙材料楔形體尖端的數(shù)值特征解與單相材料楔形體是不同的。

圖5 一般雙材料楔形體Fig.5 Definition of a general bimaterial wedge

值得注意的是， 楔形角尖端單元滿足LBB 準(zhǔn)則：應(yīng)力參數(shù)個數(shù)，也就是待定系數(shù)βn的的個數(shù)應(yīng)大于或等于所有節(jié)點(diǎn)自由度總數(shù)減去混合單元的剛體模式[44]，二維問題的剛體模式為3。 圖6 為建立的9 節(jié)點(diǎn)雙材料超級角尖端單元。 對圖7 所示的單軸拉伸界面上帶有楔形缺口的雙材料面板進(jìn)行處理， 楔形體尖端域采用一個雙材料超級角尖端單元，周圍采用常規(guī)單元相結(jié)合。 分析結(jié)果精確，開發(fā)的雙材料超級角尖端單元在處理楔形問題方面具有通用性和適用性[34]。

圖6 9 節(jié)點(diǎn)雙材料超級角尖端單元Fig.6 Definition of nine-node bimaterial super wedge-tip element

圖7 雙材料楔形角問題的幾何尺寸和網(wǎng)格劃分Fig.7 Geometry and mesh division for the bimaterial wedge problems

1.3 壓電復(fù)合材料奇異單元

壓電材料是具有壓電特性的智能材料，例如一般電子監(jiān)控設(shè)備中的傳感器材料中就包含壓電材料。 自Parton[45]以來，壓電斷裂力學(xué)的理論有一定發(fā)展，然而在這些研究中有一個缺陷：假定裂紋內(nèi)部介質(zhì)的介電常數(shù)為零，而實(shí)際情況是電場可以自由導(dǎo)通。 Scherzer 和Kuna[46]分析了嵌入式智能復(fù)合材料的裂紋結(jié)構(gòu)。 在他們的研究中，裂紋面間的電邊界條件仍然局限于非導(dǎo)通條件。 一般情況下，對于裂紋力學(xué)與電學(xué)不同狀態(tài)，主要存在不導(dǎo)通、導(dǎo)通和有限導(dǎo)通三種邊界條件。 不導(dǎo)通或?qū)ǖ倪吔鐥l件可以看作是部分導(dǎo)通的極限情況。 對于楔形裂紋，發(fā)現(xiàn)如果開口楔形角大于3.5°，則空氣的介電常數(shù)可以忽略，也就是說，不導(dǎo)通條件的影響與導(dǎo)通邊界條件的影響相同。

根據(jù)線彈性理論的應(yīng)變位移關(guān)系和電場電勢關(guān)系，可得到線性壓電材料的本構(gòu)方程[47]

式中：σ（x，y）為應(yīng)力張量；ε（x，y）為應(yīng)變張量；D（x，y）為電位移張量；C 為壓電材料彈性常數(shù)矩陣；E（x，y）為電場張量；e 為材料壓電常數(shù)矩陣。

通過應(yīng)力平衡方程、電荷守恒方程、應(yīng)力自然及本質(zhì)邊界條件、電荷自然及本質(zhì)邊界條件求解平面壓電材料的廣義H-R 變分泛函為

針對材料裂紋尖端奇異性，應(yīng)用超級角尖端混合元模型[34]，通過分部積分和散度定理處理H-R變分泛函， 采用拉格朗日線性插值方法建立相鄰節(jié)點(diǎn)間的形函數(shù)， 使單元間位移和電勢自動滿足相容性。 通過泛函∏的駐值條件，得到壓電材料超級角尖端單元的剛度矩陣Kw。 建立包含部分角尖端的整個壓電域的有限元代數(shù)方程， 確定整體位移U（x，y）和電勢Φ（x，y），將超級角尖端單元進(jìn)一步推廣到壓電斷裂力學(xué)中[48-50]。

圖8 為一個含9 節(jié)點(diǎn)的超級角尖端單元，用以分析楔形體頂點(diǎn)處的應(yīng)力場強(qiáng)度和電彈場強(qiáng)度。 如圖9 所示為遠(yuǎn)場拉伸載荷σy∞和電位移載荷Dy∞下含中心裂紋試件，該裂紋為部分導(dǎo)通。 由于對稱，只考慮右半面板分析。 單元劃分采用傳統(tǒng)單元和超級角尖端單元組合形式。 沿正向x 軸（即x＞0）的節(jié)點(diǎn)電勢設(shè)為0。 通過奇異電場的數(shù)值解得到能量釋放率， 結(jié)果精確， 并避免了在裂紋尖端劃分高密度網(wǎng)格，提高了計(jì)算效率；通過裂紋尖端的應(yīng)力強(qiáng)度因子和電位移強(qiáng)度因子結(jié)果發(fā)現(xiàn)， 裂紋面電邊界條件對KⅠ和KⅡ影響不大，而對電位移強(qiáng)度因子KⅣ有較大影響，導(dǎo)通裂紋時，對應(yīng)于Ⅳ型斷裂模式的特征值λⅣ不再呈現(xiàn)-0.5 奇異性，表明電場奇異性很小[51]。

圖8 9 節(jié)點(diǎn)超級角尖端單元Fig.8 A nine-node super corner-tip element

圖9 壓電平板中心裂紋問題的幾何模型和網(wǎng)格劃分Fig.9 Geometry and mesh division for the central crack problem in a piezoelectric pane

1.4 平面夾雜角尖端奇異單元

復(fù)合材料中，夾雜物的形狀可以是圓形、橢圓形或其它不規(guī)則形。 在制造或使用過程中，夾雜界面處往往會出現(xiàn)應(yīng)力集中， 裂紋可以在夾雜-基體界面或附近形核，影響耐久性。

為求解彈性材料中夾雜角的奇異彈性場，建立如圖10 所示的含部分夾雜角的超級n 邊形單元。將原始問題（圖10（a））分解為：在夾雜域Ω2中的混合邊界條件為C2，ΓA和ΓB（圖10（b））；在基體域Ω1中的混合邊界條件為C1，ΓA和ΓB（圖10（c））。同樣，根據(jù)H-R 原理， 利用散度定理定義以下兩個單獨(dú)的混合泛函[52]

圖10 單元分解Fig.10 Element decomposition

為了從兩個分解問題的解中恢復(fù)原問題的解，需要在邊界ΓA和ΓB施加牽引交換條件和位移協(xié)調(diào)條件

利用彈性材料夾角附近漸進(jìn)位移和應(yīng)力的一般表達(dá)式進(jìn)而得到超級夾雜角尖端單元的剛度矩陣。 超級夾雜角尖端單元用于近場區(qū)域的建模，并與遠(yuǎn)場區(qū)域的傳統(tǒng)4 節(jié)點(diǎn)單元結(jié)合，最終得到夾雜角附近的奇異性應(yīng)力場[53-54]。

圖11 為在平面應(yīng)力條件下， 受拉伸和剪切荷載下包含矩形夾雜物的各向同性無限平板。 圖12為矩形夾雜角附近網(wǎng)格劃分的構(gòu)型。 由于幾何形狀和載荷的對稱性，有限元網(wǎng)格劃分只需要幾何形狀的1/4。 在數(shù)值計(jì)算中，可以結(jié)合使用一個超級夾雜角尖端單元和285 個常規(guī)4 節(jié)點(diǎn)四邊形單元。 相對于傳統(tǒng)有限元，超級夾雜角尖端單元法能夠節(jié)省夾雜尖端的網(wǎng)格數(shù)量，方便建立模型，提高計(jì)算效率。研究得到的無量綱應(yīng)力強(qiáng)度因子數(shù)值結(jié)果具有精確性[55]；對于單矩形夾雜物問題，當(dāng)l/b＞10 時，可以忽略l/b 對廣義強(qiáng)度因子的影響， 將矩形夾雜物視為纖維。

圖11 夾雜物形態(tài)Fig.11 Configuration of inclusions

圖12 矩形夾雜角周圍的網(wǎng)格劃分Fig.12 Configuration of mesh division around a rectangular corner

1.5 熱-機(jī)載荷下的奇異單元

熱-機(jī)載荷下由于材料彈性特性和熱膨脹系數(shù)不匹配，在制造或使用期間材料界面處通常會出現(xiàn)高度集中的應(yīng)力[56]。

1.5.1 熱載荷下的等效奇異性應(yīng)力場

對于由于溫度變化而引起的奇異熱應(yīng)力場的變化問題，首先需要從物理角度出發(fā)，分析機(jī)械載荷與熱載荷的關(guān)系，然后計(jì)算機(jī)械載荷作用下奇異性應(yīng)力場，以及等效熱載荷下的奇異熱應(yīng)力場。

如圖13 所示： ①機(jī)械載荷σ0和溫度變化△T下的無限基體平板與夾雜平板（圖13（a））；②切割某形狀夾雜， 嵌入去除同形狀基體平板的孔洞中，保持滿足力學(xué)平衡條件和位移協(xié)調(diào)條件 （圖13（b））；③反向施加載荷σ0（圖13（c））。

圖13 ΔT 下的熱應(yīng)力場變化Fig.13 Change of thermal stress field under ΔT

經(jīng)分析可知，在邊界自由、界面全結(jié)合、溫度均勻變化的條件下，夾雜角尖端奇異性應(yīng)力場為

式中：（r，θ） 為以夾雜角頂點(diǎn)o 為原點(diǎn)的極坐標(biāo)；σ0為機(jī)械載荷； KⅠ，λ1｜σ0和KⅡ，λ2｜σ0為σ0產(chǎn)生的Ⅰ型和Ⅱ型奇異性應(yīng)力場強(qiáng)度因子；fⅠrr（θ），fⅡrr（θ）等分別為對應(yīng)的Ⅰ型和Ⅱ型應(yīng)力角分布函數(shù)；λ1，λ2分別為對應(yīng)的Ⅰ型和Ⅱ型應(yīng)力奇異指數(shù)。

由溫度變化帶來的奇異熱應(yīng)力場變化值等于遠(yuǎn)場機(jī)械載荷σ0導(dǎo)致的奇性應(yīng)力和非奇性項(xiàng)相加，σrθ沒有附加常數(shù)項(xiàng)。 溫度變化產(chǎn)生熱應(yīng)力場的等效遠(yuǎn)場機(jī)械載荷為[57]

式中：μ1，α1，k1為基體材料參數(shù)， 包括剪切模量，熱膨脹系數(shù)， 楊氏彈性模量；μ2，α2，k2為夾雜材料參數(shù)；ΔT 為變化溫度。

對于夾雜尖端部奇異性熱應(yīng)力場分析，可以直接利用超級夾雜角尖端奇異單元法分析的機(jī)械載荷作用下奇異性應(yīng)力場數(shù)值結(jié)果，計(jì)算奇異性熱應(yīng)力場數(shù)值解[57]。

1.5.2 熱-機(jī)載荷下的超級楔形角尖端單元

在線性彈性理論中，溫度變化下的本構(gòu)方程與純機(jī)械載荷下的本構(gòu)方程不同。 在熱機(jī)械載荷下，雙材料楔形體中的應(yīng)力和位移場可以寫為

式中：qs和qc分別為超級楔角尖端單元和常規(guī)單元的公共節(jié)點(diǎn)位移向量；Ks和Kc分別為超級楔角尖端單元和常規(guī)單元的單元剛度矩陣；tFc為熱載荷引起的節(jié)點(diǎn)反作用力向量；mFc為常規(guī)單元的邊界牽引力向量；H，f為自定義矩陣[59]；∏0為熱負(fù)荷下的與位移和應(yīng)力相關(guān)的總電勢。

為了組裝超級楔形角尖端單元和常規(guī)單元，令邊界線處的節(jié)點(diǎn)對應(yīng)的位移具有連續(xù)性。 對泛函∏t進(jìn)行變分導(dǎo)出下式， 進(jìn)而求解節(jié)點(diǎn)位移qc和熱-機(jī)載荷下的奇異性應(yīng)力場。

圖14 超級楔形角尖端單元Fig.14 The super wedge-tip element

圖15 具有90°～150°粘合楔形角結(jié)構(gòu)Fig.15 The structure with a 90°～150° bonded wedge

1.5.3 熱-機(jī)載荷下的超級夾雜角尖端單元

熱-機(jī)械載荷作用下夾雜角尖端附近的奇異性應(yīng)力場可以通過使用超級夾雜角尖端單元和混合應(yīng)力單元來獲得。 如圖16 為夾雜附近的局部坐標(biāo)系，x 軸位于夾雜角分線上，Ω1為基體域，Ω2為夾雜域。 熱-機(jī)載荷下夾雜角尖端附近奇異性應(yīng)力場的特征解與應(yīng)力平衡條件、牽引交換條件和相容條件有關(guān)。

圖16 夾雜附近的局部坐標(biāo)系Fig.16 Local coordinate system near the inclusion

應(yīng)力平衡條件為

牽引交換條件為

從某種意義上說，上式中第2 項(xiàng)和第3 項(xiàng)是由位移而不是規(guī)定的牽引邊界條件導(dǎo)出，它與傳統(tǒng)的虛功原理是不同的。 如圖17 為一個n 邊型超級夾雜角尖端單元。 根據(jù)熱-機(jī)載荷下夾雜角特征解建立插值函數(shù)。 建立的超級夾雜角尖端單元，能夠保證超級單元和常規(guī)單元間公共節(jié)點(diǎn)的位移連續(xù)性。圖18 為熱載荷下含矩形夾雜物的矩形單胞， 均勻變化的溫度為ΔT=100 ℃。 在夾雜角尖端采用一個超級夾雜角尖端單元，周圍采用傳統(tǒng)單元。 如圖所示利用2 種超級單元來分析奇異性應(yīng)力場：8 節(jié)點(diǎn)單元和16 節(jié)點(diǎn)單元。 2 種超級單元尺寸hx=hy=0.05l1；hx=hy=0.1l1。分析的奇異性應(yīng)力場對超級夾雜角尖端單元的節(jié)點(diǎn)數(shù)量和單元尺寸不敏感，并且在粗網(wǎng)格的情況下能夠提供令人滿意的結(jié)果[60-62]。

圖17 n 邊型超級夾雜角尖端單元Fig.17 Super n-sided polygonal inclusion corner-tip element

圖18 熱載荷下含矩形夾雜物的矩形單胞Fig.18 Rectangular unit cell containing a rectangular inclusion for pure thermal loads

2 三維楔形角尖端奇異單元

三維角部尖端問題所帶來的數(shù)學(xué)困難遠(yuǎn)遠(yuǎn)大于二維角部尖端問題。 與二維問題不同的是，三維角部尖端附近的漸近場包含了3 種斷裂模式。 目前三維問題研究的目標(biāo)之一是解決角部前沿奇異性應(yīng)力場。 為此，基于H-R 變分原理，提出一種新的包含部分角尖端的超級奇異三維單元。 由于解析本征解通常應(yīng)用于簡單構(gòu)型， 數(shù)學(xué)推導(dǎo)過程較復(fù)雜，可以采用廣義平面應(yīng)變問題的一維有限元公式求解三維位移場和應(yīng)力場強(qiáng)度。

2.1 包含直線角線的三維角構(gòu)型奇異單元

如圖19（a）所示為包含恒定二面角直邊的三維均勻線彈性實(shí)體。 在這種情況下，一個沿直線角前沿的特定直線坐標(biāo)z 軸可以代替一般曲線角問題的曲線坐標(biāo)。 依據(jù)Sze 和Wang[63]，任何法向平面的漸近位移和應(yīng)力場可以由二維廣義平面應(yīng)變解推導(dǎo)而得。 根據(jù)控制方程的弱形式，由有限元特征分析方法來確定特征值[64]。

式中：u（r，θ）包含的位移分量為ur，uθ，uz；σ（r，θ）包含的應(yīng)力分量為σr，σθ，σz，σrθ，σrz，σθz；?e{ }表示取復(fù)數(shù)的實(shí)部；λn為特征值；N 為復(fù)數(shù)特征值個數(shù)，M為實(shí)數(shù)特征值個數(shù)；βn為未知應(yīng)力強(qiáng)度系數(shù)， 其與遠(yuǎn)場邊界和載荷條件有關(guān)；un（θ），σˉn（θ）分別為位移和應(yīng)力角分布函數(shù)；Un是位移級數(shù)向量表達(dá)式的第n 項(xiàng)。

如圖19（b）所示，為確定彈性材料中未知系數(shù)βn和計(jì)算圍繞角尖端奇異彈性場，在笛卡爾坐標(biāo)系中采用包含部分角尖端域的超級n 邊單元。 根據(jù)Zhang 和Katsube[20]定義角尖端單元域問題的混合泛函形式

圖19 包含直線角線的三維角構(gòu)型奇異單元Fig.19 Three-dimensional angular configuration singular element containing straight corners

為求解泛函的駐值，構(gòu)造應(yīng)力和位移場，通過散度定理進(jìn)行泛函簡化[64]。在三維問題中，βn隨z 值的變化而變化。 引入新的自然坐標(biāo)建立三維角尖端等參單元。 通過極坐標(biāo)和笛卡爾坐標(biāo)的位移和應(yīng)力轉(zhuǎn)換，利用δπ=0 得到公式

圖20 所示為沿厚度方向有階躍變化的復(fù)合材料層壓板。 圖21 展示了兩種具有不同三維角構(gòu)型奇異單元的網(wǎng)格劃分形式。 圖21（a）和圖21（b）分別包含一個14 節(jié)點(diǎn)和一個26 節(jié)點(diǎn)的奇異單元。 兩種模式的單元尺寸不同， 但分析結(jié)果令人滿意，而相同區(qū)域內(nèi)，Pageau 和Biggers[7]則需要一個包含182個集中單元才能獲得滿意的結(jié)果。 三維角構(gòu)型奇異單元比集中單元的數(shù)量要少，具有高效性。 應(yīng)用建立的三維角構(gòu)型奇異單元，還可以分析三維貫穿中心裂紋問題、三維單邊裂紋問題和三維對接接頭問題的奇異性應(yīng)力場[64]。

圖20 厚度有階梯變化的復(fù)合材料層壓板Fig.20 Composite laminate with step change in thickness

圖21 兩種不同的超級三維角尖端單元網(wǎng)格Fig.21 Meshes with two difference kinds of 3D corner-tip elements

2.2 三維曲線裂紋前沿奇異單元

在疲勞裂紋擴(kuò)展過程中，裂紋前沿通常為曲線形狀（如圓形裂紋或橢圓裂紋），必須考慮裂紋問題的三維特征[65]。 Ayhan 等[66]采用三維富集有限元計(jì)算了曲線裂紋的應(yīng)力強(qiáng)度因子， 模擬了疲勞裂紋擴(kuò)展，而利用富集有限元必須面臨3 個困難：①數(shù)值結(jié)果仍然取決于特殊單元的大??； ②由于節(jié)點(diǎn)處的自由度不同， 不能滿足特殊單元與常規(guī)單元之間的單元間兼容性， 通常需要過渡單元將特殊元素與常規(guī)單元結(jié)合起來； ③級數(shù)解的階項(xiàng)不包含剛體運(yùn)動模式， 在常規(guī)有限元方法中必須與位移插值一起使用。

基于奇異性位移場和應(yīng)力場的數(shù)值特征解和H-R 變分原理， 建立包含部分三維曲線裂紋前沿的奇異單元，可以求解包括應(yīng)力奇點(diǎn)階數(shù)、位移角和應(yīng)力角變化在內(nèi)的數(shù)值系列特征解， 用于三維實(shí)體中的硬幣形裂紋、 圓柱體中的周向裂紋和無限體中的橢圓裂紋等問題分析[67]。在曲線裂紋前沿附加的笛卡爾（x，y，z）和圓柱坐標(biāo)（r，θ，z）中建立局部曲線坐標(biāo)系（ρ，?，θ），可以方便地描述裂紋前沿附近的局部變形行為。 對于裂紋尖端的彈性漸近場可表示為

一般來說，沿任意曲線裂紋前沿不同法向面上的奇異性位移場和應(yīng)力場的強(qiáng)度不是恒定的。 為了解決三維效應(yīng)的影響， 建立如圖22 所示的三維曲線裂紋前沿奇異單元。 單元域由兩個端面Γe和一定數(shù)量的周向面Γc包圍。 在三維曲線裂紋前沿奇異單元中，需要引入另一個自然坐標(biāo)η 來表示裂紋前沿法平面ρ-?。此外，還需定義用于離散化裂紋前沿位移和應(yīng)力角變化的自然坐標(biāo)ξ。 因而系數(shù)向量β 可以被認(rèn)為是η 的線性函數(shù)。 在局部三維坐標(biāo)系（ρ，?，θ）定義位移和應(yīng)力函數(shù)

圖22 包含部分曲線裂紋前沿的超級單元Fig.22 A super element containing a part of the curved crack front

為了獲得三維曲線裂紋前沿奇異單元的剛度矩陣，需要利用H-R 變分的駐值和散度定理，通過位移邊界條件求解[67]。圖23 為承受扭轉(zhuǎn)載荷下的含硬幣形裂紋。 為了清楚地說明裂紋前沿附近的網(wǎng)格劃分，圖24 中繪制了模型的四分之一，在這種情況下， 沿裂紋前沿分布一定數(shù)量的超級18 節(jié)點(diǎn)三維曲線裂紋前沿奇異單元，其周圍是傳統(tǒng)的三維8 節(jié)點(diǎn)單元。 當(dāng)單元尺寸小于或等于0.02a×0.02a 且單元數(shù)大于30 時，數(shù)值結(jié)果可以趨于收斂。 分析周向裂紋問題的應(yīng)力強(qiáng)度因子表明， 當(dāng)使用超過30 個單元時， 即使應(yīng)用尺寸為0.1a×0.1a 的奇異單元也足以得到收斂解[67]。

圖23 扭轉(zhuǎn)作用下的幣狀裂紋Fig.23 A penny-shaped crack problem subjected to torsion loading

圖24 幣狀裂紋的網(wǎng)格劃分Fig.24 Mesh division for the penny-shaped crack problem

2.3 三維V 型缺口前沿奇異單元

大多數(shù)情況下，缺陷端部通常不是圓形或橢圓形狀。在孔邊可能出現(xiàn)V 型角引起的局部奇異應(yīng)力導(dǎo)致嚴(yán)重的損傷或疲勞斷裂。三維V 型角前沿奇異場的研究與三維裂紋奇異場研究原理相似， 但與Griffith 裂紋問題不同的是，V 型角的應(yīng)力奇點(diǎn)數(shù)值特征解一般不等于-0.5， 奇異性應(yīng)力場也隨夾角的變化而變化。如圖25（a）為三維類鉆石型缺陷。為了求解三維V 型角附近的三維奇異性應(yīng)力場和相應(yīng)的高階項(xiàng)級數(shù)解，采用廣義平面應(yīng)變問題的特殊有限元方法[63]，在局部坐標(biāo)系中確定數(shù)值特征解。通過該方法可以得到應(yīng)力奇異性的階數(shù)λn，以及位移和應(yīng)變的角變函數(shù)un（?，θ）和εn（?，θ）。

局部坐標(biāo)系中V 型角前沿的位移場u（ρ，?，θ）和應(yīng)力場σ（ρ，?，θ）的一般表達(dá)式為

式中：βn為應(yīng)力強(qiáng)度參數(shù)；D 為線彈性材料的彈性常數(shù)矩陣；uˉn（?，θ），εn（?，θ）分別為位移和應(yīng)變角分布函數(shù)；un（ρ，?，θ），εn（ρ，?，θ）是位移級數(shù)向量表達(dá)式和應(yīng)變級數(shù)向量表達(dá)式的第n 項(xiàng)。

在任意外載荷作用下，V 型角前沿的應(yīng)力強(qiáng)度不是恒定的。如圖25（b）所示可以建立三維V 型缺口前沿奇異單元來解決目前的三維V 型缺口問題。 利用局部坐標(biāo)系下的位移和應(yīng)力函數(shù)， 得到H-R 變分泛函和三維V 型缺口前沿奇異單元的剛度矩陣[68]。

圖25 三維V 型缺口問題Fig.25 Three-dimensional V-notch problem

為了將傳統(tǒng)單元與三維V 型缺口前沿奇異單元組合起來，在公共節(jié)點(diǎn)上使傳統(tǒng)單元與奇異單元的剛度矩陣進(jìn)行組合。

式中： 下標(biāo)i 和j 分別表示公共節(jié)點(diǎn)域和非公共節(jié)點(diǎn)域；ui（x，y，z）和uj（x，y，z）分別表示公共節(jié)點(diǎn)和非公共節(jié)點(diǎn)的位移；Ff是等效節(jié)點(diǎn)力向量；Ks為奇異單 元 的 剛 度矩陣；Kii，Kij，Kji，Kjj分 別 是傳統(tǒng) 單 元 剛度矩陣的分量。

采用新型有限元法計(jì)算圓桿中V 型缺口廣義應(yīng)力強(qiáng)度因子，如圖26 所示，沿四分之一周角前沿分布40 個超級單元，周圍采用傳統(tǒng)單元。 在向量β維數(shù)小于24 時，數(shù)值結(jié)果的計(jì)算精度令人滿意，這意味著三維V 型缺口前沿奇異單元的收斂結(jié)果很容易得到。 應(yīng)用此方法可以進(jìn)一步分析嵌入類金剛石缺陷、接近自由表面缺陷、表面缺陷、缺陷干涉等模型的應(yīng)力強(qiáng)度問題[68]。

圖26 含周向V 型缺口的圓形桿件網(wǎng)格劃分Fig.26 Mesh division for a round bar with a circumferential a V-shaped notch

2.4 三維曲線型界面角前沿奇異單元

孔洞周圍的應(yīng)力集中可能導(dǎo)致雙材料界面的退化和損傷，這是雙材料界面結(jié)合強(qiáng)度減弱的重要原因，有必要建立相應(yīng)的方法求解孔洞與界面交線附近的力學(xué)行為，以確定界面角尖端處的應(yīng)力奇異強(qiáng)度。 雙材料界面與三維軸對稱孔洞的相互作用將形成一個周向的界面角，界面角可以是尖角或不同角度的鈍角。 三維孔洞幾何形式的復(fù)雜性增加了求解漸近應(yīng)力場和邊值問題的難度。 如圖27，以雙材料圓周界面線上點(diǎn)O 為中心建立局部曲線坐標(biāo)系（ρ，?，θ）， 孔洞與ρ-? 平面交點(diǎn)的曲率半徑記為rI?；诰植壳€坐標(biāo)系建立三維曲線型界面角前沿奇異單元，單元包含部分圓周界面角線，由兩種不同體積的彈性材料V1和V2組成。 為了建立奇異單元的有限元方程，用界面邊緣奇異性應(yīng)力場的漸近解來離散位移和應(yīng)力，引入一個自然坐標(biāo)η 來描述角坐標(biāo)θ 的位置。

圖27 包含部分周向界面交線的奇異界面單元Fig.27 A singular interface edge element containing a part of the circumferential interface intersection line

三維曲線型界面角前沿奇異單元域的泛函形式為

式中：N 為4 節(jié)點(diǎn)四邊形單元的形狀函數(shù)矩陣；n 為邊界Sm外法線向量矩陣；βc為應(yīng)力強(qiáng)度參數(shù)；Ecm和Ucm分別為自定義應(yīng)力和位移矩陣[69]；dc為含4 節(jié)點(diǎn)四邊形單元所有節(jié)點(diǎn)分量的位移向量。

應(yīng)用泛函駐值條件δΠH-R=0，導(dǎo)出奇異單元的剛度矩陣

式中：Γc為奇異界面單元的周向界面，如圖27 所示。

對于圖28 所示的界面孔洞問題， 幾何模型可以看作是包含有限界面孔洞缺陷的無限固體問題。為了避免當(dāng)前雙材料界面與三維軸對稱孔洞交線附近的網(wǎng)格細(xì)化問題，在三維曲線型界面角前沿采用奇異單元，周圍使用傳統(tǒng)三維單元。 圖29 顯示了新型有限元法與傳統(tǒng)有限元法在界面角尖端附近的網(wǎng)格劃分區(qū)別。 如圖29（a）所示，傳統(tǒng)有限元在圓周線附近，使用高度細(xì)化的網(wǎng)格。 如圖29（b）所示， 沿雙材料界面和孔洞的周向交線上使用了30 個奇異單元， 新型有限元方法數(shù)值分析結(jié)果具有精確性[69]。

圖28 雙材料界面三維孔洞Fig.28 3D defects at the interface of a bimaterial

圖29 含三維孔洞問題的網(wǎng)絡(luò)劃分Fig.29 Mesh division for problems with 3D voids

3 結(jié)束語

本文主要綜述了新型超級奇異單元在二維或三維中的裂紋尖端、夾雜角尖端、V 型缺口角尖端、孔洞角尖端、多相材料界面角尖端的應(yīng)用情況。 與其它奇異單元相比，新型超級奇異單元具有以下特點(diǎn)：

1） 奇異單元和常規(guī)單元之間節(jié)點(diǎn)處自由度一致，無需任何過渡單元即可與常規(guī)單元連接，單元相容性好。

2） 超級奇異單元在相對較低的計(jì)算成本下可以獲得優(yōu)異的結(jié)果，是分析二維或三維中的裂紋尖端、夾雜角尖端、V 型缺口角尖端、孔洞角尖端、多相材料界面角尖端的奇異性應(yīng)力場的有效方法，適用性廣。

3） 可應(yīng)用于任意楔形角和任意材料組合的奇異場問題，具有良好的通用性和工程應(yīng)用前景。

4） 新型超級奇異單元分析結(jié)果與單元尺寸無關(guān)，具有較高的穩(wěn)定性。

5） 可以用于分析各向同性、各向異性和壓電材料奇異電彈性場，可以解決熱-機(jī)耦合載荷和機(jī)-電耦合載荷作用下的奇異性應(yīng)力場。

6） 在三維應(yīng)用領(lǐng)域，新型超級奇異單元只能對角部尖端線為光滑連續(xù)曲線或直線的奇異性應(yīng)力場進(jìn)行計(jì)算，對幾何不連續(xù)角前沿線奇異性應(yīng)力場的分析目前還有一定的局限性，需要進(jìn)一步研究發(fā)展。