王 楠,李恩普,湯忠斌,李發(fā)東
(1.西北工業(yè)大學(xué) 理學(xué)院 陜西省光信息技術(shù)重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710129;2.西北工業(yè)大學(xué) 航空學(xué)院,陜西 西安 710072)
數(shù)字圖像相關(guān)(digital image correlation,DIC)方法是一種非接觸、全場變形測量技術(shù),因其具有設(shè)備簡單、對環(huán)境要求低、測量精度高等優(yōu)點(diǎn),已被廣泛用于材料的力學(xué)性能測試中[1,2]。DIC方法分二維(2-D)和三維(3-D)兩種:3-D DIC需要兩臺攝像機(jī),實(shí)驗(yàn)較為繁瑣;2-D DIC僅用一臺攝像機(jī),且不需要相機(jī)標(biāo)定。盡管2-D DIC精度略低于3-D DIC[3],但是通過修正系統(tǒng)誤差,仍可達(dá)到相當(dāng)高的精度,因此仍然具有相當(dāng)高的實(shí)用價值。
研究人員對2-D DIC的誤差做了大量的研究,且將誤差源分為硬件和軟件兩類[4]:(1)硬件方面主要包括散斑圖質(zhì)量[5]、鏡頭畸變[6]、離面位移[3]等;(2)軟件方面主要包括子區(qū)域大?。?]、相關(guān)函數(shù)[8]、亞像素插值[9]等。這些工作從原理上發(fā)展和完善了DIC方法,使其在實(shí)際應(yīng)用上又向前邁進(jìn)了一步。
力學(xué)實(shí)驗(yàn)中最常見是拉伸實(shí)驗(yàn),將2-D DIC方法應(yīng)用于拉伸實(shí)驗(yàn)已有很多報道[1,2],但卻沒有規(guī)范的實(shí)驗(yàn)方法,研究人員多是憑自身理解和經(jīng)驗(yàn)進(jìn)行實(shí)驗(yàn)。DIC作為一種圖像測量方法,由硬件及實(shí)驗(yàn)方法引入的系統(tǒng)誤差和偶然誤差對測量精度和可靠性都有很大的影響,因此為提高實(shí)驗(yàn)的可靠性,對實(shí)驗(yàn)誤差分析具有重要意義。
文中主要分析了實(shí)驗(yàn)條件對2-D DIC影響,對多晶銅試樣進(jìn)行拉伸實(shí)驗(yàn),以應(yīng)變片的測量結(jié)果為基準(zhǔn),將2-D DIC的測量結(jié)果與之比較,驗(yàn)證2-D DIC的測量精度,針對出現(xiàn)的誤差,尋找誤差源,并進(jìn)行系統(tǒng)修正。
DIC方法是通過處理變形前后被測物體表面的圖像獲得位移和應(yīng)變場信息的測量方法。將變形前后的圖像分別稱為“參考圖像”和“變形后圖像”,利用灰度分布的相關(guān)性求形變量。首先在參考圖像中定義計算區(qū)域(region of interst,ROI),一般為矩形。計算區(qū)域進(jìn)一步被均分為虛擬網(wǎng)格,通過計算每個網(wǎng)格節(jié)點(diǎn)的位移得到全場位移信息。2-D DIC方法的基本原理在于對變形前后兩幅圖像中的相同像素點(diǎn)進(jìn)行追蹤或匹配,如圖1所示,為計算P點(diǎn)的位移,在參考圖像的計算區(qū)域內(nèi)選擇一個以P(x0,y0)為中心的含(2 M+1)×(2 M+1)個像素的正方形參考子區(qū),在變形后圖像中通過一定的搜索方法,按預(yù)先定義的互相關(guān)函數(shù)進(jìn)行相關(guān)計算,尋找與參考圖像子區(qū)的互相關(guān)系數(shù)最大或最?。ㄈQ于所選擇的相關(guān)函數(shù))的以P′(x′0,y′0)為中心的目標(biāo)圖像子區(qū),從而確定P(x0,y0)點(diǎn)在X、Y 方向的位移分量U、V。
為了衡量參考子區(qū)與目標(biāo)子區(qū)的相似程度,必須預(yù)先定義互相關(guān)函數(shù)作為評價標(biāo)準(zhǔn),根據(jù)文獻(xiàn)[4]的研究,選擇歸一化的最小平方距離函數(shù)(zero-normalized sum of squared differences,ZNSSD):
式(1)中,f(xi,yi)為參考圖像中坐標(biāo)(xi,yi)點(diǎn)的灰度,g(x′i,y′i)為變形后圖像中坐標(biāo)(x′i,y′i)點(diǎn)的灰度,fm、gm分別為變形前后圖像子區(qū)的平均灰度。與其它相關(guān)函數(shù)相比,該函數(shù)具有抗干擾性強(qiáng),對目標(biāo)圖像子區(qū)灰度的線性變換不敏感,相關(guān)系數(shù)峰值全場唯一且尖銳等特點(diǎn),因而能更準(zhǔn)確地尋找到整個搜索區(qū)域的相關(guān)系數(shù)極值。
圖1 正方形參考子區(qū)變形前和變形后的示意圖Fig.1 Schematic of a reference square subset before and after deformation
圖2 實(shí)驗(yàn)裝置示意圖Fig.2 Schematic of the experimental set-up
實(shí)驗(yàn)裝置簡圖如圖2所示,采用Instron 5848試驗(yàn)機(jī)進(jìn)行單軸拉伸加載。試驗(yàn)機(jī)載荷傳感器分辨力為0.00001N,最大載荷2kN,位移傳感器分辨力0.00001mm。圖像傳感器為一臺大恒DHHV1303UM CMOS攝像機(jī),分辨力為1280×1024pixel,鏡頭為Computar MLM-3XMP變焦鏡頭。實(shí)驗(yàn)過程中利用磁性底座把攝像機(jī)固定在鋼鐵基座上,以保證攝像機(jī)穩(wěn)定且光軸與試樣表面垂直,拍攝時用冷光源照明試樣。
試樣材料是牌號為T2的紫銅,幾何尺寸如圖3所示。由于相關(guān)運(yùn)算的精度與散斑質(zhì)量關(guān)系密切,因此為了增加散斑圖的平均灰度梯度[5],實(shí)驗(yàn)中的散斑圖是在白漆基底上噴涂直徑約為0.5μm的霧化黑色碳素墨水顆粒得到的,如圖4所示,白框?yàn)樗x計算區(qū)域。
圖3 試樣的幾何尺寸(厚度:0.5mm)Fig.3 Geometrical dimensions of the specimens(thickness:0.5mm)
圖4 散斑圖樣Fig.4 Speckle pattern
現(xiàn)從軟件和硬件兩方面分析了實(shí)驗(yàn)條件及設(shè)備可能引入的誤差,確定最佳拍攝條件,并進(jìn)行了拉伸實(shí)驗(yàn)。
2.2.1 軟件計算誤差
DIC方法是先計算位移場,然后再通過位移場計算應(yīng)變場,先利用雙線性插值法對散斑圖像進(jìn)行灰度的插值,然后利用式(1)計算插值后散斑圖相關(guān)區(qū)域的相關(guān)系數(shù),從而得到亞像素位移,再通過逐點(diǎn)局部最小二乘法[4]來計算位移的導(dǎo)數(shù),即應(yīng)變。
由于位移的誤差會導(dǎo)致應(yīng)變計算不準(zhǔn),因此為確定軟件對實(shí)驗(yàn)圖像的位移計算精度,選取一幅實(shí)驗(yàn)圖像為參考圖像,對其施加0.01~1pixel的模擬位移,比較計算得到位移和虛擬位移之間的差別。
2.2.2 硬件誤差實(shí)驗(yàn)
在保證散斑圖質(zhì)量、光照的均勻、穩(wěn)定及實(shí)驗(yàn)臺隔振的情況下,2-D DIC的硬件誤差主要由以下幾方面引入,因此需逐個分析:
(1)拍攝條件的影響
影響圖像拍攝的主要因素有:光圈、焦距、物距、像距、快門速度(也叫曝光時間)等。而圖像的質(zhì)量直接影響DIC計算的結(jié)果,現(xiàn)通過剛體平移和零位移實(shí)驗(yàn)來檢驗(yàn)拍攝狀況。
剛體平移由于不包含任何變形,所以DIC計算區(qū)域內(nèi)的位移值應(yīng)該是相同的,位移場應(yīng)為一平面;零位移實(shí)驗(yàn)是對靜止的試樣表面連續(xù)拍照,然后對圖像進(jìn)行DIC計算,所得位移場應(yīng)是全為零的平面分布??紤]到軟件存在計算精度,因此若計算得到結(jié)果在軟件計算精度范圍內(nèi)波動,則說明攝像機(jī)的拍攝狀況比較理想。
(a)放大倍數(shù)影響
將鏡頭放大倍率調(diào)至0.3×和1.0×,各做一組零位移和剛體平移實(shí)驗(yàn)。剛體平移是樣品在試驗(yàn)機(jī)上沿豎直方向平移,以0.03mm為步長,平移0.3mm,依次采集10幅散斑圖像。
(b)快門速度的影響
快門速度需配合光源設(shè)置,設(shè)置不當(dāng)也會影響成像,實(shí)驗(yàn)所用的攝像機(jī)快門速度可在1μs~1s范圍內(nèi)調(diào)節(jié),但為了配合白光冷光源,快門速度必須設(shè)為10ms的整數(shù)倍。實(shí)驗(yàn)中,把鏡頭放大倍率設(shè)置在1.0×,在不同快門速度下拍攝零位移圖像。
(2)離面位移實(shí)驗(yàn)
要成功地應(yīng)用2-D DIC實(shí)驗(yàn),要求試樣表面應(yīng)足夠平,且與攝像機(jī)光軸盡可能垂直。然而實(shí)際應(yīng)用中,因?yàn)榧虞d裝置的缺陷,以及材料的泊松效應(yīng)[3,4],故試樣表面會偏離理想平面,離面位移很難避免。為降低離面位移的影響,主要有兩個辦法:一是采用遠(yuǎn)心鏡頭,二是盡可能地將攝像機(jī)放置在遠(yuǎn)離試樣表面的地方,近似形成一個遠(yuǎn)心成像系統(tǒng)[3,4]。由于實(shí)驗(yàn)中使用的鏡頭屬于微距鏡頭,物距較短,不得不考慮離面位移的影響,因此通過數(shù)值計算,討論了實(shí)驗(yàn)中可能出現(xiàn)的離面位移與測量應(yīng)變的關(guān)系。
2.2.3 拉伸實(shí)驗(yàn)
為了驗(yàn)證DIC系統(tǒng)應(yīng)變測量的精度,將多晶銅大試樣單軸拉伸變形的DIC與應(yīng)變片測量的結(jié)果進(jìn)行比較。實(shí)驗(yàn)中,試驗(yàn)機(jī)每拉伸100~200μm記錄一次載荷和應(yīng)變儀讀數(shù),同時采集圖像,直至試樣拉斷實(shí)驗(yàn)停止,每次實(shí)驗(yàn)記錄40~50幅圖像。
總共進(jìn)行了7組實(shí)驗(yàn),編號為A1~A7,為了對比和驗(yàn)證誤差源分析的正確性,各組實(shí)驗(yàn)的拍攝條件不盡相同,其中,A5~A7實(shí)驗(yàn)條件相同,相應(yīng)的攝像機(jī)參數(shù)如表1所示:
表1 A1~A7實(shí)驗(yàn)中攝像機(jī)參數(shù)Tab.1 Camera parameters of A1~A7
施加0.01~1pixel虛擬位移后的散斑圖DIC計算結(jié)果如表2所示。可以看到,劃線處的兩個相對誤差值差別較大,因而斷定實(shí)驗(yàn)所用散斑圖的DIC軟件計算精度大約在0.04~0.05pixel之間,完全可以滿足DIC實(shí)驗(yàn)要求。
表2 不同虛擬位移與DIC計算結(jié)果的比較Tab.2 Comparison of DIC results with virtual displacement
通過軟件計算可知,DIC軟件的位移識別精度滿足實(shí)驗(yàn)要求,那么誤差主要是由外部實(shí)驗(yàn)條件所引入的。
(1)拍攝條件的影響
(a)放大倍數(shù)影響
鏡頭0.3×放大倍率時的零位移實(shí)驗(yàn)X、Y方向的位移場U和V分別如圖5(a)和5(b)所示,可以看到兩個位移場均非理想中零附近的平面,而是呈波浪形分布,在1.0×放大倍率時所得結(jié)果與之類似。
由于剛體平移不包含任何變形,因此計算區(qū)域內(nèi)的位移值應(yīng)該是相同的。然而如圖6所示,在鏡頭0.3×放大倍率時,平移量為0.15mm時的位移場V形狀為一斜面,且其余平移位置的計算結(jié)果與之類似。以Vmax-Vmin來衡量整個計算區(qū)域位移偏差的程度,從圖7中可以看到,在位移量比較大的情況下,偏差有增大的趨勢,而且相比于0.3×放大倍率,1.0×放大倍率時增大得更快,說明鏡頭放大倍率越大,對測量結(jié)果的影響越大。
圖5 0.3×放大倍率,零位移實(shí)驗(yàn)Fig.5 0.3× magnification,null displacement experiment
圖6 0.3×放大率,剛體平移y=0.15mm時的位移場VFig.6 Displacement field along the loading direction(0.3×magnification,y=0.15mm)
圖7 剛體平移實(shí)驗(yàn)位移場Vmax-Vmin比較Fig.7 Comparison of Vmax-Vminin rigid body motion displacement field
零位移和剛體平移計算結(jié)果均說明攝像機(jī)拍攝狀況不理想,導(dǎo)致采集到的圖像不能準(zhǔn)確反映試樣表面變化,從而引入了大的實(shí)驗(yàn)誤差。
(b)快門速度的影響
先前實(shí)驗(yàn)的快門速度均為默認(rèn)的60ms。通過實(shí)驗(yàn)比較,發(fā)現(xiàn)當(dāng)快門設(shè)為400~600ms時,計算出的位移場數(shù)值最為均勻。如圖8所示,快門速度600ms時的零位移實(shí)驗(yàn)的位移場U和V整體呈水平面分布,這一分布特征與理論相符。圖9對比了不同快門速度下零位移實(shí)驗(yàn)的Vmax-Vmin值,可以看到快門速度設(shè)為400~600ms以后,位移場偏差已經(jīng)完全都在0.1pixel以內(nèi),因此可以斷定,在此條件下,實(shí)驗(yàn)誤差最小。
圖8 快門為600ms時的零位移實(shí)驗(yàn)位移場Fig.8 Displacement field of null displacement experiment(shutter speed:600ms)
(2)離面位移
2D-DIC實(shí)驗(yàn)通常是憑經(jīng)驗(yàn)判斷試樣表面與攝像機(jī)光軸是否垂直,因此偏差總是存在的。實(shí)驗(yàn)中使用的微距鏡頭因物距很短,對離面位移非常敏感,故必須具體分析其對測量結(jié)果的影響。
實(shí)驗(yàn)中的離面運(yùn)動通常是離面平移和離面轉(zhuǎn)動的疊加。假定試樣表面與豎直方向夾角為θ,拉伸位移沿著試樣表面向上,且在拉伸過程中角θ不變,仿照文獻(xiàn)[3]建立離面位移模型,如圖10所示。
圖9 不同快門速度條件下,零位移實(shí)驗(yàn)Vmax-VminFig.9 Value of Vmax-Vminin null displacement experiment under different shutter speed
圖10 試樣表面和攝像機(jī)光軸不垂直對位移和應(yīng)變測量的影響Fig.10 Effect of misalignment of the optical camera axis relative to the object plane
若試樣表面和攝像機(jī)光軸之間是垂直的,試樣表面沿著Y方向運(yùn)動,有:
當(dāng)試樣表面向攝像機(jī)方向傾斜θ角時,試樣表面沿著Y′方向運(yùn)動,根據(jù)三角形幾何關(guān)系,Y方向位移為:
取泰勒展式的一階分量,得到:
令Y′方向位移為vY′(Y′)=Y(jié)′2-Y′1,Y′1=Y(jié)′代表試樣上任意的位置,則有:
綜上可得:
此式即在試樣表面任意一個位置Y′,DIC計算得到的應(yīng)變εy與試樣表面的實(shí)際應(yīng)變εY′之間的關(guān)系,可以看到,應(yīng)變與像距L無關(guān)。當(dāng)Y′=0時,可簡化為:
對于實(shí)驗(yàn)所用拍攝系統(tǒng),利用式(8),可算得試樣表面與攝像機(jī)的夾角θ,如圖11所示,它由兩部分組成,首先試驗(yàn)機(jī)上下夾具不對中,存在一定的傾角θ1;其次由于鏡頭自重使攝像機(jī)下傾,及攝像機(jī)基座不水平產(chǎn)生的傾角θ2。若要克服這個夾角θ,可以把攝像機(jī)安裝在一個五維調(diào)整基座上進(jìn)行調(diào)整,直到成像平面與試樣表面平行。
圖11 實(shí)驗(yàn)系統(tǒng)的離面位移Fig.11 Out-of-plane displacement
圖12所示為應(yīng)變片測量和DIC分析所得A2試樣的應(yīng)力應(yīng)變曲線比較。可以看到應(yīng)變片和DIC分析所得結(jié)果之間存在明顯差別,而且其余六組實(shí)驗(yàn)結(jié)果相同。
圖12 A2試樣的應(yīng)力應(yīng)變曲線Fig.12 Stress-strain curve of A2
圖13 A2試樣DIC應(yīng)變的相對誤差隨應(yīng)變的變化趨勢Fig.13 Variation tendencies of strain relative error of A2obtained by DIC method
為了比較DIC與應(yīng)變片測量應(yīng)變的差別,現(xiàn)以相對誤差εi來衡量應(yīng)變誤差
式(9)中,i代表應(yīng)變測量數(shù)據(jù)點(diǎn)序號。以橫坐標(biāo)為應(yīng)變片測量的應(yīng)變,縱坐標(biāo)為DIC應(yīng)變計算的相對誤差,計算發(fā)現(xiàn)五組實(shí)驗(yàn)DIC應(yīng)變測量的相對誤差變化趨勢相同,當(dāng)εy<0.02%時呈散亂分布,而當(dāng)εy>0.2%以后分布逐漸趨于穩(wěn)定,圖13所示為A2的計算結(jié)果。
分析εy<0.2%時相對誤差的散亂分布,主要是因?yàn)閷?shí)驗(yàn)開始階段試樣表面粗糙不平造成的,隨著實(shí)驗(yàn)進(jìn)行,試樣逐漸被拉緊,其表面開始變得平整,這樣的變化在一定程度上改變了試樣與攝像機(jī)之間的距離[10],從而引入離面位移,導(dǎo)致測量的不準(zhǔn)確。當(dāng)εy>0.2%時,試樣僅存在面內(nèi)位移,因而七組實(shí)驗(yàn)的相對誤差都逐漸趨于穩(wěn)定。
為了橫向比較每組實(shí)驗(yàn)間的相對誤差的大小,取平均相對誤差
式(10)中,N為每組實(shí)驗(yàn)測量應(yīng)變數(shù)據(jù)點(diǎn)的數(shù)目,以ε平均為基準(zhǔn)比較實(shí)驗(yàn)結(jié)果,如表3所示:
表3 A1~A7DIC應(yīng)變相對誤差比較Tab.3 Strain average relative error of A1~A7under DIC method
對比發(fā)現(xiàn),相對于其它組實(shí)驗(yàn),當(dāng)快門速度設(shè)為600ms后,A3和A4的誤差大幅減小,這與前面分析吻合,然而和應(yīng)變片測量結(jié)果相比,DIC測量誤差仍然較大,可以斷定,這些誤差主要是由攝像機(jī)的噪聲和實(shí)驗(yàn)中的離面位移引起的。離面位移的影響前文已做理論推導(dǎo),但要具體計算,則需要精確的角度測量才能完成,由于受實(shí)驗(yàn)條件限制,要到以后才完成。
通過前面的論述,可以看到每組實(shí)驗(yàn)中DIC方法和應(yīng)變片測量得到的應(yīng)變的相對誤差大致是穩(wěn)定的,只是不同組實(shí)驗(yàn)的相對誤差不同而已。如果實(shí)驗(yàn)的條件不變的話,其測量的誤差也是恒定的,因此就可以對誤差進(jìn)行修正。
結(jié)合式(6)、式(7),令ε平均=εi,則有
取k=1/(1+ε平均)作為誤差修正系數(shù),其值與拍攝條件有關(guān),因此,對應(yīng)每種拍攝條件,需要事先標(biāo)定一個修正系數(shù)k,才能對此條件下實(shí)驗(yàn)結(jié)果進(jìn)行準(zhǔn)確修正。計算得A1~A7試樣的誤差修正系數(shù)如表4所示。
表4 A1~A7各組實(shí)驗(yàn)的誤差修正系數(shù)Tab.4 Error correction coefficient of A1~A7
可以看到,A5~A7三組實(shí)驗(yàn)因?qū)嶒?yàn)條件相同,系數(shù)幾乎相同,從而可以認(rèn)定在固定了拍攝條件后就有固定修正系數(shù)。其他幾組系數(shù)的稀疏變化則是由于拍攝條件變化引起的。
根據(jù)式(11),現(xiàn)對A5~A7進(jìn)行誤差修正,表4中三個k值均可作為該條件下修正系數(shù),取三者平均值對結(jié)果進(jìn)行修正,如圖14,可以看到經(jīng)過修正的DIC曲線和應(yīng)變片測量結(jié)果非常吻合,可見通過此法可以得到滿意實(shí)驗(yàn)結(jié)果。
圖14 A5~A7試樣經(jīng)過DIC應(yīng)變修正后的應(yīng)力應(yīng)變曲線Fig.14 Modified stress-strain curve of A5~A7
DIC方法是一種很具吸引力的位移和應(yīng)變場測量技術(shù),可應(yīng)用于不同領(lǐng)域中。然而,由于實(shí)驗(yàn)方法的不當(dāng),通常會引入一些誤差。文中對2-D DIC的實(shí)驗(yàn)方法進(jìn)行了研究,分析了多種拍攝條件下產(chǎn)生的實(shí)驗(yàn)現(xiàn)象及誤差,并提出了相應(yīng)的消除和抑制誤差的措施,同時對DIC方法中因離面位移引起的誤差進(jìn)行了理論分析,得出了實(shí)際應(yīng)變計算公式。最后利用2-D DIC方法測量了多晶銅試樣的單軸拉伸應(yīng)變,將結(jié)果與應(yīng)變片測量結(jié)果進(jìn)行比較,檢驗(yàn)文中所用測量系統(tǒng)的精度,針對出現(xiàn)的誤差,認(rèn)為其主要來源于攝像機(jī)噪聲及測試系統(tǒng)的離面位移,鑒于這兩種誤差源比較難以消除,探索了一種利用修正系數(shù)k對誤差修正的方法,通過該方法對系統(tǒng)誤差進(jìn)行修正并得到滿意的測量結(jié)果。
[1]TUNG S,SHIH M,KUO J.Application of digital image correlation for anisotropic plastic deformation during tension testing[J].Optics and Lasers in Engineering,2010,48(5):636-641.
[2]GRYTTEN F,DAIYAN H,POLANCO-LORIA M,et al.Use of digital image correlation to measure large-strain tensile properties of ductile thermoplastics[J].Polymer Testing,2009,28(6):653-660.
[3]SUTTON M A,YAN J H,TIWARI V,et al.The effect of out-of-plane motion on 2Dand 3Ddigital image correlation measurements[J].Optics and Lasers in Engineering,2008,46(10):746-757.
[4]PAN B,QIAN K M,XIE H M,et al.Two-dimensional digital image correlation for in-plane displacement and strain measurement:A review[J].Measurement Science and Technology,2009,20(6):62001-62005.
[5]PAN B,LU Z,XIE H M.Mean intensity gradient:an effective global parameter for quality assessment of the speckle patterns used in digital image correlation[J].Optics and Lasers in Engineering,2010,48(4):469-477.
[6]潘 兵,謝惠民,陳鵬萬,等.數(shù)字圖像相關(guān)測量中鏡頭成像畸變的估計和校正[J].計量學(xué)報,2009,30(1):62-67.
[7]PAN B,XIE H,WANG Z,et al.Study on subset size selection in digital image correlation for speckle patterns[J].Optics Express,2008,16(10):7037-7048.
[8]TONG W.An evaluation of digital image correlation criteria for strain mapping applications[J].Strain,2005,41(4):167-175.
[9]SCHREIER H W,BRAASCH J R,SUTTON M A.Systematic errors in digital image correlation caused by intensity interpolation[J].Optical Engineering,2000,39(11):2915-2921.
[10]CHENG P,SUTTON M A,SCHREIER H W,et al.Full-field speckle pattern image correlation with B-spline deformation function[J].Experimental Mechanics,2002,42(3):344-352.