基于等幾何分析的自由曲面殼體包容式節點力學分析方法*

孫遠韜 王斌賀 張 氫 秦仙蓉

同濟大學機械與能源工程學院 上海 201804

0 引言

桁架結構由于其結構簡單、承載能力強、加工方便且自重輕的優勢在起重機械中得到了廣泛的應用。節點的力學性能對桁架結構的承載能力和受力起到了至關重要的作用。但目前常見的節點形式存在應力集中，抗疲勞能力差等問題。故一種新型管相貫節點——包容式節點被提出應用[1]，如圖1所示。

圖1 包容式節點

目前包容式節點已在某風電安裝平臺、起重機、RMG等項目中得到初步應用，取得良好效果，展示出較好的力學承載能力和抗疲勞性能。但目前包容式節點的設計仍處于經驗設計階段，其曲面形狀具有較大隨機性，限制了包容式節點的進一步發展。為進一步發揮包容式節點的優勢，將曲面形狀采用自由曲面的形式進行描述并進行靈敏度分析、形狀優化等研究必不可少。

常規的有限元手段在進行自由曲面殼體形狀優化時存在容易出現網格畸變，在優化過程中需不斷重復劃分網格，無法獲取靈敏度的解析解等問題，限制有限元方法在包容式節點曲面優化設計中的應用。而由Hughes等于2005年提出的等幾何分析（isogeometric analysis，IGA）[2]是一種采用幾何模型自身的NURBS函數作為形函數的有限元分析方法。其將CAD、CAE一體化，在優化過程中無需劃分網格，不因對幾何模型的離散造成精度損失，能給出自由曲面靈敏度的全解析解，且優化結果較為光滑，適用于進一步對包容式節點自由曲面的優化設計。

因此，本文基于對包容式節點曲面優化設計的需求，研究等幾何分析在殼體自由曲面分析中的應用，并建立包容式節點的曲面模型，采用等幾何分析方法對不同厚度的包容式節點自由曲面進行力學分析，研究不同厚徑比對等幾何分析結果的影響，并采用有限元仿真結果進行驗證，為接下來對包容式節點曲面的優化設計提供力學分析手段。

目前，國內外對于等幾何分析在殼體分析領域的研究包括：Kiendl.J等[3]基于Kirchhoff-Love假設建立了三自由度的等幾何殼單元，對相關理論進行推導后建立一個基于CAD程序的等幾何分析程序驗證該方案的有效性。Nguyen[4]將T樣條曲線應用到等幾何分析方法中，并采用實例說明該方案允許局部細化。張升剛等[5,6]研究了Reissner-Mindlin退化殼理論在等幾何分析中的應用，并研究了基于殼體等幾何形狀靈敏度的曲面優化。李新康等[7]進一步研究了基于Kirchhoff-Love理論的等幾何分析在薄殼靜態問題中應用，并對NURBS單元的階次與數量與該方法精度之間的聯系進行了分析。常峰[8]分析不同罰函數對Kirchhoff-Love板殼等幾何分析方法的精度影響，并將等幾何分析理論應用到斷裂力學當中。袁沛等[9]基于Kirchhoff-Love理論建立了非結構化T樣條薄殼單元，并對其在模態計算和彈性變形分析中的應用進行了研究。

1 NURBS方法基本原理

NURBS方法是一種特殊的樣條曲線，起源于航天領域對于曲面設計的要求。其從數學形式上可寫為控制點的插值分段有理多項式形式，與常見的拉格朗日等插值方式不同的是，NURBS曲線并非經過控制點，其控制點落在擬合曲線外側，即凸包性。同時，一段NURBS曲線形狀僅與k+1個控制點相關(k為NURBS曲線的階次)，即局部性。幾何不變性、參數連續性等也是NURBS方法被廣泛應用于描述曲線曲面形狀的重要原因[10]。

NURBS曲線通常可寫為以下分段有理多項式矢函數形式

式中：Ni,k(u)為k次B樣條規范基，其具體數值一般根據節點矢量U=[u0，u1，…un+k+1]由de-door遞推公式獲得為

此外，di(i=0，1，…n)為控制頂點，wi(i=0，1，…n)為控制點頂點對應的權因子。

