分層加密網格的懸鏈面膜結構數值找形優化

于 琪，張少欽，熊 磊，李云鵬

（南昌航空大學 土木建筑學院，南昌 330063）

引 言

膜結構作為一種柔性結構，其膜材不能抗壓、抗彎，故需要在符合給定邊界條件的情況下，通過施加一定的初始預應力使膜面產生只有張力作用的自平衡曲面來維持結構形狀和抵抗外荷載，這種由應力分布狀態求解結構形狀的過程常被稱為找形（form finding）。從膜結構分析的三個過程（初始形態分析、荷載效應分析和裁剪分析）來看，初始形態分析（找形）是膜結構分析過程中最為活躍、最為關鍵的部分，結構初始形態的好壞直接影響后續分析過程[1]。

從最開始的物理找形到近年來的數值找形，眾多學者廣泛探索有效的找形方法及不斷加以改進。用于膜結構初始形態確定的數值方法主要有：力密度法(Force Density Method)[2]、動力松弛法(Dynamic Relaxation Method)[3]和非線性有限元法(Nonlinear Finite Element Method)[4]。在計算機軟硬件飛速發展的今天，非線性有限元法應用更為廣泛[5]。利用非線性有限元法確定膜結構初始形態的計算方法主要包括節點平衡法和支座移動法[6]。衛東等[6]在對具體薄膜結構算例進行初始形態確定的過程中，通過對各種傳統方法的分析比較，提出了綜合節點平衡法。葉繼紅等[7]根據建筑師對膜面造型的變化要求，針對動力松弛法中阻尼參數的取值，對常規動力松弛法進行了改進。周樹路等[8]針對力密度法進行膜結構找形時總需要多次試算才能確定合適的力密度取值的相關問題，提出了改進找形過程的措施。魏德敏等[5]通過計算發現，膜的預張力、單元的劃分、邊索和脊索的預張拉力值以及初始彈性模量等參數的選取都會影響找形結果?？v智育等[9]根據幾何非線性有限元理論，提出運用新的8結點曲面四邊形等參單元進行張力膜結構找形分析的具體做法。袁駟等[10]提出一種基于EEP（單元能量投影）法的二維非線性有限元自適應方法，成功將之應用于膜結構的找形分析。

針對結構找形分析過程中，采用傳統網格加密方法能提高結構找形精度但計算量較大且效率低下的問題，本文從調整有限元網格的角度，提出了一種在保證單元總量不變的情況下，通過分層加密網格的方式來實現膜結構找形優化的方法，并結合具體的懸鏈面膜結構算例驗證其有效性。

1 懸鏈面膜結構數值找形優化的實現

1.1 數值找形的實現

懸鏈面（圖1）是旋轉面中唯一有解析解的最小曲面，其曲面方程為：

圖1 懸鏈面結構模型示意圖[6]

其中，a表示懸鏈面內圓半徑，h表示懸鏈面高度，，其含義為結構內任一高度上圓的半徑，(X，Y，Z)為懸鏈面內節點坐標。h的求解方法如下式所示。

其中，b表示懸鏈面外圓半徑。

基于ANSYS軟件采用節點平衡法對懸鏈面膜結構進行數值找形，具體流程如下：

采用SHELL41三維殼體單元模擬懸鏈面膜面，建立初始形狀（可通過曲面擬合[6]等方式得到）模型，設定很小的彈性模量（如取實際彈性模量的10?3~10?6），利用式(3)[11]將膜面初始預應力轉換為下降溫度值施加于膜單元以確定其平衡曲面（稱為第一次找形）；進行模型修正，直到膜面應力趨于相等或滿足一定的應力誤差為止，得到近似的最小曲面（稱為第二次找形）[12]。

其中，ΔT為溫度下降值；σ為膜面初始預應力；E為膜材彈性模量；α為膨脹系數；t為膜材厚度。

1.2 分層加密網格的實現

分層加密網格分兩步進行。第一步，根據初始條件進行結構有限元建模及分析，求出正常網格劃分情況下的膜面數值解與解析解誤差。第二步，重新建立模型，對原結構誤差數值相近的區域進行歸類，將結構劃分為不同的網格劃分區域，在單元總量保持不變的情況下，對誤差較小區域采用較稀疏的網格，對誤差較大區域則將網格劃分的密集。這就達到了分析計算量增加不大的同時使結構找形結果得到優化的目的。

首先，全面收集區域概況資料、土壤污染源資料、土壤環境和農產品資料及其對應的圖件資料等，并通過分析，研究有關信息劃分區域內的地理單元，為污染風險評價提供基礎支撐；其次，通過評估土壤及農產品調查點位的土壤重金屬污染情況、農產品重金屬污染情況，反映該點位土壤和農產品重金屬污染狀況，劃分點位土壤重金屬風險等級，為區域風險評價提供基礎依據；最后，通過評價單元內土壤重金屬污染風險、農產品重金屬危害風險、重金屬生物可利用性等因素的情況，反映區域內耕地土壤和農產品重金屬污染情況，劃分區域內土壤重金屬風險等級，為重金屬污染耕地分級管理提供依據。

