裂紋問題超奇異積分方程法的程序?qū)崿F(xiàn)

杜云海，王 志，候建華

(鄭州大學(xué)力學(xué)與工程科學(xué)學(xué)院，河南鄭州450001)

0 引言

超奇異積分方程方法近些年來引起許多斷裂力學(xué)工作者的關(guān)注并開展深入地研究［1-4］，在研究獲得描述裂紋問題的超奇異積分方程組的基礎(chǔ)上，建立計(jì)算應(yīng)力強(qiáng)度因子的數(shù)值方法.由于利用了奇異積分的有限部積分，有效地解決了對奇異積分方程組求解的困難，使數(shù)值計(jì)算容易實(shí)現(xiàn)，計(jì)算的收斂速度也比較快捷，便于工程應(yīng)用.由于目前尚無該方法在編程方面的相關(guān)報(bào)道，筆者根據(jù)近年來在相關(guān)研究［5-6］中對于數(shù)值方法的編程實(shí)現(xiàn)所積累的經(jīng)驗(yàn)，就如何提高計(jì)算速度和輸出滿足精度要求的計(jì)算結(jié)果作以簡要論述.

1 求解問題的線性方程組

作者的最新研究發(fā)現(xiàn)，雙材料平面內(nèi)一條任意曲線裂紋(長度為2a)問題均可以用一個(gè)以裂紋岸位移間斷為基本未知量的超奇異積分方程組來描述，通過適當(dāng)?shù)淖兞颗c函數(shù)替換，對方程組進(jìn)行歸一化處理，并將未知的裂紋間斷函數(shù)用第二類契比謝夫多項(xiàng)式表示，取該多項(xiàng)式在區(qū)間(-1，1)內(nèi)的N個(gè)零點(diǎn)作為配置點(diǎn)，最終可把問題歸結(jié)為對一個(gè)有2N個(gè)方程、2N個(gè)未知量的線性方程組的求解.這個(gè)線性方程組可用張量形式表示為

其基本未知量變成未知的第二類契比謝夫多項(xiàng)式系數(shù)akj.方程組(1)的系數(shù)Gikj(τm)與裂紋形狀有關(guān)，以一條圓弧裂紋(圖1)為例，方程組(1)的系數(shù)為

圖1 雙材料平面中的圓弧裂紋Fig.1 A arc crack in plane bimaterial

這里Hik(ξ，τ)為正常積分核函數(shù)，而Gik為奇異積分的有限部積分，具體形式與結(jié)果為

其中:Tj+1(τ)和Uj(τ)分別為第一類和第二類契比謝夫多項(xiàng)式，λ為圓弧半徑比.方程組(1)中的gi(τm)決定于裂紋岸的分布載荷qi.ζ，τ表示裂紋弧長坐標(biāo).第二類契比謝夫多項(xiàng)式在區(qū)間(-1，1)上的N個(gè)零點(diǎn)為

2 應(yīng)力強(qiáng)度因子計(jì)算公式

由斷裂力學(xué)，裂紋尖端 A，B兩點(diǎn)的應(yīng)力強(qiáng)度因子可表達(dá)為

這里Fk(±1)由第二類契比謝夫多項(xiàng)式在區(qū)間端點(diǎn)的極限值來確定，可利用方程(1)的解進(jìn)行計(jì)算，即

3 程序設(shè)計(jì)問題及其例程

在通過編程實(shí)現(xiàn)數(shù)值方法時(shí)，程序須按上述方法建立的關(guān)系來編寫，這里不必贅述.但編程時(shí)須注意如何控制計(jì)算過程使之盡可能快速輸出滿足精度要求的應(yīng)力強(qiáng)度因子.這樣就必須把握好數(shù)值積分方法、方程組求解方法、以及適當(dāng)配置方程個(gè)數(shù)以滿足精度要求等幾個(gè)環(huán)節(jié).下面就這些關(guān)鍵問題逐一進(jìn)行闡述.

3.1 關(guān)于數(shù)值積分的方法

數(shù)值積分是超奇異積分方程方法計(jì)算方程組系數(shù)的必需步驟，由于積分函數(shù)涉及契比謝夫多項(xiàng)式和復(fù)雜的積分核函數(shù)，實(shí)踐表明，裂紋的位置越靠近雙材料界面，采用一般數(shù)值積分法的收斂速度就越慢，如梯形積分法等會(huì)達(dá)到讓人無法忍耐的緩慢程度.高斯積分法是精度相對較高而且較快的方法，但它只能配置固定的積分點(diǎn)，具有相同的計(jì)算時(shí)間.由于裂紋位置不同時(shí)，數(shù)值積分的收斂速度是不同的，這樣就不利于控制在不同情況下獲得統(tǒng)一的計(jì)算精度，或者不利于根據(jù)精度要求適當(dāng)伸縮計(jì)算時(shí)間.相對而言，逐次分半加速法是一種比較適用的積分方法，既可以避免重復(fù)計(jì)算，又可根據(jù)設(shè)定的積分精度控制計(jì)算過程.假定設(shè)置精度為δ1(如取為1×10-4)，則當(dāng)前、后兩次計(jì)算結(jié)果相對誤差小于δ1時(shí)便可返回計(jì)算結(jié)果，否則繼續(xù)計(jì)算.

