郝文峰， 陳新文， 鄧立偉， 王 翔， 姚學鋒

(1.北京航空材料研究院航空材料檢測與評價北京市重點實驗室，北京100095;2.清華大學工程力學系，北京100084)

數(shù)字圖像相關方法(digital image correlation，DIC)是一種非條紋光測方法，它是現(xiàn)代先進的光電子技術、數(shù)字圖像處理技術與計算機技術相結合的產(chǎn)物，它直接對物體變形前、后的兩幅散斑圖像進行相關運算處理，以提取散斑圖中所攜帶的變形信息。散斑是變形信息的載體，可以通過人工制斑獲得或者就直接以試件表面的自然紋理作為標記。

在二維數(shù)字圖像相關方法的應用方面，國內(nèi)外都開展了大量的工作。Corr等［1］利用數(shù)字圖像相關方法研究了纖維復合材料與混凝土界面的粘接強度，得到較為精確的表面變形場。Fedele等［2］利用數(shù)字圖像相關方法研究GLARE層板層間粘接層材料參數(shù)，將材料參數(shù)用于有限元模擬，并將模擬結果與試驗進行對比。Rjafiallah等［3］提出一種三相模型來分析玻璃微珠增強聚合物基復合材料界面相影響，并用數(shù)字圖像相關方法進行試驗驗證。Desai等［4］利用數(shù)字圖像相關方法測得的位移場計算得到雙材料界面裂紋應力強度因子，并將試驗結果與理論和有限元結果進行對比。Carloni等［5］利用數(shù)字圖像相關方法研究疲勞載荷作用下復合材料-混凝土界面裂紋擴展過程。姚學鋒及其合作者對數(shù)字圖像相關方法開展系統(tǒng)研究［6］，并且將數(shù)字圖像相關方法成功應用于復合材料壓力容器測量［7］、人體頸椎生物力學［8］、有機玻璃黏接力學［9］、圓管殘余應力測試［9］和編織復合材料力學性能表征［10～12］等一系列新型材料和新型問題的力學測試。

本工作將數(shù)字圖像相關方法用于芳綸纖維復合材料I型斷裂問題研究，通過最小二乘擬合法提取裂紋尖端應力強度因子，分析數(shù)字圖像相關方法中的子區(qū)域和步長大小對計算結果的影響。

1 基于數(shù)字圖像相關方法的應力強度因子測量

1.1 二維數(shù)字圖像相關方法的基本原理

對比變形前、后散斑圖像的灰度特征是數(shù)字圖像相關方法的基本原理。被測試樣表面需要有足夠多的可識別標志點，一般通過在試樣表面噴黑白相間的漆斑作為散斑標志點，并且假設這些散斑點的變形就是試樣本身的變形。

進行圖像前處理以及選擇合理的匹配原則來定義相關度是數(shù)字圖像相關方法算法的關鍵。合理的對應關系，通常是基于如下兩個前提條件［13］:(1)物體表面上的同一個點在變形前后圖像上的灰度保持不變;(2)隨機分布的散斑使得圖像上的任一個包含有足夠多的像素點的子集在灰度分布上具有唯一性。

在二維數(shù)字圖像相關方法中通常將待計算區(qū)域劃分為許多小的方形小區(qū)域，一般稱為相關計算子區(qū)域。子區(qū)域中像素點的數(shù)量決定了子區(qū)域的大小，對于一個M×M的方形區(qū)域，表示其在橫向和縱向各包含M個像素點。相鄰兩個子區(qū)域之間的距離稱為步長，步長大小決定了相關計算所得位移場的分辨率。

目前常用的各種二維數(shù)字圖像相關方法的算法，基本上都是以位移為參數(shù)，建立變形前后圖像上的像素點間的映射關系，如圖1所示，變形前圖像上的像素點(x，y)的位移為(u，v)，對應于變形后圖像上的點(x'，y')。

圖1 數(shù)字圖像相關方法原理Fig.1 Principle of digital image correlation method