以懸鏈面1/4結構平面圖模型為例，分層加密網格法的實現過程如下：

首先建立正常網格劃分方法下的結構模型，將模型周向圓弧邊及徑向直線邊劃分為12等份，得到如圖2所示網格尺寸為12×12的結構模型。

圖2 懸鏈面1/4結構(平面展開)網格

按照圖2所示網格劃分方式求得結構數值解與解析解誤差后，依據誤差數值相近進行歸類的原則，對重新建立的結構模型進行結構網格區域劃分，如圖3a所示。其中，A、B、C、D表示不同節點數值解與解析解誤差對應區域。根據分析結果，在誤差的數量關系上從區域A至區域D誤差逐級遞增，實際結構劃分區域數視誤差情況而定。將圖3a中四個區域從A至D逐級加密網格。具體做法為：周向圓弧邊上各區域網格劃分份數相同，徑向直線邊上每個區域按等分方式劃分網格，由于每個區域徑向寬度不同（越靠近小圓，徑向寬度越?。?，此法就實現了從外圓至內圓分層加密網格的構想，網格劃分結果如圖3b所示。

圖3 分層加密網格法示意圖

2 計算結果與分析

算例取內圓半徑a=8 m，外圓半徑b=40 m，高度h=18.339 5 m。結構上下邊界均固定，初始預應力σ=2.5 MPa，彈性模量E=250 MPa，厚度t=1 mm，泊松比μ=0.34，膨脹系數α=0.01。采用兩種不同的網格尺寸（1/4結構網格尺寸分別為策略一：12 ×12和策略二：16×16），對正常網格劃分以及采用分層加密網格法下（網格調整前后）結構變形圖、應力云圖及節點數值解與解析解誤差進行對比分析。

2.1 策略一結果對比

圖4中，網格調整前后結構單元及節點總量保持不變，其中單元有1152個，節點624個?；趫D4可以看出，采用分層加密網格法下的結構網格分布更加均勻，特別是靠近小圓位置，結構曲面較前者更為平滑。

圖4 策略一結構變形云圖

從圖5可以得出，網格調整前結構應力在4.27~4.68 MPa之間，應力差值為0.41 MPa，網格調整后的結構應力分布在4.38~4.72 MPa之間，應力差值為0.34 MPa。整體應力方面，網格調整后增加了約2.5%，應力稍有提高；應力差值方面，網格調整后降低了約20.5%，結構應力分布較網格調整前提高明顯。

圖5 策略一結構應力云圖(俯視圖)

由于懸鏈面是旋轉面[13]，其在各個方向上變化相同，為了方便計算，在判斷節點數值解與解析解誤差時，文章選取的節點均在XZ軸平面內，故Y坐標取為0。此外，結構邊界上的節點未產生位移，即不存在誤差，故未在表中列出。

表1 策略一網格調整前節點數值解與解析解誤差

表2 策略一網格調整后節點數值解與解析解誤差

表3 策略一坐標相近節點誤差對比

2.2 策略二結果對比

圖6中，網格調整前后結構單元及節點個數相同，其中單元個數為2048個，節點1088個。由圖6可以看出，新方法下的結構網格分布較前者更加均勻，曲面也更為平滑。

圖6 策略二結構變形云圖

依據圖7，網格調整前結構應力在4.28~4.65 MPa之間，應力差值為0.37 MPa，網格調整后的結構應力分布在4.38~4.72 MPa之間，應力差值為0.34 MPa。整體應力方面，網格調整后增加了約2.3%，應力稍有提高；應力差值方面，網格調整后減小了約9%，結構應力分布較網格調整前也有一定減小。

圖7 策略二結構應力云圖(俯視圖)

從表4、表5、表6中可以看出，采用新方法進行結構找形的節點數值解與解析解誤差較正常網格劃分方式下的誤差仍有顯著降低。

表4 策略二網格調整前節點數值解與解析解誤差

表5 策略二網格調整后節點數值解與解析解誤差

表6 策略二坐標相近節點誤差對比

對比圖5a、圖7a、表1及表4還可以發現，分析時采用16×16的網格比12×12的網格應力分布更均勻，數值解與解析解誤差也更小。

3 結 論

1) 結構變形方面，分層加密網格法比常規方法求得的結構變形更協調；應力分布方面，采用分層加密網格法計算出的結構應力分布比常規方法更加均勻；數值解與解析解誤差方面，采用分層加密網格法計算出的誤差相比于常規方法，其計算結果有明顯改善。

2) 進行結構分析時采用較密集的網格能提高分析結果的準確度。