3.2 關(guān)于方程組求解的方法

對于線性方程組的求解，可以采用不同的方法，如主元素消去法、賽德爾迭代法等都可以快速求解.這里推薦比較簡單的主元素消去法，它可以有效抑制系數(shù)矩陣變換過程所產(chǎn)生的計(jì)算誤差.

3.3 關(guān)于應(yīng)力強(qiáng)度因子精度的控制

在編寫計(jì)算程序時(shí)，須注意計(jì)算結(jié)果應(yīng)能輸出滿足精度要求的應(yīng)力強(qiáng)度因子.這個(gè)問題可通過適當(dāng)選擇配置點(diǎn)數(shù)(N階第二類契比謝夫多項(xiàng)式的N個(gè)零點(diǎn))來解決.假如應(yīng)力強(qiáng)度因子的計(jì)算精度要求為δ2=1×10-3，在裂紋位置遠(yuǎn)離雙材料界面時(shí)，配置20幾個(gè)點(diǎn)就足以滿足精度要求;但在裂紋靠近界面時(shí)，就需要較多的配置點(diǎn)才能滿足精度要求.為了得到剛好滿足精度的計(jì)算結(jié)果，在編程時(shí)可以先設(shè)定一個(gè)配置點(diǎn)數(shù)的初值N0，通過計(jì)算求得應(yīng)力強(qiáng)度因子后先用變量存儲起來，然后再將配置點(diǎn)數(shù)增加一個(gè)增量ΔN，用點(diǎn)數(shù)N1=N0+ΔN重新計(jì)算，計(jì)算結(jié)果與上次結(jié)果比較，如滿足精度要求則停止計(jì)算而輸處結(jié)果，如不滿足要求則再用點(diǎn)數(shù)N2=N1+ΔN進(jìn)行計(jì)算，如結(jié)果與點(diǎn)數(shù)為N1的結(jié)果相比滿足精度要求，則輸出結(jié)果，否則再增加點(diǎn)數(shù)ΔN進(jìn)行計(jì)算，如此類推，直至當(dāng)下點(diǎn)數(shù)的計(jì)算結(jié)果與上次結(jié)果相比滿足精度要求為止，這樣最終就得到穩(wěn)定的應(yīng)力強(qiáng)度因子值.

采用上述方法，建立和求解方程組的次數(shù)至少需要兩次，一般3到4次計(jì)算就可獲得理想的計(jì)算結(jié)果.當(dāng)然，循環(huán)次數(shù)肯定與最初設(shè)置的配置點(diǎn)數(shù)N0有關(guān).若初始點(diǎn)數(shù)小，數(shù)值積分速度快，但在控制應(yīng)力強(qiáng)度因子精度的情況下可能會(huì)造成更多的循環(huán)計(jì)算次數(shù);若初始點(diǎn)數(shù)大，數(shù)值積分速度慢，但在控制精度情況下的循環(huán)計(jì)算次數(shù)反而會(huì)少.故選擇適當(dāng)?shù)某跏键c(diǎn)數(shù)可有效減少計(jì)算所需的總機(jī)時(shí).以下是圓弧裂紋問題采用這種算法的Fortran9.0主程序，其中配置點(diǎn)數(shù)的初始值選為N0=n1=38.每次循環(huán)計(jì)算增加點(diǎn)數(shù)選為ΔN=4，應(yīng)力強(qiáng)度因子控制精度為1×10-3.

program sclyhlw

double precision ka(2)，kb(2)，k1a(2)，

*k1b(2)，tb，rb，xob，xta0，xta1，gs1(2，0:100)，

*gs(2，2，0:100，0:100)，fy(2) ，a(2，0:100)，

*g(2，0:100，0:202)，s(0:100)，g1(0:201，0:202)，

*x(202)，rs(0:100，0:100)，fa(2)，fb(2)，

*g11(202，203)，ss1，dta，dtb

character*20 jsjg

write(* ，*)＇輸入結(jié)果文件名:'

read(*，*)jsjg

open(6，file=jsjg，status='new')

write(* ，*)'輸入?yún)?shù):tb，rb，xob，xta0，fy(1)

*，fy(2)'

read(* ，*)tb，rb，xob，xta0，fy(1)，fy(2)

