周 陽(yáng)，黃維平

（中國(guó)海洋大學(xué)山東省海洋工程重點(diǎn)實(shí)驗(yàn)室，山東 青島266100）

對(duì)渦激振動(dòng)的研究由來(lái)已久，主要針對(duì)物體在一定速度的流體中，其尾流產(chǎn)生交替的渦泄，導(dǎo)致物體兩側(cè)壓強(qiáng)的交替變化，從而引起結(jié)構(gòu)的振動(dòng)，對(duì)于Spar平臺(tái)等大直徑浮式結(jié)構(gòu)，尾流的渦旋脫落周期較長(zhǎng)，且其自身剛度遠(yuǎn)遠(yuǎn)大于系泊系統(tǒng)的剛度，發(fā)生渦激振動(dòng)時(shí)，其運(yùn)動(dòng)幅度以殼體的整體運(yùn)動(dòng)為主，這種振動(dòng)稱為渦激運(yùn)動(dòng)（Vortex－Induced Motion，VIM）。目前海洋工程方面，主要是針對(duì)立管和海底管線等大長(zhǎng)細(xì)比的柔性結(jié)構(gòu)進(jìn)行渦激振動(dòng)研究［1－2］。但是，對(duì)于Spar平臺(tái)這種大型新型海洋結(jié)構(gòu)物的渦激運(yùn)動(dòng)研究才剛起步［3－7］。國(guó)際上已有一些學(xué)者對(duì)渦激運(yùn)動(dòng)進(jìn)行了運(yùn)動(dòng)特性研究［8－10］，其中包括側(cè)重于工程應(yīng)用角度的對(duì)海洋平臺(tái)等結(jié)構(gòu)物進(jìn)行的渦激運(yùn)動(dòng)整體預(yù)報(bào)的相關(guān)研究［8－9］，還包括側(cè)重于考察非流線型柱狀主體周圍的復(fù)雜流動(dòng)現(xiàn)象、渦的生成和脫落特征方面的機(jī)理性研究［10］。

渦激運(yùn)動(dòng)包括橫流向和順流向的兩向渦激運(yùn)動(dòng)。大部分軟件對(duì)渦激運(yùn)動(dòng)進(jìn)行計(jì)算時(shí)，會(huì)忽略流固耦合的作用，即使考慮流固耦合，也只是考慮結(jié)構(gòu)位移對(duì)流場(chǎng)變化帶來(lái)的影響，而不考慮結(jié)構(gòu)速度帶來(lái)的影響，即非線性阻尼帶來(lái)的影響。為了更好的研究流固耦合對(duì)于Spar平臺(tái)兩向渦激運(yùn)動(dòng)的影響，本文運(yùn)用Matlab程序分別對(duì)Spar平臺(tái)順流向和橫流向渦激運(yùn)動(dòng)進(jìn)行了流固耦合數(shù)值模擬，探討結(jié)構(gòu)速度和兩向渦激運(yùn)動(dòng)之間的關(guān)系，并對(duì)結(jié)構(gòu)非線性阻尼的影響進(jìn)行了探討。

1 模型介紹

本文主要針對(duì)一硬艙外徑為28m，長(zhǎng)為70m，連接段外徑為14m，長(zhǎng)為84m，軟艙外徑為28m，長(zhǎng)為15m的簡(jiǎn)易Spar平臺(tái)模型進(jìn)行分析，其中干舷為16m（見(jiàn)圖1）。系泊系統(tǒng)由9根系泊纜組成，9根系泊纜分成3組對(duì)稱布置（見(jiàn)圖2），其水平剛度為93.745kN／m。進(jìn)行了流速為0.2～2m／s的情況下的兩向渦激運(yùn)動(dòng)的數(shù)值模擬。

圖1 Spar平臺(tái)模型Fig.1 Model of Spar platform

2 數(shù)值模擬方法

本文將Spar平臺(tái)看做剛體進(jìn)行整體分析，運(yùn)用Matlab軟件根據(jù)Newmark－β法編程完成平臺(tái)的兩向渦激運(yùn)動(dòng)響應(yīng)計(jì)算［11］。

圖2 Spar平臺(tái)系泊系統(tǒng)平面布置Fig.2 Floorplan of Spar mooring system

2.1 數(shù)模建立

