半潛式平臺系泊運動數學模型及其仿真視景構建

張磊， 賈鵬， 蘇峰， 王恒， 姚廣旭

（1.海洋石油工程股份有限公司，天津 300451；2.哈爾濱工程大學機電工程學院，哈爾濱 150001）

0 引言

隨著近海油氣的開采日漸減少，對深海油氣資源的開發已成必然的趨勢，海洋平臺相關技術成為海洋工程界的研究熱點[1]，半潛式采油平臺的應用越來越廣泛。因此使用數學方法對半潛式平臺錨泊運動進行研究，對于平臺的安全性評估、保證平臺正常運行與生產，以及對遠洋深海工程錨泊定位技術的發展具有重大意義。

現階段隨著越來越多的深海浮式生產系統的服役，對各種浮式生產系統的運動方面的仿真也越來越多。對于油氣田浮式生產系統的仿真，基本上分為對各種類型的半潛式平臺的錨泊運動的仿真以及對與其相連的錨泊鎖鏈的仿真。其中，由多點系泊錨泊鎖鏈組成的錨泊系統，在環境的作用下為平臺提供了巨大的回復力，使平臺固定在海面的指定范圍內[2]，半潛式平臺在海洋環境中，在環境風、浪、流載荷的作用以及錨泊系統的回復力作用下，能夠產生6自由度方向上的運動，基于以上原理，平臺和錨泊系統組成了一個具有復雜運動的耦合運動系統[3-6]。

目前，對于半潛式錨泊系統的仿真技術已經相當成熟，C.M.Leech使用繩子模擬，在理論和實際方面分別對懸鏈線方法做了分析，能夠分析出懸鏈線扭曲的狀態，也能夠對輸水軟管做出模擬；Joel Brown、Jean Claude、Latombel等基于繩子提出了一種計算方法，能夠模擬出繩子在運動狀態下打結的效果；Gunnar Teichelmann、Mike Schaub、Bernd Simeon等基于鐵路電纜線的研究，在穩定受力平衡和動態振蕩的效果下，研究了對懸鏈線仿真方法。由于半潛式平臺結構相對復雜且種類較多，現階段并沒有普遍的數學模型和經驗公式對其6自由度運動進行模擬，對其運動的模擬，更多的是采用有限元分析軟件，基于有限元模型對其頻率特性、運動特性等做出分析[7-9]。

為了模擬復雜海況中各種環境載荷對平臺錨泊系統所產生的影響，本文根據所給定的平臺尺寸構建合適的半潛式平臺動力學模型，然后針對所得到的數學模型，研究系統錨泊運動特性，在Unity3d中開發虛擬海洋環境，使用3dsmax建立半潛式平臺的三維模型，并將三維模型導入到Unity3d的虛擬場景中，實現對半潛式平臺仿真實景的構建。

1 半潛式平臺錨泊運動方程

在實際海洋油田環境中，半潛式平臺與眾多錨泊線連接組成耦合系統[10]。對半潛式平臺的運動分析需建立2個方程：一為半潛式平臺本身的6自由度運動方程，另一為錨泊線拉力方程。

1.1 半潛式平臺運動方程

建立半潛式平臺運動方程，將平臺考慮為剛體，以平臺自身重心為原點、以平臺側弦方向為X軸、以頭部方向為Y軸、以垂直于甲板方向為Z軸建立慣性的空間笛卡爾坐標系。在慣性空間坐標系中應用牛頓第二運動定律建立平臺運動方程式如下：

式中：Mb為平臺質量矩陣；Cb為平臺輻射阻尼矩陣；Kb為平臺剛度矩陣；U為平臺的位移矢量；為平臺的速度矢量為平臺的加速度矢量；Fb為平臺受到的環境載荷。

1.2 錨泊線方程

對于錨泊系統的分析，需分析自然狀態下的單根錨泊線對平臺的拉力，然后基于懸鏈線方程建立單根錨泊線數學模型，懸鏈線方程是一根固定于兩點間的自然懸垂狀態下的曲線。在笛卡爾坐標系中，懸鏈線方程為：

方程三式分別代表在笛卡爾坐標系中，懸鏈線方程縱坐標y、長度L、拉力T隨橫坐標x的變化，其中a=H/（ρg）；H為懸鏈線橫向拉力；ρ為懸鏈線在水中單位重量；θ為懸鏈線初始張角。

