基于脈動風(fēng)作用下的大跨度橋梁拉索風(fēng)場數(shù)值模擬

孫亭亭,楊吉新,黎建華,周興宇

?

基于脈動風(fēng)作用下的大跨度橋梁拉索風(fēng)場數(shù)值模擬

孫亭亭1,2,楊吉新2*,黎建華2,周興宇3

1. 安慶職業(yè)技術(shù)學(xué)院 建筑工程系, 安徽 安慶 246003 2. 武漢理工大學(xué) 交通學(xué)院, 湖北武漢 430063 3. 中國市政工程中南設(shè)計研究總院有限公司, 湖北 武漢 430010

采用Kaimal脈動風(fēng)速譜，利用MATLAB模擬脈動風(fēng)的瞬時風(fēng)速時程并導(dǎo)入計算流體力學(xué)軟件中，選擇RNG-模型，考慮在高雷諾數(shù)的作用下，對大跨度橋梁的拉索風(fēng)場數(shù)值模擬。采用結(jié)構(gòu)性網(wǎng)格建立流場模型，分別選取30 m和80 m高度處分析流場和風(fēng)壓力系數(shù)的變化；選取50 m高度處分析速度矢量；選取80 m高度處分析流速和跡線變化；選取拉索中心位置、拉索尾端分析拉索周圍壓力分布，并進行數(shù)值模擬。計算結(jié)果表明，對比文獻中風(fēng)洞試驗的結(jié)果，流場中流體的分布特征符合類似研究的結(jié)論，可以得到合理精度的數(shù)值解；在拉索周圍上、下表面形成渦街，渦流沿中軸線運動并呈現(xiàn)對稱的形狀；拉索表面壓力分布和文獻中風(fēng)洞試驗結(jié)果較吻合，并且符合工程實際。

Kaimal脈動風(fēng)速譜; 三維風(fēng)場; RNG-模型; 數(shù)值模擬

在大跨度橋梁中斜拉橋是常用的橋型，斜拉索是斜拉橋中重要的結(jié)構(gòu)組成，斜拉索的風(fēng)荷載在全橋風(fēng)荷載中占得比重比較大，所以在風(fēng)荷載作用下斜拉索周圍流場及其響應(yīng)的問題研究就顯得尤為重要[1]。從以往的研究結(jié)果來看，大跨度橋梁風(fēng)荷載作用主要采用風(fēng)洞試驗的方法來檢驗斜拉索和橋梁的抗風(fēng)性能[2,3]。比如在蘇通長江公路大橋的建設(shè)中，全橋有272根斜拉索，其中最長的斜拉索長度是577 m[4]。研究表明在橫橋向風(fēng)的作用下，斜拉索產(chǎn)生的位移和內(nèi)力對于主梁位移來說占整個風(fēng)荷載的60%~70%[5,6]。目前，關(guān)于斜拉索風(fēng)荷載的計算依據(jù)主要有歐洲規(guī)范、日本規(guī)范和我國的JTG/T D60-01-2004《公路橋梁抗風(fēng)設(shè)計規(guī)范》，同時這些規(guī)范也被作為風(fēng)洞試驗的依據(jù)[7]，本文在計算中的參數(shù)主要采用我國的JTG/T D60-01-2004《公路橋梁抗風(fēng)設(shè)計規(guī)范》[8]。風(fēng)洞試驗對于研究斜拉索以及斜拉橋抗風(fēng)具有重要意義，然而風(fēng)洞試驗造價昂貴，時間周期較長[9,10]。本文采用計算流體力學(xué)CFD數(shù)值風(fēng)洞的方法來對斜拉索進行研究并結(jié)合相關(guān)文獻的風(fēng)洞試驗結(jié)果進行分析[11]。

