葉瀟瀟 游景皓 宋金寶

實驗室造波條件對內(nèi)孤立波發(fā)展影響的直接數(shù)值模擬*

葉瀟瀟1游景皓2宋金寶1①

(1. 浙江大學(xué)海洋學(xué)院物理海洋與遙感研究所 舟山 316021; 2. 四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室 成都 610065)

內(nèi)孤立波具有振幅尺度大、能量集中的特點, 其引起流場和密度場的迅速變化可能對海洋工程結(jié)構(gòu)物以及水下潛體造成嚴重威脅。因此研究不同造波條件下生成的內(nèi)孤立波運動的流場特征具有重要的學(xué)術(shù)意義和實際應(yīng)用價值。采用直接數(shù)值模擬方法和給定的初始密度場密度躍遷函數(shù), 對重力塌陷激發(fā)內(nèi)孤立波的運動過程進行研究, 探討了不同造波條件下, 激發(fā)產(chǎn)生的內(nèi)孤立波波型、渦度、振幅和水平速度等流場特征。結(jié)果表明: (1)直接模擬數(shù)值方法能夠模擬內(nèi)孤立波傳播過程中的密度界面波型反轉(zhuǎn)現(xiàn)象; (2)從定性和定量的角度, 證實了不穩(wěn)定內(nèi)孤立波傳播過程中存在能量的向后傳遞; (3) 對于相同的臺階深度(水閘兩側(cè)初始密度界面的高度差), 初始渦流保持相同, 但是隨著上下層水深比的減小, 其強度下降顯著; (4)臺階深度對初始渦流的垂直結(jié)構(gòu)的影響要大于上下層水深比, 且臺階深度對內(nèi)孤立波的振幅、水平速度的影響顯著。

內(nèi)孤立波; 初始條件; 重力塌陷; 直接數(shù)值模擬

內(nèi)波是發(fā)生在密度穩(wěn)定分層的流體內(nèi)部的一種波動, 而內(nèi)孤立波是一種特殊的內(nèi)波。內(nèi)孤立波頻繁地發(fā)生在海洋中, 在世界海洋的許多大陸邊緣斜坡上都曾觀測到內(nèi)孤立波痕跡(Filonov, 2000; Azevedo, 2006; Vázquez, 2008)。內(nèi)孤立波具有強垂向流和水平流, 受大陸或近岸斜坡影響的內(nèi)孤立波會淺化、破碎, 并造成強烈的湍流混合(Forgia, 2018), 進而影響溫度(Leichter, 1996; Davis, 2011; Walter, 2012)、氧氣和營養(yǎng)物質(zhì)的空間分布(Leichter, 1996; Omand, 2011; Walter, 2014)。

內(nèi)孤立波對于結(jié)構(gòu)物具有潛在的破壞性。在非線性效應(yīng)和頻散效應(yīng)的平衡作用下, 在傳播過程中, 內(nèi)孤立波的波形和速度可以保持數(shù)百公里不變(Grue, 2000; Grimshaw, 2002)。在中國南海(South China Sea, SCS)曾記錄到振幅最大為170 m的內(nèi)孤立波, 且上、下層之間的水平速度差超過3.4 m/s (Warn- Varnas, 2010)。正是由于內(nèi)孤立波具有振幅尺度大、能量集中的特點, 其引起流場和密度場的迅速變化可能對海洋工程的結(jié)構(gòu)物以及水下潛體造成嚴重威脅(鄒麗等, 2020)。由此可見, 內(nèi)孤立波研究在海洋工程、海洋觀測等方面具有十分重要的應(yīng)用意義。因此研究不同造波條件下生成的內(nèi)孤立波運動的流場特征具有重要的學(xué)術(shù)意義和實際應(yīng)用價值。

在Chen等(2007)的實驗室實驗中, 內(nèi)孤立波是在水閘兩側(cè)流體界面高度不同的坍塌機制下生成的(由水閘劃分為左側(cè)環(huán)境區(qū)和右側(cè)擾動區(qū))。目前在內(nèi)孤立波的實驗室研究過程中, 在分層水體制取方面, 要配制出質(zhì)量高的穩(wěn)定的分層水體, 需保證流入水槽流體運動的平穩(wěn)性及均勻性, 使得配制過程要耗用很長的時間, 而且試驗成功率低, 需經(jīng)多次重新配制才能達到要求(徐鑫哲, 2012)。此外, 由實驗儀器收集到的信息也往往比較有限, 即使使用了高頻相機和超聲波測試儀, 依然難以得出流場內(nèi)部精確的密度、渦度的時空變化。

