劉 超，詹海鵬

（同濟大學(xué)土木工程學(xué)院，上海 200092）

索網(wǎng)是由相互正交、曲率相反的2組拉索直接疊交而形成的一種負高斯曲率曲面結(jié)構(gòu)［1］。索網(wǎng)的拉索一般采用由高強鋼絲組成的鋼絞線，通過拉索的軸向力抵抗外荷載作用，當(dāng)拉索采用高強度材料時可以大幅減輕自重。索網(wǎng)結(jié)構(gòu)能夠比較經(jīng)濟地實現(xiàn)較大的跨度，目前已經(jīng)成為大跨建筑的主要結(jié)構(gòu)形式之一［2］。

索網(wǎng)結(jié)構(gòu)不同于一般的剛性結(jié)構(gòu)，施加預(yù)應(yīng)力之前，結(jié)構(gòu)形狀是不穩(wěn)定的，只有在適當(dāng)?shù)念A(yù)應(yīng)力作用下，結(jié)構(gòu)才具有一定的剛度。剛性結(jié)構(gòu)的一般分析方法是在已知結(jié)構(gòu)形狀的基礎(chǔ)上進行的，而索網(wǎng)這種柔性結(jié)構(gòu)，在對其結(jié)構(gòu)受力分析之前，必須確定一定的邊界和預(yù)應(yīng)力條件下對應(yīng)的結(jié)構(gòu)形狀，如何確定索網(wǎng)的平衡形狀就是索網(wǎng)結(jié)構(gòu)的找形問題。在設(shè)計實踐中，結(jié)構(gòu)的邊界條件通常是已知的，需要通過特定的找形方法獲得結(jié)構(gòu)形狀。早期，通過懸掛相應(yīng)的結(jié)構(gòu)模型或制造模型等物理方法尋找空間結(jié)構(gòu)的形狀［3-5］。隨著計算機技術(shù)的發(fā)展，針對形狀優(yōu)化的研究主要集中于數(shù)值計算方法，力密度法［6-12］、動態(tài)松弛法［13-15］和有限元法［16］已經(jīng)廣泛應(yīng)用于空間索網(wǎng)結(jié)構(gòu)的找形。

1960年，Clough［16］提出了有限元法，而后各國學(xué)者均對有限元法進行了完善，目前該方法已經(jīng)成為結(jié)構(gòu)分析的有力工具。在索網(wǎng)結(jié)構(gòu)找形問題中，有限元法應(yīng)用較為廣泛。有限元法大多直接采用已有的有限單元進行整體結(jié)構(gòu)分析，包括兩節(jié)點直線單元法、內(nèi)插多節(jié)點單元法、樣條函數(shù)單元法和兩節(jié)點曲線單元法等［17-18］。有限元法的主要思路是根據(jù)平衡方程或能量原理推導(dǎo)剛度矩陣（彈性剛度矩陣與幾何剛度矩陣），再計算力與位移的關(guān)系，最后集成整體平衡方程進行迭代計算。由于拉索的大位移特性，有限元法得到的平衡方程是非線性方程，迭代過程可能出現(xiàn)不收斂問題，因此需要不斷調(diào)整初始線形和迭代參數(shù)。1965年動態(tài)松弛法被提出，該方法可用于索網(wǎng)找形分析。利用虛設(shè)的節(jié)點質(zhì)量和運動阻尼，對節(jié)點施加激振力使之圍繞平衡位置振動，整個體系經(jīng)過平衡位置的瞬間動能最大、勢能最小。因此，將平衡位置初速度設(shè)為零，這樣整個體系的節(jié)點逐漸趨向于平衡位置，直到體系的動能小到可以忽略為止。動態(tài)松弛法的本質(zhì)是利用特殊條件下的動力方法解決靜力問題，與有限元法類似，在找形過程中同樣會出現(xiàn)不收斂的情況。Schek［6］于1974年提出了用于索網(wǎng)找形分析的力密度法，并證明了力密度法可以求解無荷載索網(wǎng)的最小總長度。力密度法由于其簡易的計算理念而得到快速發(fā)展，廣泛應(yīng)用于索網(wǎng)和膜結(jié)構(gòu)設(shè)計中。王春江等［19］指出了力密度法矩陣組裝的特點，編寫了一維變帶寬存儲規(guī)則及其實現(xiàn)算法，并用C++語言編制了相應(yīng)的計算程序。廖理等［20］提出了一種膜曲面離散的方法并將其作為索網(wǎng)劃分初始網(wǎng)格的方法。韓大建等［21］構(gòu)造了張拉索膜結(jié)構(gòu)中的T單元來協(xié)調(diào)索和膜的變形，完善了索膜結(jié)構(gòu)找形分析的力密度法。力密度法將幾何非線性問題轉(zhuǎn)化為線性問題，求解的是線性方程組，避免了非線性方程坐標初值不理想而導(dǎo)致迭代不收斂的問題，是求解空間結(jié)構(gòu)找形問題的有效方法。