現(xiàn)在流體力學(xué)研究方法主要包括理論分析、數(shù)值計算、試驗研究和現(xiàn)場實測四個方面[12]，這些研究方法針對不同的角度相互補充。理論分析研究主要是從各參數(shù)影響出發(fā)，為數(shù)值計算和實驗研究提供有效的理論指導(dǎo)[13]；試驗研究是驗證客觀科學(xué)規(guī)律的一種有效手段，用來驗證理論分析和數(shù)值計算的正確性；現(xiàn)場實測是有效收集現(xiàn)場數(shù)據(jù)的一種方法，能夠反映自然界中的真實情況，但是由于各種環(huán)境因素的影響，測量往往有困難，甚至是無法實施[14]；數(shù)值計算主要通過計算流體力學(xué)的理論來模擬實際的情況，包括實際環(huán)境中的復(fù)雜情況，從而補充理論及試驗的空缺，更重要的是提供了廉價的模擬、設(shè)計和優(yōu)化的工具，以及提供了分析三維復(fù)雜流動的工具[13,15]。從上世紀80年代以來，計算流體力學(xué)(Computational ?uid dynamics CFD)發(fā)展異常迅速，結(jié)構(gòu)風(fēng)工程的研究除了繼續(xù)采用風(fēng)洞試驗的傳統(tǒng)方法以外，研究者開始采用計算流體力學(xué)的技術(shù)并使用數(shù)值模擬的方法研究結(jié)構(gòu)在流場中的作用[16,17]。在數(shù)值風(fēng)洞分析中高雷諾數(shù)下的流場變化無論對于理論研究還是工程實踐都具有重要意義[18]，同時也是風(fēng)工程中需要解決的關(guān)鍵問題之一[19,20]。本文利用CFD數(shù)值模擬的方法來研究斜拉橋拉索在高雷諾數(shù)作用下采用RNG-模型拉索表面壓力分布和不同高度處流場的分布。

1 計算理論分析

1.1 脈動風(fēng)速譜

在工程實際應(yīng)用中，通常認為風(fēng)速包含平均風(fēng)和脈動風(fēng)兩部分，風(fēng)荷載的模擬主要是針對脈動風(fēng)而言[21]。脈動風(fēng)有兩個主要概率特性：脈動風(fēng)的風(fēng)速譜和相干函數(shù)[11]。本文在風(fēng)荷載計算中主要采用Kaimal脈動風(fēng)速譜[21]。

在橋梁風(fēng)工程中，氣流相對于橋梁而言是湍流充分發(fā)展的狀態(tài)，而Kaimal脈動風(fēng)速譜考慮了大氣湍流運動中的湍流功率譜隨高度的變化，所以采用Kaimal脈動風(fēng)速譜能夠更好的反應(yīng)出風(fēng)場的真實情況。

1.2 脈動風(fēng)荷載的模擬

模型參數(shù)，設(shè)[()]={()]=[1(),2(),3(),…,u()]，u()為第點處的風(fēng)速序列為時間序列，可以知道，根據(jù)自回歸模型特點，第個變量不僅與本身前一個或幾個時間階段的風(fēng)速u(-D)有關(guān)，而且與其他位置處同時段的風(fēng)速u()和前一個或幾個時段的風(fēng)速u(-D)有關(guān)，此外還受到隨機因素的影響[22]。寫成數(shù)學(xué)表達式如下：

式中[(-D)]=[u(-D)]，是與時刻相關(guān)的個風(fēng)速，是自回歸模型的階數(shù)。

采用Kaimal脈動風(fēng)速譜，對所模擬出的風(fēng)速時程離散值做數(shù)理統(tǒng)計，得到的模擬功率譜與目標功率譜吻合[23]。模型高度共分10個部分，選取第3部分和第8部分作為研究對象，通過MATLAB軟件基于AR模型方法模擬拉索多變量互相關(guān)水平脈動風(fēng)速時程[20]，對此斜拉索進行結(jié)構(gòu)風(fēng)場模擬，提取脈動風(fēng)場湍流強度。將此拉索在豎直方向上分為10段，每10 m一段，斜拉索的10個點編號從下往上依次為1~10。

圖1為3號節(jié)點和8號節(jié)點瞬時風(fēng)速時程的模擬結(jié)果，模擬時間步長為0.1 s，共設(shè)置2048步。通過圖1觀察可以看出：脈動風(fēng)為完全不規(guī)則的隨機脈動，且同一時刻拉索不同位置的風(fēng)速也不盡相同。脈動風(fēng)速在-10 m/s~10 m/s之間浮動，節(jié)點瞬時風(fēng)速在30 m/s~50 m/s，體現(xiàn)風(fēng)速在時間及空間位置上的隨機性。

圖 1 3、8號節(jié)點瞬時風(fēng)速時程

圖 2 節(jié)點脈動風(fēng)速模擬譜與目標譜

