將同幾何分析(IGA)引入到有限元分析(FEA)的框架中，目的是使數值分析模型與計算機輔助設計(CAD)的幾何模型相同。與標準的、低階的有限元單元相比，許多研究論文已經證明了使用更高階和更高連續性的基函數是有益和優越的分析特性。B-樣條曲線(B-splines)和非均勻有理B樣條曲線 (NURBS)是CAD中使用最廣泛的幾何描述，近幾年基于NURBS的有限元技術快速發展并被應用到LS-DYNA中。同時，LS-DYNA也開發了一些適合分析的非結構樣條，以滿足更復雜的工業幾何形狀和精細化分析的需要。本文將主要介紹：

• LS-DYNA中IGA方法概述

• IGA的基本概念

• 如何使用IGA

• IGA在LS-DYNA中的一些功能（下一期分享）

• 演示IGA在碰撞模擬、鈑金成形到頻域分析中的應用等（下一期分享）

基本概念

IGA的提出基于一項現實研究，美國德克薩斯大學計算與應用數學系主任、航空航天工程與工程力學教授 T. J. R. Hughes 在論文中引用了桑迪亞國家實驗室在研究中發現，工程計算從模型建模、劃分網格再到計算結束整個過程中，劃分網格工作花費了大量時間，占據了整個流程近 70-80% 時間，而如果能找到一種解決辦法減少網格劃分的工作量，計算效率將大大提高，這是最初提出IGA想法的原因。

傳統有限元最初由CAD幾何模型進行網格劃分，然后進行計算，這里需要大量的時間在網格劃分上。對于曲面復雜的模型有限元分析雖然細化了網格，但與初始CAD模型只是近似，并不完全一樣，在計算過程中可能出現結果不準確的問題。同幾何分析是利用CAD相同的幾何， CAD幾何模型與用于計算上的幾何模型是同一個幾何模型。

傳統有限元中用的是同參數法（ISOPARAMETRIC），即用來計算的幾何模型和形變幾何模型用相同的參數近似，都是用低階拉格朗日多項式（比如LS-DYNA線性單元一樣）。而同幾何法（ISOGEOMERTRIC）的區別在于，用于計算的幾何模型和CAD模型一致，好處是計算完成之后根據計算結果，如果需要反饋回去修改原來的幾何模型，可以直接在計算模型上修改，而有限元方法可能需要用計算模型近似或mapping回CAD 模型等的操作。總結來說，IGA是用FEA的方式去計算，使用準確的幾何模型，以簡化建模的過程，提高計算效率，同時計算用的幾何與CAD設計用的是同一個幾何，這樣設計和分析、模型修改的數據是互通，使得雙方反饋更流暢，流程更簡易。

IGA中最常用的樣條曲線為NURBS（Non-Uniform Rational B-splines）是一種B-樣條曲線，此外還有T-splines，LR-splines，LRB-splines，HR-splines, THR-splines等曲線，后四種可以做局部的網格細化。本文將主要介紹NURBS，是用得比較多的，也是在LS-DYNA中最早采用的樣條曲線。

基礎知識

隨著越來越多的專家學者開始研究IGA，產生了各種各樣不同的樣條曲線，但并非所有的樣條曲線都適合分析。圖中展示了部分適合用來進行工程分析的樣條曲線，這些曲線未來也將會添加到LS-DYNA中，里面包括：

• Unstructured Splines，可以模擬復雜的形狀；

• Hierarchical splines可進行局部網格細化；

• B-splines是一種比較規則的曲線；

• Bernstein-Bézier forms等

B-樣條曲線（B-splines）和非均勻有理B樣條曲線 (NURBS)是CAD中使用最廣泛的幾何描述，下文將就此進行介紹。

圖中的公式為B-splines的基函數，i代表基函數序號，p為階數，ξ為節點向量（是一系列參數坐標），B-樣條曲線的特征包括：

• 不減的

• 遞歸方式計算

• 數值大于等于零的，沒有負值

• 基于節點向量和階數

• 單位分解

B樣條曲線與普通有限元參數坐標形式對比。

B-樣條曲線的表達類似傳統有限元，選取一些控制點，前面乘上基函數，然后相加。

如何定義B樣條曲線單元？包含兩個概念控制點Control point（紅色圓點）以及節點向量Knot vector（藍色菱形點）。每一個單元是一段非零的節點向量，如上圖，一共有5個單元。

與有限元的網格細分功能一樣，B樣條曲線也可以做網格細分，包含h-refinement以及p-refinement， h-refinement是指element數量會增加，但曲線的連續性不變，通過插入knot value（如0與1之間插入0.5）增加單元數來實現。同時還可以保證curve形狀不變。h-refinement的Control point會增加。

