葛仁余，張佳宸，馬國強(qiáng)，劉小雙，牛忠榮

（1.安徽工程大學(xué) 力學(xué)重點(diǎn)實(shí)驗(yàn)室，安徽 蕪湖 241000；2.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，合肥 230009）

引 言

不同材料的連接是機(jī)械工程和材料科學(xué)中經(jīng)常遇到的一種情況.由于組成材料的彈性性質(zhì)不同，連接點(diǎn)是產(chǎn)生結(jié)構(gòu)破壞失效的薄弱部位和損傷源.結(jié)構(gòu)在這些薄弱部位產(chǎn)生強(qiáng)應(yīng)力集中，以至于在彈性力學(xué)意義上應(yīng)力趨于無窮大，這種彈性力學(xué)范圍內(nèi)應(yīng)力趨于無窮大的特性稱為應(yīng)力奇異性.當(dāng)外載荷作用時(shí)，在結(jié)構(gòu)的這些薄弱部位會(huì)產(chǎn)生很高的應(yīng)力集中，會(huì)導(dǎo)致結(jié)構(gòu)在V 形切口或接頭連接點(diǎn)處出現(xiàn)龜裂等多種破壞形式.由于接頭應(yīng)力場的奇性指數(shù)對(duì)結(jié)構(gòu)強(qiáng)度有決定性的影響，因此，獲取奇性指數(shù)對(duì)確保結(jié)構(gòu)服役的可靠性和安全性具有重要意義[1-2].

接頭連接點(diǎn)處應(yīng)力奇點(diǎn)研究方法有若干種，每種方法都有各自的特點(diǎn)，例如，Mellin 變換法是一種計(jì)算應(yīng)力奇性指數(shù)簡便而十分有效的方法[3-5].Williams[6]（1952年）用直接的Airy 應(yīng)力函數(shù)法研究了冪應(yīng)力奇異性.England[7]（1971年），Stern 和Soni[8]（1976年），以及Carpenter 和Byers[9]（1987年）使用了簡單的復(fù)勢(shì)方法來研究冪應(yīng)力奇異性.Paggi 和Carpinteri[10]對(duì)含切口結(jié)構(gòu)和多材料接頭應(yīng)力奇異性進(jìn)行了廣泛的研究，從數(shù)學(xué)上證明了特征函數(shù)展開法、復(fù)變函數(shù)法和Mellin 變換法在求解應(yīng)力奇性指數(shù)方面的等價(jià)性.

為了確定和降低材料微觀結(jié)構(gòu)中的應(yīng)力奇性指數(shù)，人們開始關(guān)注不同材料組成的接頭結(jié)構(gòu)的應(yīng)力奇異性問題[11-13].例如，張金輪和葛仁余等[14]運(yùn)用插值矩陣法分析了平面接頭與界面裂紋的應(yīng)力奇性指數(shù)；Sator等[15]從特征方程出發(fā)，獲得了平面接頭應(yīng)力奇性指數(shù)實(shí)數(shù)解的精確值，用Newton 法獲得了應(yīng)力奇性指數(shù)復(fù)數(shù)解的數(shù)值解；在假設(shè)特征值為復(fù)數(shù)的前提下，Carpinteri 等[16]運(yùn)用特征函數(shù)展開的方法研究了多材料連接點(diǎn)處的應(yīng)力奇異性問題；Cho 等[17]利用復(fù)勢(shì)方法以及復(fù)根的概念，研究了雙材料V 形切口裂紋的冪對(duì)數(shù)應(yīng)力奇異性問題.

1971年，Bellman 和Casti[18]提出了微分求積法（DQM）基本理論，在求解微分方程特征值和邊值問題時(shí)，由于該理論不依賴于變分原理，具有公式簡單、編程方便、高效、精度高等特點(diǎn)，是一種有吸引力的直接求解微分方程的數(shù)值計(jì)算方法，并且能以較少的網(wǎng)格點(diǎn)求得微分方程的高精度數(shù)值解，所以其在振動(dòng)工程領(lǐng)域得到了廣泛應(yīng)用，主要用來求解工程結(jié)構(gòu)的自由振動(dòng)和強(qiáng)迫振動(dòng)問題.本文創(chuàng)新性地將該方法運(yùn)用到工程斷裂力學(xué)的研究領(lǐng)域.根據(jù)彈性力學(xué)基本理論，將異質(zhì)雙材料平面接頭連接點(diǎn)處應(yīng)力奇性指數(shù)的計(jì)算轉(zhuǎn)化為一組常微分方程組（ODEs）的特征值問題，再由DQM 基本理論將其轉(zhuǎn)化為標(biāo)準(zhǔn)型廣義代數(shù)方程組特征值問題，最后，一次性地計(jì)算出異質(zhì)材料平面接頭連接點(diǎn)處應(yīng)力奇性指數(shù)及其相應(yīng)的位移和應(yīng)力特征函數(shù).

