多介質(zhì)ALE方法流固耦合影響因素及參數(shù)分析

汪 振，吳茂林，孫玉松

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

1 引言

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

隨著當(dāng)今計(jì)算機(jī)技術(shù)的成熟以及數(shù)值分析方法的不斷完善，可以使得流體與固體單獨(dú)建立數(shù)值模型并將流體模型與結(jié)構(gòu)模型進(jìn)行耦合。使得解決入水流固耦合問題成為可能。

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

2 影響因素及參數(shù)分析

2.1 流固耦合方法介紹

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

1)共節(jié)點(diǎn)算法。該方法在網(wǎng)格劃分時使流體與結(jié)構(gòu)在兩者界面處的網(wǎng)格密度相同，節(jié)點(diǎn)數(shù)及位置相同，即使流體和結(jié)構(gòu)共界面共節(jié)點(diǎn)。計(jì)算時流體和結(jié)構(gòu)網(wǎng)格一起變形。

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

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

F=ki·d

(1)

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

2.2 有限元模型的建立

有限元模型的網(wǎng)格采用八節(jié)點(diǎn)六面體單元，算法上計(jì)算域采用LS-DYNA中多介質(zhì)ALE算法，流體材料水和空氣選擇本構(gòu)模型*MAT_NULL和Gruneisen狀態(tài)方程來描述

(γ0+αμ)E

(2)

(3)

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

(4)

計(jì)算過程中所使用的數(shù)據(jù)如表1所示。

表1 水與空氣的狀態(tài)方程參數(shù)

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

表2 柱體模型參數(shù)

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

2.3 關(guān)鍵字說明

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

*CONTROL_ALE是ALE算法控制關(guān)鍵字，參數(shù)設(shè)置上面，選擇介質(zhì)數(shù)值方法為ALE方法，兩次對流間循環(huán)1次，對流方法采用Van Leer+半漂移指數(shù)的對流算法，關(guān)閉光滑選項(xiàng)。

*CONTROL_BULK_VISCOSITY是設(shè)置體積粘度系數(shù)的關(guān)鍵字。總體調(diào)整模型的體積粘度可以減少沙漏變形，人工粘度最初用于處理應(yīng)力波問題，因?yàn)樵诳焖僮冃蔚倪^程中，結(jié)構(gòu)內(nèi)部產(chǎn)生應(yīng)力波，形成壓力、密度、質(zhì)點(diǎn)加速度和能量的跳躍。為了求解穩(wěn)定性，必須引入人工壓力Q，可以將沖擊波間斷面模糊化處理為快速變化但保持連續(xù)傳遞的區(qū)域。如果無耗散，在沖擊波后會出現(xiàn)速度的亂真震蕩。采用該方法可使數(shù)值解不受沖擊波的干擾。LS-DYNA包含有作為壓力附加項(xiàng)的耗散項(xiàng)。體積粘度包括速度散度中的線性項(xiàng)和二次項(xiàng)。

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

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

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

在仿真過程中一些敏感參數(shù)，如網(wǎng)格密度，接觸剛度參數(shù)PFAC(文中簡稱為P值)等對最終計(jì)算的結(jié)果有較大影響，不合理的設(shè)置甚至?xí)霈F(xiàn)穿透，負(fù)體積等問題導(dǎo)致計(jì)算中斷或不合理的結(jié)果出現(xiàn)。為保證數(shù)值計(jì)算的準(zhǔn)確性，得到穩(wěn)定收斂的結(jié)果，本文以大口徑彈體為研究對象，對相關(guān)的參數(shù)設(shè)置進(jìn)行了反復(fù)計(jì)算試驗(yàn)。