NURBS曲線的表現形式如圖2所示。

圖2 NURBS曲線圖

但在程序中采用de-door計算NURBS的基函數需要大量運算，因此，在實際運算過程中可以通過公式推導將NURBS曲線寫為矩陣表示形式。同樣條件下，采用矩陣表示形式可以提高運算速率20倍以上。

NURBS曲面一般可寫為

式(3)、式(4)中：dij(i=0，1，…m；j=0，1，…n)為拓撲矩形陣列，形成一個控制網格。wij(i=0，1，…m；j=0，1，…n)為其對應的權因子。Ni,k(u)(i=0，1，…m)和Nj,l(u)(j=0，1，…m)分別為u向k次B樣條規范基和v向l次B樣條規范基，分別由節點矢量U=[u0，u1，…um+k+1]，V=[v0，v1，…vn+l+1]經de-door遞推公式獲得。

NURBS曲面的表現形式如圖3所示。綠色曲面為NURBS曲面，紅色星形點為NURBS曲面對應的控制點。事實上，由于幾何連續性以及精度方面的優勢，工程中自由曲面一般采用NURBS方法進行描述，并于1991年NURBS方法被國際標準化組織定義為描述工業產品幾何形狀的唯一數學方法。

圖3 雙三次5×5控制點曲面

2 基于Reissner-Mindlin退化殼理論的等幾何分析

2.1 基礎理論

Reissner-Mindlin退化殼單元具有5個節點自由度，u,v和w以及繞u,v兩軸的轉動自由度，相較于Reissner-Mindlin單元去除了uv平面內的旋轉自由度，Beason[11]采用了大量實例證明了將Reissner-Mindlin退化殼單元應用于等幾何分析的可靠性。曲面殼體坐標系的定義如圖4所示。

圖4 曲面殼體坐標系

等幾何分析采用描述曲面自身的NURBS函數作為插值形函數，故有殼體上任意一點坐標的表示為

式中：xi(u,v,ζ)為殼體內一點坐標的第i個維度，為控制點對應的法向量，Hk(ζ)為殼體法向插值形函數，Rk(u,v)為NURBS曲面的k次有理基函數。

任意一點的位移為

式中：kδ為第k個控制點對應的位移，為控制點對應的切向量，則為控制點繞單位向量的轉角。

由此可得，等幾何分析對應的應力應變矩陣為

2.2 等幾何分析剛度矩陣的計算

有限元分析的單元剛度矩陣一般為

等幾何分析方法采用等參元思想，其Reissner-Mindlin殼單元的單元剛度矩陣為

式中：[T]為轉置矩陣，將單元剛度矩陣從局部坐標系轉置到總體剛度矩陣的坐標系下。

其轉置矩陣T為

式中：lij為積分點對應的切向量或法向量Vi相對于整體坐標系下的3個坐標軸的方向余弦。

3 等幾何分析在包容式節點中的應用

3.1 等幾何分析的程序實現

基于Reissner-Mindlin退化殼理論的包容式節點自由曲面等幾何分析由以下幾個步驟組成：

1）基于NURBS方法建立包容式節點的自由曲面模型。一般可通過逼近和插值獲取包容式節點對應的NURBS曲面模型，也可直接設定NURBS參數建模對應自由曲面模型。

2）輸入NURBS曲面模型數據 不同于有限元方法的需要根據曲面坐標通過劃分網格來劃分單元，等幾何分析可直接采用幾何模型自身NURBS網絡進行單元劃分，無需在幾何形狀改變后重新劃分網格，也不會因CAD、CAE模型之間的轉化造成誤差。它一般通過NURBS曲面自身的節點矢量劃分單元，并對單元進行編號，記錄每個單元對應的控制點編號。計算控制點對應的曲面最近點(u,v)坐標，計算該坐標的法向量作為控制點對應的法向量。賦予所有單元計算單元剛度矩陣所需數據。如圖5所示等幾何分析單元1以及其對應的控制點。

圖5 等幾何分析單元與控制點對應圖

3）遍歷所有單元 在單元參數域中取高斯積分點并計算對應的R基函數的值以及其對于參數u，v的導數、雅克比矩陣、應力應變矩陣、彈性矩陣、轉置矩陣等，通過高斯積分公式進行積分運算獲取單元剛度矩陣。

