縱搖容器中液體晃蕩的非線性數(shù)值模擬

寧德志，宋偉華，2，滕 斌

（1.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024；2.中交第四航務工程勘察設計院有限公司，廣州 510230）

縱搖容器中液體晃蕩的非線性數(shù)值模擬

寧德志1，宋偉華1，2，滕 斌1

（1.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024；2.中交第四航務工程勘察設計院有限公司，廣州 510230）

在海洋工程領(lǐng)域，液體晃蕩是一種普遍存在的物理現(xiàn)象。對于船舶而言，轉(zhuǎn)動比平動有著更重要的影響。該文針對縱搖容器中的液體晃蕩問題，采用高階邊界元方法建立自由水面滿足完全非線性邊界條件的時域數(shù)學模型。通過大地坐標系和隨體坐標系之間的坐標變換，使得計算域僅控制在隨體坐標系內(nèi)。求解中采用半混合歐拉—拉格朗日方法追蹤流體瞬時水面，運用四階龍格庫塔方法更新下一時間步的波面和速度勢。通過與已發(fā)表試驗和數(shù)值結(jié)果的對比，驗證了建立模型的準確性。進而開展大量數(shù)值試驗研究容器縱搖運動頻率、縱搖轉(zhuǎn)動中心和容器中布置一豎直隔板對晃動波面與荷載的影響。

液體晃蕩；縱搖；完全非線性；高階邊界元；隔板

0 引 言

在海洋工程領(lǐng)域中，液體晃蕩是一種普遍存在的物理現(xiàn)象。對于海上航行的船舶而言，縱搖運動對其穩(wěn)定性的影響遠遠大于平動運動。尤其是當運動幅值很大或外界運動頻率接近液艙的固有頻率時，液體晃蕩便會趨于劇烈，在內(nèi)壁上產(chǎn)生極大的作用力，從而對結(jié)構(gòu)造成危害，影響船舶的安全性和穩(wěn)定性。

對于縱搖液艙中的液體晃蕩，國內(nèi)外許多學者已經(jīng)進行了研究。Faltinsen[1]早在1974年就對二維矩形液艙作縱搖運動時的液體晃蕩問題進行了非線性理論分析。此非線性理論解析解為后來學者進行數(shù)值模擬提供了很好的參照。Nakayama和Washizu[2]采用數(shù)值方法對二維矩形液艙縱搖運動的液體晃蕩進行了模擬研究，他們利用有限元方法對空間進行網(wǎng)格劃分，采用有限差分方法來計算時間遞進。隨后他們又運用邊界元方法對二維縱搖矩形容器中的液體晃蕩進行了研究[3]。同時，Chen[4-8]對于液艙縱搖運動的液體晃蕩也做了大量工作。他針對理想流體的液體晃蕩問題展開研究，分析了運動頻率等因素的影響作用；同時還建立了粘性流模型，并將粘性流結(jié)果與理想流體結(jié)果進行對比分析。此外，他還針對復合運動液艙內(nèi)的液體晃蕩進行了研究，即液艙作平動和縱搖聯(lián)合運動。

近幾年來，Akyildiz和Unal[9-11]又針對矩形縱搖運動液艙中的液體晃蕩進行了研究。他們分別采取了物理實驗和數(shù)值模擬兩種方法對液體晃蕩產(chǎn)生的動水壓強進行了研究對比，并且考慮了縱搖的頻率和角度幅值對液體晃蕩的影響。除此之外，Armenio[12]、Bae[13]、Celebi[14]和Xue[15]等學者也在這方面做了大量的研究工作。同時，國內(nèi)學者肖龍飛等[16]、尹立中等[17]、洪亮等[18]、徐國徽等[19]也在這個領(lǐng)域有著較為深入的研究。然而，與平動運動液艙中液體晃蕩的研究相比，人們對縱搖運動液艙中液體晃蕩的研究還是相對較少，尤其是在減晃運動研究方面。

本文通過高階邊界元方法建立自由水面滿足完全非線性邊界條件的時域數(shù)學模型，針對液艙在不同轉(zhuǎn)動中心和運動頻率情況下的液體晃蕩進行了數(shù)值試驗，并分析了各因素的影響。在此基礎(chǔ)上，考慮到對液體晃蕩的抑制，本文還針對液艙中安置垂直隔板的情況進行了研究，給出了不同隔板位置、隔板長度情況下液體晃蕩的波面幅值和液艙所受壓力的變化趨勢，可為實際工程提供一定的參考借鑒。

