基于ANSYS 的水電站廠房剛度中心數值分析

趙廷紅，陳 健

（蘭州理工大學能源與動力工程學院，甘肅 蘭州 730050）

0 引 言

在結構動力分析的過程中，剛度中心是一項極為重要的參數，大量研究表明［1，2］它與質量中心的相對位置決定了結構在地震等動力作用下發生平扭耦聯等不良現象的劇烈程度，隨著二者偏心距離的增大，結構的動位移與動應力也隨之增大，控制偏心距可有效避免結構在地震作用下因扭轉響應過大而發生扭轉破壞。對此，我國相關規范［3］也對結構的剛度分布進行了一定的約束，從而保證結構在設計層面的合理性，同時防止因建筑物的剛度中心與質量中心之間的距離過大而出現的動力響應異常增大、乃至危及結構體系安全的問題。

對于單一單層建筑結構，通常直接采用靜力學法來求解其剛度中心，即先分析整體結構的剛度分布情況，而后直接應用靜力學公式求解出該結構抗側反力的合力點［4］，即得出其剛度中心。該方法也可用于求解大型多層復雜結構剛度中心的計算，但其求解過程將非常復雜，而且也很難得到理想的結果。為了克服這種現象，研究者［5，6］提出了一些求解多層復雜結構剛度中心的方法，這些方法雖然各有特色，但它們的主要思路都是：當整體結構只在某層剛度中心作用水平力時，該層結構只有平動而不發生扭轉，而其他層可以有扭轉和平動，而后通過觀測結構各部分在荷載作用下的響應去求解剛度中心。按照該思路，該層的剛度中心既反映了本層構件的剛度分布又反映了整體結構中其他各層構件的剛度分布情況。

目前，水電站廠房抗震研究已有大量成果并逐漸聚焦于流固耦合模擬以及非線性分析等方面，但從整體來看，其地震響應研究是滯后的［7］，關于水電站廠房平扭耦聯效應的分析以及相應的剛度中心計算求解尚未獲得關注。隨著計算機技術的飛速發展，結構仿真模擬軟件也應運而生，結合有限元技術，研究人員可以非常方便地通過數值計算近似地得到結構的剛度中心，但現階段具備剛度中心自動求解功能的工程軟件局限性太大，不適用于除土木工程外其他行業的應用，無法很好利用有限元方法的優勢。水電站廠房與土木工程中的框架結構類似，但因其內部設有廊道、流道、孔洞等功能區域，而這些因素使得它的幾何外形較土木工程中的框架結構體系要復雜得多，若將其內部的各部分構件按梁、板、柱等形狀規則的構件進行簡化再用土木行業的數值模擬軟件去求解剛度中心，必然使得數值計算的準確性大打折扣，基于此，有必要將剛度中心的求解方法進行擴展。

ANSYS 是目前應用最為廣泛的一款有限元分析軟件，可以發揮出有限元方法最大的優勢，任意形狀的結構均可在精確三維建模后，通過增大結構的離散程度獲取可接受的數值結果。本文擬應用ANSYS 軟件結合現有求解土木工程結構剛度中心的方法來求解水電站廠房的剛度中心；此外，因考慮到有限元網格的劃分是影響數值計算結果精度的主要因素之一，本文擬將研究網格的劃分對剛度中心計算精度的影響，以期能為水工結構求解剛度中心時提供參考。

1 剛度中心的計算原理

1.1 基本假設

在求解建筑結構的剛度中心時，必須遵循兩條基本假設［8］：①假定樓板在平面內為無限剛性，即受到平面內作用力時僅有剛體位移產生而無彈性變形；②剛性樓板上存在這樣一個點，當平面內作用力的合力通過此點時，樓板只有平動位移而不發生扭轉位移，否則反之。

1.2 理論公式

根據相關文獻［6］，基于結構整體各層相互影響下剛度中心求解方法的理論推導過程如下：如圖1，設某層剛度中心的坐標值為R（xr，yr），在質心或任意非剛度中心點O（x，y）處施加平行于x軸方向的單位力（無量綱數1），記為fx，則剛度中心R（xr，yr）處將隨之產生抗側合力，記為p，則fx與p二者同時形成一組力偶mx，樓板會由此產生扭轉位移δx，若在P（x，y）處作用反向力偶mz，則樓板由此所產生扭轉位移與δx方向相反，二者可相互抵消，當mz=mx時，該層樓板的扭轉位移歸零。

圖1 計算原理Fig.1 Theory of calculation

