趙廷紅，陳 健

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

0 引 言

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

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

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

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

1 剛度中心的計算原理

1.1 基本假設(shè)

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

1.2 理論公式

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

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

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

也即：

其中，［δ］為結(jié)構(gòu)的柔度矩陣，δ11至δ33為對應單位荷載作用下相應的位移，即柔度系數(shù)。

由前述有：

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

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

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

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

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

1.3 剛性樓板的處理

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

在幾何建模階段將樓板單獨剖分出來進行建模，計算時將其剛度行為設(shè)置為剛性，則內(nèi)部程序?qū)跇前宓馁|(zhì)心處創(chuàng)建一個包含其整體質(zhì)量、質(zhì)心坐標值以及轉(zhuǎn)動慣量等動力特性的質(zhì)量單元，并作為導向節(jié)點來表達樓板的剛體運動，故而不必再劃分應力應變的網(wǎng)格。但為了讓樓板在后處理中也能顯示出來，通過設(shè)置讓程序劃分出表示樓板幾何形狀的網(wǎng)格。剛體與柔性體（樓板與結(jié)構(gòu)主體）之間必須設(shè)置接觸來進行連接，二者不能合并在一個整體（part）內(nèi)，而剩余的柔性體（結(jié)構(gòu)主體與地基）將合并為一個整體，使所有柔性體連續(xù)、共用節(jié)點。

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

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

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

2 工程實例

2.1 工程概況

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

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

2.2 坐標系設(shè)定

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

圖3 坐標系設(shè)定Fig.3 Coordinate setting

2.3 計算模型

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

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

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

表1 網(wǎng)格參數(shù)Tab.1 Mesh parameter

表2 各組荷載作用下方案3的結(jié)構(gòu)扭轉(zhuǎn)位移Tab.2 The structural torsional displacement of scheme 3 under each group of loads

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

2.4 剛度中心的求解及驗證

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

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

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

圖6 結(jié)構(gòu)變形Fig.6 Structural deformation

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

表3 各方案的結(jié)構(gòu)扭轉(zhuǎn)位移Tab.3 Structural torsional displacement of each scheme

按照剛度中心的定義，若在求解得出的剛度中心上施加水平力，結(jié)構(gòu)應再無扭轉(zhuǎn)位移產(chǎn)生，但由于數(shù)值計算存在誤差，仍需通過結(jié)構(gòu)的響應來驗證其精確性，故在所求解得出的剛度中心上施加平行于x與y軸的作用力來進行驗證，記為Fxr與Fyr。在ANSYS中施加荷載需要具體的幾何元素或網(wǎng)格節(jié)點，而在求解得出的剛度中心處沒有相應的加載點?？紤]到樓板為剛體，根據(jù)力的平移定理，在剛度中心（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

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

圖8 精確性驗證中的結(jié)構(gòu)變形Fig.8 Structural deformation in accuracy verification

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

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

圖9 網(wǎng)格加密對扭轉(zhuǎn)位移的影響Fig.9 Influence of mesh refinement on torsional displacement

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

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

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

2.5 關(guān)于剛度中心數(shù)值精度的討論

在本文的有限元模型中，當作用力大小、方向均不變時，隨著網(wǎng)格加密，扭轉(zhuǎn)位移θx、θy與θz的變化是不一致的，結(jié)合理論知識，這種現(xiàn)象說明結(jié)構(gòu)本身的剛度分布發(fā)生了變化。從x與y方向上的平動以及繞z軸的扭轉(zhuǎn)這3 個自由度建立如式（11）所示的受力平衡方程［17］，由此可知結(jié)構(gòu)在x、y方向平動位移值的大小以及繞z軸扭轉(zhuǎn)位移值的大小最終取決于剛度矩陣中x、y方向的抗側(cè)剛度的分布。根據(jù)有限元方法的特點，產(chǎn)生該現(xiàn)象是因為x與y方向上的網(wǎng)格加密不均勻?qū)е掠嬎隳Ｐ驮诰W(wǎng)格加密的過程中x與y方向上的抗側(cè)剛度變化程度不一樣的緣故。

其中：

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

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

表5 各方案的結(jié)構(gòu)位移Tab.5 Structural displacement of each scheme

圖10 網(wǎng)格加密對結(jié)構(gòu)位移的影響Fig.10 Influence of mesh refinement on structural displacement

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

3 結(jié) 論

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

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

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