基于瑞利分布的平面磨削溫度場的仿真研究

王　艷　謝建華　劉建國　張　省

上海理工大學，上海，200093

基于瑞利分布的平面磨削溫度場的仿真研究

王艷謝建華劉建國張省

上海理工大學，上海，200093

考慮磨粒在砂輪表面隨機分布對接觸區熱流密度分布的影響，基于未變形磨屑厚度的瑞利分布理論，假定流入工件的熱流密度呈瑞利分布，建立了瑞利分布熱源模型。采用基于有限元法的瑞利分布熱源模型對典型磨削工況進行仿真計算,將計算結果與矩形分布及三角形分布仿真結果進行對比，系統地分析了熱源模型、磨削液對溫度場的影響規律。結合磨削實驗測量值，發現基于瑞利分布的熱源模型仿真結果與實驗測量值吻合較好。

熱源模型；磨削溫度場；有限元仿真；磨削液

0　引言

相對于其他加工方法,磨削去除單位體積材料需要更高的能量輸入，進而導致磨削區溫度急劇升高，高溫會對工件造成不同程度的損傷，如燒傷、相變[1]，因此為提高磨削加工質量，有必要對工件磨削溫度分布進行深入研究。

磨削工件表面溫度分布受熱源模型及佩克萊(Peclet)數的影響非常大[2]。針對磨削區的熱源模型，許多學者做了大量的研究工作。Fang等[3]基于三角形分布熱源模型分析了間歇砂輪平面磨削的溫度分布；毛聰等[4]考慮了未變形磨屑厚度的非均勻變化，研究了基于拋物線分布的熱源模型對平面磨削溫度場的影響；Jaeger[5]則考慮磨削區的熱源分布為矩形分布，計算研究了工件表面溫度的分布情況。未變形磨屑厚度影響著熱流密度在磨削區的分布情況，熱流密度三角形分布考慮了未變形磨屑厚度由零逐漸增至最大值的變化過程，將熱流密度近似為線性變化；熱流密度拋物線分布基于接觸區工件與砂輪的圓弧幾何接觸形態，考慮了未變形磨屑厚度的非線性變化過程,圓弧熱源模型即建立在熱流密度拋物線分布的基礎上。早期的矩形分布忽略了未變形磨屑厚度的非均勻性變化，其計算結果與實際測量值相差較大。

綜上所述，許多學者利用各種熱源模型對磨削工況進行溫度場的數值模擬仿真，并取得了一定的研究成果，但都是對磨削實況經過不同程度簡化后的結果，未考慮磨粒分布的隨機性對未變形磨屑厚度及熱流密度的影響，這與磨削實況相差較大。本文充分考慮磨粒分布的隨機性，基于接觸區未變形磨屑厚度的瑞利分布特性，建立了瑞利分布熱源模型。分析了熱源模型及磨削液對接觸區磨削溫度分布的影響，并結合磨削實驗結果及理論公式分析了基于有限元法建立的瑞利分布熱源模型仿真結果的可靠性。

1　磨削區熱源模型的建立

1.1磨削溫度場數學模型

為了估計磨削區的溫度分布情況及討論有關磨削參數對磨削溫度影響的規律，必須建立一種可以用數學計算而又可以模擬磨削實況的理論模型[6]。磨削時，由于磨削深度較小，接觸弧長相對工件尺寸也很小，因此可以將磨削的熱問題視作帶狀熱源在半無限體表面上移動的情況來考慮。

基于上述假設及傳熱學理論，采用有限元方法建立磨削溫度場的理論模型，且整個磨削溫度場遵守能量守恒定律，得到最終磨削溫度場的熱傳導方程[7](求解域Ω內)為

(1)

式中，Q為熱量；ρ為工件材料密度；c為材料的質量熱容；t為時間；kx、ky、kz為材料沿x、y、z方向的熱導率；T為磨削工件溫度。

此外，在求解域Ω內，應滿足三類邊界條件。

(1)在Γ1邊界上:

(2)

(2)在Γ2邊界上:

(3)

(3)在Γ3邊界上:

(4)

求解域Ω的全部邊界條件滿足：Γ1+Γ2+Γ3=Γ。其中，Γ1為強制邊界條件，為給定的邊界溫度；Γ2、Γ3屬于自然邊界條件，分別表示邊界上的熱流密度條件及對流傳熱條件。

