基于各向異性分形幾何理論的摩擦非線性數學模型

田紅亮，趙春華，方子帆，朱大林，秦紅玲，劉芙蓉，鐘先友

將式(21)、(42)、(29)代入式(45)得:

摩擦是研究相互作用表面及其有關理論與實踐的科學，涉及物理、化學及材料等學科，摩擦的應用覆蓋機械、材料、化工、石油及冶金等領域。我國的空間固體潤滑材料、微機械制造摩擦學、納米顆粒潤滑材料與技術等成果已應用于我國空間技術領域，解決了航天、航空和船舶等領域多項潤滑技術難題。摩擦是一種復雜的、非線性的、具有不確定性的自然現象，人類目前對于摩擦的物理過程的了解還只停留在定性認識階段，無法通過數學方法對摩擦過程給出精確描述［1－2］，所以準確預測靜摩擦特性會變得困難。傳統的Coulomb摩擦第1定律認為:摩擦力與接觸面所承受的法向載荷成正比例，靜摩擦系數與法向載荷無關，也不決定于接觸面所承受的切向載荷，只隨接觸材料的變化而不同。但隨著近代摩擦學的發展，發現①Coulomb摩擦第1定律中摩擦力與法向載荷不成正比例，而呈現出非線性關系，這是因為法向載荷對摩擦系數的影響與材料之間的真實接觸面積有關，根據粘著理論，只有當材料之間真實接觸面積與法向載荷成正比時，材料摩擦系數才與法向載荷無關，這種情況只發生在微突起主要成塑性變形或微突起成彈性變形，而其高度成指數分布的狀態，在一般狀態下，微突起處于彈塑性變形狀態，故摩擦系數將與法向載荷有關［3］;②靜摩擦系數對于同種材料不變的結論不能很好地反映工程實踐;③在一定條件下，靜摩擦系數隨著法向載荷的變大而變大。

自1508年達芬奇開始把摩擦學引入理論研究的途徑以來，諸多學者一直開展摩擦特性的工作。Chang等［4］依照經典的Greenwood-Williamson綜合彈性接觸理論［5］，涉及粘著力的作用，架構金屬粗糙表面間的靜摩擦系數計算模型。尤晉閩等［6］綜合探究微突起的完全彈性、彈－塑性及完全塑性等不同變形機制，構建干摩擦的靜摩擦系數的統計計算模型，該模型考慮完全彈性區、彈－塑性Ⅰ區的接觸點承擔切向載荷的水平，得到“靜摩擦系數隨著法向載荷的增大而減小”的結論，支持古典摩擦結論［7－9］。表面粗糙度對摩擦系數有較大影響［10］，盛選禹等［11］沿襲 Majumdar-Bhushan 彈塑性接觸分形模型［12］，綜合考慮分形維數D、分形粗糙度G、材料性能參數φ、相關因子K的作用，預測靜摩擦系數的變化規律。但文獻［12］在推導微突起曲率半徑R時存在1個計算原理性瑕疵［13］。另外，文獻［11］在計算方法上具有一些不足之處:①法向總載荷、總摩擦力皆應該是有條件等式［14］，這是因為接觸面在法向載荷作用下，處于塑性變形的微突起由于受到局部接觸載荷作用而發生塑性流動，將不能繼續承受切向載荷，因此在計算靜摩擦力時只包含完全彈性區的微突起;②按文獻［12］的計算方法，引起一些公式錯誤;③給出“f隨φ的增加而降低”的觀點與實際相反;④靜摩擦系數的數字仿真數據在0.001～1 000之間，與工程實際表面靜摩擦系數的分布區間0～1相隔較遠;⑤假定粗糙表面的微觀形貌各向同性。

鑒于以上原因，本文應用各向異性分形幾何理論，探索金屬接觸面最大靜摩擦力的來源機理，構建法向載荷、最大靜摩擦力、靜摩擦系數的非線性數學模型，使用數值模擬，預測摩擦特性，本文分析結果與經典摩擦學論點有偏離，探討導致這些偏離的原因。

1 各向異性分形幾何理論

配對摩擦副實質上是由兩個粗糙表面組成的，在微觀狀態下可以形象化地把兩個粗糙表面的接觸看作是一系列微突起的接觸，當表面受到擠壓，微突起將發生彈性或者塑性變形。