基于上述思路，以O（x，y）為原點，列出結構由x、y及θ三個方向為分量且以柔度矩陣（即剛度矩陣的逆矩陣）表示的受力平衡方程，如下式（扭轉位移與力矩以逆時針為正）：

也即：

其中，［δ］為結構的柔度矩陣，δ11至δ33為對應單位荷載作用下相應的位移，即柔度系數。

由前述有：

求解θ=0，可得下式（力偶與扭轉位移以逆時針為正）：

此時注意到：δ31=δx，即上述推導中施加單位力后樓板的扭轉位移；而δ33為單位扭矩作用下樓板的扭轉位移，記為δz。所以有：

上式中僅有剛度中心的坐標值yr未知，故可求得。同理，可求得剛度中心的另一個坐標值xr。經整理，考慮結構整體各層相互影響時剛度中心的求解公式如下：

由于在ANSYS中無法添加單位荷載，所以柔度系數并不能夠直接從柔度矩陣中獲取，而且對于節點數量龐大的有限元模型而言，直接分析其剛度矩陣的組成也是不現實的。根據柔度系數的定義，δx、δy與δz的值可由點O處作用實際的荷載Fx、Fy及Mz后產生實際扭轉位移θx、θy及θz求解其比值獲得，所以可將計算公式可變換為：

式中：xir、yir為第i層樓板剛度中心的坐標值；xi、yi為第i層樓板荷載作用點的坐標值；θix為第i層樓板荷載作用點沿x軸正方向施加作用力時樓板的扭轉位移；θiy為在第i層樓板荷載作用點沿y軸正方向施加作用力時樓板的扭轉位移；θiz為第i層樓板荷載作用點沿逆時針方向施加力矩時樓板的扭轉位移；Miz為在第i層樓板荷載作用點沿逆時針方向作用的力矩；Fix為第i層樓板荷載作用點沿x正方向作用的水平力；Fiy為第i層樓板荷載作用點沿y正方向作用的水平力。

1.3 剛性樓板的處理

在數值計算中，不同構件之間剛度差異太大將使得計算結果誤差異常增大，甚至引起求解無法收斂的情況［9］。利用公式（8）、（9）求解結構剛度中心需要將樓板視為剛體，而建筑結構的其余部分皆為柔性體，為克服上述困難，在計算時采用ANSYS中的剛柔耦合仿真分析方法來進行處理，即通過接觸算法將剛體與柔性體接觸面上的節點進行耦合，從而將二者連接在一起，該方法廣泛應用于機械制造、軌道交通等研究領域，是多體系統動力學與有限元分析方法相結合的表現［10］。結構剛度中心其實就是約束剛性樓板位移的柔性體所產生抗側反力的合力點，所以在剛柔耦合分析過程中并不需要進行瞬態動力學分析，這類問題只進行靜力學分析便能滿足要求［11］。

在幾何建模階段將樓板單獨剖分出來進行建模，計算時將其剛度行為設置為剛性，則內部程序將會在樓板的質心處創建一個包含其整體質量、質心坐標值以及轉動慣量等動力特性的質量單元，并作為導向節點來表達樓板的剛體運動，故而不必再劃分應力應變的網格。但為了讓樓板在后處理中也能顯示出來，通過設置讓程序劃分出表示樓板幾何形狀的網格。剛體與柔性體（樓板與結構主體）之間必須設置接觸來進行連接，二者不能合并在一個整體（part）內，而剩余的柔性體（結構主體與地基）將合并為一個整體，使所有柔性體連續、共用節點。

由于本研究并不需要進行碰撞、摩擦等高度非線性的接觸分析，樓板與結構主體在分析過程中始終保持緊密連接，所以將剛體與柔性體的接觸類型設置為綁定連接。剛體與柔性體的接觸算法采用多點約束法，該方法是ANSYS眾多連接方法中適用性較為廣泛、計算精度較高的一種方法，多點約束算法不僅適用于剛柔耦合分析，也可以將不同自由度的單元通過自動創建約束方程聯系起來，所以用多點約束法對水電站廠房這類復雜建筑結構進行模擬前處理，將會大大減小數值模擬前處理的工作量［12］。多點約束法的本質是通過一個或多個節點的自由度為標準值，令其余節點的自由度與之通過數學方程建立耦合關系［12，13］，如式（10）所示，使用該方法將樓板處理為剛體后如圖2所示。

圖2 基于多點約束法的剛性樓板Fig.2 Rigid floor based on multi-point constraint method

式中：j為主節點的節點編號；i為從節點的節點編號；U為自由度，如位移、溫度等，在本文中為位移值；Ci為加權系數；C0為常數項。

2 工程實例