1.2磨削區能量分配比

要得到磨削區溫度場的分布情況，就必須確定流入工件的熱量。磨削熱來源于磨削功率，總的熱流密度可通過下式計算得到：

q=Ftvs/(lgb)

(5)

式中，Ft為切向磨削力；vs為砂輪線速度；lg為磨削弧長；b為磨削寬度。

磨削時磨削液很難進入磨削區，因此忽略磨削液帶走的熱量；而磨屑帶走的熱量占總能量的比例也很小，本文假設磨屑帶走的熱流密度為零。因此有

q=qw+qs

(6)

基于Hahn模型[8]定義的熱量分配比Rws如下：

(7)

式中，βw為工件熱接觸系數；kg為磨粒的熱導率；kw為工件材料熱導率；γ0為磨粒的有效接觸半徑，本文取其典型計算值γ0=15 μm[9]。

1.3熱源模型的建立

磨削區未變形磨屑厚度并不均勻，而是沿著接觸弧長方向由零逐漸增至最大值，如圖1所示。考慮到磨粒在砂輪工作表面的隨機性分布，文獻[10-11]用瑞利概率密度函數描述未變形磨屑厚度的變化情況，如下式：

(8)

式中，β為待定參數，由磨削條件確定，如砂輪組織結構、工件材料性質等；h為未變形磨屑厚度。

圖1　未變形磨屑厚度三維幾何形貌

磨削深度較小時，通常假設流入工件的熱流密度在磨削區呈帶狀均勻分布，但用均布熱源模型計算所得溫度與實際測量結果相差較大，這是因為熱流密度跟未變形磨屑厚度有關。三角形熱源模型將流入工件的熱流密度分布簡化為線性變化；毛聰等[4]考慮了熱流密度在接觸區的非線性分布，建立了圓弧熱源模型；三角形熱源模型與圓弧熱源模型在考慮未變形磨屑厚度對熱流密度的影響時忽略了磨粒的隨機性分布對熱流密度的影響，這與磨削實況相差較大。因此本文充分考慮磨粒的隨機性分布，研究其對熱流密度的影響。由式(8)可知，未變形磨屑厚度在磨削區服從瑞利概率密度函數，因此本文假定磨削區內的熱流密度分布呈現瑞利分布，即

（3）急性壞疽穿孔性闌尾炎：本組7例，表現為：闌尾腫大明顯，性狀不規則、位置較固定、輪廓模糊，邊界不清；闌尾壁可見連續性中斷回聲，壁增厚，各層次不清；闌尾腔擴張無回聲或呈低回聲。超聲可見闌尾區混合性包塊，且存在實性低回聲，強弱不等、回聲雜亂，如可見強回聲光斑于闌尾腔內，有結石，后方有聲影，有明顯腹腔腸間積液，且腸蠕動減弱。

(9)

式中，l為磨削區內距離磨削起點為l的長度值，0≤l≤lg；lg為磨削弧長。

由瑞利概率密度函數可知，當l/β≥4時，其概率F(l/β)≈1；如前所述，磨削熱來源于磨削功率，因此無論采用何種熱源模型，流入工件的熱量應該是相等的；得到瑞利分布函數F(l)應滿足下述條件：

lg/β=4

(10)

(11)

式中，A為修正系數，A值由式(8)～式(11)聯立確定。

由瑞利分布特點可知，其自變量取值區間為0～+∞，而磨削區內弧長僅為lg，不能完整地描述其自變量的變化情況。本文借助待定系數β進行修正，β由式(10)確定，令x=l/β，當x≥4時，f(x)逼近1，如圖2所示，則相應的概率密度函數在0≤x≤4范圍內，其概率也為1，從而保證了自

圖2　熱流密度的瑞利分布函數

變量分布的完整性。圖3顯示了磨削過程中流入工件的熱流密度分布，圖中oxy為工件坐標系，o為磨削起點，o′x′y′為移動熱源動態坐標系，o′為熱源中心，x′為動態坐標系距離熱源中心的坐標值;y′為相應的熱流密度(J/m2·s);x為工件坐標系距離磨削起點的坐標值;y為相應的熱流密度。最終建立磨削區熱流密度模型如下：

