葛仁余,張佳宸,馬國(guó)強(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)力場(chǎ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ù)簡(jiǎn)便而十分有效的方法[3-5].Williams[6](1952年)用直接的Airy 應(yīng)力函數(shù)法研究了冪應(yīng)力奇異性.England[7](1971年),Stern 和Soni[8](1976年),以及Carpenter 和Byers[9](1987年)使用了簡(jiǎn)單的復(fù)勢(shì)方法來研究?jī)鐟?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í),由于該理論不依賴于變分原理,具有公式簡(jiǎn)單、編程方便、高效、精度高等特點(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為各向同性雙材料平面接頭模型,Γ1和 Γ3為 兩自由楔邊,Γ2為兩種不同材料的交界邊.在接頭連接點(diǎn)O處同時(shí)定義一個(gè)直角坐標(biāo)系xOy和 一個(gè)極坐標(biāo)系rOθ,交界 Γ2與x軸重合,θ逆時(shí)針方向?yàn)檎?,順時(shí)針方向?yàn)樨?fù).對(duì)于一個(gè)平面接頭構(gòu)件,基于Williams 漸近展開理論[6],將接頭連接點(diǎn)處的位移場(chǎng)表達(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)分布方式采用等步長(zhǎng)均勻網(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 語(yǔ)言編寫了通用程序DQMEI 用于DQM 計(jì)算平面接頭應(yīng)力奇性指數(shù) λk及其相應(yīng)的位移特征函數(shù)(θ)和 應(yīng)力特征函數(shù)(θ).
圖2為平面接頭模型1,α=180°是 一定值,β角是變化的.考慮平面應(yīng)變情況,并且假定兩種材料的Poisson 比ν(1)=ν(2)=0.2 ;離散單元總數(shù)分別取N=4,6,8,10,12,離散節(jié)點(diǎn)分布形式采用等步長(zhǎng)均勻網(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
圖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)力場(chǎ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
圖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°
DQM 在振動(dòng)工程領(lǐng)域應(yīng)用比較廣泛,主要用來求解結(jié)構(gòu)的自然頻率及振型,本文創(chuàng)新性地將該法運(yùn)用到彈性力學(xué)和工程斷裂力學(xué)的研究領(lǐng)域.基于平面接頭連接點(diǎn)處位移場(chǎng)和應(yīng)力場(chǎ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)力場(chǎng)和應(yīng)力強(qiáng)度因子時(shí)是不可或缺的物理量.文中給出數(shù)值算例的DQM 計(jì)算值與已有文獻(xiàn)結(jié)果完全一致,證明了求解一般平面接頭連接點(diǎn)處應(yīng)力奇性指數(shù)時(shí),DQM 是一種有效且準(zhǔn)確的手段.