將懸鏈線初始張角θ近似取零，則可得到簡化后的懸鏈線錨泊方程[11]：

2 環境載荷

對于海洋上的受約束低速運動結構，其受到的環境載荷包括環境風、浪和海流的作用。對于海洋波浪力載荷的確定，一般基于線性波理論，確定波浪力數學模型。風載荷和流載荷作用對平臺微幅運動的效果影響并不明顯，在實際計算中，一般將風載荷和流載荷考慮為線性作用[12]。

2.1 波浪干擾力數學模型建立

海洋中的波浪具有高度的隨機性，海浪并沒有固定的參數、形狀等供人們研究。實際工程中，波浪一般被分為規則波和不規則波，規則波與不規則波又根據其對漂浮物的作用特性，被分為主體力與橫漂力。本文分別建立規則波和不規則波的主體力及橫漂力數學模型，來描述整體的波浪干擾力。

2.1.1 規則波浪干擾力數學模型

對于被視作剛體的海洋漂浮物，其受到的規則波主干擾力可被分解為作用在其6個自由度上對其結構體積的三重積分式：

式中：k為2π內波的個數，k=2π/λ；λ為波浪平均波長；a為波浪運動幅度；c為波浪的傳播速度，其與頻率有關系，c=ω/k；ψ為平臺艏搖角。

以上對規則波浪主干擾力的計算公式，是基于對半潛式平臺浸水體積的積分。在計算時，將需要計算的部分考慮為若干立方體的組合以簡化計算。

2.1.2 規則波浪橫漂力數學模型

在波浪力中，橫漂力是作用在水平面方向上的作用力。實際計算時，通過已知條件波浪橫漂力與波幅的平方成正比來簡化計算[13]。規則波浪橫漂力標準形式如下式：

式中：ρ為海水密度；L為半潛式平臺長度；a為平均波幅；CXwD、CYwD、CNwD分別為縱向、橫向以及波浪力矩飄移系數，χ為波向角。

其中，CXwD、CYwD、CNwD可通過下式計算：

2.1.3 不規則波主干擾力數學模型

對于不規則波的數值計算，一般采用統計學對其進行分析。從統計理論方面考慮，不規則波可以視為是對若干個參數不同的單位規則波的疊加，基于以上理論，海浪特征可用波能譜的概念進行描述。Jonswap波能譜函數為

式中：ωP為峰值頻率，rad/s，參數β和σ一般取：

γ為譜峰提升因子，一般取平均值3.3。α取下值：

其中Hs為有義波高，一般取平均波高的1.6倍。基于波能譜函數，不規則波主干擾力數學模型如下式所示：

對于不規則波主干擾力的數學模型，也是6自由度方向上對不同頻率的規則波在平臺浸水體積的三重積分的疊加[4]，為簡化計算，在波能譜函數中，僅取影響較為明顯的規則波疊加到數學模型中。

2.1.4 不規則波橫漂力數學模型

不規則波浪中波浪飄移力的計算，仍然考慮為對各種參數的規則波波浪飄移力的疊加[14]，不規則波二階波浪干擾力如下式：

2.2 風干擾力數學模型

對于半潛式平臺而言，受風部分主要是平臺上層建筑以及水線面以上的平臺浮樁，風力大小受到受風面積大小、受風點高度、受風面形狀以及風速的影響。半潛式平臺受到的簡化水平風力如下式：

式中：Cw為半潛式平臺風力形狀系數，對于平臺的表面結構取值為1；AT為平臺橫向受風面面積；AL為平臺縱向受風面面積；Lpp為平臺兩立柱間距離；ρw為空氣密度；VwR為平面上10 m處相對風速。

2.3 海流干擾力數學模型

海洋中海流流速一般隨海水深度的增大而減小，對于半潛式平臺，其整個水下部分基本處于淺水區，垂直方向流速變化不明顯，故將流速做常數值處理。對于半潛式平臺，海流力可通過如下公式計算：

式中：Cc為半潛式平臺流力系數，對于平臺的表面結構取值為1.5；Ap為平臺橫向浸水面面積；Aq為平臺縱向浸水面面積；Lpp為平臺兩立柱間距離；ρc為海水密度，VcR為平均相對流速。

3 半潛式平臺6自由度方程水動力系數分析