1 平面接頭應(yīng)力奇性指數(shù)及其DQM 計(jì)算

圖1為各向同性雙材料平面接頭模型，Γ1和 Γ3為 兩自由楔邊，Γ2為兩種不同材料的交界邊.在接頭連接點(diǎn)O處同時(shí)定義一個(gè)直角坐標(biāo)系xOy和 一個(gè)極坐標(biāo)系rOθ，交界 Γ2與x軸重合，θ逆時(shí)針方向?yàn)檎槙r(shí)針方向?yàn)樨?fù).對(duì)于一個(gè)平面接頭構(gòu)件，基于Williams 漸近展開理論[6]，將接頭連接點(diǎn)處的位移場表達(dá)成關(guān)于徑向r的級(jí)數(shù)漸近展開形式：

圖1 兩相材料平面接頭模型Fig.1 The 2-phase isotropic multi-material junction model

關(guān)于平面應(yīng)力問題，將式（1）代入各向同性材料本構(gòu)方程，獲得接頭連接點(diǎn)處應(yīng)力分量如下：

式（1）和（2）中，λk為應(yīng)力奇性指數(shù)，Ak為位移幅值系數(shù)(k=1,2,3,···,M),M是截取的級(jí)數(shù)項(xiàng)數(shù)；上標(biāo)m=1,2分別表示材料域1和材料域2，(θ)(i=r,θ)為相應(yīng)材料域的徑向和周向位移特征函數(shù)，(θ)(ij=rr,θθ,rθ)為相應(yīng)材料域的應(yīng)力特征函數(shù)，(θ)為(θ)的線性組合，即

式中，E(m)和 ν(m)分 別為相應(yīng)材料域的彈性模量和Poisson 比.考慮平面應(yīng)變問題時(shí)，只需將式（3）中的E(m)用替換，用替換.在極坐標(biāo)系下，無體力的彈性力學(xué)平面問題平衡方程為

為了便于編程計(jì)算，引入歸一化變量ξ，0 ≤ξ≤1.在區(qū)間[θ1,θ2]中，歸一化變量ξ 為

在區(qū)間[θ2,θ3]中，歸一化變量ξ 為

將式（2）、（5）代入平面問題的平衡方程（4）中，得

其中，矩陣列向量L,M和N為(ξ)，(ξ)及其第1 階、第2 階導(dǎo)數(shù)的線性組合（m=1,2），即

圖1平面接頭兩自由楔邊 Γ1和 Γ3上的面力為零，則邊界條件可表示為

圖1平面接頭交界 Γ2上位移連續(xù)、面力相等，則交界條件可由下式表示：

至此，平面接頭應(yīng)力奇性指數(shù)的分析變成了求解在邊界條件（7）和交界條件（8）下，ODEs(6)的特征值問題.這里，運(yùn)用DQM 計(jì)算平面接頭結(jié)構(gòu)的應(yīng)力奇性指數(shù).圖1平面接頭的兩個(gè)楔形圓弧區(qū)間θ1≤θ≤θ2和θ2≤θ≤θ3上的離散單元總數(shù)皆為N，離散節(jié)點(diǎn)總數(shù)皆為N+1，即，其中 ξ0=0，ξN=1.離散節(jié)點(diǎn)分布方式采用等步長均勻網(wǎng)格模式，即

DQM 的基本思想是：將歸一化的位移特征函數(shù)(ξ)，(ξ)（0 ≤ξ≤1）及其r階導(dǎo)數(shù)在其求解域中任意離散節(jié)點(diǎn)上的近似值，表示為所有離散節(jié)點(diǎn)上值的線性加權(quán)和.即

當(dāng)i=j時(shí)，有

式中，R(1)為

