多介質ALE方法流固耦合影響因素及參數分析

汪 振，吳茂林，孫玉松

(海軍工程大學兵器工程學院，湖北 武漢 430033)

1 引言

對于入水問題的研究[1][2]，自1900年采用閃光攝影技術開始。從早期飛機水上迫降到載人航天工程返回艙的海上降落，魚水雷入水以及當今迅速發展的UUV等入水問題上。入水問題漸漸成為大量學者研究的重要方向。對于入水問題的研究，早期主要采用理論和實驗的方法。理論研究自Von Karman[3]以及Wagner[4]提出并完善楔形體入水的基本理論開始。實驗研究主要通過高速相機觀察結構體入水初期空泡現象以及對入水結構體搭載傳感器進行。

隨著當今計算機技術的成熟以及數值分析方法的不斷完善，可以使得流體與固體單獨建立數值模型并將流體模型與結構模型進行耦合。使得解決入水流固耦合問題成為可能。

大口徑彈體入水問題屬于流固耦合問題，當彈體速度達到200m/s時，在彈體入水的瞬間工況較為復雜，從理論和實驗角度研究存在較大的難度。本文應用LS-DYNA軟件對大口徑彈體高速入水進行數值模擬，主要研究LS-DYNA軟件流固耦合關鍵字參數PFAC以及網格等因素對流固耦合模型的影響，并對數值結果進行分析。

2 影響因素及參數分析

2.1 流固耦合方法介紹

LS-DYNA程序[5][6]是求解流固耦合問題的非線性動力分析有限元軟件，是顯示有限元理論和程序的鼻祖。該程序中的流固耦合算法可以采用以下幾種方法處理：

1)共節點算法。該方法在網格劃分時使流體與結構在兩者界面處的網格密度相同，節點數及位置相同，即使流體和結構共界面共節點。計算時流體和結構網格一起變形。

2)接觸算法。該方法結構和流體在界面處有各自的邊界面，并且可以劃分為不同密度的網格。流體-結構接觸算法在結構與流體間需要定義接觸滑移面，同時定義“光滑約束”條件保證流體網格在變形過程中避免發生畸變。

3)歐拉-拉格朗日耦合算法。該耦合方法與前兩種不同，其主要特點是：分別建立結構與流體模型，同時有限元劃分的網格可以重疊在一起。計算時通過一定的約束方法將結構與流體耦合在一起來實現力學參量的傳遞。LS-DYNA程序中對拉格朗日結構進行約束實現力學參量傳遞的流固耦合算法的約束方法有：加速度約束、加速度和速度約束、罰函數約束等。對于彈體入水問題，需要采用罰函數約束方法。罰函數耦合法是通過追蹤拉格朗日節點(結構，即從物質)和歐拉流體(主物質)物質位置間的相對位移d。檢查每一個節點對主物質表面的貫穿，界面力F就會分布到歐拉流體的節點上。界面力的大小與發生的貫穿數量成正比

F=ki·d

(1)

式中ki表示基于主從節點質量模型特性的剛度系數。由于在每一個時間積分步上都要對等式中的界面力進行求解，所以可以認為F是等式中的一個外力。因此，每一時間積分上都可以對總節點力Fn進行求解，有總節點力才引起了結構加速度、速率和位移的變化。在LS-DYNA程序通過添加關鍵字的方式來進行計算。

2.2 有限元模型的建立

有限元模型的網格采用八節點六面體單元，算法上計算域采用LS-DYNA中多介質ALE算法，流體材料水和空氣選擇本構模型*MAT_NULL和Gruneisen狀態方程來描述

(γ0+αμ)E

(2)

(3)

式(2)、(3)中：P表示壓力；V表示相對體積；E表示單位體積內能；C、S1、S2、S3、γ0為材料常數，其中：C是沖擊波速度us—質點速度up，曲線的截距；ρ為波前介質密度；γ0為Gruneisen常數；S1、S2和S3是us—up曲線斜率的系數。α是對常數γ0的一階體積修正。質點速度up采用式(4)與沖擊波速度us相關聯

(4)

計算過程中所使用的數據如表1所示。

表1 水與空氣的狀態方程參數

柱體的單元類型SOLID164，算法上采用Lagrange算法，選擇模型*MAT_PLASTIC_KINEMATIC來定義材料，計算過程中所使用的數據如表2所示。

表2 柱體模型參數

為了節約計算資源，建立1/2模型，對對稱面進行約束；設置工況為彈體以200m/s的速度垂直入水；空氣域及水域外圍設置無反射邊界條件。

2.3 關鍵字說明

在LS-DYNA前處理模塊中生成關鍵字文件后，對關鍵字文件進行編輯和修改，現將關鍵字文件中個別較為重要的關鍵字[7][8]進行說明，并對其中對模型影響較大的關鍵字參數進行研究。

*CONTROL_ALE是ALE算法控制關鍵字，參數設置上面，選擇介質數值方法為ALE方法，兩次對流間循環1次，對流方法采用Van Leer+半漂移指數的對流算法，關閉光滑選項。