(12)

圖3　流入工件的熱流密度

1.4磨削區溫度場的解析解

文獻[12]根據移動熱源理論，得出磨削區工件表面最高溫度如下：

(13)

其中，C為與Pe[2]有關的熱導率,其取值如表1所示；vw為工件進給速度；Pe可由下式得到：

(14)

式中，K為工件材料的熱擴散系數；l′為磨削弧長的一半。

表1　不同Pe情況下的C值

由式(13)可知，在理論溫度值已知的情況下，可求得流入工件的熱流密度：

(15)

2　基于瑞利分布的算例

2.1磨削條件

本文對廣泛應用于模具、量具中的GCr15軸承鋼進行磨削溫度仿真，表2列出了GCr15工件不同溫度下的物理性能，工件尺寸為15mm×12.7mm×5mm。采用Al2O3砂輪，砂輪直徑ds為355.6mm，其熱物理性能參數見表3，具體磨削條件及相關參數見表4。

表2　不同溫度下工件的物理性能[9]

表3　砂輪的物理性能參數[13]

表4　磨削條件及相關參數

2.2基于有限元法的數值模型

2.2.1有限元模型的建立

磨削過程中的熱傳導在空間上為三維傳熱，因此有限元模型采用體單元SOLID90。磨削過程中接觸區內砂輪與工件的接觸時間極短，磨削熱往往來不及傳遞到工件底層而在工件表層聚集，導致工件表面溫度梯度大，故對工件表面網格進行加密處理，從而在保證計算精度的前提下縮短計算時間，幾何實體離散后即得到有限元模型。

2.2.2邊界條件的確定

如前所述，磨削溫度場的全部邊界條件滿足：Γ1+Γ2+Γ3=Γ，其中Γ1為強制邊界條件，用來給定環境溫度邊界條件；Γ2、Γ3分別用來給定熱流密度及對流傳熱邊界條件，即模擬接觸區內的熱傳導過程及工件表面的冷卻液，Γ2、Γ3的加載由移動熱源的位置決定，未定義的表面為絕熱邊界。干磨時，由于空氣的導熱性能差，一般情況下，忽略工件與空氣的對流換熱，即Γ3為零；濕磨時，由于冷卻液很難進入接觸區，因此只在工件的側面及除接觸區之外的工件頂部施加熱導率，其中，水基乳化液的熱導率為82 000 W/(m2·K)，油基冷卻液的熱導率為32 800 W/(m2·K)[14]。

熱流密度加載的位置隨移動熱源的位置變化，本文采用瑞利分布熱源模型確定接觸區的熱流密度的變化規律。磨削過程是個連續不斷的過程，即磨削溫度場的場函數是時間域與空間域的耦合函數，在仿真計算過程中，需要處理兩個問題：一是對移動熱源的連續變化過程進行離散；二是由式(11)可知，描述熱流密度的瑞利分布函數不能進行定積分運算，因此需對瑞利分布函數進行離散處理。熱源的移動過程即是工件材料連續不斷的去除過程，因此將熱載荷的加載進行離散，可以采用“載荷步-ITS”來處理，即每經過一個載荷步，即去除工件上的一層材料，顯然時間取得越短，精度越高，但所需計算時間也越長，經過多次計算，本文載荷步的時間步長定為0.0024 s，每個載荷步分為6個子步；對于隨瑞利分布變化的熱流密度的加載仍可借助積分理論進行離散，離散后的熱流密度分布如圖4所示。

x′為動態坐標系距離熱源中心的坐標值，y′為相應的熱流密度(J/(m2·s))圖4　基于瑞利分布的熱流密度的離散加載

2.2.3仿真結果分析

圖5所示為第10個載荷步接觸區溫度場的分布情況，顯示了此時接觸區內的最高溫度及其所在位置。由圖5可知，該載荷步接觸區內最高溫度為633.85 ℃。圖6顯示的是不同載荷步接觸區的溫度場:在磨削初期，磨削溫度較低，隨著磨削過程的進行，磨削溫度逐漸升高，并在第5個載荷步時最高溫度達到穩定狀態。

