三角形單元的案例
巖土問題二維有限元實體單元類型(element types)的選擇
1 引言
有限元分析的單元類型以及單元的質量影響著數值模擬的結果，在論壇中有人問到了這個問題，因此本文僅對二維有限元實體單元類型作簡要的討論，討論的條件包括: (1) 巖土工程問題使用的主要單元類型；(2) 不包括結構元、界面元【板樁墻(Sheet Pile Wall)模擬---FEM中的界面元】以及一些特殊類型的單元，例如不連續元【使用Plaxis的不連續單元(Discontinuity Elements)】。
2 單元類型
(1) 3節點的三角形單元T3
三節點的三角形單元(3 Noded Trangles)是有限元分析中最簡單的實體單元，在實踐中已經很少使用，典型的應用包括RS2, ADONIS【二維FEM分析軟件ADONIS(V3.50.1)命令流】，EnFEM【EnFEM---一個小型的有限元分析教學軟件】等。
(2) 6節點的三角形單元T6
六節點的三角形單元(6 Noded Trangles)是有限元分析中經常使用的實體單元，6節點的三角形對位移進行了二階插值，數值積分涉及三個高斯點。如果劃分足夠多的單元，那么T6在標準的平面應變分析中可以得出非常精確的結果，不過在軸對稱應用或在使用強度折減法(phi-c )時應小心使用，它很可能高估其破壞載荷和安全系數。因此最好使用下面的T15單元。T6典型的應用包括Plaxis, RS2, ADONIS, EnFEM等。
(3) 15節點的三角形單元T15
15節點的三角形單元(15 Noded Trangles)能夠更精確地表示應力應變關系，特別是對于軸對稱問題，目前只有Plaxis提供這種單元。T5對位移進行了四階插值，對12個高斯點(應力點)進行數值積分。T15單元能夠精確計算一些巖土工程特殊問題的應力，例如不可壓縮土層中的塌陷計算。
展開 三角形常應變單元解平面問題實施步驟與注意事項
利用上面討論的三角形常應變單元解平面問題，其具體步驟可歸納如下:
1)將要計算的彈性體劃分成三角形單元。對結點進行編號，列出結點坐標作為輸入信息。
(2)對單元進行編號，列出單元三個結點的號碼作為輸入信息。
(3)計算載荷的等效結點力，把等效結點力作為輸入信息。
(4)按照(6)式計算各單元的常數bi、ci、bj、cj、bm、cm，再按照(4)計算2A。
(5)按照(35)式計算各單元的剛度矩陣。
(6)形成整體剛度矩陣。
(7)處理約束及消除剛體位移。
(8)解線性方程組(32)式，求結點位移。
(9)按照(20)式計算應力矩陣，再按(18)式計算單元應力。根據需要計算主應力和主方向。
通常步驟(4)至(9)均由計算機來完成，而步驟(1)至(3)可以用手工完成，也可由計算機來完成。在實現以上各步驟時，為了達到一定的計算精度，節約計算機存儲量，縮短計算機運行時間等目的，還需要注意下列事項。
1、利用對稱性
在劃分單元前要研究一下，計算對象是否有對稱變形或反對稱變形存在，從而確定是否需要取整個物體，還是取部分物體作為計算模型。例如圖8a所示受純彎曲的梁，它對于x,y軸都對稱，而載荷對于y軸對稱，對于x軸反對稱。可見，應力和應變亦將具有同樣的對稱和反對稱特性，所以我們只需計算1/4梁就行了。分離體如圖8b所示。對于刪去部分結構的影響可以這樣考慮：對于處于y軸對稱面內各結點的x方向位移和y方向分布力都應等于零，而對于處在x軸反對稱面上的各結點的x方向位移和y方向分布力亦都應等于零。這些條件相當于安置如圖8b中的約束。圖中o點上安置y方向的約束是為了消除剛體位移而設置的。又例如在分析圖9中所承受均勻壓力的厚壁圓筒時，根據結構和載荷軸對稱的性質，我們可以取出一個小扇形(圖中陰影部分)進行計算。
展開 有限元法講解及運用常應變三角形單元解彈性力學平面問題（FORTRAN語言編寫有限元法程序算例）
元素（單元）的形狀原則上是任意的。二維問題一般采用三角形單元或矩形單元，三維空間可采用四面體或多面體等，每個單元的頂點稱為節點（或結點）。
步驟2：單元分析
進行分片插值，即將分割單元中任意點的未知函數用該分割單元中形狀函數及離散網格點上的函數值展開，即建立一個線性插值函數。
步驟3：求解近似變分方程
用有限個單元將連續體離散化，通過對有限個單元作分片插值求解各種力學、物理問題的一種數值方法。有限元法把連續體離散成有限個單元：桿系結構的單元是每一個桿件；連續體的單元是各種形狀（如三角形、四邊形、六面體等）的單元體。每個單元的場函數是只包含有限個待定節點參量的簡單場函數，這些單元場函數的集合就能近似代表整個連續體的場函數。根據能量方程或加權殘量方程可建立有限個待定參量的代數方程組，求解此離散方程組就得到有限元法的數值解。
有限元法已被用于求解線性和非線性問題，并建立了各種有限元模型，如協調、不協調、混合、雜交、擬協調元等。有限元法十分有效、通用性強、應用廣泛，已有許多大型或專用程序系統供工程設計使用。結合計算機輔助設計技術，有限元法也被用于計算機輔助制造中。
展開 【JY】 ABAQUS子程序UEL的有限元原理與應用
算例中采用的三角形單元為ABAQUS單元庫中的二維三角形單元CPS3。
算例一：平面三角形單元
一、三角形單元有限元計算參數
圖3 平面三角形單元示意圖
三角形單元的相關計算參數如表2所示：
表2 三角形單元計算參數
該三角形單元共有三節點（如圖1所示），每個節點有2個自由度
二、邊界條件及載荷設置
將節點1與節點3的位移自由度施加約束，在節點2的水平x方向施加100N的節點力。
基于Matlab的有限元網格自動生成算法 | CST、LST單元網格
單元網格的形成實際上屬于有限元計算中的前處理部分，即確定單元節點信息，當模型較為復雜時，用戶可在Abaqus、Ansys等大型商業有限元軟件中進行建模，導出網格信息。
當模型較為簡單時，如二維平面板模型，用戶可基于一些較為基礎的網格生成算法，在自己的程序中通過控制模型長、寬等信息，即可生成有限元網格。
看似應用有限，但是在一些比較復雜的領域內，往往需要先在簡單的模型中得到理論驗證，如此以來，有利于自編程代碼的完整性，即前處理、內核計算、后處理于一體。
本篇推文，木木就帶著大家，學習一下CST、LST單元網格的自動生成。
CST單元網格
單元自動網格劃分
如下圖所示，為3節點三角形單元網格生成示意圖，圖中NXE和NYE分別是模型橫向和縱向單元個數，dhx和dhy分別是單元的橫向、縱向長度。
展開 有限元分析的底層邏輯是什么？
剛度矩陣：將上述應變能公式轉化為剛度矩陣形式，得到單元剛度矩陣：
??=??????(1?1?11)
其中，??
是楊氏模量, ??
是截面面積，??
是單元的長度。
2.2 2D 單元（如三角形單元）
三角形單元（線性單元）
假設一個簡單的三角形單元有三個節點，節點1、節點2和節點3。其剛度矩陣的推導過程也采用能量原理。 步驟： - 位移場：假設每個節點的位移是線性插值的。位移場可以寫作：??(??,??)=??1(??,??)??1+??2(??,??)??2+??3(??,??)??3
其中??1(??,??)
，??2(??,??)
，??3(??,??)
是形函數。 - 應變能：應變能是通過單元內的應變能密度（與應力和應變的關系）計算的。由于三角形單元涉及到二維應變，通常通過應變-位移矩陣來計算。 - 剛度矩陣：經過推導后，三角形單元的剛度矩陣為：??=∫??????????,????
其中，??
是應變-位移矩陣，??
是材料的剛度矩陣，??
是單元的面積。
三角形單元剛度矩陣（簡化版）
對于簡單的線性三角形單元，剛度矩陣一般可以通過以下積分公式推導：
??=??4(120?12012?12?12?1212)
其中??
是三角形單元的面積。
2.3 3D 單元（如四面體單元和六面體單元）
在三維問題中，常用的單元有四面體單元和六面體單元。
四面體單元 四面體單元的推導較為復雜，通常需要通過形函數和積分來獲得剛度矩陣。這里給出大致的推導思路。 步驟：
- 位移場：四面體單元的位移場是通過節點的位移進行插值，通常是多項式形式的。
- 應變能：通過應變-位移關系，計算應變能。應變-位移矩陣 ??
是由形函數的導數組成的。
展開 基于python分析中心孔的均勻薄板受到單軸壓力將有限元的近似解與基于彈性力學理論的精確解進行對比 ￥59.9
首先對此結構進行單元剖分，確定單元按照有限元的分析流程，要先對此結構進行單元剖分，確定單元與結點編號、以及單元的自由度編號。因為這里是平面應力問題，所以可以采用常應變三角形單元進行網格劃分，并且采用的是非結構化的網格。
在該結構中采用的是常應變三角形單元，在整體坐標系中單元剛度矩陣均用下式進行計算：
對于三個節點，有
還可以得到
代入可以得到：
其中A為三角形單元面積，給定三角形節點坐標，可以通過下式得到：
三角形單元應變能為
在三角形單元中，外力做功為
單元總勢能為
單元處于平衡狀態，總勢能最小，有
代入上式得到
其中k即為單元剛度矩陣。
展開 ABAQUS 單元-節點排布順序解析（重點講解分析方法）
進行ABAQUS UEL二次開發、或者研究界面問題的時候，比如cohsive單元界面問題，會涉及到單元-節點的排布順序。ABAQUS inp文件中的單元-節點順序排布很有講究，不能搞錯，這是因為節點的排布順序與內部程序有關聯，內部的程序我們不得而知，但是節點順序的規律可以從inp文件中看到，再對比cae節點中的節點編號可以總結規律。以下內容介紹一些經典的實體單元-節點排布順序。
這部分內容其實挺好玩兒的！！！
（一）三角形單元
以一個10x10（無單位）的挖孔板為例，劃分的網格為：（注意到被單元為30、138和23包圍的三角形沒有編號，這是因為那個三角形不是單元，我在那里挖了個孔，因為背景顯示的問題，網格顏色和背景一樣，故沒有區分。作此說明）
這里我們首先查看節點的排布方式。分別以單元1、12、123為例。
編號為“1”的單元，在inp文件中，單元-節點的編號為：
1, 61, 1, 2
在cae界面中，我們只顯示編號為“1”的單元為：
發現節點排布的方式是“逆時針”順序。
編號為“12”的單元，在inp文件中，單元-節點的編號為：
12, 30, 9, 10
在cae界面中，我們只顯示編號為“12”的單元為：
發現節點排布的方式是“逆時針”順序。
編號為“123”的單元，在inp文件中，單元-節點的編號為：
123, 83, 34, 33
在cae界面中，我們只顯示編號為“12”的單元為：
發現節點排布的方式是“逆時針”順序。
至此我們可以得出結論，在一階三角形單元中，單元-節點的索引方式遵循“逆時針”順序，即按照節點順序彎曲右手四個手指，大拇指指向屏幕外。
展開 ABAQUS二維裂紋擴展模擬詳解
需要特別注意的是，由于裂紋尖端的網格本質上是退化的四邊形單元，因此不能劃分為三角形單元，因此在Mesh Control中需要將該區域設置為四邊形單元為主（Quad-dominated）的掃略（Sweep）劃分方式，不能指定為三角形單元（Tri）。
圖10 裂紋尖端的網格控制
對于裂紋尖端區域外圍的環形區域，主要起到單元過渡的作用，因此劃分方式選擇結構化（Structured）的四邊形（Quad）單元。對于其他區域，可采用四邊形單元或以四邊形單元為主的自由劃分方式，最終得到裂紋區域的網格控制如圖11所示。
圖11 裂紋區域網格控制
在ABAQUS的Mesh模塊中，不同顏色的區域代表不同的網格劃分方式，黃色代表掃略劃分方式，綠色代表結構化的映射劃分方式，粉色代表自由劃分方式。為了保證裂紋區域的網格正常劃分，還需要通過Seed Edges工具調整裂紋區域的網格數量，如圖12所示。
圖12 裂紋區域網格布種
最終可以得到裂紋區域的網格如圖13所示。
圖13 裂紋區域網格
再次注意，雖然通過前面的劃分方法最終在裂紋尖端區域得到了三角形的網格，但其實際上是退化的四邊形單元，三角形單元是不能用于圍線積分計算的。這一點從單元類型的選擇也可以看出，本分析中假設中心裂紋板處于平面應力狀態，采用帶有縮減積分的8節點四邊形單元CPS8R進行計算，雖然此時ABAQUS自動檢測處模型中存在三角形單元，并推薦使用6節點的平面應力單元CPS6M，但從后面的分析可知，該單元類型實際上并沒有被用于裂紋尖端區域。
展開 【JY】Abaqus殼單元概述與應用（一）
數值積分在沿殼厚度方向的一系列積分點上進行，這些積分點就是單元變量可以被輸出的位置。最外層的積分點位于殼單元的表面。
復合材料鋪層定義：對于復合材料層合殼，殼截面定義變得更加復雜，需要定義各鋪層的材料、厚度、方向和順序。在 Abaqus/CAE 中，可以通過 "Composite Layup" 工具方便地定義復合材料鋪層。常規殼方法將層合板建模為一系列薄層，同時忽略厚度方向的復雜剪切變形；連續殼方法則使用三維網格離散化，能夠更準確地描述厚度方向的剪切變形。
4.2 殼單元網格劃分最佳實踐
網格密度控制：殼單元的網格密度應根據結構的幾何特征、載荷條件和分析目標確定。在應力集中區域、幾何突變處或載荷復雜的區域，應使用更精細的網格。對于以彎曲為主的結構，沿厚度方向應至少劃分 2-3 個單元，以捕捉厚度方向的應力梯度。對于高精度分析，可能需要更多的單元。
單元形狀優化：殼單元應盡量保持規則的形狀，如接近正方形的四邊形單元。應避免使用高度扭曲或長寬比過大的單元，這些單元可能導致計算精度下降或收斂困難。對于三角形單元，應盡量使其接近等邊三角形，避免出現尖銳角度。
混合單元類型的使用：在同一模型中可以混合使用四邊形和三角形單元，但應注意避免在關鍵區域使用三角形單元，因為它們通常精度較低。在幾何復雜區域，可以使用三角形單元過渡，但應盡量限制其使用范圍。
結構化網格的優勢：對于規則幾何區域，應優先使用結構化網格劃分技術，生成高質量的四邊形網格。結構化網格不僅質量高，而且計算效率更高，能夠提供更準確的結果。在可能的情況下，應避免使用自由網格劃分，特別是對于高精度分析。
4.3 殼單元結果評估與解釋
應力結果位置：殼單元的應力結果通常在殼的頂面、中面和底面的積分點處計算，然后外推到節點。
展開 有限單元法重要知識點 附有限單元法原理與應用朱伯芳第三版下載
三角形單元中，三結點的常應變單元?其單元內應力是常量，它是一種簡單但精度低的單元：六結點的二次三角形單元精度高但不能適應曲線邊界。而矩形單元，其精度雖比相應的三角 形單元高，但不易改變單元尺寸，以及不能適應曲線邊界和非直角的直線邊界。平面等參數 單元適應了曲線邊界和非直角的直線邊界。
四、 有限元網格劃分的過程中應注意哪些問題？
1. 網格數目
網格數目的多少將影響計算結果的精度和計算規模的人小。一般來講，網格數目增加，計算 精度會仃所進步，但同時計算規模也會增加。實際應用時可以比較兩種網格劃分的計算結 果，假如兩次計算結果相差較大，可以繼續增加網格，相反則停止計算。
2. 網怡疏密
網格疏密是指在結構不同部位采用人小不同的網格，這是為了適應計算數據的分布特點。在 計算數據變化梯度較人的部位（如應力集中處），為了較好地反映數據變化規律，需要采用比 較密集的網格。而在計算數據變化梯度較小的部位，為減小模型規模，則應劃分相對稀疏的 網格。
3. 單元階次
選用高階單元可進步計算精度，所以當結構外形不規則、應力分布或變形很復雜時可以選 用高階單元.但高階單元的節點數較多，在網格數目相同的情況下由高階單元組成的模型 規模要大得多，因此在使用時應權衡考虎計算精度和時間。
4.
展開 
學習有限元需了解的知識點
26、 假設平面三節點三角形單元的的位移模式為：
U=a1x2+a2xy+a3y2
V=a4x2+a5xy+a6y2
試計算該單元的形函數矩陣、單元剛度矩陣，并討論該單元的特性。
27、 平面矩形單元的位移、應力在單元內、單元邊界上有何特性？試說明矩形單元剛度矩陣的計算與坐標原點位置無關。
答：常數項和線性項的系數反映了單元的剛體位移和常應變，滿足收斂性的必要條件；在單元邊界上，由于u，v分別僅為x或y的線性函數，則這樣的單元的位移函數是雙線性函數，這說明單元邊界上的兩點能唯一確定變形后的邊界，而對于相鄰的單元公共邊界，它們具有公共節點，則不論按哪個單元確定公共邊界上的位移，都能保證公共邊界上具有相同的位移，即單元邊界處位移具有連續性，滿足協調性要求。
28、 何謂面積坐標？其特點是什么？
答：Li=Ai/A；Lj=Aj/A;Lm=Am/A特點：只有兩個坐標是獨立的：Ai+Aj+Am=1
29、 試分析以下幾種平面單元的位移在單元公共邊界上的連續性：1）常應變三角形單元；2）四節點矩形單元；3）六節點三角形單元；4）四節點直線邊界四邊形等參單元；
5）八節點曲線邊界四邊形等參單元。
答：常應變三角形單元：形函數只與節點坐標有關；單元應變分量均為常量；
收斂性：位移函數含單元常量應變；反應單元剛體位移；單元內部位移連續；相鄰公共邊界連續協調。
四節點矩形單元：位移函數滿足收斂性條件，為協調單元；較常應變單元有更高的計算精度。
六節點三角形單元：比常應變三角形單元精度高
30、 非節點載荷等效的基本原則是什么？
答：能量等效原則和圣維南原理。
31、 試計算三節點三角形邊界上不同線性分布載荷的等效節點載荷。
展開 設計仿真 | 聯合仿真助力灣流航空機翼建模求解時間提升50%
灣流工程師最初使用三角形單元用于過渡語的網格，單元的拆分需要大量的手動調節，并且由于三角形單元本身的局限性也降低了求解的準確性。
圖 1: 通過三角形單元由粗四邊形網格向精細四邊形網格過渡
新機型的設計上每個控制面上都有大量的機構運動副構成，通常是方向舵、左右舷升降舵、左右舷副翼，以及所有六個擾流板。方向舵有四個鉸鏈位置和兩個驅動器，副翼有兩個驅動器和四個鉸鏈，升降舵有兩個驅動器和五個鉸鏈，擾流板有一個驅動器和四個鉸鏈。這些構件的建模和分析需要在一到兩個月內完成，以滿足典型飛機計劃的調度要求。
灣流公司工程師Brian Smith表示:“建模過程中涉及的大量手工調整，需要一個由六人組成的團隊來完成大部分工作。前一種方法的另一個問題是三角形單元的使用。通常，四邊形單元比三角形單元更受歡迎。三角形單元通常不如四邊形單元準確。多數情況下使用四邊形單元居多。
解決方案
灣流現在使用Patran中嵌入的Rutman實用程序對緊固件進行自動化建模。
展開 結構有限元計算中的Locking和F-Bar單元
比如說在做熱-變形耦合分析時，溫度（在三角形單元中線性分布）和應變（在三角形單元中均勻分布）計算精度帶來的誤差。
對于該問題的進一步解釋建議參看文獻[2]，在哪里給出了通俗和準確的解釋。
1.2 如何減輕鎖死
可以通過單元技術來減輕鎖死的影響。如Selective reduce integration, Incompatible，Assumed strain，Mixed formulation，MITC等等。對于近似不可壓縮材料的體積鎖死，最常見的對應方法是采用B-bar單元，這種單元的新式變種是F-bar單元，這是本文下面要講的內容。
2. B-bar單元[3,4]
一般的，單元內任意一點的應變ε表為單元節點i的位移u的函數。在B-bar單元中，將應變分解，其中為單元的平均體積應變或單元中心的體積應變值為偏差應變值。相應地，單元內任意一點的應變用下式計算
這里將體積應變做平均計算的方法緩和了體積鎖死。這種方法實現起來非常簡單，計算費用不高，效果也不錯，因此得到非常廣泛的應用。但是這種方法視乎不能滿足LBB條件[5]。
3. F-bar單元[6,7]
相應于B-bar單元對應變的加法分解，單元中的任意一點的變形梯度（deformation gradient）可以分解為體積變形梯度和偏差變形梯度兩個部分 F=Fs Fv。 理論上來說這種分解方式是嚴密的，而B-bar單元中的應變加法分解是一種近似，只是在應變很小時與變形梯度的乘法分解相近。因此把B-bar單元推廣到 F-bar單元是非常自然的，在變形較大時采用F-bar單元精度更好。
在F-bar單元中，單元中的任意一點的變形梯度表述為，J=detF. 其中
可以是單元中Jocobian J的平均值 ，也可以指單元中心的Jocobian值。
展開 設計仿真 | 聯合仿真助力灣流航空機翼建模求解時間提升50%
灣流工程師最初使用三角形單元用于過渡語的網格，單元的拆分需要大量的手動調節，并且由于三角形單元本身的局限性也降低了求解的準確性。
圖 1: 通過三角形單元由粗四邊形網格向精細四邊形網格過渡
新機型的設計上每個控制面上都有大量的機構運動副構成，通常是方向舵、左右舷升降舵、左右舷副翼，以及所有六個擾流板。方向舵有四個鉸鏈位置和兩個驅動器，副翼有兩個驅動器和四個鉸鏈，升降舵有兩個驅動器和五個鉸鏈，擾流板有一個驅動器和四個鉸鏈。這些構件的建模和分析需要在一到兩個月內完成，以滿足典型飛機計劃的調度要求。
灣流公司工程師Brian Smith表示:“建模過程中涉及的大量手工調整，需要一個由六人組成的團隊來完成大部分工作。前一種方法的另一個問題是三角形單元的使用。通常，四邊形單元比三角形單元更受歡迎。三角形單元通常不如四邊形單元準確。多數情況下使用四邊形單元居多。
解決方案
灣流現在使用Patran中嵌入的Rutman實用程序對緊固件進行自動化建模。
展開 