亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        考慮土-結(jié)構(gòu)動(dòng)力相互作用的冷卻塔地震響應(yīng)分析

        2017-01-10 08:14:38張俊發(fā)陳厚群
        振動(dòng)與沖擊 2016年23期
        關(guān)鍵詞:結(jié)構(gòu)模型

        陶 磊, 張俊發(fā), 陳厚群

        (1.西安理工大學(xué) 水利水電學(xué)院,西安 710048;2.西安理工大學(xué) 土木建筑工程學(xué)院,西安 710048;3.中國(guó)水利水電科學(xué)研究院 工程抗震研究中心,北京 100048)

        考慮土-結(jié)構(gòu)動(dòng)力相互作用的冷卻塔地震響應(yīng)分析

        陶 磊1, 張俊發(fā)2, 陳厚群3

        (1.西安理工大學(xué) 水利水電學(xué)院,西安 710048;2.西安理工大學(xué) 土木建筑工程學(xué)院,西安 710048;3.中國(guó)水利水電科學(xué)研究院 工程抗震研究中心,北京 100048)

        大型高聳鋼筋混凝土冷卻塔屬于典型的風(fēng)敏感型結(jié)構(gòu),近年來在地震、飛機(jī)撞擊、爆破等極端外部作用下動(dòng)力響應(yīng)研究也得到了工程界的廣泛關(guān)注。土-結(jié)構(gòu)動(dòng)力相互作用效應(yīng)(Soil Structure Interaction, SSI)對(duì)于大壩、橋梁等類型工程結(jié)構(gòu)地震響應(yīng)的影響研究成果較為豐富,而對(duì)冷卻塔結(jié)構(gòu)體系的影響程度方面很少涉及。依據(jù)彈性波動(dòng)理論結(jié)合有限單元法,建立和推導(dǎo)了考慮土-結(jié)構(gòu)動(dòng)力相互作用的三維黏彈性人工邊界模型和公式,并通過半空間自由場(chǎng)模型驗(yàn)證該地震動(dòng)輸入方法的準(zhǔn)確性。以國(guó)內(nèi)實(shí)際工程項(xiàng)目——某火電廠大型冷卻塔為研究背景,以通用有限元程序ANSYS為平臺(tái),分別建立了冷卻塔剛性地基模型、無質(zhì)量地基模型和黏彈性人工邊界模型,開展模態(tài)分析及彈性時(shí)程分析,研究不同計(jì)算模型相應(yīng)的動(dòng)力特性及內(nèi)力變化,探討了土-結(jié)構(gòu)動(dòng)力相互作用的影響規(guī)律。研究結(jié)果表明:考慮剛性地基,冷卻塔結(jié)構(gòu)體系自振頻率分布十分密集,絕大多數(shù)振型為環(huán)向諧波和子午向諧波組合的局部振型;考慮彈性地基后,結(jié)構(gòu)的自振頻率略有降低,整體振型較早出現(xiàn)。通過時(shí)程分析可知,采用黏彈性人工邊界模型,考慮無限地基輻射阻尼效應(yīng),與剛性地基模型相比,塔筒的絕對(duì)加速度最大值降低43.4%,塔筒沿子午向彎矩軸力幅值降低約50%,而環(huán)向內(nèi)力卻均顯著提高,X支柱內(nèi)力幅值降低約20%~50%.因此,在進(jìn)行冷卻塔地震響應(yīng)分析時(shí),土-結(jié)構(gòu)動(dòng)力相互作用的影響不可忽視。

        冷卻塔結(jié)構(gòu);地震響應(yīng)分析;土-結(jié)構(gòu)動(dòng)力相互作用;黏彈性人工邊界;自由場(chǎng)分析;地震動(dòng)輸入

        冷卻塔是電廠中用于散熱冷卻的重要電力構(gòu)筑物,屬火電站或核電站的重要組成部分,常見結(jié)構(gòu)型式為一種典型的高聳軸對(duì)稱旋轉(zhuǎn)雙曲線空間薄殼結(jié)構(gòu)。由于冷卻塔體型龐大、經(jīng)受荷載復(fù)雜及其使用功能在電廠中的重要地位,一直是工程界關(guān)注和重點(diǎn)研究的對(duì)象。近年來,隨著電力行業(yè)需求的不斷增加,火電、核電的單機(jī)容量隨之增加,對(duì)冷卻塔提出了更高的要求,其發(fā)展呈現(xiàn)高度高(目前國(guó)內(nèi)最高的冷卻塔已超過200 m)、跨度大、下部支撐結(jié)構(gòu)高的趨勢(shì)[1],這些超高的冷卻塔的塔高已經(jīng)超過相關(guān)規(guī)范規(guī)定的限值。對(duì)建造于高烈度地區(qū)的大型冷卻塔在地震作用下的動(dòng)力分析和抗震設(shè)計(jì)已經(jīng)得到了設(shè)計(jì)人員的廣泛關(guān)注,《構(gòu)筑物抗震設(shè)計(jì)規(guī)范》(GB50191—2012)中規(guī)定[2]:“冷卻塔系在地震時(shí)使用功能不能中斷或需盡快恢復(fù)的構(gòu)筑物,按其使用功能的重要性分類,應(yīng)屬乙類抗震設(shè)防類別的構(gòu)筑物?!币虼藢?duì)地震等極端外部荷載作用下,冷卻塔結(jié)構(gòu)的強(qiáng)度和穩(wěn)定性開展深入研究具有重要意義。

        國(guó)內(nèi)外學(xué)者對(duì)冷卻塔地震反應(yīng)分析開展了一系列研究,取得了較為豐富的成果。GUPTA等[3]認(rèn)為對(duì)于冷卻塔抗震設(shè)計(jì),只取水平向運(yùn)動(dòng)即可。WOLF[4]將冷卻塔離散為有限元模型,環(huán)梁、支柱與軸對(duì)稱殼單元合成為動(dòng)態(tài)剛度矩陣,地基作為獨(dú)立頻率參數(shù)的彈簧和阻尼,進(jìn)行了時(shí)域時(shí)程分析。NASIR等[5]采用有限元法分析雙曲線殼體的自由振動(dòng)和受地震激勵(lì)的響應(yīng),研究了雙曲線殼體的厚度、高度和曲率對(duì)動(dòng)態(tài)響應(yīng)的影響。SABOURI GHOMI等[6]研究了冷卻塔受地震激勵(lì)時(shí)的動(dòng)態(tài)特性,通過實(shí)測(cè)水平向和豎向加速度進(jìn)行線性和非線性時(shí)程分析,評(píng)估了支柱中塑性鉸的位置對(duì)冷卻塔穩(wěn)定性的影響。柯世堂首次對(duì)冷卻塔進(jìn)行了考慮行波效應(yīng)的多點(diǎn)激勵(lì)分析,開展了模態(tài)分析、振型分解反應(yīng)譜和彈性及考慮材料和幾何非線性的動(dòng)力彈塑性時(shí)程分析。葉志浩等[7]建立基于混凝土損傷塑性模型,利用ABAQUS建立了自然通風(fēng)冷卻塔的三維有限元模型,通過薄弱部位支柱部分節(jié)點(diǎn)拉伸損傷因子的輸出統(tǒng)計(jì)進(jìn)行分析,研究不同跨間距和不同支柱對(duì)數(shù)的人字形支柱對(duì)冷卻塔結(jié)構(gòu)抗震性能的影響。綜上所述,以上學(xué)者關(guān)于冷卻塔結(jié)構(gòu)地震反應(yīng)分析研究的側(cè)重點(diǎn)在于抗震分析方法、分析模型的建立手段、結(jié)構(gòu)體型的影響以及結(jié)構(gòu)抗震性能的評(píng)估等。以往的文獻(xiàn)中,開展冷卻塔結(jié)構(gòu)地震響應(yīng)分析時(shí),往往僅研究冷卻塔環(huán)基、支撐及塔筒結(jié)構(gòu),采用剛性地基假定,未考慮基礎(chǔ)變形對(duì)結(jié)構(gòu)地震反應(yīng)的影響;或用彈簧單元模擬地基的彈性作用,作為封閉系統(tǒng)的振動(dòng)問題來求解。從波動(dòng)理論上分析,考慮地基土與冷卻塔結(jié)構(gòu)動(dòng)力相互作用的時(shí)候,應(yīng)作為開放系統(tǒng)的波動(dòng)問題進(jìn)行研究:來自地基無限域的入射波,經(jīng)過地基和幾何邊界反射和折射產(chǎn)生相對(duì)于人工邊界的外行波,最終要靠人工邊界吸收截?cái)噙吔缟贤庑胁??!稑?gòu)筑物抗震設(shè)計(jì)規(guī)范》GB50191—2012亦給出了下列規(guī)定:“塔筒的地震作用計(jì)算宜計(jì)及地基與上部結(jié)構(gòu)的相互作用,計(jì)算時(shí)應(yīng)采用土的動(dòng)力參數(shù)。”,但對(duì)于如何建立考慮地基-結(jié)構(gòu)的相互作用的計(jì)算模型等一系列問題未給出詳細(xì)說明。因此,對(duì)大型冷卻塔結(jié)構(gòu)地震動(dòng)輸入機(jī)制及土-結(jié)構(gòu)動(dòng)力相互作用影響效果做進(jìn)一步深入分析是十分必要的。

        也有少量文獻(xiàn)開展了考慮土-結(jié)構(gòu)動(dòng)力相互作用的冷卻塔地震響應(yīng)分析研究。其中,HORR等[8]根據(jù)復(fù)雜空間單元法和分?jǐn)?shù)階微積分理論提出了一種混合復(fù)雜解耦的空間方法用于分析冷卻塔的土-結(jié)構(gòu)動(dòng)力相互作用。高標(biāo)等[9]考慮了冷卻塔與地基的動(dòng)力相互作用,分析了樁基、地基土阻尼比對(duì)冷卻塔結(jié)構(gòu)抗震性能的影響。李輝等[10]建立了考慮土-結(jié)構(gòu)動(dòng)力相互作用的有限元模型,開展動(dòng)力特性及反應(yīng)譜分析,并與剛性地基結(jié)果進(jìn)行了對(duì)比。以上學(xué)者雖探討了土-結(jié)構(gòu)動(dòng)力相互作用的影響,但未能從波動(dòng)角度出發(fā),對(duì)土-結(jié)構(gòu)動(dòng)力相互作用效應(yīng)對(duì)冷卻塔的地震響應(yīng)影響規(guī)律進(jìn)行深入分析。

        鑒于大型冷卻塔結(jié)構(gòu)體系考慮土-結(jié)構(gòu)動(dòng)力相互作用的地震反應(yīng)分析方面研究尚不深入,本文依據(jù)波動(dòng)理論結(jié)合有限元法,建立三維黏彈性人工邊界的冷卻塔有限元模型,模擬無限地基輻射阻尼效應(yīng),以考慮土-結(jié)構(gòu)動(dòng)力相互作用的影響。在人工邊界上將波動(dòng)分解為自由波和散射波,并將輸入地震波動(dòng)轉(zhuǎn)化為作用于人工邊界上的等效荷載以實(shí)現(xiàn)波動(dòng)輸入。以國(guó)內(nèi)某火電廠采用間接空冷系統(tǒng)的大型冷卻塔為研究背景,建立冷卻塔剛性地基模型、無質(zhì)量地基模型和黏彈性人工邊界模型,通過模態(tài)分析及動(dòng)力時(shí)程分析,探討考慮土-結(jié)構(gòu)動(dòng)力相互作用對(duì)冷卻塔結(jié)構(gòu)動(dòng)力響應(yīng)的影響規(guī)律。

        1 土-結(jié)構(gòu)動(dòng)力相互作用理論

        土-結(jié)構(gòu)動(dòng)力相互作用(Soil Structure Interaction, SSI)是指地震波通過場(chǎng)地土介質(zhì)傳播結(jié)構(gòu)底部,激起結(jié)構(gòu)發(fā)生振動(dòng),結(jié)構(gòu)體系產(chǎn)生的慣性力如同新的震源反過來作用與場(chǎng)地土,這種現(xiàn)象稱為土-結(jié)構(gòu)動(dòng)力相互作用。從波動(dòng)理論上講,土-結(jié)構(gòu)動(dòng)力相互作用就是波動(dòng)在結(jié)構(gòu)-地基系統(tǒng)內(nèi)進(jìn)行傳播時(shí)引起的結(jié)構(gòu)和地基的動(dòng)力反應(yīng)問題。

        土-結(jié)構(gòu)動(dòng)力相互作用的理論分析方法有多種劃分方法[11-12],按求解域可分為時(shí)域法和頻域法;按求解方法可分為解析法、半解析法、數(shù)值法和等效法;按結(jié)構(gòu)系統(tǒng)劃分為子結(jié)構(gòu)法和整體法。子結(jié)構(gòu)法中將結(jié)構(gòu)與近域地基用有限元進(jìn)行離散,遠(yuǎn)場(chǎng)地基則通過近場(chǎng)和遠(yuǎn)場(chǎng)交界面上,用彈簧和阻尼器單元代替地基土的影響,而計(jì)算參數(shù)通過地基的動(dòng)力阻抗函數(shù)G(ia0)來確定。G(ia0)取決于基礎(chǔ)的尺寸、形狀、地基介質(zhì)的力學(xué)參數(shù)及強(qiáng)迫振動(dòng)頻率等。對(duì)于彈性半空間上剛性無質(zhì)量圓板的地基阻抗函數(shù)形式如下:

        G(ia0)=GR(a0,ν)+iGI(a0,ν)

        (1)

        式中:上標(biāo)R和I分別表示實(shí)部和虛部,a0是無量綱頻率,定義為

        (2)

        式中:R代表圓板半徑,Vs代表均勻半空間材料的剪切波速,ν代表泊松比。

        地基阻抗是一個(gè)復(fù)數(shù),實(shí)部體現(xiàn)了體系的剛度特征和慣性性質(zhì),而虛部則代表了體系的能量損失,并且實(shí)部和虛部均是頻率相關(guān)的,即地基阻抗具有頻率相關(guān)性。對(duì)于半空間上的剛性無質(zhì)量圓板,VELETSOS等[13],LUCO等[14]給出了確定阻抗函數(shù)的經(jīng)驗(yàn)公式。賀廣零等[15]通過在基礎(chǔ)與土體的交界面上設(shè)置豎向、水平、搖擺、扭轉(zhuǎn)等方向的彈簧-阻尼器單元實(shí)現(xiàn)了考慮土-結(jié)構(gòu)動(dòng)力相互作用情況下的風(fēng)力發(fā)電機(jī)結(jié)構(gòu)系統(tǒng)的動(dòng)力反應(yīng)分析。這一模型并不適用于模擬大型冷卻塔結(jié)構(gòu)考慮土-結(jié)構(gòu)動(dòng)力相互作用情況下的動(dòng)力響應(yīng)分析問題。一方面,冷卻塔結(jié)構(gòu)本身的動(dòng)力特性較風(fēng)力發(fā)電機(jī)結(jié)構(gòu)更為復(fù)雜;另一方面,實(shí)際的地基土層分布情況很復(fù)雜,很少能直接簡(jiǎn)化為彈性半空間。

        早期研究大多將結(jié)構(gòu)簡(jiǎn)化為置于彈性半空間的表面之上,采用WINKLER模型,簡(jiǎn)化的兩彈簧地基模型以及SR模型等。隨著計(jì)算理論及計(jì)算技術(shù)的發(fā)展,開展更為復(fù)雜、真實(shí)的考慮土-結(jié)構(gòu)動(dòng)力相互作用的三維有限元分析已經(jīng)成為可能。整體法將廣義結(jié)構(gòu)和地基組成的系統(tǒng)直接進(jìn)行整體有限元分析的方法,對(duì)于無限地基,需要引入人工邊界條件來近似模擬無限地基輻射阻尼效應(yīng)。目前工程抗震分析常用的人工邊界有透射邊界和黏彈性人工邊界。透射邊界精度較高,但存在高頻振蕩和漂移失穩(wěn)問題,且在有限元程序中處理不便。黏彈性人工邊界是在人工截?cái)嗟倪吔缟显O(shè)置彈簧和阻尼器,以吸收外行波能量,且能模擬半無限地基的彈性恢復(fù)能力。

        2 三維黏彈性人工邊界地震動(dòng)輸入

        三維黏彈性人工邊界(3D Viscous-Spring Artificial Boundary)是實(shí)現(xiàn)土-結(jié)構(gòu)動(dòng)力相互作用分析的一種直接有限元方法,其主要核心思想是將波動(dòng)的散射問題轉(zhuǎn)化為波源問題的方法來實(shí)現(xiàn),即通過在人工邊界上施加等效荷載實(shí)現(xiàn)波動(dòng)輸入,其本質(zhì)是一種應(yīng)力型局部人工邊界,基于波場(chǎng)分離技術(shù),將總波場(chǎng)分解為自由波場(chǎng)和散射波場(chǎng)。為了實(shí)現(xiàn)波動(dòng)輸入,首先需通過SHAKE91等計(jì)算程序計(jì)算得到整個(gè)自由場(chǎng)的應(yīng)力、速度和位移時(shí)程。何建濤等[18]對(duì)地震動(dòng)輸入公式進(jìn)一步推導(dǎo),將自由場(chǎng)應(yīng)力求解轉(zhuǎn)化為自由場(chǎng)速度求解,簡(jiǎn)化了黏彈性人工邊界輸入過程。

        2.1 三維黏彈性人工邊界推導(dǎo)

        根據(jù)彈性波動(dòng)理論,球坐標(biāo)系中球面膨脹波(P波)的波動(dòng)方程為:

        (3)

        式中:φ為位移勢(shì)函數(shù),cp為介質(zhì)的P波波速,R為徑向坐標(biāo)。上式的通解為:

        (4)

        式中:f(·)和g(·)為任意函數(shù),分別表示外行波和內(nèi)聚波。

        人工邊界處只考慮外行波,則垂直邊界處的法向應(yīng)力和位移滿足:

        (5)

        式中:G為剪切模量,ρ為介質(zhì)密度。

        為了建立人工邊界,將無限域連續(xù)介質(zhì)截?cái)?,在截?cái)嗵幵O(shè)置彈簧-阻尼器-集中質(zhì)量系統(tǒng),如圖1所示,通過推導(dǎo)可知,在人工邊界結(jié)點(diǎn)應(yīng)力與位移滿足微分方程:

        (6)

        式中:M,C,K分別為集中質(zhì)量,阻尼器的阻尼系數(shù),彈簧的剛度系數(shù)。比較(5)、(6)兩式可知,各物理元件的參數(shù)分別為:

        (7)

        同樣可以推導(dǎo)對(duì)于剪切波(S波)彈簧、阻尼器物理元件參數(shù)為:

        (8)

        然而式(7)、式(8)并非黏彈性人工邊界參數(shù)值的唯一形式。由于構(gòu)建思想的存在一定區(qū)別,眾多學(xué)者給出幾種彈簧、阻尼器不同參數(shù)形式,經(jīng)過大量數(shù)值計(jì)算,黏彈性人工邊界具有良好的魯棒性,都能近似處理一般的近場(chǎng)波動(dòng)問題。

        圖1 法向人工邊界等效系統(tǒng)

        2.2 地震動(dòng)輸入方法

        已知入射波在人工邊界結(jié)點(diǎn)B(xB,yB,zB)產(chǎn)生的三向?yàn)槲灰品謩e為u0(xB,yB,zB,t),v0(xB,yB,zB,t),w0(xB,yB,zB,t),法向和切向應(yīng)力分別為σ0(xB,yB,zB,t)和τ0(xB,yB,zB,t)。黏彈性人工邊界上的彈簧和阻尼器產(chǎn)生的法向和切向應(yīng)力分別為fBN和fBT,設(shè)人工邊界結(jié)點(diǎn)B的法向和切向等效應(yīng)力分別為FBN和FBT,等效結(jié)點(diǎn)荷載分別為PBN和PBT。

        三向位移和應(yīng)力必須滿足以下關(guān)系:

        (9)

        (10)

        人工邊界結(jié)點(diǎn)B的應(yīng)力可以表示為:

        (11)

        而彈簧阻尼器產(chǎn)生的應(yīng)力為:

        (12)

        施加在人工邊界結(jié)點(diǎn)B的等效結(jié)點(diǎn)荷載為:

        (13)

        綜上所述,通過彈性波動(dòng)理論結(jié)合有限元法,推導(dǎo)了三維黏彈性人工邊界彈簧阻尼器法向和切向的剛度系數(shù)和阻尼系數(shù)公式;并給出了用于自由場(chǎng)地震動(dòng)輸入的等效結(jié)點(diǎn)荷載公式。

        3 黏彈性人工邊界地震動(dòng)輸入方法正確性驗(yàn)證

        以有限元程序ANSYS為平臺(tái),通過APDL語言編制命令實(shí)現(xiàn)無限地基輻射阻尼效應(yīng),用以驗(yàn)證黏彈性人工邊界地震動(dòng)輸入模型及公式的正確性。

        算例:從三維半無限空間中取邊長(zhǎng)為50 m的立方體,建立整體笛卡爾坐標(biāo)系,模型示意圖如圖2所示,自由場(chǎng)表面中心為坐標(biāo)原點(diǎn),取此點(diǎn)為散射源O。散射源O到人工邊界的距離rb,選取散射源到人工邊界各結(jié)點(diǎn)的直線距離,采用八結(jié)點(diǎn)六面體塊體單元SOLID185離散,網(wǎng)格尺寸為1 m。假定材料為均質(zhì)彈性各向同性材料,彈性模量E=24 MPa,泊松比μ=0.2,剪切模量G=10 MPa,質(zhì)量密度ρ=1 000 kg/m3,剪切波速Cs=100 m/s,縱波波速Cp=163.3 m/s,材料的阻尼系數(shù)取零。在模型底部輸入三條單位脈沖位移波,位移、速度時(shí)程的表達(dá)式如式(14)、式(15)所示,時(shí)間步取0.001 s,總時(shí)長(zhǎng)為2.0 s,頻率為f=4.0 Hz。位移波、速度波波形如圖3所示。取散射源O點(diǎn)(0,0,0)及底部點(diǎn)A(0,0,-50)和中部點(diǎn)B(0,0,-25)為代表點(diǎn)。

        圖2 算例模型有限元網(wǎng)格

        (1)位移函數(shù)

        (14)

        (2)速度函數(shù)

        (15)

        圖3 位移、速度時(shí)間歷程曲線

        圖4給出了代表點(diǎn)的水平、豎向位移時(shí)程曲線。由圖可見,波動(dòng)從立方體底部黏彈性邊界向上傳播,0.25 s時(shí),到達(dá)塊體中部后,在0.5 s后到達(dá)立方體頂部,在頂部位移波幅值放大為原入射位移波的2倍左右,數(shù)值解結(jié)果與解析解十分接近。波動(dòng)在自由表面反射后,向下繼續(xù)傳播,且傳播回底邊界的波動(dòng)不再向上反射,黏彈性邊界具有較好的吸能效果。因此經(jīng)過自由場(chǎng)模型的驗(yàn)證,三維黏彈性邊界地震動(dòng)輸入方法是合理的。

        圖4 代表點(diǎn)位移時(shí)程曲線

        4 SSI效應(yīng)對(duì)冷卻塔地震響應(yīng)的影響分析

        4.1 工程概況

        本文研究對(duì)象為國(guó)內(nèi)某設(shè)計(jì)容量為2×350 MW熱電廠一座大型雙曲冷卻塔。空冷系統(tǒng)采用表面式凝汽間接空冷系統(tǒng),采用兩機(jī)一塔方案。冷卻塔整體高度為165.0 m,進(jìn)風(fēng)口高度27.5 m,喉部高度為132.0 m;±0.000位置直徑為155.8 m,塔頂直徑為108.20 m。塔筒由進(jìn)風(fēng)口至塔頂共劃分107節(jié)模板,殼體厚度通過塔筒殼體整體和局部穩(wěn)定驗(yàn)算確定,最大和最小厚度分別為1.70 m和0.29 m。斜支柱采用40對(duì)矩形截面的X型支柱,截面尺寸為1.0 m(環(huán)向)×1.7 m(徑向)。環(huán)基為矩形截面,尺寸為9.0 m(徑向)×2.0 m(豎向),中面半徑為79.442 m,環(huán)基頂面標(biāo)高為-4.6 m。塔筒、斜支柱和支墩采用C40混凝土,環(huán)基采用C30混凝土,鋼筋采用HRB400級(jí)鋼筋。該冷卻塔結(jié)構(gòu)特征尺寸示意圖如圖5所示。

        圖5 冷卻塔結(jié)構(gòu)特征尺寸(m)

        4.2 冷卻塔數(shù)值計(jì)算模型的建立

        利用ANSYS程序建立該冷卻塔的有限元模型,水平向分別為X、Y向,豎向?yàn)閆向,以±0.000處冷卻塔軸對(duì)稱中心為坐標(biāo)原點(diǎn),建立整體笛卡爾坐標(biāo)系。

        模型1:剛性地基模型

        冷卻塔塔筒采用三維空間殼體單元SHELL63離散,X支柱、環(huán)形基礎(chǔ)、支墩采用三維空間梁?jiǎn)卧狟EAM188離散,環(huán)基上表面的覆土及設(shè)備采用附加質(zhì)量單元MASS21離散,共劃分28 320個(gè)結(jié)點(diǎn),26 880個(gè)單元。建立的模型如圖6(a)所示。

        模型2:無質(zhì)量地基模型

        冷卻塔部分與剛性地基模型一致,近域地基土范圍:300 m(水平向)×300 m(水平向)×150 m(豎向),采用三維空間八結(jié)點(diǎn)塊體單元SOLID185離散,共劃分42 910個(gè)結(jié)點(diǎn),41 080個(gè)單元。

        模型3:黏彈性人工邊界模型

        通過APDL語言編制程序,采用三維彈簧單元COMBIN14建立三維黏彈性人工邊界,共劃分49 063個(gè)結(jié)點(diǎn),47 233個(gè)單元。地基—冷卻塔結(jié)構(gòu)體系的有限元模型,如圖6(b)所示。

        地基采用彈性各向同性體,彈性模量E取48 MPa,泊松比μ取0.2,質(zhì)量密度ρ取2 000 kg/m3。

        4.3 結(jié)構(gòu)體系動(dòng)力特性分析

        大型冷卻塔屬于典型的高聳薄壁旋轉(zhuǎn)空間軸對(duì)稱殼體結(jié)構(gòu),其動(dòng)力特性與房屋建筑、大壩、橋梁等常規(guī)結(jié)構(gòu)明顯不同?!稑?gòu)筑物抗震設(shè)計(jì)規(guī)范》:GB50191—2012規(guī)定“當(dāng)按冷卻塔專用有限元程序計(jì)算時(shí),建議每階諧波宜取不少于5個(gè)振型;按通用程序計(jì)算時(shí),建議宜取不少于300個(gè)振型?!北疚奶崛傂缘鼗P秃蜔o質(zhì)量地基模型前500階頻率分布,如圖7所示。由圖7可見,冷卻塔頻率密集前500階頻率未超過10 Hz。由于結(jié)構(gòu)的軸對(duì)稱性,其雙向諧波耦合振型為冷卻塔主振型,數(shù)量超過90%,側(cè)彎傾覆振型、豎向伸縮振型、豎向扭轉(zhuǎn)振型等結(jié)構(gòu)整體振型屬于結(jié)構(gòu)高階振型,且分布較為分散。由于側(cè)彎傾覆振型、豎向伸縮振型對(duì)冷卻塔水平、豎向地震響應(yīng)存在很大貢獻(xiàn),為滿足“計(jì)算振型數(shù)應(yīng)使振型參與質(zhì)量不小于總質(zhì)量的90%”的要求,保證計(jì)算精度,振型分解反應(yīng)譜法計(jì)算需取足夠多振型。

        表1、表2分別列出了兩模型前10階及整體振型情況下的動(dòng)力特性。對(duì)于剛性地基模型,前22階都是局部振型,結(jié)構(gòu)的基頻為環(huán)向4個(gè)諧波,子午向2個(gè)豎向諧波,第23階為水平整體傾覆振型,45階為塔筒整體繞Z軸扭轉(zhuǎn)振型,205階為整體豎向振型。而對(duì)于無質(zhì)量地基模型,前6階為局部振型,結(jié)構(gòu)基頻為環(huán)向3個(gè)諧波,子午向2個(gè)豎向諧波,第7階為水平整體傾覆振型,18階為塔筒整體繞Z軸扭轉(zhuǎn)振型,15階為整體豎向振型。

        (a) 冷卻塔模型

        (b) 地基-冷卻塔有限元模型

        圖7 冷卻塔前500階自振頻率分布Fig.7 First 500 natural frequency distribution for the cooling tower

        表1 結(jié)構(gòu)前10階動(dòng)力特性(模型1:剛性地基模型)

        表2 結(jié)構(gòu)前10階動(dòng)力特性(模型2:無質(zhì)量地基模型)

        從以上剛性地基模型和無質(zhì)量地基模型的模態(tài)分析的結(jié)果來看,①冷卻塔結(jié)構(gòu)體系自振頻率分布非常密集,前500階頻率均未超過10 Hz;②由于冷卻塔結(jié)構(gòu)具有軸對(duì)稱性,其低階局部振型成對(duì)出現(xiàn),頻率相同,低階模態(tài)中各階振型均為不同數(shù)量環(huán)向諧波與子午向諧波的組合;③考慮地基的彈性作用,結(jié)構(gòu)的各階自振頻率有所降低,整體振型可以較早出現(xiàn)。

        4.4 彈性地震響應(yīng)分析

        4.4.1 地震動(dòng)輸入

        在開展水工混凝土大壩等建于高山峽谷中的大型水工結(jié)構(gòu)的地震反應(yīng)分析時(shí),工程界對(duì)于選取自由場(chǎng)表面具有不同的看法;而冷卻塔結(jié)構(gòu)的建設(shè)場(chǎng)地平坦,選定自由場(chǎng)平面的位置一般不存在爭(zhēng)議。

        地震動(dòng)輸入采用El Centro南北向記錄,如圖8(a)所示,步長(zhǎng)為0.02 s,共計(jì)30 s,按8度多遇地震進(jìn)行彈性時(shí)程分析,將加速度幅值調(diào)整為70 cm/s2,圖8(b)、圖8(c)分別給出了利用作者自編程序所消除基線漂移后速度、位移時(shí)程曲線。無質(zhì)量地基模型的地震動(dòng)輸入,直接按動(dòng)力學(xué)基本方程式(16)采用加速度記錄輸入;而黏彈性人工邊界模型采用前文提到的方法進(jìn)行輸入。根據(jù)國(guó)內(nèi)外學(xué)者實(shí)測(cè)的結(jié)果[20],阻尼ζ比取0.025,計(jì)算RAYLEIGH阻尼的質(zhì)量和剛度阻尼系數(shù)α和β,感興趣的頻率段按結(jié)構(gòu)基頻和整體傾覆振型一階頻率選取。α和β計(jì)算公式按式采用式(17)、式(18)進(jìn)行計(jì)算。

        圖8 El Centro記錄加速度、速度、位移時(shí)程曲線

        (16)

        [C]=α[M]+β[K]

        (17)

        (18)

        4.4.2 單元坐標(biāo)系及內(nèi)力符號(hào)規(guī)定

        對(duì)于具有軸對(duì)稱性的冷卻塔結(jié)構(gòu),配筋設(shè)計(jì)時(shí)一般僅需關(guān)注各荷載組合作用下,沿結(jié)構(gòu)高度處的內(nèi)力的環(huán)向的包絡(luò)值,而對(duì)內(nèi)力環(huán)向分布特征的關(guān)注程度則較少。為了便于描述塔筒殼體、X支柱的內(nèi)力最大值沿環(huán)向分布特點(diǎn),以水平向地震激勵(lì)正方向與塔筒環(huán)向交點(diǎn)分別為θ=0°和θ=180°,其中,θ以俯視圖逆時(shí)針方向?yàn)檎?,如圖9(a)所示。

        塔筒殼體內(nèi)力主要考察沿子午向、環(huán)向的單寬軸力(Tx,Ty)、單寬彎矩(Mx,My)。對(duì)于X支柱,主要考察軸力剪力(Fx,Fy,Fz)和扭矩彎矩(Mx,My,Mz)。殼單元,梁?jiǎn)卧鴺?biāo)系及內(nèi)力正負(fù)號(hào)規(guī)定如圖9(b)所示,圖中所示內(nèi)力符號(hào)均為正值。

        (a) 觀測(cè)點(diǎn)位置示意圖

        (b) 殼單元、梁?jiǎn)卧鴺?biāo)系及內(nèi)力符號(hào)

        4.4.3 塔筒的地震響應(yīng)分析結(jié)果

        (1) 塔筒絕對(duì)加速度最大值沿高度變化規(guī)律

        絕對(duì)加速度大小是結(jié)構(gòu)受地震作用響應(yīng)大小的直觀體現(xiàn),對(duì)于剛性地基模型和無質(zhì)量地基模型,直接按式(16)的動(dòng)力方程求解,其絕對(duì)加速度需考慮地面運(yùn)動(dòng)加速度,按式(19)確定;黏彈性人工邊界模型則直接提取結(jié)構(gòu)上的絕對(duì)加速度。

        (19)

        圖10給出了三個(gè)模型在冷卻塔環(huán)向幾個(gè)特征角度塔筒子午線上各點(diǎn)的絕對(duì)加速度最大值沿其高度的變化規(guī)律。由圖10(a)可見應(yīng)用傳統(tǒng)剛性地基模型:①沿YOZ平面兩側(cè)對(duì)稱的子午線上同一高度處的絕對(duì)加速度最大值基本一致,即0°和180°,45°和135°結(jié)果一致,這與其結(jié)構(gòu)軸對(duì)稱性的因素是分不開的;②塔筒底部的絕對(duì)加速度最大值較大,沿高度先增加,后逐漸變小,塔頂?shù)慕^對(duì)加速度最大值為1.483 m/s2;由圖10(b)可見,對(duì)于無質(zhì)量地基模型而言,塔筒絕對(duì)加速度最大值喉部以下變化幅度不大,喉部以上急劇增加,塔筒頂部最大值達(dá)到1.917m/s2;由圖10(c)可見,對(duì)于考慮土-結(jié)構(gòu)動(dòng)力相互作用的黏彈性人工邊界模型而言,絕對(duì)加速度最大值沿高度變化,先減小而后呈現(xiàn)線性增加,頂部最大為0.840 m/s2。由圖10(d)可見,對(duì)于同一條地震波輸入的時(shí)程分析,體系的絕對(duì)加速度最大值沿其高度呈現(xiàn)完全不同的發(fā)展趨勢(shì),但模型3,考慮土-結(jié)構(gòu)動(dòng)力相互作用的黏彈性人工邊界模型的塔筒的絕對(duì)加速度最大值降低43.4%。

        圖10 塔筒絕對(duì)加速度最大值沿其高度變化

        (2) 塔筒的內(nèi)力響應(yīng)分析

        為了研究塔筒內(nèi)力沿環(huán)向變化規(guī)律,選取了最底層(即1號(hào)模板)分析其內(nèi)力的環(huán)向變化規(guī)律,如圖11(c)和圖11(d)所示;圖11(e)給出了黏彈性人工邊界模型與剛性地基模型相應(yīng)最大值的比值。黏彈性人工邊界模型的子午向彎矩Mx和Tx僅為剛性地基模型結(jié)果的一半;而環(huán)向彎矩My和Ty變化幅度較大,最大分別為剛性地基結(jié)果的3.05倍和1.59倍。

        (a) 塔筒子午向內(nèi)力最大值分布 (剛性地基)(b) 塔筒子午向內(nèi)力最大值分布 (黏彈性邊界)(c) 1號(hào)模板環(huán)向內(nèi)力分布 (剛性地基)

        (d) 1號(hào)模板環(huán)向內(nèi)力分布(黏彈性邊界)(e) 1號(hào)模板環(huán)向各內(nèi)力對(duì)比(黏彈性邊界/剛性地基)圖11 塔筒內(nèi)力環(huán)向分布特征

        Fig.11 Distribution characteristic of internal force for the cooling tower in the ring direction

        表3給出了兩種模型塔筒內(nèi)力幅值及變化率。綜上所述,黏彈性人工邊界模型的塔筒內(nèi)力相對(duì)于剛性地基模型,子午向內(nèi)力(Mx,Tx)有所降低,而環(huán)向內(nèi)力(My,Ty)則有所增加。這主要是由于剛性地基模型各點(diǎn)地震動(dòng)輸入一致,地震動(dòng)響應(yīng)主要為整體傾覆振型貢獻(xiàn);而黏彈性人工邊界模型考慮了波動(dòng)過程,因此各輸入點(diǎn)引起的運(yùn)動(dòng)并不同步,提高了其環(huán)向內(nèi)力。

        表3 塔筒內(nèi)力幅值

        注:變化率η=(剛性模型-黏性模型)/黏性模型×100%,下同。

        (3) X支柱的地震響應(yīng)分析

        大型冷卻塔的X支柱是抗震薄弱部位,相比于塔筒,其內(nèi)力分布更為復(fù)雜。圖12分別給出了兩種模型情況下,X支柱軸力Fx、剪力Fy以及繞環(huán)向彎矩My、繞徑向彎矩Mz沿塔筒環(huán)向的變化規(guī)律。表4給出了兩種模型X支柱的內(nèi)力幅值及變化率。

        表4 X支柱內(nèi)力幅值

        綜上所述,無限地基輻射阻尼可大幅降低冷卻塔X支柱的內(nèi)力幅值,軸力降低達(dá)20.5%、剪力36.8%,柱底彎矩約50%。

        (a) X支柱軸力、剪力(剛性地基)(b) X支柱彎矩(剛性地基)

        (c) X支柱軸力、剪力(黏彈性邊界)(d) X支柱彎矩(黏彈性邊界)圖12 X支柱內(nèi)力環(huán)向分布特征

        Fig.12 Distribution characteristic of internal force for X columns in the ring direction

        5 結(jié) 論

        本文探討土-結(jié)構(gòu)動(dòng)力相互作用對(duì)于冷卻塔結(jié)構(gòu)體系動(dòng)力響應(yīng)的影響規(guī)律。主要結(jié)論如下:

        (1)以有限元軟件為平臺(tái),進(jìn)行二次開發(fā),驗(yàn)證了黏彈性人工邊界地震動(dòng)輸入模型及輸入公式的正確性,結(jié)果表明,黏彈性人工邊界具有較好的吸能效果。

        (2)通過對(duì)冷卻塔計(jì)算模型進(jìn)行模態(tài)分析可知,冷卻塔自振頻率分布十分密集,絕大多數(shù)振型為環(huán)向諧波和子午向諧波組合的局部振型;考慮地基彈性后,體系的自振頻率略有降低,整體傾覆振型較早出現(xiàn)。

        (3)從塔筒的絕對(duì)加速度最大值分布規(guī)律可知,輻射阻尼效應(yīng)顯著降低結(jié)構(gòu)體系的動(dòng)力響應(yīng),與剛性地基模型相比,塔筒絕對(duì)加速度最大值降低43.4%。

        (4)考慮地基土-冷卻塔動(dòng)力相互作用,計(jì)入無限地基輻射阻尼效應(yīng)對(duì)冷卻塔塔筒及X支柱內(nèi)力影響較大,塔筒子午向內(nèi)力有所降低,而環(huán)向內(nèi)力則顯著提高。對(duì)X支柱的內(nèi)力最大值則降低20%~50%。顯然,考慮土-結(jié)構(gòu)動(dòng)力相互作用效應(yīng)可使冷卻塔結(jié)構(gòu)設(shè)計(jì)更加經(jīng)濟(jì)合理,避免浪費(fèi);然而不考慮土-結(jié)構(gòu)動(dòng)力相互作用的計(jì)算模型的地震響應(yīng)并非都是保守的,也會(huì)低估塔筒的子午向內(nèi)力值。

        [ 1 ] 柯世堂,陳少林,葛耀君,等.超大型冷卻塔隨機(jī)地震響應(yīng)分析[J].地震工程與工程振動(dòng),2012,32(6):159-165. KE Shitang, CHEN Shaolin, GE Yaojun, et al. Earthquake stochastic response analysis of super large cooling towers[J].Journal of earthquake Engineering and Engineering Vibration, 2012, 32(6): 159-165.

        [ 2 ] 構(gòu)筑物抗震設(shè)計(jì)規(guī)范:GB 50191—2012[S].北京:中國(guó)計(jì)劃出版社, 2012.

        [ 3 ] GUPTA A K, SCHNOBRICH W C. Seismic analysis and design of hyperbolic cooling towers[J].Nuclear Engineering and design, 1976, 36(2): 251-260.

        [ 4 ] WOLF J P. Seismic analysis of cooling towers[J].Engineering Structures, 1986, 8(3): 191-198.

        [ 5 ] NASIR A M, THAMBIRATNAM D P, BUTLER D, et al. Dynamics of axisymmetric hyperbolic shell structures[J].Thin-Walled Structures, 2002, 40(7/8): 665-690.

        [ 6 ] SABOURI-GHOMI S, NIK F A, ROUFEGARINEJAD A, et al. Numerical study of the nonlinear dynamic behavior of reinforced concrete cooling towers under earthquake excitation[J].Advances in Structural Engineering, 2006, 9(3): 433-442.

        [ 7 ] 葉浩,李華峰,黃志龍.基于混凝土損傷塑性模型的自然通風(fēng)冷卻塔非線性抗震分析[J].電力建設(shè),2014,35(9):88-92. YE Hao, LI Huafeng, HUANG Zhilong. Nonlinear seismic response of natural draft cooling tower based on concrete damaged plasticity model[J].Electric Power Construction, 2014, 35(9): 88-92.

        [ 8 ] HORR A M, SAFI M. Full dynamic analysis of large concrete cooling towers: soil structure interaction[J].International Journal of Space Structures, 2002, 17 (4): 301-312.

        [ 9 ] 高標(biāo),盧紅前.SSI效應(yīng)對(duì)大型雙曲線冷卻塔結(jié)構(gòu)抗震性能的影響[J].武漢大學(xué)學(xué)報(bào)(工學(xué)版),2009,42(增刊1):427-431. GAO Biao, LU Hongqian. Influence of SSI effect on aseismic performance of hyperbolic cooling tower[J].Engineering Journal of Wuhan University Engineering Science, 2009, 42(Sup1): 427-431.

        [10] 李輝,張俊發(fā).超大型冷卻塔考慮土-結(jié)構(gòu)相互作用的三維有限元抗震分析[J].長(zhǎng)春工程學(xué)院學(xué)報(bào)(自然科學(xué)版), 2012,13(4):19-21. LI Hui, ZHANG Junfa. The Three-dimensional finite element seismic analysis to super large cooling tower considering soil-structure interaction[J].Journal of Changchun Institute Technology(Natural Science Edition), 2012, 13(4): 19-21.

        [11] 房營(yíng)光.巖土介質(zhì)與結(jié)構(gòu)動(dòng)力相互作用理論及其應(yīng)用[M].北京:科學(xué)出版社,2005.

        [12] 薛素鐸,劉毅,李雄彥.土-結(jié)構(gòu)動(dòng)力相互作用研究若干問題綜述[J].世界地震工程,2013,29(2):1-9. XUE Suduo, LIU Yi, LI Xiongyan. Review of some problems about research on soil-structure dynamic interaction[J].World Information on Earthquake Engineering,2013,29(2): 1-9.

        [13] VELETSOS A S, WEI Y T. Lateral and Rocking vibrations of footings[J].Journal of the Soil Mechanism and Foundations Division, ASCE, 1971, 97(1):1227-1248.

        [14] LUCO J E, WESTMAN R A. Dynamic response of circular footings[J].Journal of the Engineering Mechanism Division, ASCE, 1971, 97(1):1381-1395.

        [15] 賀廣零.考慮土-結(jié)構(gòu)相互作用的風(fēng)力發(fā)電高塔系統(tǒng)地震動(dòng)力響應(yīng)分析[J].機(jī)械工程學(xué)報(bào),2009,45(7):87-94. HE Guangling. Seismic response analysis of wind turbine tower systems considering soil-structure interaction[J].Journal of Mechanical Engineering,2009,45(7): 87-94.

        [16] 柯世堂,王同光,曹九發(fā),等.考慮土-結(jié)相互作用大型風(fēng)力發(fā)電結(jié)構(gòu)風(fēng)致響應(yīng)分析[J].土木工程學(xué)報(bào),2015,48(2):18-25. KE Shitang, WANG Tongguang, CAO Jiufa, et al. Analysis on wind-induced responses of large wind power structures considering soil-structure interaction[J].China Civil Engineering Journal, 2015, 48(2): 18-25.

        [17] 曹青,張豪.考慮土-結(jié)構(gòu)相互作用的風(fēng)力發(fā)電機(jī)塔架地震響應(yīng)分析[J].西北地震學(xué)報(bào),2011,33(1):62-66. CAO Qing, ZHANG Hao. Seismic response analysis of wind turbine tower with soil-structure interaction[J].Northwestern Seismological Journal, 2011, 33(1): 62-66.

        [18] 何建濤,馬懷發(fā),張伯艷,等.黏彈性人工邊界地震動(dòng)輸入方法及實(shí)現(xiàn)[J].水利學(xué)報(bào),2010,41(8):960-969. HE Jiantao, MA Huaifa, ZHANG Boyan, et al. Method and realization of seismic motion input of viscous-spring boundary[J]. Journal of Hydraulic Engineering, 2010, 41(8): 960-969.

        [19] 劉晶波,王振宇,杜修力,等.波動(dòng)問題中的三維時(shí)域黏彈性人工邊界[J].工程力學(xué),2005,22(6):46-51. LIU Jingbo, WANG Zhenyu, DU Xiuli, et al. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems[J].Engineering Mechanics, 2005, 22(6): 46-51.

        [20] RADWANSKA M, WASZCZYSZYN Z.Buckling analysis of a cooling tower shell with measured and theoretically-modeled imperfections[J].Thin-Walled Structures, 1995, 23(1): 107-121.

        Seismic response analysis for a cooling tower structure considering soil-structure dynamic interaction

        TAO Lei1, ZHANG Junfa2, CHEN Houqun3

        (1. College of Water Resources and Hydroelectric Engineering, Xi’an University of Technology, Xi’an 710048, China; 2. College of Civil Engineering and Architecture, Xi’an University of Technology, Xi’an 710048, China; 3. Engineering Earthquake Research Center, China Institute of Water Resources and Hydropower Research, Beijing 100048, China)

        A large-scale high-rise reinforced concrete cooling tower belongs to a typical wind-sensitive structure. In recent years, its dynamic response study under extreme external actions, such as, earthquake, jet-crash, blasting et. al receives extensive attention in engineering field. The study achienements about the effects of soil-structure interaction (SSI) on seismic responses of dams and bridges are relatively abundant, however, they are seldom involved in cooling towers. Here, based on the elastic wave motion theory combined with the finite element method, a 3D visco-elastic artificial boundary model and its formulas were built and deduced considering soil-structure dynamic interaction, and the correctness of the earthquake input method was verified with the half space free field model. A large cooling tower in a certain coal-fired power plant, one of the domestic engineering projects was taken as a study background, based on the FE software ANSYS, the cooling tower’s rigid foundation model, its massless foundation model and its viscoelastic artificial boundary model were modeled, respectively. The modal analysis and the elastic time history analysis were performed in order to study dynamic characteristics and changes of internal force of these models. The influence laws of soil-structure dynamic interaction on different models were explored. The research results showed that the cooling tower structure system natural vibration frequency distribution considering rigid foundation is very dense, and the majority of vibration modes are local coupling vibration modes in ring and meridian directions; the natural vibration frequencies are slightly reduced considering elastic foundation, and the whole body vibration modes appear earlier; with the time history analysis, using the visco elastic artificial boundary model, considering the infinite foundation radiation damping effect, compared with the rigid foundation model, the maximum absolute acceleration of the tower tube reduces 43.4%, the moment and axial force in meridian direction reduce 50%, and they in ring direction increase significantly; internal force amplitude values of X column reduce about 20%~50%; therefore, when performing the seismic response analysis of the cooling tower, the influences of soil-structure interaction can not be neglected.

        cooling tower structure; seismic response analysis; soil-structure dynamic interaction; visco-elastic artificial boundary; free field analysis; earthquake input

        國(guó)家自然科學(xué)基金(51179162)

        2016-02-29 修改稿收到日期:2016-07-05

        陶磊 男,博士生,1982年生

        張俊發(fā) 男,博士,教授,1961年生 E-mail:zhangjf-4314@163.com

        TU279.7+41;TU33+2;P315.92

        猜你喜歡
        結(jié)構(gòu)模型
        一半模型
        《形而上學(xué)》△卷的結(jié)構(gòu)和位置
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        論結(jié)構(gòu)
        中華詩詞(2019年7期)2019-11-25 01:43:04
        新型平衡塊結(jié)構(gòu)的應(yīng)用
        模具制造(2019年3期)2019-06-06 02:10:54
        論《日出》的結(jié)構(gòu)
        3D打印中的模型分割與打包
        FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
        創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長(zhǎng)
        国产精品麻花传媒二三区别 | 免费中文熟妇在线影片| 伊人影院在线观看不卡| 蜜桃视频在线免费视频| 日韩精品中文字幕人妻中出| 成人久久精品人妻一区二区三区| 免费国产在线精品一区二区三区免| 中文字幕亚洲精品一区二区三区| 亚洲啪av永久无码精品放毛片| 99久久久无码国产精品9| 8av国产精品爽爽ⅴa在线观看| 一本色道久久综合亚州精品| 亚洲熟女熟妇另类中文| 国产aⅴ无码专区亚洲av| 亚洲色欲色欲www在线播放| 国产av综合一区二区三区最新| 美女福利一区二区三区在线观看| 国产成人国产三级国产精品| а√天堂资源官网在线资源| 人体内射精一区二区三区| 亚洲精品二区在线观看| 国产成人综合精品一区二区| 4hu四虎永久免费地址ww416| 亚洲国产综合精品 在线 一区 | 熟女中文字幕一区二区三区| 国产成人精品综合在线观看| 日本午夜国产精彩| 亚洲人成在线播放a偷伦| 99久久精品一区二区国产| 亚洲精品v欧洲精品v日韩精品 | 天堂69亚洲精品中文字幕| 日韩av在线手机免费观看| 精品亚洲一区二区三区四区五区| 色八区人妻在线视频免费| 思思久久99er热只有频精品66| 国产码欧美日韩高清综合一区| 手机在线免费观看av不卡网站| 狠狠色噜噜狠狠狠777米奇小说| 香蕉视频一级| 亚洲精品中文字幕尤物综合| 国产av剧情一区二区三区|