微分求積法分析平面接頭應力奇異性*

葛仁余，張佳宸，馬國強，劉小雙，牛忠榮

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

引 言

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

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

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

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

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

圖1為各向同性雙材料平面接頭模型，Γ1和 Γ3為 兩自由楔邊，Γ2為兩種不同材料的交界邊.在接頭連接點O處同時定義一個直角坐標系xOy和 一個極坐標系rOθ，交界 Γ2與x軸重合，θ逆時針方向為正，順時針方向為負.對于一個平面接頭構件，基于Williams 漸近展開理論[6]，將接頭連接點處的位移場表達成關于徑向r的級數(shù)漸近展開形式：

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

關于平面應力問題，將式（1）代入各向同性材料本構方程，獲得接頭連接點處應力分量如下：

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

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

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

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

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

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

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

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

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

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

當i=j時，有

式中，R(1)為

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

其中

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

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

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

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

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

表1 η = 3.0 時，平面接頭模型1 的第1 階應力奇性指數(shù)λ1 計算值隨離散單元數(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

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

表2 η = 4.0 時，平面接頭模型1 的第1 階應力奇性指數(shù)λ1 計算值隨離散單元數(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 時，平面接頭模型1 的第1 階應力奇性指數(shù)λ1 計算值隨離散單元數(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 當β = 90°時，平面接頭模型1 的應力奇性指數(shù)Fig.3 The singular index of stress in plane joint model 1 for β = 90°

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

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

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

圖6 E(2)/E(1) = 0.1 時，平面接頭模型1 第1 階應力奇性指數(shù)λ1 對應的位移和應力特征函數(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 當β = 180°時，平面接頭模型2 應力奇性指數(shù)Fig.8 The singular index of stress in plane joint model 2 for β = 180°

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

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

圖10 E(2)/E(1) = 3 時，平面接頭模型2 第1 階應力奇性指數(shù)λ1 對應的位移和應力特征函數(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 時，平面接頭模型2 第2 階應力奇性指數(shù)λ2 對應的位移和應力特征函數(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 為兩個材料完全粘接在一起的平面接頭模型3，其中材料1 和材料2 均為各向同性材料，材料Poisson 比為 ν(1)=ν(2)=0.2，α=270°，β=90°.在每個材料域上離散單元總數(shù)皆取N=12 時，圖13 給出了η=E(2)/E(1)從0.001 至10 000 時，DQM 計算的應力奇性指數(shù)的變化曲線，表明了兩相材料屬性關系對奇異性應力狀況的影響情況，本文計算值與文獻[16]的結果完全一致.

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

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

3 結 論

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