1 數(shù)學模型

1.1 控制方程和邊界條件

如圖1所示，一個長度為L、水深為H的矩形容器在外界激勵作用下繞O點進行縱搖運動。定義兩個笛卡爾坐標系：一個是固定大地坐標系O0X0Y0；另一個是運動的隨體坐標系oxy。其中，x軸平行于容器靜止時的自由水面且向右為正方向，y軸和容器的中心線重合且向上為正方向。當容器靜止時，兩個坐標系是重合的。同時，定義靜止水面和x軸間的距離為e，且當x軸在靜止水面以下時e值取為正。

假設流體是無粘、不可壓縮且運動無旋的理想流體，則流體運動滿足Laplace方程，

式中：φ表示速度勢。

運動角速度ω（）t定義為

式中：θ（）

t表示液艙偏離靜止位置的角度，且順時針偏離為正值。

圖1 示意圖Fig.1 Definition sketch

同時假設：一是自由表面上的壓強等于大氣壓強的自由表面邊界條件；二是自由表面上的流體質(zhì)點將始終保持在自由表面上。本文中將大氣壓強設定為0?；诘谝粭l假設，自由表面的動力學邊界條件可寫為[2]

式中：η（x， ）

t是自由表面偏離靜止水面的位移，g是重力加速度?；诘诙l假設，自由表面的運動學邊界條件可寫為[2]

式中：nx和ny是自由表面上的外法向量與x、y軸夾角的余弦值。?/?n表示在外法向上的偏導數(shù)。

在固壁上，邊界條件可寫為如下形式[2]：

當t＝0時，容器靜止不動，相應的初始條件可寫為

此外，作用在液艙上的動水壓力可由Bernoulli方程求得：

進而通過在容器瞬時濕表面上進行壓強積分求出作用在容器上的荷載。

1.2 數(shù)值求解

為求解上節(jié)中的數(shù)學方程，在整個流域內(nèi)對速度勢應用格林第二定理，可得到如下邊界積分方程[20]：

式中：p為源點，q為場點，C為固角系數(shù)，S為流域邊界，包括自由水面邊界和固體邊界。本文用三節(jié)點高階邊界元離散計算邊界成一些曲線單元，單元內(nèi)任一點的幾何坐標和速度勢等物理量可以用形狀函數(shù)插值得到，這與曲面單元內(nèi)的做法是相同的。G是簡單格林函數(shù)，可寫為

式中：r1=，H為水深。

接下來，積分方程（8）經(jīng)過三節(jié)點曲線高階邊界單元離散后，便可以通過求解線性方程組得到相關(guān)未知量的值。計算中認為當前時刻物面Sn上的速度勢法向?qū)?shù)和自由水面Sf上的速度勢是已知的，根據(jù)積分方程計算當前時刻物面Sn上的速度勢和自由水面Sf上的速度勢法向?qū)?shù)，然后應用四階Runga-Kutta法，根據(jù)自由水面條件式（3）、（4）計算下一時刻的水質(zhì)點位置和自由水面Sf上的速度勢，再對自由水面重新劃分網(wǎng)格，重新應用積分方程計算下一時刻物面上的速度勢和自由水面上的速度勢法向?qū)?shù)。這樣計算周而復始，直到計算結(jié)束[21]。

2 模型驗證

作為算例，首先考慮一個長度為0.9m、水深為0.6m的矩形容器，運動方程為θ=θ0cosωt。其中，θ0是角度幅值，ω是運動頻率。選取θ0和ω為0.8度和5.5 rad/s。轉(zhuǎn)動中心設定在靜止水面上，即e＝0。由運動方程可知，當t＝0時，容器有一個0.8度的傾角。初始波面設為靜止不動，即初始波面與軸也有一個0.8度的夾角。通過數(shù)值收斂性分析，選取時間步長和網(wǎng)格大小為t＝0.01 s、△x＝0.045m、△z＝0.05m。圖2給出了容器右側(cè)壁處的波面時間歷程圖及本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果[2]的對比。由圖可知，本文數(shù)值結(jié)果與Nakayama和Washizu的數(shù)值結(jié)果吻合得非常好，且波面變化趨于共振。

