地震作用下貯液結(jié)構(gòu)象足位置樣條插值分析方法

胡明祎 ,　張令心,　婁　宇 ,　陳　騮 ,　黃　偉

(1.中國電子工程設(shè)計院, 北京　100142; 2.北京市微振動環(huán)境控制工程技術(shù)研究中心, 北京　100086; 3.清華大學(xué) 土木水利工程學(xué)院,北京　100084; 4.中國地震局工程力學(xué)研究所,黑龍江 哈爾濱　150080; 5.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥　230009; 6.土木工程結(jié)構(gòu)與材料安徽省重點實驗室,安徽 合肥　230009)

地震作用下貯液結(jié)構(gòu)象足位置樣條插值分析方法

胡明祎1,2,3,張令心1,4,婁宇1,2,陳騮1,2,黃偉5,6

(1.中國電子工程設(shè)計院, 北京100142; 2.北京市微振動環(huán)境控制工程技術(shù)研究中心, 北京100086; 3.清華大學(xué) 土木水利工程學(xué)院,北京100084; 4.中國地震局工程力學(xué)研究所,黑龍江 哈爾濱150080; 5.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥230009; 6.土木工程結(jié)構(gòu)與材料安徽省重點實驗室,安徽 合肥230009)

摘要:貯液結(jié)構(gòu)的象足式破壞是一種典型的震害現(xiàn)象,為了有效地減輕或避免象足式震害,文章采用有限元方法(FEM)和數(shù)學(xué)插值方法對地震中貯液結(jié)構(gòu)產(chǎn)生象足破壞的位置進行了研究。首先將數(shù)值計算結(jié)果進行指數(shù)公式擬合,確定出快速判斷象足破壞位置的方法,然后對象足部位進行合理的結(jié)構(gòu)優(yōu)化設(shè)計以提高貯液結(jié)構(gòu)抗震性能。在分析過程中,利用FEM方法對貯液結(jié)構(gòu)地震響應(yīng)進行數(shù)值模擬,考慮了液固耦合效應(yīng),同時根據(jù)已知結(jié)構(gòu)響應(yīng)在某一區(qū)域收斂于某個極值的條件,采用樣條插值方法確定象足產(chǎn)生的位置,提出在有限計算資源情況下,利用樣條函數(shù)插值分析低精度FEM數(shù)值計算結(jié)果以獲取高精度結(jié)構(gòu)響應(yīng)的方法,并通過實例驗證了該方法的實用性。

關(guān)鍵詞:貯液結(jié)構(gòu);地震響應(yīng);有限元方法;插值

貯液結(jié)構(gòu)的地震破壞一般包括如下3種形式:① 傾覆彎矩過大引起結(jié)構(gòu)整體傾覆; ② 基底剪力過大引起結(jié)構(gòu)底部開裂、材料失效; ③環(huán)向力和軸向力共同作用引起結(jié)構(gòu)下部材料屈曲,亦稱象足屈曲。本文主要針對象足屈曲震害現(xiàn)象進行貯液結(jié)構(gòu)地震響應(yīng)優(yōu)化應(yīng)用方法研究。研究內(nèi)容是在有限計算資源條件下,尋求一種簡易計算方法,可以快速確定象足位置,據(jù)此對結(jié)構(gòu)進行優(yōu)化設(shè)計以提高貯液結(jié)構(gòu)抵抗象足屈曲的能力。

利用FEM進行數(shù)值計算時,保證有限元計算精度的常用措施是增加網(wǎng)格劃分密度,但在計算資源有限的情況下,特別是對于考慮液固耦合狀態(tài)非線性的貯液結(jié)構(gòu)地震反應(yīng)問題,增加網(wǎng)格密度會因計算機運算能力所限而降低求解效率。但是,如果在低密度網(wǎng)格情況下,已經(jīng)準(zhǔn)確獲悉結(jié)構(gòu)響應(yīng)在某一區(qū)域內(nèi)收斂于極大值或極小值,那么通過適當(dāng)?shù)靥岣呔W(wǎng)格劃分密度,采用數(shù)學(xué)插值的手段對該區(qū)域響應(yīng)值進行高精度插值,便可以有效地提高結(jié)構(gòu)有限元計算精度。因此,本文的技術(shù)問題便歸結(jié)于利用數(shù)學(xué)插值方法對一定精度條件下的有限元計算結(jié)果進行插值分析,從而獲取高精度結(jié)構(gòu)響應(yīng)結(jié)果。

