高懷玉，張 峰，秦忠國(guó)

(1.河海大學(xué)力學(xué)與材料學(xué)院，江蘇 南京 210098;2.南京市水利規(guī)劃設(shè)計(jì)院有限責(zé)任公司工程分院，江蘇南京 210006)

隨著工程技術(shù)領(lǐng)域的不斷開拓和發(fā)展，對(duì)大型復(fù)雜結(jié)構(gòu)的分析要借助高階數(shù)值分析模型［1］和大規(guī)模數(shù)學(xué)計(jì)算，以保證數(shù)值解的精度［2］。在普通的計(jì)算設(shè)備上采用傳統(tǒng)的串行計(jì)算面臨著因計(jì)算資源不足而無法計(jì)算以及計(jì)算時(shí)間過長(zhǎng)、工作效率低的問題。使用高性能計(jì)算機(jī)進(jìn)行并行有限元分析是解決上述難題的一個(gè)方案［3-4］。自從美國(guó)國(guó)家航空航天局NASA的Noor等［4］在1975年發(fā)表第1篇有限元并行計(jì)算的論文以來，有限元并行算法的研究越來越具體和廣泛，適用于不同類型計(jì)算機(jī)的有限元并行算法也相繼產(chǎn)生。20世紀(jì)90年代以來，單處理機(jī)的計(jì)算速度不斷提高，并行計(jì)算機(jī)的體系結(jié)構(gòu)趨于成熟，數(shù)據(jù)傳輸網(wǎng)絡(luò)的標(biāo)準(zhǔn)化和傳輸速率的大幅度提升，使得并行計(jì)算機(jī)的研制周期從幾年縮短至幾個(gè)月，為研制并行計(jì)算機(jī)系統(tǒng)創(chuàng)造了有利條件［5-7］。高性能計(jì)算機(jī)的發(fā)展勢(shì)頭異常迅猛，擁有并行計(jì)算功能的有限元軟件［8］也相繼出現(xiàn)。李根國(guó)等［9］將NASTRAN在“神威I”超級(jí)計(jì)算機(jī)系統(tǒng)上進(jìn)行移植和并行計(jì)算功能的二次開發(fā)，獲得了成功。一些優(yōu)秀的并行數(shù)值算法包也日益引起關(guān)注。借助并行算法包開發(fā)有限元并行程序無需考慮復(fù)雜的數(shù)據(jù)分布和通信，大大降低了有限元并行計(jì)算的難度和成本，縮短了開發(fā)周期［10］。

可移植可擴(kuò)展科學(xué)計(jì)算工具箱PETSc(portable，extensible toolkit for scientific computation)由一套數(shù)據(jù)結(jié)構(gòu)和程序組成，由美國(guó)Argonne國(guó)家實(shí)驗(yàn)室承擔(dān)開發(fā)并獲得成功，主要用于分布式存儲(chǔ)環(huán)境下高效求解偏微分方程組及相關(guān)問題的求解［3］。PETSc在軟件實(shí)現(xiàn)上采用ANSI C和消息傳遞標(biāo)準(zhǔn)接口MPI［11-12］。針對(duì)單處理機(jī)多核計(jì)算機(jī)，筆者利用PETSc提供的高層應(yīng)用程序開發(fā)平臺(tái)開發(fā)出一種線性并行求解方法，并應(yīng)用編寫的相應(yīng)求解程序?qū)η蠼馄鬟M(jìn)行了測(cè)試。

1 PETSc的數(shù)據(jù)結(jié)構(gòu)

PETSc的方程組求解器組件定義在各種數(shù)據(jù)結(jié)構(gòu)組件之上，包括線性方程組求解器KSP、非線性方程組求解器SNES、時(shí)間依賴方程求解器TS。PETSc基于基礎(chǔ)線性代數(shù)子程序庫BLAS、線性代數(shù)包LAPACK、MPI等實(shí)現(xiàn)計(jì)算功能。PETSc為用戶提供了豐富的Krylov子空間迭代方法和預(yù)條件子［13-14］，并具有錯(cuò)誤檢測(cè)、性能統(tǒng)計(jì)和圖形打印等功能?；赑ETSc提供的大量對(duì)象和求解庫，用戶可以靈活地開發(fā)自己的應(yīng)用程序，還可隨意添加和完善某些功能。