基于索網(wǎng)結(jié)構(gòu)幾何形狀由力所決定的特點，提出了歐拉坐標找形法。利用靜力平衡和幾何協(xié)調(diào)關(guān)系，可以得到整個系統(tǒng)的平衡方程，然后建立力與節(jié)點坐標的線性關(guān)系。既可以將力作為初始參數(shù)求解相應(yīng)的平衡狀態(tài)，也可以將力密度作為描述力學(xué)平衡狀態(tài)的參數(shù)，不同的力密度對應(yīng)著不同的平衡狀態(tài)。與力密度法相比，歐拉坐標找形法建立了更簡單、更方便的拓撲關(guān)系，易于理解和計算。傳統(tǒng)的力密度法只能針對固定節(jié)點的邊界條件，而不能選擇性地控制節(jié)點的三維坐標。歐拉坐標找形法能夠在索網(wǎng)結(jié)構(gòu)平面投影固定的條件下，無需復(fù)雜迭代就可以完成找形計算，有利于索網(wǎng)結(jié)構(gòu)的外形控制和現(xiàn)場施工。

1 歐拉坐標找形法

1.1 計算假定

作出如下計算假定：

（1）忽略桿件抗彎剛度的影響，視作桿單元。

（2）單個桿件的自重均布荷載按照有限元基本理論等分到桿單元兩端節(jié)點，形成等效節(jié)點力。

（3）每個桿件作為獨立單元，截面面積保持不變。

1.2 基本原理

空間歐拉坐標系桿單元如圖1所示?；谏鲜黾僭O(shè)，結(jié)構(gòu)的合理受力狀態(tài)需要滿足力學(xué)平衡的2個條件：①每個桿單元的橫向、豎向長度與桿端力成比例，滿足桿單元平衡；②每個節(jié)點各個方向的桿端力的代數(shù)和為零，即滿足節(jié)點力平衡。

圖1 空間桿單元示意圖Fig.1 Schematic diagram of spatial rod element

在桿長和軸力已知的情況下，若要保持桿單元平衡，桿端力之間則必然符合一定的三角函數(shù)關(guān)系。此外，在迭代過程中若要滿足正確的收斂趨勢，則各力之間符號的相對關(guān)系是不變的。

該桿單元向量可表示為：

式中：lij為桿單元向量；i、j、k分別為3個方向的單位向量。

設(shè)定軸力符號與桿單元i端力的符號保持一致，以受拉為正。根據(jù)桿單元平衡可以得到力與坐標的關(guān)系式，如下所示：

式中：Nij為軸力；lij為桿長度。通過整理可以得到以下桿單元平衡迭代方程：

在桿長和軸力已知的情況下，可以直接建立單元桿端力與節(jié)點坐標的桿單元平衡迭代矩陣。令

可得

式中：qij為力密度，定義為單位長度的軸力，即Nij/lij=qij。

1.3 桿件總長度最小狀態(tài)方程

