本發(fā)明屬于有限元結(jié)構(gòu)靜力學(xué)分析中的數(shù)值求解算法領(lǐng)域，具體涉及一種基于稀疏存儲(chǔ)格式的三維結(jié)構(gòu)靜力學(xué)仿真模擬方法。

背景技術(shù)：

目前，有限元靜力學(xué)分析的過(guò)程一般包括進(jìn)行有限元建模、計(jì)算單元?jiǎng)偠染仃?、集成總體剛度矩陣、計(jì)算載荷列陣、處理邊界條件、求解矩陣方程組和繪制場(chǎng)圖。這個(gè)分析過(guò)程的缺點(diǎn)在于總體剛度矩陣的集成和矩陣方程的求解將會(huì)產(chǎn)生大量的內(nèi)存消耗，同時(shí)時(shí)間性能較差。其實(shí)，在總體剛度矩陣中，大多數(shù)元素都是零元素，非零元素占有很小的比例，而在計(jì)算中有作用的正是這些非零元素，以及這些非零元素所處的位置。雖然現(xiàn)在有些有限元問(wèn)題的求解都采用了稀疏存儲(chǔ)，但是所采用的稀疏存儲(chǔ)的方式就是先生成總體剛度矩陣，然后從總體剛度矩陣中生成稀疏格式，是一種間接的稀疏求解方式，這樣做雖然能夠節(jié)省一定的內(nèi)存空間，但是并沒(méi)有從本質(zhì)上解決問(wèn)題。在文獻(xiàn)“稀疏矩陣存儲(chǔ)技術(shù)”[作者：張永杰，孫秦]中，描述了很多種稀疏格式存儲(chǔ)方法，但是并沒(méi)有講解如何從最初的模型數(shù)據(jù)直接到稀疏格式。在實(shí)際的結(jié)構(gòu)靜力學(xué)問(wèn)題的求解中，需要做到直接從有限元網(wǎng)格數(shù)據(jù)到稀疏存儲(chǔ)，才能全面提高時(shí)空效率，發(fā)揮有限元在求解結(jié)構(gòu)靜力學(xué)問(wèn)題中的優(yōu)勢(shì)。

技術(shù)實(shí)現(xiàn)要素：

針對(duì)背景技術(shù)中結(jié)構(gòu)靜力學(xué)問(wèn)題分析出現(xiàn)的不足和缺陷，本發(fā)明直接從靜力學(xué)模型出發(fā)，綜合考慮時(shí)間和空間因素，跳過(guò)總剛矩陣的集成，直接生成稀疏存儲(chǔ)格式。從而獲得一種高性能的結(jié)構(gòu)靜力學(xué)的仿真模擬方法，能夠從根本上解決內(nèi)存消耗大，同時(shí)時(shí)間性能較差的問(wèn)題。

本發(fā)明提供一種基于稀疏存儲(chǔ)格式的三維結(jié)構(gòu)靜力學(xué)仿真模擬方法，包含如下步驟：

a.處理有限元模型，得到編程所需要的網(wǎng)格數(shù)據(jù)，包括節(jié)點(diǎn)坐標(biāo)信息和節(jié)點(diǎn)編號(hào)信息。

b.對(duì)網(wǎng)格數(shù)據(jù)進(jìn)行計(jì)算處理，獲取該模型的棱邊信息。

c.劃定該模型的零位移邊界條件。

d.計(jì)算每個(gè)單元的單元?jiǎng)偠染仃?，并將這些結(jié)果保存起來(lái)。

e.判斷總體剛度矩陣的每一個(gè)節(jié)點(diǎn)，若該節(jié)點(diǎn)為0，則跳過(guò)；若不為0，則判斷哪些單元對(duì)該節(jié)點(diǎn)產(chǎn)生作用，然后就把相應(yīng)單元疊加到該節(jié)點(diǎn)，同時(shí)將非零元素以及位置信息寫入稀疏格式中。

f.按照稀疏格式求解矩陣，繪制場(chǎng)圖。

相比于傳統(tǒng)的有限元處理方式，本發(fā)明能夠節(jié)省大量的內(nèi)存，同時(shí)做到時(shí)間性能的優(yōu)化，由于跳過(guò)了集成總體剛度矩陣的過(guò)程，直接對(duì)稀疏格式進(jìn)行矩陣求解，使得計(jì)算規(guī)?？梢詷O大的增加，求解矩陣的規(guī)?？梢赃_(dá)到上百萬(wàn)階，能夠滿足工程中的絕大多數(shù)問(wèn)題。