數(shù)值模擬是一種可以獲得精確的流場內(nèi)部變化信息的有效方法。對內(nèi)孤立波的數(shù)值模擬研究基本都是從數(shù)值求解Navier-Stokes (N-S)方程入手, 根據(jù)求解N-S方程的方法可進一步分為: (1)直接求解N-S方程(直接數(shù)值模擬方法, direct numerical simulation, DNS); (2)求解雷諾平均的N-S方程(reynolds-averaged Navier-Stokes simulation, RANS) (Vlasenko, 2002; Lin, 2012); (3)求解空間濾波的N-S方程(大渦數(shù)值模擬法, large-eddy simulation, LES) (沈治等, 2009)。然而, RANS的模擬方法只能提供湍流的平均信息, 在模擬結(jié)果中, 上凸型內(nèi)孤立波基本不呈現(xiàn)渦旋結(jié)構(gòu), 這與實驗室實驗結(jié)果不符。LES的模擬方法適用于研究運動尺度較大的湍流。目前有關(guān)直接數(shù)值模擬方法的內(nèi)孤立波研究相對較少, 理論上來說, 直接數(shù)值模擬方法能夠更好地模擬出內(nèi)孤立波發(fā)展和傳播過程中的細微結(jié)構(gòu), 如因斷裂形成的湍流渦旋結(jié)構(gòu)。

本文參考Chen等(2007)重力塌陷激發(fā)內(nèi)孤立波的水槽實驗, 構(gòu)建了重力塌陷激發(fā)內(nèi)孤立波的數(shù)學(xué)模型。根據(jù)Druzhinin等(2015)的密度躍遷公式計算初始密度場。以Boussinesq假定下的N-S方程組及鹽分對流擴散方程作為控制方程, 并采用有限差分法進行離散, 實現(xiàn)了重力塌陷造波方式下形成的內(nèi)孤立波運動的數(shù)值模擬。通過與Chen等(2007)的實驗結(jié)果對比, 驗證了本文方法的有效性, 而后本文進一步與Lin等(2012)對比不同數(shù)值模擬方法下呈現(xiàn)的渦旋結(jié)構(gòu)。之后, 通過改變造波參數(shù), 即流體密度差、上下層水深比和臺階深度(水閘兩側(cè)初始密度界面的高度差), 對一系列算例進行數(shù)值模擬, 得到流場的密度界面、渦旋結(jié)構(gòu)、渦度等的時空演化結(jié)果, 比較不同造波條件對內(nèi)波后續(xù)發(fā)展過程中的波型、渦度、水平速度、流場動能和特征波長的影響。

1 數(shù)學(xué)模型

1.1 控制方程

本文以不可壓縮N-S方程組(連續(xù)性方程及動量守恒方程)和對流擴散方程為控制方程, 三個方程分別為:

圖1 模型示意圖

Fig.1 Schematic of the model

注:: 水箱高度;: 水箱長度;s: 斜坡長度;p: 擾動區(qū)長度;1: 左側(cè)環(huán)境區(qū)域的上層流體水深;2: 左側(cè)環(huán)境區(qū)域的下層流體水深;1: 右側(cè)擾動區(qū)域的上層流體水深;2: 右側(cè)擾動區(qū)域的下層流體水深;: 波高;w: 波長;1: 上層流體密度;2: 下層流體密度;: 流體運動黏度;0: 臺階深度

表1 各算例主要參數(shù)表

Tab.1 Main parameters in the simulation cases

控制方程中的鹽分對流擴散方程用迎風(fēng)緊致差分格式(upwinding combined compact difference, UCCD)離散(Yu, 2019)。N-S方程, 利用二階顯式Adams-Bashforth法離散時間項, 二階迎風(fēng)格式離散對流項, 以二階中心差分格式離散黏性擴散項, 壓力項通過超松弛迭代法(successive over relaxation method, SOR)求解(Yu, 2019)。

表2 控制方程中參數(shù)的有量綱形式和無量綱形式

Tab.2 Dimensional and dimensionless forms of parameters in the governing equations

