廖 斌,陳善群,任超洋
(安徽工程大學 建筑工程學院,安徽 蕪湖 241000)
潛堤上波浪傳播過程數(shù)值研究
廖 斌,陳善群,任超洋
(安徽工程大學 建筑工程學院,安徽 蕪湖 241000)
以RANS方程結(jié)合Level-Set方法為研究基礎(chǔ),采用五階WENO有限差分格式進行空間離散、三階TVD Runge-Kutta格式進行時間推進、Level-Set方法追蹤波浪與空氣間自由面以及解析松弛法來實現(xiàn)數(shù)值水槽中的造波消波,建立了一種求解潛堤上波浪傳播問題的數(shù)值計算模型。選取經(jīng)典的潛堤上波浪傳播物理試驗?zāi)P蛯?shù)值模型進行了驗證,水位計算結(jié)果與試驗值吻合較好。進一步研究了波高、潛堤頂部淹沒深度以及潛堤坡度等參數(shù)對潛堤上波浪傳播過程的影響。結(jié)果表明:波高越高、淹沒深度越淺、潛堤向波坡度越小,波浪受淺水作用越明顯;潛堤背波坡度越小,波浪受反淺水作用稍大,但并不明顯。
潛堤;Level-Set方法;波浪傳播;數(shù)值研究;淺水作用
潛堤作為一種常見的護岸構(gòu)筑物在防止波浪對海岸(尤其是開敞式海岸)的侵蝕過程中發(fā)揮著重要作用。研究潛堤上波浪傳播過程,特別是波浪與空氣間自由面的演化過程對波浪傳播物理機制的深入理解以及臨岸構(gòu)筑物設(shè)計等方面具有非常重要的理論和工程實際意義。
相對于試驗研究手段,數(shù)值模擬具有快捷、經(jīng)濟、不受計算范圍及尺度約束、重復(fù)性好、在保證精度的條件下容易得到一般規(guī)律性結(jié)論的特點,使得數(shù)值模擬在波浪傳播研究領(lǐng)域越來越受研究者的青睞。以往研究采用的數(shù)值模擬方法主要包括Boussinesq方程數(shù)值模型[1-3]、淺水方程數(shù)值模型[4-5]以及雷諾時均Navier-Stokes(簡稱RANS)方程結(jié)合界面追蹤數(shù)值模型[6-8]3種。其中,Boussinesq方程數(shù)值模型和淺水方程數(shù)值模型均沒有考慮波浪和空氣之間的相互作用,只適用于求解淺水效應(yīng)引起的波浪變形問題,不適用于碎波條件下波浪水位問題的求解。劉忠波等[9]將破碎項引入Boussinesq方程仍不能很好地解決碎波條件下波浪水位問題;采用RANS方程結(jié)合界面追蹤數(shù)值模型考慮了自由面演化過程中波浪和空氣之間的相互作用以及流體黏度和湍流流動的影響,具有解決碎波條件下波浪水位問題的能力,是一種較為直接的求解模型,但目前使用這種求解模型的方法大多基于商業(yè)軟件。商業(yè)軟件中關(guān)于空間及時間離散精度以及邊界條件設(shè)定的限制使得潛堤上波浪傳播過程的數(shù)值模擬精度受到一定影響。為了克服商業(yè)軟件的缺陷,本文擬建立一種求解潛堤上波浪傳播問題的二維數(shù)值計算模型,并進行驗證計算。用該模型進一步研究波高、潛堤頂部淹沒深度以及潛堤坡度等參數(shù)對潛堤上波浪傳播過程的影響。
2.1 Level-Set界面追蹤方法
Level-Set方法由Osher等[10]于1988年提出,其中心思想是將隨時間運動的2種物質(zhì)間界面Γ(t)看作某個函數(shù)φ(x,t)的零等值面(線),即構(gòu)造函數(shù)φ(x,t)使得在任意時刻t,運動界面Γ(t)恰好是φ(x,t)的零等值面(線),滿足關(guān)系式
(1)
式中:Ω表示整個計算域;φ(x,t)又被稱為Level-Set函數(shù)。
φ(x,t)的初值需滿足在Γ(t)附近法向單調(diào),在Γ(t)上始終為0。一般可取φ(x,t)為點x到t=0時刻界面Γ(0)的符號距離,關(guān)系式為
(2)
式中:d(x,Γ(0))表示點x到界面Γ(0)的最小距離;Ω1和Ω2分別表示第1種物質(zhì)和第2種物質(zhì)所處的計算域。
在任意時刻t,對于運動界面Γ(t)上的任意點x,都要保證φ(x,t)=0,從而函數(shù)φ(x,t)需滿足關(guān)系式
(3)
式(3)被稱為Level-Set方程,求解該方程即能得到任意時刻φ(x,t)的函數(shù)值,從而確定運動界面Γ(t)的位置。
由于在數(shù)值計算過程中容易產(chǎn)生誤差,根據(jù)式(3)計算得到的φ(x,t)可能不再是t時刻點x到界面的距離,這將失去φ(x,t)的零等值面就是運動界面的意義。因此,為了保持φ(x,t)的符號距離性質(zhì),需對φ(x,t)作重新初始化,使其重新成為點x到界面Γ(t)的符號距離。重新初始化過程通過求解式(4)的初值問題來實現(xiàn)。
(4)
2.2 控制方程
本文選取不可壓縮RANS方程,結(jié)合Level-Set方程求解潛堤上波浪傳播問題,二維控制方程表達式為
(5)
式中:ui,uj均為時均速度;ρ為波浪密度;p為壓強;ν為運動黏度;νT為渦黏系數(shù);g為重力加速度。壓強p和時均速度ui的耦合迭代求解采用預(yù)處理BiCGStab算法[11],湍流模型選用k-ω模型[12]。
2.3 空間離散和時間推進格式
本文選取的空間離散格式為Jiang等[13]提出的五階WENO有限差分格式,簡稱WENO5格式。WENO5格式將用于對流項中時均速度ui,k-ω湍流模型中的湍動能k和湍流耗散率ω,以及Level-Set函數(shù)φ的離散迭代求解。下面用Level-Set重新初始化過程為例說明WENO5格式的離散求解過程。
將式(4)進行有限差分,二維差分格式為
(6)
其中:
(7)
(8)
(9)
式中si,j為符號函數(shù)sign(φ)的簡寫。
式(7)中WENO5的定義式為
(10)
其中:
(11)
式中:IS1,IS2,IS3被稱為光滑因子;ε表示一不為0的小量。時間推進格式為Shu等[14]提出的三階TVDRunge-Kutta格式,表達式為
(12)
其中:
2.4 解析松弛法
解析松弛法最早由Mayer等[15]提出,并在消波中使用,后經(jīng)由學者研究發(fā)現(xiàn)該方法可同時用于造波和消波。其基本原理是在數(shù)值水槽的前端和后端分別設(shè)立一個造波區(qū)和消波區(qū),相當于在計算區(qū)域添加進口和出口邊界條件,進口邊界條件只負責造波,而出口邊界條件只負責消波。造波區(qū)和消波區(qū)內(nèi)每一時刻對速度進行解析松弛處理,保證波浪在數(shù)值水槽中傳播時不受二次反射波以及透射波的影響。本文在對速度解析松弛處理基礎(chǔ)上,對Level-Set函數(shù)φ也作同樣處理。
解析松弛處理的關(guān)系式為
(13)
式中下標r,a,c分別代表解析松弛值、理論值以及當前計算值。
式(13)中的R(x)為解析松弛函數(shù),本文選取Jacobsen提出的松弛函數(shù),表達式為
(14)
式中x′為造波區(qū)或消波區(qū)內(nèi)水平方向坐標與造波區(qū)或消波區(qū)長度的比值。一般在造波區(qū)內(nèi)解析松弛函數(shù)取R(x′),消波區(qū)內(nèi)取R(1-x′)。
圖1 潛堤上波浪傳播數(shù)值水槽模型示意圖Fig.1 Schematic diagram of numerical flumemodel for simulating wave propagation overa submerged dike
為驗證第2節(jié)數(shù)值模型對潛堤上波浪傳播過程數(shù)值模擬的精確性,本文選取Beji等[16]建立的潛堤上波浪傳播過程物理試驗?zāi)P妥鳛樗憷?,圖1為物理試驗?zāi)P秃喕傻臐摰躺喜ɡ藗鞑?shù)值水槽模型。
整個數(shù)值水槽長38.0m,高0.80m,水槽中水位為0.40m。造波區(qū)和消波區(qū)分別位于數(shù)值水槽的前后兩端,梯形潛堤位于數(shù)值水槽中部的工作區(qū)。波浪由造波區(qū)生成后從左向右傳播至工作區(qū),經(jīng)過潛堤后最終進入消波區(qū)。梯形潛堤的前沿距造波區(qū)6.0m,向波坡面坡度1∶20,背波坡面坡度1∶10,潛堤高度0.30m,數(shù)值水槽中水位距潛堤頂部0.10m。水槽工作區(qū)內(nèi)共設(shè)置8個水位儀用于水位監(jiān)測,分別為WG1—WG8,以造波區(qū)為起點,對應(yīng)位置分別為6.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0m。數(shù)值模擬網(wǎng)格大小dx=dz=0.01 m,整個數(shù)值水槽網(wǎng)格單元數(shù)為304 000。為模擬試驗中的真實波浪,本文選取周期T=2.525 s,波高H=0.022 m的二階Stokes波為生成波。時間步大小滿足CFL條件:CFL數(shù)=0.1。
將水位儀(WG1—WG8)所處位置的波浪水位數(shù)值計算結(jié)果與Beji等[16]的試驗結(jié)果進行對比,結(jié)果如圖2所示。
圖2 潛堤上波浪傳播過程中水位數(shù)值計算與試驗結(jié)果的對比Fig.2 Comparison between numerical results and experimental data of water level in the process of wave propagation over a submerged dike
圖2表明:
(1) 圖2(a)中水位儀WG1監(jiān)測的是波浪傳播至潛堤向波坡面起點的水位。由于潛堤向波坡面的存在,波浪從深水區(qū)往淺水區(qū)傳播過程中受到淺水作用而出現(xiàn)變形。WG1數(shù)值計算水位波形水平方向不對稱并且波峰高度比波谷深度要稍微大一點。
(2) 圖2(b)中水位儀WG2監(jiān)測的是波浪傳播至潛堤向波坡面上的水位。WG2計算水位波形與WG1類似,但波峰高度和波谷深度都有增加,其中波峰高度增加更為明顯。
(3) 圖2(c)中水位儀WG3監(jiān)測的是波浪剛剛傳播至潛堤頂部平臺時的水位。WG3數(shù)值計算水位波形較WG2而言,波峰高度和波谷深度稍微增加,同樣波峰高度增加更加明顯。
(4) 圖2(d)中水位儀WG4監(jiān)測的是波浪沿潛堤平臺傳播時的水位。WG4數(shù)值計算水位波形的主波峰之后出現(xiàn)了2次波峰,并且主波峰高度較WG3進一步增加。
(5) 圖2(e)中水位儀WG5監(jiān)測的是波浪剛剛從潛堤頂部平臺傳播至背波坡面位置時的水位。WG5數(shù)值計算水位較WG4變化明顯,出現(xiàn)明顯的2次波峰且主波峰高度達到最大值。
(6) 圖2(f)和圖2(g)中水位儀WG6和WG7監(jiān)測的是波浪傳播至背波坡面上的水位。值得注意的是,由于潛堤背波坡面的存在,此后波浪從淺水區(qū)往深水區(qū)傳播并受到反淺水作用,使得WG6和WG7數(shù)值計算水位波形出現(xiàn)了3次波峰,WG7的3次波峰高度較WG6要高。
(7) 圖2(h)中水位儀WG8監(jiān)測的是波浪傳播至潛堤背波坡面終點時的水位。WG8數(shù)值計算水位波形中3次波峰高度較WG7進一步增加。
整個數(shù)值計算的各個水位監(jiān)測點的水位波形與Beji等[16]的試驗結(jié)果非常相符,表明本文建立的數(shù)值模型對潛堤上波浪傳播過程的模擬具較高精度。
為進一步研究潛堤上波浪的傳播過程,通過調(diào)整波浪波高、潛堤上淹沒深度、潛堤向波坡面坡度以及背波坡面坡度等參數(shù),深入分析各種參數(shù)對波浪傳播的影響。
4.1 波浪波高
圖3給出了波高H=0.011,0.022,0.044 m條件(其余計算條件與第3節(jié)一致)下潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果。
圖3 不同波高條件下潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果Fig.3 Numerical results of water level in the process of wave propagation over a submerged dike in the presence of varying wave height
從圖3中可以看出,不同波高條件下WG2,WG3,WG5,WG7數(shù)值計算水位波形明顯不同,具體表現(xiàn)為:①由于從深水區(qū)往淺水區(qū)傳播過程中波浪受到淺水作用,波高H=0.044 m的波浪受淺水作用影響較大,WG2位置波高H=0.044 m的水位波形水平方向不對稱較為明顯,且波峰高度和波谷深度都近似與H=0.011,0.022 m的數(shù)值計算結(jié)果成倍數(shù)關(guān)系。②WG3位置波高H=0.044 m的水位波形波峰高度進一步增加,且波形有2次波峰出現(xiàn),相比之下波高H=0.011,0.022 m波峰高度增加得并不明顯,波形變化較小。③WG5位置波高H=0.044 m的水位波形已出現(xiàn)明顯的2次波峰,波浪主波峰高度明顯下降。波高H=0.022 m的水位波形出現(xiàn)明顯的2次波峰,H=0.011 m水位波形有2次波峰出現(xiàn)的趨勢。④WG7位置波浪受反淺水作用,波高H=0.044 m水位波形明顯出現(xiàn)3次甚至4次波峰,H=0.022 m水位波形出現(xiàn)明顯的3次波峰,H=0.011 m水位波形有3次波峰出現(xiàn)的趨勢。
圖4 不同淹沒深度條件下潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果Fig.4 Numerical results of water level in the process of wave propagation over a submerged dike in the presence of varying submergence depth
4.2 潛堤頂部淹沒深度
圖4給出了潛堤頂部淹沒深度d=0.05,0.1,0.2 m條件(其余計算條件與第3節(jié)一致)下潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果。
從圖4中可以看出,不同淹沒深度條件下WG2,WG3,WG4,WG6數(shù)值計算水位波形明顯不同,具體表現(xiàn)為:①由于從深水區(qū)往淺水區(qū)傳播過程中波浪受到淺水作用,淹沒深度d=0.05 m時波浪受淺水作用影響較大,WG2位置淹沒深度d=0.05 m的水位波形水平方向不對稱較為明顯,且波峰高度和波谷深度均比d=0.1,0.2 m的數(shù)值計算結(jié)果要大。值得注意的是,淹沒深度d=0.05 m條件下波浪傳播的水平速度要更大。②WG3位置淹沒深度d=0.05 m的水位波形波峰高度進一步增加,且波形有2次波峰出現(xiàn),相比之下d=0.1,0.2 m波峰高度增加的并不明顯,波形變化較小。③WG4位置淹沒深度d=0.05 m的水位波形已出現(xiàn)明顯的2次波峰,波浪主波峰高度明顯下降。淹沒深度d=0.1 m的水位波形也有明顯的2次波峰出現(xiàn),d=0.2 m水位波形變化不大。④WG6位置淹沒深度d=0.05 m的水位波形受反淺水作用,出現(xiàn)3次甚至4次波峰。淹沒深度d=0.1 m水位波形的2次波峰高度明顯增加并出現(xiàn)3次波峰,d=0.2 m水位波形有2次波峰出現(xiàn)趨勢。
4.3 潛堤向波坡度
圖5給出了潛堤向波坡度1∶5,1∶10,1∶20條件下(其余計算條件與第3節(jié)一致)潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果。
圖5 不同向波坡度條件下潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果Fig.5 Numerical results of water level in the process of wave propagation over a submerged dike with varying seaward slope gradient
從圖5中可以看出,不同向波坡度條件下WG2,WG3數(shù)值計算水位波形明顯不同,而WG5,WG6水位波形近乎相同。其中,從深水區(qū)往淺水區(qū)傳播過程中波浪受到淺水作用,向波坡度1∶20時波浪受淺水作用影響較大,WG2與WG3位置水位波形水平方向不對稱較為明顯,波峰高度增加較多。值得注意的是,淹沒深度d=0.05 m條件下波浪傳播的水平速度要更大;相比之下坡度1∶5,1∶10的水位波形變化不明顯,波峰高度有些許增加。由于后續(xù)潛堤的水平平臺尺寸以及背波坡度都一致,波浪的后續(xù)傳播過程中外部條件影響相同,導(dǎo)致后續(xù)水位波形趨于相同。
4.4 潛堤背波坡度
圖6給出了潛堤背波坡度1∶5,1∶10,1∶20條件下(其余計算條件與第3節(jié)一致)潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果。
圖6 不同背波坡度條件下潛堤上波浪傳播過程中水位的數(shù)值計算結(jié)果Fig.6 Numerical results of water level in the process of wave propagation over a submerged dike with varying leeward slope gradient
從圖6可以看出,不同背波坡度條件下WG2,WG4數(shù)值計算水位波形幾乎一致,而WG6,WG7水位波形有些許不同。其中,由于潛堤的向波坡度以及水平平臺尺寸都一致,波浪的前期傳播過程中外部條件影響相同,導(dǎo)致WG2,WG4數(shù)值計算水位波形幾乎一致。后續(xù)波浪從淺水區(qū)往深水區(qū)傳播過程中波浪受到反淺水作用,從圖6(c),圖6(d)中可以看出,背波坡度1∶20條件下受反淺水作用稍大,WG6,WG7水位波形水平方向不對稱稍微明顯。但總體來說,3種背波坡度條件下波浪受反淺水作用影響較為接近,WG6,WG7位置水位波形主波峰、2次及3次波峰高度相差不大。由于背波坡度1∶5條件下,WG7位置不在背波坡面上,所以背波坡度1∶5的WG7水位波形不參與比較。
本文從進一步提升潛堤上波浪傳播數(shù)值模擬的精度出發(fā),利用五階WENO有限差分格式以及三階TVD Runge-Kutta格式離散求解RANS方程,采用Level-Set方法追蹤波浪與空氣間自由面并采用解析松弛法實現(xiàn)造波消波,建立了一種求解潛堤上波浪傳播問題的數(shù)值模型。采用經(jīng)典的梯形潛堤試驗驗證數(shù)值模型精度,并進一步研究了波高、潛堤頂部淹沒深度以及潛堤坡度等參數(shù)對潛堤上波浪傳播過程的影響,得到以下結(jié)論:
(1) 數(shù)值模型計算的各個水位監(jiān)測點的水位波形與Beji等[16]的試驗結(jié)果非常相符,表明數(shù)值模型對潛堤上波浪傳播過程的模擬具有較高的精度,數(shù)值模擬結(jié)果具有較高的可靠性。
(2) 波高越高時波浪受淺水作用越明顯,波浪變形較為明顯,水位波形較早出現(xiàn)次級波峰;波高H=0.044 m條件下,水位波形甚至出現(xiàn)了4次波峰。
(3) 潛堤的淹沒深度越小波浪受淺水作用越明顯,此時波浪變形非常明顯;淹沒深度d=0.05 m條件下,水位波形同樣出現(xiàn)了4次波峰。
(4) 潛堤的向波坡面越小,淺水作用越明顯;而背波坡面坡度越小,受反淺水作用稍大,但并不明顯。
[1] BROCCHINI M, DRAGO M, IOVENITTI L. The Modelling of Short Waves in Shallow Waters. Comparison of Numerical Methods based on Boussinesq and Serre Equations[C]∥American Society of Civil Engineers. Proceedings of the 23rd International Conference on Coastal Engineering. Venice, Italy. October 4-9, 1992: 76-88.
[2] BEJI S, BATTJES J A. Numerical Simulation of Nonlinear Wave Propagation over a Bar[J]. Coastal Engineering, 1994, 23(1): 1-16.
[3] 孫 先, 楊文俊, 王國棟. 基于 坐標變換和水位函數(shù)法的立面二維水動力學模型[J]. 長江科學院院報, 2012, 29(7): 6-9, 26.
[4] KOBAYASHI N, OTTA A, ROY I. Wave Reflection and Run-up on Rough Slopes[J]. Journal of Waterway, Port, Coastal and Ocean Engineering, 1987, 113(3): 282-298.
[5] 陳善群, 廖 斌, 李海峰. 用最小二乘無網(wǎng)格有限差分方法求解二維淺水方程[J]. 長江科學院院報, 2012, 29(6): 36-39, 57.
[6] 李昌良, 謝媛媛. 不同結(jié)構(gòu)形式潛堤上的隨機波浪運動[J]. 海岸工程, 2008, 27(1): 1-9.
[7] 蔣昌波, 曹永港, 陳 杰, 等. 斜坡上梯形潛堤附近流動結(jié)構(gòu)的數(shù)值研究[J]. 水動力學研究與進展A輯, 2009, 24(3): 286-295.
[8]JACOBSEN N G, FUHRMAN D R, FREDS?E J. A Wave Generation Toolbox for the Open-source CFD Library: OpenFOAM[J]. International Journal for Numerical Methods in Fluids, 2011, 70(9): 1073-1088.
[9] 劉忠波, 于德海, 孫昭晨. 潛堤上破碎波浪傳播變形的數(shù)值模型及其驗證[J]. 海洋通報, 2011, 30(6): 633-636.[10]OSHER S,SETHIAN J A.Fronts Propagating with Curvature-dependent Speed:Algorithms based on Hamil- ton-Jacobi Formulations[J].Journal of Computational Physics, 1988, 79(1): 12-49.
[11]VAN DER VORST H A. Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems[J]. SIAM Journal on Scientific and Statistical Computing, 1992, 13(2): 631-644.
[12]WILCOX D C. Turbulence Modeling for CFD[M]. Los Angeles: DCW Industries Inc., 1994.
[13]JIANG G S, SHU C W. Efficient Implementation of Weighted ENO Schemes[J]. Journal of Computational Physics, 1996, 126(1): 202-228.
[14]SHU C W, OSHER S. Efficient Implementation of EssentiallyNon-oscillatoryShockCapturingSchemes[J].Journal of Computational Physics, 1988, 77(2): 439-471.
[15]MAYER S,GARAPON A,S?RENSEN L S.A Fractional StepMethodforUnsteadyFreeSurfaceFlowwithAppli- cationstoNon-linearWaveDynamics[J].International Journal for Numerical Methods in Fluids, 1998, 28(2): 293-315.
[16]BEJI S, BATTJES J A. Experimental Investigation of Wave Propagation over a Bar[J]. Coastal Engineering, 1993, 19(1/2): 151-162.
(編輯:黃 玲)
Numerical Investigation of Wave Propagation overa Submerged Dike
LIAO Bin, CHEN Shan-qun, REN Chao-yang
(College of Civil Engineering and Architecture, Anhui Polytechnic University, Wuhu 241000, China)
In this study, a computational model is established for solving wave propagation over a submerged dike based on RANS equations combined with Level-Set method. The fifth-order finite difference WENO scheme is used for spatial discretization, and a TVD third-order Runge-Kutta explicit time scheme is employed for time discretization in the model. The Level-Set method is used for tracking the free interface between the air and water phases, and a relaxation method is employed for wave generation and absorption. In order to validate the accuracy and applicability of the model, numerical investigation of the wave propagation over a submerged dike is conducted. The numerical results show a good agreement with experimental data. Further studies are carried out to investigate the influence of physical parameters, such as wave height, submerged depth, seaward and leeward slope gradients, on the process of wave propagation over a submerged dike. Results reveal that when the wave height is higher, submerged depth smaller, and seaward slope flatter, the effect of shoaling is more obvious; when leeward slope is flatter, the effect of shoaling on wave is slightly larger, but not obvious.
submerged dike; Level-Set method; wave propagation; numerical investigation;shoaling
2016-04-25;
2016-06-16
國家自然科學基金項目(51409001);安徽省自然科學基金項目(1508085QE100);安徽高校優(yōu)秀青年人才支持計劃項目(gxyq2017015)
廖 斌(1985-),男,江西撫州人,高級實驗師,碩士,研究方向為流體力學,(電話)13515534139(電子信箱)liaobinfluid@126.com。
10.11988/ckyyb.20160393
2017,34(8):52-58
TV131.2
A
1001-5485(2017)08-0052-07