1貯液結(jié)構(gòu)FEM數(shù)值模擬

貯液結(jié)構(gòu)有限元數(shù)值模擬工作已經(jīng)開展多年,目前有多種數(shù)值簡化計算模型。文獻[1]在前人研究的基礎(chǔ)上提出了質(zhì)量彈簧系統(tǒng)的儲液罐剛性簡化模型;文獻[2]深入研究了貯液結(jié)構(gòu)彈性變形對液固耦合效應(yīng)影響的重要性,推廣應(yīng)用液固耦合半解析半數(shù)值計算模型;文獻[3]利用有限元方法研究液固耦合振動問題;文獻[4] 提出關(guān)于彈性圓柱形儲液罐地震反應(yīng)解析解分析模型;文獻[5]采用附加質(zhì)量模型,建立了液固耦合邊界處的附加質(zhì)量矩陣以確保液固耦合方程有解條件的連續(xù)性;文獻[6]采用Sanders殼體有限元理論將液固耦合振動方程按線性降階處理建立了FEM模型;文獻[7]提出了考慮彈性變形的柔性罐壁單自由度質(zhì)點模型;文獻[8]提出了Haroun-Housner三分量模型?？傊?在材料模型可靠的條件下,FEM模型要比解析解模型能更加真實地模擬貯液結(jié)構(gòu)振動問題,隨著有限元理論的飛速發(fā)展,FEM模型在貯液結(jié)構(gòu)地震響應(yīng)研究中也得到了很好的應(yīng)用[9-10]。

貯液結(jié)構(gòu)地震響應(yīng)數(shù)值模擬的難點包括結(jié)構(gòu)材料理論屬性、液固耦合相互作用、結(jié)構(gòu)彈性變形對計算結(jié)果產(chǎn)生的影響等。本文采用FEM方法建立貯液結(jié)構(gòu)數(shù)值模型時,假設(shè)貯液振動為小幅振動,貯液為無旋無黏性流體,有限元模型中,貯液結(jié)構(gòu)采用shell單元,具有6個自由度,材料為雙線性彈塑性模型,貯液采用勢流體單元,具有3個平動自由度。貯液結(jié)構(gòu)有限元模型剖面圖如圖1a所示,采用的地震荷載為ELcentro地震動,加速度峰值為341 cm/s2,加載時間為20 s,步長取為0.02 s。圖1b所示為貯液結(jié)構(gòu)地震響應(yīng)象足效應(yīng)圖。

(a)　有限元模型剖面圖　　(b)　地震響應(yīng)象足效應(yīng)圖

2樣條函數(shù)插值分析

樣條函數(shù)插值已經(jīng)被廣泛地應(yīng)用到流線型形狀設(shè)計問題中,其應(yīng)用的目標(biāo)狀態(tài)是保證有限插值控制的精確性和收斂性。相對Lagrange插值函數(shù)和三次Hermite插值而言,樣條函數(shù)插值優(yōu)勢是既能有效地避免產(chǎn)生邊界Runge現(xiàn)象,也能成功解決平滑化問題。樣條插值的函數(shù)基為三次多項式函數(shù),而且其一階、二階導(dǎo)數(shù)在插值控制點都是連續(xù)的。

假設(shè)F″(xj)=Mj,j=1,2,3,…,n,由于F(x)是二階光滑的分段三次多項式,于是F″(x)是分段線性連續(xù)函數(shù),可以在x∈[xj,xj+1]上假設(shè)F″(x)函數(shù)形式為:

(1)

其中,x∈[xj,xj+1]。對(1)式進行兩次積分可得:

(2)

建立插值條件(3)式、(4)式,可以求解出系數(shù)Aj和Bj,即

(3)

(4)

(5)

(6)

假設(shè)插值函數(shù)一階導(dǎo)數(shù)連續(xù),F′(xj+0)=F′(xj-0),同時采用第2類邊界條件F″(a)=F″(b)=0,聯(lián)立可以建立如下線性方程組:

(7)

其中,λj=(xj+1-xj)/(xj+1-xj-1);μj=1-λj;f[xj-1xjxj+1]為點xj-1、xj、xj+1的二階差商。求解(7)式進而可以得到Mj值,返回帶入(5)式、(6)式可以求解系數(shù)Aj和Bj,再帶入(2)式最后進行插值求解。

需要注意的是插值求解過程中,也可以根據(jù)情況調(diào)整邊界條件,但是通過多次調(diào)試對比,筆者認(rèn)為文中給出插值邊界條件較為穩(wěn)定,與精細網(wǎng)格有限元模型計算結(jié)果相比象足位置計算結(jié)果誤差低于10-5,因此本文的邊界條件采用第2類邊界條件,具體取值規(guī)定為二階導(dǎo)數(shù)為零。

從理論角度而言,如果所有插值區(qū)域具有同樣的屬性,并且嚴(yán)格在插值點具有收斂特征,那么樣條函數(shù)插值是一種很好的數(shù)值逼近方法。本文中,假設(shè)研究的對象材料分布連續(xù)并且均勻,其地震響應(yīng)沿結(jié)構(gòu)壁高度分布也具有連續(xù)并且收斂的特征,因此可以采用樣條函數(shù)插值方法對結(jié)構(gòu)最大響應(yīng)位置進行插值逼近求解。

圖2所示為基于FEM計算結(jié)果的插值分析程序流程圖。

圖2　基于FEM計算結(jié)果的插值分析程序流程圖

根據(jù)上述方法,按照圖2中的流程圖編制插值分析程序,可以結(jié)合貯液結(jié)構(gòu)地震反應(yīng)FEM數(shù)值計算結(jié)果進行插值分析。分析之前要確定插值分析的數(shù)據(jù)使用條件,即保證插值控制點的有效性。通過前面FEM數(shù)值模擬,本文選用的模型網(wǎng)格劃分密度為0.05 m,所以結(jié)構(gòu)壁最大反應(yīng)及象足高度位置的計算精度也為0.05 m。為了說明該精度對于插值控制點的有效性,通過細化網(wǎng)格密度范圍為0.05～0.001 m,重新計算分析,可以證實網(wǎng)格密度為0.05 m和低于0.05 m的結(jié)構(gòu)響應(yīng)結(jié)果誤差基本在1.5%以內(nèi),表1所列為不同網(wǎng)格密度對應(yīng)象足位置計算精度比較,如果按照允許誤差1%考慮,在0.05 m的網(wǎng)格劃分密度條件下,FEM數(shù)值計算結(jié)果針對于象足高度位置的精度可以達到0.01 m,這種結(jié)果同時也表明按0.05 m的網(wǎng)格密度計算出的象足位置可以保證插值控制點的有效性。

表1　不同網(wǎng)格密度對應(yīng)象足位置計算精度比較

注:斜線后為貯液結(jié)構(gòu)地震響應(yīng)計算結(jié)果相對網(wǎng)格密度為0.05 m情況下的誤差率;黑體為有效精度確定位置。