2.1 工程概況

某水電站是河段上游開發規劃中的第六個梯級電站，該水電站以發電為主，同時兼顧防洪、航運等功能。工程規模等級為Ⅱ等大（2）級，按百年一遇的洪水進行設計，千年一遇的洪水進行校核，對應的洪峰流量為30 000和38 100 m3/s，壩址處的多年平均流量732 m3/s，正常蓄水位對應的水庫庫容為1.72 億m3。廠內溢流式水電站廠房與混凝土重力壩在結構形式、應力應變行為等方面十分相近［14］，后文以壩體、壩段進行表述。

從右岸至左岸分別為右副壩、航運壩段、泄洪閘壩段、電站廠房壩段、安裝間壩段以及左副壩。其中，電站廠房壩段每兩個發電機組由兩側的橫縫分割為一個機組壩段，每個機組壩段的壩體下部是發電機組的過流通道，壩體上部是表孔溢洪通道，發電機組與表孔泄洪的流道交疊布置，水電站廠房在正常工作階段有蓄水功能，當洪水來臨時啟用表孔進行泄洪，過流形式為廠內溢流式，結構體系的質量與剛度分布錯綜復雜。

2.2 坐標系設定

ANSYS 中的坐標系遵循右手定則，全局坐標系原點位于壩體上游迎水面與地基以及左岸相交處，左岸至右岸的垂直方向為x軸正方向、上游至下游的垂直方向為z軸正方向。在全局坐標系中，整體模型x方向寬35 m、y方向高225 m、z方向長369 m；壩體x方向寬35 m、y方向高75 m、z方向總長71.3 m。為便于求解剛度中心，特建立如圖3所示的局部坐標系，其原點位于全局坐標系的（35，70，-4.5）處，右岸至左岸的垂直方向平行于x軸正方向、上游至下游的垂直方向平行于y軸正方向。本文以下所有計算及分析過程均基于該局部坐標系。

圖3 坐標系設定Fig.3 Coordinate setting

2.3 計算模型

在以往與本文研究對象相似的動靜力研究中［15，16］，取典型壩段進行分析是普遍做法，本文基于這種常規的處理方法去求解結構的剛度中心。以一段壩段為研究對象，將壩體頂部的樓板按第一節的方法設置為剛性，采用無質量地基法來考慮結構和壩基的相互作用，壩基沿上、下游及深度方向延伸1.5 至2 倍壩高，壩體結構密度為2.5×103kg/m3、彈性模量為2.8×1010Pa、泊松比為0.167；壩基的彈性模量為3×109Pa、泊松比為0.3。壩基底部為全約束，壩基四周為法相約束，由于壩段兩側由橫縫分割，本文視壩段間無相互作用，各個壩段獨立承受荷載作用［15］，壩體兩側視為自由面。需要注意壩基的屬性尺寸、邊界條件的設定等各方面因素均對結構剛度中心的求解結果產生影響。

有限元網格的類型對計算結果的精度有所影響，為得到質量較好的有限元網格，將幾何模型切分為各個可以掃掠的體，劃分為八節點六面體以及六節點五面體的網格單元。遠離壩體一定范圍外的壩基應力應變程度很小，將此部分網格尺寸增大以減輕計算成本。

在實際工程中，對于同一結構體系，數值計算往往需要劃分不同的網格方案以獲得具有相當精度的計算結果，而在網格整體加密的過程中，結構的應力應變將逐步收斂于真實值，即結構的剛度（柔度）發生了變化，所以剛度中心的坐標值也將隨之改變，為分析剛度中心的求解在網格加密過程中的變化，對有限元模型進行整體加密，改變網格尺寸大小以及掃掠路徑的單次長短，使壩體與壩基同步加密。本文共設置五組不同的網格方案進行對比分析，各方案的網格參數如表1所示，因篇幅所限，在此只列出方案3的網格劃分圖（圖4）。

表1 網格參數Tab.1 Mesh parameter

表2 各組荷載作用下方案3的結構扭轉位移Tab.2 The structural torsional displacement of scheme 3 under each group of loads

圖4 有限元模型Fig.4 Finite element model

2.4 剛度中心的求解及驗證

因篇幅所限，在此僅以方案3 為例敘述剛度中心的求解以及驗證過程。將荷載作用點設為（35，0，0），作用平行于x軸的作用力Fx=1 N、平行于y軸的作用力Fy=1 N、繞z 軸的力矩Mz=1 N·m，邊界條件的設定與荷載的施加如圖5所示。

圖5 邊界條件及荷載施加Fig.5 Boundary conditions and load application