附圖說(shuō)明

圖1基于稀疏格式的三維結(jié)構(gòu)靜力學(xué)仿真模擬方法流程圖。

圖2四面體網(wǎng)格剖分模型效果。

圖3二次函數(shù)四面體單元模型。

圖4點(diǎn)基函數(shù)示意圖。

圖5邊基函數(shù)示意圖。

圖6三維空間問(wèn)題示例模型。

具體實(shí)施方式

下面結(jié)合附圖和具體的實(shí)例來(lái)說(shuō)明本發(fā)明的技術(shù)方案。

a.處理有限元模型，得到編程所需要的網(wǎng)格相關(guān)數(shù)據(jù)，包括節(jié)點(diǎn)坐標(biāo)信息和節(jié)點(diǎn)編號(hào)信息。

對(duì)目標(biāo)模型做如圖2所示的四面體網(wǎng)格剖分，一個(gè)四面體為一個(gè)單元，產(chǎn)生了大量的節(jié)點(diǎn)，通過(guò)計(jì)算得到這些節(jié)點(diǎn)的空間坐標(biāo)數(shù)據(jù)，并對(duì)這些節(jié)點(diǎn)進(jìn)行順序編號(hào)處理。

在函數(shù)的階次上，選擇二階函數(shù)，二階函數(shù)四面體單元的模型如圖3所示。對(duì)函數(shù)的階次做進(jìn)一步說(shuō)明，一階函數(shù)四面體單元只有四個(gè)角節(jié)點(diǎn)，二次函數(shù)四面體單元有十個(gè)節(jié)點(diǎn)，除了角上的四個(gè)節(jié)點(diǎn)，每條邊的中點(diǎn)也有一個(gè)節(jié)點(diǎn)，在二次函數(shù)四面體單元中，共有三十個(gè)節(jié)點(diǎn)位移自由度，設(shè)其節(jié)點(diǎn)位移列陣為

q^e＝(u1υ1w1,.....,u10υ10w10)^t(1)

其中uk,υk,wk,(k＝1,2,3,4,5,6,7,8,9,10)為x,y,z三個(gè)方向上的位移。

由于是二階函數(shù)，則取單元的位移模式為

其中ak,(k＝1,2,3,4,5,6,7,8,9,10)為該函數(shù)的系數(shù)，此位移模式為(x,y,z)的二次函數(shù)。

b.對(duì)網(wǎng)格數(shù)據(jù)進(jìn)行計(jì)算處理，獲取該模型的棱邊信息。

對(duì)步驟a中得到的數(shù)據(jù)進(jìn)行計(jì)算，獲取棱邊的信息，并對(duì)棱邊順序編號(hào)處理。

在步驟b中，需要獲取模型的棱邊信息，是因?yàn)楸痉椒ú捎玫氖嵌A函數(shù)四面體單元，不需要用到邊節(jié)點(diǎn)(在邊上的節(jié)點(diǎn))的坐標(biāo)數(shù)據(jù)，僅需要編號(hào)。介于上述原因，該方法統(tǒng)計(jì)了模型的棱邊編號(hào)信息。這樣做的好處是減小了剖分模型的計(jì)算壓力，提高了計(jì)算效率。

c.劃定該模型的零位移邊界條件。

在對(duì)實(shí)際問(wèn)題進(jìn)行靜力學(xué)分析時(shí)，模型的某些位置總是會(huì)處于固定狀態(tài)，處在固定位置上的點(diǎn)必須設(shè)置其位移為0。

d.計(jì)算每個(gè)單元的單元?jiǎng)偠染仃?，并將這些計(jì)算結(jié)果保存起來(lái)。

在步驟d中，方法采取的是將所有單元的單元矩陣先進(jìn)行計(jì)算并保存，而不是在步驟e中進(jìn)行計(jì)算，如果采取臨時(shí)計(jì)算單元矩陣的方式，一個(gè)單元的單元矩陣可能會(huì)多次重復(fù)計(jì)算，不利于提高時(shí)間效率，先將結(jié)果進(jìn)行保存的方式雖然浪費(fèi)了小部分內(nèi)存，但是換取了良好的時(shí)間效率。