基于上述條件,對貯液結(jié)構(gòu)的FEM計算結(jié)果進行插值分析。插值結(jié)果檢驗?zāi)繕?biāo)有2項,即插值結(jié)果穩(wěn)定性和插值曲線平滑化。為了說明樣條函數(shù)插值結(jié)果的優(yōu)勢,本文采用Lagrange高階插值函數(shù)、分段Lagrange拋物線插值函數(shù)和三次樣條插值函數(shù)對FEM計算結(jié)果進行插值結(jié)果對比,這里給出其中4種工況的計算分析結(jié)果,如圖3所示。圖3中D為貯液結(jié)構(gòu)直徑;H為貯液結(jié)構(gòu)高度;Hw為貯液高度；φ為有效精度；A、B、C、D、E、F為插值控制點。由圖3a、圖3b可以清晰看到對于不同貯液量的細高貯液結(jié)構(gòu)而言,Lagrange高階插值函數(shù)在邊界處產(chǎn)生顯著的Runge現(xiàn)象,引起邊界處插值不穩(wěn)定性,分段Lagrange拋物線插值會有效降低Runge現(xiàn)象程度,但是整體曲線平滑性較差,三次樣條函數(shù)插值結(jié)果很好地避免了上述問題,并且有效地保證了精度水平為0.01 m。

圖3c、圖3d研究對象為相同貯液量不同高度的矮胖貯液結(jié)構(gòu),可以看到這種結(jié)構(gòu)地震響應(yīng)Lagrange高階插值和三次樣條插值結(jié)果基本一致,而且Runge現(xiàn)象不顯著。此外,通過比較,樣條函數(shù)插值最終的結(jié)果精度保證在0.01 m水平,圖3中有效精度誤差分別為0.14、0.11、0.08、0.13 m,最大相對精度誤差是圖3a，為28%,這意味著通過插值分析確定的象足位置高度值精度比FEM直接計算值提高28%。這種計算精度的有效提升對于象足位置的地震響應(yīng)優(yōu)化設(shè)計而言具有重要意義。

圖3　插值比較曲線

3插值結(jié)果擬合分析

通過前面針對FEM數(shù)值計算結(jié)果的插值精度分析,可以確定樣條函數(shù)插值獲取的象足位置插值結(jié)果是有效的,下面對不同工況下的插值結(jié)果進行統(tǒng)計分析,確定象足位置與D/H值和Hw/H值的關(guān)系,從而為后面的優(yōu)化工作提供依據(jù)。圖4所示為D/H值取0.5和1.0對應(yīng)的象足高度值Hmax隨Hw/H的等效坐標(biāo)值的變化曲線

通過圖4,根據(jù)插值結(jié)果Hmax曲線指數(shù)衰減關(guān)系,可以進行指數(shù)函數(shù)擬合,圖4a、圖4b對應(yīng)函數(shù)形式分別為(8)式、(9)式:

(8)

(9)

由圖5可以發(fā)現(xiàn)隨著D/H值和Hw/H不斷增大,Hmax值也在增大,這個規(guī)律在貯液結(jié)構(gòu)抗震設(shè)計中非常重要,通過這種插值擬合的方法將不同設(shè)計工況下的Hmax取值量化有利于實際應(yīng)用。

采用這種方式,可以將Hmax值隨D/H值和Hw/H值同時變化的曲線繪制在一起,可以得到圖6所示的有水剛結(jié)頂蓋Hmax隨D/H變化曲線和有水浮頂蓋Hmax隨D/H變化曲線,兩者的工況區(qū)別為貯液結(jié)構(gòu)蓋頂構(gòu)筑形式不同。

圖4　Hmax隨Hw/H變化曲線

圖5　Hmax隨D/H變化曲線

圖6　Hw/H=0.9時Hmax隨D/H變化曲線

由圖6可以發(fā)現(xiàn)Hmax隨D/H取值在貯液量較少時變化較小,貯液量較多時變化較大,也說明了液固耦合作用對結(jié)構(gòu)的響應(yīng)隨貯液增加而變大,同時也發(fā)現(xiàn)蓋頂形式對象足位置的高度影響不大。

(10)

注:D/H值的范圍為0.001～0.01,由于篇幅限制這里只給出部分值;實際應(yīng)用中A、B、C值為擬合系數(shù),可以按照表2進行線性插值取值。

4優(yōu)化應(yīng)用分析