PETSc的數(shù)據(jù)結(jié)構(gòu)包括向量Vec、矩陣Mat、網(wǎng)格管理對(duì)象、索引與排序、結(jié)構(gòu)網(wǎng)格DA以及無結(jié)構(gòu)網(wǎng)格的IS對(duì)象。用戶在PETSc環(huán)境下編寫串行或并行應(yīng)用程序時(shí)，要按照固定的框架結(jié)構(gòu)來編寫，即對(duì)于一般的應(yīng)用對(duì)象均需要聲明、創(chuàng)建、設(shè)置、使用和銷毀。

2 基于PETSc的有限元求解器概況

基于PETSc的有限元并行求解器(以下簡(jiǎn)稱求解器)以C++類［15］的形式存在，求解器的類封裝了一般線性有限元分析所必需的數(shù)據(jù)及操作，用戶按照給定的格式提供數(shù)據(jù)，依次調(diào)用成員函數(shù)，即可完成一個(gè)有限元的并行分析，得到位移解。使用求解器時(shí)，可以通過改變變量nfreedom值來改變有限元分析的類別，如nfreedom值設(shè)為1時(shí)，對(duì)應(yīng)的是一維有限元問題，nfreedom值為2和3則分別對(duì)應(yīng)平面和空間有限元問題。自定義的單元結(jié)構(gòu)體StructElement，對(duì)外提供一個(gè)單元?jiǎng)偠染仃囉?jì)算程序接口，用戶可以自由地將自定義單元接入程序。求解器的類能方便地?cái)U(kuò)展單元類，具有一定的通用性和可擴(kuò)展性，它繼承了PETSc工具箱所提供的運(yùn)行時(shí)選項(xiàng)(runtime option)、性能分析等功能，即在求解時(shí)可以改變求解程序的迭代方法、預(yù)條件方法等設(shè)置，這些功能為測(cè)試程序的效率提供了方便。

2.1 求解器結(jié)構(gòu)

求解器在耗時(shí)較多的分析部分進(jìn)行并行處理，具體流程如圖1所示。當(dāng)程序運(yùn)行到插入節(jié)點(diǎn)力、形成荷載列陣和引入邊界條件、修改方程組模塊式時(shí)，進(jìn)程0要向后層其他各進(jìn)程廣播這些信息。當(dāng)程序運(yùn)行到計(jì)算矩陣預(yù)分配參數(shù)模塊時(shí)，各進(jìn)程之間要保持相互通信。最后，各進(jìn)程把各自計(jì)算出的結(jié)果返回到進(jìn)程0，并保存結(jié)果，結(jié)束計(jì)算。并行分析部分包括預(yù)分配、形成總剛度矩陣及求解方程組，這些部分消耗的時(shí)間一般情況下占整個(gè)過程的90%以上，而且這一比例會(huì)隨著計(jì)算規(guī)模的增大而增大，用如上并行方式求解線彈性問題，理論上能取得很好的效果。

2.2 求解器調(diào)用

求解器提供的可選設(shè)置選項(xiàng)說明如下:

a.SetElement(intType，intmaxnode)與 SetElement(int Type，proc－ESElementStiffer)是提供接口給用戶接入自己的單元類型。其中參數(shù)Type是單元類型號(hào)，求解器默認(rèn)類型為八節(jié)點(diǎn)六面體空間等參元;參數(shù)maxnode為該類型單元的節(jié)點(diǎn)數(shù);自定義指針函數(shù)類型 proc－ElementStiffer則為指向該單元類型的單剛求解程序。通常這2條設(shè)置應(yīng)該設(shè)置在初始化函數(shù)SolverInitialize()之后。

b.SetKSPType(char*ktype)與 SetPCType(char*ptype)語句用于設(shè)置子空間迭代法的類型，和PETSc工具箱的命名一樣。子空間迭代法的類型 ktype可以是KSPCG，KSPGMRES，KSPTCQMR，KSPBCGS，KSPTFQMR，KSPCR，KSPLSOR，KSPBICG中的任何一種。預(yù)條件法的類型 ktype包括 PCJACOBI，PCBJACOBI，PCSOR，PCEISENSTAT。

圖1 求解器流程Fig.1 Flow chart of solver

c.SetTolerances(double abstol，double rtol，double dtol，int maxits)用于設(shè)置迭代收斂參數(shù)，前3個(gè)參數(shù)是基于－范數(shù) ( (∑)1/2)的剩余范數(shù)相對(duì)下降量rtol，剩余范數(shù)的絕對(duì)值atol及剩余中的相對(duì)增加量dtol。Maxits則為最大迭代次數(shù)。