半潛式平臺6自由度方程水動力參數包括質量矩陣系數Mb、輻射阻尼矩陣系數Cb與剛度矩陣系數Kb，其中質量矩陣包含半潛式平臺主體結構質量與附加質量，附加質量在方程中代表與平臺接觸的水體速度變化時對平臺的作用力；輻射阻尼矩陣為與平臺接觸的水體對平臺的黏性阻力作用；剛度矩陣是由平臺結構及錨泊系統對平臺產生的回復力作用。

3.1 附加質量與輻射阻尼的計算

3.1.1 基于有限元分析的平臺水動力參數計算

使用AQWA對半潛式平臺的附加質量及輻射阻尼做出分析。對于海洋浮體的水動力分析，一般只考慮浮體的外殼受到的影響而忽略浮體內部結構。AQWA默認的撓射單元為處于水線面以下的部分，而水線面以上的部分是不需要做水動力分析的。如圖1所示為平臺主體的有限元模型。

經計算后，在AQWA中打開平臺有限元模型如圖2所示。

AQWA會生成平臺的附加質量矩陣以及輻射阻尼矩陣中的全部參數，其中數量較大的為矩陣右對角線的6個參數，矩陣中其余參數全部視作0。

3.1.2 半潛式平臺附加質量矩陣

圖1 平臺有限元模型

圖2 平臺AQWA模型

在平臺的二階運動微分方程中，Mb是由平臺主體質量矩陣和附加質量矩陣組成的，其中附加質量矩陣為6×6的矩陣，在AQWA計算得到矩陣右斜對角線上的參數如圖3所示。

所以平臺的質量矩陣為：

其中，M為平臺主體質量矩陣。

3.1.3 半潛式平臺輻射阻尼矩陣

同理，平臺的二階運動微分方程中的阻尼矩陣Cb與Mb有同樣的形式，在AQWA計算得矩陣右斜對角線上的參數如圖4所示。

所以平臺的阻尼矩陣為

3.2 半潛式平臺剛度計算

半潛式平臺的剛度受到錨泊系統以及平臺自身結構的影響，其中錨泊系統提供平臺在橫蕩、縱蕩及艏搖方向上的剛度，平臺自身結構提供平臺在橫搖、縱搖及垂蕩方向的剛度。

3.2.1 錨泊系統回復力計算

平臺的錨泊系統由11根錨鏈組成，如圖5所示。求解11根錨泊線在水平方向的和作用，通過分析平臺在水平面上受到的錨泊作用隨平臺水平面上位移的變化，求解錨泊系統在橫蕩、縱蕩和艏搖方向上的剛度系數。

對于錨鏈在水平方向上拉力的計算：根據其在水平面的俯視長度和錨泊線長度方程，解得參量a；根據橫向拉力公式H=ρga解得單根錨鏈橫向拉力；再通過向量加法計算合力及合力矩，如圖5所示。經計算，平臺在水平面上指定位置受到的錨泊作用如表1所示。

圖3 平臺附加質量系數

圖4 平臺阻尼系數

圖5 錨泊系統剛性計算

以上數據點對半潛式平臺在橫蕩、縱蕩和艏搖方向上受到錨泊作用進行線性擬合，得到圖6～圖8。

3.2.2 半潛式平臺自身回復力計算

在半潛式平臺的錨泊運動中，平臺自身回復力出現在橫搖、縱搖與升沉3個方向上，且主要受平臺結構的影響。以下分別對3個方向的回復力作出分析計算：

1）升沉方向回復力計算。半潛式平臺漂浮在海面上時，其在垂直方向上具有上浮和下沉的效果，平臺的上浮和下沉均會引起平臺排水體積的變化，也就會引起其在垂直方向上受力的變化。

圖6 平臺橫蕩方向上錨泊作用擬合圖

由于浮筒在水線面上的運動面積是不發生變化的，如圖9所示，所以對于平臺在升沉方向的剛度可做簡單的線性處理，如下所示：

Ka=ρgAσ。

其中Aσ為水線面的面積，式中的算術值代表平臺在升沉方向上具有單位位移時，其受到浮力的變化。

表1 平臺在橫蕩方向上受錨泊作用點數據

表2 平臺在縱蕩方向上受錨泊作用點數據

表3 平臺在艏搖方向上受錨泊作用點數據