1.1 彈性接觸時一個微突起的實際微接觸面積

由于粗糙表面輪廓對工程研究可認為是連續的，而且隨著放大倍數的增加，在任何一點上都會呈現出更多的具有統計自相似自仿射性的粗糙度細節。Ausloos-Berman函數能模擬一個各向異性分形表面輪廓的幾何形貌，該函數［15］為:

表面輪廓的高度能用式(1)的實部導出:

將式(1)中的一些注釋代入式(2)得:

式中nmax為與測量長度Ls相關的指數，且:

這里fix為向零方向取整。

若 M =10、γ =1.5、L=0.4 μm、G=1.36 ×10－5μm時，式(3)仿真不同分形維數的各向異性表面輪廓高度見圖1。從圖1中可以看出，在x→0，y→0附近，表面輪廓高度急劇降低，故式(3)可仿真在微觀上有缺陷的接觸面，即向材料內凹陷的裂紋、點蝕坑，例如灰鑄鐵的組織結構就如同在鋼的基體上嵌入大量石墨片，石墨的存在相當于在基體上分布大量孔洞，特別是片狀石墨的尖角引起應力集中;當x→0，y→0時，表面輪廓高度很低，幾乎趨于零，事實上，按照式(3)易得二元函數的極限

若M=1，則m=1時，式(3)退變為:

有趣的是，式(5)與y無關。

設橫截微接觸點的半徑為r'，微突起波形式的最長波長為2r'，取Ls=2r'。按式(4)可得基波長微突起對應的頻率指數滿足:

上式在式(5)中對應的一項為:

選擇合適的兩自變量φ1，n0，x使z0(x)最大。如將φ1，n0=0，x=r'代入式(7)，可得微突起變形量為:

當表面受到擠壓時，一個球體狀微突起將發生彈性變形見圖2，其中Pe為一個微突起承擔的法向彈性載荷;R為一個微突起頂端的曲率半徑;r為彈性接觸時一個微突起實際微接觸面積的半徑;E'為兩接觸物體的綜合彈性模量。

圖1 不同分形維數的表面輪廓高度Fig.1 Surface profile height with different fractal dimensions

圖2 球體狀微突起與剛性光滑平面的彈性接觸Fig.2 Elastic contact of spherical asperity and rigid smooth surface

對圖2的直角三角形ocd使用勾股定理得:

化簡上式得:

將式(8)代入式(11)得:

進一步可寫為:

這里a'為一個微突起的橫截微接觸面積，a'=πr'2。

彈性接觸時一個微突起實際微接觸面積的半徑［19］為:

式中:E1、E2分別為兩個接觸物體的彈性模量;μ1、μ2分別為兩個接觸物體的泊松比。

值得指出的是，式(15)與文獻［19］的式(2－13)不同。

一個微突起承擔的法向彈性載荷［4］為:

將式(16)代入式(14)得:

將式(11)代入式(17)得:

將式(18)代入式(14)得:

接著式(19)可變形為:

將式(13)代入式(20)得:

由式(13)的注釋、式(18)可得彈性接觸時一個微突起的實際微接觸面積為:

1.2 完全塑性接觸時一個微突起的實際微接觸面積

一個微突起承擔的法向塑性載荷［12］為:

式中:K為相關因子;σy為較軟材料的屈服強度。

一個塑性變形微突起的實際微接觸面積［4］為:

這里δc為一個微突起的臨界變形量。

當2δ?δc時，完全塑性接觸時一個微突起的實際微接觸面積為:

運用式(13)的注釋、式(18)得:

將式(17)代入式(26)得:

使用式(25)暨(27)，可得完全塑性接觸時一個微突起的實際微接觸面積為:

將式(28)代入式(23)得:

1.3 接觸面承擔的法向總載荷

一個微突起的臨界變形量［4］為:

這里b為待求解的無量綱正常數。

將式(13)代入式(30)得:

式(8)可改寫為:

式(32)除以式(31)得:

當δ=δc時，可得惟一解a'，即為劃分彈塑性區域的臨界微接觸截面積:

將式(34)代入式(33)得:

由于 1≤D≤2，則 a'≥a'c時，可得 δ≤δc，接觸微突起發生彈性變形;而則a'＜a'c時，可得δ＞δc，接觸微突起發生塑性變形。