由于數(shù)值模擬的需要，需要對(duì)模型進(jìn)行簡(jiǎn)化，并引入一些假定。

平臺(tái)所受渦激升力包括兩部分，一部分為硬艙和軟艙渦泄時(shí)產(chǎn)生的兩向渦激升力，由于直徑相同，因此升力頻率一樣，故對(duì)于平臺(tái)的受力可以疊加在一起進(jìn)行。另一部分為連接段渦泄時(shí)產(chǎn)生的兩向渦激升力，由于直徑與硬艙和軟艙不同，因此其渦泄升力頻率不一樣，其作用力要單獨(dú)加在平臺(tái)上。

由于研究平臺(tái)渦激運(yùn)動(dòng)時(shí)，只考慮水平面內(nèi)的運(yùn)動(dòng)，因此系泊系統(tǒng)的剛度只考慮其水平剛度，由于位移比較小，同時(shí)忽略了平臺(tái)產(chǎn)生位移時(shí)對(duì)系泊系統(tǒng)水平剛度的影響，將水平剛度視為一個(gè)定值。

由于Spar平臺(tái)自上而下不是均勻的圓筒形，其渦激升力沿高度方向也不統(tǒng)一，因此，實(shí)際過(guò)程中，平臺(tái)會(huì)發(fā)生一定的橫搖和縱搖，為了計(jì)算需要，本文同時(shí)忽略了平臺(tái)發(fā)生橫搖和縱搖時(shí)對(duì)此次數(shù)值模擬的影響。

2.2 Newmark－β法參數(shù)選定

Newmark－β法中，根據(jù)拉格朗日中值定理，t＋Δt時(shí)刻的速度表示成

于是t＋Δt時(shí)刻的速度

由泰勒級(jí)數(shù)導(dǎo)出

利用式（3）和（4），寫(xiě)出t＋Δt時(shí)刻的速度和加速度的表達(dá)式

把式（6）代入（5）并消去lui＋1和¨ui＋1得到關(guān)于ui＋1的求解方程

其中

求解方程（7）得到ui＋1，然后根據(jù)公式（6）可求得

根據(jù)上面建立的迭代公式，可歸納出Newmark－β法的計(jì)算過(guò)程：第一步形成K＊，第二步根據(jù)F＊i＋1求解方程（7）得到位移響應(yīng)ui＋1，并按照式（6）求速度和加速度響應(yīng)。逐步進(jìn)行迭代即可得到整個(gè)時(shí)間段的響應(yīng)情況，這是一種數(shù)值解的方法。在此次計(jì)算中，取α＝

3 流固耦合時(shí)兩向渦激運(yùn)動(dòng)的數(shù)值模擬

3.1 考慮非線性阻尼時(shí)順流向渦激運(yùn)動(dòng)響應(yīng)

平臺(tái)的順流向渦激運(yùn)動(dòng)是指渦旋脫落時(shí)對(duì)平臺(tái)引起沿著來(lái)流方向的曳力，導(dǎo)致平臺(tái)產(chǎn)生沿流方向的往復(fù)運(yùn)動(dòng)，相關(guān)圓柱形平臺(tái)的運(yùn)動(dòng)方程如下：

式中：CLI為脈動(dòng)拖曳力系數(shù)，與Re有關(guān)；u為來(lái)流速度；x為立管順流向振動(dòng)位移；lx為立管順流向振動(dòng)速度為立管順流向振動(dòng)加速度；wS為渦泄頻率；M為單位長(zhǎng)度質(zhì)量；K為單位長(zhǎng)度錨鏈水平剛度；C為結(jié)構(gòu)阻尼；Cm為附加質(zhì)量系數(shù)；CD為拖曳力系數(shù)。其中根據(jù)文獻(xiàn)［3，6］，選定CLI為0.05，CD為1.4，Cm為1。渦泄圓頻率其中St為斯托哈爾數(shù)，取0.25。

本文提及的非線性阻尼項(xiàng)指的是由于平臺(tái)和流體之間的相對(duì)運(yùn)動(dòng)而引起的拖曳力項(xiàng)，對(duì)于順流向的情況，非線性阻尼項(xiàng)為