圖2表示根據(jù)選取上述節(jié)點模擬結(jié)果，利用快速傅里葉變換(FFT)求得模擬風(fēng)速和Kaimal脈動風(fēng)速譜的對比。這里采用雙對數(shù)坐標軸形式來表示，可見兩者相當(dāng)吻合[22,23]。

2 氣動力和流場特性分析

2.1 相關(guān)參數(shù)分析

本算例以圓柱體作為計算模型，模型本身定義的坐標系為體軸坐標系，根據(jù)風(fēng)向坐標定義的坐標系為風(fēng)軸坐標系。采用RNG-模型進行數(shù)值模擬，RNG-模型來源于嚴格的統(tǒng)計技術(shù)。它和標準-模型很相似，但是在計算中有效改善了計算精度，在湍流漩渦中精度比標準-模型高，為湍流Prandtl數(shù)提供了一個解析公式[24]。這些特點使得RNG-模型比標準-模型在更廣泛的流動中有更高的可信度和精度。所以本文選擇RNG-模型對拉索的風(fēng)場進行分析，使得數(shù)值計算的結(jié)果和實際工程更加符合。

由于缺乏風(fēng)洞試驗條件，為驗證數(shù)值風(fēng)洞模型計算準確性，采取該建模思想對同濟大學(xué)土木實驗室的試驗?zāi)Ｐ瓦M行建模，選取拉索傾角=90°，風(fēng)速為35 m/s的工況，求解此工況下的拉索風(fēng)阻系數(shù)，與物理風(fēng)洞試驗進行對比，對比結(jié)果如表1所示[25]。將計算結(jié)果與試驗結(jié)果對比發(fā)現(xiàn)，采用此計算模型對流場部分進行模擬是可行的[25]。

表1 數(shù)值計算與試驗結(jié)果對比

2.2 控制方程

計算風(fēng)工程中，鈍體繞流問題的控制方程是粘性不可壓Navier-Stokes方程，基于雷諾平均的控制方程可寫為[18-20]：

3 幾何建模與網(wǎng)絡(luò)劃分

拉索尺寸高度100 m直徑D是0.15 m，在建模時采用三維仿真實體建模，模型尺寸3.75 m×3 m×100 m。計算區(qū)域為一個長方體：上游來流為10 D，下游15 D，寬20 D，高度為100 m，其中拉索直徑D為0.15 m，其中拉索的材料密度為8400 kg/m3（考慮防護材料質(zhì)量），彈性模量為195 GPa。

在數(shù)值計算中網(wǎng)格的質(zhì)量直接影響到計算結(jié)果的精確性和計算時間的長短，同時網(wǎng)格的質(zhì)量也是作為工程應(yīng)用的有效工具所面臨的關(guān)鍵問題之一[14]。在計算敏感區(qū)(壁面附近、尾流區(qū)、外形曲率大的表面處)參數(shù)變化梯度大的區(qū)域，對這些區(qū)域需特別注意，選取較密網(wǎng)格。按照本文模型的特點，采用前處理軟件Icem中的O-Block型劃分，并對拉索周圍流場進行加密，流場采用結(jié)構(gòu)化網(wǎng)格，計算分析結(jié)構(gòu)化網(wǎng)格比非結(jié)構(gòu)化網(wǎng)格計算精度較高，共建立1210824個單元，經(jīng)檢查網(wǎng)格質(zhì)量較好。網(wǎng)格劃分如圖3（a，b，c）所示。

圖 3 網(wǎng)格劃分示意圖

a:總體模型，上頂面和下底面 Overall model, top and bottom surface；b:整體網(wǎng)格劃分Overall mesh generation；c:拉索周圍加密的網(wǎng)格區(qū)域 Encrypted mesh area around the cable

4 CFD數(shù)值結(jié)果討論與分析

4.1 風(fēng)壓分析

以風(fēng)壓系數(shù)表示拉索周圍的壓力分布。

為相對壓力（風(fēng)壓—標準大氣壓力），正值為壓力，負值為吸力。

本文計算考慮在高雷諾數(shù)2.4×105下的風(fēng)場的變化，計算采用RNG-模型。分別考慮在拉索30 m和80 m高度處建立切面，研究在各個不同高度處拉索表面風(fēng)壓的變化和流場的壓力分布如圖4和5所示：

