金曉清,呂鼎,張向?qū)?,李璞,周青華,胡玉梅
(1.重慶大學 機械傳動國家重點實驗室, 重慶400044; 2.四川大學 空天科學與工程學院, 四川 成都610065)
斷裂或接觸力學問題中第二類柯西奇異積分方程的一種解析方法
金曉清1,呂鼎1,張向?qū)?,李璞1,周青華2,胡玉梅1
(1.重慶大學 機械傳動國家重點實驗室, 重慶400044; 2.四川大學 空天科學與工程學院, 四川 成都610065)
第二類柯西奇異積分方程因涉及復奇異因子往往造成求解困難,而適用第一類奇異積分方程的高效數(shù)值方法并不能推廣至第二類奇異積分方程,即便是第二類奇異積分方程,其數(shù)值解法仍是一個難題. 為此提出了構造第二類奇異積分方程解析解的一種新方法. 通過分解柯西奇異項,并利用雅克比多項式的正交性,推導針對右端載荷項為單項式(monomial)的遞推解析解,進而借助級數(shù)展開的方法推廣至一般的載荷問題. 提出的基于遞推的解析解構造方案,能完美地結合maple軟件編程,從而提供一種方便、快捷、有效的算法. 由給出的算例可見,本方法適用于處理界面斷裂或接觸分析問題中含復數(shù)奇異因子的復雜情形,從而為研究該類典型力學問題提供了一種可供選擇的方法.
第二類奇異積分方程;柯西主值積分;復數(shù)奇異因子;界面裂紋
(1)
(2)
其中L為已知給定量.
在已有文獻中,第一類奇異積分方程求解方法相對成熟[7-10],第二類柯西奇異積分方程求解較難,要獲取其封閉解析解尤為困難[11]. 現(xiàn)有研究大多局限于奇異積分方程的數(shù)值解法,典型成果有:ERDOGAN等[7]和KRENK[8]采用的正交多項式法和GERASOULIS等[12]采用的分段多項式方法. MILLER等[13]采用分段二階多項式,以數(shù)值求解第二類奇異積分方程. JIN等[14-15]討論了上述各方法的利弊,并提出一種針對第二類奇異積分方程的有效數(shù)值解法. 周薇等[6]和周躍亭等[16]對第二類奇異積分方程的配置法、內(nèi)插型求積公式法和機械求積法等數(shù)值解法進行了論述.
對于第一類奇異積分方程,正交多項式法[7-10]求解效率高,并且易于編程. 然而,此方法很難推廣到第二類奇異積分方程[7],其所探討的求解方案[4]計算復雜、編程困難且效率不高.分段多項式法(piecewise polynomial approach)[12-13, 17-20]允許積分點和配置點(collocation point)任意分布,但當近似多項式的階次相對較小時,解的收斂速度和精度[18, 21]難以達到最優(yōu).雖然通過提高近似多項式的階次理論上可使函數(shù)解更精確,但其積分式將變得繁復難解[13, 17, 20]. 由此可見,第二類柯西奇異積分方程的數(shù)值解法較為煩瑣,至今仍缺少一種針對性強且高效的方法. 筆者借助maple軟件,通過分解柯西奇異項,利用雅克比多項式的正交性,提出一種實用的求解第二類奇異積分方程的解析新方法.
根據(jù)Muskhelishvili指數(shù)理論(index theory)[2],式(1)中φ(x)可以表示成新待求函數(shù)g(x)和反映問題物理奇異性態(tài)的基函數(shù)(fundamental function)ω(x)的乘積:
φ(x)=g(x)ω(x),
(3)
其中g(x)滿足H?lder連續(xù)條件,而
ω(x)=(1-x)α(1+x)β,
(4)
將式(3)(4)代入式(1),式(1)左端含柯西核的第2項可改寫為
(5)
(1-x)α(1+x)βcot(απ)+I(x),
(6)
其中,
(7)
(8)
在界面裂紋問題中,兩端裂尖具有物理奇異是一種典型的情況,此時α和β的實部取值滿足-1 當I(x)消失時,式(1)可化為 |a+bcot(απ)|ω(x)g(x)+ (9) 通過適當選取α值,式(9)中的第1項可化為0,即求: a+bcot(απ)=0, (10) 當系數(shù)a和b為復數(shù)時,α和β的值可由下式確定: (11) 其中整數(shù)N的取值需考慮問題于區(qū)間左端點-1處的物理奇異特性,即滿足-1 (12) 下面先討論載荷為單項式(monomial)即f(x)=xn的情況. 由式(3)(4)(11)可知,求解式(1)中的φ(x)只需求解式(12)中的g(x).聯(lián)立式(4)~(12),此時式(1)變換為如下形式: (13) 因載荷為n階高次項,可設g(t)為n+1階多項式:g(t)=c0+c1t+c2t2+c3t3+…+cn+1tn+1, (14) 故求解φ(x)即為求解g(t)中待定系數(shù)cn(n=0,1,2…)的值. 將多項式(14)代入,式(13)中積分項關鍵部分可變換為 (15) 其中di(t)表示式(15)中xi項的系數(shù),即 di(t)=ci+1+ci+2t+…+cn+1tn-i. (16) (17) (18) 雅克比多項式與權函數(shù)(1-t)α(1+t)β相關,存在如下正交特性: (19) 特別地,0階雅克比多項式等于1.由式(19)可知,非0階雅克比多項式存在如下性質(zhì): (20) (21) 將式(15)~(17)代入式(13),利用式(21),積分方程(13)可轉(zhuǎn)化為如下代數(shù)關系: (22) (23.0) (23.1) (23.i) (23.n) (24) 式(23)對應的關系轉(zhuǎn)變?yōu)?/p> (25.0) (25.1) … (25.i) … (25.n) g(t)=c0+c1t+c2t2+c3t3+c4t4. (26) (27) 即 φ(x)=ω(x)g(x)=-sinαπ(1-x)α(1+x)β× (28) 式(1)~(12)給出了一種求解第二類奇異積分方程的新方法.上述變換主要有2個目的: 回顧文獻[24]中界面裂紋的例子,其控制奇異積分方程為 (29) 滿足 (30) 式(29)中實常數(shù)γ與兩各向同性彈性介質(zhì)材料性質(zhì)有關,等式右邊的f(x)與裂紋面上的法向和切向載荷的綜合作用有關.指數(shù)α和β的值由式(11)確定: (31) 在此情況下,利用附錄A中給出的maple程序,得到邊界未知函數(shù)的封閉解: g(x)=-sin(πα)[3x5+(6α+4)x4+ (6α2+8α+1)x3+4α(α+1)2x2+c1x+c0], (32) 其中, (33) 表1為利用本方法求得的精確解與文獻[15]給出的數(shù)值計算結果的比較,兩者符合良好,在計算機的舍入誤差范圍內(nèi),驗證了本文給出的載荷f(x)為單項式xn的解析解. 且利用單項式載荷解析解構造多項式載荷解析解的方法可行.相應的maple程序見附錄. 第(ii)種載荷條件是泰勒級數(shù)展成多項式的組合. 將右端指數(shù)項進行泰勒展開,得到x的多項式,利用本方法近似求解. 隨著泰勒級數(shù)的增加,計算結果會逐漸收斂于精確解. 利用maple軟件并逐次增加多項式項數(shù),得到裂紋右尖端的復應力強度因子的計算值如表2所示,有效數(shù)收斂至小數(shù)點后10位. 從數(shù)值結果中可看出,收斂速度令人滿意. 通過分解柯西奇異項,消除奇異積分方程中的柯西奇異性,以便解析求解第二類奇異積分方程. 本解析方法可以用來處理界面裂紋中的物理奇異性問題,結合maple軟件編程,方便易行,對復數(shù)奇異因子同樣適用.從而,為獲取第二類柯西奇異積分方程的封閉解析解提供了一種切實可行的新途徑. 表1 本文精確解與文獻[15]數(shù)值解對比 表2 載荷(ii)情況下,例2中右裂尖應力強度因子(SIF)的數(shù)值收斂解 [1]GAKHOVFD.BoundaryValueProblems[M]. New York: Dover,1990. [2] MUSKHELISHVILI N I. Singular integral equations: Boundary problems of function theory and their application to mathematical physics[J].PNoordhoffN,2008,24(2):256-291. [3] HILLS D A, KELLY P A, DAI D N, et al. Solution of crack problems: The distributed dislocation technique[J].JournalofAppliedMechanics,1998,65(2):548. [4] 張耀明,孫翠蓮,谷巖.邊界積分方程中近奇異積分計算的一種變量替換法[J].力學學報,2008,40(2):207-214. ZHANG Y M, SUN C L, GU Y. The evaluation of nearly singular integrals in the boundary integral equations with variable transformation[J].ChineseJournalofTheoreticalandAppliedMechanics,2008,40(2):207-214. [5] 劉俊俏,段惠琴,李星.SH 波在壓電材料條中垂直界面裂紋處的散射[J].固體力學學報,2010,31(4),385-391. LIU J Q, DUAN H Q, LI X. The scattering of sh wave on a vertical crack in a coated piezoelectric strip[J] .ChineseJournalofSolidMechanics,2010,31(4),385-391. [6] 周薇.在接觸力學中的奇異積分方程的高精度數(shù)值解法[D].成都:電子科技大學,2011. ZHOU W.High-AccuracyNumericalSolutionforSingularIntegrationEquationsinContactMechanics[D]. Chengdu:University of Electronic Science and Technology of China,2011. [7] ERDOGAN F, GUPTA G D, COOK T.NumericalSolutionofSingularIntegralEquations[M].Berlin / Netherlands: Springer,1973,1(6):368-425. [8] KRENK S. On quadrature formulas for singular integral equations of the first and the second kind[J].QuarterlyofAppliedMathematics,1975,33(3): 225-232. [9] THEOCARIS P, IOAKIMIDIS N. Numerical integration methods for the solution of singular integral equations(for crack tip stress intensity factor evaluation in elastic media)[J].QuarterlyofAppliedMathematics,1977,35(1): 173-183. [10] ERDOGAN F, GUPTA G. On the numerical solution of singular integral equations[J].QuarterlyofAppliedMathematics,1972,29(4): 525-534. [11] 李星.積分方程[M].北京:科學出版社,2008. LI X.IntegralEquation[M]. Beijing: Science Press,2008. [12] GERASOULIS A, SRIVASTAV R. A method for the numerical solution of singular integral equations with a principal value integral[J].InternationalJournalofEngineeringScience,1981,19(9): 1293-1298. [13] MILLER G R, KEER L M. A numerical technique for the solution of singular integral equations of the second kind[J].QuarterlyofAppliedMathematics,1985,42(4): 455-465. [14] JIN X.AnalysisofSomeTwoDimensionalProblemsContainingCracksandHoles[D]. Chengdu: Northwestern University,2006. [15] JIN X, KEER L M, WANG Q. A practical method for singular integral equations of the second kind[J].EngineeringFractureMechanics,2008,75(5): 1005-1014. [16] 周躍亭,李星.具周期裂紋的半平面周期接觸問題的奇異積分方程數(shù)值解法[J].固體力學學報,2005,26(2): 167-174. ZHOU Y T, LI X. Singular integral equation method for periodic contact problem of an elastic half-plane with periodic cracks[J].ChineseJournalofSolidMechanics,2005,26(2): 167-174. [17] KIM P, LEE S. A piecewise linear quadrature of Cauchy singular integrals[J].JournalofComputationalandAppliedMathematics,1998,95(1/2): 101-115. [18] GERASOULIS A. Piecewise-polynomial quadratures for Cauchy singular integrals[J].SiamJournalonNumericalAnalysis,1986,23(4): 891-902. [19] GERASOULIS A. The use of piecewise quadratic polynomials for the solution of singular integral equations of Cauchy type[J].Computers&MathematicswithApplications,1982,8(1): 15-22. [20] KURTZ R D, FARRIS T N, SUN C. The numerical solution of Cauchy singular integral equations with application to fracture[J].InternationalJournalofFracture,1994,66(2): 139-154. [21] RABINOWITZ P. Convergence results for piecewise linear quadratures for Cauchy principal value integrals[J].MathematicsofComputation,1988,51(184): 741-747. [22] IOAKIMIDIS N I. On the numerical evaluation of derivatives of Cauchy principal value integrals[J].Computing,1981,27(1): 81-88. [23] ERDELYI A, MAGNUS W, OBERHETTINGER F, et al .TablesofIntegralTransforms:VolII[M]. New York / Toronto / London: McGraw-Hill,1954. [24] THEOCARIS P S, IOAKIMIDIS N I. On the numerical solution of Cauchy type singular integral equations and the determination of stress intensity factors in case of complex singularities[J].ZeitschriftFürAngewandteMathematikundPhysik,1977,28(6):1085-1098. 以本文的界面裂紋為例,說明maple程序的編制及調(diào)用. 當載荷為n-1階單項式時,可利用SIEslover子程序得到g′π(x)的解析式,其中g′π(x)=g(x)/cn. 將多項式載荷視為多個單項式的疊加,根據(jù)力學疊加原理,利用本程序可求解載荷條件(i)的情況. 同時,本程序也能求解載荷條件(ii)下α為復數(shù)的情況. >SIEslover:= proc (n,x, alpha) localpx,rx,i,j,Pj,cPj,cRj,c; c[n]:= 1; px:=c[n]; foritondo px:=px*x; rx:=px; forjfromiby -1 to 1 do Pj:= normal(simplify( JacobiP(j, alpha, -1-alpha,x), 'JacobiP')); cPj:= coeff(Pj,x^j); cRj:= coeff(rx,x^j); rx:= normal(simplify(rx-cRj*Pj/cPj)) od: c[n-i]:= -rx; px:=px+c[n-i] od: px:= collect(px, [x], factor) end: 下面為本求解器的調(diào)用示例: 載荷條件(i) >f5:= SIEslover(5,x, alpha); f4:=SIEslover (4,x, alpha); f1:=SIEslover (1,x, alpha); f1:=2α+x+1; 載荷條件(ii) >alpha:= -1/2-1/10*′I; >res2:= SIEslover (1,x, alpha); >forito 35 do c[i]:=coeff(taylor(2*exp(x*x),x= 0, 40),x^i); res2:=simplify(res2+c[i]* SIEslover (i+1,x, alpha)); sif:=evalf(subs(x= 1, res2), 20); k1:=Re(sif); k2:=Im(sif); if modp(i,5)=0 then printf("i=%3d,…k1 +Ik2 =%+15.10f,+I%+15.10f
",i,k1,k2) end if od: JIN Xiaoqing1, LYU Ding1, ZHANG Xiangning1, LI Pu1, ZHOU Qinghua2, HU Yumei1 (1.State Key Laboratory of Mechanical Transmission, Chongqing University, Chongqing 400044, China;2.School of Aeronautics and Astronautics, Sichuan University, Chengdu 610065, China) Due to the presence of complex singularity, solutions to the singular integration equation (SIE) of the second kind are still under development. As a matter of fact, numerical methods for SIE of the first kind are hardly applicable to SIE of the second kind. With the assistance of maple programming, this paper presents a novel approach to formulate an analytical solution to a typical SIE of the second kind. By splitting the Cauchy kernel, and taking advantage of the orthogonality of Jacobi polynomials, we derive an analytical solution corresponding to the monomial loading case. Furthermore, the solution to a general loading case may be obtained via series expansion. The present method appears efficient and convenient, providing an effective tool for treating tangentially loaded contact analyses and interface crack problems. SIE of the second kind;Cauchy principal value integration;complex singularity; interface crack 10.3785/j.issn.1008-9497.2017.05.009 O 343.3; O 346.1 :A :1008-9497(2017)05-548-07 2016-04-13. 國家自然科學基金資助項目 (51475057);中央高?;究蒲袠I(yè)務費專項(106112017CDJQJ328839). 金曉清(1974-),ORCID:http://orcid.org/0000-0002-8836-3505,博士,研究員,博士生導師,主要從事斷裂疲勞、細觀力學、摩擦學等研究,E-mail:jinxq@cqu.edu.cn. AnanalyticalmethodforsolvingCauchysingularintegralequationsofthesecondkindwithapplicationstofractureandcontactanalyses.Journal of Zhejiang University (Science Edition), 2017, 44(5 ): 548-5542 算例討論
3 結束語
附錄:maple程序