由于Spar平臺(tái)中間段，硬艙和軟艙的直徑不同，計(jì)算時(shí)將中間段的渦激升力以及硬艙和軟艙的渦激升力根據(jù)式（9）沿長(zhǎng)度方向積分分別加載于平臺(tái)之上，附加質(zhì)量力也如此處理。而與平臺(tái)質(zhì)量相關(guān)的慣性力，總平臺(tái)結(jié)構(gòu)阻尼以及總錨鏈水平剛度則以實(shí)際Spar平臺(tái)的參數(shù)為準(zhǔn)進(jìn)行加載，其中Spar平臺(tái)總質(zhì)量為38 000t，總結(jié)構(gòu)阻尼為4.18×105kg·s－1，總錨鏈水平剛度為93.745kN／m。此文對(duì)于Spar平臺(tái)的渦激運(yùn)動(dòng)數(shù)值模擬均采用這種加載方式進(jìn)行。

根據(jù)Newmark－β法對(duì)平臺(tái)的運(yùn)動(dòng)方程進(jìn)行編程計(jì)算，分析Spar平臺(tái)考慮流固耦合時(shí)在流速0.2～2m／s工況中，由于渦泄引起的順流向振動(dòng)響應(yīng)。平臺(tái)的位移以幅值比的形式給出，即位移比硬艙直徑的形式，同時(shí)計(jì)算流速度以約化速度的形式給出（見(jiàn)圖3）。

圖3 順流向渦激運(yùn)動(dòng)幅值比隨約化速度的變化Fig.3 The relationship between in－line amplitude ratio of VIM and reduced velocity

從圖中可以看到，順流向的幅值比隨著約化速度的增加而增加，是一個(gè)遞增曲線，可見(jiàn)流速的大小對(duì)于平臺(tái)順流向的運(yùn)動(dòng)幅值起著關(guān)鍵的作用。但是順流向的運(yùn)動(dòng)響應(yīng)很小，即使在約化速度為15的情況下，其幅值比僅約為0.006，即位移為0.168m。所以目前對(duì)于Spar平臺(tái)的順流向響應(yīng)基本不予以考慮。

圖4是1m／s流速時(shí)順流向渦激運(yùn)動(dòng)位移響應(yīng)情況，可見(jiàn)Spar平臺(tái)在順流向產(chǎn)生了恒定的剛體位移23.78m，并在這個(gè)位置附近來(lái)回振蕩。

運(yùn)用Matlab對(duì)數(shù)據(jù)采用快速傅里葉變換（FFT）進(jìn)行譜分析，得位移周期見(jiàn)圖5，清晰可見(jiàn)與渦激升力一致的運(yùn)動(dòng)周期56s。通過(guò)圖5可知，連接段處渦旋脫落對(duì)平臺(tái)整體順流向的運(yùn)動(dòng)響應(yīng)影響比較微弱，因?yàn)橥ㄟ^(guò)計(jì)算可知連接段處渦旋脫落所引起的順流向振動(dòng)周期為28s，而這個(gè)數(shù)據(jù)在周期圖上并不是很明顯。

圖4 1m／s時(shí)平臺(tái)順流向位移圖Fig.4 Displacement time history in the in－line direction when the velocity is 1m／s

圖5 1m／s時(shí)的周期圖Fig.5 Periodogram when the velocity is 1m／s

3.2 考慮非線性阻尼時(shí)橫流向渦激運(yùn)動(dòng)響應(yīng)

平臺(tái)的橫流向渦激運(yùn)動(dòng)是指渦旋脫落時(shí)對(duì)平臺(tái)引起垂直于來(lái)流方向的升力，導(dǎo)致平臺(tái)產(chǎn)生垂直于來(lái)流方向的往復(fù)運(yùn)動(dòng)，相關(guān)圓柱形平臺(tái)的運(yùn)動(dòng)方程如下：

式中：CL為升力系數(shù)，與Re有關(guān)；y為立管橫流向振動(dòng)位移；ly為立管橫流向振動(dòng)速度為立管橫流向振動(dòng)即為非線性阻尼項(xiàng)。

運(yùn)用Matlab對(duì)平臺(tái)的運(yùn)動(dòng)方程進(jìn)行編程計(jì)算，分析Spar平臺(tái)考慮流固耦合時(shí)在流速0.2～2m／s工況中，由于渦泄引起的橫流向振動(dòng)響應(yīng)（見(jiàn)圖6）。

圖6 橫流向渦激運(yùn)動(dòng)幅值比隨約化速度的變化Fig.6 The relationship between crosss－flow amplitude ratio of VIM and reduced velocity