建立映射關系后，再選擇一個相關公式，計算變形前后圖像子區(qū)間的相關系數(shù)，當相關系數(shù)取極值時，認為子區(qū)間的匹配是最佳的，然后可以得到位移等變形量。目前，相關運算方法歸納起來有以下幾個主要途徑［13］:(1)雙參數(shù)法;(2)粗細搜索法;(3) Newton-Raphson方法;(4)十字搜索法;(5)爬山搜索法;(6)相關系數(shù)梯度法;(7)遺傳算法;(8)分形算法。

1.2 最小二乘擬合提取應力強度因子

根據(jù)文獻［19］，I型裂紋尖端x和y方向位移場可以表示為:

式中，r和θ是裂尖任意點的極坐標，κ為彈性常數(shù)，平面應力狀態(tài)下為，平面應變狀態(tài)下為3-4ν，ν為泊松比。由于數(shù)字圖像相關數(shù)據(jù)由表面位移計算所得，故需考慮平面應力狀態(tài)。μ為剪切模量。對于靠近裂紋尖端區(qū)域式(3)、式(4)第一項起主導作用，而對于遠離裂尖區(qū)域必須考慮高階項的影響才能精確得到裂尖位移場。

對于式(3)、式(4)，在二維數(shù)字圖像相關方法試驗中，必須考慮由于剛體位移和平面內(nèi)轉(zhuǎn)動引起的附加項。因此，對于由數(shù)字圖像相關方法試驗數(shù)據(jù)提取的位移場公式應修訂為:

式中，上標j表示M個數(shù)據(jù)點中的第j個數(shù)據(jù)點，其坐標為(xj，yj)。uj和vj為點(xj，yj)在x方向和y方向的位移值。u*，v*分別為x方向和y方向的剛體位移，α為平面內(nèi)轉(zhuǎn)角。裂尖極坐標和直角坐標的轉(zhuǎn)換關系為

其中(x0，y0)為裂尖處坐標值。由斷裂力學獲得式(9):

本工作考慮拉伸載荷作用下的正交異性單邊裂紋芳綸纖維復合材料板，由最小二乘擬合法來確定式(5)、式(6)中的未知系數(shù)。因此，其表達式為［15］:

其中N為式(5)的項數(shù)，三個未知數(shù)為兩個方向剛體位移和轉(zhuǎn)動角。未知參數(shù){χ}，系數(shù)矩陣［R］以及向量{Λ}分別為:

式中［A］為M×N矩陣，包含了u，v表達式中的aN。

本研究拉伸載荷作用下的單邊I型裂紋裂尖應力強度因子為:

其中P為拉伸載荷，B為試驗件厚度，W為試驗件寬度，a為裂紋長度。函數(shù)是與試驗件尺寸相關的函數(shù)，可按文獻［21］中的曲線取值，對于本研究試驗件幾何尺寸=1.25。

2 試驗件及試驗裝置

試驗件為3233/023A芳綸纖維復合材料，如圖2所示。試驗件尺寸如圖2a所示，其中試驗件寬度W=20mm，長度H=200mm，厚度B=3mm，裂紋長度a=5mm。在試驗件正面首先噴上一薄層白漆作為底色，然后再噴上黑漆斑作為數(shù)字圖像相關的特征斑點，如圖2b所示。數(shù)字圖像相關測試平臺及試驗裝置如圖3所示，主要包括試驗加載系統(tǒng)、冷光源、CCD相機以及裝有數(shù)字圖像相關算法軟件的計算機。試驗加載速率為0.5mm/min，試驗載荷為0至5000N間隔500N采集一幅圖像，CCD相機和冷光源放置于隔震臺上以減少試驗噪聲干擾。試驗所用CCD相機分辨率為1280×1024，圖像中1mm長度對應于46.5個像素點。為了研究數(shù)字圖像相關方法中子區(qū)域和步長大小對計算結果的影響，本計算中選取五種子區(qū)域，其大小分別為21×21，41 ×41，61×61，81×81，101×101，五種步長取子區(qū)域長度的0.1至0.5倍，間隔0.1。

圖2 試驗件 (a)試驗件示意圖;(b)噴斑后試驗件Fig.2 Specimen (a)schematic diagram of specimen; (b)specimens after speckled