1.2 初始及邊界條件

本文采用了密躍層假設(shè)(Druzhinin, 2015), 認為流體間的密度是按雙曲正切函數(shù)變化。相對于簡單的線性分層, 雙曲正切函數(shù)具有階躍突變的特點, 這與實際海洋中的密躍情況相接近。本文設(shè)置計算流體的初始鹽度分布如下:

對于速度邊界條件, 在頂?shù)准皟蓚?cè)界面處均使用了自由滑移邊界條件(free-slip boundary condition)。以頂部為例:

對于鹽度邊界條件, 在界面處均使用了無通量條件。

2 模型檢驗

為檢驗本文數(shù)值方法的有效性, 將本文模擬結(jié)果和前人的實驗(Chen, 2007)及模擬結(jié)果(Lin, 2012)進行比較。

本文與Chen等(2007)實驗結(jié)果進行了兩方面的比對, 分別為界面波形的演化過程和定點波動時間序列。圖2為界面波形的演化過程, 圖2a為Chen等(2007)實驗的結(jié)果, 圖2b為本文的模擬結(jié)果。由圖2可以看出, 實驗與模擬方式之間最明顯的差異是: 在實驗水槽內(nèi)填充階段引入了厚度較小的水閘, 以便在左右兩側(cè)流體區(qū)域之間產(chǎn)生不同的分層。實驗中水閘的去除會引起垂直剪切應(yīng)力, 并將力傳遞到與其直接接觸的流體。但考慮到水閘的厚度較小, 其所產(chǎn)生的干擾應(yīng)是相當有限的。圖2證實, 在去除水閘不久后, 擾動區(qū)附近模擬得到的密度界面便與通過實驗觀察到的染料的空間分布非常吻合: 在無量綱時間=3.5, 最初放置在擾動區(qū)內(nèi)的較重的流體開始呈現(xiàn)出內(nèi)孤立波的典型形狀; 在=6.22至=10.5, 內(nèi)孤立波向左傳播, 同時保留其形狀; 在=6.22至=15, 在實驗和模擬中, 均存在界面尾隨內(nèi)波。

圖2 界面波形的演化過程

注: a: Chen等(2007)實驗結(jié)果; b: 本文模擬結(jié)果(表1中工況Chen1, 實際距離10.75—12 m),表示無量綱時間

總體而言, 本文的直接數(shù)值模擬結(jié)果與Chen等(2007)實驗結(jié)果所得的流場界面波形具有較好的一致性。

圖3 定點波動時間序列對比圖(表1中工況Chen2)

本文與Lin等(2012)的RANS模擬結(jié)果進行了界面渦旋結(jié)構(gòu)對比。如圖4所示, 在擾動區(qū)水體剛開始下落時, 本文與Lin等(2012)模擬結(jié)果中的密度界面位移變化基本一致。隨著時間的推移, 兩種模擬方法開始有所區(qū)別: 在Lin等(2012)模擬結(jié)果中, 上凸型內(nèi)波基本不呈現(xiàn)渦旋結(jié)構(gòu), 這與實驗室實驗結(jié)果不符; 本文模擬結(jié)果則呈現(xiàn)了因斷裂形成湍流渦旋結(jié)構(gòu)的詳細過程, 這為后續(xù)對比內(nèi)孤立波的起渦能力提供了充分條件。

圖4 本文模擬結(jié)果(表1中工況Lin)與Lin等(2012)模擬結(jié)果的界面渦旋結(jié)構(gòu)對比

注: a:=0.14; b:=0.27; c:=0.41; d:=0.54;表示無量綱高度,表示無量綱水平距離

綜上, 對驗證算例的模擬結(jié)果表明, 本文所采用的直接模擬數(shù)值模型有效地模擬了重力塌陷式內(nèi)孤立波的波動變化過程, 相比于Lin等(2012)的RANS模擬結(jié)果呈現(xiàn)出更為明顯的渦旋結(jié)構(gòu)。

3 模擬結(jié)果分析

