基于EDEM的月球極區采冰鉆具參數的研究

湯厚鑫，紀玉杰

（沈陽理工大學 機械工程學院，沈陽 110159）

0 引言

水資源是各個國家在月球進行探索的最重要物質及能源原料。1961年，由美國科學家Watson等第一次提出在月球南極和北極撞擊坑的永久陰影區中存在大量水冰資源[1]，他們認為在太陽赤道平面與月球赤道平面夾角小于1.6°的情況下，位于月球南極和北極的某些撞擊坑底部可能永久地處于溫度為40～50 K的陰影之下，在這樣的低溫下，位于月球撞擊坑永久陰影區水資源會以冰的形式存在。自Watson等首次提出在月球兩極永久陰影區存在水冰資源以來，各個國家通過遙感探測、探地雷達等多種手段對月球極區水冰資源開展了探測活動[2-3]。根據探測結果，科學家們推測在月球南極的某些區域會有水冰的存在。

本文以冰為離散元模型，選取鉆頭的回轉轉速、進給速度等重點鉆具參數進行對比分析，確定各鉆具參數對鉆具的受力及對冰的采樣量的影響趨勢，為工程設計提供一定的理論依據。

1 冰離散元模型的建立

本文為了在離散元軟件EDEM中模擬冰的動力作用，考慮將冰離散為具有一定質量和大小的球形顆粒單元，并且單元間具有相應的黏結作用。因此本文考慮建立彈性接觸力模型及平行黏結模型來分別描述冰顆粒單元間的動力作用和黏結作用。

1.1 冰顆粒離散單元間的彈性接觸力模型

冰顆粒的離散元接觸力模型表達了每個顆粒單元間力與位移的主要關系；在顆粒間相互作用的過程中，顆粒單元間會產生彈性作用力和剪切力[4]；其中彈性作用力是由顆粒單元間的相對速度及彈性變形導致的，然而在計算剪切力時需要滿足摩爾-庫倫（Mohr-coulomb）摩擦定律；彈性作用力又可以分為線性和非線性兩種類型；在理論上可以采用非線性力學模型（Hertz-Mindlin）來表示顆粒碰撞的過程[4]；但是在分析顆粒單元的動力學特性時，線性彈性接觸力模型在應用時更加方便和快捷，并且可較為準確地計算出顆粒單元的動力學行為[4]。如圖1所示，建立冰顆粒單元間的線性彈性接觸力模型，其中：顆粒A與顆粒B的質量分別為MA和MB；法向和切向剛度系數分別為Kn和Ks；法向和切向阻尼系數分別為Cn和Cs；μ為摩擦因數。

圖1 冰顆粒單元間的線性彈性接觸力模型

顆粒在法線方向上的相互作用力是由摩爾-庫倫（Mohr-coulomb）摩擦定律來決定，其中法向力Fn是由彈性力Fe和黏性力Fv兩部分組成[4-5]，因此法向力Fn的計算公式為

式中：Fe=KnXn，Fv=-CnX˙n。因此法向力分別為顆粒的法向變形和變形率）。

顆粒在切線方向上，基于Hertz-Mindlin理論和摩爾-庫倫（Mohr-coulomb）的摩擦定律[6]可知，顆粒間的切向作用力計算公式為：

式中，Xs和X˙s分別為切向變形和變形率。

兩個顆粒單元A和B碰撞的法向有效剛度系數計算公式[7]為

而法向阻尼系數按下式計算：

式中:Cn為法向阻尼系數，ζn為無量綱法向阻尼系數，e為回彈系數，M為顆粒單元A或B的有效質量。法向和切向剛度阻尼系數有如下對應關系[7]：Ks=αKn，Cs=βCn。其中：α=0.5，β=0。

在線性接觸模型當中，顆粒單元A和顆粒單元B從碰撞到分離的接觸時間[6]可被定義為

在計算過程中，計算步長為接觸時間的1/50；在線彈性的模型當中，接觸時間是一個與顆粒大小和材料相關的參數。

1.2 冰顆粒離散單元間的平行黏結模型

由于冰顆粒單元間具有凍結作用，所以考慮建立冰顆粒單元間的黏結模型來表達冰顆粒間的凍結作用；黏結模型有2種類型，即接觸黏結模型和平行黏結模型[6]；顆粒單元間的接觸黏結模型只傳遞力的作用，且接觸黏結的黏結作用只發生在接觸點上；然而平行黏結模型不僅可以將2個冰顆粒單元膠結在一起，而且還可以傳遞力和力矩的作用；在本文中主要考慮建立冰顆粒間的平行黏結模型。

