專利名稱：用于獲取工程結(jié)構(gòu)靜動(dòng)態(tài)力學(xué)特性的數(shù)值復(fù)合單元方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種用于獲取工程結(jié)構(gòu)靜動(dòng)態(tài)力學(xué)特性的數(shù)值方法，屬于結(jié)構(gòu)力學(xué)技術(shù)領(lǐng)域。
現(xiàn)有的用于工程結(jié)構(gòu)靜動(dòng)態(tài)力學(xué)特性分析的技術(shù)主要為有限元方法，它由于具有很好的邊界適應(yīng)性而得到大量和廣泛的應(yīng)用，但目前大都采用低階多項(xiàng)式形狀函數(shù)，因而計(jì)算精度不高。當(dāng)然可以通過兩種途徑來提高現(xiàn)行有限元方法的計(jì)算精度，第一是通過密化有限元網(wǎng)格(即h-version方法)，第二是采用高階次的多項(xiàng)式形狀函數(shù)(即p-version方法和等級(jí)有限元方法hierarchical FEM)，近年來，發(fā)展了綜合前兩種方法的自適應(yīng)h-p version方法，這些方法雖然可以有效地提高計(jì)算精度，但卻使計(jì)算量大大地增加，由于所構(gòu)造的高階形狀函數(shù)仍然是采用多項(xiàng)式(如Legendre正交多項(xiàng)式函數(shù))，因而構(gòu)造過程復(fù)雜、效率很低，同時(shí)，也存在數(shù)值穩(wěn)定性方面的問題。在結(jié)構(gòu)的動(dòng)態(tài)特性分析過程中，這些缺點(diǎn)則更加突出。
在結(jié)構(gòu)分析中，對(duì)于一些情形簡(jiǎn)單而幾何形狀又規(guī)整的構(gòu)件，也可采用經(jīng)典力學(xué)方法，即通過求解微分方程來得到相應(yīng)的解析解或精度很高的近似解，這些簡(jiǎn)單構(gòu)件包括桿的拉伸、桿的扭轉(zhuǎn)、梁的彎曲、方或圓板的平面變形和彎曲變形等，邊界條件包括簡(jiǎn)支、固支，由于幾何形狀及邊界條件的嚴(yán)格限制，經(jīng)典力學(xué)方法在實(shí)際工程結(jié)構(gòu)中的應(yīng)用極為有限。
本發(fā)明之目的在于1.提出復(fù)合單元方法的概念。它可使傳統(tǒng)的有限元方法和經(jīng)典力學(xué)解析解進(jìn)行復(fù)合，從而繼承各自方法的優(yōu)點(diǎn)，具有傳統(tǒng)有限元方法優(yōu)良的離散特性，即可對(duì)任意幾何形狀的結(jié)構(gòu)進(jìn)行離散化，又可繼承經(jīng)典力學(xué)解析解的高精度，因而，既具有很好的適應(yīng)性，又能達(dá)到很高的精度，從而實(shí)現(xiàn)高效率。
2.定義兩組自由度座標(biāo)體系來描述復(fù)合單元。第一組為節(jié)點(diǎn)自由度座標(biāo)體系，它主要在于運(yùn)用傳統(tǒng)有限元方法，目的在于對(duì)復(fù)雜幾何形狀進(jìn)行離散化，第二組為單元的場(chǎng)自由度座標(biāo)體系(叫作c-自由度)，它主要在于運(yùn)用經(jīng)典力學(xué)的解析解，其目的在于提高計(jì)算精度，從而實(shí)現(xiàn)高效率。
3.提出“零邊界條件”下解析解概念。即由經(jīng)典力學(xué)在零邊界條件下求得解析解，把它作為在復(fù)合單元方法中所使用的場(chǎng)自由度座標(biāo)體系的基底函數(shù)，它是兩組自由度座標(biāo)體系容合的關(guān)鍵，也是傳統(tǒng)有限元法與經(jīng)典力學(xué)解析解復(fù)合的基礎(chǔ)，本發(fā)明申請(qǐng)中可提供適用于各種邊界性要求(如C0、C1問題)的“零邊界條件”下解析解。
4.基于所提出的兩組自由度座標(biāo)體系，定義新的形狀函數(shù)描述(叫作復(fù)合形狀函數(shù))，它由兩部分組成，其一為傳統(tǒng)有限元方法中的形狀函數(shù)，其二為經(jīng)典力學(xué)在“零邊界條件”下解析解得到的形狀函數(shù)，同樣，定義新的復(fù)合自由度，它也是由節(jié)點(diǎn)自由度和場(chǎng)自由度組成。
5.在建立好復(fù)合單元的復(fù)合形狀函數(shù)后，可完全按照傳統(tǒng)有限元方法中的過程推導(dǎo)出各種單元(拉伸桿單元、扭轉(zhuǎn)桿單元、梁?jiǎn)卧?、平面三?jié)點(diǎn)單元、平面四節(jié)點(diǎn)單元、空間四節(jié)點(diǎn)四面體單元、空間八節(jié)點(diǎn)六面體單元、三節(jié)點(diǎn)彎曲板單元、四節(jié)點(diǎn)彎曲板單元......等)的剛度矩陣及質(zhì)量矩陣。
6.復(fù)合單元方法的實(shí)施與傳統(tǒng)有限元方法類似，如果將其c-自由度忽略或取為零，則復(fù)合單元方法將退化為傳統(tǒng)的有限元方法，也就是說，傳統(tǒng)有限方法是復(fù)合單元方法的一種特殊情形。
7.在復(fù)合單元方法中，可以通過兩種途徑來提高計(jì)算精度，其一為與傳統(tǒng)有限元方法類似的密化離散網(wǎng)格，(即h-version方法)，其二為增加c-自由度(叫作c-version方法)，第二種途徑是實(shí)現(xiàn)復(fù)合單元高精度的主要手段。
本發(fā)明用于獲取工程結(jié)構(gòu)靜動(dòng)態(tài)力學(xué)特性的數(shù)值復(fù)合單元方法，包括下述步驟(1)在對(duì)結(jié)構(gòu)進(jìn)行離散化得到各種單元后，定義兩組自由度座標(biāo)體系來描述單元的位移場(chǎng)，一組為節(jié)點(diǎn)自由度座標(biāo)體系，另一組為場(chǎng)自由度座標(biāo)體系(叫作c-自由度)。(2)基于節(jié)點(diǎn)自由度座標(biāo)體系中的節(jié)點(diǎn)自由度q，構(gòu)造與傳統(tǒng)有限元方法相同的位移場(chǎng)函數(shù)UFEM(ξ)及相應(yīng)的形狀函數(shù)N(ξ)UFEM(ξ)＝N(ξ)q(3)基于場(chǎng)自由度座標(biāo)體系中的場(chǎng)自由度c，構(gòu)造由經(jīng)典力學(xué)理論在“零邊界條件”下得到的解析解作為基底的位移場(chǎng)函數(shù)UCT(ξ)及相應(yīng)的形狀函數(shù)φ(ξ)UCT(ξ)＝φ(ξ)c(4)由步驟(2)及(3)得到的位移場(chǎng)函數(shù)復(fù)合相加而形成復(fù)合單元的復(fù)合位移場(chǎng)函數(shù)U(ξ)、復(fù)合形狀函數(shù)S(ξ)及復(fù)合自由度δ
U(ξ)＝UFEM(ξ)＋UCT(ξ)＝N(ξ)q＋φ(ξ)c＝S(ξ)δ其中S(ξ)＝[N(ξ)φ(ξ)]，δ＝[qc]T(5)在得到復(fù)合單元的復(fù)合形狀函數(shù)S(ξ)后，可計(jì)算與應(yīng)變?chǔ)?ξ)對(duì)應(yīng)的幾何矩陣B(ξ)ε(ξ)＝B(ξ)δ(6)然后計(jì)算復(fù)合單元的剛度矩陣ke及質(zhì)量矩陣meKe=∫VBTDBdV]]>me=∫VρSTSdV]]>其中D為彈性矩陣，ρ為密度。(7)對(duì)結(jié)構(gòu)各單元的剛度矩陣和質(zhì)量矩陣進(jìn)行從局部座標(biāo)到總體座標(biāo)的轉(zhuǎn)換，然后再進(jìn)行組裝以得到整體剛度矩陣和整體質(zhì)量矩陣K=Σe=1nTeTkeTe]]>M=Σe=1nTeTmeTe]]>其中Te為座標(biāo)轉(zhuǎn)換矩陣。(8)有了復(fù)合單元的整體剛度矩陣和質(zhì)量矩陣后，可對(duì)工程結(jié)構(gòu)進(jìn)行靜力學(xué)分析Kδ＝P和動(dòng)力學(xué)分析Kδ＝ω2Mδ其中P為外載荷向量，ω為結(jié)構(gòu)的共振自然頻率。(9)在復(fù)合單元方法中，有兩種途徑來提高數(shù)值分析精度，一種是通過與傳統(tǒng)有限元方法類似的密化離散網(wǎng)格，叫作h-version方法，另一種是增加復(fù)合單元的場(chǎng)自由度，即增加c-自由度，叫作c-version方法，后一種方法的效率遠(yuǎn)高于前一種。(10)傳統(tǒng)有限元方法可看作為本發(fā)明申請(qǐng)的復(fù)合單元方法的一種特殊情況，即將復(fù)合單元方法中的場(chǎng)自由度忽略或取為零，則復(fù)合單元方法退化為傳統(tǒng)有限元方法。
本發(fā)明的方法可用于各種結(jié)構(gòu)的靜動(dòng)力學(xué)分析、穩(wěn)定性分析以及相應(yīng)的CAD、FEM軟件中。采用復(fù)合單元方法可使數(shù)值分析效率大大提高，只需很少的幾何離散單元通過增加c-自由度就能獲得很高的計(jì)算精度并達(dá)到一種超級(jí)收斂。