圖 4 30 m處拉索周圍壓力分布和30 m高度處流場中壓力分布

Fig.4 Pressure distribution around the cable and the flow field at 30 m

圖 5 80 m處拉索周圍壓力分布和80 m高度處流場中壓力分布

圖4中(a)顯示在30 m高度時的拉索表面風(fēng)壓分布，區(qū)間在（-2.121~2.406），(b)表示30 m高度處流場壓力的分布；圖5中(a)顯示在80 m高度時的拉索表面風(fēng)壓分布，區(qū)間在（-2.745~1.768），(b)表示80 m高度處流場壓力的分布。

對比不同高度處的風(fēng)壓力分布可以看出隨著拉索高度的不斷變化在拉索迎風(fēng)面壓力不斷減小，背風(fēng)面壓力不斷增加（圖6）。

從壓力系數(shù)和壓力云圖中可以看出，拉索周圍壓力場分布的情況，拉索迎風(fēng)面為正壓區(qū)最大值1500 KN，兩側(cè)為負壓區(qū)最大值2000 KN，拉索末端壓力呈現(xiàn)-500 KN~0 KN的變化，這是由于流體繞物體流動產(chǎn)生邊界層分離，形成渦旋，渦旋消耗來流能量，因而尾流中壓力降低所引起的。

圖 6 不同高度處壓力系數(shù)分布

4.2 流場分析

對于圓柱體，由于空氣的粘性效應(yīng)，在表面會形成兩個區(qū)域：一個是離表面很近的邊界層區(qū)，另一個是主流區(qū)，其空氣粘性摩擦力很小，基本可以忽略不計。當(dāng)主流區(qū)內(nèi)存在逆壓梯度使氣流加速或減速時，這種壓力差也會影響邊界層區(qū)流體的流動。

圖 7 Z=50 m時速度矢量圖

圖 8 Z=80時流速分布圖和跡線圖

選取50 m高度處的速度矢量圖進行分析，如圖7所示，由速度矢量圖可見，在拉索周圍上、下表面脫落的渦沿尾流區(qū)的中軸線運動，并且在中軸線附近存在著一個上、下交界區(qū)，并保持著一定的慣性持續(xù)向流場下流運動，形成一個完整的下流渦街。

選取80 m處的流場切片顯示索端部流場，如圖8所示。在邊界層中的流體將會受到更為強烈的阻滯，當(dāng)能量消耗完后會被迫折回發(fā)生倒流，而主流能量相對較大會繼續(xù)前進。于是越來越多的流體在邊界區(qū)和主流區(qū)堆積，形成一個或多個旋渦。由流線圖可見，因為在計算中設(shè)置的計算域關(guān)于x軸對稱，所以在渦流的形狀、拉索表面的壓力分布、拉索端部和拉索尾部的壓力分布也呈現(xiàn)對稱的形狀。

依照建模坐標選取=0、=0.1、=0.025和=0.068處的位置作為切面，包含了拉索的中心位置和尾端，通過數(shù)值計算分析各表面的壓力分布，如圖9所示。從計算結(jié)果中可以看出各表面的壓力等值線與文獻[26]風(fēng)洞試驗結(jié)果比較吻合，F(xiàn)luent中RNG-模型對鈍體結(jié)構(gòu)迎風(fēng)面的模擬結(jié)果較為理想。

圖9 拉索周圍壓力分布

5 結(jié)論

本文主要通過數(shù)值模擬的方法分析大跨度橋梁拉索在脈動風(fēng)作用下三維流場，得到以下結(jié)論：

（1）采用Kaimal脈動風(fēng)速譜進行計算，考慮在平均風(fēng)速的作用下，利用MATLAB模擬脈動風(fēng)在參考點的瞬時風(fēng)速時程，并采用快速傅里葉變換(FFT)模擬風(fēng)速和Kaimal脈動風(fēng)速譜的對比，使得模擬結(jié)果更符合工程實際；

（2）為了提高計算精度，對構(gòu)建的三維流場進行結(jié)構(gòu)性網(wǎng)格的劃分，并采用RNG-湍流模型進行數(shù)值模擬，計算結(jié)果和風(fēng)洞試驗的結(jié)果吻合并能得到合理精度的數(shù)值解，說明選擇RNG-湍流模型能夠反應(yīng)出工程實際的風(fēng)場狀態(tài)；