圖5　第10個載荷步的接觸區溫度分布圖

圖6　不同載荷步接觸區的溫度場

圖7顯示了第10個載荷步磨削接觸區的溫度分布曲線圖，o′為熱源中心所處位置，由圖7可知，接觸區最高溫度并未處于熱源中心位置。由圖5可讀出沿磨削方向接觸區內任意位置的溫度值，得出最高溫度距離磨削起點5.9 mm，而此時熱源中心距磨削起點為6.5 mm，記移動坐標系內最高溫度所處位置為x′，則磨削區最高溫度位于x′/lg=-0.2處。

圖7　接觸區溫度分布曲線圖

3　熱源模型對溫度場的影響

磨削區熱流密度分布狀況跟未變形磨屑厚度有關，不同的熱源模型，熱流密度分布不同。如圖8所示，采用矩形分布和三角形分布熱源模型計算得到的最高溫度分別為667.5 ℃和648.79 ℃，基于瑞利分布的熱源模型計算得到的最高溫度則為633.85 ℃。本文所用磨削條件及相關參數與文獻[9]中的實驗條件一致，相關參數及磨削條件見表4，實驗通過熱電偶測得工件表面最高溫度為606.15 ℃。各種熱源模型的計算誤差如表5所示，可見不同熱源模型得到的最高溫度相差不大。

圖8　不同熱源模型下接觸區溫度分布曲線圖

熱源模型矩形分布三角形分布瑞利分布計算溫度(℃)667.50648.79633.85實測溫度(℃)606.15606.15606.15誤差(%)10.174.6

如圖8所示，不同熱源模型其最高溫度所處位置不同，矩形熱源模型最高溫度處于A，靠近磨削區的出口，即x′/lg=-0.5；三角形分布最高溫度位置為C，靠近熱源中心；而瑞利分布熱源模型其最高溫度出現在B點，處于x′/lg=-0.2附近。

圖9為磨削區的“磨削溫度-時間”變化曲線。如前所述，熱源模型不同，磨削最高溫度出現的位置有差異，導致三者達到最高溫度的時間不同。如圖9所示，瑞利分布熱源模型溫升至最高溫度時間點介于矩形分布與三角形分布之間；此外，接觸區內矩形分布熱源模型溫升相對比較緩慢，三角形分布溫升速度最快，而基于瑞利分布的熱源模型的溫升速度介于前兩者之間；在溫升至最高溫度之后，三種熱源模型冷卻速度也有差異。

圖9　不同熱源模型磨削溫度-時間變化曲線

可見，熱源模型對磨削區的溫度分布有重大影響。熱源模型影響磨削區溫度最高出現的位置；在磨削區溫升階段，熱源模型影響工件的溫升速度；在磨削溫度冷卻階段，熱源模型不同，冷卻速度也存在差異。此外，雖然不同的熱源模型其最高溫度相差不大，但由于瑞利分布熱源模型考慮了磨粒的隨機性分布，比較接近磨削實況，其計算誤差更小(表5)。

4　磨削液對溫度場的影響分析

本文用對流傳熱邊界條件來模擬磨削液對工件的冷卻作用，并采用瑞利分布熱源模型進行了仿真計算，以此考察磨削液對接觸區溫度場的影響。圖10為濕磨環境下第10個載荷步接觸區內磨削溫度分布，由于工件表面與磨削液的對流傳熱，磨削區沿著寬度方向存在著溫度梯度，工件兩側面的溫度較中心區溫度低150 ℃，而接觸區中心溫度受磨削液的冷卻作用不大，此時接觸區最高溫度為630.03 ℃，與干磨時最高溫度相差0.6%。

圖10　濕磨環境下第10個載荷步接觸區溫度分布圖

圖11顯示了不同熱源模型磨削區溫度沿寬度方向的變化曲線，圖中A、B部分為磨削液對工件兩側面的冷卻作用區域，其影響區域在寬度方向尺寸約為0.7 mm。由圖11可知，磨削液在寬度方向上的影響區域及中心區域磨削最高溫度值受熱源模型的影響不大。

圖11　濕磨下磨削表面寬度方向溫度分布