圖3 試驗裝置Fig.3 Experimental set-up

3 試驗結果及分析

試驗件變形前后的散斑場如圖4所示，其中圖4a為變形前散斑圖，圖4b為載荷1000N時的散斑圖，計算區(qū)域以及裂紋位置見圖4a，其中黑線位置為裂紋，紅色方框內(nèi)為計算區(qū)域。由圖4a和圖4b無法簡單的辨別出差異，但通過數(shù)字圖像相關方法計算后可以得到試件在x方向和y方向的位移場，如圖5a和圖5b所示。

圖5為數(shù)字圖像相關方法計算所得全場位移，其中包括剛體位移和面內(nèi)轉(zhuǎn)動。由圖可見，裂紋尖端全場位移分布具有較好的對稱性。對圖中裂尖部分取半圓按最小二乘擬合方法計算裂紋尖端應力強度因子，對于不同擬合項數(shù)和載荷大小計算所得的應力強度因子結果如圖6所示。圖6a為載荷2000N時不同擬合項次下的應力強度因子值，由圖6a可見，在載荷為2000N時，隨著最小二乘擬合項數(shù)的增加，擬合所得應力強度因子結果與式(14)計算結果的誤差越來越小。當只取一項忽略高階項時，擬合結果與計算結果誤差為22.8%，而當取N=5時所得應力強度因子誤差減小為4.1%。為減小計算量并保證一定的精度，其他計算均取N=5。圖6b為應力強度因子與施加載荷曲線，其中由式(14)計算所得應力強度因子隨施加載荷的增加而線性增大。數(shù)字圖像相關方法計算位移場擬合所得應力強度因子與式(14)計算結果偏差較小。

圖4 DIC計算散斑圖 (a)變形前散斑圖(0N);(b)變形后散斑圖(1000N)Fig.4 Speckle field using calculation (a)before deformation(0N);(b)after deformation(1000N)

圖5 DIC計算位移場 (a)U場，1000N;(b)V場，1000N;(c)U場，2000N; (d)V場，2000N;(e)U場，5000N;(f)V場，5000NFig.5 Displacement field calculated from DIC (a)U field，1000N;(b)V field，1000N; (c)U field，2000N;(d)V field，2000N;(e)U field，5000N;(f)V field，5000N

圖7為載荷2000N時，數(shù)字圖像相關計算時子區(qū)域大小和步長大小對最小二乘擬合所得應力強度因子的影響，由圖7a可見子區(qū)域大小為41時，擬合結果與式(14)計算結果最為接近，這主要是因為如果區(qū)域太小，包含的信息量不夠，難以準確匹配二幅散斑圖的相關區(qū)域;區(qū)域太大，則會造成比較嚴重的平均效果，同樣亦難以準確匹配二幅散斑圖的相關區(qū)域，這兩種情況都會降低測量精度。子區(qū)域大小為101時誤差為6.0%，子區(qū)域大小為21時誤差為4.7%，均大于子區(qū)域大小為41時的4.1%。由圖7b可見，子區(qū)域大小為41時，數(shù)字圖像相關計算的步長越小，其結果越接近于式(14)計算結果，五種步長大小計算所得結果偏差均在5%以內(nèi)。數(shù)字圖像相關計算時步長大小對結果的影響小于子區(qū)域大小的影響。

圖6 DIC應力強度因子 (a)擬合項數(shù)影響;(b)應力強度因子-載荷曲線Fig.6 Stress intensity factors (a)effect of least squares fitting terms;(b)stress intensity factors-load curve

圖7 DIC計算參數(shù)影響 (a)子區(qū)域大小影響;(b)步長大小影響Fig.7 Effect of DIC parameters (a)effect of subset size;(b)effect of subset spacing

4 結論

(1)通過最小二乘擬合法建立了考慮剛體位移的數(shù)字圖像相關方法計算所得位移場與I型裂紋尖端應力強度因子之間的關系。