4）組裝總體剛度矩陣 遍歷所有單元，根據控制單在單元剛度矩陣以及總體剛度矩陣中的對應位置，將單元剛度矩陣組裝到總體剛度矩陣之中。

5）施加載荷 根據載荷的類型及大小生成對應的總體載荷矩陣。

6）施加約束 根據所約束控制點的自由度對總體剛度矩陣進行對角元素置一法處理。

7）采用稀疏矩陣運算庫對KU=F進行求解，所獲得的位移向量U即為各個控制點對應的位移值。通過NURBS函數的R基函數進行插值處理便可獲得曲面上任意一點的位移值。

3.2 包容式節點的等幾何分析

包容式節點作為一種曲面殼體結構，一般由凸包和半圓柱面組成，如圖6所示。

圖6 包容式節點實物圖片

采用曲面擬合的方式建立起包容式節點的NURBS曲面模型，為方便展示結果，采用Python對NURBS曲面表面均布采點可視化，如圖7所示。

圖7 包容式節點NURBS曲面模型

所采用包容式節點NURBS曲面模型參數為：長度L=8 m，半圓柱面半徑R=2 m，凸包高度H=1 m，厚度=0.2 m。控制點數目20×20，且曲面為雙三次NURBS曲面。力學參數設置為：彈性模量E=2.06×1011Pa，泊松比ε=0.3。本質邊界條件為約束其4條邊的所有自由度，自然邊界條件為在包容式節點中部施加F=3 000 N的集中載荷。采用等幾何分析對該包容式節點模型進行分析，對NURBS曲面模型表面局部采點，并記錄位移值，為便于與有限元分析結果進行比較，采用Python將其可視化，如圖8所示。

圖8 等幾何分析仿真結果圖

采用有限元方法在相同條件進行力學分析可獲得如圖9所示結果。

圖9 有限元S仿真結果圖

如圖8、圖9所示，兩者位移分布基本一致，采用等幾何分析算出最大位移值為5.62×10-7m，采用Abaqus分析最大位移值為6.28×10-7m。兩者結果基本一致。

采用相同雙三次NURBS曲面模型，厚度分別取為0.01 m、0.1 m、0.15 m、0.2 m、0.25 m、0.3 m、0.5 m的包容式節點自由曲面模型進行試驗。通過等幾何分析以及有限元方法對其進行力學求解，可獲得其z方向最大位移值如表1所示。

表1 IGA與Abaqus結果對照表 m

由表1可知，當包容式節點厚度較大時，IGA結果與Abaqus分析結果較為接近。其與有限元分析結果的誤差如圖10所示。

圖10 誤差隨包容式節點厚度分布圖

圖10中紅色折線圖為Abaqus仿真結果位移值，藍色折線圖為等幾何分析仿真結果，黃色柱狀圖為誤差值。由圖10可知，采用基于Reissner-Mindlin退化殼理論的等幾何分析方法對半徑L=1的包容式節點進行分析時，當其厚徑比小于0.005時其結果誤差較大。而當其厚徑比大于0.05時結果與有限元方法分析結果較為接近，當其厚徑比大于0.15時，采用該方法的分析結果基本與有限元方法的分析結果一致。目前，主要應用的包容式節點厚徑比為0.2左右，故該方法適用于包容式節點的力學分析。

4 結論

本文簡要地介紹了包容式節點在工程中的應用價值以及其進一步發展所遇到的問題。說明了對包容式節點自由曲面模型進行等幾何分析方法研究的必要性，并通過建立基于Reissner-Mindlin退化殼理論的等幾何分析殼單元模型，給出了相應的單元剛度矩陣等的計算方式，將其應用到包容式節點的分析之中。通過實現等幾何分析程序并對包容式節點曲面進行力學分析，與有限元方法仿真結果進行對比驗證了該方法的有效性。并對不同厚度的包容式節點采用基于Reissner-Mindlin退化殼理論的等幾何分析方法進行分析實驗，發現在包容式節點厚徑比較小時分析結果嚴重失真，而當包容式節點厚徑比較大時則能取得良好的精度。目前所使用的包容式節點厚徑比滿足該方法的要求。