圖1表示本發(fā)明的實(shí)施流程框圖。圖2表示復(fù)合單元方法中的軸向單元。圖3表示復(fù)合單元方法中的扭轉(zhuǎn)單元。圖4表示復(fù)合單元法中的彎曲單元。圖5A表示一由7根軸向桿組成的桁架。圖5B表示一與x軸成γ角的任意桿單元。圖6表示一右端固定的懸臂梁。圖7表示各種方法用于計(jì)算左端固支軸向拉伸桿第4階特征頻率的相對(duì)誤差的比較。圖8表示各種方法用于計(jì)算左端固支懸臂梁第4階特征頻率的相對(duì)誤差的比較。
下面結(jié)合附圖詳細(xì)介紹本發(fā)明的內(nèi)容。
1.對(duì)工程結(jié)構(gòu)(或連續(xù)體)進(jìn)行幾何離散，如離散為桿單元、梁?jiǎn)卧鍐卧?、平面三?jié)點(diǎn)單元、空間八節(jié)點(diǎn)單元……等。
2.位移場(chǎng)函數(shù)的表達(dá)，在復(fù)合單元方法中，可定義兩組自由度座標(biāo)體系來描述單元的位移場(chǎng)函數(shù)，即U(ξ)＝UFEM(ξ)＋UCT(ξ)其中U(ξ)為單元的復(fù)合位移場(chǎng)函數(shù)，UFEM(ξ)為基于節(jié)點(diǎn)自由度座標(biāo)體系的位移場(chǎng)函數(shù)，UCT(ξ)為基于單元場(chǎng)自由度座標(biāo)體系的位移場(chǎng)函數(shù)，ξ為無量綱幾何位置座標(biāo)。如果將UFEM(ξ)取為傳統(tǒng)有限元方法中的節(jié)點(diǎn)位移場(chǎng)函數(shù)，UCT(ξ)取為由“零邊界條件”下解析解得到的位移場(chǎng)函數(shù)，則有UFEM(ξ)＝N(ξ)qUCT(ξ)＝φ(ξ)c其中N(ξ)為傳統(tǒng)有限元方法中的節(jié)點(diǎn)形狀函數(shù)，q為節(jié)點(diǎn)自由度，φ(ξ)為由經(jīng)典力學(xué)的“零邊界下解析解”，c為場(chǎng)自由度(即c-自由度)。所以，復(fù)合單元的位移場(chǎng)函數(shù)可表達(dá)為
U(ξ)＝N(ξ)q＋φ(ξ)c＝S(ξ)δ其中，S(ξ)＝[N(ξ)φ(ξ)]叫作復(fù)合單元的復(fù)合形狀函數(shù)，δ＝[qc]T叫作復(fù)合單元的復(fù)合自由度。可以看出S(ξ)和δ都是由兩個(gè)部分組成，即基于節(jié)點(diǎn)自由度座標(biāo)體系的傳統(tǒng)有限元部分和基于場(chǎng)自由度座標(biāo)體系的經(jīng)典解析解部分。值得注意的是，這里的經(jīng)典解析解是基于“零邊界條件”下得到的，復(fù)合單元方法中所說的“零邊界條件”具體可表達(dá)為對(duì)于C0連續(xù)問題，“零邊界條件”為φr(ξ)|ξ＝0＝0，φr(ξ)|ξ＝1＝0r＝1，2，3，…其中φr(ξ)為場(chǎng)函數(shù)(在工程結(jié)構(gòu)中，φr(ξ)為位移場(chǎng)函數(shù))。對(duì)于C1連續(xù)問題，“零邊界條件”為φr(ξ)|ξ＝0＝0，φr(ξ)|ξ＝1＝0dφr(ξ)dξ|ξ=0=0,dφr(ξ)dξ|ξ=1=0]]>r＝1，2，3，…對(duì)于任意的Ck連續(xù)問題，“零邊界條件”為φr(ξ)|ξ＝0＝0，φr(ξ)|ξ＝1＝0dφr(ξ)dξ|ξ=0=0,dφr(ξ)dξ|ξ=1=0]]>dkφr(ξ)dξk|ξ=0=0,dkφr(ξ)dξk|ξ=1=0]]>r＝1，2，3，…以上的“零邊界條件”將是獲取經(jīng)典解析解的關(guān)鍵，它將作為基于場(chǎng)自由度座標(biāo)體系的形狀函數(shù)，在大多數(shù)情形下，各種單元的φr(ξ)都可以通過在“零邊界條件”下求解微分方程而得到解析表達(dá)式。3.單元?jiǎng)偠燃百|(zhì)量矩陣的求取，這一步驟與傳統(tǒng)的有限元方法相同，由于復(fù)合單元方法中的形狀函數(shù)包括有兩個(gè)部分，因此所得到的單元?jiǎng)偠燃百|(zhì)量矩陣與傳統(tǒng)有限元方法中的單元?jiǎng)偠燃百|(zhì)量矩陣不同，后面將具體給出幾種類型的單元?jiǎng)偠燃百|(zhì)量矩陣。
4.將各單元的剛度及質(zhì)量矩陣進(jìn)行裝配可形成總體的剛度及質(zhì)量矩陣，然后按傳統(tǒng)有限元的步驟進(jìn)行方程求解則可得到結(jié)構(gòu)的位移、應(yīng)變、應(yīng)力、自然頻率、振型……等。
下面進(jìn)一步詳述復(fù)合單元方法中各種類型單元的位移函數(shù)、形狀函數(shù)、剛度矩陣以及質(zhì)量矩陣的表達(dá)。(1)軸向桿單元圖2為一具有軸向位移的桿單元，這是一個(gè)C0連續(xù)問題，即只要求位移函數(shù)連續(xù)，q1和q2表示節(jié)點(diǎn)自由度座標(biāo)體系中的節(jié)點(diǎn)自由度，c1c2......cr表示場(chǎng)自由度座標(biāo)體系中的場(chǎng)自由度(叫作c-自由度)，l為單元的長(zhǎng)度，ξ＝x/l為局部的無量綱幾何座標(biāo)，基于節(jié)點(diǎn)自由度座標(biāo)體系，軸向位移函數(shù)可表達(dá)為UFEM(ξ)＝q1(1－ξ)＋q2ξ＝N(ξ)q其中 N(ξ)＝[(1－ξ)ξ]，qT＝[q1q2]基于場(chǎng)自由度座標(biāo)體系的位移函數(shù)可表達(dá)為UCT(ξ)＝c1φ1(ξ)＋c2φ2(ξ)＋…＋crφ(ξ)＝φ(ξ)cφi(ξ)，i＝1，2，......r為基底函數(shù)，由經(jīng)典力學(xué)理論就桿的拉壓?jiǎn)栴}在“零邊界條件”下所求解而得到，即φ1(ξ)＝sin(πξ)φ2(ξ)＝sin(2πξ)φr(ξ)＝sin(rπξ)所以φ(ξ)＝[φ1(ξ)φ2(ξ)…φr(ξ)]＝[sin(πξ)sin(2πξ)…sin(rπξ)]cT＝[c1c2…cr]那么，桿單元的復(fù)合位移場(chǎng)為U(ξ)＝UFEM(ξ)＋UCT(ξ)＝q1(1－ξ)＋q2ξ＋c1sin(πξ)＋c2sin(2πξ)＋…＋crsin(rπξ)＝S(ξ)δ其中S(ξ)＝[N(ξ)φ(ξ)]＝[(1－ξ)ξsin(πξ)sin(2πξ)…sin(rπξ)]叫作桿單元的復(fù)合形狀函數(shù)。
δ＝[qTcT]＝[q1q2c1c2…cr]T叫作桿單元的復(fù)合自由度。
桿單元的軸向應(yīng)變?yōu)?amp;epsiv;=∂U(x)∂x=1l∂U(ξ)∂ξ]]>=1l∂S(ξ)∂ξ•δ=B(ξ)•δ]]>其中B(ξ)=1l∂S(ξ)∂ξ]]>=1l[-11πcos(πξ)2πcos(2πξ)…rπcos(rπξ)]]]>桿單元的剛度矩陣為KLe=∫VBTEBdV]]> q1q2c1c2… … cr E是楊氏模量，A為桿的橫截面積，V為單元的體積。桿單元的質(zhì)量矩陣為mLe=∫VρSTSdV]]> q1q2c1c2… … cr (2)扭轉(zhuǎn)桿單元圖3為一具有扭轉(zhuǎn)位移的桿單元，這是一個(gè)C0連續(xù)問題，即只要求位移函數(shù)連續(xù)，w1及w2表示節(jié)點(diǎn)自由度座標(biāo)體系中的節(jié)點(diǎn)自由度，c1、c2、....cr表示場(chǎng)自由度座標(biāo)體系中的場(chǎng)自由度(叫作c-自由度)，l為單元的長(zhǎng)度，ξ＝x/l為局部的無量綱幾何坐標(biāo)?；诠?jié)點(diǎn)自由度座標(biāo)體系的扭轉(zhuǎn)位移函數(shù)可表達(dá)為ΘFEM(ξ)＝w1(1－ξ)＋w2ξ＝N(ξ)q其中N(ξ)＝[(1－ξ)ξ]，q＝[w1w2]T基于場(chǎng)自由度座標(biāo)體系的位移函數(shù)可表達(dá)為ΘCT(ξ)＝c1φ1(ξ)＋c2φ2(ξ)＋…＋crφr(ξ)＝φ(ξ)c其中φi(ξ)，i＝1，2，.....r為基底函數(shù)，它由經(jīng)典力學(xué)理論就桿的扭轉(zhuǎn)問題在“零邊界條件”下求解而得到，即φ1(ξ)＝sin(πξ)φ2(ξ)＝sin(2πξ)φr(ξ)＝sin(rπξ)所以φ(ξ)＝[sin(πξ)sin(2πξ)…sin(rπξ)]c＝[c1c2…cr]那么，桿單元的復(fù)合扭轉(zhuǎn)位移為Θ(x)＝ΘFEM(x)＋ΘCT(x)＝w1(1－ξ)＋w2ξ＋c1sin(πξ)＋c2sin(2πξ)＋…＋crsin(rπξ)＝N(ξ)q＋φ(ξ)c＝S(ξ)·δ其中S(ξ)＝[N(ξ)φ(ξ)]＝[(1－ξ)ξsin(πξ)sin(2πξ)…sin(rπξ)]叫作扭轉(zhuǎn)單元的復(fù)合形狀函數(shù)。
δ＝[q c]T＝[w1w2c1c2…cn]T叫作扭轉(zhuǎn)桿單元的復(fù)合自由度。
扭轉(zhuǎn)桿單元的剪切應(yīng)變?yōu)閞θx=y∂Θ(x)∂x]]>=yl∂S(ξ)∂ξ•δ=B(ξ)]]>其中y為距中性軸的距離，B(ξ)為B(ξ)=yl∂S(ξ)∂ξ]]>=yl[-11πcos(πξ)2πcos(2πξ)…rπcos(rπξ)]]]>扭轉(zhuǎn)桿單元的剛度矩陣為KTe=∫VBTGBdV]]> w1w2c1c2… … cr G為彈性剪切模量，J為扭轉(zhuǎn)桿橫截面的極慣性矩，V為單元的體積。
扭轉(zhuǎn)桿單元的質(zhì)量矩陣為mTe=∫VρSTSdV]]>
w1w2c1c2… …cr (3)彎曲梁?jiǎn)卧獔D4為一具有橫向位移的彎曲梁?jiǎn)卧?，這是一個(gè)C1連續(xù)問題，即要求撓度的一階導(dǎo)數(shù)連續(xù)，v1和v2表示節(jié)點(diǎn)自由度座標(biāo)體系中的節(jié)點(diǎn)橫向撓度位移，θ1及θ2表示節(jié)點(diǎn)自由度座標(biāo)體系中的節(jié)點(diǎn)彎曲轉(zhuǎn)角，c1，c2，....，cr表示場(chǎng)自由度座標(biāo)體系中的場(chǎng)自由度(叫做c-自由度)，l為單元的長(zhǎng)度，ξ＝x/l為局部的無量綱幾何座標(biāo)，基于節(jié)點(diǎn)自由度座標(biāo)體系的彎曲，撓度函數(shù)可表達(dá)為WFEM(ξ)＝v1(1－3ξ2＋2ξ3)＋θ1l(ξ－2ξ2＋ξ3)＋v2(3ξ2－2ξ3)＋θ2l(ξ3－ξ2)＝N(ξ)q其中N(ξ)＝[(1－3ξ2＋2ξ3)l(ξ－2ξ2＋ξ3)(3ξ2－2ξ3)l(ξ3－ξ2)]q＝[v1θ1v2θ2]T基于場(chǎng)自由度座標(biāo)體系，彎曲撓度位移函數(shù)可表達(dá)為WCT(ξ)＝c1φ1(ξ)＋c2φ2(ξ)＋…＋crφr(ξ)＝φ(ξ)cφi(ξ)，i＝1，2，.....r為基底函數(shù)，它由經(jīng)典力學(xué)理論就梁的彎曲問題在“零邊界條件”下求解而得到，即φ1(ξ)=F1(λ1*,ξ)]]>φ2(ξ)=F2(λ2*,ξ)]]>φr(ξ)=Fr(λr*,ξ)]]>其中Fi(λi*,ξ)=sinλi*ξ-sinhλi*ξ-(sinλi*-sinhλi*cosλi*-coshλi*)(cosλi*ξ-cosλi*ξ)]]>i＝1，2，3，…并且λi*滿足關(guān)系cosλi*•coshλi*=1]]>i＝1，2，3，…即λ1*=40730041]]>λ2*=7.853205]]>λ3*=10.995608]]>λ4*=14.137165]]>λr*=(r+0.5)π,]]>r＞4所以φ(ξ)=[F1(λ1*,ξ)F2(λ2*,ξ)…Fr(λr*,ξ)]]]>c=[c1c2…cr]T那么，梁?jiǎn)卧膹?fù)合撓度位移為W(ξ)＝WFEM(ξ)+WCT(ξ)＝v1(1－3ξ2＋2ξ3)+θ1l(ξ－2ξ2+ξ3)＋v2(3ξ2－2ξ3)＋θ2l(ξ3－ξ2)＋c1F1(λ1*,ξ)+c2F2(λ2*,ξ)+…+crFr(λr*,ξ)]]>＝N(ξ)q＋φ(ξ)c＝S(ξ)·δ其中S(ξ)＝[N(ξ)φ(ξ)]＝[(1－3ξ2＋2ξ3)l(ξ－2ξ2＋ξ3)(3ξ2－2ξ3)l(ξ3－ξ2)F1(λ1*,ξ)F2(λ2*,ξ)…Fr(λr*,ξ)]]]>叫作梁?jiǎn)卧膹?fù)合形狀函數(shù)。
δ＝[qTcT]＝[v1θ1v2θ2c1c2…cr]T叫作梁?jiǎn)卧膹?fù)合自由度。
梁?jiǎn)卧膽?yīng)變?yōu)?amp;epsiv;=-y-∂2W(x)∂x2]]>=-y-l2∂2S(ξ)∂ξ2•δ=B(ξ•δ)]]>其中y為距中性軸距離，B(ξ)為B(ξ)=-y-[1l2(12ξ-6)1l(6ξ-4)-1l2(12ξ-6)1l(6ξ-2)]]>F1''(λ1*,ξ)F2''(λ2*,ξ)…Fr''(λr*,ξ)]]]>Fr''(λr*,ξ)=]]>-λr*2[sinλr*ξ+sinhλr*ξ-(sinλr*-sinhλr*cosλr*-coshλr*)cosλr*ξ+coshλr*ξ]]>r＝1，2，3，…ξ=xl]]>λ1*=4.730041]]>λ2*=7.853205]]>λ3*=10.995608]]>λ4*=14.137165]]>λr*=(r+0.5)π,]]>r＞4彎曲梁?jiǎn)卧膭偠染仃嚍閗Be=∫VBTDBdV]]>=E∫A∫0lBTBdxdA]]>
v1θ1v2θ2c 其中，I為梁?jiǎn)卧獧M截面的慣性矩，V為單元體積，子矩陣kBcc為c1c2c3c4c5… cr 彎曲梁?jiǎn)卧馁|(zhì)量矩陣為mBe=∫VρSTSdV]]>=ρA∫0lSTSdx]]>v1θ1v2θ2c ·V為單元體積，子矩陣mBqc為v1θ1v2θ2 子矩陣mBcc為c1c2c3c4c5… cr (4)其它的復(fù)合單元對(duì)于平面三節(jié)點(diǎn)單元、平面四節(jié)點(diǎn)單元、空間四節(jié)點(diǎn)四面體單元、空間八節(jié)點(diǎn)六面體單元、三節(jié)點(diǎn)彎曲板單元、四節(jié)點(diǎn)彎曲板單元等各種類型的復(fù)合單元以及對(duì)應(yīng)的等參復(fù)合單元，完全可以按照前面所述的方法構(gòu)造其相應(yīng)的基于節(jié)點(diǎn)自由度座標(biāo)體系和場(chǎng)自由度座標(biāo)體系的復(fù)合位移函數(shù)、復(fù)合形狀函數(shù)、應(yīng)變、剛度矩陣以及質(zhì)量矩陣。(5)應(yīng)用實(shí)例圖5A表示一個(gè)由7根桿組成的桁架?，F(xiàn)分析它的各階共振頻率和模態(tài)，圖中的參數(shù)如下L＝2l＝4m，h＝2m，橫截面積A＝0.001m2，桿的質(zhì)量密度ρ＝8000Ns2/m4，楊氏模量E＝2.1*1011N/m2。
圖5B為一與x軸成γ角的任意桿單元，[qiqj]為節(jié)點(diǎn)的局部軸向位移，[uiviujvj]為整體節(jié)點(diǎn)位移，它們之間的轉(zhuǎn)換關(guān)系為qi＝uicosγ＋visinγqj＝ujcosγ＋vjsinγ取任意復(fù)合桿單元的局部“復(fù)合自由度”δe為δe＝[qiqjc1c2…cr]T取任意復(fù)合桿單元的整體“復(fù)合自由度”δe為δe＝[uiviujvjc1c2…cr]T那么δe與δe之間的轉(zhuǎn)換關(guān)系為δe＝Teδe其中uiviujvjc1c2…… cr 那么，在整體座標(biāo)體系中，桿單元的剛度矩陣為Ke＝TeTkeTe
uiviujvjc1…cr 桿單元的質(zhì)量矩陣為Me＝TeTmeTeuiviujvjc1… cr 分別對(duì)圖5A中桁架中的7個(gè)桿單元形成各自的基于整體座標(biāo)的單元?jiǎng)偠染仃嚭唾|(zhì)量矩陣，然后將它們組裝，形成總剛度矩陣和質(zhì)量矩陣K=Σe=17Ke]]>M=Σe=17Me]]>再求解以下振動(dòng)方程Kδ＝ω2Mδ如果每個(gè)桿單元只取2個(gè)場(chǎng)自由度(c-自由度)，則求解的前10階自然頻率為ω1＝1648.26Hz，ω2＝1741.32Hz，ω3＝3113.83Hz，ω4＝4567.69Hz，ω5＝4829.70Hz，ω6＝7379.96Hz，ω7＝7532.30Hz，ω8＝8047.93Hz，ω9＝9997.48Hz，ω10＝10567.43Hz，如果用傳統(tǒng)的有限元方法求解此問題，則得到的前6階自然頻率為ω1＝1683.52Hz，ω2＝1776.28Hz，ω3＝3341.37Hz，ω4＝5174.35Hz，ω5＝5678.1 8Hz，ω6＝8315.40Hz，下面再給出一個(gè)梁?jiǎn)卧獞?yīng)用的例子。
圖6所示為一右端固定的懸臂梁，L為梁的長(zhǎng)度，ρ為密度，E為楊氏模量，I為橫截面的慣性矩，現(xiàn)分析它的各階共振頻率和模態(tài)。
如果只用一個(gè)復(fù)合梁?jiǎn)卧獊硖幚磉@一問題，取4個(gè)場(chǎng)自由度(c-自由度)，它的總剛度矩陣K為v1θ1c1c2c3c4 總質(zhì)量矩陣為v1θ1c1c2c3c4 令λi4=ρL4EIωi2]]>i＝1，2，…其中ωi為共振自然頻率。
由復(fù)合單元方法計(jì)算出的各階λi如下λ1＝1.875109，λ2＝4.694419，λ3＝7.857543，λ4＝11.00451，由傳統(tǒng)有限元方法取3個(gè)單元計(jì)算出的各階λi如下λ1＝1.875199，λ2＝4.701793，λ3＝7.903542，λ4＝11.86048，由解析解得到的各階精確的λi如下λ1＝1.875104，λ2＝4.694091，λ3＝7.854757，λ4＝10.99554，圖7為各種方法用于計(jì)算左端固定拉伸桿的第4階特征頻率(λ4)的誤差比較，圖8為各種方法用于計(jì)算左端固定懸臂梁的第4階特征頻率(λ4)的誤差比較，可以看出，在相同計(jì)算量的前提下，復(fù)合單元方法的計(jì)算結(jié)果的誤差大大低于常規(guī)有限元方法的誤差，其計(jì)算效率高出10倍左右，取較多場(chǎng)自由度(c-自由度)的結(jié)果明顯優(yōu)于取較少場(chǎng)自由度的結(jié)果。
權(quán)利要求
1.一種用于獲取工程結(jié)構(gòu)靜動(dòng)態(tài)力學(xué)特性的數(shù)值復(fù)合單元方法，其特征在于該方法包括下述步驟(1)在對(duì)結(jié)構(gòu)進(jìn)行離散化得到各種單元后，定義兩組自由度座標(biāo)體系來描述單元的位移場(chǎng)，一組為節(jié)點(diǎn)自由度座標(biāo)體系，另一組為場(chǎng)自由度座標(biāo)體系(叫作c-自由度)；(2)基于節(jié)點(diǎn)自由度座標(biāo)體系中的節(jié)點(diǎn)自由度q，構(gòu)造與傳統(tǒng)有限元方法相同的位移場(chǎng)函數(shù)UFEM(ξ)及相應(yīng)的形狀函數(shù)N(ξ)UFEM(ξ)＝N(ξ)q(3)基于場(chǎng)自由度座標(biāo)體系中的場(chǎng)自由度c，構(gòu)造由經(jīng)典力學(xué)理論在“零邊界條件”下得到的解析解作為基底的位移場(chǎng)函數(shù)UCT(ξ)及相應(yīng)的形狀函數(shù)φ(ξ)UCT(ξ)＝φ(ξ)c(4)由步驟(2)及(3)得到的位移場(chǎng)函數(shù)復(fù)合相加而形成復(fù)合單元的復(fù)合位移場(chǎng)函數(shù)U(ξ)、復(fù)合形狀函數(shù)S(ξ)及復(fù)合自由度δU(ξ)＝UFEM(ξ)＋UCT(ξ)＝N(ξ)q＋φ(ξ)c＝S(ξ)δ其中S(ξ)＝[N(ξ)φ(ξ)]，δ＝[qc]T(5)在得到復(fù)合單元的復(fù)合形狀函數(shù)S(ξ)后，可計(jì)算與應(yīng)變?chǔ)?ξ)對(duì)應(yīng)的幾何矩陣B(ξ)ε(ξ)＝B(ξ)δ(6)然后計(jì)算復(fù)合單元的剛度矩陣ke及質(zhì)量矩陣meKe=∫VBTDBdV]]>me=∫VρSTSdV]]>其中D為彈性矩陣，ρ為密度；(7)對(duì)結(jié)構(gòu)各單元的剛度矩陣和質(zhì)量矩陣進(jìn)行從局部座標(biāo)到總體座標(biāo)的轉(zhuǎn)換，然后再進(jìn)行組裝以得到整體剛度矩陣和整體質(zhì)量矩陣K=Σe=1nTeTkeTe]]>M=Σe=1nTeTmeTe]]>其中Te為座標(biāo)轉(zhuǎn)換矩陣；(8)有了復(fù)合單元的整體剛度矩陣和質(zhì)量矩陣后，可對(duì)工程結(jié)構(gòu)進(jìn)行靜力學(xué)分析Kδ＝P和動(dòng)力學(xué)分析Kδ＝ω2Mδ其中P為外載荷向量，ω為結(jié)構(gòu)的共振自然頻率；(9)在復(fù)合單元方法中，有兩種途徑來提高數(shù)值分析精度，一種是通過與傳統(tǒng)有限元方法類似的密化離散網(wǎng)格，叫作h-version方法，另一種是增加復(fù)合單元的場(chǎng)自由度，即增加c-自由度，叫作c-version方法，后一種方法的效率遠(yuǎn)高于前一種；(10)傳統(tǒng)有限元方法可看作為本發(fā)明申請(qǐng)的復(fù)合單元方法的一種特殊情況，即將復(fù)合單元方法中的場(chǎng)自由度忽略或取為零，則復(fù)合單元方法退化為傳統(tǒng)有限元方法。
全文摘要
本發(fā)明涉及一種用于獲取工程結(jié)構(gòu)靜動(dòng)態(tài)力學(xué)特性的數(shù)值方法。該方法是在對(duì)結(jié)構(gòu)進(jìn)行離散后，定義兩組自由度坐標(biāo)體系來描述單元的位移場(chǎng)，基于節(jié)點(diǎn)坐標(biāo)體系，應(yīng)用常規(guī)插值多項(xiàng)式構(gòu)造出位移場(chǎng)函數(shù)U
文檔編號(hào)G06F17/00GK1139242SQ9610046
公開日1997年1月1日 申請(qǐng)日期1996年2月2日 優(yōu)先權(quán)日1996年2月2日
發(fā)明者曾攀 申請(qǐng)人:曾攀