根據(jù)上述基于對貯液結(jié)構(gòu)地震反應(yīng)FEM計算結(jié)果插值分析而擬合獲取的象足位置計算公式及系數(shù)表格的確定方法,將該表格擴展到不同高度、厚度、形狀或者其他影響因素,可以得到關(guān)于象足高度快速確定的設(shè)計參照系列表格,利用這些表格可以快速確定對應(yīng)象足部位計算函數(shù)中的系數(shù)A、B、C值,從而獲得不同抗震設(shè)計方案下的象足高度位置,進而對原結(jié)構(gòu)形式進行優(yōu)化設(shè)計。

具體的優(yōu)化設(shè)計思路是在象足部位對貯液結(jié)構(gòu)增加環(huán)形肋梁,通過這種實際的應(yīng)用分析,可以快速對原設(shè)計方案進行抗震性能驗算和確定對象足部位進行結(jié)構(gòu)加固和功能補足。對于具體實例,圖7a所示為FEM優(yōu)化模型,圖7b所示為優(yōu)化后貯液結(jié)構(gòu)地震響應(yīng)象足效應(yīng)圖。在實際情況下,具體的設(shè)計方式也可以采用內(nèi)部加環(huán)形肋梁,或者外部增加環(huán)向加勁肋,其增加的位置都是通過(10)式和設(shè)計參照表格得到的。

(a)　有限元模型剖面圖　　　　　　(b)　地震響應(yīng)示意圖

圖8所示為優(yōu)化前和優(yōu)化后的結(jié)構(gòu)側(cè)向最大位移隨D/H變化的曲線對比圖。

圖8　優(yōu)化設(shè)計前、后Hmax隨D/H變化曲線

通過圖8可以看出,優(yōu)化后的FEM計算結(jié)果要明顯小于優(yōu)化前計算結(jié)果,其中對于Hw/H=0.8的工況下,結(jié)構(gòu)最大側(cè)向位移減小了40%。

5結(jié)論

本文針對考慮液固耦合相互作用的圓柱形貯液結(jié)構(gòu)的象足破壞問題,在計算資源有限的情況下,利用FEM方法對貯液結(jié)構(gòu)地震響應(yīng)進行數(shù)值模擬,并利用樣條插值手段對低精度FEM計算結(jié)果進行分析和比較,提出了利用樣條函數(shù)插值分析低精度FEM數(shù)值計算結(jié)果獲取高精度結(jié)構(gòu)響應(yīng)的方法,同時采用該方法計算結(jié)果,對貯液結(jié)構(gòu)象足出現(xiàn)位置進行了指數(shù)公式擬合,從而對結(jié)構(gòu)進行優(yōu)化設(shè)計,改善其地震響應(yīng)效果,最后對優(yōu)化前后的貯液結(jié)構(gòu)地震響應(yīng)進行了對比,驗證了該方法的實用性。具體結(jié)論如下:

(1) 在進行插值分析時,提出采用第2類插值邊界條件F″(a)=F″(b)=0的方式進行,該邊界條件可以很好地保證插值邊界連續(xù)且平滑。

(2) 將FEM方法和插值方法相結(jié)合,編制了插值分析程序,通過確定FEM網(wǎng)格密度對貯液結(jié)構(gòu)地震響應(yīng)的影響,保證了FEM計算結(jié)果可用來進行插值分析的條件。

(3) 對象足部位確定公式,采用基于FEM和插值結(jié)果的指數(shù)衰減公式,繪制指數(shù)參數(shù)設(shè)計參照表格,這種方法可以應(yīng)用于結(jié)構(gòu)抗震驗算中,并且地震響應(yīng)優(yōu)化效果較好。

雖然該方法可以在一定條件下有效地提高低精度FEM計算精度,但是在有限計算資源條件下,優(yōu)化結(jié)構(gòu)計算的精度劃分區(qū)域,如局部采用FEM子結(jié)構(gòu)或者局部進行插值分析的方法可以更好地提高FEM計算結(jié)果的精度,本文為此提供了有效的理論分析依據(jù)。

[參考文獻]

[1]Housner G W. Dynamic pressures on accelerated fluid containers[J]. Bulletin of the SeismologicalSociety of America 1957,47:15-35.