在采用了式(2)位移模式之后，通過(guò)體積坐標(biāo)來(lái)構(gòu)建形函數(shù)，將位移模式寫為

uk,υk,wk,(k＝1,2,3,4,5,6,7,8,9,10)為x,y,z三個(gè)方向上的位移，用矩陣可以表示為

u＝nq^e(4)

根據(jù)形函數(shù)具有的性質(zhì)，形函數(shù)nk(x,y,z)＝δij,(i,j＝1,2,3,4)，可以得到式(5)中形函數(shù)

其中角節(jié)點(diǎn)上的形函數(shù)為

ni＝(2li-1)li,i＝1,2,3,4(6)

邊節(jié)點(diǎn)上的形函數(shù)為

式中

其中，v表示當(dāng)前四面體單元的體積

在式(9)和式(10)中，xi,yi,zi(i＝1,2,3,4)分別表示了四面體單元四個(gè)角節(jié)點(diǎn)的空間坐標(biāo)。

最后通過(guò)彈性力學(xué)的平衡方程，幾何方程，物理方程來(lái)構(gòu)建單元矩陣，這個(gè)過(guò)程為本領(lǐng)域人員公知常識(shí)，不在此詳述。

e.判斷總體剛度矩陣的每一個(gè)節(jié)點(diǎn)，如果該節(jié)點(diǎn)為0，則不考慮該節(jié)點(diǎn)，如果不為0，則判斷哪些單元對(duì)該節(jié)點(diǎn)產(chǎn)生作用，然后就把相應(yīng)單元疊加到該節(jié)點(diǎn)。

步驟a,b,c,d為本領(lǐng)域公知常識(shí)，本發(fā)明的主要?jiǎng)?chuàng)新點(diǎn)在于步驟e。通過(guò)步驟d，可以計(jì)算任意一個(gè)單元的單元?jiǎng)偠染仃嚕∈璐鎯?chǔ)需要知道兩個(gè)重要的因素，一個(gè)是總剛矩陣的每一個(gè)節(jié)點(diǎn)是否為0，第二個(gè)是如果不為0，哪些單元會(huì)對(duì)該節(jié)點(diǎn)產(chǎn)生貢獻(xiàn)。為了方便敘述，以一個(gè)二維模型為例說(shuō)明，因?yàn)槿S空間問(wèn)題完全一致，首先考慮角節(jié)點(diǎn)。有限元模型如圖4所示。

得到這個(gè)模型后，可以非常輕易的知道此模型的節(jié)點(diǎn)數(shù)目和單元數(shù)目，從而可以知道總體剛度矩陣的規(guī)模，本方法只需要知道總剛矩陣的規(guī)模，而不需要去計(jì)算它。按照傳統(tǒng)的有限元總剛矩陣的生成方法，在先考慮角節(jié)點(diǎn)的情況下，單元矩陣為3*3，則可以得到這個(gè)模型的總體剛度矩陣如下

對(duì)于該總剛矩陣的(1,1)位置，從圖3可以發(fā)現(xiàn)，第一個(gè)節(jié)點(diǎn)僅處在第一個(gè)單元的范圍內(nèi)，就是說(shuō)只有第一個(gè)單元對(duì)該節(jié)點(diǎn)有貢獻(xiàn)，所以可以確定的是只有第一個(gè)單元的單元?jiǎng)偠染仃嚂?huì)疊加到該位置，現(xiàn)在需要進(jìn)一步確認(rèn)是第一個(gè)單元?jiǎng)偠染仃嚨哪囊粋€(gè)元素會(huì)加進(jìn)來(lái)，將四個(gè)單元的單元矩陣和全局編號(hào)羅列出來(lái)。

在總體剛度矩陣的(1,1)位置，按照第一個(gè)單元所有節(jié)點(diǎn)全局編號(hào)的順序，將第一個(gè)單元的(1,1)位置加到總剛矩陣的(1,1)位置，值得說(shuō)明的是，單元?jiǎng)偠染仃嚨?1,1)位置是按照全局編號(hào)來(lái)的，而不是按照行列標(biāo)來(lái)的。