2）橫搖與縱搖方向回復力計算。當平臺具有橫搖或縱搖運動時，其受到的重力和浮力的合作用會產生反方向的傾覆力矩，則平臺在橫搖和縱搖上的剛度，可用下式表示：

其中zb為平臺重心與浮心間的垂直距離，上式代表平臺具有單位角位移時，其在相應方向上傾覆力矩的變化，由于平臺的搖動角度相對較小，所以忽略θ角的影響。

4 半潛式平臺錨泊運動仿真視景構建

半潛式平臺錨泊運動的仿真視景，一方面是支撐三維視景系統運行的硬件部分；另一方面為視景軟件部分，包括虛擬海洋場景、三維模型以及視景的后臺程序。

圖7 平臺縱蕩方向上錨泊作用擬合圖

圖8 平臺艏搖方向上錨泊作用擬合圖

4.1 視景硬件系統

視景系統的硬件部分包括信號采集設備及視景顯示設備兩部分。信號采集設備用以探測海空的各類環境參數，并將其處理成計算機能夠識別的信息，并在視景機端通過數學模型解算，將解算的結果以視景的形式在視景顯示設備上顯示出來。

4.2 半潛式平臺視景構建

圖9 平臺自身回復力計算原理圖

在Unity3d中開發虛擬海洋環境，使用3dsmax建立半潛式平臺的三維模型，并將半潛式平臺的三維模型的.FBX格式導入到Unity3d的虛擬場景中，組成半潛式平臺的虛擬視景。

4.3 差分法解平臺運動微分方程組

在視景系統中，使用VS2012作為開發環境，在VS2012中，使用C#語言開發腳本，編輯差分法解算平臺的運動微分方程組，并將解算得到的平臺位移通過腳本實時賦予視景中的平臺，從而賦予平臺相應的運動效果。

4.3.1 差分法解微分方程組

半潛式平臺的運動微分方程是六元二階微分方程。求解二階微分方程最終目的是求解半潛式平臺6自由度方向上的實時位移，在二階微分方程中可看出，方程結果受到海況參數和環境參數影響。系統中，海況參數和環境參數通過通信接口實時獲取，在具備這些參數的情況下，半潛式平臺6個自由度方向上的位移也就是關于時間的函數。根據以上分析，將半潛式平臺運動微分方程簡化成下式：

其中，p（t）、q（t）、f（t）均為[a,b]上的連續函數，則根據二階微分方程解的條件，可得微分方程存在唯一解。

二階微分方程求解往往伴隨著初值問題，方程的初值［u（a）,u′（a）］代表半潛式平臺在某一方向上的位移與速度，如下式所示：

圖10 視景硬件系統

差分法求解微分方程是基于連續方程離散化的原理，首先將區間a,［］b離散為N等份，則有下式：

圖11 半潛式平臺錨泊運動仿真視景

用ui表示u在點xi處的取值，對連續可導函數u（t），由泰勒級數展開式可得：

則將式（1）、式（4）與式（5）聯合，省略局部截斷誤差，可得二階微分方程三點差分格式式：

上式可改寫為：

其中：ai=2/h-pi,bi=-4/h+2hqi,ci=2/h+pi,di=2hfi。

以上即為二階微分方程差分離散形式。

由數值微分公式u′（a）≈u1-u0/h以及初值條件u′（a）=β可得：

則根據式（7）可依次求解各差分點值。

本文需求解的二階微分方程如下：

則對于平臺的單一自由度，其標準差分形式如下式：

式中的各參量中，僅有Fb為隨時間變化，其余參量均與平臺水動力系數有關，半潛式平臺在一定運動頻率范圍內時，取平臺水動力系數為常數進行迭代計算。

4.3.2 差分法解微分方程組腳本編輯

在視景系統中，使用C#語言編輯差分法解算平臺的運動微分方程組。半潛式平臺的6自由度微分方程經處理后，其最后能夠簡化成6個平臺位移、平臺受風向角、計算風速、波浪周期、波浪平均波長、波浪運動頻率、波浪傳播速度、洋流速度、平均浪高以及時間有關的6個二階微分方程，通過通信接口或者系統數據中的環境參數和海況參數，并通過計算獲得波浪運動頻率，則平臺的二階運動微分方程化簡為6個位移隨時間變化的二階微分方程。