將式（10）中表示的DQM 規(guī)則應(yīng)用于ODEs（6），從而得到關(guān)于平面接頭應(yīng)力奇性指數(shù) λk為特征值的一組代數(shù)方程組為

其中

最后，再將邊界條件（7）和交界條件（8）中關(guān)于E(m)和 ν(m)的常系數(shù)表達(dá)式替換式（15）矩陣，和中相應(yīng)的首行和末行.根據(jù)以上公式，筆者用FORTRAN 語言編寫了通用程序DQMEI 用于DQM 計(jì)算平面接頭應(yīng)力奇性指數(shù) λk及其相應(yīng)的位移特征函數(shù)(θ)和 應(yīng)力特征函數(shù)(θ).

2 數(shù)值算例及問題討論

2.1 算例1：平面接頭模型1

圖2為平面接頭模型1，α=180°是 一定值，β角是變化的.考慮平面應(yīng)變情況，并且假定兩種材料的Poisson 比ν(1)=ν(2)=0.2 ；離散單元總數(shù)分別取N=4，6，8，10，12，離散節(jié)點(diǎn)分布形式采用等步長均勻網(wǎng)格模式，運(yùn)用DQM 求解ODEs（6）、邊界條件（7）和交界條件（8），獲得兩相材料接頭連接點(diǎn)處應(yīng)力奇性指數(shù)隨材料的彈性模量比值 η =E(2)/E(1)以及楔角 β 變化的計(jì)算值，如表1～3 所示.從表1～3 可知，在離散單元總數(shù)N=4 和N=6 時(shí)，DQM 計(jì)算值與實(shí)際值誤差較大，隨著N逐漸增加，DQM 計(jì)算值加速收斂.當(dāng)N=12 時(shí)，DQM 計(jì)算值與文獻(xiàn)[14]的計(jì)算結(jié)果至少有5 位有效數(shù)字相同，表明了本文數(shù)值方法計(jì)算平面接頭應(yīng)力奇性指數(shù)的有效性和精確性.

圖2 平面接頭模型1Fig.2 Plane junction model 1

表1 η = 3.0 時(shí)，平面接頭模型1 的第1 階應(yīng)力奇性指數(shù)λ1 計(jì)算值隨離散單元數(shù)N 的變化Table 1 Variation of the 1st-order stress singularity index of plane joint model 1 with number of discrete elements N for η = 3.0

當(dāng) β=90°和 β=15°時(shí)，平面接頭模型1 連接點(diǎn)處應(yīng)力奇性指數(shù)隨材料的彈性模量比值η =E(2)/E(1)的變化曲線如圖3、4 所示.本文DQM 計(jì)算值與文獻(xiàn)[15]的解析解及Newton 迭代法計(jì)算結(jié)果一致，同時(shí)，DQM 計(jì)算結(jié)果還表明：圖3中，β=90°，0.001≤η≤10 000 時(shí)，接頭僅存在一個(gè)實(shí)數(shù)奇性指數(shù)；圖4中，β=15°，0.001<η≤1時(shí)，R e λ=0 平面接頭模型1 連接點(diǎn)處無應(yīng)力奇異性，1 <η≤511.9時(shí)，連接點(diǎn)處存在1 個(gè)實(shí)數(shù)奇性指數(shù)，5 11.9<η≤851.3時(shí)，連接點(diǎn)處存在2 個(gè)實(shí)數(shù)奇性指數(shù)，尤其當(dāng) 8 51.3<η≤10 000時(shí)，由于材料的相互不匹配性，平面接頭模型1 連接點(diǎn)附近區(qū)域應(yīng)力奇性指數(shù)是復(fù)數(shù)，而不是實(shí)數(shù)，產(chǎn)生振蕩應(yīng)力奇異性[23].

表2 η = 4.0 時(shí)，平面接頭模型1 的第1 階應(yīng)力奇性指數(shù)λ1 計(jì)算值隨離散單元數(shù)N 的變化Table 2 Variation of the 1st-order stress singularity index of plane joint model 1 with number of discrete elements N for η = 4.0

表3 η = 5.0 時(shí)，平面接頭模型1 的第1 階應(yīng)力奇性指數(shù)λ1 計(jì)算值隨離散單元數(shù)N 的變化Table 3 Variation of the 1st-order stress singularity index of plane joint model 1 with number of discrete elements N for η = 5.0