從圖中可以看到，隨著約化速度的增加，橫流向的幅值比先增加后減小，在約化速度為5的附近出現(xiàn)一個(gè)峰值，此時(shí)硬艙渦泄引起的橫流向升力頻率接近平臺(tái)的固有頻率，結(jié)構(gòu)發(fā)生共振。從圖形可見(jiàn)，這與順流向的運(yùn)動(dòng)響應(yīng)有一定的不同，制約橫流向運(yùn)動(dòng)響應(yīng)幅值的不再是單一的流速，而是與結(jié)構(gòu)的固有頻率也有一定的關(guān)系，來(lái)流引起的升力頻率越遠(yuǎn)離固有頻率，響應(yīng)則越小。同時(shí)又可從圖中看到，橫流向的響應(yīng)幅值比明顯比順流向的大得多，在共振時(shí)，平臺(tái)的幅值比達(dá)到了0.11附近，即位移為3.08m，這使得橫流向渦激運(yùn)動(dòng)的研究成為Spar平臺(tái)兩向渦激運(yùn)動(dòng)的重點(diǎn)。

圖7是1m／s流速時(shí)橫流向渦激運(yùn)動(dòng)位移響應(yīng)情況，可見(jiàn)Spar平臺(tái)在位移為0處來(lái)回振蕩，做簡(jiǎn)諧運(yùn)動(dòng)，幅值接近1m。

圖7 1m／s時(shí)平臺(tái)橫流向位移圖Fig.7 Displacement time history in the cross－flow direction when the velocity is 1m／s

運(yùn)用Matlab對(duì)數(shù)據(jù)采用快速傅里葉變換（FFT）進(jìn)行譜分析，得位移周期圖8，清晰可見(jiàn)與渦激升力一致的運(yùn)動(dòng)周期56，112s和平臺(tái)的固有周期200s。其中周期56s是連接段渦激升力產(chǎn)生的，周期112s是硬艙和軟艙渦激升力產(chǎn)生的。從圖8可見(jiàn)，渦激升力周期與平臺(tái)固有周期比較接近，隨著流速的降低，平臺(tái)將會(huì)和渦激產(chǎn)生的升力產(chǎn)生共振，這和圖6顯示的結(jié)果一致。

圖8 1m／s時(shí)的周期圖Fig.8 Periodogram when the velocity is 1m／s

3.3 非線性阻尼對(duì)順流向渦激運(yùn)動(dòng)響應(yīng)的影響

為了進(jìn)一步了解流固耦合中非線性阻尼力對(duì)結(jié)構(gòu)運(yùn)動(dòng)的影響，不考慮非線性阻尼力時(shí)的結(jié)構(gòu)運(yùn)動(dòng)模型同樣通過(guò)Matlab建立起來(lái)，這里先討論非線性阻尼力對(duì)順流向渦激運(yùn)動(dòng)響應(yīng)的影響。

Spar平臺(tái)在海流作用下產(chǎn)生渦激運(yùn)動(dòng)，順流向的運(yùn)動(dòng)較橫流向的運(yùn)動(dòng)幅度小得多。通過(guò)Newmark－β法同樣建立不考慮非線性阻尼時(shí)的順流向模型，方程如下，可與公式（9）進(jìn)行對(duì)比：

運(yùn)用方程（11）計(jì)算得平臺(tái)此時(shí)的運(yùn)動(dòng)響應(yīng)，同樣平臺(tái)的位移以幅值比的形式給出。為了方便比較，同時(shí)將考慮非線性阻尼時(shí)，平臺(tái)順流向的響應(yīng)情況在圖中給出（見(jiàn)圖9）。

圖9 順流向渦激運(yùn)動(dòng)幅值比隨約化速度的變化Fig.9 The relationship between in－line amplitude ratio of VIM and reduced velocity

通過(guò)上圖比較發(fā)現(xiàn)平臺(tái)考慮非線性阻尼時(shí)的順流向響應(yīng)比不考慮時(shí)來(lái)得小，而且響應(yīng)均隨著約化速度的增加而增加，在約化速度為2.5附近都有一定的突變?cè)黾?，這是由于在這個(gè)約化速度時(shí)，平臺(tái)上部硬艙的順流向渦激升力頻率與平臺(tái)固有頻率接近。然而在共振時(shí)，顯然非線性阻尼的作用更為明顯，它能有效的降低平臺(tái)順流向共振時(shí)位移響應(yīng)的計(jì)算結(jié)果，這對(duì)于Spar平臺(tái)將來(lái)可能的順流向位移計(jì)算有一定的參考價(jià)值。

