詳解LS-DYNA爆炸仿真計算的模型與算法
炸藥的起爆和爆炸過程是一種快速的化學反應過程，對于該過程的描述，主要存在 CJ 模型和 ZND 模型兩種。LS-DYNA 中包含兩種炸藥爆轟模型：高能燃燒模型和點火生成模型，前者屬于 CJ 模型，后者屬于 ZND 模型。點火生成模型需要輸入炸藥反應率方程參數和未反應炸藥的 JWL 方程參數，對于大多數炸藥來說這些參數缺乏足夠的試驗數據支撐，而工程中使用 CJ 模型就足夠了，因此下面主要介紹高能燃燒模型。
高能燃燒模型根據炸藥上各點距起爆點的距離和炸藥爆速來確定每點的起爆時間，如某個炸藥單元中心離起爆點位置的距離為ri ，炸藥爆速為D ，則該單元起爆時間ti=ri/D ，如果存在多個起爆點，各單元起爆時間按照最近起爆點距離計算。該模型定義爆炸產物壓力 P：
其中s P 是依據產物狀態方程計算得到的壓力， F 為燃燒系數：
其中 Lmin為單元最小特征尺寸。當上式計算的F值大于1時，將F值設置為1，使0≤F≤1 。從式(0.2)可以看出，爆轟波到達前F為零，炸藥單元壓力為零；當爆轟波經過很短時間后，F迅速增大為1，爆炸產物壓力幾乎瞬間從零增大到 Ps 。
流體運動方程的描述，按照所采用的坐標系可以分為拉格朗日方法和歐拉方法兩大類。拉格朗日方法在物質域內求解流體運動方程，坐標系固定在物質上并跟隨物質一起運動和變形，因此也被稱為物質描述；歐拉方法在空間域內求解流體運動方程，坐標系固定在空間不動，物質在計算網格之間流動，因此也稱為空間描述。
LS-DYNA 中采用任意拉格朗日歐拉方法（ALE）來描述流體運動。該方法在拉格朗日坐標和歐拉坐標之外引入一個可以任意運動的參考坐標系，計算域基于參考域，可以在空間中以任意形式運動。采用ALE 算法的網格同時具有歐拉網格和拉格朗日網格的優點，網格可以隨物質一起運動，也可以固定在空間不動，甚至可以在一個方向上隨物質運動，而另一個方向上固定不動。
任意拉格朗日歐拉方法描述下的物質導數為：
其中f為物理量，為質點X 的速度，
為參考點ξ 的速度，也即計算網格運動速度。當
時，計算網格在空間固定不動，退化為歐拉描述；當
時，計算網格隨同物體一起運動，退化為拉格朗日描述；當
時，計算網格在空間中獨立運動，對應于一般的ALE 描述。
由于爆炸產物和空氣的粘性很小，而且爆炸流場運動被視為絕熱等熵運動，一般都采用無粘性可壓縮流體運動方程來描述爆炸流場。通過式(0.3)將歐拉方法描述的無粘性可壓縮流體力學方程變換得到ALE 方法描述的控制方程：
式(0.4)－(0.6)結合空氣和爆炸產物的狀態方程可以構成封閉的控制方程組。
ALE 方法引入了運動網格，通過在移動邊界法向上采用拉格朗日方法，可以很簡單地描述邊界運動，解決了歐拉方法中移動邊界描述困難的問題，給計算帶來了很大的方便，但計算過程中需要確定網格的位置。
LS-DYNA 程序中提供了簡單平均算法、體積加權算法、等參算法、等勢算法以及混合算法等用于ALE 運動網格位置的確定。但由于爆炸流場計算過程中，爆炸產物和空氣界面存在很大的壓力和密度梯度，采用以上任何一種算法都會產生異常小的界面網格，從而導致計算無法正常進行。因此爆炸流場計算中一般僅在邊界上采用物質描述，使邊界節點速度與界面法向運動速度相等，對于除邊界節點外的網格要關閉程序中的網格運動算法，使內部網格退化為空間描述。
當需要考慮殼體影響時，殼體和流場邊界可通過共用節點聯結，殼體為爆炸流場提供運動邊界條件，爆炸流場為殼體施加壓力載荷條件，在每一個時間步分步求解即可實現爆炸流場和殼體結構的流固耦合；而當采用剛性壁面假設之后，ALE 方法進一步退化為純粹的歐拉方法。
炸藥爆炸后，爆炸容器內存在爆炸氣體產物和空氣兩種物質，這兩種氣體的流場都可以通過上述無粘性可壓縮方程描述，但計算過程中需要區分兩種不同介質，并捕捉兩種物質的界面。LS-DYNA 中通過定義物質組號來區分不同介質，采用楊氏流體體積法（Yong’s VOF）來捕捉兩種物質的運動界面，具體過程是：采用多物質單元劃分爆炸流場計算域，一個單元中允許同時存在多種物質；先假設界面沿單元邊界，根據與節點相鄰的所有單元中存在的物質計算各個節點上的物質體積分數；由同一單元網格各節點的物質體積分數梯度確定界面的法向，并構造該單元內界面；計算每個時間步內通過四周流到相鄰單元的流體體積，修改網格單元和相鄰單元中的流體體積分數；由各個邊界單元內界面組成整個物質界面。
歐拉和ALE 描述的運動方程求解一般有兩種方式：全耦合求解和算子分裂法，前者是在整個計算域同時求解運動方程，后者將每個時間步分為拉格朗日階段和對流階段依次計算。全耦合求解一個計算單元只能存在一種物質，不適合求解多物質問題。
LS-DYNA 程序中采用算子分裂法求解。拉格朗日階段采用有限元方法計算由于外力和內部應力產生的速度、壓力和內能變化以及現時密度，單元采用單點積分并通過沙漏粘性控制零能模式，引入人工粘性以捕捉沖擊波，時間推進采用二階精度的中心差分法；求得這些參數后再進行界面捕捉，構造多物質內界面。對流階段采用有限體積法計算通過單元邊界的質量、動量和能量通量，通量的計算可以采用一階精度的迎風算法或者二階精度的Van Leer 對流算法；該階段時間步不發生變化，保持與拉格朗日階段一致。爆炸流場計算一般采用Van Leer對流算法，因為這種算法不僅具有二階精度，而且具備總變差遞減（TVD）性質。
其中P 為空氣壓力；γ 為多方指數；ρ 為空氣現時密度，ρ0為初始密度；E為空氣能量密度。理想氣體狀態方程可以通過設置LS-DYNA 程序中預定義的線性多項式狀態方程的相關常數得到。
目前存在的爆炸產物狀態方程，按照是否顯含產物組分可以分為兩類，一類是顯含產物組分的狀態方程，如BKW 方程、KHT 方程等，另一類是不顯含產物組分的狀態方程，如多方方程、JWL 方程等，其中JWL 方程是數值計算中應用最為廣泛的狀態方程。
JWL 方程是1968 年由美國LLNL 的E. L. Lee 在H. Jones 和M. L. Wilkins 的工作基礎上提出的半經驗狀態方程。JWL 等熵方程形式如下：
其中為相對比容，即現時比容
與初始比容
之比，A 、B 、C 、1 R 、2 R 、ω 為六個參數。
根據熱力學第一定律和等熵條件，比內能增量，將(0.8)式代入并積分可得比內能
。
根據式(0.8)和(0.9)，可消去參數C ，得到
其中為以初始體積表示的能量密度，初始能量密度
。式 (0.10)為 LS-DYNA、AUTODYN、CTH、MESA 等爆炸動力學計算軟件中所普遍采用 JWL 方程形式。該形式只要求輸入
五個參數值和初始能量密度
，結合炸藥高能燃燒模型要求輸入的炸藥密度
、爆速 D 和 CJ 爆壓
，可以得到爆炸時刻的
等，然后通過增加相對體積微量
，根據式 (0.10) 和原先的
值，計算得到新的
值，新的
值調整為
，如此反復一直計算。
炸藥 JWL 方程參數的確定需要通過圓筒試驗和二維流體彈塑性數值計算相結合的方法確定：先進行圓筒試驗，將待測炸藥裝入紫銅管，一端起爆，用高速攝像儀記錄下銅管外徑的運動軌跡；設定一組值，根據式(2.11)、(2.12)和(2.13)求得A 、B 、C ，利用假設的JWL 方程通過二維流體彈塑性程序數值模擬炸藥驅動圓管的外徑膨脹軌跡，如果數值計算結果與和試驗結果相對誤差小于1%，則假設參數即為真實JWL 方程參數，如果不滿足相對誤差要求，則繼續調整系數，直到和試驗結果相對誤差小于1%為止。
式中為炸藥等容爆熱，
為炸藥CJ 狀態時的比容，
文章來源：CAE仿真之家
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
