王忠雷，趙國群，馬新武

(1.山東大學材料液固結(jié)構(gòu)演變與加工教育部重點實驗室，濟南250061;2.山東建筑大學機電工程學院，濟南250101)

三維有限元剛度矩陣的壓縮存儲算法

王忠雷1，2，趙國群1，馬新武1

(1.山東大學材料液固結(jié)構(gòu)演變與加工教育部重點實驗室，濟南250061;2.山東建筑大學機電工程學院，濟南250101)

為提高有限元分析效率、減少存儲空間消耗，對剛度矩陣的壓縮存儲算法進行了研究.研究了“廣義相鄰節(jié)點對”與剛度矩陣中非零子矩陣的關(guān)系，確定了剛度矩陣中非零子矩陣的分布規(guī)律;提出了一種新的剛度矩陣壓縮存儲方法—“改進的CSR存儲方法”，給出了基于壓縮存儲的剛度矩陣的生成過程以及線性方程組迭代解法方法，并將提出的算法應(yīng)用于三維體積成形有限元分析軟件.有限元分析實例表明，該算法可以有效地減少存儲空間，提高計算效率.

三維有限元方法;剛度矩陣;壓縮存儲;六面體網(wǎng)格

隨著計算機技術(shù)和有限元技術(shù)的發(fā)展，有限元模擬在金屬塑性加工中獲得了廣泛應(yīng)用，并發(fā)揮了不可替代的作用［1-3］.近年來，所要分析的工程問題日益復雜、數(shù)據(jù)量大，尤其是在大型三維體積成形過程數(shù)值分析中，工件形狀的復雜性導致單元模型節(jié)點數(shù)量的激增，從而使得有限元軟件在分析該類問題時計算效率低下，甚至在分析大型、超大型問題時出現(xiàn)存儲空間不足的情況，這大大限制了有限元軟件在大型體積成形分析中的應(yīng)用.

工程問題的有限元分析最終都歸結(jié)為求解大型線性方程組，有限元的求解效率在很大程度上取決于線性方程組的解法及剛度矩陣的存儲結(jié)構(gòu).有限元剛度矩陣一般為大型、帶狀稀疏矩陣，矩陣的非零元素一般集中在對角線區(qū)域，為此出現(xiàn)了很多節(jié)約存儲空間的存儲方法，如半帶存儲、變帶寬存儲等.這些方法可以減少剛度矩陣的存儲空間，但都不可避免地存儲了部分非零元素，浪費了儲存空間，而且受到節(jié)點編號的影響，需要對節(jié)點編號進行優(yōu)化，浪費大量的計算時間.為了進一步節(jié)約存儲空間，鮑益東等［4］針對板料成形問題提出了一種適合迭代算法的改進一維變帶寬壓縮存儲方法，該方法采用一維數(shù)組按行的順序保存剛度矩陣的非零元素，同時采用2個輔助數(shù)組記錄列號和每列起始元素在剛度矩陣數(shù)組的位置.該方法節(jié)約了大量存儲空間，也避免了節(jié)點編號的優(yōu)化難題，但該方法采用的輔助數(shù)組與剛度矩陣存儲數(shù)組具有同樣長度，仍存在消耗部分存儲空間的問題.

本文根據(jù)剛度矩陣的組成特點，提出了一種新的剛度矩陣壓縮存儲方法，該方法在輔助數(shù)組中只記錄非零子矩陣的位置，間接記錄剛度矩陣非零元素的位置，進一步節(jié)約了剛度矩陣的存儲空間，尤其適合于求解大型、復雜的三維工程問題.

1 剛度矩陣非零元素確定

在有限元網(wǎng)格模型中，同一單元內(nèi)的任意2個節(jié)點構(gòu)成的節(jié)點對(包括節(jié)點自身與自身構(gòu)成的節(jié)點對)稱為“廣義相鄰節(jié)點對”.根據(jù)剛度矩陣的生成原理，有限元剛度矩陣中的非零元素與“廣義相鄰節(jié)點對”有著密切的關(guān)系，每組“廣義相鄰節(jié)點對”對應(yīng)一個與所分析問題維數(shù)相同的非零子矩陣，非零子矩陣的位置由廣義節(jié)點對的全局節(jié)點序號決定.例如同一個單元中的2個節(jié)點的全局節(jié)點編號為(I，J)，構(gòu)成一個“廣義相鄰節(jié)點對”Bij，該“廣義相鄰節(jié)點對”產(chǎn)生的非零子矩陣在廣義矩陣中的位置為I行J列.

