裂紋問題超奇異積分方程法的程序實現

杜云海，王 志，候建華

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

0 引言

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

1 求解問題的線性方程組

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

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

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

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

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

2 應力強度因子計算公式

由斷裂力學，裂紋尖端 A，B兩點的應力強度因子可表達為

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

3 程序設計問題及其例程

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

3.1 關于數值積分的方法

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

3.2 關于方程組求解的方法

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

3.3 關于應力強度因子精度的控制

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

采用上述方法，建立和求解方程組的次數至少需要兩次，一般3到4次計算就可獲得理想的計算結果.當然，循環次數肯定與最初設置的配置點數N0有關.若初始點數小，數值積分速度快，但在控制應力強度因子精度的情況下可能會造成更多的循環計算次數;若初始點數大，數值積分速度慢，但在控制精度情況下的循環計算次數反而會少.故選擇適當的初始點數可有效減少計算所需的總機時.以下是圓弧裂紋問題采用這種算法的Fortran9.0主程序，其中配置點數的初始值選為N0=n1=38.每次循環計算增加點數選為ΔN=4，應力強度因子控制精度為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(* ，*)＇輸入結果文件名:'

read(*，*)jsjg

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

write(* ，*)'輸入參數: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)計算方程組系數

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

調用數值積分子程序szjf，并將結果存于四維臨時數組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

整理方程組，將系數賦給一個二維數組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)

調用主元素法子程序zyss求解方程組g11，并將結果賦給二維數組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)= ，F12.6)')n，a(1，n)，n，a(2，n)

290 continue

利用數組a計算應力強度因子，并與上次結果比較，在誤差dta、dtb不滿足精度要求的條件下，增加配置點數，暫存當前結果，并返回到2 0號語句重新計算

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

輸出滿足精度的強度因子當前值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

該程序將計算結果寫入開始運行時文本變量“jsjg”要求賦值的數據文件.

3.4 其他問題

在程序中還需要正常積分核函數，這需要根據正常積分核函數的具體形式編制函數子程序，如上述例程中的“sjf”，由于積分核函數的關系復雜，這里不再給出相應的例程，具體函數的形式可關注作者近期的相關文獻.

3.5 算例

啟動程序后，按要求須首先輸入計算結果的存儲文件名(如“JG”)，然后再輸入彈性模量比、圓弧半徑比、圓弧中點坐標比、圓弧中點半徑與x軸的夾角、以及裂紋岸載荷等所需參數，輸入后即可自動實施計算，并將輸入數據與運行結果寫入開始定義的存儲文件“JG”.用Windows系統附件的記事本或寫字板工具可打開瀏覽結果文件“JG”，顯示內容如下所示:

輸入參數: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

這一計算過程所用時間一般不超過5 min.

4 結論

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

［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］.機械強度，2004，26(3):326-331.

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