結(jié)構(gòu)桿件的總長度是衡量成本的一項重要指標。在靜力平衡狀態(tài)下，結(jié)構(gòu)內(nèi)力與桿件長度達到比例平衡，通過控制兩者的比例關(guān)系，理論上能達到結(jié)構(gòu)桿件總長度最小的目標?？梢园l(fā)現(xiàn)，矩陣K與坐標向量H的乘積表示相鄰節(jié)點之間在3個坐標軸上的投影距離。

設(shè)桿件的長度向量為l，則長度平方加權(quán)總和為

對稱矩陣K由三部分組成，即Kx、Ky、Kz，如下所示：

對于lTl取得最小值的問題，可以通過對x、y、z求導(dǎo)來解決，如下所示：

合并為

比較式（4）與式（13）可以發(fā)現(xiàn)，桿件最小長度的狀態(tài)方程與平衡方程相似，只需將式（4）中qij=Nij/lij=1，并且荷載向量P=0，兩式就完全等價，即：當(dāng)結(jié)構(gòu)荷載為零時，力密度取值為1，最小長度問題就可以使用歐拉坐標找形法的平衡方程來求解，此時得到的是桿件總平方長度最小的形狀。由Schek［6］提出的力密度法推導(dǎo)得到的結(jié)論與此相同，說明了歐拉坐標找形法的真實性和有效性。

對于更一般的情況，可以將長度平方加權(quán)累加。設(shè)置特殊的權(quán)重gij=Nconst/lij，其中Nconst為結(jié)構(gòu)桿件內(nèi)力，可以通過以下計算式得到總長度最小的狀態(tài)方程：

式中：G為權(quán)重向量，由各單元的權(quán)重組成。

1.4 迭代格式

①②2個桿單元總共3個節(jié)點，2號節(jié)點為共同節(jié)點。按照式（3），根據(jù)平衡條件與比例關(guān)系，建立格式相同的平衡方程。可以發(fā)現(xiàn)，2個桿單元可以依照組集規(guī)則，構(gòu)建總體平衡迭代矩陣。為方便組集規(guī)則演示，可以表示為與其中q1表示①單元的力密度，q2表示②單元的力密度，k11、k12、k21、k22、k23、k32、k33為剛度矩陣元素。

以共節(jié)點為例，組集總體平衡迭代矩陣如下所示：

類似地，組集所有的節(jié)點單元平衡迭代矩陣可以形成總體平衡迭代矩陣Ktotal。迭代過程中，平衡迭代矩陣是關(guān)于上一階段（第（n-1）階段）桿件軸力和長度的函數(shù)，可以表示為Ktotal(Nn-1，ln-1)，其中Nn-1為第(n-1)階段的軸力向量，ln-1為第(n-1)階段的長度向量。依照組集規(guī)則，可以推導(dǎo)出第n階段整個節(jié)點的總體平衡迭代方程式，如下所示：

式中：Ptotal為結(jié)構(gòu)的外荷載向量；Hn為第n階段節(jié)點坐標向量。軸力向量的元素與長度向量的元素分別為

式（17）和式（18）中：Fix，n-1、Fiy，n-1、Fiz，n-1分別為第(n-1)階段i端的桿端力；xi，n-1、yi，n-1、zi，n-1、xj，n-1、yj，n-1、zj，n-1分 別 為 第(n-1)階 段i、j點 的 坐 標；lij，n-1、Nij，n-1分別表示迭代第(n-1)階段所確定的單元桿長與軸力。

1.5 平衡方程的迭代求解條件

（1）初始參數(shù)的設(shè)定

歐拉坐標找形法初始參數(shù)設(shè)置可分為2種：①確定桿件恒定軸力向量N，在迭代過程中通過設(shè)置Nn=Nn-1ln/ln-1達到索網(wǎng)軸力均勻的目的；②根據(jù)經(jīng)驗調(diào)整桿件力密度，直接得到索網(wǎng)找形后的節(jié)點坐標。

（2）邊界條件的處理

因為總體平衡方程的未知數(shù)是坐標，而不是有限元中的位移，所以在考慮邊界條件時，對已知的節(jié)點坐標類似于強迫位移進行處理。采用主元充大數(shù)法，對于總體平衡迭代矩陣的第i個自由度點坐標Ui，原第i個平衡方程為