如果考慮“廣義相鄰節(jié)點對”的順序，即Bij和Bji作為2個不同的節(jié)點對，它們分別對應(yīng)第I行、第J列和第J行、第I列的非零子矩陣，有序“廣義相鄰節(jié)點對”與整個剛度矩陣的非零子矩陣數(shù)量相同，且一一對應(yīng).由于剛度矩陣為對稱矩陣，為節(jié)約存儲空間，可只存儲上三角矩陣，因此只需考慮Bij(I≤J)的節(jié)點對對應(yīng)的非零子矩陣，不需考慮節(jié)點對的順序問題.Bij和Bji作為同一個節(jié)點對，無序廣義節(jié)點對與上三角矩陣存儲的剛度矩陣的非零子矩陣數(shù)量相同，一一對應(yīng)，因此，完全可以通過無序廣義節(jié)點對及其數(shù)量來確定剛度矩陣中非零子矩陣的數(shù)量和位置.

對于如圖1所示的八節(jié)點六面體單元，其“廣義相鄰節(jié)點對”分為3個類型:一是邊相鄰節(jié)點對，例如節(jié)點對12、14、15;二是面相鄰節(jié)點對，例如節(jié)點對16、18、13;三是體相鄰節(jié)點對，例如節(jié)點對17.從整個單元模型來看，邊相鄰節(jié)點對由單元內(nèi)相鄰節(jié)點連線的數(shù)量決定，每條連線對應(yīng)1個無序節(jié)點對;面相鄰節(jié)點對的數(shù)量由單元面的數(shù)量決定，每個面有2條對角線，每條面對角線對應(yīng)1個無序節(jié)點對;體相鄰節(jié)點對的數(shù)量由單元的數(shù)量決定，每個單元有4條體對角線，每條體對角線對應(yīng)1個無序節(jié)點對.

圖1所示單元模型的無序廣義節(jié)點對的數(shù)量即非零子矩陣的數(shù)量，可由式(1)確定:

式中:NS為網(wǎng)格模型無序節(jié)點對數(shù)量;N1為節(jié)點數(shù)量，表示節(jié)點與自身生成的節(jié)點對;N2為同一單元內(nèi)相鄰節(jié)點連線的數(shù)量，表示邊相鄰節(jié)點對的數(shù)量;N3為單元的面數(shù)，2N3表示面相鄰節(jié)點對的數(shù)量;N4為單元數(shù)量，4N4表示體相鄰節(jié)點對的數(shù)量.

圖1 八節(jié)點六面體單元模型

對角線非零子矩陣的數(shù)量也為N1，每個對角線非零子矩陣只需存儲上三角子矩陣和對角線元素，若設(shè)N為所求解問題的維數(shù)，則所需存儲的非零元素為(N2+N)/2.對于三維問題N=3，則每個對角線非零子矩陣只需存儲6個數(shù)據(jù)，其他非零子矩陣需要存儲N2=9個數(shù)據(jù)，因此，總共需要保存的非零剛度矩陣元素數(shù)量為

式中，NZ為網(wǎng)格模型非零剛度矩陣元素數(shù)量，N1、N2、N3、N4與式(1)中的意義相同.

同理，對于如圖2所示的由2個單元構(gòu)成的簡單六面體網(wǎng)格，N1=12、N2=20、N3=11、N4= 2.采用上三角矩陣存儲時，其非零子矩陣數(shù)量為

總共需要保存的剛度矩陣非零元素數(shù)量為

一旦確定了非零子矩陣和非零元素的數(shù)量，即可準確地分配剛度矩陣的存儲空間，這對于剛度矩陣存儲算法的實現(xiàn)具有重要意義.