對(duì)于總剛矩陣的(1,2)位置，雖然節(jié)點(diǎn)2被多個(gè)單元共用，但是節(jié)點(diǎn)1只處于單元1中，所以總剛矩陣的(1,2)位置也僅僅只有單元1有貢獻(xiàn)，所以同樣按照全局編號(hào)的方式，將單元1的(1,2)位置疊加到總剛矩陣的(1,2)位置。

對(duì)于總剛矩陣的(1,4)位置，節(jié)點(diǎn)1處于單元1，而節(jié)點(diǎn)4處于單元2，這兩個(gè)節(jié)點(diǎn)都沒(méi)有同時(shí)被某個(gè)單元共用，所以該位置必然為0，所以在這個(gè)位置就不用計(jì)算。

在總剛矩陣的(2,2)位置，這個(gè)位置只考察節(jié)點(diǎn)2，節(jié)點(diǎn)2被單元1、單元2，單元3三個(gè)單元共用，所以必須按照全局編號(hào)的方式，分別將單元1、單元2、單元3三個(gè)矩陣的(2,2)位置疊加到總剛矩陣的(2,2)位置。

在總剛矩陣的(2,3)位置，節(jié)點(diǎn)2和節(jié)點(diǎn)3同時(shí)處于單元1和單元3中，所以在該位置單元1和單元3都有貢獻(xiàn)，按照全局編號(hào)的方式將單元1和單元3的(2,3)位置上的元素疊加進(jìn)來(lái)。

在進(jìn)行單元判斷的過(guò)程中，還需要處理位移邊界條件，本方法采取了對(duì)角線元素置1法來(lái)處理邊界，如果發(fā)現(xiàn)此節(jié)點(diǎn)位于邊界上，則將對(duì)應(yīng)的位置置為1。

所以，按照上述方法就可以知道總剛矩陣的哪些位置為零，哪些是非零，哪些單元會(huì)對(duì)該位置產(chǎn)生貢獻(xiàn)。方法最終的目的是生成標(biāo)準(zhǔn)稀疏格式，標(biāo)準(zhǔn)稀疏格式只存儲(chǔ)矩陣的非零元素。一般用一個(gè)一維數(shù)組存儲(chǔ)矩陣的非零元素，用另一個(gè)索引數(shù)組存儲(chǔ)這些非零元素在原矩陣中的行和列值。例如

通過(guò)對(duì)上面的矩陣進(jìn)行計(jì)算，可以得到標(biāo)準(zhǔn)稀疏格式如下

a＝[1.02.03.04.05.06.07.08.09.010.011.0]

ia＝[11223344555]

ja＝[12243524145]

標(biāo)準(zhǔn)稀疏格式中，a表示的是非零元素，ia表示的每一個(gè)非零元素的行標(biāo)，ja表示的非零元素的列標(biāo)。所以在對(duì)每一個(gè)位置進(jìn)行判斷計(jì)算的過(guò)程中，將非零的元素及其行、列標(biāo)等信息存入稀疏格式中。即可完成稀疏格式的高效生成。

考慮邊節(jié)點(diǎn)的情況，在圖5二維的情況下，邊節(jié)點(diǎn)上除了邊界上的點(diǎn)外，其他的始終都被兩個(gè)單元共用，只需要按照同樣的方式將單元矩陣的相應(yīng)節(jié)點(diǎn)疊加上來(lái)。在三維空間問(wèn)題中，如圖5、6所示，處理的方式方法和二維情況下完全一樣，只是單元?jiǎng)偠染仃囈?guī)模大了很多，但是疊加方式不變。

f.按照稀疏格式求解矩陣，繪制場(chǎng)圖。

通過(guò)一般的矩陣求解方法來(lái)求解上述矩陣，便可以得到位移解，這里不在詳述，然后根據(jù)求解結(jié)果畫出場(chǎng)圖分布。

按照本方法來(lái)處理一個(gè)力學(xué)模型，假設(shè)模型規(guī)模為3141個(gè)單元，761個(gè)點(diǎn)，4339條邊，總共存在66378個(gè)非零元素，如果按照常規(guī)方式來(lái)處理，會(huì)產(chǎn)生1.74gb的內(nèi)存消耗，而采用了本方法做稀疏存儲(chǔ)之后，僅僅會(huì)產(chǎn)生800kb的內(nèi)存消耗，兩者相差幾千倍，相比于傳統(tǒng)的靜力學(xué)仿真方法，本發(fā)明具有明顯的優(yōu)勢(shì)。