d.SetTimeInfo()函數(shù)用于保存求解時(shí)間明細(xì)，包括各主要部分的時(shí)間消耗，以此來分析求解效率。SaveBinary()與SaveInfo()是將剛度矩陣和等效節(jié)點(diǎn)力列陣分別保存為文本形式和二進(jìn)制形式，但只在單進(jìn)程下有效。其中二進(jìn)制保存的格式為PETSc規(guī)定的格式。使用求解器的類，即為依次調(diào)用類中的成員函數(shù)來實(shí)現(xiàn)有限元分析的過程。求解器的一個(gè)完整的使用示例如下:

SOLVER S;聲明一個(gè)求解器對(duì)象

S.Solver Initialize(argc，args);求解器初始化

if(S.rank=0)

S.Read Information(Char*path);由根進(jìn)程讀入數(shù)據(jù)

Set Element(int Type，int maxnode);設(shè)置單元節(jié)點(diǎn)數(shù)，Type為單元類型編號(hào)

Set Element(int Type，proc－ESElement Stiffer);設(shè)置單剛求解程序

S.Preallocation();創(chuàng)建矩陣、向量結(jié)構(gòu)、矩陣預(yù)分配

S.Add KERE();計(jì)算單剛、組集總剛和節(jié)點(diǎn)力列陣

Set KSP Type(char*ktype);設(shè)置求解迭代方法

Set PC Type(char*ptype);設(shè)置預(yù)條件方法

S.Set Bound Condition();根據(jù)邊界條件修正支配方程

Set Tolerances(double abstol，double rtol，double dtol，int maxits);設(shè)置收斂誤差

Save Binary();保存系數(shù)矩陣和右端頂，僅單進(jìn)程運(yùn)行時(shí)有效

腎功能指標(biāo)水平：取上述血清3 ml，用肌氨酸氧化酶檢測(cè)血肌酐(Scr)法酶偶聯(lián)速率法檢測(cè)血尿素氮(BUN)；采用核素動(dòng)態(tài)顯像法測(cè)定腎小球過濾率(GFR)。

Save ASCII();文本保存系數(shù)矩陣和右端項(xiàng)，僅單進(jìn)程運(yùn)行時(shí)有效

S.SOLVE();求解方程組

S.Get Answer();保存結(jié)果

Set Time Info();保存時(shí)間耗用明細(xì)

S.Finalize();結(jié)束

3 求解器的輸入輸出

求解器的輸入目前提供了2種形式:讀入文件和傳入網(wǎng)格信息結(jié)構(gòu)體指針。用戶在調(diào)用Solver Initialize()函數(shù)進(jìn)行初始化以后，即可在進(jìn)程0里接受該結(jié)構(gòu)體指針或執(zhí)行讀入文件函數(shù)。如以讀入文件的方式輸入，需要讀入節(jié)點(diǎn)信息、單元信息、材料信息、約束信息、節(jié)點(diǎn)力信息。

求解器的輸出包括3個(gè)部分:(a)方程構(gòu)成。求解器能夠在一個(gè)進(jìn)程運(yùn)行時(shí)輸出方程的系數(shù)矩陣和右端項(xiàng)，默認(rèn)的狀態(tài)是不輸出，用戶需要在求解方程之前將其設(shè)置成輸出。輸出文件的格式有文本文檔和二進(jìn)制格式。文本文檔可供查閱，保存時(shí)間較長(zhǎng)，存儲(chǔ)較大;二進(jìn)制格式存儲(chǔ)小，速度快，可由PETSc程序讀入。(b)求解信息。求解網(wǎng)格信息的簡(jiǎn)要說明，求解設(shè)置，以及可選的時(shí)間信息，供性能分析使用。(c)方程的解。包括方程的解和節(jié)點(diǎn)的位移矢量大小，在文件的開始則為位移的最大值及其節(jié)點(diǎn)號(hào)。

4 算例分析