圖2 八節(jié)點六面體單元網(wǎng)格模型

2 剛度矩陣的存儲方案

稀疏矩陣有多種存儲方式，在此采用了比較常用的 CSR格式(Compressed Sparse Row format).該方法采用3個一維數(shù)組保存稀疏矩陣，浮點型數(shù)組A按行順序存儲矩陣中的各非零元素;整形數(shù)組J中記錄對應(yīng)元素所在的列號;整形數(shù)組K中記錄每行第1個元素在A中的位置.其中A和J的長度相同，即為稀疏矩陣中非零元素的數(shù)量，而整形數(shù)組K的長度則為稀疏矩陣的維數(shù).例如，對于圖3所示的稀疏矩陣，采用CSR格式的存儲格式如下所示:

由于剛度矩陣為對稱矩陣，可以只存儲其上三角矩陣，因此K數(shù)組應(yīng)保存對角線元素在數(shù)組A中的位置.CSR方法可以很大程度上減少剛度矩陣的存儲空間，而且按行引用，其效率較高，比較適合存儲有限元的剛度方程.但是整型數(shù)組J與浮點數(shù)組A長度一致，占用了大約1/3的存儲空間，造成了存儲效率的降低.由前面的分析可以看出，有限元剛度矩陣中非零元素的分布是有規(guī)律的，只分布于與“廣義相鄰節(jié)點對”對應(yīng)的非零子矩陣中，因此，只要記錄了非零子矩陣的位置，即可間接地確定非零元素的位置，該方法稱為改進的CSR方法.對于三維問題，非零子矩陣有9個元素，只記錄非零子矩陣的位置，就可以減少數(shù)組J的存儲空間約90%，進而使整個剛度矩陣的存儲空間減少約30%.在此方法中，整形數(shù)組J中記錄對應(yīng)非零子矩陣所在的列號;整形數(shù)組K中記錄廣義矩陣中每行第一個非零子矩陣在J中的位置.

下面以如圖2所示的單元模型的剛度矩陣為例進行分析，采用CSR方法的J數(shù)組數(shù)據(jù)長度為所有非零元素的數(shù)量，即NZ=522，K數(shù)組數(shù)據(jù)長度為剛度方程的維數(shù)，即為12×3=36，因此，總數(shù)據(jù)量為558;而采用改進的CSR方法后，J數(shù)組數(shù)據(jù)長度則為剛度矩陣非零子矩陣數(shù)量，即NS= 62，K數(shù)組數(shù)據(jù)長度為單元數(shù)量12，則總數(shù)據(jù)量減少為74.如果只比較2個輔助數(shù)組，則改進的CSR方法所節(jié)約存儲空間為如果浮點型數(shù)組A采用單精度浮點數(shù)保存，則改進的CSR方法后，整個剛度矩陣所節(jié)約的存儲空間為

顯而易見，改進的CSR方法進一步明顯地節(jié)約了剛度矩陣的存儲空間，但該方法間接表示非零元素的位置，增加了訪問剛度矩陣的計算時間消耗，對計算效率有一定的影響，屬于以時間換空間的方法.因此建議在單元數(shù)量較多時采用改進的CSR方法，以減少內(nèi)存消耗;而單元數(shù)量較少時，則采用前一種方法，提高對剛度矩陣的訪問效率，達到最好的綜合效益.對于三維體積成形問題，內(nèi)存不足是分析的主要問題，因此，宜采用能夠節(jié)約內(nèi)存的改進的CSR方法.

3 剛度矩陣的生成

在生成剛度矩陣時，首先根據(jù)單元模型的拓撲結(jié)構(gòu)來確定非零子矩陣的數(shù)量和位置，生成索引數(shù)組J、K;然后根據(jù)剛度矩陣的生成規(guī)則和2個索引數(shù)組的數(shù)據(jù)，將單元剛度矩陣的數(shù)據(jù)對號入座地進行組合.對于改進的CSR存儲方法，還需要對索引數(shù)據(jù)進行轉(zhuǎn)換.下面通過一個實例說明剛度矩陣的生成過程.