*CONTROL_BULK_VISCOSITY是設置體積粘度系數的關鍵字。總體調整模型的體積粘度可以減少沙漏變形，人工粘度最初用于處理應力波問題，因為在快速變形的過程中，結構內部產生應力波，形成壓力、密度、質點加速度和能量的跳躍。為了求解穩定性，必須引入人工壓力Q，可以將沖擊波間斷面模糊化處理為快速變化但保持連續傳遞的區域。如果無耗散，在沖擊波后會出現速度的亂真震蕩。采用該方法可使數值解不受沖擊波的干擾。LS-DYNA包含有作為壓力附加項的耗散項。體積粘度包括速度散度中的線性項和二次項。

如果div(u)<0，材料處于壓縮狀態，則有Q=C0ρ△x2div(u)2-C1ρcdiv(u)△x，

如果div(u)>0，材料處于拉伸狀態，則有Q=0，

*CONSTRAINED_LARGRANGE_IN_SOLID是實現拉格朗日幾何實體與ALE或者歐拉幾何實體耦合較為重要的關鍵字。設置ALE結構為主實體，在本文中為空氣域和水域；拉格朗日結構為從實體，在本文中為彈體結構；設置NQUAD耦合積分點數為3，即在每個拉格朗日段的表面上分布3*3個耦合點；對于剛性體，設置CTYPE耦合類型為4，即無侵蝕罰函數耦合；設置DIRECT耦合方向為2，僅壓縮方向耦合，該方法具有較好的穩定性及抗干擾性；設置FRCMIN耦合時流體材料的百分比設置為0.3。

在仿真過程中一些敏感參數，如網格密度，接觸剛度參數PFAC(文中簡稱為P值)等對最終計算的結果有較大影響，不合理的設置甚至會出現穿透，負體積等問題導致計算中斷或不合理的結果出現。為保證數值計算的準確性，得到穩定收斂的結果，本文以大口徑彈體為研究對象，對相關的參數設置進行了反復計算試驗。

本文主要通過對能量檢驗的角度來判斷數值模型是否合理。LS-DYNA理論手冊[9]中說明對內能的主要影響因素有以下三點：界面力、摩擦耗散的能量、接觸阻尼耗散的能量。在一個定義良好的模型中，滑動能量應該是正的，因為摩擦等其他因素將導致正的接觸能。如果沒有設置接觸阻尼和接觸摩擦系數，數值計算的結果可能得到的接觸能為零或一個很小的數值。這里所說的小根據以下判斷：在沒有接觸摩擦系數時，接觸能為模型峰值內能的10%內可以被認為是較為合理的。但是對于入水侵徹模型來說只有阻尼沒有摩擦，能量的消耗發生在罰函數耦合時滲透發生的過程中。對于該類問題，接觸能為模型內能的10%甚至更小；滑動能量應該穩定零附近或者處于一個較小的正值，可以認為模型是較為合理的。然而，數值計算結果出現的一個常見的數值錯誤是產生負滑動能量，這對該模型計算的一個報錯信號。在數值上，造成負滑動能量出現的一般原因如下：

1)部件之間互相滑動(沒有摩擦)，由罰函數法進行耦合時固液表面沒有良好分離即發生了滲透現象，這跟摩擦是沒有關系的，這里說的負接觸能是由法向接觸力和法向穿透產生。當一個穿透的節點從它原來的主面滑動到臨近的沒有連接的主面時，如果穿透突然被檢測到，那么就會產生負的接觸能。

2)粗糙的、較大的、不匹配的等不合理的網格，以及長時間計算步長等因素會導致耦合效果不好，界面力無法良好的傳遞，主從界面不能良好的分離，會導致某些方向上數值計算的丟失。

2.4 網格大小對平頭彈體入水數值模擬的影響

首先探究計算域網格大小對平頭彈體入水模型數值模擬的影響。建立直徑為30cm，長度為120cm，頭部帶有5cm的圓倒角的彈體模型，并對其進行網格劃分，如圖1所示。

圖1 彈體網格圖

對水域和空氣域進行網格劃分，設置其為ALE算法，允許網格中存在水和空氣兩種介質，空氣域與水域網格的比例為1：1。將計算域劃分為5種不同的網格類型，不同計算域模型的網格數目如表3所示。

表3 網格類型表

改變計算域域的ALE網格大小，保持彈體的網格不變，設置罰函數耦合剛度系數PFAC值為0.1，探究計算域網格大小對數值模擬結果的影響。計算結果如圖2所示。

圖2 泄漏能曲線圖

圖2為泄漏能變化曲線，從中可以看出：隨著水域及空氣域網格尺寸的增大，泄漏能不斷增大。當網格尺寸較小如類型1、2、3時，泄漏能較小，并且改變網格尺寸對泄漏能的影響不大，模型從泄漏量角度來看較為合理。當網格為類型5的時候，泄漏能明顯增大，說明此時網格模型已經不合理。不同網格模型的滑移能變化曲線如圖3所示。

圖3 滑移能曲線圖

對于流固耦合類問題，滑移能應該穩定在零附近并保持一個較小的正值較為合理。從圖中也可以看出，當網格尺寸較為合理時滑移能變化較為合理。當網格尺寸為類型5時，滑移能較大且開始出現負的滑移能，也說明此時網格尺寸已經不合理。