式中：m和n為剛度矩陣維度；Pi為荷載向量的元素。將矩陣對角元中充一大數(shù)D，并采取修改荷載項的方式進行處理，平衡方程變?yōu)?/p>

由于大數(shù)D數(shù)值很大，可得到Um≈U*i，從而將限制坐標體現(xiàn)在平衡迭代方程中。

（3）收斂條件的確定

第1種初始參數(shù)設(shè)置方法需要增加收斂條件，相較于上一迭代階段的3個方向坐標差平方和epx小于10-6，計算式如下所示：

滿足條件后，說明迭代過程中節(jié)點坐標不再發(fā)生變化。

1.6 歐拉坐標找形法的流程

歐拉坐標找形法的核心思想是在整體坐標系內(nèi)，按照桿單元的力學(xué)特性直接建立節(jié)點坐標和桿端力的平衡迭代矩陣，將“坐標-桿端力”的平衡迭代矩陣比擬成傳統(tǒng)“位移-桿端力”的有限元剛度矩陣。將坐標約束條件比擬成強迫位移條件，按照一定規(guī)則集合平衡迭代矩陣，從而建立總體平衡方程的迭代格式。依據(jù)前述原理利用Matlab軟件編制求解程序，該程序的流程如圖2所示。

圖2 歐拉坐標找形法流程Fig.2 Flow chart of Euler coordinates form-finding method

2 應(yīng)用

歐拉坐標找形法是從桿件無彎矩的假定出發(fā)，將初始結(jié)構(gòu)中的節(jié)點坐標與桿單元軸力作為初始參數(shù)，并建立迭代格式，從而求解出無彎矩狀態(tài)下的最終結(jié)構(gòu)形態(tài)。從力學(xué)角度來看，最終得到的結(jié)構(gòu)形態(tài)是最為合理的一種受力形狀。本節(jié)中主要對坐標迭代法進行驗證，對邊界條件、網(wǎng)格類型、初始參數(shù)等參數(shù)進行研究，以說明歐拉坐標找形法的適用性。

2.1 無荷載結(jié)構(gòu)的最小長度狀態(tài)找形

由第1.3節(jié)的推導(dǎo)可知，當(dāng)外荷載為零、網(wǎng)殼的桿件力相同時，桿件總長度最小。在數(shù)學(xué)領(lǐng)域，普拉托提出過極小曲面的概念，是由肥皂泡引發(fā)的猜想，即利用一根鐵絲，就會形成一個處于平衡狀態(tài)的彩色薄膜。如果忽略混合液體自身的質(zhì)量，也不考慮風(fēng)力等外部干擾因素，薄膜的勢能在表面張力的作用下就會達到最小值，而肥皂膜所呈現(xiàn)的曲面形狀必然具有最小的面積。目前具有解析式的常見極小曲面有懸鏈線曲面、螺旋線曲面和Scherk曲面，利用本方法可以便捷地形成上述3種曲面。

圖3 懸鏈線曲面Fig.3 Catenary curved surface

圖4 螺旋線曲面Fig.4 Spiral surface

圖5 Scherk曲面Fig.5 Scherk’s surface

2.2 雙曲拋物面找形

為了進一步驗證本方法的正確性和適用性，選取建筑設(shè)計中最為常見的雙曲拋物面作為實例。雙曲拋物面網(wǎng)格參數(shù)與迭代條件按照文獻［22］設(shè)置，正方形平面邊長為10，立面高度為5。初始網(wǎng)格如圖6a所示，其平面投影如圖6b直線規(guī)則部分所示。圖2中流程1的迭代結(jié)果如圖6c所示，所有桿件軸力為1，網(wǎng)格分布均勻，也驗證了此類型網(wǎng)格適用于流程1。從圖6b平面投影變化可以看出，網(wǎng)格向內(nèi)收縮，對角線網(wǎng)格位置不變，其他網(wǎng)格向?qū)蔷€偏移靠近，越遠離邊界的網(wǎng)格偏移的幅度越大。