圖12、圖13是在濕磨及干磨兩種磨削環境下溫度對時間的響應曲線。通過對比發現，不同的熱源模型，在干磨及濕磨狀態下，磨削溫度的變化趨勢基本一致，都呈現出磨削溫度急速上升又迅速下降的變化趨勢。在磨削溫度上升階段，溫升速度基本一致；但在濕磨環境下，由于磨削液的冷卻作用，工件在冷卻階段其速度要大于干磨狀態下的冷卻速度。

圖12　不同磨削環境下基于瑞利分布的溫度-時間曲線圖

圖13　不同磨削環境下基于三角形分布的溫度-時間曲線圖

5　瑞利分布仿真結果的驗證分析

圖14為文獻[9]的實驗裝置圖，磨削條件及相關磨削參數見表4。實驗通過測力儀獲得磨削過程中的切向磨削力Fx，并由式(5)、式(7)、式(13)、式(14)計算出磨削溫度理論值Tmax，Tmax取值見表6。磨削方式為逆磨，磨削行程為152.4 mm，而由仿真過程可知，在第5個載荷步時最高溫度已達到穩定狀態，因此仿真所選用的工件長度足夠；實驗分別采用CBN砂輪及Al2O3砂輪對磨削深度分別為0.0254 mm、0.0381 mm、0.0305 mm、0.0127 mm時濕磨及干磨下的磨削溫度值進行測量，得到共16組數據；其中，對磨削深度為0.0254 mm下的磨削溫度測量值與理論計算值進行標定，得到一個確定的標定因子，其他的溫度測量值乘以這個確定的標定因子，本文對采用Al2O3砂輪的磨削過程進行研究。

圖14　實驗裝置[9]

仿真值(℃)549.73633.85735.63865.77Tmax值(℃)526.64606.15777.71921.40實測值(℃)514.62606.15683.86777.12仿真值與實測值的誤差(%)6.84.67.611.4Tmax與仿真值的誤差(%)4.44.65.46.0

表7為文獻[9]與本文在相同磨削條件及相關參數下所做磨削實驗的實測溫度值。圖15顯示了不同磨削條件下的實測磨削溫度值、瑞利分布仿真值以及Tmax解析解。由圖15可知，隨著磨削深度的增大，上述三者均呈現逐漸上升的趨勢，而且實測磨削溫度值低于瑞利分布仿真值及Tmax解析解，這可能是因為：①實驗過程中的熱電偶的靈敏度不夠理想；②溫度標定存在誤差；③瑞利分布仿真計算過程及Tmax解析解忽略了磨屑帶走的熱量；④測量本身存在誤差。

表7　不同磨削深度下工件表面的磨削溫度

圖15　干磨狀態下的磨削溫度

表6列出了上述三種溫度值之間的誤差分析數據，通過對比發現，基于瑞利分布的磨削溫度仿真值與實測值最大誤差為11.4%，仿真值與Tmax解析解的最大誤差僅為6.0%，表明仿真結果真實可信。

6　結論

(1)熱源模型影響磨削區最高溫度值的位置，矩形分布熱源模型最高溫度位于磨削區出口附近，三角形分布熱源模型最高溫度位于熱源中心附近，而基于瑞利分布的熱源模型最高溫度出現在x′/lg=-0.2附近。

(2)熱源模型影響磨削區工件的溫升及冷卻速度，瑞利分布熱源模型溫升速度介于矩形分布及三角形分布之間，三角形分布溫升速度最快，此外在工件溫度冷卻階段，熱源模型不同，冷卻速度也有差異。

(3)熱源模型不同，磨削區最高溫度值相差不大。

(4)濕磨環境下，沿磨削寬度方向，磨削中心區域受磨削液影響不大，磨削液主要影響工件兩側的溫度。

(5)濕磨環境下，在磨削寬度方向上磨削液對工件兩側面的冷卻作用區域范圍受熱源模型的影響不大。

(6)考慮了磨粒的隨機性對熱流密度分布的影響，建立了瑞利分布熱源模型并進行了仿真計算，仿真結果與實驗測量值及磨削溫度理論值相差不大，因此可為進一步研究磨削溫度場提供理論依據。

[1]Tahvilian A M, Liu Zhaoheng, Champliaud H, et al. Experimental and Finite Element Analysis of Temperature and Energy Partition to the Workpiece While Grinding with a Flexible Robot[J]. Journal of Materials Processing Technology, 2013, 213(12):2292-2303.