在系統中，定義6個自定義函數，使用差分法分別對平臺6個自由度方向的運動微分方程進行求解，并將解得的運動數據賦予平臺。在單個自定義函數中，賦予平臺運動的流程如圖12所示。

圖12 賦予平臺運動的流程

圖中的程序中，每隔0.5 s取得平臺6自由度中單一方向上的位移量，并賦予平臺驅動平臺產生運動。其中子函數實現了對平臺運動微分方程的解算，其流程如圖13所示。

分別在上文的6個自定義函數中編輯圖的程序，同時計算出平臺6自由度的實時位移，圖中設置迭代步長 為 0.5，也就是系統每間隔0.5 s輸出一個位置點，根據精度的要求，可縮小迭代步長。

圖13 平臺運動微分方程的解算流程

5 半潛式平臺錨泊運動仿真視景系統測試

使用Matlab對半潛式平臺運動微分方程組加以解算，得到平臺的位移隨時間變化的關系；并在視景系統中記錄平臺位移隨時間變化的關系，驗證視景系統性能，并驗證視景系統對平臺運動的準確性。

5.1 運動微分方程組解算

取平臺6自由度上的初始位移與初始速度均為0，使用表4的海況數據及表5的平臺結構尺寸，分析平臺在6自由度方向上的運動，如圖14所示。

使用表4海況數據對平臺6自由度方向上運動的解算結果如圖14（a）～（f）。

本文對于微分方程的求解，賦予方程的初始值為[0,0]，其物理意義代表平臺在單一自由度方向上的位移和速度均為0，在宏觀上，其相當于在三維空間的（0,0,0）位置，賦予平臺初始速度為0，讓其在海水作用、環境載荷和錨泊作用下做進一步的運動。對以上各圖中平臺的運動結果進行分析，在海水作用、環境載荷以及錨泊作用下，由于三維空間中的[0,0,0]位置并不是平臺受力平衡的位置，在平臺的橫蕩和縱蕩以及艏搖方向上，平臺會在風力、海流力以及二階波浪力的作用下產生偏移，偏移達到一定程度時，會因為錨泊系統的回復力作用，將平臺拉回反方向，如此反復直到平臺到達平衡位置，在運動圖像中，其反應為平臺的振蕩運動，在垂蕩、橫搖和縱搖方向上，由于平臺自身的回復力作用，其具備相同的振蕩效果。當各自由度方向振蕩結束后，在規則及不規則波浪干擾力的作用下，平臺會在6自由度方向上產生具有一定規則性的小周期運動，這部分運動即是平臺的主體運動。通過對以上各圖的分析，在海況較為不穩定的情況下，平臺在6自由度方向上的小周期振動也會比較強烈。

表4 計算用海況數據

表5 平臺主要尺寸表 m

圖14 平臺歷時運動效果圖

圖15 平臺運動驗證

為驗證半潛式平臺錨泊運動的準確性，需要將計算結果與實測數據進行比較。本文在計算數據與實測數據的比較上，忽略掉其運動的振蕩環節，取平臺運動規律性較強的平穩運動時段進行比較與分析。圖15（a）～（f）為對平臺運動計算數據與實測數據的比較。

通過以上的分析，應用數學模型計算的結果，與實際海洋中的大型漂浮物的錨泊運動規律是大體類似的。將平臺在6自由度方向上運動的計算結果同實測數據進行比較，可發現數學模型計算結果同實測數據是基本吻合的，其中若干數據點的運動數據脫離數學模型的控制，經分析是由于海洋海況不穩定性造成的，如果不考慮這些不穩定性因素的影響，使用本文中的平臺運動微分方程，是可以在很大程度上描述出平臺的運動規律的。并能夠應用于三維視景系統，描述平臺的6自由度運動。

5.2 視景系統性能測試

在視景系統軟件端，是采用差分法[15]，每隔0.5 s計算一個平臺位置點，從而實現對平臺位移的解算。其結果受差分間隔和視景機性能的影響。

圖16 平臺運動視景效果圖

表6 初始時刻平臺垂蕩方向運動數據 m

表7 系統運行20 min后平臺垂蕩方向運動數據 m

表8 系統運行40 min后平臺垂蕩方向運動數據 m

表9 系統運行60 min后平臺垂蕩方向運動數據 m

表10 系統運行80 min后平臺垂蕩方向運動數據 m