接下來考慮另一個算例，矩形容器長度L＝0.92m、水深H＝0.31m。運動方程為θ=θ0sinωt。其中，θ0是角度幅值，ω是運動頻率。此例中，選取θ0和ω為8.0度和2.0 rad/s。轉(zhuǎn)動中心設定在靜止水面上，即e＝0。初始波面是靜止不動的。通過數(shù)值收斂性分析，選取時間步長和網(wǎng)格大小為t＝0.01 s，△x＝0.046m，△z＝0.052m。圖3給出了容器右側(cè)壁處的波面時間歷程圖及本文數(shù)值結(jié)果與Akyildiz實驗數(shù)據(jù)[11]、Akyildiz的CFD數(shù)值結(jié)果[11]的對比。由圖可知，本文數(shù)值結(jié)果比Akyildiz的數(shù)值結(jié)果更接近實驗數(shù)據(jù)。

圖2 余弦縱搖運動容器右側(cè)壁處的波面時間歷程Fig.2 Time series of surface elevation at the right wall subject to cosine pitch excitation

圖3 正弦縱搖運動容器右側(cè)壁處的波面時間歷程Fig.3 Time series of surface elevation at the right end subject to sine pitch excitation

圖4 右側(cè)壁底部上0.075m處的壓強時間歷程Fig.4 Time series of pressure ata point0.075m from the rightwall-bottom

圖5 t=5.52 s時右側(cè)壁上的壓強分布Fig.5 Distribution of pressure along the right wall at t=5.25 s

在上述波面驗證的基礎(chǔ)上，進一步開展數(shù)值模型的壓強驗證。選取容器長度L為1.0 m，水深H為0.35 m，運動方程為θ=θ0cosωt，θ0和ω為5.0度和5.19 rad/s。轉(zhuǎn)動中心在靜止水面以下0.1 m，即e＝0.1m。圖4給出了容器右側(cè)壁上一點的壓強時間歷程圖及本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果[3]的對比，這點距離液艙底0.075m。由圖可知，本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果吻合良好。同時，圖5給出了t＝5.52 s時液艙右側(cè)壁上的壓強分布圖及本文結(jié)果與Nakayama和Washizu數(shù)值結(jié)果[3]、Higuchi實驗數(shù)據(jù)[22]的對比。由圖可知，本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果吻合很好，與實驗數(shù)據(jù)對比也吻合良好。

通過上述算例計算，驗證了本文數(shù)值模型能夠準確模擬縱搖容器中液體晃蕩的波面變化和壓強變化。

3 計算及分析

3.1 運動頻率對液體晃蕩的影響

眾所周知，運動頻率在很大程度上決定著液體晃蕩的劇烈程度。當運動頻率接近固有頻率時，液體晃蕩將會趨于共振，在結(jié)構(gòu)上產(chǎn)生極大的壓力。因此，本文將針對運動頻率對液體晃蕩的影響進行模擬研究。一般情況下，液體晃蕩的強度主要通過波面幅值和作用在結(jié)構(gòu)上的壓力來反映。本文中，容器左側(cè)壁處的波面幅值和作用在液艙上的水平力將被用來衡量液體晃蕩強度。其中，水平力方向與x軸平行，且指向x軸正方向為正。為方便起見，本文統(tǒng)一規(guī)定：液艙尺寸L＝1.0 m，H＝0.5m，運動方程為θ=θ0sinωt，初始波面為靜止水面。通過數(shù)值收斂性驗證，時間步長和網(wǎng)格分別為△t=0.01 s和△x=△z= 0.05m。

對于L=1.0 m、H＝0.5 m的矩形容器，其固有頻率可由求得。其中，n表示固有頻率的階數(shù)。各階固有頻率分別為ω1＝5.32 rad/s、ω2＝7.84 rad/s、ω3＝9.61 rad/s、ω4＝11.10 rad/s、ω5＝12.41 rad/s和ω6＝13.60 rad/s。下面將針對運動頻率是固有頻率的情況進行模擬分析。選取了6個固有頻率作為運動頻率，分別為ω1、ω2、ω3、ω4、ω5和ω6。轉(zhuǎn)角幅值統(tǒng)一為1.0度。轉(zhuǎn)動中心固定在靜止水面上，即e＝0。