（3）對拉索進行分段研究，分別計算在拉索30 m和80 m高度處建立切面，比較各個高度的風(fēng)壓力系數(shù)。拉索迎風(fēng)面正壓力最大值1500 KN，兩側(cè)為負壓力最大值2000 KN，拉索末端壓力呈現(xiàn)-500 KN~0 KN的變化；

（4）通過數(shù)值模擬出的速度矢量圖和流線分布圖可以看出，在拉索周圍上、下表面脫落的渦沿尾流區(qū)的中軸線運動。在氣流經(jīng)過拉索壁面時，在壁面摩擦力和壓力梯度的共同作用下，使得拉索壁面的空氣質(zhì)點的運動受阻，質(zhì)點的動量不斷消耗，不能維持原來的運動，一些質(zhì)點無法沿壁面向前運動，出現(xiàn)回流；

（5）選擇拉索的中心位置（關(guān)于軸對稱）進行研究，通過數(shù)值計算分析各表面的壓力分布，從計算結(jié)果中可以看出渦流和壓力呈現(xiàn)對稱的分布，是因為拉索在豎向平面內(nèi)周期性運動，拉索壁面上下分離的剪切層之間的相互作用形成渦街；

（6）對拉索的尾端以及下游區(qū)域進行研究，從計算結(jié)果中可以看出隨著拉索高度的不斷增加，拉索表面及周圍的壓力梯度不斷增大，旋渦的尺寸也發(fā)生變化，邊界層出現(xiàn)層流分離，其尾流逐漸變成紊流渦街。

[1] 張志田,張顯雄,陳政清.橋梁氣動力CFD模擬中湍流模型的應(yīng)用現(xiàn)狀[J].工程力學(xué),2016,33(6):1-8

[2] 陶天友,王浩.基于Hermite插值的簡化風(fēng)場模擬[J].工程力學(xué),2017,34(3):182-188

[3] 曾要爭,喻梅.大跨度獨塔連續(xù)剛構(gòu)斜拉橋抗風(fēng)性能研究[J].公路交通技術(shù),2015(6):55-59

[4] 張掙鑫.斜拉橋拉索參數(shù)振動的半主動控制及MR阻尼器優(yōu)化布置研究[D].長沙:中南大學(xué),2013

[5] Zheng Z, Lei J. Application of the-Re θ, Transition Model to Simulations of the Flow Past a Circular Cylinder[J]. Flow, Turbulence and Combustion, 2016,97(2):401-426

[6] 呂東陽.大跨度屋蓋結(jié)構(gòu)表面風(fēng)壓的數(shù)值模擬研究[D].西安:長安大學(xué),2015

[7] 張建,李天,楊慶山.雷諾數(shù)對大跨屋蓋結(jié)構(gòu)表面風(fēng)壓的影響[J].哈爾濱工業(yè)大學(xué)學(xué)報,2016,48(6):38-42

[8] 何艾狄.基于模態(tài)參數(shù)識別的高聳結(jié)構(gòu)五點等效風(fēng)荷載識別技術(shù)[D].武漢:武漢理工大學(xué),2013

[9] Sarkar S, Ganguly S, Biswas G,. Effect of cylinder rotation during mixed convective flow of nanofluids past a circular cylinder[J]. Computers & Fluids, 2016,127(C):47-64

[10] 方治華,李晨.斜拉索上瞬時風(fēng)速與風(fēng)荷載的數(shù)值模擬[J].振動與沖擊,2010,29(7):210-212

[11] 徐豐,詹昊,梁琛.基于自回歸模型的橋梁脈動風(fēng)場模擬[J].武漢工程大學(xué)學(xué)報,2015,37(3):25-28

[12] 周彬彬,蔡建國,馮健.基于ARMA模型和時空Kriging插值聯(lián)合模擬大跨結(jié)構(gòu)的脈動風(fēng)速時程[J].振動與沖擊,2014,33(3):29-34

[13] Zhang H, Yang JM, Xiao LF,. Large-eddy simulation of the flow past both finite and infinite circular cylinders at=3900[J]. Journal of Hydrodynamics Series B, 2015,27(2):195-203

[14] 汪權(quán),王建國,張鳴祥.高層建筑結(jié)構(gòu)隨機風(fēng)速場的數(shù)值模擬及風(fēng)振控制[J].建筑科學(xué)與工程學(xué)報,2011,28(2):32-37