圖3 當(dāng)β = 90°時(shí)，平面接頭模型1 的應(yīng)力奇性指數(shù)Fig.3 The singular index of stress in plane joint model 1 for β = 90°

圖4 當(dāng)β = 15°時(shí)，平面接頭模型1 的應(yīng)力奇性指數(shù)Fig.4 The singular index of stress in plane joint model 1 for β = 15°

圖2中，在平面應(yīng)變條件下，當(dāng) β=90°、離散單元總數(shù)N=12 時(shí)，在η =E(2)/E(1)=0.01和 η =E(2)/E(1)=0.1兩種情況下，由本文DQM 獲得連接點(diǎn)處的第1 階應(yīng)力奇性指數(shù)分別為 λ1=?0.209 247 37、 λ1=?0.141 009 68，它們對(duì)應(yīng)的位移特征函數(shù)和應(yīng)力特征函數(shù)分布曲線如圖5、6 所示.從計(jì)算結(jié)果中看出，位移分量和在粘接界面上雖然連續(xù)，位移特征函數(shù)曲線卻出現(xiàn)轉(zhuǎn)折，兩種材料的彈性模量相差越大，轉(zhuǎn)折越明顯，而且，徑向應(yīng)力特征函數(shù)在粘接界面上出現(xiàn)突變，彈性模量相差越大突變?cè)郊ち?，這就是接頭連接點(diǎn)處發(fā)生開裂破壞的 主要原因.但是，另外兩個(gè)應(yīng)力特征函數(shù)和在粘接界面上連續(xù)，無明顯的轉(zhuǎn)折.

圖5 E(2)/E(1) = 0.01 時(shí)，平面接頭模型1 第1 階應(yīng)力奇性指數(shù)λ1 對(duì)應(yīng)的位移和應(yīng)力特征函數(shù)曲線圖Fig.5 Displacement and stress characteristic function curves corresponding to 1st-order stress singular index λ1 of plane joint model 1 for E(2)/E(1) = 0.01

2.2 算例2：平面接頭模型2

圖7為平面接頭模型2，設(shè)交界面兩側(cè)分屬不同的均質(zhì)彈性體材料，材料Poisson 比為ν(1)=ν(2)=0.2；α=180°是一定值，β角是變化的.在平面應(yīng)變條件下，兩個(gè)區(qū)間 [?α,0°]和[0°,β]上離散單元總數(shù)取N=12時(shí)，由DQM 計(jì)算獲得以下結(jié)論：① 如圖8所示，當(dāng) β=180°時(shí)，平面接頭模型2 退化為雙材料界面裂紋模型，當(dāng)0.001≤η≤10 000 時(shí)，界面裂紋尖端應(yīng)力奇性指數(shù)是一對(duì)復(fù)數(shù)，實(shí)部值 Re λ=?0.5，虛部值 Im λ很小.② 如圖9所示，當(dāng) β=135°和 0 .17<η<2.82 時(shí)，平面接頭連接點(diǎn)處應(yīng)力場存在2 個(gè)實(shí)數(shù)奇性指數(shù)；當(dāng) 0 .001≤η≤0.17和2.82≤η≤10 000時(shí)，平面接頭連接點(diǎn)處為2 個(gè)復(fù)數(shù)奇性指數(shù)，產(chǎn)生振蕩應(yīng)力奇異性[23].根據(jù)圖8、9 可知，本文DQM 計(jì)算值與文獻(xiàn)[15]計(jì)算結(jié)果一致，再次證明本文DQM 對(duì)分析一般平面接頭應(yīng)力奇性指數(shù)是一種有效且準(zhǔn)確的手段.

圖6 E(2)/E(1) = 0.1 時(shí)，平面接頭模型1 第1 階應(yīng)力奇性指數(shù)λ1 對(duì)應(yīng)的位移和應(yīng)力特征函數(shù)曲線圖Fig.6 Displacement and stress characteristic function curves corresponding to the 1st-order stress singular index λ1 of plane joint model 1 for E(2)/E(1) = 0.1

圖7 平面接頭模型2Fig.7 Plane junction model 2

圖8 當(dāng)β = 180°時(shí)，平面接頭模型2 應(yīng)力奇性指數(shù)Fig.8 The singular index of stress in plane joint model 2 for β = 180°