圖6給出了6個運動頻率下液艙左側(cè)壁處波面時間歷程圖。在圖6（a）、（c）、（e）中，波面幅值隨著時間逐漸增大；同時，運動頻率越大，波面幅值增長的速率越小，波面幅值與時間的線性關(guān)系愈不明顯。在圖6（b）、（d）、（f）中，波面幅值隨著時間沒有增大的趨勢，而是保持在一個穩(wěn)定的數(shù)值，亦即沒有發(fā)生共振現(xiàn)象。因此，相對于偶數(shù)階固有頻率來說，液體晃蕩對奇數(shù)階固有頻率更為敏感并且當運動頻率是奇數(shù)階固有頻率時，頻率越小，其對液體晃蕩的影響越大，共振現(xiàn)象越明顯。

圖6 不同運動頻率情況下的容器左側(cè)壁處波面時間歷程Fig.6 Time series of surface elevation at the leftwall at differentexcitation frequency

3.2 轉(zhuǎn)動中心對液體晃蕩的影響

選定三種有代表性的轉(zhuǎn)動中心位置e分別為-0.2 m、0、0.2 m。運動頻率ω和轉(zhuǎn)角幅值分別為ω1（容器的一階固有頻率為5.316 rad/s）和1.0度。

圖7 容器左側(cè)壁處的波面時間歷程 Fig.7 Time series of surface elevation at the leftwall

圖8 作用在容器上的水平力時間歷程Fig.8 Time series of horizontal force acted on the container

圖7給出了三種情況下的容器左側(cè)壁處波面時間歷程。由圖可知，當轉(zhuǎn)動中心在靜止水面（e＝0）時，容器左側(cè)壁處的波面幅值最小。當轉(zhuǎn)動中心偏離靜水面，轉(zhuǎn)動中心在靜止水面以上（e＝-0.2m）情況下的容器左側(cè)壁處的波面幅值要遠遠大于轉(zhuǎn)動中心在靜止水面下面（e＝0.2m）情況下的。同時，圖8給出了作用在容器上水平力的時間歷程圖。由圖可知，轉(zhuǎn)動中心在靜止水面（e＝0）時的水平力最小。當轉(zhuǎn)動中心距離靜止水面0.2m時，轉(zhuǎn)動中心在靜止水面以上（e＝-0.2m）情況下的水平力要明顯大于轉(zhuǎn)動中心在靜止水面下面（e＝0.2m）情況下的，這同圖7反映了同樣的實質(zhì)，即轉(zhuǎn)動中心在靜止水面以上更能激起容器中的液體晃蕩。

3.3 隔板對液體晃蕩的影響

對于在海上航行的液貨船，縱搖是影響其內(nèi)液體晃蕩強度的重要因素。通常，人們在船艙內(nèi)部加上不同形式的隔板來抑制液體的劇烈晃蕩，從而保證船體的穩(wěn)定性。本文中主要考慮垂直隔板對縱搖容器中液體晃蕩的影響?？紤]液艙具有一個長為0.3m、厚為0.01m的隔板。隔板是剛性的，固定在容器底部的左側(cè)。隔板中心到容器中心的距離定義為D。運動頻率ω為無隔板液艙的一階固有頻率ω1，轉(zhuǎn)角幅值為1.0度。隔板上的網(wǎng)格大小為△x=0.01m和△z=0.02m。圖9給出了隔板在不同位置D情況下左右側(cè)壁處的最大波高值。由圖可知，左側(cè)壁處的波高值要大于右側(cè)壁處的。兩側(cè)壁處的波高值都隨著D的增大而逐漸增大。當隔板位置由0移動到0.3m時，波高值只增大很小的幅度。而當隔板位置由0.3m增大到0.4m時，波高值卻急劇增大。圖10給出了隔板在不同位置D情況下作用在液艙上的最大水平力。由圖可知，最大水平力隨著D的增大而逐漸增大。當隔板位置由0移動到0.3m時，最大水平力增大的幅值非常小。當隔板位置從0.3m增大到0.4m時，最大水平力卻增幅明顯。這與圖9中波面值與D的關(guān)系是一致的。 因此，當隔板位置在區(qū)間（0，0.3L）內(nèi)時，其位置對液體晃蕩的影響不太明顯。當隔板位置大于0.3L時，液體晃蕩將對其位置將非常敏感。