假設(shè)一個三維八節(jié)點六面體單元，其單元節(jié)點的整體編號為I［8］，單元剛度矩陣為U［8× 3］［8×3］，單元剛度矩陣的任意元素可描述為

式中:i、j=1、2、Λ、8為節(jié)點局部編號;k、l為物理量維數(shù)編號，k、l=1、2、3.其在整體剛度矩陣數(shù)組中的序號按如下步驟確定:

1)確定Uijkl的整體節(jié)點編號為

式中，I、J為Uijkl所在的非零子矩陣AIJ中的廣義行號和列號，非零子矩陣的位置如圖3所示.

圖3 非零子矩陣分布示意圖

2)確定非零子矩陣AII之前非零子矩陣的數(shù)量NI:

3)確定非零子矩陣AII之前非零元素的數(shù)量NNI:式中，(I-1)×3是由于對角非零子矩陣引起的非零元素的減少項.如圖4顯示的AII子矩陣只存儲了6個非零元素，與其他子矩陣相比少存了3個.

圖4 非零元素位置示意圖

4)確定第I行非零子矩陣中第k行之前非零元素的數(shù)量NNIK(未考慮對角非零矩陣引起的非零元素的減少項):

5)搜索數(shù)組J從標號K［I］到K［I+1］中J［II］=J的數(shù)組標號II，則第I行中，第J列之前的非零子矩陣數(shù)量NJ為

確定NJ的列查詢是整個剛度矩陣生成過程中消耗時間最多的步驟，而在此步驟當中，普通的CSR方法按元素搜索 ，而改進CSR方法按照子矩陣進行搜索.在三維問題中，普通CSR方法的搜索時間為該進CSR方法搜索時間的3倍，因此，改進的CSR方法在此步驟中較常規(guī)CSR方法更節(jié)約時間.

6)確定 Uijkl在剛度矩陣中所在行中的編號NNJl:

式中M為由于對角非零子矩陣引起的非零元素的減少項，k=1、2、3時M=0、1、3.

7)確定 Uijkl在剛度矩陣數(shù)組 A中的編號NNIJkl:

一旦確定了Uijkl在剛度矩陣數(shù)組A中的位置后，即可根據(jù)剛度矩陣的生成規(guī)則進行累加，從而獲得總體剛度矩陣.

4 線性方程組的求解

大型有限元方程組的求解方法一般分為直接法和迭代法2種.所謂直接法是指在不考慮舍入誤差存在的情況下，經(jīng)過有限步運算就可以求得方程組的精確解的方法，此方法也稱為精確法.而迭代法是一種逐次逼近的方法，它從一個初始向量出發(fā)，通過一個給定的計算格式，構(gòu)造出一個無窮向量列，此向量列即是方程組的近似解，其極限為方程組的精確解.直接法求解精度高、穩(wěn)定性好，但是在方程數(shù)量較多時，求解效率低，因此，在求解大型線性方程組時，往往采用迭代法［5-8］.

線性方程組的迭代方法有很多，其中松弛預處理共額梯度法(SSOR-PCG)是求解對稱正定系數(shù)矩陣最受歡迎的方法之一.SSOR-PCG方法的基本原理如下［9］.

設(shè)線性方程組(11)的系數(shù)矩陣A為n階對稱正定矩陣，

設(shè)M=(STS)-1為對稱正定矩陣，則上式等價于

式中A'=SA ST，b'=Sb，u'=S-Tu.

M通常取為對稱逐步超松弛迭代法(SSOR法)的分裂矩陣，即

式中:0＜ω＜2，為松弛因子;D為A的對角陣;L為A的嚴格下三角陣，即

在改進的SSOR-PCG方法中，不必對M求逆，也不必另開辟數(shù)組存放M矩陣，計算中僅對A的非零元素進行運算，所以僅需存儲A的下三角陣(或上三角陣)中的非零元素，因此，該方法特別適合于本文中只存儲非零子矩陣的剛度方程.

5 算例分析與討論