平行黏結模型可以被想象為一對法向和切向剛度恒定的彈簧設置在兩顆粒接觸平面的圓盤上[6]，并且兩顆粒間的圓盤在接觸平面的正中央；圖2所示為兩冰顆粒單元間的平行黏結模型；兩顆粒間的力和力矩可轉化為作用在黏結材料邊界上的最大法向和切向應力；其中任意一個最大的應力超過了對應的連接強度，冰顆粒間的平行黏結模型就會發生斷裂[6]；其中：XA和XB分別為冰顆粒單元A和冰顆粒單元B的位置矢量，n為接觸方向的矢量，在兩顆粒接觸點的圓盤上可以傳遞冰顆粒單元間的作用力及力矩，即顆粒單元間的拉力、剪力、彎矩和轉矩；與平行黏結有關的總的力和力矩可以用Fb和Mb表示；每一個向量都可以分解為法向分量和切向分量：

圖2 兩冰顆粒單元間的平行黏結模型

式中：Fb和Mb分別為力和力矩的矢量，分別為力和力矩的法向和切向分量。

在平行黏結模型當中，根據梁的拉伸、扭轉和彎曲理論，黏結圓盤上最大拉應力和最大剪應力的計算公式為：

式中，A、J、I分別為兩顆粒間圓盤的面積、極慣性矩和慣性矩。通常有A=πR2，J=πR4/4，I=πR4/4，其中R為兩顆粒間圓盤的半徑。當兩顆粒圓盤間的最大拉應力及剪應力超過其拉伸強度及剪切強度時，兩顆粒之間的平行黏結模型就會發生斷裂。

2 基于EDEM鉆具鉆取冰塊的仿真試驗設計

為了研究鉆具的進給速度、轉動速度對鉆具的受力和冰采樣量的影響，需要制定出不同的仿真工況進行對比分析。

2.1 鉆具設計

如圖3所示，本文主要設計一種鉆具，該鉆具的內徑為40 mm，外徑為60 mm，鉆具的高度為275 mm，螺旋葉片的螺距為35 mm。

圖3 內徑為40 mm鉆具的結構

2.2 鉆具運動參數的設計

衡量鉆具的回轉鉆進運動的主要運動參數有回轉轉速和進給速度；主要根據所查文獻設計，對10、15、20 rad/s等3種轉動速度進行對比分析[8]。對于進給速度，本文主要設計出一種進給速度為0.15 m/s進行對比分析。因此本文主要針對轉動速度設計3組仿真工況（如表1），對鉆具的受力及冰的采樣量進行對比分析。

表1 仿真工況設計

2.3 冰的離散元參數以及模型的建立

在離散元仿真軟件EDEM中對顆粒模型設置需要進行兩個方面的處理：1）定義冰顆粒單元的幾何模型；2）定義冰顆粒單元的計算參數及在前處理中生成冰顆粒單元。其生成的冰顆粒離散元模型如圖4所示，主要的計算參數如表2所示[7,9]。

表2 冰顆粒離散元模型計算參數

圖4 冰顆粒離散元模型

2.4 鉆具工作過程的仿真

鉆具與冰的相互作用過程中，鉆具以一定的回轉速度和進給速度運動，而冰以剪切破碎為主。本文主要截取T=0.9 s時鉆具的位置，如圖5所示。

圖5 T=0.9 s時鉆具的位置

3 仿真計算結果對比分析

取表1所示的第1、2、3組仿真數據進行對比分析；3組數據分別對應鉆具的3種不同的轉動速度；對于鉆具的受力狀況，分別記錄鉆具所受的力矩。對于冰顆粒的采樣量，在離散元軟件EDEM后處理中測量出鉆具內部冰顆粒的高度。再從離散元軟件的后處理中記錄數據進行對比分析，研究鉆具的轉動速度對鉆具的受力及冰采樣量的影響。不同轉動速度的鉆具對冰的采樣量如圖6所示，不同轉動速度的鉆具力矩如圖7所示。

圖6 不同轉動速度的鉆具對冰的采樣量

圖7 不同轉動速度的鉆具力矩圖

從冰顆粒的采樣量分析，對于不同轉動速度的鉆具，鉆具的轉動速度分別為10、15、20 rad/s，可測量出鉆桿內部冰樣品的高度分別為79.05、86.01、101.11 mm。由此可以看出，隨著轉動速度的增大，鉆桿內部冰樣品的高度逐漸增加。

從鉆具的受力情況分析，對于不同轉動速度的鉆具，鉆具的轉動速度分別為10、15、20 rad/s。如圖7所示，隨著轉動速度提高，鉆具所受的力矩呈現逐漸減小的趨勢。

4 結語

本文首先建立了冰顆粒離散元間彈性接觸力模型及冰顆粒間的平行黏結模型。通過查閱文獻得到冰顆粒離散元模型的計算參數，采用此參數在離散元仿真軟件EDEM中生成了冰顆粒的離散元模型，利用此模型完成了3種工況下月球極區冰鉆取過程的離散元仿真。通過記錄3種工況下鉆具的受力情況及冰的采樣量，分析鉆具的轉動速度對回轉鉆進的影響，并得到以下結論：隨著鉆具轉動速度增大，冰的采樣量逐漸增大，鉆具所受的力矩逐漸減小。