圖9 不同D情況下的左右側(cè)壁處最大波面高度Fig.9 Maximum surface-elevation at leftand rightwalls vs. distance between the baffle and bottom center

圖10 不同D情況下的容器所受的最大水平力Fig.10 Maximum horizontal-loads on the container vs. distance between the baffle and bottom center

4 結(jié) 論

針對縱搖液艙中液體晃蕩問題，本文利用高階邊界元方法建立模擬帶有簡晃隔板縱搖容器液體晃蕩的完全非線性時域數(shù)學模型，通過坐標轉(zhuǎn)換，使得計算域僅僅控制在隨體坐標系內(nèi)。通過縱搖容器中波面變化和作用在容器上荷載變化的模擬，并與已發(fā)表試驗和數(shù)值結(jié)果的對比，驗證了本文模型的準確性。進一步開展數(shù)值模擬研究發(fā)現(xiàn)，相對于容器的偶數(shù)階固有頻率來說，液體晃蕩對奇數(shù)階固有頻率更為敏感，并且當運動頻率是奇數(shù)階固有頻率時，頻率越小，其對液體晃蕩的影響越大；同時，轉(zhuǎn)動中心在靜水面以上容器中的液體晃蕩效果更加顯著，而轉(zhuǎn)動中心在靜水面上液體晃動效果最弱。此外，本文還研究了安置垂直隔板液艙中的液體晃蕩。對于隔板長度一定的情況，隔板位置D在區(qū)間（0，0.3L）內(nèi)時，其位置對液體晃蕩的影響不太明顯；當D大于0.3L時，液體晃蕩將對其位置將非常敏感。

參 考 文 獻：

[1]Faltisen OM.A nonlinear theory of sloshing in rectangular containers[J].Journal of Ship Research,1974,18:224-241.

[2]Nakayama T,Washizu K.Nonlinear analysis of liquid motion in a container subjected to forced pitching oscillation[J].International Journal for NumericalMethods in Engineering,1980,15:1207-1220.

[3]Nakayama T,Washizu K.The boundary elementmethod applied to the analysis of two-dimensional nonlinear sloshing problems[J].International Journal for NumericalMethods in Engineering,1981,17:1631-1646.

[4]Chen B F,Chiang H W.Complete 2D and fully nonlinear analysis of ideal fluid in containers[J].Journal of Engineering Mechanics,1999,125(1):70-78.

[5]Chen B F,Chiang HW.Complete two-dimensional analysis of sea-wave-induced fully non-linear sloshing fluid in a rigid floating container[J].Ocean Engineering,2000,27:953-977.

[6]Chen B F,Nokes R.Time-independent finite difference analysis of fully non-linear and viscous fluid sloshing in a rectangular container[J].Journal of Computational Physics,2005,209:47-81.

[7]Chen B F,Asce M.Viscous fluid in container under coupled surge,heave and pitchmotions[J].Journal ofWaterway, Port,Coastal and Ocean Engineering,2005,131(5):239-256.

[8]Chen B F,Wu CH.Effects of excitation angle and coupled heave-surge-swaymotion on fluid sloshing in a three-dimensional container[J].Journal ofMarine Science and Technology,2011,16:22-50.

[9]Akyildiz H,Unal E.Experimental investigation of pressure distribution on a rectangular container due to the liquid sloshing[J].Ocean Engineering,2005,32:1503-1516.

[10]Akildiz H,Unal E.Sloshing in a three-dimensional rectangular container:numerical simulation and experimental validation[J].Ocean Engineering,2006,33:2135-2149.

[11]Akildiz H.A numerical study of the effects of the vertical baffle on liquid sloshing in two-dimensional rectangular container[J].Journal of Sound and Vibration,2012,332:41-52.

[12]Armenio V,Rocca M L.On the analysis of sloshing ofwater in rectangular containers:numerical study and experimental validation[J].Ocean Engineering,1996,23(8):705-739.

[13]Bae H S,Yun YW,Park M K.Flow response ofmagnetic fluid surface by pitchingmotion[J].Journal of Mechanical Science and Technology,2010,24(2):583-592.