另一方面，不同網格引起的加速度變化如圖4所示。類型5的網格尺寸數值模型加速度計算結果較為異常。對于類型1、2、3的網格尺寸的數值模型加速度曲線基本重合，但不同的尺寸對入水沖擊加速度的第一次峰值數值影響較大。分析原因為：入水沖擊一瞬間由于類似于平板接觸，導致出現多個節點工況相同的情況發生。在該工況下，采用罰函數耦合方式追蹤拉格朗日節點貫穿主物質的節點的位移受到網格尺寸的影響較大，造成網格尺寸對入水加速度的第一次沖擊峰值影響較大。對于平板類問題入水，應當考慮入水瞬間多節點工況相同的情況發生，需要對數值結果進行一定的理論計算或試驗驗證。

圖4 加速度曲線圖

所以，尺寸較小的網格計算結果相對較為合理。但由于計算采用單點積分，時間步長由最小網格尺寸決定，過小的網格會對計算機的要求較高，并且需要大量的計算時間。另一方面，較為合理的網格尺寸即可以得到收斂性較好的結果。同時數值模擬結果也需要一定的理論計算或試驗驗證進行比較分析以及優化。

針對網格尺寸較為合理的類型2，改變罰函數耦合剛度PFAC數值，來研究PFAC數值對數值模擬結果的影響。分別建立PFAC數值為0.05、0.1、0.2和0.3的模型。能量變化如圖5所示，結果顯示：當網格尺寸較為合理時，不同模型產生的泄漏能均較小。但是從滑移能曲線中可以看出，較大的P值會產生負的滑移能。另一方面當P值過小時會產生較大的正滑移能，說明此時耦合力不足。所以PFAC數值在0.1左右較為合理。

圖5 能量變化圖

同時，由PFAC數值值引起的加速度變化如圖6所示。該數值對入水沖擊加速度第一次峰值數值模擬結果影響較大，P值增加一倍，加速度增量在50%左右。原因仍是由于類似平板入水造成的入水瞬間多節點工況相同的情況發生。導致耦合時，總結節點力Fn的求解受到罰函數耦合剛度PFAC數值的影響較大。

圖6 加速度曲線圖

2.5 PFAC數值對圓頭彈體數值模擬的影響

在確定網格尺寸以及PFAC參數大小對平頭彈體入水數值模型的影響后。為了避免在模擬彈體入水時出現大量節點的工況一致的情況，建立直徑為30cm，總長度為120cm，頭部直徑30cm的球頭彈體模型，如圖7所示。并對其進行合理的網格劃分與參數設置。

圖7 球頭彈體模型

改變PFAC數值的大小來探究該數值對模型的影響，確定較為優化的P值參數。設置PFAC數值分別為0.05、0.1、0.2進行比較計算。數值模擬的結果如圖8所示，從泄漏能曲線來看，與平頭類似，對于合適的網格P值對泄漏能的影響較小。并且泄漏能的變化并不是與P值的增量呈線性關系。另一方面P值對滑移能的影響較大，較高的P值(P=0.2)會產生明顯較大的負滑移能。在P值為0.3時，計算過程報錯并出現負體積，所以在圖8中未列出P=0.3的能量變化曲線。觀察能量信息發現：當P=0.2時滑移能表現為很大的負值。P值為0.05時，出現了一定的滲透，滑移能在上升區間出現較大的正峰值。

圖8 能量變化圖

當彈體為球頭時，由P值引起的加速度變化如圖9所示。由于球頭較平頭入水避免了入水瞬間多節點工況相同的情況發生。該數值對入水沖擊加速度第一次峰值數值模擬結果影響較小，P值增加一倍，加速度增量在20%左右。

圖9 加速度曲線圖

3 結論

綜合以上數值模擬測試結果得出以下結論：

1)針對網格研究得出：網格達到一定的密度時，即可以得到較為收斂的結果，同時節約計算機資源。不合理的網格會產生較大的泄漏能以及負的滑移能。

2)在合理的網格尺寸及其它參數下，極大的p值計算時會出現負體積并報錯，滑移能表現為很大的負值。p值較小時會有滲透出現，滑移能表現為很大的正值。合理的p值能使得滑移能穩定在零附近或者取較小的正值。

3)對于平頭彈體，在入水的一瞬間由于多個節點受到的工況相同，在采用罰函數耦合的方法時，導致彈體入水瞬間第一次沖擊加速度受P值影響較大。當第一次沖擊過后，結構體多節點受力工況不同，此時P值對沖擊加速度的影響慢慢趨于穩定。當鈍頭彈體入水，本文研究球頭彈體入水，彈體入水瞬間節點工況不同，彈體入水時的沖擊加速度受罰函數P值的影響較小，沖擊加速度在一個穩定的范圍內。對于平頭物體在入水第一次沖擊，由于多節點工況相同(類似平板)。當這一特殊情況發生時，需要進行能量分析數值模型是否合理，并結合相關理論與實驗進行驗證。