筆者以計(jì)算某軟土地基上水閘在自重荷載下所產(chǎn)生的位移為算例，對(duì)本文所開發(fā)的并行求解器進(jìn)行測(cè)試。水閘結(jié)構(gòu)由3部分組成:閘身、樁基、軟土地基，材料參數(shù)見表1。本文僅作線彈性分析，故將地基和樁基密度置為零，閘身密度為2500kg/m3，以計(jì)算由水閘自身重力產(chǎn)生的位移。因水閘結(jié)構(gòu)對(duì)稱分布，因此取其一半結(jié)構(gòu)進(jìn)行建模［16］計(jì)算(圖2)。從求解器求解規(guī)模、并行效率、對(duì)于病態(tài)方程組的迭代收斂情況來測(cè)試本文開發(fā)的并行求解器。求解環(huán)境為一臺(tái)四核I5、主頻為2.8 GHz的32位臺(tái)式機(jī)，可使用內(nèi)存為3 G。

表1 水閘材料參數(shù)Table 1 Material parameters of water lock

4.1 求解規(guī)模

按網(wǎng)格疏密將水閘結(jié)構(gòu)離散為八節(jié)點(diǎn)六面體單元網(wǎng)格，表2為6組有限元模型的網(wǎng)格信息，圖2為網(wǎng)格單元數(shù)為51112的水閘網(wǎng)格圖。對(duì)每一組有限元模型都進(jìn)行四進(jìn)程并行求解，計(jì)算的時(shí)間如表2所示。

第6組模型在計(jì)算時(shí)由于每個(gè)進(jìn)程的內(nèi)存使用都超過600 M，系統(tǒng)需要使用虛擬內(nèi)存，導(dǎo)致計(jì)算時(shí)間達(dá)到2 409.33 s。在同樣的計(jì)算環(huán)境下，單進(jìn)程能計(jì)算節(jié)點(diǎn)數(shù)不超過40萬的網(wǎng)格，Ansys軟件能求解節(jié)點(diǎn)數(shù)不超過20萬的線彈性問題，可見本文求解器能在普通的多核計(jì)算機(jī)下處理更大規(guī)模的網(wǎng)格。

取節(jié)點(diǎn)數(shù)為58 622的網(wǎng)格模型比較本文求解器和Ansys軟件計(jì)算得到的節(jié)點(diǎn)位移，兩者的x，y，z方向最大位移及總位移間的誤差分別為0，0，0.007%，0.001%，最大位移值的誤差不超過1/10000(表3)。將節(jié)點(diǎn)位移導(dǎo)入軟件Tecplot中，繪制位移云圖。本文求解器得到的位移云圖與Ansys結(jié)果分布規(guī)律基本一致。求解器在提高收斂精度下能得到更精確的解，但因迭代次數(shù)增加而消耗更多的時(shí)間。

圖2 水閘結(jié)構(gòu)及整體網(wǎng)格Fig.2 Structure of water locks and overall grid

表2 不同規(guī)模的網(wǎng)格信息Table 2 Grid information with different scales

表3 本文求解器和Ansys軟件計(jì)算的部分節(jié)點(diǎn)的結(jié)果對(duì)比Table 3 Comparison of node calculation results obtained by solver developed in this study and ANSYS software

不同規(guī)模網(wǎng)格計(jì)算得出的最大節(jié)點(diǎn)位移值如圖3所示，可看出隨著單元總數(shù)增加，節(jié)點(diǎn)最大位移值也逐漸增大，此現(xiàn)象符合有限單元法中位移近似解是從下方逼近精確解的理論，因而水閘結(jié)構(gòu)的最大位移應(yīng)該在0.067 m附近。

4.2 并行加速效果

為驗(yàn)證并行加速效果，選取節(jié)點(diǎn)總數(shù)為178 007的網(wǎng)格數(shù)據(jù)，分別進(jìn)行1～4個(gè)進(jìn)程的并行計(jì)算。迭代方法為共軛梯度法(CG)［17］，預(yù)條件子為塊雅可比方法(bjacobi)。圖4為程序在單機(jī)多進(jìn)程并行求解時(shí)的并行加速比、并行效率與進(jìn)程數(shù)相關(guān)關(guān)系曲線。

程序的并行分析部分主要包括矩陣預(yù)分配參數(shù)計(jì)算、形成總剛度矩陣和迭代求解方程組。表4列出了程序主要并行部分的并行加速比和并行效率。由表4可知，求解器用于計(jì)算求解矩陣分配參數(shù)的時(shí)間較少，對(duì)于一個(gè)節(jié)點(diǎn)數(shù)為17萬的網(wǎng)格，計(jì)算時(shí)間不到1 s。計(jì)算單元?jiǎng)偠染仃?、組集總剛度矩陣時(shí)，由于并行程度較高，其加速效果比較理想。另外，求解的迭代次數(shù)隨著進(jìn)程數(shù)目的增加而變大。

