趙廷紅，陳 健

（蘭州理工大學(xué)能源與動(dòng)力工程學(xué)院，甘肅 蘭州 730050）

0 引 言

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

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

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

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

1 剛度中心的計(jì)算原理

1.1 基本假設(shè)

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

1.2 理論公式

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

圖1 計(jì)算原理Fig.1 Theory of calculation

基于上述思路，以O(shè)（x，y）為原點(diǎn)，列出結(jié)構(gòu)由x、y及θ三個(gè)方向?yàn)榉至壳乙匀岫染仃嚕磩偠染仃嚨哪婢仃嚕┍硎镜氖芰ζ胶夥匠?，如下式（扭轉(zhuǎn)位移與力矩以逆時(shí)針為正）：

也即：

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

由前述有：

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

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

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

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

式中：xir、yir為第i層樓板剛度中心的坐標(biāo)值；xi、yi為第i層樓板荷載作用點(diǎn)的坐標(biāo)值；θix為第i層樓板荷載作用點(diǎn)沿x軸正方向施加作用力時(shí)樓板的扭轉(zhuǎn)位移；θiy為在第i層樓板荷載作用點(diǎn)沿y軸正方向施加作用力時(shí)樓板的扭轉(zhuǎn)位移；θiz為第i層樓板荷載作用點(diǎn)沿逆時(shí)針?lè)较蚴┘恿貢r(shí)樓板的扭轉(zhuǎn)位移；Miz為在第i層樓板荷載作用點(diǎn)沿逆時(shí)針?lè)较蜃饔玫牧?；Fix為第i層樓板荷載作用點(diǎn)沿x正方向作用的水平力；Fiy為第i層樓板荷載作用點(diǎn)沿y正方向作用的水平力。

1.3 剛性樓板的處理

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

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

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

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

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

2 工程實(shí)例

2.1 工程概況

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

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

2.2 坐標(biāo)系設(shè)定

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

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

2.3 計(jì)算模型

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

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

在實(shí)際工程中，對(duì)于同一結(jié)構(gòu)體系，數(shù)值計(jì)算往往需要?jiǎng)澐植煌木W(wǎng)格方案以獲得具有相當(dāng)精度的計(jì)算結(jié)果，而在網(wǎng)格整體加密的過(guò)程中，結(jié)構(gòu)的應(yīng)力應(yīng)變將逐步收斂于真實(shí)值，即結(jié)構(gòu)的剛度（柔度）發(fā)生了變化，所以剛度中心的坐標(biāo)值也將隨之改變，為分析剛度中心的求解在網(wǎng)格加密過(guò)程中的變化，對(duì)有限元模型進(jìn)行整體加密，改變網(wǎng)格尺寸大小以及掃掠路徑的單次長(zhǎng)短，使壩體與壩基同步加密。本文共設(shè)置五組不同的網(wǎng)格方案進(jìn)行對(duì)比分析，各方案的網(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 剛度中心的求解及驗(yàn)證

因篇幅所限，在此僅以方案3 為例敘述剛度中心的求解以及驗(yàn)證過(guò)程。將荷載作用點(diǎn)設(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)位移強(qiáng)烈；在Fy作用下結(jié)構(gòu)的變形位移規(guī)律與Fx作用時(shí)基本一致，結(jié)構(gòu)整體向y軸正方向彎曲，受荷載直接作用的左岸側(cè)的總體位移量比右岸側(cè)大，約為1.3至1.5倍，扭轉(zhuǎn)位移明顯；在Mz作用下結(jié)構(gòu)變形以繞z軸的扭轉(zhuǎn)變形為主，從結(jié)構(gòu)中心部位至最外圍的位移變形逐漸增大，壩體沿河流方向三面墻的變形位移程度并無(wú)太大差異。結(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ù)計(jì)算得出的扭轉(zhuǎn)位移數(shù)值，結(jié)合剛度中心求解公式（8）、（9），可得出方案3 剛度中心的x和y坐標(biāo)值分別為17.522 503 4與42.452 870 24。將上述過(guò)程中荷載的量級(jí)放大，令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)位移以及對(duì)應(yīng)的剛度中心如表3所示，可以看出，荷載大小與結(jié)構(gòu)的響應(yīng)呈嚴(yán)格線性關(guān)系，每組剛度中心的坐標(biāo)值恒定不變，結(jié)合剛度中心計(jì)算公式的推導(dǎo)過(guò)程可知，當(dāng)支撐樓板的結(jié)構(gòu)處于線彈性階段時(shí)，柔度系數(shù)的值不會(huì)隨著荷載大小發(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)應(yīng)再無(wú)扭轉(zhuǎn)位移產(chǎn)生，但由于數(shù)值計(jì)算存在誤差，仍需通過(guò)結(jié)構(gòu)的響應(yīng)來(lái)驗(yàn)證其精確性，故在所求解得出的剛度中心上施加平行于x與y軸的作用力來(lái)進(jìn)行驗(yàn)證，記為Fxr與Fyr。在ANSYS中施加荷載需要具體的幾何元素或網(wǎng)格節(jié)點(diǎn)，而在求解得出的剛度中心處沒(méi)有相應(yīng)的加載點(diǎn)。考慮到樓板為剛體，根據(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 精確性驗(yàn)證中的荷載施加Fig.7 Load application in accuracy verification

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

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

以在方案3 結(jié)構(gòu)的剛度中心上施加荷載前后的響應(yīng)為例，若求解的剛度中心準(zhǔn)確無(wú)誤，則在剛度中心上施加作用力時(shí)應(yīng)無(wú)扭轉(zhuǎn)產(chǎn)生，但實(shí)際上仍然有-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%以下，認(rèn)為求解得出的剛度中心滿足要求。

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

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

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

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

表4 各方案的剛度中心坐標(biāo)及精確性驗(yàn)證Tab.4 Stiffness center coordinates of each scheme and accuracy verification

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

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

其中：

式中：［K］為剛度矩陣；xi、yi為第i個(gè)抗側(cè)構(gòu)件以荷載作用點(diǎn)為原點(diǎn)的坐標(biāo)值；∑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作用下的響應(yīng)進(jìn)行說(shuō)明。令該點(diǎn)在Fx作用下x方向的位移值為α、y方向的位移值為β，該點(diǎn)在Fy作用下x方向的位移值為ζ、y方向的位移值為ξ，各位移值的變化曲線按等比例繪制于圖10，各方案中各位移值的具體數(shù)據(jù)如表5所示。

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

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

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

3 結(jié) 論

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

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

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