為了驗證本文所提出的算法在節(jié)約內(nèi)存、提高計算效率方面的效果，將剛度矩陣壓縮存儲方法應(yīng)用于自行開發(fā)的三維體積成形有限元分析程序CasForm中［10］.以如圖5所示的道釘成形為例，從內(nèi)存占用和計算效率兩個方面，將本文提出的壓縮存儲、迭代求解的算法與半帶存儲、直接求解的算法進行了比較.

圖5 道釘成形模擬結(jié)果示意圖

在內(nèi)存占用方面，以八節(jié)點六面體單元模型為例，按照剛度方程采用單精度浮點數(shù)保存，即每個數(shù)據(jù)占用4個字節(jié);輔助數(shù)組采用整形數(shù)據(jù)保存，即每個數(shù)據(jù)占用2個字節(jié).在不同的節(jié)點和單元數(shù)量下，對采用不同方式存儲的剛度矩陣所占存儲空間進行比較，結(jié)果見表1.

表1 不同存儲方法存儲空間的比較

從表1可以看出，壓縮存儲方法所占的存儲空間遠遠小于半帶存儲方法的存儲空間，存儲空間節(jié)約超過90%，而且隨著節(jié)點數(shù)量的增加，空間節(jié)約的比例逐漸增大.這主要是因為壓縮存儲方法只保存了剛度矩陣的非零元素，剛度矩陣中每一行非零元素的數(shù)量僅與該行對應(yīng)節(jié)點的廣義相鄰節(jié)點的數(shù)量有關(guān)，而與整體模型節(jié)點的數(shù)量無關(guān).換言之，對于壓縮存儲方法，存儲空間與節(jié)點數(shù)量的正比關(guān)系僅與剛度矩陣的行數(shù)增加有關(guān)，呈線性變化趨勢.而對于半帶存儲方法，每一行的存儲空間由剛度矩陣的半帶寬決定，隨節(jié)點數(shù)量增加，半帶寬一般也隨之增加，因此，隨節(jié)點數(shù)量增加，不僅剛度矩陣的行數(shù)增加，而且每行需要存儲的數(shù)據(jù)也有所增加，因而呈現(xiàn)拋物線變化趨勢.由此可知，隨著節(jié)點數(shù)量的增加，半帶存儲方法所需的存儲空間呈二次冪指數(shù)函數(shù)增加，而采用壓縮存儲方法，其節(jié)約空間的比例更加顯著.

在計算效率方面，以求解單個方程的計算效率為例，對半帶存儲直接求解和壓縮存儲迭代求解兩種方法的進行了比較，結(jié)果見表2.從表2可以看出，壓縮存儲迭代求解的效率遠遠高于半帶存儲直接求解的效率，而且節(jié)點數(shù)量越多，效率提高越明顯.

表2 不同單元模型不同存儲方法求解效率的比較

6 結(jié)論

1)依據(jù)由網(wǎng)格拓撲結(jié)構(gòu)決定的“廣義相鄰節(jié)點對”，可以確定有限元剛度矩陣中非零子矩陣的分布和數(shù)量.

2)采用改進的CSR存儲方法保存有限元剛度矩陣，可顯著減少剛度矩陣的存儲空間，而且節(jié)點數(shù)量越多，空間節(jié)約的比例越大.

3)采用基于壓縮矩陣的SSOR-PCG迭代方法求解有限元線性方程組，可顯著提高方程組的求解效率，而且節(jié)點數(shù)量越多，求解效率提高的幅度越大.

［1］ GHAEI A，MOVAHHEDY R M.Die design for the radial forging process using 3D FEM［J］.Journal of Materials Processing Technology，2007，182:534-539.

［2］ BEHRENS B A.Finite element analysis of die wear in hot forging processes［J］.CIRP Annals-Manufacturing Technology，2008，57:305-308.

［3］ LU C，ZHANG L W，MU Z J，et al.3D FEM simulation of the multi-stage forging process of a gas turbine compressor blade［J］.Journal of Materials Processing Technology，2008，198:463-470.

［4］ 鮑益東，杜國康，陳文亮，等.沖壓成形模擬中有限元方程組求解算法［J］.南京航空航天大學學報，2009，41(5):606-610. BAO Yi-dong，DU Guo-kang，CHEN Wen-liang，et al. Solving algorithm for finite element equations in sheet metal forming simulation［J］.Journal of Nanjing University of Aero nautics＆Astronautics，2009，41(5): 606-610.