圖3 不同網(wǎng)格規(guī)模下最大節(jié)點(diǎn)位移Fig.3 Maximum displacement of node for different scales of grid

求解器由于其總剛度矩陣的計(jì)算采取了精確預(yù)分配空間的方法，只存非零元素，故程序運(yùn)行時(shí)占用內(nèi)存較小。多個(gè)進(jìn)程運(yùn)行時(shí)，每個(gè)進(jìn)程占用內(nèi)存會(huì)明顯減小。

4.3 對(duì)于病態(tài)方程組的迭代收斂

在有限元中，若結(jié)構(gòu)材料參數(shù)相差過大會(huì)導(dǎo)致支配方程呈病態(tài)，即出現(xiàn)迭代次數(shù)過多或者迭代不收斂的問題。PETSc提供的迭代法的預(yù)條件子，能有效地解決或緩解這一問題。不考慮實(shí)際情況，將閘身的彈性模量作表5所示改變，即改變閘身彈性模量使之分別相差10倍、100倍至100萬倍。將這7組模型分別進(jìn)行4個(gè)進(jìn)程的并行迭代求解，得出的位移與Ansys軟件計(jì)算出的位移進(jìn)行對(duì)比(表6)。由表6可得，采用預(yù)條件子技術(shù)后，求解程序的迭代方法對(duì)因材料參數(shù)相差過大(最大相差10萬倍)而產(chǎn)生的病態(tài)方程組具有良好的適應(yīng)能力。

圖4 1～4個(gè)進(jìn)程計(jì)算時(shí)的并行加速比與并行效率Fig.4 Parallel speedup ratio and efficiency for one to four processes

表4 并行分析各部分的并行加速比和并行效率Table 4 Parallel speedup ratio and efficiency for various parts

表5 閘身彈性模量變化Table 5 Lock body elastic modulus change Pa

表6 修改材料參數(shù)后求解器的計(jì)算結(jié)果Table 6 Calculated results of solver after material parameters were modified

5 結(jié) 語

基于PETSc的數(shù)據(jù)結(jié)構(gòu)(矩陣Mat和向量Vec)和線性求解器(KSP)組件開發(fā)了有限元求解器，并以水閘在自重荷載下所產(chǎn)生的位移為算例，對(duì)該求解器進(jìn)行了測(cè)試。測(cè)試結(jié)果表明，該求解器能在普通的內(nèi)存為3G的四核計(jì)算機(jī)上，計(jì)算高達(dá)80萬節(jié)點(diǎn)的線彈性有限元問題，并行加速效果也較為明顯(4個(gè)進(jìn)程計(jì)算時(shí)加速比能達(dá)到3.24)。另外由于使用了PETSc提供的較為成熟的預(yù)條件子技術(shù)，求解器對(duì)病態(tài)方程組也有良好的適應(yīng)能力，為求解由于材料參數(shù)相差10萬倍而產(chǎn)生的病態(tài)方程組提供了新方法。由于現(xiàn)階段各方面的限制，本文開發(fā)的求解器僅僅使用了PETSc工具箱的很小一部分功能，程序也有很多地方需要改進(jìn)?；赑ETSc的有限元高性能求解方法的研究和開發(fā)還可以進(jìn)行以下幾方面的深入研究:

a.對(duì)程序進(jìn)一步完善擴(kuò)展，能夠?qū)€性有限元問題作進(jìn)一步的求解，如求出應(yīng)力、應(yīng)變。

b.將PETSc中的非線性求解器SNES和時(shí)間步進(jìn)器TS對(duì)應(yīng)到非線性有限元、時(shí)間依賴的動(dòng)力學(xué)問題的求解上。

c.本文所開發(fā)的求解器是一個(gè)沒有界面的程序，使用起來不方便，若能將其與圖像處理軟件結(jié)合起來或自己形成界面將會(huì)使該求解器更接近實(shí)際應(yīng)用。

［1］范穎，王磊，章青.多尺度有限元法及其應(yīng)用研究進(jìn)展［J］.水利水電科技進(jìn)展，2012，32(3):90-94.(FAN Ying，WANG Lei，ZHANG Qing.Research progress and application of multiscale finite element method［J］.Advances in Science and Technology of Water Resources，2012，32(3):90-94.(in Chinese))