p-refinement代表階數continuity增加，但element數量不變。實現方法：通過加入相同的knot({0,0,0}變成{0,0,0,0}，此時element數量不變但階數提高了1，由于knot vector變化它的Control point數量和原來相比也相應增加（位置可能也不一樣），這是為了保證curve在refine之后形狀保持不變。

與有限元網格細化不同的是，有限元中沒有k-refinement概念（k-refinement為同時進行h-refinement和p-refinement），有限元中沒有對應的細分網格方法。IGA可進行k-refinement，在增加element數量的同時，又能提高階數。

以上圖右上角案例為例，原來的knot vector是{0,0,1,1}，階數p為1。第一種方案先插入knot，使得參數坐標變成{0,0,1/3,2/3,1,1}，3個element，階數p仍然為1。然后提高階數，通過插入重復的knot，參數坐標變成{0,0,0,1/3,1/3,2/3,2/3,1,1,1}，此時階數p為2。

若反過來采用先提高階數的第二種方案，{0,0,1,1}變為{0,0,0,1,1,1}，隨后再插入knot，最后得到結果完全不一樣的參數坐標{0,0,0,1/3,2/3,1,1,1}，最后也是得到3個單元，但他們的基函數是不一樣的。因此從這里可以得出h-refinement和p-refinement順序不能隨意改變，現實中通常采用第二種方案（先提高階數，再增加knot）。k-refinement結合了h-refinement和p-refinement以確保C（p-1）的光滑性以及連續性。第二種方案單元和單元之間可以提高階數（p等于2），但是左邊的第一種方案單元和單元之間的連接處，p仍然是1并沒有變化(光滑性沒有變化)。因此，若需要將整個區域進行全部的Pure k-refinement，除非只有一個單元。

為什么我們需要NURBS呢？如果是一個圓，我們用B-樣條曲線就可以。如果是橢圓形或拋物形，雙曲線形，B-樣條曲線就無法準確的來描述這類圖形，這里就需要用到NURBS。從名字可以看出為非均勻的，每個Control point貢獻的權重可能不一樣。

所以在NURBS里面會引進一個權重參數w，NURBS curve形式與B-樣條曲線類似，只是多考慮了權重參數的線性組合。右邊展示了權重不同導致某個曲線變化的案例，若改變第9個Control point權重，黃色曲線代表 w9 等于0的情況，紫色曲線代表 w9 等于10的情況，也說明Control point的權重不一樣，那么畫出來的曲線是不一樣的。

前面是介紹曲線，這里再介紹曲面和體。NURBS surfaces也是線性組合（linear combination），不同點在于它是二維的，正如上圖所示。公式中二位的基函數是兩個B-Spline曲線考慮權重后基函數的積：

以右圖為例，它展示了R方向與 S方向分別有各自B-Spline的基函數，將他們相乘最后可以得到曲面的基函數。整個曲面呈現向量積的形狀，也比較容易得到所有的Control points數量。假設一個方向Control points數量為n個，另一個方向是m個，總數量就是n*m。總共有兩個互相獨立的Knot vector，以及兩個互相獨立的基函數。

NURBS Solids概念相同，區別在于它是三維（R、S、T三個方向都分別有各自的基函數N、M、L，然后利用向量積得到NURBS Solids的基函數）。NURBS Solids特性與NURBS surfaces基本一致。

上面介紹的B-Spline和NURBS是CAD中最常用的樣條曲線。由于我們在unstructured splines時會通過Bézier extraction轉化成B-樣條曲線，所以這里介紹一下Bernstein-Bézier forms。圖中展示了Bernstein-Bézier forms的基函數表達式，與B-樣條曲線類似，也是通過遞歸方式計算得出基函數。右側圖表分別展示了階數為1、2、3、4時基函數的表現形式。基函數的個數是它的階數加上1。同樣地，Bernstein-Bézier forms基函數與B-樣條曲線的基函數有著類似的特性，均為遞歸方式得到，均為非負值，滿足單位分解（partition of unity和總為1），線性獨立，continuity連續性p的值可以無窮大，p值越大曲線越平滑。

Bézier曲線是P+1個 Bernstein多項式基函數的線性組合，右圖展示的曲線包含4個Control point（p為3），起始點P1和P4在曲線之上，而P2和P3在曲線之外，P1 P2為曲線的切線，P3 P4為曲線的切線（element boundary）。

這里對比B-樣條曲線與Bézier曲線，研究它們是否可以相互轉化？上圖分別展示了B-樣條曲線和Bézier曲線的形狀，同樣為二階函數，knot vector相同均為{0,0,0,1,1,2,2,3,3,3}，可以看到B-樣條曲線一共有5個基函數，5個Control point，而Bézier曲線卻有7個基函數。盡管B-樣條曲線與Bézier曲線的基函數和Control point數量不同，但卻可以畫出形狀相同的曲線（右側紅色曲線），可以說明它們之間存在轉化關系。

B-樣條曲線橫軸第一個element通過P1 P2 P3 三個ctronl point來控制，對應的基函數為N1、N2、N3。Bézier曲線第一個element由B1、B2、B3構成，B-樣條曲線與Bézier曲線之間的轉換通過C^e完成，C^e為Bézier extraction operator。圖中最后一行為C^e的關系式。第二第三個等都有各自不同的C^e關系式，可以將Bézier unstructured splines轉化成B-樣條曲線，方便后續計算。

Unstructured splines。右側案例可以看到，大部分的地方還是tensor product形式，其中紅點Extraordinary point（EP）和T-junctions的地方不是tensor product，導致整個Patch并不是tensor product，但是每個局部的element還是可以寫成tensor product形式，利用C^e將Bézier extraction形式轉換成B-樣條曲線，方便后續計算。對于Unstructure splines，三維的案例也是目前比較活躍的研究領域，因為二維案例中想要改變任何一個形狀，可以對它進行trimming（LS-DYNA中比較成熟的方法之一），但是對于solid，如果用同樣的算法的來trimming會相當復雜，這里并不推薦，但如果用Unstructure splines，在三維中也可以模擬任何形狀，所以在三維中，使用Unstructure splines是比較常見的。

LS-DYNA中的同幾何分析

IGA的概念由Dave Benson于2008年開始引進LS-DYNA，最初稱為GENERALIZED EMEMENTS，直到2010年通過 *ELEMENT_SHELL_NURBS_PATCH 開發NURBS Shells，2015年李利萍博士開始通過*ELEMENT_SOLID_NURBS_PATCH開發NURBS Solids，2019年引入新的關鍵字命名*IGA（從R13版本開始）。這個keyword可以和CAD design的data Structure聯系起來，如幾何結構信息，拓撲結構信息等，而且也比較容易加邊界條件，所以引入這個新的關鍵字*IGA。

LS-DYNA中的IGA經過多年發展，目前已可以支持SMP/MPP，Contact以及各種boundary condition，trimmed shells等。

NURBS Patch，也有element概念（與有限元對應更方便理解）。例如圖中的藍色部分為NURBS element從方便理解上，可以對應有限元中的element。假設有限元中某個殼單元有4個節點來定義，藍色區域代表一個NURBS element，而這里的節點就對應NURBS中的Control point（圖中藍色點），Control point數量多少由階數決定，圖中案例階數為二，說明一個方向有3個Control point，總共就有3*3，9個Control points貢獻到這個element上。右上圖是這個NURBS element的基函數。

其他的element也是由9個control point，單元之間幾個共同的control point，使得element之間有一種overlap重疊(使其具有higher continuity高連續性)，重疊的大小取決于階數。

LS-DYNA *IGA keyword可以與標準CAD data Structure兼容，會包括幾何信息和拓撲信息等。幾何信息包括定義NURBS Patch， Control point，knot vector等；拓撲信息包括Face，edge，point等（更易于施加邊界）。*IGA 關鍵字目前包括殼單元，固體單元，將來也會包括trimmed solid，同時也包括Bézier format 的殼單元和固體單元。

有限元關鍵詞與IGA關鍵詞信息對比，以Shell單元為例。

有限元關鍵詞與IGA關鍵詞信息對比，以Solid單元為例。

（NISR/NISS定義Number of interpolation elements per NURBS-Element(r-/s-direction)，Shell例子（Solid會多一個NIST）在R direction或S direction，藍色區域為NURBS Element，綠色區域為interpolation element。這里定義NISR和NISS均為2，代表R方向與S方向分別使用2個element，來interpolate 這個IGA element。使用interpolation element的原因在于，post processing目前無法直接在IGA上看結果，但是可以將IGA上的結果map 到interpolation element上，且因為interpolation element是Finite element mesh, 也可以用有限元的方式來讀取結果。

此外，interpolation element還可以應用于接觸分析， IGA Surface可以在接觸中用來檢測計算兩個接觸的物體還有多遠的距離開始接觸，interpolation element也同樣可以用來檢索，且方法更簡便。兩個不同的部件之間連接也能應用interpolation element，具體的用法下文有一個例子介紹。）

*IGA_INCLUDE關鍵字可以引入unstructured splines，視頻中案例展示了幾個solid的模型，由于solid模型有一些extraordinary point和孔洞，這幾個模型無法用NURBS或只用一個Patch模擬出完整的形狀。而IGA中的unstructured splines卻可以模擬這些復雜的形狀。右側為左邊模型的模態計算結果，以及另外兩個案例的模態分析結果。

---

更多內容將在下一期繼續分享，歡迎關注我們！

文章來源：2022 LS-DYNA網絡研討會，作者：李利萍博士，Ansys高級研發工程師

視頻鏈接：LS-DYNA IGA等幾何分析介紹

技術校對：王強， Ansys高級應用工程師；整理編輯：俞琴