(2)試驗研究了含單邊I型裂紋芳綸纖維復合材料在拉伸載荷作用下裂紋尖端位移場演化規(guī)律。提取裂紋尖端應力強度因子，并將試驗結果與理論公式計算結果進行對比，結果表明試驗所得應力強度因子與理論結果偏差在10%以內(nèi)。

(3)分析最小二乘擬合項數(shù)、數(shù)字圖像相關計算子區(qū)域大小和步長大小對試驗提取應力強度因子的影響，結果表明增加最小二乘擬合項數(shù)可以提高應力強度因子測試精度，而數(shù)字圖像相關計算子區(qū)域大小和步長大小的影響較為復雜，需要通過經(jīng)驗或者試算確定合理的區(qū)域和步長大小來提高計算精度。

［1］CORR D，ACCARDI M，GRAHAM-BRADY L，et al. Digital image correlation analysis of interfacial debonding properties and fracture behavior in concrete［J］.Engineering Fracture Mechanics，2007，74(1/2):109-121.

［2］FEDELER，RAKA B，HILD F，etal.Identification of adhesive properties in GLARE assemblies using digital image correlation［J］.JMech Phys Solids，2009，57(7):1003 -1016.

［3］RJAFIALLAH S，GUESSASMA S.Three-phasemodel and digital image correlation to assess the interphase effect on the elasticity of carbohydrate polymer-based composites reinforced with glass-silica beads［J］.Carbohyd Polym，2011，83(1):246-256.

［4］DESAIC K，BASU S，PARAMESWARAN V.Determination of complex stress intensity factor for a crack in a bimaterial interface using digital image correlation［J］.Opt Laser Eng，2012，50(10):1423-1430.

［5］CARLONIC，SUBRAMANIAM K V.Investigation of subcritical fatigue crack growth in FRP/concrete cohesive interface using digital image analysis［J］.Composites Part B:Engineering，2013，51:35-43.

［6］MENG L B，JIN G C，YAO X F.Application of iteration and finite element smoothing technique for displacement and strain measurement of digital speckle correlation［J］. Optics and Lasers in Engineering，2007，45(1):57-63.

［7］YAO X F，MENG L B，JIN JC，etal.Full-field deformation measurement of fiber composite pressure vessel using digital speckle correlation method［J］.Polymer Testing，2005，24(2):245-251

［8］MA Y J，XIONG T C，YAO X F.Experimental investigation of interface curing stresses between PMMA［J］.Theoretical＆ Applied Mechanics Letters，2011，1(5): 051003.

［9］MA Y J，YAO X F，ZHANG D W.Axially symmetrical stresses measurement in the cylindrical tube using DIC with hole-drilling［J］.Optics and Lasers in Engineering，2015，66:174-180

［10］MA Y J，YAO X F，ZHENGQ S，etal.Carbon nanotube films change Poisson's ratios from negative to positive［J］. Applied Physics Letters，2010，97(6):061909.

［11］MA Y J，YAO X F，WANG D.Experimental investigation onmechanical properties of CNT film using digital speckle correlation method［J］.Optics and Lasers in Engineering，2012，50(11):1575-1581.

［12］HAOW F，GED Y，MA Y J，etal.Experimental investigation on deformation and strength of carbon/epoxy laminated curved beams［J］.Polymer Testing，2012，31(4): 520-526.

［13］孟利波.數(shù)字散斑相關方法的研究和應用［D］.北京:清華大學，2005. (MENG L B.Study and Application of Digital Speckle Correlation Method［D］.Beijing:Tsinghua University，2005.)

［14］YONEYAMA S，MORIMOTO Y，TAKASHIM.Automatic evaluation of mixed-mode stress intensity factors utilizing digital image correla-tion［J］.Strain，2006，42:21-29.

［15］ESHRAGHI I，YADEGARI DEHNAVIM R，SOLTANI N.Effect of subset parameters selection on the estimation ofmode-I stress intensity factor in a cracked PMMA specimen using digital image correlation［J］.Polymer Testing，2014，37:193-200.

［16］CRAVERO S，RUGGIERIC.Estimation procedure of J-resistance curves for SE(T)fracture specimens using unloading compliance［J］.Engineering Fracture Mechanics，2007，74:2735-2757.