文章快速检索     高级检索
  浙江大学学报(理学版)  2017, Vol. 44 Issue (5): 548-554  DOI:10.3785/j.issn.1008-9497.2017.05.009
0

引用本文 [复制中英文]

金晓清, 吕鼎, 张向宁, 李璞, 周青华, 胡玉梅. 断裂或接触力学问题中第二类柯西奇异积分方程的一种解析方法[J]. 浙江大学学报(理学版), 2017, 44(5): 548-554. DOI: 10.3785/j.issn.1008-9497.2017.05.009.
[复制中文]
JIN Xiaoqing, LYU Ding, ZHANG Xiangning, LI Pu, ZHOU Qinghua, HU Yumei. An analytical method for solving Cauchy singular integral equations of the second kind with applications to fracture and contact analyses[J]. Journal of Zhejiang University(Science Edition), 2017, 44(5): 548-554. DOI: 10.3785/j.issn.1008-9497.2017.05.009.
[复制英文]

基金项目

国家自然科学基金资助项目(51475057);中央高校基本科研业务费专项(106112017CDJQJ328839)

作者简介

金晓清(1974-), ORCID:http://orcid.org/0000-0002-8836-3505, 博士, 研究员, 博士生导师, 主要从事断裂疲劳、细观力学、摩擦学等研究, E-mail: jinxq@cqu.edu.cn

文章历史

收稿日期:2016-04-13
断裂或接触力学问题中第二类柯西奇异积分方程的一种解析方法
金晓清1 , 吕鼎1 , 张向宁1 , 李璞1 , 周青华2 , 胡玉梅1     
1. 重庆大学 机械传动国家重点实验室, 重庆 400044;
2. 四川大学 空天科学与工程学院, 四川 成都 610065
摘要: 第二类柯西奇异积分方程因涉及复奇异因子往往造成求解困难,而适用第一类奇异积分方程的高效数值方法并不能推广至第二类奇异积分方程,即便是第二类奇异积分方程,其数值解法仍是一个难题.为此提出了构造第二类奇异积分方程解析解的一种新方法.通过分解柯西奇异项,并利用雅克比多项式的正交性,推导针对右端载荷项为单项式(monomial)的递推解析解,进而借助级数展开的方法推广至一般的载荷问题.提出的基于递推的解析解构造方案,能完美地结合maple软件编程,从而提供一种方便、快捷、有效的算法.由给出的算例可见,本方法适用于处理界面断裂或接触分析问题中含复数奇异因子的复杂情形,从而为研究该类典型力学问题提供了一种可供选择的方法.
关键词: 第二类奇异积分方程    柯西主值积分    复数奇异因子    界面裂纹    
An analytical method for solving Cauchy singular integral equations of the second kind with applications to fracture and contact analyses
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
Abstract: 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.
Key words: SIE of the second kind    Cauchy principal value integration    complex singularity    interface crack    
0 引言

在接触力学、断裂力学和量子理论学等问题中,含柯西核(Cauchy kernel)的奇异积分方程十分常见[1-6].典型的第二类柯西奇异积分方程:

$ a\varphi \left( x \right) + \frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {\frac{{\varphi \left( t \right)}}{{t - x}}{\rm{d}}t} = f\left( x \right), - 1 < x < 1, $ (1)

式(1) 左端积分项所包含的$ \frac{1}{{t - x}} $称为柯西核,而该积分项需要在“柯西主值(Cauchy principal value)” [3]意义下计算.右端项f(x)为已知函数且通常假定满足Hölder连续条件[2]; 系数ab为给定的复常数.如果a=0,未知函数φ(x)只出现在积分项中,式(1) 退化为第一类奇异积分方程.假定与式(1) 相关联的边界条件待求函数φ(x)满足

$ \int_{ - 1}^1 {\varphi \left( x \right){\rm{d}}x} = L, $ (2)

其中L为已知给定量.

在已有文献中,第一类奇异积分方程求解方法相对成熟[7-10], 第二类柯西奇异积分方程求解较难,要获取其封闭解析解尤为困难[11].现有研究大多局限于奇异积分方程的数值解法,典型成果有:ERDOGAN等[7]和KRENK[8]采用的正交多项式法和GERASOULIS等[12]采用的分段多项式方法. MILLER等[13]采用分段二阶多项式,以数值求解第二类奇异积分方程. JIN等[14-15]讨论了上述各方法的利弊,并提出一种针对第二类奇异积分方程的有效数值解法.周薇等[6]和周跃亭等[16]对第二类奇异积分方程的配置法、内插型求积公式法和机械求积法等数值解法进行了论述.