[2]Marinescu I D, Rowe W B, Dimitrov B，et al.Tribology of Abrasive Machining Processes[M]. Norwich：William Andrew Publishing, 2004.

[3]Fang Congfu,Xu Xipeng.Analysis of Temperature Distributions in Surface Grinding with Intermittent Wheels[J]. The International Journal of Advanced Manufacturing Technology,2014,71(1/4):23-31.

[4]毛聰,周志雄,周德旺,等. 平面磨削溫度場三維數值仿真的研究[J]. 中國機械工程,2009,20(5):589-595.

Mao Cong, Zhou Zhixiong, Zhou Dewang,et al. Research on Three Dimensional Finite Element Simulation of the Temperature Field in Surface Grinding[J].China Mechanical Engineering, 2009,20(5):589-595.

[5]Jaeger J C.Moving Sources of Heat and the Temperature at Sliding Contacts[J]. Proc of the Royal Society of New South Wales, 1942,76:203-224.

[6]李伯民,趙波. 現代磨削技術[M]. 北京：機械工業出版社,2003.

[7]楊世銘,陶文銓. 傳熱學[M]. 北京：高等教育出版社,2006.

[8]Hahn R S.On the Nature of the Grinding Process[C]//Proceedings of 3rd Internatinal MachineTool Design and Research Conference.Birmingham,England, 1962:129-154.

[9]Pavel R, Srivastava A. An Experimental Investigation of Temperatures during Conventional and CBN Grinding[J]. The International Journal of Advanced Manufacturing Technology,2007,33(3):412-418.

[10]Agarwal S, Rao P V. Predictive Modeling of Force and Power Based on a New Analytical Undeformed Chip Thickness Model in Ceramic Grinding[J]. International Journal of Machine Tools and Manufacture,2013,65:68-78.

[11]Agarwal S, Rao P V. Predictive Modeling of Undeformed Chip Thickness in Ceramic Grinding[J]. International Journal of Machine Tools and Manufacture,2012,56:59-68.

[12]Rowe W B. Thermal Analysis of High Efficiency Deep Grinding[J]. International Journal of Machine Tools and Manufacture,2001,41(1):1-19.

[13]JinTan,Rowe W B,McCormack D. Temperatures in Deep Grinding of Finite Workpieces[J]. International Journal of Machine Tools and Manufacture,2001,42(1):53-59.

[14]Jin Tan, Stephenson D J, Rowe W B. Estimation of the Convection Heat Transfer Coefficient of Coolant Within the Grinding Zone[J]. Journal of Engineering Manufacture, 2003,217:397-407.

(編輯王艷麗)

Research on Finite Element Simulation of Temperature Field Based on Rayleigh Distribution in Surface Grinding

Wang YanXie JianhuaLiu JianguoZhang Sheng

University of Shanghai for Science and Technology,Shanghai，200093

On the basis of Rayleigh’s probability density function of undeformed chip thickness and considering the influences of grain tips’randomness on the heat flux,the heat flux entering the workpiece was assumed to have a Rayleigh distribution,a Rayleigh distribution heat source model was developed to evaluate the temperature field in surface grinding. The FEM simulation for some typical grinding conditions were carried out by using the Rayleigh distribution heat source model,the simulation results were compared with those obtained from triangle heat source model and rectangle heat source model, the influences of heat source model and grinding fluid on the temperature field were investigated. By comparing with experimental results, a good agreement is found between the FEM results obtained from Rayleigh distribution heat source model and experimental results.

heat source model; grinding temperature field; finite element simulation; grinding fluid

2014-06-23

上海市研究生創新基金資助項目(JWCXSL1402)

TG580.14DOI：10.3969/j.issn.1004-132X.2015.04.011

王艷，女，1969年生。上海理工大學機械工程學院教授、博士。研究方向為精密加工與特種加工。發表論文60余篇。謝建華(通信作者)，男，1990年生。上海理工大學機械工程學院碩士研究生。劉建國，男，1962年生。上海理工大學機械工程學院高級工程師。張省，男，1988年生。上海理工大學機械工程學院碩士研究生。