［2］張汝清.概說并行計(jì)算結(jié)構(gòu)力學(xué)［J］.計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用，1995，12(4):447-484.(ZHANG Ruqing.Introduction to parallel computational structural mechanics［J］.Computational Structural Mechanics and Applications，1995，12(4):447-484.(in Chinese))

［3］張健飛.采用PETSc的有限元并行計(jì)算實(shí)現(xiàn)與優(yōu)化［J］.計(jì)算機(jī)工程與應(yīng)用，2010，46(10):57-59.(ZHANG Jianfei.Parallel implementation and tuning of finite element computing with PETSc［J］.Computer Engineering and Applications，2010，46(10):57-59.(in Chinese))

［4］姜弘道.水利高性能計(jì)算的進(jìn)展［J］.水利水電科技進(jìn)展，2006，26(2):70-76.(JIANG Hongdao.Progress in highperformance computing in water science［J］.Advances in Sciences and Technology of Water Resources，2006，26(2):70-76.(in Chinese))

［5］NOOR A K，VOIGT S J.Hypermatrix scheme for finite element system on CDC STAR-100 computer［J］.Computers and Structures，1975，5(5):287-296.

［6］ DIMITRIS Z.Parallel processing techniques for FE analysis:stiffness，loads and stress evaluation［J］.Computers and Structures，1988，28(2):247-260.

［7］ BITZARAKIS S，PAPADRAKAKIS M.Parallel solution techniques in computational structural mechanics［J］.Computer Methods in Applied Mechanics and Engineering，1997，148:75-104.

［8］CHIEN L S，SUNCT.Parallel processing techniques for finite element analysis of nonliear large truss structures［J］.Computers and Structures，1989，31(6):1023-1029.

［9］李麗君，金先龍，李淵印，等.一種新型并行化有限元結(jié)構(gòu)模態(tài)分析集成系統(tǒng)［J］.計(jì)算力學(xué)學(xué)報(bào)，2004，21(5):546-550.(LI Lijun，JIN Xianlong，LI Yuanyin，et al.A new parallel integration system of finite element model analysis［J］.Chinese Journal of Computational Mechanics，2004，21(5):546-550.(in Chinese))

［10］張林波，遲學(xué)斌，莫?jiǎng)t堯，等.并行計(jì)算導(dǎo)論［M］.北京:清華大學(xué)出版社，2006:1-487.

［11］Michael JQ.MPI與OpenMP并行程序設(shè)計(jì)［M］.北京:清華大學(xué)出版社，2004:34-40.

［12］都志輝.高性能計(jì)算并行編程技術(shù)MPI并行程序設(shè)計(jì)［M］.北京:清華大學(xué)出版社，2001:1-336.

［13］KRYLOV A N.On the numerical solution of the equation by which the frequency of small oscillations is determined in technical problems［J］.Izv Akad NaukAzerb Ser Fiz-Tekh Mat Nauk，1931，4:491-539.

［14］李曉梅，吳建平.Krylov子空間方法及其并行計(jì)算［J］.計(jì)算機(jī)科學(xué)，2005，32(1):19-20.(LI Xiaomei，WU Jianping.Krylovsubspace methods and parallel computation［J］.Computer Science，2005，32(1):19-20.(in Chinese))

［15］付成華，周洪波.面向?qū)ο蟮蔫F路路基加固設(shè)計(jì)系統(tǒng)集成與實(shí)現(xiàn)［J］.水利水電科技進(jìn)展，2012，32(4):70-73.(FU Chenghua，ZHOU Hongbo.Development and implementation of object-oriented railway subgrade strengthening design system［J］.Advances in Science and Technology of Water Resources，2012，32(4):70-73.(in Chinese))

［16］林江，母德偉，余葵，等.萬州長(zhǎng)江公路大橋防撞裝置結(jié)構(gòu)穩(wěn)定性及通航凈空尺度［J］.水利水電科技進(jìn)展，2013，33(2):59-62.(LIN Jiang，MU Dewei，YU Kui，et al.Structural stability and navigable clearance of anti-collision device of Wanzhou Yangtze River Highway Bridge［J］.Advances in Science and Technology of Water Resources，2013，33(2):59-62.(in Chinese))

［17］HESTENSEM R，STIEFEL E.Methods of conjugate gradients for solving linear systerms［J］.Tournal of Reaserch of the Natural Bureau of Standards，1952，49(6):409-436.