write(6，*)tb，rb，xob，xta0，fy(1)，fy(2)計(jì)算方程組系數(shù)

n1=38

xta1=xta0*3.1415927D-00/180.00

do 10 k=1，2

ka(k)=0.0

kb(k)=0.0

10 continue

20 do 100 i=1，2

do 90 k=1，2

do 80 j=0，n1

s(j)=dcos((j+1)*3.1415927d-00/

* (n1+2.0))

gs1(i，j)=-fy(i)*3.1415927

do 70 n=0，n1

if(i.eq.1.and.k.eq.1)then

if(s(j).gt.1.0-5.0d-05)then

rs(j，n)=-3.1415927*(n+1)＊＊2

else if(s(j).lt.-1.0+5.0d-05)then

rs(j，n)=-3.1415927*(n+1)＊＊2*(-1)＊＊n

else

rs(j，n)=-3.1415927*(n+1)* dsin((n+1)*

＊dacos(s(j)))/dsin(dacos(s(j)))

end if

else if(i.eq.1.and.k.eq.2)then

if(s(j).gt.1.0-5.0d-05)then

rs(j，n)=3.1415927/rb

else if(s(j).lt.-1.0+5.0d-05)then

rs(j，n)=3.1415927*(-1)＊＊(n+1)/rb

else

rs(j，n)=3.1415927*dcos((n+1)*

＊dacos(s(j)))/rb

end if

else if(i.eq.2.and.k.eq.1)then

if(s(j).gt.1.0-5.0d-05)then

rs(j，n)=-3.1415927/rb

else if(s(j).lt.-1.0+5.0d-05)then

rs(j，n)=-3.1415927*(-1)＊＊(n+1)/rb

else

rs(j，n)=-3.1415927*dcos((n+1)*

＊dacos(s(j)))/rb

end if

else

if(s(j).gt.1.0-5.0d-05)then

rs(j，n)=-3.1415927*(n+1)＊＊2

else if(s(j).lt.-1.0+5.0d-05)then

rs(j，n)=-3.1415927*(n+1)＊＊2*(-1)＊＊n

else

rs(j，n)=-3.1415927*(n+1)* dsin((n+1)*

＊dacos(s(j)))/dsin(dacos(s(j)))

end if

end if

調(diào)用數(shù)值積分子程序szjf，并將結(jié)果存于四維臨時(shí)數(shù)組gs

call szjf(i，k，n，s(j)，tb，ss1，rb，xob，xta1)

gs(i，k，j，n)=rs(j，n)+ss1

write(* ，'(5x，3hg( ，i1，1h，，i1，1h，，i2，1h，，i2，

*4h)= ，f18.8)')i，k，j，n，gs(i，k，j，n)

70 continue

80 continue

90 continue

100 continue

整理方程組，將系數(shù)賦給一個(gè)二維數(shù)組g11

write(*，*)'solve equtions:'

do 200 i=1，2

do 190 j=0，n1

do 180 n=0，n1

g(i，j，n)=gs(i，1，j，n)

g(i，j，n1+1+n)=gs(i，2，j，n)

180 continue

g(i，j，2*n1+2)=gs1(i，j)

190 continue

200 continue

do 270 j=0，n1

do 260 n=0，2*n1+2

g1(j，n)=g(1，j，n)

g1(j+n1+1，n)=g(2，j，n)

260 continue

270 continue

do 280 i=1，2*n1+2

do 275 j=1，2*n1+3

g11(i，j)=g1(i-1，j-1)

275 continue

280 continue

write(* ，285)((g11(i，j)，j=1，2*n1+3)，i=1，2*n1+2)

285 format(1x，11f7.4)

調(diào)用主元素法子程序zyss求解方程組g11，并將結(jié)果賦給二維數(shù)組a

call zyss(g11，x，2*n1+2，2*n1+3)

write(*，*)'outo zyss'

do 290 n=0，n1

a(1，n)=x(n+1)

a(2，n)=x(n+n1+2)

write(* ，'(5x，4ha(1，，i2，2h)= ，f12.6，10x，

*4ha(2，，i2，2h)= ，F(xiàn)12.6)')n，a(1，n)，n，a(2，n)

290 continue

利用數(shù)組a計(jì)算應(yīng)力強(qiáng)度因子，并與上次結(jié)果比較，在誤差dta、dtb不滿足精度要求的條件下，增加配置點(diǎn)數(shù)，暫存當(dāng)前結(jié)果，并返回到2 0號語句重新計(jì)算

do 600 k=1，2 fa(k)=0.0 fb(k)=0.0

do 500 n=0，n1

fa(k)=fa(k)+a(k，n)*(n+1)*(-1)＊＊n