[2]Veletsos A S. Seismic effects in flexible liquid storage tanks[C]//Proc 5th World Conference on Earthquake Engineering, Rome, Italy, Vol 1,1973:630-639.

[3]Edward N W.A procedure for dynamic analysis of thin walled cylindrical liquid storage tanks subjected to lateral ground motions[D]. Univ of Michigan, 1969.

[4]Tedesco J W, Kostem C N,Kalnins A. Free vibration analysis of cylindrical liquid storage tanks[J]. Comput Struct,1987,26 (6):957-964.

[5]Rajasankar J, Iyer N R, Appa R T V S R.A new 3-D finite element model to evaluate added mass for analysis of fluid structure interaction problems[J]. Int J Numer Methods Eng,1993, 36:997-1012.

[6]Nash W A. Response of liquid storage tanks to seismic motion[C]// Koiter W T,Mikhailov G K. Theory of Shells. North-Holland Publishing Company, 1980:393-403.

[7]Veletsos A S, Yang J Y. Earthquake response of liquid storage tanks[C]// Advances in Civil Engineering through Engineering Mechanics:Proceedings Annual EMD Specialty Conference ASCE, Raleigh, N C, USA,1977:1-24.

[8]Haroun M A, Housner G W. Seismic design of liquid storage tanks[J]. ASCE Journal of Technical Councils,1981,107 (1):191-207.

[9]Hu Mingyi,Zhang Lingxin,Liu Jieping,et al. Structural frequency topology optimization in seismic design based on FEM method[C]//The 14th WCEE,Beijing,China,2008:12-17.

[10]Hu Mingyi,Li Wei, Hou Yuxin.Numerical study on seismic response analysis of reinforced concrete aqueduct strucure[C]//International Symposium on Innovation & Sustainability of Structures in Civil Engineering,Shanghai, China,2007:28-30.

(責(zé)任編輯張镅)

Spline interpolation methods to confirm elephant-foot position of maximum deformation of liquid container during earthquake

HU Ming-yi1, 2, 3, ZHANG Ling-xin1, 4, LOU Yu1, 2, CHEN Liu1, 2, HUANG Wei5, 6

(1.China Electronics Engineering Design Institute, Beijing 100142, China; 2.Micro-vibration Environmental Control Engineering Technology Research Center of Beijing, Beijing 100086, China; 3.School of Civil Engineering, Tsinghua University, Beijing 100084, China; 4.Institute of Engineering Mechanics, China Earthquake Administration, Harbin 150080, China; 5.School of Civil and Hydraulic Engineering, Hefei University of Technology, Hefei 230009, China; 6.Anhui Key Laboratory of Structure and Materials in Civil Engineering, Hefei 230009, China)

Abstract:Damage of liquid container may be great due to elephant-foot-type buckling during earthquake in general. In this paper,the finite element method(FEM) and spline interpolation method are used to analyze the elephant-foot position of the container structure during earthquake.The numerical calculation results are fitted by exponential formula,and the method for rapidly evaluating the elephant-foot position is obtained.On this basis, the reasonable optimal design of the elephant-foot position is maded to improve the structure’s seismic behavior.In the whole analysis, the FEM is used to simulate the structure’s seismic response with consideration of fluid-structure interaction, and based on the finite element analysis result that the structure’s responses converge to an extreme value within certain region, the elephant-foot position is determined by spline interpolation method. Although the computing resource is limited, high precision structural response can be gotten by the presented method.The practicality of this method is validated by an example.

Key words:liquid container; seismic response; finite element method(FEM); interpolation

收稿日期：2015-08-03；修回日期：2015-09-30

基金項目：北京市科學(xué)技術(shù)委員會資助項目(Z131101002813083)

作者簡介：胡明祎(1982-),男,河南固始人,博士,中國電子工程設(shè)計院高級工程師.

doi:10.3969/j.issn.1003-5060.2016.05.020

中圖分類號:TU352

文獻標(biāo)識碼：A

文章編號：1003-5060(2016)05-0677-07