考慮式(21)與(29)，劃分彈塑性區域的臨界微接觸截面積a'c匹配的一個微突起承擔的法向彈性載荷、法向塑性載荷分別為:

兩接觸粗糙表面在逐漸變大的法向總載荷作用下，并列分布的一系列微突起承擔法向載荷。如果一個微突起承擔的法向載荷函數在點劃分彈塑性區域的臨界微接觸截面積ac'連續，即 Pe(ac')=Pp(ac')，或:

將式(34)代入式(39)得:

需要說明的是式(40)不同于文獻［15］。

將式(40)代入式(34)得:

微突起微接觸截面積的大小分布函數［20］為:n(a')=0.5Dψ1－0.5Da'L0.5Da'－1－0.5D，0 ＜ a'≤ a'L(42)這里ψ為域擴展因子。

整個接觸面的實際接觸面積為:

需要強調式(43)異于文獻［15］的式(8)。

將式(22)、(42)代入式(43)得:

整個接觸面承受的法向總載荷為:

將式(21)、(42)、(29)代入式(45)得:

值得指出式(46)異于文獻［15］的式(13)、(16)。

仿照式(22)得:

將式(47)代入式(44)得:

將式(47)、(48)代入式(46)得:

這里Arc為實際臨界接觸面積，且:

依靠式(49)得:

式中:Aa為整個接觸面的宏觀接觸面積為整個接觸面的接觸率。

將式(52)代入式(50)得:

2 接觸面的最大靜摩擦力與靜摩擦系數

對整個接觸面承受法向載荷和切向力，應力場目前還沒有完整的解析解。現有的有關文獻都根據一定的假設進行了簡化。田紅亮等［21］指出，如果兩個力T(一個微突起承擔的切向力)、Pe(一個微突起承擔的法向彈性載荷)的比例小于0.3，屈服將發生在表層以下，如果比例大于0.3，屈服將發生在接觸點的邊緣。為使應力表達式簡化，假設屈服將發生在接觸點的邊緣。此時接觸點邊緣的接觸面三個主應力分別為:

這里μ為較軟材料的泊松比。

需要說明式(57)不同于文獻［6］的式(2)。

解釋屈服失效的第3強度理論［22］為:

將式(57)、(59)代入式(60)中，有:

將式(22)代入式(61)得:

假設已經發生塑性變形的微突起，由于局部的接觸載荷已經使其發生塑性流動，將不能繼續承受切向載荷，因此在計算最大靜摩擦力時不包括已經發生塑性變形的微突起。沒有達到彈性極限的微突起能夠繼續承受切向載荷，對最大靜摩擦力有貢獻。當切向力逐漸增大時，只有彈性變形的微突起最終將達到屈服，此時的切向力就是最大靜摩擦力。所以最大靜摩擦力為:

將式(42)、(21)代入式(63)中，有:

需要說明式(67)不同于文獻［15］的式(25)。

下面將對本文所提出的理論公式進行數字仿真。在接下來的仿真研究中，先假設某個分形參數確定，然后讓另外一個或二個分形參數變化，研究靜摩擦非線性特征的變化規律，但分形參數他們之間不是相互獨立的，各參數需滿足不等式獻［12］，近似選擇接觸率的上限為0.1，即0.1≥A*r，故適合不等式劃分分形維數D、域擴展因子ψ、相關因子K、材料性能參數φ、表面輪廓頻譜密度相關的參數γ、無量綱分形粗糙度G*有效范圍的依據，且0.5≤K≤3。因此后續的6個仿真參數 D、ψ、K、φ、γ、G*的選擇絕不是簡單意義上的純數學游戲，否則不考慮不等式，就會出現文獻［23］中的3個弊端(圖3、4、6的縱坐標接觸率A*r在0刻度以下都有坐標)。

3 法向總載荷對實際接觸面積的影響

在理論上如何精確計算整個接觸面的接觸率A*

r，是傳統摩擦學的一個懸而未決難點。依據式(53)，整個接觸面承受的無量綱法向總載荷P*是整個接觸面的接觸率A*r的顯函數。在機械工程實際的螺栓結合部中，啞鈴狀測試試件具有三種接觸面:①零件1、零件2間的平面接觸面;②螺栓與零件2的螺紋接觸面;③螺栓頭與零件1階梯光孔的平面接觸面。其中后兩種接觸面與螺栓的共同作用是將數顯扳手控制的擰緊力矩轉化為作用在整個接觸面上的法向總載荷。所以按照式(53)、給定工程預緊力［24］P*，可間接計算接觸率

若 K =0.5、φ = 0.01、γ =1.5 時，選擇不同分形維數D及無量綱分形粗糙度G*，因變量P*對自變量的影響見圖3。

圖3 接觸率隨無量綱法向總載荷的變化規律Fig.3 Varying law of contact ratio with dimensionless normal total load

4 靜摩擦非線性特征的預測

4.1 分形維數D對靜摩擦的影響

若 K=2.8、φ =1、G*=10－10、μ =0.3、γ =1.5 時，分形維數D對靜摩擦特征的影響見圖4。

圖4 分形維數對靜摩擦特征的影響Fig.4 Effect of fractal dimension to static friction features

4.2 無量綱分形粗糙度G*對靜摩擦的影響

若 K=2.8、φ =1、μ =0.3、γ =1.5 時，無量綱分形粗糙度G*對靜摩擦特征的影響見圖5。

4.3 材料性能參數φ對靜摩擦的影響

若 K=2.8、D=1.7、G*=10－10、μ =0.3、γ =1.5時，材料性能參數φ對靜摩擦特征的影響見圖6，其中圖6(a)不同于文獻［15］的圖4(3條曲線沿順時針方向分別對應 φ =0.01、φ =1、φ =0.1，φ =0.01 代表的曲線幾乎與y軸重合)。

圖5 無量綱分形粗糙度對靜摩擦特征的影響Fig.5 Influence of dimensionless fractal roughness to static friction performances

圖6 材料性能參數對靜摩擦特征的影響Fig.6 Impact of material property parameter to static friction traits

通過數值模擬的計算結果可以看出:

(1)依據圖3，接觸率隨著無量綱法向總載荷的增大而增大，且表面越平整(即分形維數越大、無量綱分形粗糙度越小)接觸率增大的相對比例越大。

(2)由圖4(a)～4(b)、圖5(a)～5(d)、圖6(a)可知，靜摩擦系數隨著無量綱法向總載荷的增大而微向上凸弧式增大，并無限接近于1，這是由于隨著法向總載荷增大，潤滑油的擠壓損失隨之加大，摩擦副表面油膜變薄，進入混合潤滑狀態，摩擦表面接觸微突起增多，法向總載荷一部分由油膜承擔，一部分由接觸微突起承擔，其摩擦學特征由潤滑油的流變學和接觸微突起的彈塑性變形共同決定，本文論斷與文獻［4］中結論“靜摩擦系數隨著無量綱法向總載荷的增大而微向上凹弧式減小”、文獻［6］的論斷“靜摩擦系數隨著無量綱法向總載荷的增大而微向上凸弧式減小”恰恰相反，與文獻［25］中有關結果“靜摩擦系數隨著無量綱法向總載荷的增大而微向上凹弧式增大”稍有不同，這是因為文獻［25］在無量綱法向總載荷表示的橫坐標中采用常用對數刻度，而本文在無量綱法向總載荷表示的橫坐標中使用絕對刻度。

(3)當φ=1時，由圖4(a)分析可知，在1.1≤D≤1.4范圍內，靜摩擦系數隨著分形維數的增大而增大;從圖4(a)～4(b)可以看出，在1.4≤D≤1.99范圍內，靜摩擦系數隨著分形維數的增大而下降，故分形維數存在某個最優數，該數讓靜摩擦系數取得極大值。

(4)傳統摩擦學觀點認同，降低表面粗糙度，可以改善摩擦性能——降低靜摩擦系數。但按照圖5(a)～5(d)分析可知，靜摩擦系數隨著無量綱分形粗糙度的減小而微向上凸弧式增加，當表面粗糙度降低到一定程度，達到超光潔狀態時，摩擦副的摩擦性能反而下降。

(5)通過圖6(a)分析可知，靜摩擦系數f隨φ的增加而微向上凸弧式增加，這和文獻［11］的言論“f隨φ的增加而降低”差異很大。

(6)理想的Coulomb干摩擦模型認定，靜摩擦系數與法向總載荷無關，法向總載荷越大，最大靜摩擦力會變大。從圖4(c)～4(d)、圖5(e)～5(h)、圖6(b)可知，法向總載荷和最大靜摩擦力近似呈現出線性正比例的關系，驗證經典摩擦學論點。

5 結論

(1)基于各向異性分形幾何理論，提出一種摩擦非線性數學模型，研究表明該模型可以用于進行摩擦非線性特征的預測。

(2)靜摩擦系數不像理想的Coulomb干摩擦模型描述那樣為恒定值，它由兩接觸材料的分形維數D、分形粗糙度G、域擴展因子ψ、相關因子K、材料性能參數φ、表面輪廓頻譜密度相關的參數γ與法向總載荷P一起決定。

(3)靜摩擦系數隨著法向總載荷、材料性能參數的增大而微向上凸弧式增大，但隨著分形粗糙度的減小而微向上凸弧式增加。

(4)當分形維數較小時，靜摩擦系數隨著分形維數的增大而增大;當分形維數較大時，靜摩擦系數隨著分形維數的增大而變小。

(5)在絕對刻度坐標系下，法向總載荷和最大靜摩擦力近似呈現出線性正比例的關系。

［1］ Peitgen H O.Benoi^t B mandelbrot(1924－2010)［J］.Science，2010，330(6006):926.

［2］ Gomory R.Benoi^t mandelbrot(1924－2010)［J］.Nature，2010，468:378.

［3］謝伯元，張 琦，魯毅飛，等.緊急制動摩擦片的摩擦因數研究［J］.農業機械學報，2006，37(12):33 －35.XIE Bo-yuan，ZHANG Qi，LU Yi-fei，et al.Research on friction coefficient of drum brake’s lining based on urgency brake course［J］.Transactions of the Chinese Society for Agricultural Machinery，2006，37(12):33 －35.

［4］ Chang W R，Etsion I，Bogy D B.Static friction coefficient model for metallic rough surfaces［J］.ASME Journal of Tribology，1988，110(1):57 －63.

［5］Greenwood JA，Williamson J B P.Contact of nominally flat surfaces［J］.Proceedings of the Royal Society of London:Series A Mathematical and Physical Sciences，1966，295(1442):300－319.

［6］尤晉閩，陳天寧.結合面靜摩擦系數的統計模型［J］.振動與沖擊，2010，29(12):26 －29，235.YOU Jin-min，CHEN Tian-ning.Statistical model of static friction coefficient between joint surfaces［J］.Journal of Vibration and Shock，2010，29(12):26 －29，235.

［7 ］ Ibrahim R A.Friction-induced vibration，chatter，squeal，and chaos partⅠ:mechanics of contact and friction［J］.Applied Mechanics Reviews，1994，47(7):209 －226.

［8 ］ Ibrahim R A.Friction-induced vibration，chatter，squeal，and chaos partⅡ:dynamics and modeling［J］.Applied Mechanics Reviews，1994，47(7):227 －253.

［9］王國強，馬若丁，劉巨元，等.金屬摩阻材料間摩擦系數與滑動速度關系的研究［J］.農業工程學報，1997，13(1):35－38.WANG Guo-qiang，MA Ruo-ding，LIU Ju-yuan，et al.Study on the coefficient of dynamic friction between frictional metal material［J］. Transactions of the Chinese Society of Agricultural Engineering，1997，13(1):35 －38.

［10］馬晨波，朱 華，孫見君.考慮粗糙度影響的表面織構最優參數設計模型［J］.華中科技大學學報(自然科學版)，2011，39(8):14 －18，35.MA Chen-bo，ZHU Hua，SUN Jian-jun.Optimal design model of surface texture with surface roughness considered［J］. Journal of Huazhong University of Science and Technology(Natural Science Edition)，2011，39(8):14 －18，35.

［11］盛選禹，雒建斌，溫詩鑄.基于分形接觸的靜摩擦系數預測［J］.中國機械工程，1998，9(7):16 －18.SHENG Xuan-yu， LUO Jian-bin， WEN Shi-zhu. Static friction coefficient model based on fractal contact［J］.China Mechanical Engineering，1998，9(7):16 －18.

［12］ Majumdar A，Bhushan B.Fractal model of elastic－plastic contact between rough surfaces［J］.ASME Journal of Tribology，1991，113(1):1 －11.

［13］田紅亮，朱大林，秦紅玲.固定接觸界面法向靜彈性剛度的改進彈簧分形模型［J］.三峽大學學報(自然科學版)，2012，34(6):83 －88.TIAN Hong-liang，ZHU Da-lin，QIN Hong-ling.Improved spring fractal model for normal static elastic stiffness of stationary contact interface［J］.Journal of China Three Gorges University(Natural Sciences)，2012，34(6):83 －88.

［14］ Tian H L，Li B，Liu H Q，et al.A new method of virtual material hypothesis-based dynamic modeling on fixed joint interface in machine tools［J］.International Journal of Machine Tools ＆ Manufacture，2011，51(3):239 －249.

［15］蘭國生，張學良，丁紅欽，等.基于分形理論的結合面靜摩擦因數改進模型［J］.農業機械學報，2012，43(1):213－218.LAN Guo-sheng，ZHANGXue-liang，DINGHong-qin，et al.Modified model of static friction coefficient of joint interfaces based on fractal theory［J］.Transactions of the Chinese Society for Agricultural Machinery，2012，43(1):213 －218.

［16］ Yan W，Komvopoulos K.Contact analysis of elastic-plastic fractal surfaces［J］.Journal of Applied Physics，1998，84(7):3617－3624.

［17］田紅亮，趙春華，朱大林，等.彈塑性三維各向異性分形表面的接觸分析［J］.三峽大學學報(自然科學版)，2012，34(1):69－73.TIAN Hong-liang，ZHAO Chun-hua，ZHU Da-lin，et al.Contact analysis of elastoplastic three-dimensional anisotropic fractal surfaces ［J］. Journal of China Three Gorges University(Natural Sciences)，2012，34(1):69 －73.

［18］黃曉樂，周正軍，許文年.植被混凝土基材2中草本植物根－土復合體抗剪強度與根系分形特征研究［J］.三峽大學學報(自然科學版)，2012，34(2):59 －62.HUANG Xiao-le， ZHOU Zheng-jun， XU Wen-nian.Research on shear strength and fractal feature of root－soil composite system by two herb plants in vegetation-growing concrete matrix ［J］. Journal of China Three Gorges University(Natural Sciences)，2012，34(2):59 －62.

［19］溫詩鑄，黃 平.摩擦學原理［M］.4版.北京:清華大學出版社，2012:30，219.

［20］Wang S，Komvopoulos K.A fractal theory of the interfacial temperature distribution in the slow sliding regime:partⅠ—elastic contact and heat transfer analysis［J］.ASME Journal of Tribology，1994，116(4):812 －823.

［21］田紅亮，朱大林，秦紅玲.結合面靜摩擦因數分形模型的建立與仿真［J］.應用力學學報，2011，28(2):158 －162.TIAN Hong-liang，ZHU Da-lin，QIN Hong-ling.Fractal model of static friction coefficient of joint interface and its simulation ［J］.Chinese Journal of Applied Mechanics，2011，28(2):158 －162.

［22］ Bhushan B.Analysis of the real area of contact between a polymeric magnetic medium and a rigid surface［J］.ASME Journal of Tribology，1984，106(1):26 －34.

［23］魏 龍，劉其和，張鵬高.基于分形理論的滑動摩擦表面接觸力學模型［J］.機械工程學報，2012，48(17):106 －113.WEI Long，LIU Qi-he，ZHANG Peng-gao.Sliding friction surface contact mechanics model based on fractal theory［J］.Journal of Mechanical Engineering，2012，48(17):106 －113.

［24］田紅亮，劉芙蓉，朱大林.4種NET歐拉－伯努利直梁的動力學特性［J］.三峽大學學報(自然科學版)，2013，35(4):85－93.TIAN Hong-liang， LIU Fu-rong， ZHU Da-lin. Dynamic properties of four NET euler-bernoulli straight beam ［J］.Journal of China Three Gorges University (Natural Sciences)，2013，35(4):85 －93.

［25］田紅亮，趙春華，方子帆，等.金屬材料表面靜摩擦學特性的預測研究——理論模型［J］.振動與沖擊，2013，32(12):40 －44，66.TIAN Hong-liang，ZHAO Chun-hua，FANG Zi-fan，et al.Predication investigation on static tribological performance of metallic material surfaces—theoretical model［J］.Journal of Vibration and Shock，2013，32(12):40 －44，66.