表11 系統運行100 min后平臺垂蕩方向運動數據 m

表12 系統運行120 min后平臺垂蕩方向運動數據 m

在視景系統接口接入表4中海況信息，從而在系統中實現以差分法對平臺運動微分方程組的解算。在系統中，平臺運動效果可通過視景輸出及數據輸出表現出來，如圖16所示。

表13 系統運行140 min后平臺垂蕩方向運動數據 m

表14 系統運行160 min后平臺垂蕩方向運動數據 m

表15 系統運行180 min后平臺垂蕩方向運動數據 m

如圖16所示，當三維視景系統中具有海況信息時，系統會根據海況數據信息，通過數學模型解算出平臺的6自由度方向上的運動規律。圖中為根據數學模型解算出的平臺運動規律在視景、圖像和數據方面的表現。本文中監控的海況波浪平均周期為8.7 s，取圖像的監控時長為60 s，監控時長可根據實際需要手動調節其大小。

由于數據量較大，本文取平臺垂蕩方向運動數據進行分析，并同第3節中對方程的計算結果進行比較，分析三維視景系統對平臺運動描述的準確性，如表6～表15所示，每隔20 min取一組一周期上的運動數據，并對系統進行連續3 h的觀測。

對表中平臺的運動信息進行分析，由于系統運行的遲滯性，當系統長時間運行后，系統對平臺運動數據的表述也會出現遲滯。在60 min以內，系統能夠對平添運動信息進行準確表述，并進一步表明了平臺運動微分組對平臺運動描述的準確性。系統運行超過60 min后，系統內計算數據與理論數據的偏差會越來越大。

6 結論

建立了六元的二階微分方程組，并分析了方程組中的質量矩陣、阻尼矩陣、剛度矩陣及環境載荷，并建立了半潛式平臺的運動仿真視景，驗證了方程組描述半潛式平臺運動的合理性。

［參考文獻］

[1] 孫濤,桂文彬,俞志剛.半潛式平臺定位系泊控制試驗系統設計與應用[J].船舶工程,2012,34(2)：84-86.

[2] 張維.深水平臺錨泊系統錨索鏈受力分析及視景仿真[D].武漢：武漢理工大學,2014.

[3] XIE S.Design of Corrosion of Marine Oil and Gas Fields Underwater Manifold Pry[J].Guangdong Chemical Industry,2016,43(7)：158-159.

[4] 韋紅術,張偉國,陳國明,et al.Position determining method and position determining device of Christmas tree of marine oiland-gas field producing platform：CN 103132991 A[P].2013-06-05.

[5] 王俊榮,謝彬.半潛式平臺水動力性能及運動響應研究綜述[J].中國造船,2009(a11)：255-261.

[6] 史琪琪,楊建民.半潛式平臺運動及系泊系統特性研究[J].海洋工程,2010,28(4)：1-8.

[7] KOO B J,KIM M H.Hydrodynamic interactions and relative motions of two floating platforms with mooring lines in sideby-side offloading operation[J].Applied Ocean Research,2005,27(6)：292-310.

[8] YANG M,TENG B,NING D,et al.Coupled dynamic analysis for wave interaction with a truss spar and its mooring line/riser system in time domain[J].Ocean Engineering,2012,39(1)：72-87.

[9] 朱忠顯,尹勇,申和龍.基于集中質量法的錨泊系統動力學仿真[J].中國航海,2015,38(4)：38-42.

[10]童波,楊建民,李欣.深水半潛式平臺系泊系統動力特性研究[J].海洋工程,2009,27(1)：1-7.

[11] 王冬姣.索-鏈-浮子/沉子組合錨泊線的靜力分析[J].中國海洋平臺,2007,22(6)：16-20.

[12] 羅強.981半潛式鉆井平臺錨泊系統研究[D].荊州：長江大學,2012.

[13]韓鑫,藍尉冰.二階波浪載荷計算方法探究——以深海船式平臺動力定位系統為例[J].欽州學院學報,2014,29(5)：1-3.

[14] 姜宗玉,崔錦,董剛,等.不規則波中動力定位半潛平臺慢漂運動研究[J].中國海洋平臺,2015,30(3)：89-94.

[15] 張文生.科學計算中的偏微分方程有限差分法[M].北京：高等教育出版社,2006.