［5］ 張永杰.大型稀疏方程組預處理迭代快速求解技術(shù)研究［D］.西安:西北工業(yè)大學，2006.

［6］ 劉艷芳，施法中，徐向陽.板料沖壓成形數(shù)值模擬中有限元方程求解算法的研究［J］.塑性工程學報，2006，13(4):15-19. LIU Yan-fang，SHI Fa-zhong，XU Xiang-yang.Study on algorithms of solving FE equations in the numerical simulation of sheet metal stamping［J］.Journalof Plasticity Engineering，2006，13(4):15-19.

［7］ 李靜，陳健，云周晶.改進的SSOR-PCG迭代法在接觸問題研究中應(yīng)用［J］.大連理工大學學報，2006，46(4):533-537. LI Jing，CHEN Jian，YUN Zhou-jing.Application of an improved SSOR-PCG iteration method to contact problems［J］.Journal of Dalian University of Technology， 2006，46(4):533-537.

［8］ 包勁青，楊強，陳英儒，等.對稱超松弛預處理共軛梯度法在高拱壩整體大規(guī)模彈塑性有限元分析中的應(yīng)用［J］.水利學報，2009，40(5):589-594. BAO Jin-qing，YANG Qiang，CHEN Ying-ru，et al.Application of symmetric successive over relaxation precondition conjugate method to large-scale elasto-plastic FEM analysis for high arch dams［J］.Journal of Hydraulic Engineering，2009，40(5):589-594.

［9］ 林紹忠.對稱逐步超松弛預處理共軛梯度法的改進迭代格式［J］.數(shù)值計算與計算機應(yīng)用，1997(4): 266-270. LIN Shao-zhong.Improved iterative format of symmetric successive over relaxation-preconditioned conjugated gradient method［J］.Numerical Methods and Computer Applications，1997(4):266-270.

［10］ 王忠雷，趙國群，黃麗麗，等.六面體網(wǎng)格體積成形有限元分析關(guān)鍵技術(shù)［J］.材料科學與工藝，2010，18(4):509-513. WANG Zhong-lei，ZHAO Guo-qun，HUANG Li-li，et al.Key technologies of 3D FEM analysis system for bulk forming with hexahedral mesh［J］.Materials Science＆Technology，2010，18(4):509-513.

Compressed storage algorithm of 3D-FEM stiffness matrix

WANG Zhong-lei1，2，ZHAO Guo-qun1，MA Xin-wu1
(1.Key Laboratory for Liquid-Solid Structural Evolution and Processing of Materials(Ministry of Education)，Shandong University，Jinan 250061，China;2.Institute of Electrical and Mechanical，Shandong Architectural University，Jinan 250101，China)

To improve the efficiency and reduce the storage space of finite element analysis，compression and storage algorithm of 3D-FEM stiffness matrix is studied.The relationship between“generalized adjacent double nodes"and the non-zero sub-matrix in stiffness matrix is researched for getting distribution of non-zero sub-matrix in stiffness matrix.A new algorithm of stiffness matrix of compressed storage-”improved CSR storage method"is proposed.Based on the algorithm，the generation process of stiffness matrix is given and iterative solution of linear equations method is proposed to improve the efficiency of solving linear equations.The algorithm is applied to the three-dimensional bulk forming finite element analysis software and the numerical results show that the algorithm can effectively decrease the storage space and improve the computation efficiency.

3D-FEM;stiffness matrix;compressed storage;hexahedral mesh

TG316.8 文獻標志碼：A 文章編號：1005-0299(2012)02-0096-05

2011-02-21.

國家自然科學基金資助課題(50875155);教育部“長江學者和創(chuàng)新團隊發(fā)展計劃”創(chuàng)新團隊資助項目(IRT0931).

王忠雷(1977-)，男，博士研究生;

趙國群(1962-)，男，教授，博士生導師.

趙國群，E-mail:zhaogq@sdu.edu.cn.

(編輯 程利冬)