对于第一类奇异积分方程,正交多项式法[7-10]求解效率高,并且易于编程.然而,此方法很难推广到第二类奇异积分方程[7],其所探讨的求解方案[4]计算复杂、编程困难且效率不高.分段多项式法(piecewise polynomial approach)[12-13, 17-20]允许积分点和配置点(collocation point)任意分布,但当近似多项式的阶次相对较小时,解的收敛速度和精度[18, 21]难以达到最优.虽然通过提高近似多项式的阶次理论上可使函数解更精确,但其积分式将变得繁复难解[13, 17, 20].由此可见,第二类柯西奇异积分方程的数值解法较为烦琐,至今仍缺少一种针对性强且高效的方法.笔者借助maple软件,通过分解柯西奇异项,利用雅克比多项式的正交性,提出一种实用的求解第二类奇异积分方程的解析新方法.

1 解析解计算方法

根据Muskhelishvili指数理论(index theory) [2],式(1) 中φ(x)可以表示成新待求函数g(x)和反映问题物理奇异性态的基函数(fundamental function)ω(x)的乘积:

$ \varphi \left( x \right) = g\left( x \right)\omega \left( x \right), $ (3)

其中g(x)满足Hölder连续条件,而

$ \omega \left( x \right) = {\left( {1 - x} \right)^\alpha }{\left( {1 - x} \right)^\beta }, $ (4)

将式(3)(4) 代入式(1),式(1) 左端含柯西核的第2项可改写为

$ \begin{array}{l} \frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {\frac{{\varphi \left( t \right)}}{{t - x}}{\rm{d}}t} = \frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {{{\left( {1 - t} \right)}^\alpha }{{\left( {1 - t} \right)}^\beta } \times } \\ \frac{{g\left( t \right) - g\left( x \right)}}{{t - x}}{\rm{d}}t + g\left( x \right)\frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {\frac{{{{\left( {1 - t} \right)}^\alpha }{{\left( {1 - t} \right)}^\beta }}}{{t - x}}{\rm{d}}t} . \end{array} $ (5)

经过上述变换,将柯西奇异项分解成两部分,其中,式(5) 右端第1项中的柯西奇异性被消除.论述如下:根据Hölder条件,假设g(x)在[-1, 1]内p+1次连续可导,由泰勒级数可知,$ \frac{{g\left( t \right) - g\left( x \right)}}{{t - x}} $xt至少p次局部连续可导[22].此外,对式(5) 右端第2项积分可求得封闭解(closed-form solution).根据已知结果[7, 23]有:当x<1, Re(α)>-1, Re(β)>-1, α≠0, 1,…,时,

$ \begin{array}{l} \frac{1}{{\rm{\pi }}}\int_{ - 1}^1 {\frac{{{{\left( {1 - t} \right)}^\alpha }{{\left( {1 - t} \right)}^\beta }}}{{t - x}}{\rm{d}}t} = \\ \;\;\;\;\;\;\;\;\;{\left( {1 - x} \right)^\alpha }{\left( {1 - x} \right)^\beta }\cot \left( {\alpha {\rm{\pi }}} \right) + I\left( x \right), \end{array} $ (6)

其中,

$ \begin{array}{l} I\left( x \right) = - \frac{{\Gamma \left( {\beta + 1} \right)\Gamma \left( \alpha \right)}}{{{\rm{\pi }}\Gamma \left( {\alpha + \beta + 1} \right)}}{2^{\alpha + \beta }} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;F\left( { - \alpha - \beta ,1;1 - \alpha ;\frac{{1 - x}}{2}} \right), \end{array} $ (7)

式(7) 中Γ是伽马函数,F是超几何函数,满足如下关系:

$ I\left( x \right) = \left\{ \begin{array}{l} 0,\;\;\;\;\alpha + \beta = - 1,\\ - {\rm{\pi csc}}\left( {{\rm{\pi }}\alpha } \right),\;\;\;\;\alpha + \beta = 0,\\ {\rm{\pi csc}}\left( {{\rm{\pi }}\alpha } \right)\left( {2\alpha - 1 - x} \right),\;\;\;\;\alpha + \beta = 1. \end{array} \right. $ (8)

在界面裂纹问题中,两端裂尖具有物理奇异是一种典型的情况,此时αβ的实部取值满足-1<Re(α)<0和-1<Re(β)<0,且α+β=-1,即对应式(8) 的第1类情况.本文以此为例,即I(x)=0,其他2种情况亦可通过类似方法求解.

I(x)消失时,式(1) 可化为

$ \begin{array}{l} \left| {a + b\cot \left( {\alpha {\rm{\pi }}} \right)} \right|\omega \left( x \right)g\left( x \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {\omega \left( t \right)\frac{{g\left( t \right) - g\left( x \right)}}{{t - x}}{\rm{d}}t} = f\left( x \right), \end{array} $ (9)

通过适当选取α值,式(9) 中的第1项可化为0,即求:

$ a + b\cot \left( {\alpha {\rm{\pi }}} \right) = 0, $ (10)

当系数ab为复数时,αβ的值可由下式确定:

$ \alpha = \frac{1}{{2{\rm{\pi i}}}}\ln \left( {\frac{{a - {\rm{i}}b}}{{a + {\rm{i}}b}}} \right) + N,\;\;\;\;\beta = - 1 - \alpha , $ (11)

其中整数N的取值需考虑问题于区间左端点-1处的物理奇异特性,即满足-1<Re(α)<0.当系数ab为实数时,由式(10) 确定的指数α也是实数,与式(11) 一致.据式(10),原第二类奇异积分方程(1) 可进一步简化为

$ \frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {\omega \left( t \right)\frac{{g\left( t \right) - g\left( x \right)}}{{t - x}}{\rm{d}}t} = f\left( x \right), - 1 < x < 1. $ (12)

下面先讨论载荷为单项式(monomial)即f(x)=xn的情况.由式(3)(4)(11) 可知,求解式(1) 中的φ(x)只需求解式(12) 中的g(x).联立式(4)~(12),此时式(1) 变换为如下形式:

$ \begin{array}{*{20}{c}} {\frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {\frac{{{{\left( {1 - t} \right)}^\alpha }{{\left( {1 - t} \right)}^\beta }\left[ {g\left( t \right) - g\left( x \right)} \right]}}{{t - x}}{\rm{d}}t} = {x^n},}\\ { - 1 < x < 1,} \end{array} $ (13)

因载荷为n阶高次项,可设g(t)为n+1阶多项式:

$ g\left( t \right) = {c_0} + {c_1}t + {c_2}{t^2} + {c_3}{t^3} + \cdots + {c_{n + 1}}{t^{n + 1}}, $ (14)

故求解φ(x)即为求解g(t)中待定系数cn(n=0, 1, 2…)的值.将多项式(14) 代入,式(13) 中积分项关键部分可变换为

$ \begin{array}{*{20}{c}} {\frac{{g\left( t \right) - g\left( x \right)}}{{t - x}} = {c_{n + 1}}{x^n} + \left( {{c_n} + {c_{n + 1}}t} \right){x^{n - 1}} + \cdots + \left( {{c_1} + } \right.}\\ {\left. {{c_2}t + \cdots + {c_{i + 1}}{t^i} + \cdots + {c_{n + 1}}{t^n}} \right) = \sum\limits_{i = 0}^n {{d_i}\left( t \right){x^i}} ,} \end{array} $ (15)

其中di(t)表示式(15) 中xi项的系数,即

$ {d_i}\left( t \right) = {c_{i + 1}} + {c_{i + 2}}t + \cdots + {c_{n + 1}}{t^{n - i}}. $ (16)

由式(16) 可知,di(t)由di+1(t)乘以t加上常数项ci+1得到.求解系数cn(n=0, 1, …, 4) 的一种简单方法是将式(15) 中每一项di(t)扩展成雅克比多项式(Jacobi polynomial)Pk(α, β)(t)的形式,即

$ \begin{array}{*{20}{c}} {{d_i}\left( t \right) = \lambda _i^0P_0^{\left( {\alpha ,\beta } \right)}\left( t \right) + \lambda _i^1P_1^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + }\\ {\lambda _i^kP_k^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + \lambda _i^{n - i}P_{n - i}^{\left( {\alpha ,\beta } \right)}\left( t \right),} \end{array} $ (17)

其中λik(k=0, 1, …, n-i)为与k阶雅克比多项式相关的展开系数.0~3阶雅克比多项式示例如下:

$ \left\{ \begin{array}{l} P_0^{\left( {\alpha , - 1 - \alpha } \right)}\left( t \right) = 1;\\ P_1^{\left( {\alpha , - 1 - \alpha } \right)}\left( t \right) = \frac{1}{2}\left( {t + 2\alpha + 1} \right);\\ \begin{array}{*{20}{c}} {P_2^{\left( {\alpha , - 1 - \alpha } \right)}\left( t \right) = \frac{1}{2}\left( {\alpha + 2} \right)\left( {\alpha + 1} \right) + }\\ {\left( {\alpha + 2} \right)\left( {t - 1} \right) + \frac{3}{4}{{\left( {t - 1} \right)}^2};} \end{array}\\ P_3^{\left( {\alpha , - 1 - \alpha } \right)}\left( t \right) = \frac{1}{6}\left( {\alpha + 1} \right)\left( {\alpha + 2} \right)\left( {\alpha + 3} \right) + \\ \;\;\;\;\;\;\frac{3}{4}\left( {\alpha + 2} \right)\left( {\alpha + 3} \right)\left( {t - 1} \right) + \\ \;\;\;\;\;\;\frac{3}{2}\left( {\alpha + 3} \right){\left( {t - 1} \right)^2} + \frac{5}{4}{\left( {t - 1} \right)^3};\\ \;\;\;\;\;\; \vdots \end{array} \right. $ (18)

雅克比多项式与权函数(1-t)α(1+t)β相关,存在如下正交特性:

$ \begin{array}{*{20}{c}} {\int_{ - 1}^1 {{{\left( {1 - t} \right)}^\alpha }{{\left( {1 + t} \right)}^\beta }P_m^{\left( {\alpha ,\beta } \right)}\left( t \right)P_n^{\left( {\alpha ,\beta } \right)}\left( t \right) = 0,} }\\ {m \ne n;\;\;\;m,n = 0,1,2, \cdots .} \end{array} $ (19)

特别地,0阶雅克比多项式等于1.由式(19) 可知,非0阶雅克比多项式存在如下性质:

$ \int_{ - 1}^1 {{{\left( {1 - t} \right)}^\alpha }{{\left( {1 + t} \right)}^\beta }P_m^{\left( {\alpha ,\beta } \right)}\left( t \right) = 0} ,\;\;\;m = 1,2, \cdots . $ (20)

利用式(20) 性质,并考虑式(15)~(17),di(t)关于雅克比权函数的积分项可化简至仅保留0阶雅克比项,得到如下重要结论:

$ \begin{array}{l} \int_{ - 1}^1 {{{\left( {1 - t} \right)}^\alpha }{{\left( {1 + t} \right)}^\beta }{d_i}\left( t \right){\rm{d}}t} = \lambda _i^0\int_{ - 1}^1 {{{\left( {1 - t} \right)}^\alpha }{{\left( {1 + t} \right)}^\beta }{\rm{d}}t} = \\ \;\;\;\;\;\lambda _i^0{2^{\alpha + \beta + 1}}\Gamma \left( {\alpha + 1} \right)\Gamma \left( {\beta + 1} \right)/\Gamma \left( {\alpha + \beta + 2} \right) = \\ \;\;\;\;\;\lambda _i^0\Gamma \left( {\alpha + 1} \right)\Gamma \left( { - \alpha } \right) = - \frac{{\rm{\pi }}}{{\sin \alpha {\rm{\pi }}}}\lambda _i^0. \end{array} $ (21)

将式(15)~(17) 代入式(13),利用式(21),积分方程(13) 可转化为如下代数关系:

$ \begin{array}{*{20}{c}} {\frac{b}{{\rm{\pi }}}\int_{ - 1}^1 {\frac{{{{\left( {1 - t} \right)}^\alpha }{{\left( {1 + t} \right)}^\beta }\left[ {g\left( t \right) - g\left( x \right)} \right]}}{{t - x}}{\rm{d}}t} = }\\ { - \frac{b}{{\sin \alpha {\rm{\pi }}}}\sum\limits_{i = 0}^n {\lambda _i^0{x^i}} = {x^n}.} \end{array} $ (22)

由式(22) 知,λn0=-(sin απ)/b,当i=0, 1, 2, …, n-1时,求和公式中的系数λi0均为0,即λi0=0(i=0, 1, 2, …, n-1).

雅克比多项式Pk(α, β)(t)可由maple的内置子程序得到,故利用maple软件能十分方便地计算此类问题.式(16) 和(17) 为关于t等价的ni次多项式,根据两式对应的系数关系及λi0=0(i=0, 1, 2…, n-1),ci与雅克比多项式存在如下关系:

$ {c_{n + 1}} = \lambda _n^0P_0^{\left( {\alpha ,\beta } \right)}\left( t \right), $ (23.0)
$ {c_n} + {c_{n + 1}}t = \lambda _{n - 1}^1P_1^{\left( {\alpha ,\beta } \right)}\left( t \right), $ (23.1)
$ \begin{array}{*{20}{c}} \cdots \\ {{c_{n + 1 - i}} + {c_{n + 2 - i}}t + \cdots + {c_{n + 1}}{t^i} = \lambda _{n - i}^1P_1^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + }\\ {\lambda _{n - i}^kP_k^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + \lambda _{n - i}^iP_i^{\left( {\alpha ,\beta } \right)}\left( t \right),}\\ \cdots \end{array} $ (23.i)
$ \begin{array}{*{20}{c}} {{c_1} + {c_2}t + {c_3}{t^2} + \cdots + {c_{n + 1}}{t^n} = \lambda _0^1P_1^{\left( {\alpha ,\beta } \right)}\left( t \right) + }\\ {\lambda _0^2P_2^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + \lambda _0^nP_n^{\left( {\alpha ,\beta } \right)}\left( t \right),} \end{array} $ (23.n)

上述关系式中,式(23.i+1) 等号左端的式子由式(23.i)等号左边的式子乘以t加上常数项cn-i得到.由式(17) 与(23) 知,等号右端为左端函数的雅克比多项式扩展形式,其中,P0(α, β)(t)=1,λn0=-(sin απ)/b,故dn(t)中cn+1=λn0=-(sin απ)/b.此时,式(23.0) 中的系数全部求得.代入式(18)P1(α.β)(t)中,根据t系数的对应关系,进而可计算得λn-11的值,对比(23.1) 的常数项可求得cn的值.至此式(23.1) 的未知系数全部求得.依此类推,通过单项式ti(i=0, 1, 2, …, n)的系数对应关系, 可以逐次求得式(23.n)中所有的未知系数ci(i=0, 1, 2, …, n),进而求得函数g(t)的解析式.已知cn+1=λn0=-(sin απ)/b,结合式(3)(4) 及(14),为方便编程,可将式(3) 中xn+1的系数提出,即

$ \begin{array}{l} \varphi \left( x \right) = \omega \left( x \right)g\left( x \right) = \frac{{ - \sin \alpha {\rm{\pi }}}}{b}{\left( {1 - x} \right)^\alpha }{\left( {1 + x} \right)^\beta } \times \\ \;\;\;\;\;\;\;\;\;\;\;\left( {{x^{n + 1}} + {{c'}_n}{x^n} + \cdots + {{c'}_i}{x^i} + \cdots + {{c'}_0}} \right), \end{array} $ (24)

其中, ci′=ci/λn0cn+1=1,ci/cn+1=ci/cn+1

式(23) 对应的关系转变为

$ {{c'}_{n + 1}} = P_0^{\left( {\alpha ,\beta } \right)}\left( t \right) = 1, $ (25.0)
$ {{c'}_n} + {{c'}_{n + 1}}t = \bar \lambda _{n - 1}^1P_1^{\left( {\alpha ,\beta } \right)}\left( t \right), $ (25.1)
$ \begin{array}{*{20}{c}} \cdots \\ {{{c'}_{n + 1 - i}} + {{c'}_{n + 2 - i}}t + \cdots + {{c'}_{n + 1}}{t^i} = \bar \lambda _{n - i}^1P_1^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + }\\ {\bar \lambda _{n - i}^kP_k^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + \bar \lambda _{n - i}^iP_i^{\left( {\alpha ,\beta } \right)}\left( t \right),}\\ \cdots \end{array} $ (25.i)
$ \begin{array}{*{20}{c}} {{{c'}_1} + {{c'}_2}t + {{c'}_3}{t^2} + \cdots + {{c'}_{n + 1}}{t^n} = \bar \lambda _0^1P_1^{\left( {\alpha ,\beta } \right)}\left( t \right) + }\\ {\bar \lambda _0^2P_2^{\left( {\alpha ,\beta } \right)}\left( t \right) + \cdots + \bar \lambda _0^nP_n^{\left( {\alpha ,\beta } \right)}\left( t \right),} \end{array} $ (25.n)

其中λik=λik/λn0,且cn+1=1.由上述对应关系,可得到所有ci(i=0, 1, 2, …, n)的值,具体可参考本文附录提供的maple程序.

2 算例讨论

a=1,b=1,f(x)=x3为例,举例验证本文提出的方法.由式(10) 可知$ \alpha = - \frac{1}{4} $,且β须满足α+β=-1.由此可设g(t)是一个4阶多项式

$ g\left( t \right) = {c_0} + {c_1}t + {c_2}{t^2} + {c_3}{t^3} + {c_4}{t^4}. $ (26)

结合式(10)~(25),通过ci(i=1, 2, 3, 4) 与λik对应的比例关系,可逐步求得系数cn(n=1, 2, 3, 4).

$ \left\{ \begin{array}{l} {c_4} = - \sin \alpha {\rm{\pi ,}}\;\;\;{c_3} = \left( {2\alpha + 1} \right){c_4},\\ {c_2} = 2\alpha \left( {\alpha + 1} \right){c_4},\;\;\;{c_1} = \frac{2}{3}\alpha \left( {\alpha + 1} \right)\left( {2\alpha + 1} \right){c_4}, \end{array} \right. $ (27)

$ \begin{array}{*{20}{c}} {\varphi \left( x \right) = \omega \left( x \right)g\left( x \right) = - \sin \alpha {\rm{\pi }}{{\left( {1 - x} \right)}^\alpha }{{\left( {1 + x} \right)}^\beta } \times }\\ {\left( {{x^4} + \frac{1}{2}{x^3} - \frac{3}{8}{x^2} - \frac{1}{{16}}x + {c_0}} \right),} \end{array} $ (28)

c0φ(x)积分的齐次条件相关,其中$ {c_0} = \frac{2}{3}\alpha \left( {\alpha + 1} \right)\left( {{\alpha ^2} + \alpha + 1} \right) - \frac{{L\sin \;\alpha {\rm{\pi }}}}{{\rm{\pi }}} $,将求得的φ(x)代入式(1),利用maple验证了所得φ(x)的正确性.

式(1)~(12) 给出了一种求解第二类奇异积分方程的新方法.上述变换主要有2个目的:

1) 将$ \frac{{g\left( t \right) - g\left( x \right)}}{{t - x}} $看作一个新的未知函数,式(1) 所示的第二类奇异积分问题可转化为式(12) 所示的第一类奇异积分,之后将$ \frac{{\left[{g\left( t \right)-g\left( x \right)} \right]}}{{t - x}} $表示成雅克比多项式的叠加,便可很容易地求解当前问题.这种求解方法与ERDOGAN等[7]求解第二类奇异积分方程的方法类似,但本方法更简单.

2) 本方法的解析推导部分也可由数值计算代替,因而可开发一种实用又无须复杂推导的算法来求解第二类奇异积分方程[15].本问题的难点之一是处理由基函数和柯西核所引起的奇异性.因$ \frac{{g\left( t \right) - g\left( x \right)}}{{t - x}} $是一个平滑函数,通过分解柯西核,式(5) 中柯西项的奇异性被消除,故求积的难度只剩下处理基函数所引起的奇异性.当本问题涉及其他指数类型的基函数时,上述消除奇异性的方法依然可用,这为求解涉及4类物理问题的第二类奇异积分方程提供了数值方法.

回顾文献[24]中界面裂纹的例子,其控制奇异积分方程为

$ - {\rm{i}}\gamma \varphi \left( x \right) + \frac{1}{{\rm{\pi }}}\int_{ - 1}^1 {\frac{{\varphi \left( t \right)}}{{t - x}}{\rm{d}}t = f\left( x \right)} , $ (29)

满足

$ \int_{ - 1}^1 {\varphi \left( t \right){\rm{d}}t = 0} . $ (30)

式(29) 中实常数γ与两各向同性弹性介质材料性质有关,等式右边的f(x)与裂纹面上的法向和切向载荷的综合作用有关.指数αβ的值由式(11) 确定:

$ \alpha = - \frac{1}{2} - {\rm{i}}\psi ,\beta = - \frac{1}{2} + {\rm{i}}\psi , $ (31)

其中$ \psi = \frac{1}{{2{\rm{\pi }}}}\ln \frac{{1 + \gamma }}{{1 - \gamma }} $,在数值计算中,ψ取为0.1,取两类载荷条件为(ⅰ)f(x)=1+x3+3x4和(ⅱ)f(x)= 2ex2.

在此情况下,利用附录A中给出的maple程序,得到边界未知函数的封闭解:

$ \begin{array}{l} g\left( x \right) = - \sin \left( {{\rm{\pi }}\alpha } \right)\left[ {3{x^5} + \left( {6\alpha + 4} \right){x^4} + } \right.\\ \;\;\;\;\;\;\left. {\left( {6{\alpha ^2} + 8\alpha + 1} \right){x^3} + 4\alpha {{\left( {\alpha + 1} \right)}^2}{x^2} + {c_1}x + {c_0}} \right], \end{array} $ (32)

其中,

$ \left\{ \begin{array}{l} {c_1} = 2{\alpha ^4} + \frac{{16}}{3}{\alpha ^3} + 6{\alpha ^2} + \frac{8}{3}\alpha + 1,\\ {c_0} = \frac{5}{4}{\alpha ^5} + \frac{8}{3}{\alpha ^4} + \frac{{16}}{3}{\alpha ^3} + \frac{{16}}{3}{\alpha ^2} + \frac{{58}}{{15}}\alpha + 1. \end{array} \right. $ (33)

表 1为利用本方法求得的精确解与文献[15]给出的数值计算结果的比较,两者符合良好,在计算机的舍入误差范围内,验证了本文给出的载荷f(x)为单项式xn的解析解.且利用单项式载荷解析解构造多项式载荷解析解的方法可行.相应的maple程序见附录.

表 1 本文精确解与文献[15]数值解对比 Table 1 Comparison of present exact values with numerical values of ref.[15]

第(ⅱ)种载荷条件是泰勒级数展成多项式的组合.将右端指数项进行泰勒展开,得到x的多项式,利用本方法近似求解.随着泰勒级数的增加,计算结果会逐渐收敛于精确解.利用maple软件并逐次增加多项式项数,得到裂纹右尖端的复应力强度因子的计算值如表 2所示,有效数收敛至小数点后10位.从数值结果中可看出,收敛速度令人满意.

表 2 载荷(ⅱ)情况下,例2中右裂尖应力强度因子(SIF)的数值收敛解 Table 2 Numerical convergent results of the normalized SIF at the right crack tip of the second example for the loading condition (ⅱ)
3 结束语

通过分解柯西奇异项,消除奇异积分方程中的柯西奇异性,以便解析求解第二类奇异积分方程.本解析方法可以用来处理界面裂纹中的物理奇异性问题,结合maple软件编程,方便易行,对复数奇异因子同样适用.从而,为获取第二类柯西奇异积分方程的封闭解析解提供了一种切实可行的新途径.

附录:maple程序

以本文的界面裂纹为例,说明maple程序的编制及调用.当载荷为n-1阶单项式时,可利用SIEslover子程序得到g′π(x)的解析式,其中g′π(x)=g(x)/cn.将多项式载荷视为多个单项式的叠加,根据力学叠加原理,利用本程序可求解载荷条件(ⅰ)的情况.同时,本程序也能求解载荷条件(ⅱ)下α为复数的情况.

>SIEslover:= proc (n, x, alpha)

  local px, rx, i, j, Pj, cPj, cRj, c;

  c[n]:= 1;

  px:= c[n];

  for i to n do

    px:= px*x;

    rx:= px;

    for j from i by -1 to 1 do

      Pj:= normal(simplify(Jacobi P(j, alpha, -1-alpha, x), 'Jacobi P'));

      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:

下面为本求解器的调用示例:

载荷条件(ⅰ)

>f5:= SIEslover(5, x, alpha);

f4:=SIEslover (4, x, alpha);

f1:=SIEslover (1, x, alpha);

$ \begin{array}{l} {\rm{f}}5: = {x^5} + \left( {2\alpha + 1} \right){x^4} + 2\alpha \left( {\alpha + 1} \right){x^3} + \frac{2}{3}\alpha \left( {\alpha + 1} \right) \times \\ \;\;\;\;\;\;\left( {2\alpha + 1} \right){x^2} + \frac{2}{3}\alpha \left( {\alpha + 1} \right)\left( {{\alpha ^2} + \alpha + 1} \right)x + \\ \;\;\;\;\;\;\frac{2}{{15}}\alpha \left( {\alpha + 1} \right)\left( {2\alpha + 1} \right)\left( {{\alpha ^2} + \alpha + 3} \right); \end{array} $
$ \begin{array}{l} {\rm{f}}4: = {x^4} + \left( {2\alpha + 1} \right){x^3} + 2\alpha \left( {\alpha + 1} \right){x^2} + \frac{2}{3}\alpha \left( {\alpha + 1} \right)\\ \;\;\;\;\;\;\left( {2\alpha + 1} \right)x + \frac{2}{3}\alpha \left( {\alpha + 1} \right)\left( {{\alpha ^2} + \alpha + 1} \right); \end{array} $

f1:=2α+x+1;

>res1:= collect(normal(f1+f4+3*f5), [x], factor);

$ \begin{array}{l} {\rm{resl}}: = 3{x^5} + \left( {4 + 6\alpha } \right){x^4} + \left( {6{\alpha ^2} + 8\alpha + 1} \right){x^3} + 4\alpha {\left( {\alpha + 1} \right)^2}{x^2} + \\ \;\;\;\;\;\;\;\;\;\left( {2{\alpha ^4} + \frac{{16}}{3}{\alpha ^3} + 6{\alpha ^2} + \frac{8}{3}\alpha + 1} \right)x + \\ \;\;\;\;\;\;\;\;\;1 + \frac{4}{5}{\alpha ^5} + \frac{8}{3}{\alpha ^4} + \frac{{16}}{3}{\alpha ^3} + \frac{{16}}{3}{\alpha ^2} + \frac{{58}}{{15}}\alpha . \end{array} $

载荷条件(ⅱ)

>alpha:= -1/2-1/10*′I;

$ \alpha :\; = - \frac{1}{2} = \frac{1}{{10}}I $

>res2:= SIEslover (1, x, alpha);

$ {\rm{res}}2:\; = - \frac{1}{5}I + x $

>for i to 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 + I k2 =%+15.10f, +I%+15.10f\n", i, k1, k2)

end if

od:

i= 5, … k1 + I k2 = +3.744 615 986 8, +I -1.211 230 446 2

i= 10, … k1 + I k2 = +3.865 123 759 7, +I -1.268 098 985 6

i= 15, … k1 + I k2 = +3.865 755 781 5, +I -1.268 449 900 5

i= 20, … k1 + I k2 = +3.865 765 299 0, +I -1.268 455 554 0

i= 25, … k1 + I k2 = +3.865 765 306 9, +I -1.268 455 559 0

i= 30, … k1 + I k2 = +3.865 765 306 9, +I -1.268 455 559 0

i= 35, … k1 + I k2 = +3.865 765 306 9, +I -1.268 455 559 0

参考文献
[1] GAKHOV F D. Boundary Value Problems[M]. New York: Dover, 1990.
[2] MUSKHELISHVILI N I. Singular integral equations:Boundary problems of function theory and their application to mathematical physics[J]. P Noordhoff N, 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]. Journal of Applied Mechanics, 1998, 65(2): 548. DOI:10.1115/1.2789094
[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]. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(2): 207–214. DOI:10.6052/0459-1879-2008-2-2007-123
[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]. Chinese Journal of Solid Mechanics, 2010, 31(4): 385–391.
[6] 周薇. 在接触力学中的奇异积分方程的高精度数值解法[D]. 成都: 电子科技大学, 2011.
ZHOU W.High-Accuracy Numerical Solution for Singular Integration Equations in Contact Mechanics[D]. Chengdu:University of Electronic Science and Technology of China, 2011.
[7] ERDOGAN F, GUPTA G D, COOK T. Numerical Solution of Singular Integral Equations[J]. 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]. Quarterly of Applied Mathematics, 1975, 33(3): 225–232. DOI:10.1090/qam/1975-33-03
[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]. Quarterly of Applied Mathematics, 1977, 35(1): 173–183. DOI:10.1090/qam/1977-35-01
[10] ERDOGAN F, GUPTA G. On the numerical solution of singular integral equations[J]. Quarterly of Applied Mathematics, 1972, 29(4): 525–534. DOI:10.1090/qam/1972-29-04
[11] 李星. 积分方程[M]. 北京: 科学出版社, 2008.
LI X. Integral Equation[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]. International Journal of Engineering Science, 1981, 19(9): 1293–1298. DOI:10.1016/0020-7225(81)90148-8
[13] MILLER G R, KEER L M. A numerical technique for the solution of singular integral equations of the second kind[J]. Quarterly of Applied Mathematics, 1985, 42(4): 455–465. DOI:10.1090/qam/1985-42-04
[14] JIN X. Analysis of Some Two Dimensional Problems Containing Cracks and Holes[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]. Engineering Fracture Mechanics, 2008, 75(5): 1005–1014. DOI:10.1016/j.engfracmech.2007.04.024
[16] 周跃亭, 李星. 具周期裂纹的半平面周期接触问题的奇异积分方程数值解法[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]. Chinese Journal of Solid Mechanics, 2005, 26(2): 167–174.
[17] KIM P, LEE S. A piecewise linear quadrature of Cauchy singular integrals[J]. Journal of Computational and Applied Mathematics, 1998, 95(1/2): 101–115.
[18] GERASOULIS A. Piecewise-polynomial quadratures for Cauchy singular integrals[J]. Siam Journal on Numerical Analysis, 1986, 23(4): 891–902. DOI:10.1137/0723057
[19] GERASOULIS A. The use of piecewise quadratic polynomials for the solution of singular integral equations of Cauchy type[J]. Computers & Mathematics with Applications, 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]. International Journal of Fracture, 1994, 66(2): 139–154. DOI:10.1007/BF00020079
[21] RABINOWITZ P. Convergence results for piecewise linear quadratures for Cauchy principal value integrals[J]. Mathematics of Computation, 1988, 51(184): 741–747. DOI:10.1090/S0025-5718-1988-0958639-3
[22] IOAKIMIDIS N I. On the numerical evaluation of derivatives of Cauchy principal value integrals[J]. Computing, 1981, 27(1): 81–88. DOI:10.1007/BF02243440
[23] ERDELYI A, MAGNUS W, OBERHETTINGER F, et al. Tables of Integral Transforms:Vol Ⅱ[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]. Zeitschrift Für Angewandte Mathematik und Physik, 1977, 28(6): 1085–1098. DOI:10.1007/BF01601675