本文主要通過對能量檢驗(yàn)的角度來判斷數(shù)值模型是否合理。LS-DYNA理論手冊[9]中說明對內(nèi)能的主要影響因素有以下三點(diǎn)：界面力、摩擦耗散的能量、接觸阻尼耗散的能量。在一個定義良好的模型中，滑動能量應(yīng)該是正的，因?yàn)槟Σ恋绕渌蛩貙?dǎo)致正的接觸能。如果沒有設(shè)置接觸阻尼和接觸摩擦系數(shù)，數(shù)值計(jì)算的結(jié)果可能得到的接觸能為零或一個很小的數(shù)值。這里所說的小根據(jù)以下判斷：在沒有接觸摩擦系數(shù)時，接觸能為模型峰值內(nèi)能的10%內(nèi)可以被認(rèn)為是較為合理的。但是對于入水侵徹模型來說只有阻尼沒有摩擦，能量的消耗發(fā)生在罰函數(shù)耦合時滲透發(fā)生的過程中。對于該類問題，接觸能為模型內(nèi)能的10%甚至更?。换瑒幽芰繎?yīng)該穩(wěn)定零附近或者處于一個較小的正值，可以認(rèn)為模型是較為合理的。然而，數(shù)值計(jì)算結(jié)果出現(xiàn)的一個常見的數(shù)值錯誤是產(chǎn)生負(fù)滑動能量，這對該模型計(jì)算的一個報(bào)錯信號。在數(shù)值上，造成負(fù)滑動能量出現(xiàn)的一般原因如下：

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

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

2.4 網(wǎng)格大小對平頭彈體入水?dāng)?shù)值模擬的影響

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

圖1 彈體網(wǎng)格圖

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

表3 網(wǎng)格類型表

改變計(jì)算域域的ALE網(wǎng)格大小，保持彈體的網(wǎng)格不變，設(shè)置罰函數(shù)耦合剛度系數(shù)PFAC值為0.1，探究計(jì)算域網(wǎng)格大小對數(shù)值模擬結(jié)果的影響。計(jì)算結(jié)果如圖2所示。

圖2 泄漏能曲線圖

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

圖3 滑移能曲線圖

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

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

圖4 加速度曲線圖

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

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

圖5 能量變化圖

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

圖6 加速度曲線圖

2.5 PFAC數(shù)值對圓頭彈體數(shù)值模擬的影響

在確定網(wǎng)格尺寸以及PFAC參數(shù)大小對平頭彈體入水?dāng)?shù)值模型的影響后。為了避免在模擬彈體入水時出現(xiàn)大量節(jié)點(diǎn)的工況一致的情況，建立直徑為30cm，總長度為120cm，頭部直徑30cm的球頭彈體模型，如圖7所示。并對其進(jìn)行合理的網(wǎng)格劃分與參數(shù)設(shè)置。

圖7 球頭彈體模型

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

圖8 能量變化圖

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

圖9 加速度曲線圖

3 結(jié)論

綜合以上數(shù)值模擬測試結(jié)果得出以下結(jié)論：

1)針對網(wǎng)格研究得出：網(wǎng)格達(dá)到一定的密度時，即可以得到較為收斂的結(jié)果，同時節(jié)約計(jì)算機(jī)資源。不合理的網(wǎng)格會產(chǎn)生較大的泄漏能以及負(fù)的滑移能。

2)在合理的網(wǎng)格尺寸及其它參數(shù)下，極大的p值計(jì)算時會出現(xiàn)負(fù)體積并報(bào)錯，滑移能表現(xiàn)為很大的負(fù)值。p值較小時會有滲透出現(xiàn)，滑移能表現(xiàn)為很大的正值。合理的p值能使得滑移能穩(wěn)定在零附近或者取較小的正值。

3)對于平頭彈體，在入水的一瞬間由于多個節(jié)點(diǎn)受到的工況相同，在采用罰函數(shù)耦合的方法時，導(dǎo)致彈體入水瞬間第一次沖擊加速度受P值影響較大。當(dāng)?shù)谝淮螞_擊過后，結(jié)構(gòu)體多節(jié)點(diǎn)受力工況不同，此時P值對沖擊加速度的影響慢慢趨于穩(wěn)定。當(dāng)鈍頭彈體入水，本文研究球頭彈體入水，彈體入水瞬間節(jié)點(diǎn)工況不同，彈體入水時的沖擊加速度受罰函數(shù)P值的影響較小，沖擊加速度在一個穩(wěn)定的范圍內(nèi)。對于平頭物體在入水第一次沖擊，由于多節(jié)點(diǎn)工況相同(類似平板)。當(dāng)這一特殊情況發(fā)生時，需要進(jìn)行能量分析數(shù)值模型是否合理，并結(jié)合相關(guān)理論與實(shí)驗(yàn)進(jìn)行驗(yàn)證。