圖6 雙曲拋物面Fig.6 Hyperbolic paraboloids surface

2.3 馬鞍形曲面找形

馬鞍形曲面平面網(wǎng)格參數(shù)與雙曲拋物面相同，正方形平面邊長為10，初始網(wǎng)格如圖7a所示。邊界條件為兩邊半徑5的圓弧，另外兩邊固定。圖2中流程1的迭代結(jié)果如圖7b所示，所有桿件軸力為1，馬鞍形曲面特征明顯。由如圖7c所示的平面投影變化可以看出，網(wǎng)格呈現(xiàn)對稱分布，網(wǎng)格均勻。對角線網(wǎng)格為適應(yīng)馬鞍形曲面的特征，其位置發(fā)生變化，網(wǎng)格向兩固定邊收縮。

圖7 馬鞍形曲面Fig.7 Saddle-shaped surface

2.4 邊界條件的影響

從雙曲拋物面和馬鞍形曲面的實例可以看出，相同的網(wǎng)格劃分會因為不同的邊界條件而產(chǎn)生不同的優(yōu)化結(jié)果。為了進一步探究邊界條件對于算法找形結(jié)果的影響，選取與馬鞍形曲面相同的平面網(wǎng)格，將圓弧邊界保留，釋放另外兩邊的約束。圖2中流程1的迭代結(jié)果如圖8a所示，所有桿件力為1。由如圖8b所示的平面投影可以看出，釋放邊界約束后，網(wǎng)格向中間猛烈收縮，導(dǎo)致中間網(wǎng)格密度較大，同時受約束的邊界網(wǎng)格則分布較為均勻。從極小曲面實例也可以看出，邊界約束會影響周圍網(wǎng)格的長度和分布均勻程度。

圖8 四點固定的馬鞍形曲面Fig.8 Four-point fixed saddle-shaped surface

歐拉坐標找形法可以在邊界條件中指定三維坐標的數(shù)值，而力密度法只能限制固定點，因此歐拉坐標找形法具有很大的優(yōu)勢。選擇了如圖6所示的雙曲拋物面，邊界和網(wǎng)格劃分條件保持不變，限制平面坐標，從而實現(xiàn)投影網(wǎng)格在優(yōu)化過程中不變的目標。圖2中流程1的迭代結(jié)果如圖9所示，三維形狀變化不大，并且平面投影網(wǎng)格更加均勻，這說明了本方法控制坐標的優(yōu)越性。

2.5 網(wǎng)格劃分的影響

選取螺旋線曲面，采用三角形網(wǎng)格劃分，其余參數(shù)和圖4a的參數(shù)相同，結(jié)果如圖10a所示。圖2中流程1的迭代結(jié)果如圖10b所示，桿件力均相同。圖10b與圖4b相同，因此網(wǎng)格劃分對于找形結(jié)果無明顯影響。

圖10 三角形網(wǎng)格的螺旋線曲面Fig.10 Spiral surface with triangle mesh

3 結(jié)論

（1）理論和實例都驗證了在無荷載條件下，最終形狀能夠滿足總長度最小或者總長度平方和最小。

（2）本方法相較于傳統(tǒng)力密度法，可以方便地改變邊界，控制節(jié)點三維坐標，并且迭代格式簡單，適應(yīng)性較強，收斂穩(wěn)定。

（3）邊界條件對最終的找形結(jié)果影響較大。相同的初始條件，改變邊界條件會產(chǎn)生完全不同的找形結(jié)果。同時，內(nèi)力一致或力密度一致時，邊界處節(jié)點坐標相較于初始狀態(tài)變化較小，遠離邊界處節(jié)點坐標變化較大，而且桿件長度更小。網(wǎng)格劃分的方式不同（如四邊形網(wǎng)格和三角形網(wǎng)格）基本不會影響最終的找形結(jié)果。

作者貢獻聲明：

劉超：提出算法，指導(dǎo)算法推導(dǎo)。

詹海鵬：算法的編程和算例的計算，論文撰寫。