為研究造波條件對重力塌陷式內(nèi)孤立波發(fā)展的影響, 本文設(shè)計了8組算例(工況1—工況8), 各算例的主要參數(shù)參見表1。由于本文主要考慮的是造波條件對內(nèi)孤立波后續(xù)發(fā)展過程中的波型、渦度、振幅和水平速度的影響, 而根據(jù)驗證算例和Chen等(2007)實驗可知, 如果斜坡與擾動區(qū)距離較遠, 在內(nèi)孤立波沒有傳遞到斜坡之前, 斜坡的影響基本可以忽略, 因此本文沒有在這些算例中加入斜坡。

3.1 流場渦度

在研究過程中發(fā)現(xiàn)流體密度差對無量綱速度幾乎沒有影響, 因此本節(jié)將不給出工況6、工況7流場渦度對比圖。圖5a、5b、5c為不同臺階深度的流場渦度時間演化的對比圖。結(jié)果表明, 臺階深度0對流場渦度的影響較為顯著。在無量綱時間=2.5時, 工況3幾乎沒有出現(xiàn)正渦度, 在之后的時間=3.5至=5.5, 工況1、工況2均表現(xiàn)出正渦度先增大后衰減的過程, 工況3的衰減過程則相對不明顯。這表明臺階深度越大, 對尾部的湍流振蕩影響越顯著, 湍流振蕩越劇烈, 且振蕩持續(xù)時間相對越長。除此之外, 臺階深度越大, 產(chǎn)生的內(nèi)孤立波波高越大, 內(nèi)孤立波傳播的平均速度越大。

圖5a、5f為不同水層條件下(分別為穩(wěn)定和不穩(wěn)定的內(nèi)孤立波)的流場渦度時間演化的對比圖。結(jié)果表明, 工況8的負渦度衰減的更快, 且正渦度相對工況1數(shù)值更大。這定性地表明了不穩(wěn)定內(nèi)孤立波在向前傳播的過程中, 存在能量的向后部傳遞。

3.2 最大水平速度

(1) 增大流體密度差, 流場最大水平速度將隨之增大;

(2) 隨著臺階深度的增大, 流場最大水平速度將隨之增大; 臺階深度較大的算例(工況3)初始會有一段最大水平速度保持不變的階段, 這可能是因為增大臺階深度會加長初始水體順時針渦流運動的持續(xù)時間;

圖5 不同算例流場渦度的時間演化

圖6 不同算例流場最大水平速度隨水平位置的變化

3.3 最大波高及振幅

圖7 不同算例最大波高Hmax的時間演化

(1) 流體密度差對最大波高基本沒有影響;

(2) 隨著臺階深度的增大, 最大波高第二階段的反復(fù)現(xiàn)象越明顯, 即振蕩越劇烈;

圖8為內(nèi)孤立波的振幅隨各項造波參數(shù)的變化曲線。結(jié)果表明: 臺階深度對內(nèi)孤立波的振幅影響顯著, 上下層水深比對振幅影響相對較小, 流體密度差對振幅基本不產(chǎn)生影響。

3.4 波形反轉(zhuǎn)

圖9為工況1流場界面內(nèi)孤立波的演化過程。穩(wěn)定的內(nèi)孤立波在傳播過程中, 其前部內(nèi)孤立波處的順時針渦流始終維持較大速度, 這表明能量集中在前部波內(nèi)。圖10為工況8流場界面內(nèi)孤立波波型反轉(zhuǎn)的演化過程。與圖9對比, 可清楚觀察到工況8的內(nèi)孤立波在傳播過程中發(fā)生波型不穩(wěn)定轉(zhuǎn)換的現(xiàn)象:

(2) 不穩(wěn)定的內(nèi)孤立波在傳播過程中, 前部上凸型內(nèi)孤立波處的順時針渦流速度逐漸減小; 反之, 后部的下凹型內(nèi)孤立波處的逆時針渦流速度逐漸增大。這證實了不穩(wěn)定內(nèi)孤立波在傳播過程中會發(fā)生內(nèi)部能量傳遞, 能量由前部上凸型內(nèi)孤立波, 逐漸傳遞到后部下凹型內(nèi)孤立波。

以上均是定性判斷了內(nèi)孤立波運動過程中的能量變化, 下面將通過給出內(nèi)孤立波動能的定量計算, 分析內(nèi)孤立波運動過程中的能量變化。

圖9 穩(wěn)定內(nèi)孤立波界面波形的演化過程(表1中工況1)

圖10 不穩(wěn)定內(nèi)孤立波界面波形的演化過程 (表1中工況8)