結構變形如圖6，在Fx作用下結構發生明顯的不均勻位移變形，結構整體向x軸正方向彎曲，受荷載直接作用的上游側總體位移量比下游側大，約為1.5 至1.8 倍，扭轉位移強烈；在Fy作用下結構的變形位移規律與Fx作用時基本一致，結構整體向y軸正方向彎曲，受荷載直接作用的左岸側的總體位移量比右岸側大，約為1.3至1.5倍，扭轉位移明顯；在Mz作用下結構變形以繞z軸的扭轉變形為主，從結構中心部位至最外圍的位移變形逐漸增大，壩體沿河流方向三面墻的變形位移程度并無太大差異。結構在各荷載作用下的θx、θy及θz為1.995 6×10-12rad、8.209 8×10-13rad、4.6971×10-14rad。

圖6 結構變形Fig.6 Structural deformation

根據計算得出的扭轉位移數值，結合剛度中心求解公式（8）、（9），可得出方案3 剛度中心的x和y坐標值分別為17.522 503 4與42.452 870 24。將上述過程中荷載的量級放大，令Fx=1 N、Fy=1 N 及Mz=1 N·m 為a 組，Fx=100 N、Fy=100 N 及Mz=100 N·m 為b 組，Fx=100 00 N、Fy=100 00 N 及Mz=100 00 N·m為c 組，各組荷載作用下方案3 的扭轉位移以及對應的剛度中心如表3所示，可以看出，荷載大小與結構的響應呈嚴格線性關系，每組剛度中心的坐標值恒定不變，結合剛度中心計算公式的推導過程可知，當支撐樓板的結構處于線彈性階段時，柔度系數的值不會隨著荷載大小發生而改變，所以yir與xir為確定值，這種情況可以解釋為結構剛度中心的位置保持不變，其位置只與結構本身的剛度（柔度）分布有關。

表3 各方案的結構扭轉位移Tab.3 Structural torsional displacement of each scheme

按照剛度中心的定義，若在求解得出的剛度中心上施加水平力，結構應再無扭轉位移產生，但由于數值計算存在誤差，仍需通過結構的響應來驗證其精確性，故在所求解得出的剛度中心上施加平行于x與y軸的作用力來進行驗證，記為Fxr與Fyr。在ANSYS中施加荷載需要具體的幾何元素或網格節點，而在求解得出的剛度中心處沒有相應的加載點。考慮到樓板為剛體，根據力的平移定理，在剛度中心（17.522 503 4，42.452 870 24，0）處施加的作用力完全可以沿其作用方向等效地平移至剛性樓板的外部面上，故在（0，42.452 870 24，0）處施加平行于x軸的作用力Fxr=1 N、在（17.522 503 4，0，0）處施加平行于y軸的作用力Fyr=1 N，如圖7。

圖7 精確性驗證中的荷載施加Fig.7 Load application in accuracy verification

施加作用后對應的結構變形如圖8 所示。由圖8 可以看出，在Fxr作用下結構位移變形均勻，結構整體向x軸正方向彎曲，扭轉位移難以直接觀測；在Fyr作用下結構的位移變形與Fxr作用時的規律基本一致，結構整體向y軸正方向彎曲，結構無明顯扭轉位移產生。樓板在Fxr、Fyr作用下對應的扭轉位移θxr、θyr分別為-2.927 7×10-17rad、8.870 2×10-19rad。分別對比Fx與Fxr、Fx與Fyr作用下結構的位移變形可知，當荷載作用在求解得出的剛度中心時，樓板的扭轉位移明顯減弱。

圖8 精確性驗證中的結構變形Fig.8 Structural deformation in accuracy verification

以在方案3 結構的剛度中心上施加荷載前后的響應為例，若求解的剛度中心準確無誤，則在剛度中心上施加作用力時應無扭轉產生，但實際上仍然有-2.927 7×10-17rad 的扭轉位移殘余，但θxr十分微小，Δxr=|-2.927 7×10-17|/1.995 6×10-12×100%≈0.001 467 078%，同理，Δyr約為0.000 108 044%，誤差在0.005%以下，認為求解得出的剛度中心滿足要求。

運用有限元方法進行數值計算的特點是數值結果會隨著結構離散化程度的增大而逐漸收斂于精確解，所以數值模擬需要進行網格無關性的驗證來說明結果的準確，而剛度中心的坐標值取決于各柔度系數的比值，對扭轉位移的變化比較敏感，若各方案的扭轉位移在網格加密過程中變化得不一樣，則各方案下計算得出的剛度中心也不盡相同。不同網格方案下計算得出的樓板扭轉位移θx、θy及θz如表3 所示，其變化趨勢如圖9所示。