3.4 非線性阻尼對(duì)橫流向渦激運(yùn)動(dòng)響應(yīng)的影響

Spar平臺(tái)渦激運(yùn)動(dòng)時(shí)橫流向的響應(yīng)是目前研究熱點(diǎn)，是Spar平臺(tái)設(shè)計(jì)的重要指標(biāo)。平臺(tái)渦激運(yùn)動(dòng)時(shí)橫流向運(yùn)動(dòng)的幅值較順流向運(yùn)動(dòng)的幅值大得多，對(duì)結(jié)構(gòu)和系泊系統(tǒng)的影響更為顯著。對(duì)于Spar平臺(tái)渦激運(yùn)動(dòng)的橫流向響應(yīng)必須控制在一定范圍內(nèi)。

通過(guò)Matlab軟件對(duì)Spar平臺(tái)不考慮非線性阻尼時(shí)的情況進(jìn)行編程，方程如（12），可與公式（10）進(jìn)行比較。

將不考慮非線性阻尼時(shí)的橫流向響應(yīng)數(shù)據(jù)同樣以幅值比的形式給出，同時(shí)引入考慮非線性阻尼時(shí)的數(shù)值模擬結(jié)果進(jìn)行比較（見(jiàn)圖10）。

圖10 橫流向渦激運(yùn)動(dòng)幅值比隨約化速度的變化Fig.10 The relationship between crosss－flow amplitude ratio of VIM and reduced velocity

通過(guò)上圖比較發(fā)現(xiàn)平臺(tái)考慮非線性阻尼時(shí)的橫流向響應(yīng)比不考慮時(shí)的響應(yīng)來(lái)得小，這與順流向的比較情況是一致的，然而橫流向響應(yīng)隨著約化速度的增加先增加后減小，在約化速度為5附近出現(xiàn)峰值，這是由于在這個(gè)約化速度時(shí)，平臺(tái)上部硬艙的橫流向渦激升力頻率與平臺(tái)固有頻率接近。在共振時(shí)，顯然非線性阻尼的作用也是更為明顯，它能有效的降低平臺(tái)橫流向共振時(shí)位移響應(yīng)的計(jì)算結(jié)果。在峰值兩邊，兩條曲線基本重合，可見(jiàn)，當(dāng)平臺(tái)不發(fā)生共振時(shí)，非線性阻尼的作用對(duì)于平臺(tái)的橫流向運(yùn)動(dòng)響應(yīng)的計(jì)算結(jié)果影響很小，這對(duì)于Spar平臺(tái)橫流向位移計(jì)算有一定的參考價(jià)值。

3.5 模型數(shù)值驗(yàn)證

為了驗(yàn)證本文數(shù)值模擬結(jié)果的正確性，本文將考慮非線性阻尼的橫流向數(shù)值模擬結(jié)果與文獻(xiàn)［6］進(jìn)行了對(duì)比（見(jiàn)圖11）。

圖11 數(shù)值模擬對(duì)比圖Fig.11 Comparison diagram of numerical simulation

從圖中可以看出，本文計(jì)算結(jié)果和文獻(xiàn)［6］有著一定的吻合性，在約化速度4～6區(qū)間內(nèi)，幅值比均有明顯上升趨勢(shì)，由于本文考慮了非線性阻尼的影響，而且本身Spar平臺(tái)的結(jié)構(gòu)阻尼較文獻(xiàn)［6］中的結(jié)構(gòu)阻尼大，因此在鎖定區(qū)本文的橫流向幅值比比文獻(xiàn)［6］中的小是合理的。在約化速度小于4或者大于6時(shí)，本文的計(jì)算數(shù)據(jù)與文獻(xiàn)［6］有著明顯的一致性，這很好的驗(yàn)證了本文理論數(shù)值計(jì)算的正確性。

4 結(jié)論

通過(guò)對(duì)Spar平臺(tái)兩向渦激運(yùn)動(dòng)的數(shù)值模擬，得到以下結(jié)果：