結(jié)果表明, 在無量綱時間=12.5至=16的這段時間內(nèi), 穩(wěn)定內(nèi)孤立波(工況1)的流場動能水平分布形態(tài)基本不發(fā)生改變, 第一個轉(zhuǎn)折點處判斷為前部波波峰位置; 不穩(wěn)定內(nèi)孤立波(工況8)的流場動能水平分布形態(tài)發(fā)生明顯改變, 對比=12.5與=16的第一個轉(zhuǎn)折點處的動能占比, 可看到隨時間的推移, 前部流體動能占比減小。這表明在不穩(wěn)定內(nèi)孤立波在向前傳播的過程中, 存在能量的向后部傳遞。除此之外, 工況8的后部動能占比明顯大于工況1, 這可能是由于工況8的上下層水深比接近1, 不穩(wěn)定內(nèi)孤立波在主波后方產(chǎn)生了連續(xù)的內(nèi)波波列, 造成能量分散。

圖11 不同時間的流場動能的水平分布

注:kr: 不同水平位置處的流體動能占整個流場動能的比例;: 無量綱水平距離;: 無量綱時間

4 結(jié)論

本文通過建立重力塌陷激發(fā)內(nèi)孤立波運動的數(shù)值模型, 采用直接模擬數(shù)值方法, 研究了不同造波條件對內(nèi)孤立波發(fā)展的影響, 結(jié)果表明:

(1) 通過直接數(shù)值模擬方法得到的流場界面呈現(xiàn)更為明顯的渦旋結(jié)構(gòu)。此外, 從定性和定量的角度, 證實了不穩(wěn)定內(nèi)孤立波傳播過程中存在能量的向后傳遞。

(2) 增大臺階深度、流體密度差、上下層水深比, 流場最大水平速度都將隨之增大, 其中臺階深度的影響最為顯著。

(3)對于相同的臺階深度, 初始渦流保持相同, 但是隨著上下層水深比的減小, 其強度明顯下降; 臺階深度對初始渦流的垂直結(jié)構(gòu)的影響要大于上下層水深比。

(4)臺階深度對內(nèi)孤立波的振幅影響顯著, 上下層水深比對振幅影響相對較小, 流體密度差對振幅基本不產(chǎn)生影響。

鄒 麗, 馬鑫宇, 李振浩等, 2020. 分層流體孤立子模型的實驗重現(xiàn)及流場分析. 哈爾濱工程大學(xué)學(xué)報, 41(2): 263—270

沈 治, 李宇鵬, 崔桂香等, 2009. 穩(wěn)定分層湍流的大渦模擬. 中國科學(xué) G輯: 物理學(xué)力學(xué)天文學(xué), 39(10): 1503—1513

王德軍, 唐 云, 于洪川等, 2003. 水平集方法與距離函數(shù). 應(yīng)用數(shù)學(xué)和力學(xué), 08: 839—848

徐鑫哲, 2012. 內(nèi)波生成機理及二維內(nèi)波數(shù)值水槽模型研究. 哈爾濱: 哈爾濱工程大學(xué)碩士學(xué)位論文

Azevedo A, da Silva J C B, New A L, 2006. On the generation and propagation of internal solitary waves in the southern Bay of Biscay. Deep Sea Research Part I: Oceanographic Research Papers, 53(6): 927—941

Chen C Y, Hsu J R C, Chen C W, 2007. Generation of internal solitary wave by gravity collapses. Journal of Marine Science and Technology, 15(1): 1—7

Davis K A, Monismith S G, 2011. The modification of bottom boundary layer turbulence and mixing by internal waves shoaling on a barrier reef. Journal of Physical Oceanography, 41(11): 2223—2241

Druzhinin O A, Ostrovsky L A, 2015. Dynamics of turbulence under the effect of stratification and internal waves. Nonlinear Processes in Geophysics, 22(3): 337—348

Filonov A E, Trasvi?a A, 2000. Internal waves on the continental shelf of the gulf of Tehuantepec, Mexico. Estuarine, Coastal and Shelf Science, 50(4): 531—548

Forgia G L, Tokyay T, Adduce C, 2018. Numerical investigation of breaking internal solitary waves. Physical Review Fluids, 3(10): 104801