圖9 當(dāng)β = 135°時(shí)，平面接頭模型2 應(yīng)力奇性指數(shù)Fig.9 The singular index of stress in plane joint model 2 for β = 135°

圖7中，在平面應(yīng)變條件下，當(dāng) β=90°，η =E(2)/E(1)=3時(shí)，由本文DQM 計(jì)算獲得平面接頭的前2 階應(yīng)力奇性指數(shù)分別為 λ1=?0.498 701 52和 λ2=?0.222 442 62，它們對(duì)應(yīng)的位移和應(yīng)力特征函數(shù)分布曲線如圖10、11 所示.從圖11 中看出，第2 階應(yīng)力奇性指數(shù) λ2對(duì) 應(yīng)的位移特征函數(shù)分布曲線和在粘接界面上連續(xù)且轉(zhuǎn)折，徑向應(yīng)力特征函數(shù)分布曲線在粘接界面上出現(xiàn)突變，與圖10 第1 階應(yīng)力奇性指數(shù) λ1對(duì)應(yīng)的特征函數(shù)曲線分布規(guī)律一致.

圖10 E(2)/E(1) = 3 時(shí)，平面接頭模型2 第1 階應(yīng)力奇性指數(shù)λ1 對(duì)應(yīng)的位移和應(yīng)力特征函數(shù)曲線圖Fig.10 Displacement and stress characteristic function curves corresponding to the 1st-order stress singular index λ1 of plane joint model 2 for E(2)/E(1) = 3

圖11 E(2)/E(1)=3 時(shí)，平面接頭模型2 第2 階應(yīng)力奇性指數(shù)λ2 對(duì)應(yīng)的位移和應(yīng)力特征函數(shù)曲線圖Fig.11 Displacement and stress characteristic function curves corresponding to the 2nd-order stress singular index λ2 of plane joint model 2 for E(2)/E(1)=3

2.3 算例3：平面接頭模型3

圖12 為兩個(gè)材料完全粘接在一起的平面接頭模型3，其中材料1 和材料2 均為各向同性材料，材料Poisson 比為 ν(1)=ν(2)=0.2，α=270°，β=90°.在每個(gè)材料域上離散單元總數(shù)皆取N=12 時(shí)，圖13 給出了η=E(2)/E(1)從0.001 至10 000 時(shí)，DQM 計(jì)算的應(yīng)力奇性指數(shù)的變化曲線，表明了兩相材料屬性關(guān)系對(duì)奇異性應(yīng)力狀況的影響情況，本文計(jì)算值與文獻(xiàn)[16]的結(jié)果完全一致.

圖12 平面接頭模型3Fig.12 Plane junction model 3

圖13 當(dāng)β = 90°時(shí)，平面接頭模型3 應(yīng)力奇性指數(shù)Fig.13 The singular index of stress in plane joint model 3 for β = 90°

3 結(jié) 論

DQM 在振動(dòng)工程領(lǐng)域應(yīng)用比較廣泛，主要用來求解結(jié)構(gòu)的自然頻率及振型，本文創(chuàng)新性地將該法運(yùn)用到彈性力學(xué)和工程斷裂力學(xué)的研究領(lǐng)域.基于平面接頭連接點(diǎn)處位移場和應(yīng)力場的Williams 漸近展開式，將其典型項(xiàng)代入平面問題平衡方程中，將平面接頭的線彈性理論微分方程轉(zhuǎn)換成一類ODEs 特征值問題，根據(jù)DQM 理論的計(jì)算列式，用 FORTRAN 編制了一個(gè)新求解器DQMEI 用于求解一般ODEs 特征值問題.最后應(yīng)用DQMEI 分析了平面接頭的應(yīng)力奇性指數(shù)，同時(shí)，也一并求出了平面接頭連接點(diǎn)處對(duì)應(yīng)的位移和應(yīng)力特征函數(shù)，這些特征函數(shù)在分析平面接頭完整應(yīng)力場和應(yīng)力強(qiáng)度因子時(shí)是不可或缺的物理量.文中給出數(shù)值算例的DQM 計(jì)算值與已有文獻(xiàn)結(jié)果完全一致，證明了求解一般平面接頭連接點(diǎn)處應(yīng)力奇性指數(shù)時(shí)，DQM 是一種有效且準(zhǔn)確的手段.