（1）得到了Spar平臺(tái)考慮流固耦合時(shí)兩向渦激運(yùn)動(dòng)位移響應(yīng)情況，并通過(guò)頻譜分析得到了相應(yīng)的周期圖。其中橫流向周期圖能較好地反映出硬艙和連接段2個(gè)不同直徑處的渦泄周期，同時(shí)還能從圖中觀察到平臺(tái)的固有周期。

（2）順流向渦激運(yùn)動(dòng)響應(yīng)隨著約化速度的增大而增大。橫流向渦激運(yùn)動(dòng)響應(yīng)隨著約化速度的增大先增大后減小，圖形曲線與順流向的完全不同，在約化速度5附近出現(xiàn)峰值，此時(shí)發(fā)生共振，由于其幅值較順流向的大得多，因此橫流向運(yùn)動(dòng)成為Spar平臺(tái)渦激運(yùn)動(dòng)考慮的重點(diǎn)。

（3）計(jì)算非線性阻尼對(duì)Spar平臺(tái)兩向渦激運(yùn)動(dòng)影響時(shí)發(fā)現(xiàn)，非線性阻尼能降低Spar平臺(tái)兩向渦激運(yùn)動(dòng)的計(jì)算結(jié)果，并且在共振時(shí)，減小的幅度更大?？梢?jiàn)在Spar平臺(tái)兩向渦激運(yùn)動(dòng)的計(jì)算中，只有在遠(yuǎn)離共振的情況下才可以忽略非線性阻尼。

（4）為了方便計(jì)算，本文忽略了渦泄引起的平臺(tái)的橫搖和縱搖問(wèn)題，有待日后繼續(xù)深入研究。

［1］ 葛斐，龍旭，王雷，等.大長(zhǎng)細(xì)比圓柱體順流向與橫向耦合渦激振動(dòng)的研究 ［J］.中國(guó)科學(xué)，2009，39（5）：752－759.

［2］ 張立武，陳偉民.深水細(xì)長(zhǎng)柔性立管渦激振動(dòng)響應(yīng)形式判定參數(shù)研究 ［J］.中國(guó)海上油氣，2010，22（3）：202－206.

［3］ 王穎，楊建明，李欣，等.均勻來(lái)流中浮式圓柱的渦激運(yùn)動(dòng)研究［J］.中國(guó)海洋平臺(tái)，2010，25（1）：8－15.

［4］ 王穎，楊建民，楊晨俊.Spar平臺(tái)渦激運(yùn)動(dòng)關(guān)鍵特性研究進(jìn)展［J］.中國(guó)海洋平臺(tái)，2008，23（3）：1－10.

［5］ 高云，宗智，于馨.Spar平臺(tái)渦激運(yùn)動(dòng)響應(yīng)分析 ［J］.中國(guó)海洋平臺(tái)，2011，26（1）：17－22.

［6］ 高云，宗智，于馨.均勻來(lái)流中Spar平臺(tái)渦激運(yùn)動(dòng)響應(yīng)研究 ［J］.中國(guó)科學(xué)，2011，41（2）：132－139.

［7］ 張蕙，楊建民，肖龍飛，等.均勻流中深水系泊Truss Spar平臺(tái)渦激運(yùn)動(dòng)試驗(yàn)研究 ［J］.海洋工程，2011，29（4）：14－20.

［8］ Sampath Atluri，John Halkyard，Senu Sirnivas.CFD simulation of truss spar vortex－induced motion ［C］.∥ Hamburg，Germany：The 25th International Conference on Offshore Mechanics and Arctic Engineering，2006.

［9］ Dominique Roddier，Tim Finnigan，Stergios Liapis.Influence of the reynolds number on spar vortex induced motions（VIM）：multiple scale model test comparisons［C］.∥Honolulu，Hawaii：Proceedings of the ASME 28th International Conference on Ocean，Offshore and Arctic Engineering，2009.

［10］ Wang Ying，Yang Jianmin，Lv Haining.Computational fluid dynamics and experimental study of lock－in phenomenon in vortexinduced motions of a cell－truss spar［J］.Journal of Shanghai Jiaotong University（Science），2009，14（6）：757－762.

［11］ 孟韜.基于 MATLAB的Newmark－β動(dòng)力反應(yīng)數(shù)值分析法的精度穩(wěn)定性分析 ［J］.金色年華，2010，7（6）：6.