Grimshaw R, Pelinovsky E, Poloukhina O, 2002. Higher-order Korteweg—de Vries models for internal solitary waves in a stratified shear flow with a free surface. Nonlinear Processes in Geophysics, 9(3—4): 221—235

Grue J, Jensen A, Rus?s P O, 2000. Breaking and broadening of internal solitary waves. Journal of Fluid Mechanics, 413: 181—127

He Z G, Zhao L, Hu P, 2018. Investigations of dynamic behaviors of lock-exchange turbidity currents down a slope based on direct numerical simulation. Advances in Water Resources, 119: 164—177

Leichter J J, Wing S R, Miller S L, 1996. Pulsed delivery of subthermocline water to Conch Reef (Florida Keys) by internal tidal bores. Limnology and Oceanography, 41(7): 1490—1501

Lin Z H, Song J B, 2012. Numerical studies of internal solitary wave generation and evolution by gravity collapse. Journal of Hydrodynamics, 24(4): 541—553

Omand M M, Leichter J J, Franks P J S, 2011. Physical and biological processes underlying the sudden surface appearance of a red tide in the nearshore. Limnology and Oceanography, 56(3): 787—801

Vázquez A, Bruno M, Izquierdo A, 2008. Meteorologically forced subinertial flows and internal wave generation at the main sill of the Strait of Gibraltar. Deep Sea Research Part I: Oceanographic Research Papers, 55(10): 1277—1283

Vlasenko V, Hutter K, 2002. Numerical Experiments on the breaking of solitary internal waves over a slope-shelf topography. Journal of Physical Oceanography, 32(6): 1779—1793

Walter R K, Woodson C B, Arthur R S, 2012. Nearshore internal bores and turbulent mixing in southern Monterey Bay. Journal of Geophysical Research: Oceans, 117(C7): C07017

Walter R K, Woodson C B, Leary P R, 2014. Connecting wind-driven upwelling and offshore stratification to nearshore internal bores and oxygen variability. Journal of Geophysical Research: Oceans, 119(6): 3517—3534

Warn-Varnas A, Hawkins J, Lamb K G, 2010. Solitary wave generation dynamics at Luzon Strait. Ocean Modelling, 31(1—2): 9—27

Yu C H, Zhao L, Wen H L, 2019. Numerical study of turbidity current over a three-dimensional seafloor. Communications in Computational Physics, 25(4): 1177—1212

DIRECT NUMERICAL SIMULATION OF THE EFFECT OF LABORATORY WAVE-MAKING CONDITIONS ON THE DEVELOPMENT INTERNAL SOLITARY WAVES

YE Xiao-Xiao1, YOU Jing-Hao2, SONG Jin-Bao1

(1. Institute of Physical Oceanography and Remote Sensing, Ocean College, Zhejiang University, Zhoushan 316021, China; 2. State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China)

Internal solitary waves are characteristic of large scale and concentrated energy, and could cause rapid changes in the flow field and density field, which poses a threat to marine engineering structures and underwater submersibles. Therefore, studying the flow field characteristics of internal solitary wave motion generated under different conditions has important academic significance and practical application value. Therefore, the direct numerical simulation method was conducted with a given density transition function of initial density field to study the motion process of gravity-collapsed internal solitary waves. In addition, the flow field characteristics of internal solitary waves generated in different wave-making conditions, such as waveform, vorticity, amplitude, and horizontal velocity, were discussed. Results show that the direct simulation numerical method used in this paper can simulate the density interface waveform inversion during internal solitary wave propagation. Meanwhile, from the qualitative and quantitative perspectives, backward energy transmission during the propagation of unstable internal solitary waves was confirmed. The initial vortex remained identical at the same step depth but the decrease in strength as the depth of upper layer increases was significant. The step depth influenced the vertical structure of initial vortex more than the depth ratio of the upper vs lower layer did on wave generation. Moreover, the step depth had a significant effect on the amplitude and horizontal velocity in wave generation.

internal solitary wave; initial conditions; gravity collapse; direct numerical simulation

* 國家重點基礎(chǔ)研究發(fā)展計劃項目, 2017YFA0604102號, 2016YFC1401404號。葉瀟瀟, 碩士研究生, E-mail: shawnye@zju. edu.cn

宋金寶, 博士生導(dǎo)師, 教授, E-mail: songjb@zju.edu.cn

2020-06-01,

2020-07-13

P733

10.11693/hyhz20200600157