[15] 賴馬樹金.大跨度懸索橋分離式雙箱梁渦激振動研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2013

[16] Zeng C, Li CW. A hybrid RANS-LES model for combining flows in open-channel T-junctions[C]. Shanghai: 9thInternational Conference on Hydrodynamics, 2010:154-159

[17] 孫芳錦,梁爽.考慮流固耦合作用的大跨度橋梁風(fēng)振響應(yīng)研究[J].公路交通科技,2015,32(7):84-91

[18] 季小勇,陳文禮.斜拉橋斜拉索二維渦激振動的數(shù)值模擬研究[J].公路工程,2012,37(3):6-11

[19] Zheng Z, Lei J, Wu X. Numerical Simulation of the Negative Magnus Effect of a Two-Dimensional Spinning Circular Cylinder[J]. Flow Turbulence & Combustion, 2016,98(1):1-22

[20] 徐闖,孫延國.重慶永川長江大橋渦激振動性能風(fēng)洞試驗研究[J].世界橋梁,2015(5):25-29

[21] 李薇.斜拉索風(fēng)雨激振試驗與數(shù)值模擬研究[D].西安:長安大學(xué),2012

[22] 劉克同,湯愛平,曹鵬.橋梁氣動導(dǎo)數(shù)的格子Boltzmann大渦模擬仿真[J].工程力學(xué),2015,32(5):111-119

[23] 喬永亮,桂洪斌,劉祥鑫.三維圓柱繞流數(shù)值模擬湍流方法的選擇[J].水利水運工程學(xué)報,2016(3):119-125

[24] 譚艷,許棟,白玉川,等.關(guān)于圓柱繞流水動力特性的大渦數(shù)值模擬研究[J].港工技術(shù),2013(3):5-8

[25] 林志興,楊立波,李文勃,等.斜拉索順橋向風(fēng)阻系數(shù)的試驗研究及數(shù)值模擬[C]//2004全國結(jié)構(gòu)風(fēng)工程實驗技術(shù)研 討會論文集.長沙:中國土木工程學(xué)會,2004:51-56

[26] 楊偉,顧明.高層建筑三維定常風(fēng)場數(shù)值模擬[J].同濟大學(xué)學(xué)報:自然科學(xué)版,2003,31(6):647-651

Numerical Simulation of Wind Field on Large Span Bridge Cable Based on Pulsating Winds

SUN Ting-ting1,2, YANG Ji-xin2*, LI Jian-hua2, ZHOU Xing-yu3

1.246003,2.430063,3.430010,

MATLAB is used to simulate the instantaneous wind speed time history of pulsating wind and import it into computational fluid dynamics software by using Kaimal fluctuating wind speed spectrum, the RNG-model is chosen to simulate the wind field of long span bridges under the action of high Reynolds number. The flow field model was established by using structural grid, and the changes of flow field and wind pressure coefficient were analyzed at the heights of 30 m and 80 m respectively. The velocity vector is analyzed at a height of 50 m, the flow rate and trace change were analyzed at the height of 80 m. The center and the end of the cable were selected to analyze the pressure distribution around the cable, and the numerical simulation was carried out. The calculation results showed that, the characteristics of fluid distribution in the flow field are in accordance with the conclusion of similar research, and a reasonable accurate numerical value can be obtained. Vortex was formed on the upper and lower sides of the cable, and the eddy current moved along the central axis and took a symmetrical shape, the pressure distribution of the cable surface was consistent with that of the wind tunnel in the literature, and in line with the engineering reality.

Kaimal pulsating wind speed spectrum; three-dimensional wind field; RNG-model; numerical simulation

O368;U441+.2; TU312+.1

A

1000-2324(2018)05-0809-06

10.3969/j.issn.1000-2324.2018.05.017

2017-06-09

2017-06-25

國家高技術(shù)研究發(fā)展計劃(八六三)項目(2007AA11Z107);安徽省高校自然科學(xué)研究項目(KJ2016A448);安徽省高校自然科學(xué)研究項目(KJ2017A787)

孫亭亭(1986-),男,博士研究生,主要從事橋梁結(jié)構(gòu)工程方面的研究工作. E-mail:kunpengting@163.com

通訊作者:Author for correspondence.E-mail:whutvses@163.com