[14]Celebi M S,Akyildiz H.Nonlinearmodeling of liquid sloshing in amoving rectangular container[J].Ocean Engineering, 2002,29:1527-1553.

[15]Xue M A,Lin P Z.Numerical study of ring baffle effects on reducing violent liquid sloshing[J].Computers and Fluids, 2011,52:116-129.

[16]肖龍飛,楊建民,郭 彬.淺水FPSO垂蕩和縱搖運動的低頻響應[J].艦船科學技術(shù),2009,31(11):120-125. Xiao L F,Yang JM,Guo B.Low frequency heave and pitch motions of FPSO in shallow water[J].Ship Science and Technology,2009,31(11):120-125.

[17]尹立中,鄒經(jīng)湘,王本利.俯仰運動圓柱貯箱中液體的非線性晃動[J].力學學報,2000,32(3):280-290. Yin L Z,Zou JX,Wang B L.Nonlinear sloshing of liquid in a circular cylindrical container under pitching excitation[J]. Chinese Journal of Theoretical and Applied Mechanics,2000,32(3):280-290.

[18]洪 亮,朱仁傳,繆國平,范 菊.液艙晃蕩對船舶橫搖運動影響的數(shù)值研究[J].中國造船,2015,1:1-10. Hong L,Zhu R Z,Miu G P,Fan J.Numerical research on effectof sloshing on ship rollmotion[J].Shipbuilding of China, 2015,1:1-10.

[19]徐國徽,胡嘉駿,顧學康.耦合運動下液貨艙晃蕩壓力預報研究[J].船舶力學,2016,20(7):884-891. Xu G H,Hu JJ,Gu X K.Prediction of sloshing pressure in coupledmotions for Liquid Cargo Carriers[J].Journal of Ship Mechanics,2016,20(7):884-891.

[20]李玉成,滕 斌.波浪對海上建筑物的作用(第2版)[M].北京:海洋出版社,2002. Li Y C,Teng B.Wave action onmarine structures(2nd Eidtion)[M].Beijing:Ocean Press,2002.

[21]寧德志.快速多極子邊界元方法在完全非線性水波問題中的應用[D].大連:大連理工大學,2005. Ning D Z.Application of the fastmultiple boundary elementmethod to the fully nonlinear water wave problems[D].Dalian: Ph.D thesis,Dalian University of Technology,2005.

[22]HiguchiM,Tanaka T,Endo S.Study on hull vibration-induced container liquid sloshing in LPG containerers[R].Nippon Kokan Technology Report,1976,72:111-122.

Nonlinear numerical simulation of liquid sloshing in a container subjected to pitch excitation

NING DE-zhi1,SONGWei-hua1,2,TENG Bin1
(1.State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China; 2.CCCC-FHDIEngineering Co.,Ltd.,Guangzhou 510230,China)

In the field ofmarine engineering,liquid sloshing is a kind of universal physical phenomenon. For the ship,the rotation has amore important influence than translation.Therefore,a time-domain numerical model is developed by higher-order boundary elementmethod to solve the liquid sloshing in a tank subjected to pitchmotion,in which the fully nonlinear boundary conditions satisfied on the free surface.With the coordinate exchanged between the global and local coordinates,the computational domain can be governed only in the local coordinates.In the solving process,a semi-mixed Eulerian-Lagrangian technique is applied to track the transient free surface and the 4th-order Runge-Kuttamethod is used to refresh wave elevation and velocity potential on the free surface at each time-step.The proposed numerical modelwas testified by comparison with the other published experimental and numerical results.On the base of validation,lots of numerical experiments are carried out to investigate the effect of pitch motion frequency,rotational center and a vertical bafflemounted on the container bottom on the sloshing surface and loads on the container.

liquid sloshing;pitchmotion;fully nonlinear;higher order boundary element;baffle

O353.2

：Adoi:10.3969/j.issn.1007-7294.2017.01.003

2016-09-01

國家自然科學基金資助項目（51679036，51490672）；教育部新世紀優(yōu)秀人才支持計劃（NCET-13-0076）

寧德志（1975-），男，教授，博士生導師，E-mail：dzning@dlut.edu.cn；

宋偉華（1987-），男，工程師，E-mail：songweihua9@163.com。

1007-7294（2017）01-0015-08