圖9 網格加密對扭轉位移的影響Fig.9 Influence of mesh refinement on torsional displacement

由表3 及圖9 可以看出：所有方案中結構扭轉位移在加密前后的誤差不大于1%，可認為本次試驗所有方案均已達到與網格密度無關的狀態，當網格數量大約高于25 萬時，結構位移變形的變化將更加平緩；各方案結構扭轉位移的變化程度并不一樣，3 組數據中，θx每次加密后的變化程度最大，θy次之，θz最小，所以后續以結構扭轉位移計算得出各方案的剛度中心也存在著偏差。

各方案剛度中心的計算值與精確性驗證如表4所示。由表可看出各方案的剛度中心坐標值在某個范圍內變化，通過對比表3 與表4 的數據，這是因為網格劃分發生變化時，結構的剛度（柔度）發生了不均勻的變化，使得每個方案求解得出的剛度中心總是不一樣的，這也導致了θxr與θyr的值并不能嚴格呈現下降趨勢；總體而言，在各方案剛度中心上施加荷載所引起的扭轉位移甚微，Δxr與Δyr皆不大于0.005%，對實際工程而言滿足精度要求。

表4 各方案的剛度中心坐標及精確性驗證Tab.4 Stiffness center coordinates of each scheme and accuracy verification

2.5 關于剛度中心數值精度的討論

在本文的有限元模型中，當作用力大小、方向均不變時，隨著網格加密，扭轉位移θx、θy與θz的變化是不一致的，結合理論知識，這種現象說明結構本身的剛度分布發生了變化。從x與y方向上的平動以及繞z軸的扭轉這3 個自由度建立如式（11）所示的受力平衡方程［17］，由此可知結構在x、y方向平動位移值的大小以及繞z軸扭轉位移值的大小最終取決于剛度矩陣中x、y方向的抗側剛度的分布。根據有限元方法的特點，產生該現象是因為x與y方向上的網格加密不均勻導致計算模型在網格加密的過程中x與y方向上的抗側剛度變化程度不一樣的緣故。

其中：

式中：［K］為剛度矩陣；xi、yi為第i個抗側構件以荷載作用點為原點的坐標值；∑k（ex，i）、∑k（ey，i）為結構x、y方向的抗側剛度之和；∑k（ex，i）yi2+∑k（ey，i）xi2為結構的抗扭剛度。

以不同方案的結構在Fx及Fy作用下的響應進行說明。令該點在Fx作用下x方向的位移值為α、y方向的位移值為β，該點在Fy作用下x方向的位移值為ζ、y方向的位移值為ξ，各位移值的變化曲線按等比例繪制于圖10，各方案中各位移值的具體數據如表5所示。

表5 各方案的結構位移Tab.5 Structural displacement of each scheme

圖10 網格加密對結構位移的影響Fig.10 Influence of mesh refinement on structural displacement

由表5及圖10可以看出：在4組結構響應曲線中，Fx作用時結構在x方向上的位移最大，Fy作用時結構在y方向上的位移次之，其余2 組的位移大小十分接近，從支撐樓板的三面墻體分析，其在y方向上的抗側剛度遠比x方向大，所以導致了α 曲線上的值比ξ 要小，計算結果符合物理規律；從4 組曲線的變化率可以看出，α 曲線遠比其他3 組曲線的變化程度要大，說明在網格加密的過程中x與y方向上的離散程度變化不均勻，使得對應方向上的抗側剛度變化不一致，最終導致剛度中心計算誤差的出現。

3 結 論

（1）使用ANSYS求解剛度中心是可行且高效的，為研究人員分析計算結構的剛度中心、掌握結構動力特性提供了新的途徑。

（2）由表3 及表5、圖9 及圖10 可知，隨著網格的不斷加密，結構的響應越來越趨近于真實值，但網格的加密導致計算模型在不同方向上的抗側剛度發生不均勻變化，從而使得剛度中心的坐標值發生偏差，為減小此過程中引起的誤差，網格的加密應盡量保持均勻；剛度中心驗證與求解過程中的計算模型需要保證是一致的，網格劃分參數不允許發生改變。

（3）對于（2）而言，在實際工程應用中，要做到整體網格完全均勻加密是不現實的。通過本試驗可看出，計算模型達到網格無關性后，繼續加密網格所引起的誤差是可以忽略的，所以，研究者可根據實際需求，綜合考慮重點關注的物理量以及剛度中心計算值的數值誤差求得剛度中心。