[發明專利]一種基于距離勢函數可變形三維凸多面體塊體離散單元法

有效

• 申請號: 201811283806.0
• 申請日: 2018-10-30
• 公開（公告）號: CN109408977B
• 公開（公告）日: 2022-07-22
• 發明（設計）人: 趙蘭浩;邵琳玉;劉勛楠;李同春;毛佳;彭大偉
• 申請（專利權）人: 河海大學
• 主分類號: G06F30/23
• 分類號: G06F30/23;G06F119/14;G06F111/10
• 代理公司: 南京縱橫知識產權代理有限公司 32224
• 代理人: 母秋松;董建林

1.一種基于距離勢函數三維可變形凸多面體塊體離散單元法，其特征在于：包括以下步驟：

步驟一，選取研究對象，建立可變形離散單元系統；

步驟二，確定可變形離散單元系統的時間步長Δt；

步驟三，在當前時刻t，采用No Binary Method接觸檢測方法對離散單元外圍一層的網格單元進行接觸檢測，將其中一個離散單元外圍一層產生接觸的網格單元定義為目標單元，與之接觸的離散單元外圍一層的網格單元則定義為接觸單元，并且根據距離勢函數的定義，計算作用于接觸單元與目標單元之間總的接觸力；

步驟四，將步驟三計算得出作用在接觸單元以及目標單元上的總的接觸力用形函數轉化成載荷的等效節點力矢量

步驟五，由步驟四計算得出的載荷的等效節點力矢量計算得出下一時刻t+Δt每個網格單元的位移，速度，以及加速度；

步驟六，根據步驟五中網格單元的位移，更新下一時刻t+Δt每個網格單元節點坐標；

所述采用下式計算：

其中，是網格單元當前時刻t載荷的等效節點力矢量，和分別是當前時刻t的網格單元體力和面力的載荷矢量，N是網格單元節點的形函數，V₀是網格單元的體積，A_0t是當前時刻t網格單元的表面積，A₀是網格單元的表面積。

2.根據權利要求1所述的一種基于距離勢函數三維可變形凸多面體塊體離散單元法，其特征在于：步驟一中所述可變形離散單元系統包括：多個離散單元，以及將離散單元剖分網格后形成的有限單元；離散單元剖分網格后的每個網格定義為網格單元，所述網格單元的節點坐標與有限單元的節點坐標相一致，所述離散單元的參數包括：離散單元的節點坐標、質量、阻尼比、剛度，所述有限單元的參數包括：有限單元的節點坐標、質量矩陣、阻尼矩陣、剛度矩陣。

3.根據權利要求1所述的一種基于距離勢函數三維可變形凸多面體塊體離散單元法，其特征在于：所述步驟二中計算時間步長Δt須滿足：

Δt＝min(Δt_D,Δt_s)

Δt_s≤L/C

其中，Δt_D為離散單元的時間步長；ξ為離散單元的阻尼比，m為離散單元的質量，c為離散單元的阻尼系數，k為離散單元的剛度系數，Δt_s為有限單元的時間步長，L為所有有限單元的最小邊長，C為常數；C的取值范圍為9000-12000。

4.根據權利要求3所述的一種基于距離勢函數三維可變形凸多面體塊體離散單元法，其特征在于：所述C的取值為10000。

5.根據權利要求1所述的一種基于距離勢函數三維可變形凸多面體塊體離散單元法，其特征在于：所述步驟五包括：根據當前時刻t可變形離散單元系統的動力控制方程求解得到當前時刻t到下一時刻t+Δt的加速度增量Δü；其中，M是有限單元的質量矩陣，D是有限單元的阻尼矩陣，K是有限單元的剛度矩陣，Δü是有限單元的加速度增量，是有限單元的速度增量，Δu是有限單元的位移增量，再由廣義Newmark法進行時間域離散，計算出每個有限單元下一時刻t+Δt的位移，速度，以及加速度，由于網格單元與有限單元的節點坐標相一致，因此計算得出下一時刻t+Δt每個網格單元的位移，速度，以及加速度。

6.根據權利要求1所述的一種基于距離勢函數三維可變形凸多面體塊體離散單元法，其特征在于：所述更新下一時刻t+Δt每個有限單元的節點坐標計算公式為：

x(t+Δt)＝x(t)+(r(t+Δt))_x

y(t+Δt)＝y(t)+(r(t+Δt))_y

z(t+Δt)＝z(t)+(r(t+Δt))_z

其中，x(t)、y(t)、z(t)分別為網格單元在當前時刻t，x,y,z方向的節點坐標，(r(t+Δt))_x、(r(t+Δt))_y、(r(t+Δt))_z分別為網格單元的位移在下一時刻t+Δt，x,y,z方向的分量。