基于真實轉角與“試探-修正\"插值的四邊形單元

中圖分類號：TH165 文獻標志碼：A

Abstract：To addressthe problem of incompatibility between diffrent elements，the paper introduces a“trialcorection”displacement interpolation scheme.This scheme is utilized to construct a four-nodequadrilateral plane elementthat incorporates adriling degree of freedom.The“trial-correction”interpolation-based four-node quadrilateral plane element takes translational displacementsand the driling degree of freedomasnodal parameters，and higher-order interpolation functions are employed to approximate the displacement field.Firstly， bi-linear interpolation isused to trail the displacementfields.And then，according to thedeviation ofthedrilling degree offreedom，the displacement fieldsare corrcted with bi-cubic interpolation.Theconvergence of the“trialcorrection”interpolation-basedfour-node quadrilateral plane element is provedbythe patch test，and its performance is furtherverified by three examples.The numerical results show that the“trial-corrction”

interpolation-based four-nodequadrilateral plane element hasa high convergence rate，a high numericalaccuracy， and can be compatible with beam elements for convenient mix-element modeling.Moreover，the“trial-corection” interpolation method is parameterizedandhas good extendibility，which laysafoundation forfuture research about other elements with true rotation angle.

Key Words：finite element method; interpolation; convergence of numerical methods;quadrilateral plane ele. ment; true rotation angle

隨著計算機技術的發展，有限元法已被廣泛應用于解決各種科學和工程問題1.對于同時具有不同結構的復雜機構，在運用有限元法進行分析的過程中，可能會遇到奇異剛度矩陣的問題，這會導致分析無法成功[2-3].而帶有旋轉自由度的平面單元可以通過與板單元結合來構造殼體單元，從而使殼體單元避免奇異剛度矩陣問題.不過必須指出的是，開發高性能的低階帶有旋轉自由度的平面單元仍然是研究人員面臨的挑戰.