fb(k)=fb(k)+a(k，n)*(n+1)500 continue

k1a(k)=fa(k)

k1b(k)=fb(k)

dta=dabs(k1a(k)-ka(k))

dtb=dabs(k1b(k)-kb(k))

if(dta.gt.1.0d-03.or.dtb.gt.1.0d-03)then

n1=n1+4

if(n1.gt.100)then

write(* ，*)'n1＞100，stop!'

stop

end if

ka(k)=k1a(k)

kb(k)=k1b(k)

goto 20

end if

輸出滿足精度的強(qiáng)度因子當(dāng)前值k1a、k1b

write(6，'(1x，1x，3htb= ，f5.2，1x，3hrb= ，f8.2，1x，4h

*xob=，f8.2 ，1x，5hxta0=，f5.2，1x，6hfy(1)=，f5.2，

*1x，6hfy(2)= ，f5.2，1x，2hk= ，i1，1x，3hka= ，f7.3，1x，

*3hkb=，f7.3)')tb，rb，xob，xta0，fy(1)，fy(2)，k，

*k1a(k)，k1b(k)600 continue

end

該程序?qū)⒂?jì)算結(jié)果寫入開始運(yùn)行時(shí)文本變量“jsjg”要求賦值的數(shù)據(jù)文件.

3.4 其他問題

在程序中還需要正常積分核函數(shù)，這需要根據(jù)正常積分核函數(shù)的具體形式編制函數(shù)子程序，如上述例程中的“sjf”，由于積分核函數(shù)的關(guān)系復(fù)雜，這里不再給出相應(yīng)的例程，具體函數(shù)的形式可關(guān)注作者近期的相關(guān)文獻(xiàn).

3.5 算例

啟動(dòng)程序后，按要求須首先輸入計(jì)算結(jié)果的存儲文件名(如“JG”)，然后再輸入彈性模量比、圓弧半徑比、圓弧中點(diǎn)坐標(biāo)比、圓弧中點(diǎn)半徑與x軸的夾角、以及裂紋岸載荷等所需參數(shù)，輸入后即可自動(dòng)實(shí)施計(jì)算，并將輸入數(shù)據(jù)與運(yùn)行結(jié)果寫入開始定義的存儲文件“JG”.用Windows系統(tǒng)附件的記事本或?qū)懽职骞ぞ呖纱蜷_瀏覽結(jié)果文件“JG”，顯示內(nèi)容如下所示:

輸入?yún)?shù):tb，rb，xob，xta0，fy(1)，fy(2)

1.000000000000000 5.000000000000000

3.000000000000000 0.000000000000000

1.000000000000000 0.000000000000000

tb=1.00 rb=5.00 xob=3.00 xta0=.00 fy(1)=1.00 fy(2)=.00 k=1 ka=.980 kb=.980 tb=1.00 rb=5.00 xob=3.00 xta0=.00 fy(1)=1.00 fy(2)=.00 k=2 ka=.098 kb=-.099

這一計(jì)算過程所用時(shí)間一般不超過5 min.

4 結(jié)論

基于雙材料平面中裂紋岸受分布載荷作用的圓弧裂紋問題的歸一化超奇異積分方程組及其數(shù)值方法，重點(diǎn)討論了用Fortran9.0編程工具進(jìn)行程序設(shè)計(jì)中需要考慮的數(shù)值積分方法、方程組求解方法、以及控制應(yīng)力強(qiáng)度因子精度的循環(huán)算法等相關(guān)問題，并給出部分例程編碼.實(shí)際例程的運(yùn)行結(jié)果表明，所采用的算法可以獲得較高運(yùn)算速度和滿足精度要求的應(yīng)力強(qiáng)度因子.

［1］ISIDA M，NOGUCHI H.An arbitrary array of cracks in bonded semi-infinite bodies under in-plane loads［J］.Trans.JSME，1983，49:36-45.

［2］DONG C Y，LEE K Y.A new integral equation formulation of two-dimensional inclusion-crack problems［J］.International Journal of Solid and Structure，2005，42:5010-5020.

［3］KAYA A C，ERDOYAN F.On the solution of integral equations with strongly singular kenels［J］.Quarterly of Applied Mathematics，1987(A5):105-122

［4］MARTIN P A，IZZO F J R.On boundary integral equations for crack problems［J］.Proc.R.Soc.Lond，1988，421:341-355.

［5］杜云海，樂金朝.雙材料平面斜裂紋問題超奇異積分方程方法［J］.機(jī)械強(qiáng)度，2004，26(3):326-331.

［6］杜云海，呂存景，董棟，等.雙材平面中環(huán)形界面環(huán)向裂紋問題的超奇異積分方程方法［J］.機(jī)械強(qiáng)度，2006，28(5):733-738.