Allman4提出對節點使用二次插值構造出帶真實轉角的三角形平面單元.在此基礎之上，Cook5證明了通過線性應變三角形單元和角節點旋轉的三角形單元之間的轉換技術可以得到Allman等效元，并在此基礎上提出了一種帶轉動自由度的四邊形單元.近十年間，帶轉角的平面單元有了多方面的完善與擴充，其中Shang等6提出了一種帶有真實旋轉自由度的非對稱膜單元，包括開發了具有真實轉動自由度的非對稱三角形單元[7和非對稱四邊形單元[8].Rezaiee-Pajand等9使用泰勒級數的方法對應變場進行插值，提出了基于應變的三角形膜單元.Bouta-gouga[]提出了一種具有真鉆井旋轉自由度的基于位移的膜單元公式.面內轉角在雜交有限元、混合有限元中也扮演重要角色，例如Daszkiewicz等[1]基于Hu-Washizu原理，提出了含鉆井旋轉的混合四節點殼單元，并將該單元用于非線性分析， .Wang^[12] 提出了一種具有真實轉動自由度的平面單元.文穎等[13]構建了一種含有面內自由度的四節點板殼單元，用于非線性分析中.Sangtarash等[14]利用Airy應力函數的解析解以及四節點單元的等參位移構造了一個高性能且帶有旋轉自由度的非對稱四邊形單元.基于修正偶應力理論， Wu 等[15]采用罰函數法約束引入的旋轉量并構造了一個八節點四邊形單元，而Long等[1通過引入薄板彎曲的假設，修改偶應力彈性的三維控制方程構造了類似于非線性有限元模型性能的Trefftz板單元.文獻[17]對平面單元的發展做了比較詳細的歸納.

與傳統等參單元相比，這些包含真實轉角的單元表現出更好的性能.然而，這些單元在應力/應變貼片測試中未能嚴格滿足要求，其收斂性和數值穩定性有待進一步探討.此外，帶轉動自由度的單元在公式推導方面過于復雜，不利于在工程實際中應用.

本文將平動位移和轉動位移作為節點參數，提出了一種“試探-校正\"插值方法，用于構建參數化位移插值.基于“試探-校正\"插值方法，推導了一個具有真實旋轉自由度的四節點四邊形平面單元（four-node-quadrilateral planeelement with trial-correctioninterpolation，plane-Q4-TC單元）.plane-Q4-TC單元具有三次位移場，展現出較高的數值精度和較快的收斂速度，可以與平面梁單元兼容，便于進行混合單元建模.同時，參數化的“試探-校正\"插值方案具有良好的擴展性，可用于構建梁、殼和實體等多種單元，從而形成一個多維度兼容的有限元系統.

本文內容安排如下：第1章為plane-Q4-TC單元的“試探-校正\"插值方案的理論公式推導，并通過分片測試證明所構建單元的收斂性.第2章用數值算例來驗證plane-Q4-TC單元的性能.第3章為結論.

1理論公式

1.1節點參數

假設一個無限小平面單元的真實轉動如圖1所示，圖1（a）和圖1（b）分別指剛性單元的真實轉角以及可變形單元的真實轉角，剛性單元的真實轉角 θ_u 和 θ_v 相同，因此可用轉角 表示：

可變形單元的真實轉角 θ_u 和 θ_v 不相同，因此真實轉角 θ_z 可定義為：

包含真實轉角的四節點四邊形平面單元模型如圖2所示，圖2（a）和圖2（b）分別指物理坐標系和自然坐標系下的平面單元.每個單元節點 p_i（i=1～4） 都包含3個參數，分別是平動位移 u_?v 以及繞 z 軸旋轉的真實轉角 θ_z ，位移場 u（ξ，η） 和 v（ξ，η） 均由以上3個節點參數插值得到.其中， （ξ，η） 為自然坐標，與物理坐標之間的變換關系為：

式中： x_e 和 y_e 分別是 x 和 y 方向節點坐標矢量；N（ξ，η） 是形函數矩陣，可表達為：

1.2位移插值

本節采用\"試探-修正\"法構造位移場，假設位移場 u（ξ，η） 和 v（ξ，η） 均由試探項位移 u^t（ξ，η） 、 和修正項位移 u^c（ξ，η），v^c（ξ，η） 組成，表達式為：

式中：試探項位移場 由雙線性插值得到，表達式為：

式中： u_e 和 v_e 分別是沿著 x 軸和 y 軸的節點位移； d_e 指節點參數的位移向量，表達式為：

形函數 和 的表達式為：

由式（2）得知，試探項構造的真實旋轉角可以計算為：

式中： 的表達式可以寫成：

B_z^t（ξ，η）=[B_z，1^tB_z，2^tB_z，3^tB_z，4^t]

結合式（3），將式（11）的偏導數寫成：

式中： J（ξ，η） 是雅可比矩陣，表達式為：

對于平面單元，節點 p_i（i=1～4） 可以滿足以下表達式：

式中： ξ_i 和 η_i 是節點 p_i 的自然坐標.修正項 u^c（ξ，η） ）v^c（ξ，η） 和 .θ_z^c（ξ，η） 需要滿足以下情況：

通過三次插值的方法構造修正項位移， u^c（ξ，η） 和 v^c（ξ，η） 可表達為：

式中： P（ξ，η） 是由12個插值基函數組成的行向量，如表1所示； α 和 β 是待定系數，可表達為：

為了求解待定系數，在節點 p_i 處建立以下方程：

式中：下標 ξ，η，x 和 y 表示對應的偏導數.假設修正項 u^c（ξ，η） 和 v^c（ξ，η） 呈現剛性節點特征，即修正項在節點處不引起應變，表達為：

將式（21）代入式（15）可得到：

式中： S_i^u 和 S_i^v 分別為：

將式（22）代人式（19）和式（20），得到：

其中，

將式（24）寫成聯立方程的形式，可表達為：

其中，

求解方程（26），求出待定系數：

將式（29）代人方程（17），得到：

u^c（ξ，η）=P（ξ，η）A^-1J_eS_ud_e

v^c（ξ，η）=P（ξ，η）A^-1J_eS_νd_e

因此面內位移的形函數為：

將式（6）式（30）式（31）代人式（5）可得到：

1.3單元剛度矩陣

根據式（32），膜應變可計算為：

式中： ε_x 和 ε_y 代指膜單元的正應變； γ_xy 代指剪應變；B_m^t 和 B_m^c 分別對應試探項和修正項的應變矩陣，可定義為：

其中，關于 x 和 y 的偏導數可以計算為：

令：

B_m^s=B_m^t+B_m^c

式中： B_m^s 指膜單元的應變矩陣，膜剛度矩陣可表達為：

式中： h 和 分別為單元的厚度和積分域； c 為本構矩陣.對于平面應力問題和平面應變問題， c 可分別表示為：

或者

式中： E 為彈性模量； u 為泊松比，

1.4分片試驗

位移場的設置如表2所示，分別產生恒定應變ε_x，ε_y 和 γ 以 ε_x 為例，將節點坐標 p_i 代入位移函數，可得到節點位移：

由式（40）可得：

因此可推導出位移場：

即在應變 ε_x 為常數的情況下，plane-Q4-TC單元等價于經典等參元，在 ε_x 或 ε_y 為常數的情況下也可得到類似結論.因此，plane-Q4-TC單元能夠通過分片試驗.

2數值算例

本章主要通過4個不同的算例：懸臂梁、Mac-Neal梁、Cook問題以及plane-Q4-TC單元在工程應用中的可行性分析，來驗證“試探-修正\"法構造plane-Q4-TC單元的精度以及收斂性.

2.1懸臂梁

懸臂梁的長寬比為10，梁的左端受到固定約束，在右端A點處施加垂直向下的拉力 F_y ，具體的受力情況以及結構參數如圖3所示.該懸臂梁使用不同的有限單元網格（網格尺寸 l_e=1 ，2.5，5，10，20，25，50， 100mm ）進行離散化，在 A 點處沿 y 方向的位移結果如表3所示.

作為比較，本節使用等參元和ABAQUSS4R單元對懸臂梁案例進行靜力學分析，分析結果如表3所示，以此來驗證本研究構建的單元收斂性.3種單元的位移結果如圖4所示.此外，將ABAQUS的高階四邊形單元作為參考解繪制在圖4中來直觀體現各種單元的收斂趨勢．

由圖4可知，在相同的單元網格劃分條件下，使用“試探-修正”法構造的plane-Q4-TC單元的數值結果與參考值更加接近，比等參元和ABAQUSS4R的收斂速度更快.此外，從表3的數值結果可以得知，plane-Q4-TC單元使用 10mm 網格離散時的計算結果優于其他兩種單元使用 2.5mm 網格劃分的結果.plane-Q4-TC單元使用 2.5mm 網格離散時求解的位移值與另外兩種單元使用 1mm 網格劃分得到的位移值基本相等.

綜合上述分析，使用\"試探-修正\"位移插值法構造的plane-Q4-TC單元具有很高的數值精度，在相同的單元尺寸下，該單元更加接近參考值，而且比等參元和ABAQUSS4R單元具有更快的收斂速度.

2.2MacNeal梁

圖5所示為MacNeal提出的細長懸臂梁模型，其通常用來檢驗單元剪切自鎖現象以及非常規網格對計算精度的影響，本次算例應用了兩種梯形網格[圖5（a）和圖5（b）]來做本次測試.懸臂梁長度 L= 6mm ，高度 h=0.2mm ，厚度 t=0.1mm ，楊氏模量 E= 107MPa ，泊松比 ν=0.3. 在載荷 M 作用下右端撓度的參考值為- -0.0054mm ;在載荷 P 作用下右端撓度的參考值為 0.1081mm^[6] 將計算結果歸一化并記錄在表4中.

根據表4可以看出等參元的精度在處理MacNeal梁時出現剪切自鎖現象，主要表現為精度欠佳，在不同測試環境下最優精度僅為 3.1% 為解決該問題，眾多學者在理論上進行了不斷完善，Shang等提出的US-Q40單元在該問題上甚至無限趨于解析解.

盡管本文所提出的plane-Q4-TC單元在精度上不及US-Q40單元、HS-A7單元等，但其精度依然超過了 86%

2.3 Cook問題

為了驗證plane-Q4-TC單元的性能，本節對Cook問題進行求解，即求解非規則幾何形狀的懸臂梁在端部受剪問題.Cook懸臂梁示意圖如圖6所示，梁的左端被固定，而右端則受到面內剪力的作用.通過3項指標對plane-Q4-TC單元的性能進行了評估： C 點撓度的絕對值 點的最大主應力 （σ_Amax） 以及 B 點的最小主應力（ （σ_Bmin） .隨后，將所得結果與其他學者的研究成果[19-20]進行了對比分析，相關數據如表5所示.因為plane-Q4-TC單元的應力和應變是不連續的，所以表5中平面單元的每種情況均有3個應變和應力結果： S₁ 指左側單元的結果， S₂ 指右側單元的結果， S_m 指平均值.為了直觀地比較不同方法的收斂性，分別繪制了相應指標的收斂曲線，分別在圖7＼～圖9中展示.

圖7為 C 點處的位移收斂曲線，可以觀察到，plane-Q4-TC單元具有良好的精度，比 Allman^[4] 、Q4S^[21] 和Pimpinelli[22]提出的單元具有更好的性能，但略遜于 GQ12M8^[16] 、Zhang等[19]和Choi等[23]提出的單元.對于 16×16 的網格，只有3個對比數據，plane一Q4-TC平面單元的性能比 Q4S^[21] 好，但不及Choi等[23]提出的單元.

圖8為A點處的最大主應力的收斂曲線，可以觀察到 S₁ 的收斂曲線呈自上而下的趨勢，而 S₂ 的收斂曲線則與其他學者提出的單元趨勢一致，均為自下而上收斂.此外， S₂ 的精度高于 S₁. 對于較粗的網格，plane-Q4-TC單元的性能大致與Allman單元和Q4S單元相同，但不如Zhang等提出的單元.對于16×16 的網格，plane-Q4-TC單元的性能與參考值一致，優于Q4S單元和Choi等提出的單元.

圖9為 B 點處的最小主應力的收斂曲線，就最小主應力而言，總體上 S_ν 展現出比 S₂ 更優的性能.對于包含 2×2 和 4×4 的網格規模，除了 GQ12M8^[16] 單元外， S₁ 比其他單元更接近參考值.對于細網格， S₁ 比Choi等[23]提出的單元和Q4S單元[21]收斂性更好，但稍遜于Zhang等[19]和Allman[4]提出的單元.

上述分析結果充分表明，即使與諸如混合有限元這些先進的單元理論相比，plane-Q4-TC單元也展現出了有競爭力的收斂速度.

2.4工程應用

剛柔耦合定位平臺物理模型如圖10所示，主要由3部分組成：剛性框架、柔性鉸鏈和工作平臺.柔性鉸鏈的兩端分別與工作平臺和剛性框架連接.該平臺借助柔性鉸鏈的彈性變形來進行微小運動，從而補償摩擦死區實現納米級定位，能夠滿足半導體行業的精密定位需求.結構設計對于平臺的精度具有顯著影響，而力學分析的作用不可或缺.因此為了提高剛柔耦合定位平臺在結構設計上的準確性以提升分析的效率，本節以剛柔耦合定位平臺為應用對象，驗證平面plane-Q4-TC單元在工程應用中的實用性.

根據剛柔耦合定位平臺的結構特點，提取其主視圖截面作為本節的簡單分析模型，其二維簡化模型如圖11所示.剛性框架的外側施加固定約束，工作臺中心施加驅動力 F_y

將剛性框架和工作臺使用不同單元尺寸（1，2，3，5，8mm ）的四邊形網格離散，然后采用plane-Q4-TC單元建立有限元模型.柔性鉸鏈使用梁單元建模，plane-Q4-TC單元能夠與梁單元直接耦合，無須額外添加約束.作為比較，用經典等參元和ABAQUSS4R單元代替plane-Q4-TC單元進行分析.為了避免出現兼容性問題，經典等參元和ABAQUSS4R單元需要在連接節點處添加一些梁單元作為額外約束，具體的添加位置如圖11所示.

為了評估3種單元的性能，用其求解平臺的載荷點豎直方向位移，結果如表6所示.此外，圖12描繪了隨著網格的細化，3種單元求解的收斂趨勢，以便更加直觀地比較不同單元的性能.同時，為確保分析結果的準確性，將整個二維模型（包括柔性鉸鏈）用 0.5mm 的四節點四邊形網格進行離散，并使用ABAQUSS4R單元進行建模，旨在避免平面單元和梁單元之間的不兼容問題.然后將分析的結果作為

從圖12中可以看出，隨著網格尺寸的減小，plane-Q4-TC單元可以穩定向參考解趨近.與此相反，另外兩種單元的計算結果對單元尺寸變化的敏感度較低，其變化趨勢呈現出不規則性，表明含轉角的平面單元具備更好的性能.此外，網格大小為1mm 時，plane-Q4-TC單元與參考解之間的誤差為0.54% ，而ABAQUSS4R單元和等參元的計算誤差達到 25.08% 和 24.57% ，表明plane-Q4-TC單元具有更高的計算精度.由于ABAQUSS4R單元和等參元需要額外添加梁單元約束才能避免兼容問題，這種處理實質上是一種粗糙的近似，因此這種近似不可避免地對分析結果造成一定的誤差.

3結論

本文基于“試探-修正\"位移插值方案推導一種plane-Q4-TC單元.通過數值算例表明，該單元具有較快的收斂速度以及較高的數值精度.該單元還能夠與平面梁單元直接耦合建模，無需額外的約束條件來確保兼容性，適用于工程實際.此外，“試探-修正”插值法可參數化，且具有良好的可擴展性，可作為推導包含真實轉角的梁單元、殼單元和實體單元的理論基礎，為相關研究提供參考.

最后，本文單元只涉及位移插值的改進，理論上邊界條件和非線性材料本構模型不影響位移插值的有效性.此外，本文單元所使用的真實轉角是基于小變形假設的位移梯度張量定義，因此不適用于處理幾何非線性問題.今后的工作主要將小變形單元與多體系統動力學的浮動坐標法結合，也可對材料非線性和復雜問題邊界進行擴展.

參考文獻

[1]HUANG YQ，YANGYF，WANG JZ，et al.Unsymmetric extensionsofWilson’sincompatible four-nodequadrilateral and eight-node hexahedral elements[J]. International Journal for Numerical Methods in Engineering，2022，123（1）：101-127.

[2]KARPIK A，COSCO F，MUNDO D.Higher-order hexahedral finite elements for structural dynamics：a comparative review[J]. Machines，2023，11（3）：326.

[3]HUANG GX，LI HX，LUYH，et al. Hexahedral solid element withrotational degrees offreedom based on a novel trailcorrectiondisplacement interpolation scheme[J].Iranian Journal ofScience and Technology，TransactionsofMechanical Engineering，2024，48（4）：1717-1730.

[4]ALLMAN DJ.A quadrilateral finite element including vertex rotations for plane elasticity analysis[J].International Journal for Numerical Methods in Engineering，1988，26（3）：717-730.

[5]COOK R D.On the allman triangle and a related quadrilateral element[J].Computersamp; Structures，1986，22（6）：1065-1067.

[6]SHANGY，OUYANG WG.4-node unsymmetric quadrilateral membrane element with drilling DOFs insensitive to severe meshdistortion［J].International Journal for Numerical Methodsin Engineering，2018，113（10）：1589-1606.

[7]SHANG Y，CEN S，QIAN Z H，et al.High-performance unsymmetric 3-node triangular membrane element with drilling DOFs can correctly undertake in-plane moments［J]. Engineering Computations，2018，35（7）：2543-2556.

[8]SHANGY，QIAN ZH，CEN S，et al.A simple unsymmetric 4- node 12-DOF membrane element for the modified couple stress theory[J].International Journal for Numerical Methodsin Engineering，2019，119（9）：807-825.

[9]REZAIEE-PAJAND M，GHARAEI-MOGHADDAM N，RAME ZANI M. Two triangular membrane elementsbased on strain[J]. International Journal ofApplied Mechanics，2O19，11（1）： 1950010.

[10]BOUTAGOUGA D.A formulation of membrane finite elements with true driling rotation[J].Engineering Computations，2019， 37（1）：203-236.

[11]DASZKIEWICZ K，WITKOWSKI W，BURZYNSKI S，et al. Robust four-nodeelements based on Hu-Washizu principle for nonlinear analysis of Cosserat shells[J].Continuum Mechanics and Thermodynamics，2019，31（6）：1757-1784.

[12]WANG AP.A quadrilateral membrane hybrid stress element with drilling degrees of freedom[J].Acta Mechanica Sinica，2012， 28（5）：1367-1373.

[13］文穎，戴公連，曾慶元.基于帶面內轉角自由度四節點平板殼 單元的板殼非線性分析[J].中南大學學報（自然科學版）， 2013，44（4）：1525-1531. WENY，DAIGL，ZENGQY.4-node flat shell element with drilling degrees of freedom for nonlinear analysisof plates and shells[J].Journal of Central South University（Science and Technology），2013，44（4）：1525-1531.（in Chinese）

[14]SANGTARASH H，ARAB HG，SOHRABI M R，et al.A highperformance four-node flat shell element with drilling degrees of freedom[J].Engineering with Computers，2021，37（4）：2837- 2852.

[15]WU H P，SHANG Y，CEN S，et al．Penalty CO 8-node quadrilateral and 2O-node hexahedral elements for consistent couple stress elasticity based on the unsymmetric finite element method[J].Engineering Analysis with Boundary Elements，2023， 147：302-319.

[16]LONG YQ，XU Y.Generalized conforming quadrilateral membrane element withvertex rigid rotational freedom[J]. Computersamp; Structures，1994，52（4）： 749-755.

[17]BOUTAGOUGA D.A review on membrane finite elements with drilling degree of freedom[J].Archives of Computational Methods in Engineering，2021，28（4）：3049-3065.

[18]REZAIEE-PAJAND M， KARKON M. An effective membrane element based on analytical solution［J].European Journal of Mechanics-A/Solids，2013，39：268-279.

[19]ZHANG H X，KUANGJS. Eight-node membrane element with drillingdegreesoffreedom foranalysisofin-plane stiffnessof thick floor plates[J]. International Journal forNumerical Methods in Engineering，2008，76（13）：2117-2136.

[20]BERGANPG，FELIPPA CA.A triangular membrane element with rotational degrees of freedom[J].ComputerMethodsin Applied Mechanicsand Engineering，1985，50（1）：25-69.

[21]MACNEAL R H，HARDER R L.A refined four-noded membrane element with rotational degrees of freedom[J].Computersamp; Structures，1988，28（1）：75-84.

[22]PIMPINELLI G.An assumed strain quadrilateral element with drillingdegrees of freedom[J].Finite Elementsin Analysis and Design，2004，41（3）：267-283.

[23]CHOI N，CHOO Y S，LEE B C.A hybrid Trefftz plane elasticity element with drilling degrees of freedom[J].Computer Methods in Applied Mechanics and Engineering，2006，195（33-36）： 4095-4105.